time | calls | line |
---|
| | 1 | function [YYt ui vi Pi rhoi] = DOorthnorm(app, dx, dy, YYt, ui, uid, vi, vid, Pi, pid, rhoi)
|
| | 2 | % function [YYt u v P rho a] = DOorthnorm(YYt, ui, uid, vi, vid, Pi, pid, rhoi)
|
| | 3 | % This function orthonormalizes the DO modes while preserving the original
|
| | 4 | % realizations (upto a correction for the stochastic energy) and the
|
| | 5 | % stochastic energy.
|
| | 6 | %
|
| | 7 | % INPUTS:
|
| | 8 | % dx: Size of x discretization
|
| | 9 | % dy: Size of y discretization
|
| | 10 | % YYt(MC,S): Stochastic coefficients
|
| | 11 | % ui(:,S): Mode shapes for u-velocity
|
| | 12 | % uid: Vector of ids for interior u-velocities
|
| | 13 | % vi(:,S): Mode shapes for v-velocity
|
| | 14 | % vid: Vector of ids for interior v-velocities
|
| | 15 | % Pi(:,S): Mode shapes for Pressure
|
| | 16 | % Pid: Vector of ids for interior Pressures
|
| | 17 | % rhoi(Optional): Mode shapes for density
|
| | 18 | %
|
| | 19 | % OUTPUTS:
|
| | 20 | % YYt(MC,S): Stochastic coefficients
|
| | 21 | % ui(:,S): Mode shapes for u-velocity
|
| | 22 | % vi(:,S): Mode shapes for v-velocity
|
| | 23 | % Pi(:,S): Mode shapes for Pressure
|
| | 24 | % rhoi(Optional): Mode shapes for density
|
| | 25 | %
|
| | 26 | % Written by: Matt Ueckermann
|
| | 27 |
|
| 20 | 28 | if ~isfield(app,'IP'),app.IP=[1,1,1];end
|
| 20 | 29 | if ~isfield(app,'Sb'), app.Sb = 0;end
|
| 20 | 30 | if ~isfield(app,'EvalTol'), app.EvalTol = 0;end
|
| | 31 |
|
| | 32 | %% Old orthonormalization procedure
|
| 20 | 33 | if isfield(app,'GSorth')
|
| | 34 | % display('Using old orthonormalization');
|
| | 35 | for i=1:app.S(1)
|
| | 36 | for j=1:app.S(1)
|
| | 37 | if i~=j
|
| | 38 | if nargin == 11
|
| | 39 | coef =app.IP(1) * InProd(ui(:,i), ui(:,j), dx, dy, uid)...
|
| | 40 | + app.IP(2) * InProd(vi(:,i), vi(:,j), dx ,dy, vid)...
|
| | 41 | + app.IP(3) * InProd(rhoi(:,i), rhoi(:,j), dx, dy, pid);
|
| | 42 | rhoi(:,i) = rhoi(:,i) - coef * rhoi(:,j);
|
| | 43 | else
|
| | 44 | coef = InProd(ui(:,i), ui(:,j), dx, dy, uid)...
|
| | 45 | + InProd(vi(:,i), vi(:,j), dx, dy, vid);
|
| | 46 | end
|
| | 47 | ui(:,i) = ui(:,i) - coef * ui(:,j);
|
| | 48 | vi(:,i) = vi(:,i) - coef * vi(:,j);
|
| | 49 | end
|
| | 50 | end
|
| | 51 | if nargin == 11
|
| | 52 | nrm = sqrt(app.IP(1)*InProd(ui(:,i), ui(:,i), dx, dy, uid)...
|
| | 53 | + app.IP(2)*InProd(vi(:,i), vi(:,i), dx ,dy, vid)...
|
| | 54 | + app.IP(3)*InProd(rhoi(:,i), rhoi(:,i), dx, dy, pid));
|
| | 55 | rhoi(:,i) = rhoi(:,i) / nrm;
|
| | 56 | else
|
| | 57 | nrm = sqrt(InProd(ui(:,i), ui(:,i), dx, dy, uid)...
|
| | 58 | + InProd(vi(:,i), vi(:,i), dx, dy, vid));
|
| | 59 | end
|
| | 60 | ui(:,i) = ui(:,i) / nrm;
|
| | 61 | vi(:,i) = vi(:,i) / nrm;
|
| | 62 | end
|
| 20 | 63 | else
|
| | 64 | %% First rotate the stochastic coefficients so that we have
|
| | 65 | %% uncorrelated samples
|
| 20 | 66 | Sintid = app.Sb+1:app.S;
|
0.03 | 20 | 67 | CYY=cov(YYt(:, Sintid));
|
< 0.01 | 20 | 68 | [VC, tmp] = eig(CYY) ; %Diagonal Covariance, and Diagonal Vectors
|
1.24 | 20 | 69 | ui(:, Sintid) = fliplr(ui(:, Sintid) * VC);
|
1.22 | 20 | 70 | vi(:, Sintid) = fliplr(vi(:, Sintid) * VC);
|
1.36 | 20 | 71 | Pi(:, Sintid) = fliplr(Pi(:, Sintid) * VC);
|
| 20 | 72 | if nargin == 11
|
2.30 | 20 | 73 | rhoi(:,Sintid) = fliplr(rhoi(:, Sintid) * VC);
|
| | 74 | else
|
| | 75 | rhoi = 0;
|
| | 76 | end
|
0.03 | 20 | 77 | YYt(:,Sintid) = fliplr(YYt(:,Sintid) * VC);
|
| | 78 |
|
| | 79 | %% Next orthonormalize the modes
|
| 20 | 80 | if nargin == 11
|
7.67 | 20 | 81 | M = (app.IP(1)* (ui(uid, :)' * ui(uid, :)) + app.IP(2)* (vi(vid, :)' * vi(vid, :)) ...
|
| | 82 | + app.IP(3) * (rhoi(pid, :)' * rhoi(pid, :)) ) * (dx * dy);
|
| | 83 | else
|
| | 84 | M = (ui(uid, :)' * ui(uid, :) + vi(vid, :)' * vi(vid, :)) * (dx * dy);
|
| | 85 | end
|
0.03 | 20 | 86 | [VC DC] = eig(M);
|
0.03 | 20 | 87 | YYt = fliplr(YYt * VC * sqrt(DC));
|
| | 88 | % YYt = fliplr(YYt * VC);
|
| 20 | 89 | if (app.EvalTol<=0)
|
| 20 | 90 | DC = diag (1 ./ diag(sqrt(DC)));
|
| | 91 | else
|
| | 92 | DC = sqrt(diag(DC));
|
| | 93 | zeroind = find (cumsum(DC) < (app.EvalTol.*sum(DC)) );
|
| | 94 | DC = 1./DC;
|
| | 95 | DC(zeroind) = 0;
|
| | 96 | DC = diag(DC);
|
| | 97 | end;
|
1.59 | 20 | 98 | ui = fliplr(ui * VC * DC);
|
1.48 | 20 | 99 | vi = fliplr(vi * VC * DC);
|
| 20 | 100 | if nargin == 11
|
3.63 | 20 | 101 | rhoi = fliplr(rhoi * VC * DC);
|
| 20 | 102 | end
|
1.52 | 20 | 103 | Pi = fliplr(Pi * VC * DC);
|
| 20 | 104 | Nbcs = length(ui) - length(uid);
|
| | 105 | %Now rotate the modes such that the boundary is orthogonal once again
|
| | 106 | %Now do the boundary-inner product to separate edge modes from non-boundary
|
| | 107 | %modes
|
| 20 | 108 | if nargin == 11
|
0.03 | 20 | 109 | [VC,DC] = eig(ui(1:Nbcs,:)'*ui(1:Nbcs,:) + vi(1:Nbcs,:)'*vi(1:Nbcs,:) + ...
|
| | 110 | rhoi(1:Nbcs,:)'*rhoi(1:Nbcs,:));
|
1.69 | 20 | 111 | rhoi = fliplr(rhoi*(VC));
|
| | 112 | else
|
| | 113 | [VC,DC] = eig(ui(1:Nbcs,:)'*ui(1:Nbcs,:) + vi(1:Nbcs,:)'*vi(1:Nbcs,:));
|
| | 114 | end
|
0.88 | 20 | 115 | ui = fliplr(ui*(VC));
|
0.86 | 20 | 116 | vi = fliplr(vi*(VC));
|
0.90 | 20 | 117 | Pi = fliplr(Pi*(VC));
|
0.04 | 20 | 118 | YYt = fliplr(YYt*(VC));
|
| | 119 |
|
| | 120 | %% Rotate the stochastic coefficients back to the uncorrelated case
|
0.03 | 20 | 121 | CYY=cov(YYt(:, Sintid));
|
0.01 | 20 | 122 | [VC,DC] = eig(CYY) ; %Diagonal Covariance, and Diagonal Vectors
|
1.23 | 20 | 123 | ui(:, Sintid) = fliplr(ui(:, Sintid) * VC);
|
1.14 | 20 | 124 | vi(:, Sintid) = fliplr(vi(:, Sintid) * VC);
|
1.33 | 20 | 125 | Pi(:, Sintid) = fliplr(Pi(:, Sintid) * VC);
|
| 20 | 126 | if nargin == 11
|
2.09 | 20 | 127 | rhoi(:, Sintid) = fliplr(rhoi(:, Sintid) * VC);
|
| 20 | 128 | end
|
| | 129 | %Also correct for, or make sure that, stochastic energy is preseved.
|
0.05 | 20 | 130 | YYt(:, Sintid) = fliplr(YYt(:, Sintid) * VC * ...
|
| | 131 | sqrt(sum(tmp(:)) / sum(DC(:))));
|
| | 132 | % YYt=fliplr(YYt * VC);
|
| 20 | 133 | end
|