| 1 |
function [ub,vb]=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 |
%rAw=dxc.*dyg; |
| 53 |
%rAs=dyc.*dxg; |
| 54 |
|
| 55 |
%--- recip_hFac & mask : |
| 56 |
mw=ceil(hw); mw=min(1,mw); |
| 57 |
ms=ceil(hs); ms=min(1,ms); |
| 58 |
|
| 59 |
%hw(find(hw==0))=Inf; |
| 60 |
%hs(find(hs==0))=Inf; |
| 61 |
%hw_recip=1./hw; %hw_recip(find(hw==0))=0; |
| 62 |
%hs_recip=1./hs; %hs_recip(find(hs==0))=0; |
| 63 |
|
| 64 |
switch GMform |
| 65 |
|
| 66 |
%%%%%%% Skew-flux form case |
| 67 |
case 'Skew' |
| 68 |
kwx_all = reshape(d.GM_Kwx,[6*nc*nc,nr,nt]); |
| 69 |
kwy_all = reshape(d.GM_Kwy,[6*nc*nc,nr,nt]); |
| 70 |
|
| 71 |
kwx_all = 0.5*kwx_all; |
| 72 |
kwy_all = 0.5*kwy_all; |
| 73 |
|
| 74 |
for it = 1:nt |
| 75 |
kwx = kwx_all(:,:,it); |
| 76 |
kwy = kwy_all(:,:,it); |
| 77 |
|
| 78 |
%-- K*ra + add 1 overlap : |
| 79 |
kwx = (ra*ones(1,nr)).*kwx; |
| 80 |
kwy = (ra*ones(1,nr)).*kwy; |
| 81 |
kwx = reshape(kwx,[6*nc,nc,nr]); |
| 82 |
kwy = reshape(kwy,[6*nc,nc,nr]); |
| 83 |
v6X = split_C_cub(kwx,1); |
| 84 |
v6Y = split_C_cub(kwy,1); |
| 85 |
k6x = v6X(:,[2:nc+1],:,:); |
| 86 |
k6y = v6Y([2:nc+1],:,:,:); |
| 87 |
|
| 88 |
%----------------- |
| 89 |
v6X = zeros(nc,nc,nr,6); |
| 90 |
v6Y = zeros(nc,nc,nr,6); |
| 91 |
|
| 92 |
v6X([1:nc],:,:,:) = k6x([2:nc+1],:,:,:) + k6x([1:nc],:,:,:); |
| 93 |
v6Y(:,[1:nc],:,:) = k6y(:,[2:nc+1],:,:) + k6y(:,[1:nc],:,:); |
| 94 |
|
| 95 |
v6X = v6X/2; |
| 96 |
v6Y = v6Y/2; |
| 97 |
|
| 98 |
psiX = zeros(6*nc,nc,nr+1); |
| 99 |
psiY = zeros(6*nc,nc,nr+1); |
| 100 |
|
| 101 |
for n = 1:6 |
| 102 |
is = 1+nc*(n-1);ie=nc*n; |
| 103 |
psiX([is:ie],[1:nc],[1:nr]) = v6X([1:nc],[1:nc],[1:nr],n); |
| 104 |
psiY([is:ie],[1:nc],[1:nr]) = v6Y([1:nc],[1:nc],[1:nr],n); |
| 105 |
end |
| 106 |
|
| 107 |
psiX = reshape(psiX,6*nc*nc,nr+1); |
| 108 |
psiY = reshape(psiY,6*nc*nc,nr+1); |
| 109 |
|
| 110 |
psiX(:,[1:nr]) = mw.*psiX(:,[1:nr]); |
| 111 |
psiY(:,[1:nr]) = ms.*psiY(:,[1:nr]); |
| 112 |
ub = psiX(:,[2:nr+1]) - psiX(:,[1:nr]); |
| 113 |
vb = psiY(:,[2:nr+1]) - psiY(:,[1:nr]); |
| 114 |
|
| 115 |
dr = reshape(dr,[1,length(dr)]); |
| 116 |
% ub = reshape(hw_recip.*ub./(rAw*dr),[6*nc,nc,nr]); |
| 117 |
ub = reshape(ub./(rAw*dr),[6*nc,nc,nr]); |
| 118 |
% vb = reshape(hs_recip.*vb./(rAs*dr),[6*nc,nc,nr]); |
| 119 |
vb = reshape(vb./(rAs*dr),[6*nc,nc,nr]); |
| 120 |
|
| 121 |
ub_all(:,:,:,it) = ub; |
| 122 |
vb_all(:,:,:,it) = vb; |
| 123 |
end |
| 124 |
|
| 125 |
%%%%%%% Advective form case |
| 126 |
case 'Advc' |
| 127 |
|
| 128 |
PsiX_all = reshape(d.GM_PsiX(1:6*nc,1:nc,1:nr,:),[6*nc*nc,nr,nt]); |
| 129 |
PsiY_all = reshape(d.GM_PsiY(1:6*nc,1:nc,1:nr,:),[6*nc*nc,nr,nt]); |
| 130 |
|
| 131 |
dr3d = ones(6*nc*nc,1)*reshape(dr,[1,length(dr)]); |
| 132 |
|
| 133 |
for it = 1:nt |
| 134 |
|
| 135 |
psiX = zeros(6*nc*nc,nr+1); |
| 136 |
psiY = zeros(6*nc*nc,nr+1); |
| 137 |
|
| 138 |
psiX(:,1:nr) = mw.*PsiX_all(:,:,it); |
| 139 |
psiY(:,1:nr) = ms.*PsiY_all(:,:,it); |
| 140 |
|
| 141 |
% psiX(:,[1:nr]) = mw.*psiX(:,[1:nr]); |
| 142 |
% psiY(:,[1:nr]) = ms.*psiY(:,[1:nr]); |
| 143 |
ub = psiX(:,[2:nr+1]) - psiX(:,[1:nr]); |
| 144 |
vb = psiY(:,[2:nr+1]) - psiY(:,[1:nr]); |
| 145 |
|
| 146 |
% ub = reshape(hw_recip.*ub./(rAw*dr),[6*nc,nc,nr]); |
| 147 |
% ub = reshape(ub./(rAw*dr),[6*nc,nc,nr]); |
| 148 |
ub = reshape(ub./dr3d,[6*nc,nc,nr]); |
| 149 |
% vb = reshape(hs_recip.*vb./(rAs*dr),[6*nc,nc,nr]); |
| 150 |
% vb = reshape(vb./(rAs*dr),[6*nc,nc,nr]); |
| 151 |
vb = reshape(vb./dr3d,[6*nc,nc,nr]); |
| 152 |
|
| 153 |
ub_all(:,:,:,it) = ub; |
| 154 |
vb_all(:,:,:,it) = vb; |
| 155 |
end |
| 156 |
|
| 157 |
otherwise |
| 158 |
disp('Ce portnawak') |
| 159 |
end |
| 160 |
|
| 161 |
ub = ub_all; |
| 162 |
vb = vb_all; |
Parent Directory
| 
Revision Log
| 
Revision Graph