| 1 |
enderton |
1.1 |
function [ub,vb]=calc_Bolus_CS(rAc,dXc,dYc,dXg,dYg,delR,kwx,kwy,hFacW,hFacS); |
| 2 |
|
|
% [ub,vb]=calc_Bolus_CS(rAc,dXc,dYc,dXg,dYg,delR,kwx,kwy,hFacW,hFacS); |
| 3 |
|
|
%-- GM-Bolus transp. : |
| 4 |
|
|
% from Kwx & Kwy (Skew-Flx => /2), |
| 5 |
|
|
% compute Volume Stream function psiX,psiY above uVel.vVel |
| 6 |
|
|
% (at interface between 2 levels), units=m^3/s : |
| 7 |
|
|
% psiX=(rAc*kwx)_i / dXc ; psiY=(rAc*kwy)_j / dYc ; |
| 8 |
|
|
% and then the bolus velocity (m/s): |
| 9 |
|
|
% ub = d_k(psiX)/dYg/drF ; vb = d_k(psiY)/dXg/drF ; |
| 10 |
|
|
%--------------------------------------------------------------------- |
| 11 |
|
|
% set_axis |
| 12 |
|
|
[ncx,nc]=size(rAc); nr=length(delR); |
| 13 |
|
|
nPg=6*nc*nc; ncp=nc+1; |
| 14 |
|
|
|
| 15 |
|
|
rAu=dXc.*dYg; rAu=reshape(rAu,nPg,1); |
| 16 |
|
|
rAv=dYc.*dXg; rAv=reshape(rAv,nPg,1); |
| 17 |
|
|
|
| 18 |
|
|
%-- K*rAc + add 1 overlap : |
| 19 |
|
|
tmp=reshape(rAc,nPg,1)*ones(1,nr); |
| 20 |
|
|
kwx=reshape(kwx,nPg,nr); |
| 21 |
|
|
kwy=reshape(kwy,nPg,nr); |
| 22 |
|
|
tmp=0.5*tmp; |
| 23 |
|
|
kwx=tmp.*kwx; |
| 24 |
|
|
kwy=tmp.*kwy; |
| 25 |
|
|
kwx=reshape(kwx,ncx,nc,nr); |
| 26 |
|
|
kwy=reshape(kwy,ncx,nc,nr); |
| 27 |
|
|
v6X=split_C_cub(kwx,1); |
| 28 |
|
|
v6Y=split_C_cub(kwy,1); |
| 29 |
|
|
k6x=v6X(:,[2:ncp],:,:); |
| 30 |
|
|
k6y=v6Y([2:ncp],:,:,:); |
| 31 |
|
|
|
| 32 |
|
|
%--- recip_hFac & mask : |
| 33 |
|
|
hFacW=reshape(hFacW,nPg,nr); |
| 34 |
|
|
hFacS=reshape(hFacS,nPg,nr); |
| 35 |
|
|
rhFcW=1./hFacW; rhFcW(find(hFacW==0))=0; |
| 36 |
|
|
rhFcS=1./hFacS; rhFcS(find(hFacS==0))=0; |
| 37 |
|
|
hFacW=ceil(hFacW); hFacW=min(1,hFacW); |
| 38 |
|
|
hFacS=ceil(hFacS); hFacS=min(1,hFacS); |
| 39 |
|
|
|
| 40 |
|
|
%----------------- |
| 41 |
|
|
|
| 42 |
|
|
v6X=zeros(nc,nc,nr,6); |
| 43 |
|
|
v6X([1:nc],:,:,:)=k6x([2:ncp],:,:,:)+k6x([1:nc],:,:,:); |
| 44 |
|
|
v6X=v6X/2; |
| 45 |
|
|
psiX=zeros(ncx,nc,nr+1); |
| 46 |
|
|
for n=1:6 |
| 47 |
|
|
is=1+nc*(n-1);ie=nc*n; |
| 48 |
|
|
psiX([is:ie],[1:nc],[1:nr])=v6X([1:nc],[1:nc],[1:nr],n); |
| 49 |
|
|
end |
| 50 |
|
|
psiX=reshape(psiX,nPg,nr+1); |
| 51 |
|
|
psiX(:,[1:nr])=hFacW.*psiX(:,[1:nr]); |
| 52 |
|
|
ub=zeros(nPg,nr); |
| 53 |
|
|
ub=psiX(:,[2:nr+1])-psiX(:,[1:nr]); |
| 54 |
|
|
delR = reshape(delR,[1,length(delR)]); |
| 55 |
|
|
tmp=rAu*delR; |
| 56 |
|
|
ub=ub./tmp; ub=rhFcW.*ub; |
| 57 |
|
|
ub=reshape(ub,ncx,nc,nr); |
| 58 |
|
|
|
| 59 |
|
|
%----------------- |
| 60 |
|
|
|
| 61 |
|
|
v6Y=zeros(nc,nc,nr,6); |
| 62 |
|
|
v6Y(:,[1:nc],:,:)=k6y(:,[2:ncp],:,:)+k6y(:,[1:nc],:,:); |
| 63 |
|
|
v6Y=v6Y/2; |
| 64 |
|
|
psiY=zeros(ncx,nc,nr+1); |
| 65 |
|
|
for n=1:6 |
| 66 |
|
|
is=1+nc*(n-1);ie=nc*n; |
| 67 |
|
|
psiY([is:ie],[1:nc],[1:nr])=v6Y([1:nc],[1:nc],[1:nr],n); |
| 68 |
|
|
end |
| 69 |
|
|
psiY=reshape(psiY,nPg,nr+1); |
| 70 |
|
|
psiY(:,[1:nr])=hFacS.*psiY(:,[1:nr]); |
| 71 |
|
|
vb=zeros(nPg,nr); |
| 72 |
|
|
vb=psiY(:,[2:nr+1])-psiY(:,[1:nr]); |
| 73 |
|
|
tmp=rAv*delR; |
| 74 |
|
|
vb=vb./tmp; vb=rhFcS.*vb; |
| 75 |
|
|
vb=reshape(vb,ncx,nc,nr); |
| 76 |
|
|
|
| 77 |
|
|
%----------------- |
| 78 |
|
|
return |
Parent Directory
| 
Revision Log
| 
Revision Graph