| 1 |
function [ub,vb,wb]=calcBolusVelCube(d,g,GMform); |
| 2 |
|
| 3 |
% [ub,vb] = calcBolusVelCube(d,g,GMform); |
| 4 |
% |
| 5 |
% Input arguments: |
| 6 |
% The incoming field data (d) and grid data (g) must be in a structured |
| 7 |
% array format (which is the format that comes from rdmnc): |
| 8 |
% d [Field data] Kwx,Kwy |
| 9 |
% g [Grid data ] drF,rA,dxC,dyC,dxG,dyG,HFacW,HFacS |
| 10 |
% GMform [string] GM form, 'Skew' or 'Advc' |
| 11 |
% |
| 12 |
% Output arguments: |
| 13 |
% ub, vb: GM-Bolus mass-weigthed velocity (i.e include |
| 14 |
% implicitly hFac factor) |
| 15 |
% |
| 16 |
% Comments: |
| 17 |
% For Skew-Flux form: uses Kwx & Kwy divided by 2 |
| 18 |
% compute Volume Stream function psiX,psiY above uVel.vVel |
| 19 |
% (at interface between 2 levels), units=m^3/s : |
| 20 |
% psiX=(rAc*kwx)_i / dXc ; psiY=(rAc*kwy)_j / dYc ; |
| 21 |
% and then the bolus velocity (m/s): |
| 22 |
% ub = d_k(psiX)/rAw/drF ; vb = d_k(psiY)/rAs/drF ; |
| 23 |
% |
| 24 |
%--------------------------------------------------------------------- |
| 25 |
|
| 26 |
|
| 27 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 28 |
% Prepare / reform incoming data % |
| 29 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 30 |
|
| 31 |
nc = size(g.XC,2); |
| 32 |
nr = length(g.drF); |
| 33 |
|
| 34 |
switch GMform |
| 35 |
case 'Skew' |
| 36 |
nt = size(d.GM_Kwx,4); |
| 37 |
case 'Advc' |
| 38 |
nt = size(d.GM_PsiX,4); |
| 39 |
end |
| 40 |
|
| 41 |
dr = g.drF; |
| 42 |
hw = reshape(g.HFacW(1:6*nc,1:nc,1:nr),[6*nc*nc,nr]); |
| 43 |
hs = reshape(g.HFacS(1:6*nc,1:nc,1:nr),[6*nc*nc,nr]); |
| 44 |
ra = reshape(g.rA(1:6*nc,1:nc) ,[6*nc*nc,1]); |
| 45 |
rAs = reshape(g.rAs(1:6*nc,1:nc) ,[6*nc*nc,1]); |
| 46 |
rAw = reshape(g.rAw(1:6*nc,1:nc) ,[6*nc*nc,1]); |
| 47 |
dxc = reshape(g.dxC(1:6*nc,1:nc),[6*nc*nc,1]); |
| 48 |
dyc = reshape(g.dyC(1:6*nc,1:nc),[6*nc*nc,1]); |
| 49 |
dxg = reshape(g.dxG(1:6*nc,1:nc),[6*nc*nc,1]); |
| 50 |
dyg = reshape(g.dyG(1:6*nc,1:nc),[6*nc*nc,1]); |
| 51 |
|
| 52 |
%--- recip_hFac & mask : |
| 53 |
mw=ceil(hw); mw=min(1,mw); |
| 54 |
ms=ceil(hs); ms=min(1,ms); |
| 55 |
|
| 56 |
hw(find(hw==0))=Inf; |
| 57 |
hs(find(hs==0))=Inf; |
| 58 |
hw_recip=1./hw; %hw_recip(find(hw==0))=0; |
| 59 |
hs_recip=1./hs; %hs_recip(find(hs==0))=0; |
| 60 |
|
| 61 |
ub_all = zeros(6*nc,nc,nr,nt); |
| 62 |
vb_all = zeros(6*nc,nc,nr,nt); |
| 63 |
wb_all = zeros(6*nc,nc,nr,nt); |
| 64 |
|
| 65 |
switch GMform |
| 66 |
|
| 67 |
%%%%%%% Skew-flux form case |
| 68 |
case 'Skew' |
| 69 |
kwx_all = reshape(d.GM_Kwx,[6*nc*nc,nr,nt]); |
| 70 |
kwy_all = reshape(d.GM_Kwy,[6*nc*nc,nr,nt]); |
| 71 |
|
| 72 |
kwx_all = 0.5*kwx_all; |
| 73 |
kwy_all = 0.5*kwy_all; |
| 74 |
|
| 75 |
for it = 1:nt |
| 76 |
kwx = kwx_all(:,:,it); |
| 77 |
kwy = kwy_all(:,:,it); |
| 78 |
|
| 79 |
%-- K*ra + add 1 overlap : |
| 80 |
kwx = (ra*ones(1,nr)).*kwx; |
| 81 |
kwy = (ra*ones(1,nr)).*kwy; |
| 82 |
kwx = reshape(kwx,[6*nc,nc,nr]); |
| 83 |
kwy = reshape(kwy,[6*nc,nc,nr]); |
| 84 |
v6X = split_C_cub(kwx,1); |
| 85 |
v6Y = split_C_cub(kwy,1); |
| 86 |
k6x = v6X(:,[2:nc+1],:,:); |
| 87 |
k6y = v6Y([2:nc+1],:,:,:); |
| 88 |
|
| 89 |
%----------------- |
| 90 |
v6X = zeros(nc,nc,nr,6); |
| 91 |
v6Y = zeros(nc,nc,nr,6); |
| 92 |
|
| 93 |
v6X([1:nc],:,:,:) = k6x([2:nc+1],:,:,:) + k6x([1:nc],:,:,:); |
| 94 |
v6Y(:,[1:nc],:,:) = k6y(:,[2:nc+1],:,:) + k6y(:,[1:nc],:,:); |
| 95 |
|
| 96 |
v6X = v6X/2; |
| 97 |
v6Y = v6Y/2; |
| 98 |
|
| 99 |
psiX = zeros(6*nc,nc,nr+1); |
| 100 |
psiY = zeros(6*nc,nc,nr+1); |
| 101 |
|
| 102 |
for n = 1:6 |
| 103 |
is = 1+nc*(n-1);ie=nc*n; |
| 104 |
psiX([is:ie],[1:nc],[1:nr]) = v6X([1:nc],[1:nc],[1:nr],n); |
| 105 |
psiY([is:ie],[1:nc],[1:nr]) = v6Y([1:nc],[1:nc],[1:nr],n); |
| 106 |
end |
| 107 |
|
| 108 |
psiX = reshape(psiX,6*nc*nc,nr+1); |
| 109 |
psiY = reshape(psiY,6*nc*nc,nr+1); |
| 110 |
|
| 111 |
psiX(:,[1:nr]) = mw.*psiX(:,[1:nr]); |
| 112 |
psiY(:,[1:nr]) = ms.*psiY(:,[1:nr]); |
| 113 |
ub = psiX(:,[2:nr+1]) - psiX(:,[1:nr]); |
| 114 |
vb = psiY(:,[2:nr+1]) - psiY(:,[1:nr]); |
| 115 |
|
| 116 |
dr = reshape(dr,[1,length(dr)]); |
| 117 |
% ub = reshape(hw_recip.*ub./(rAw*dr),[6*nc,nc,nr]); |
| 118 |
ub = reshape(ub./(rAw*dr),[6*nc,nc,nr]); |
| 119 |
% vb = reshape(hs_recip.*vb./(rAs*dr),[6*nc,nc,nr]); |
| 120 |
vb = reshape(vb./(rAs*dr),[6*nc,nc,nr]); |
| 121 |
|
| 122 |
ub_all(:,:,:,it) = ub; |
| 123 |
vb_all(:,:,:,it) = vb; |
| 124 |
|
| 125 |
%%%%%%%%%%%%% |
| 126 |
[u6t,v6t] =split_UV_cub(ub,vb,0,1); |
| 127 |
[hw6t,hs6t]=split_UV_cub(reshape(hw,6*nc,nc,nr),reshape(hs,6*nc,nc,nr),0,1); |
| 128 |
[dy6t,dx6t]=split_UV_cub(reshape(dyg,6*nc,nc),reshape(dxg,6*nc,nc),0,1); |
| 129 |
|
| 130 |
%F6tX = u6t.*hw6t.*permute(repmat(dy6t,[1 1 1 nr]),[1 2 4 3]); |
| 131 |
%F6tY = v6t.*hs6t.*permute(repmat(dx6t,[1 1 1 nr]),[1 2 4 3]); |
| 132 |
F6tX = u6t.*permute(repmat(dy6t,[1 1 1 nr]),[1 2 4 3]); |
| 133 |
F6tY = v6t.*permute(repmat(dx6t,[1 1 1 nr]),[1 2 4 3]); |
| 134 |
|
| 135 |
Hdiv = zeros(nc,nc,nr,6); |
| 136 |
Hdiv = F6tX([2:nc+1],:,:,:) - F6tX([1:nc],:,:,:) ... |
| 137 |
+ F6tY(:,[2:nc+1],:,:) - F6tY(:,[1:nc],:,:); |
| 138 |
for k=1:nr |
| 139 |
Hdiv(:,:,k,:) = -Hdiv(:,:,k,:)*dr(k); |
| 140 |
end |
| 141 |
|
| 142 |
%psiX = zeros(6*nc,nc,nr); |
| 143 |
%for n = 1:6 |
| 144 |
% is = 1+nc*(n-1);ie=nc*n; |
| 145 |
% psiX([is:ie],[1:nc],[1:nr]) = Hdiv([1:nc],[1:nc],[1:nr],n); |
| 146 |
%end |
| 147 |
psiX = reshape(permute(Hdiv,[1 4 2 3]),6*nc,nc,nr); |
| 148 |
%wb = psiX ./ repmat(reshape(ra,6*nc,nc),[1 1 nr]); |
| 149 |
|
| 150 |
wb = zeros(6*nc,nc,nr+1); |
| 151 |
for k=nr:-1:1 |
| 152 |
wb(:,:,k)= psiX(:,:,k) ./ reshape(ra,6*nc,nc) + wb(:,:,k+1); |
| 153 |
end |
| 154 |
|
| 155 |
wb_all(:,:,:,it) = wb(:,:,1:30); |
| 156 |
%%%%%%%%%%%% |
| 157 |
end |
| 158 |
|
| 159 |
%%%%%%% Advective form case |
| 160 |
case 'Advc' |
| 161 |
|
| 162 |
PsiX_all = reshape(d.GM_PsiX(1:6*nc,1:nc,1:nr,:),[6*nc*nc,nr,nt]); |
| 163 |
PsiY_all = reshape(d.GM_PsiY(1:6*nc,1:nc,1:nr,:),[6*nc*nc,nr,nt]); |
| 164 |
|
| 165 |
dr3d = ones(6*nc*nc,1)*reshape(dr,[1,length(dr)]); |
| 166 |
|
| 167 |
for it = 1:nt |
| 168 |
|
| 169 |
psiX = zeros(6*nc*nc,nr+1); |
| 170 |
psiY = zeros(6*nc*nc,nr+1); |
| 171 |
|
| 172 |
psiX(:,1:nr) = mw.*PsiX_all(:,:,it); |
| 173 |
psiY(:,1:nr) = ms.*PsiY_all(:,:,it); |
| 174 |
|
| 175 |
% psiX(:,[1:nr]) = mw.*psiX(:,[1:nr]); |
| 176 |
% psiY(:,[1:nr]) = ms.*psiY(:,[1:nr]); |
| 177 |
ub = psiX(:,[2:nr+1]) - psiX(:,[1:nr]); |
| 178 |
vb = psiY(:,[2:nr+1]) - psiY(:,[1:nr]); |
| 179 |
|
| 180 |
% ub = reshape(hw_recip.*ub./(rAw*dr),[6*nc,nc,nr]); |
| 181 |
% ub = reshape(ub./(rAw*dr),[6*nc,nc,nr]); |
| 182 |
ub = reshape(ub./dr3d,[6*nc,nc,nr]); |
| 183 |
% vb = reshape(hs_recip.*vb./(rAs*dr),[6*nc,nc,nr]); |
| 184 |
% vb = reshape(vb./(rAs*dr),[6*nc,nc,nr]); |
| 185 |
vb = reshape(vb./dr3d,[6*nc,nc,nr]); |
| 186 |
|
| 187 |
ub_all(:,:,:,it) = ub; |
| 188 |
vb_all(:,:,:,it) = vb; |
| 189 |
end |
| 190 |
|
| 191 |
otherwise |
| 192 |
disp('Ce portnawak') |
| 193 |
end |
| 194 |
|
| 195 |
ub = ub_all; |
| 196 |
vb = vb_all; |
| 197 |
wb = wb_all; |
Parent Directory
| 
Revision Log
| 
Revision Graph