| 1 |
function [fldUdiv,fldVdiv,fldDivPot]=diffsmooth2D_div_inv(fldU,fldV); |
| 2 |
|
| 3 |
%object: compute the divergent part of a tranport vector field |
| 4 |
% by solving a Laplace equation with Neumann B.C. |
| 5 |
% |
| 6 |
%input: fldU,fldV transport vectors (masked with NaNs) |
| 7 |
%output:fldUdiv,fldVdiv divergent part |
| 8 |
% fldDivPot Potential for the divergent flow (optional) |
| 9 |
% |
| 10 |
%note 1: here we do everything in terms of transports, so we omit all |
| 11 |
% grid factors (i.e. with dx=dy=1), and the resulting potential |
| 12 |
% has usits of transports. |
| 13 |
%note 2: this routine was derived from diffsmooth2D_extrap_inv.m but |
| 14 |
% it differs in several aspects. It : 1) omits grid factors; |
| 15 |
% 2) uses of a neumann, rather than dirichlet, boundary condition; |
| 16 |
% 3) gets the mask directly from fld. |
| 17 |
|
| 18 |
global mygrid; |
| 19 |
|
| 20 |
warning('off','MATLAB:divideByZero'); |
| 21 |
warning('off','MATLAB:nearlySingularMatrix'); |
| 22 |
|
| 23 |
fld=calc_UV_conv(fldU,fldV,{});%no grid scaling factor in div. or transports |
| 24 |
fld=fld.*mygrid.mskC(:,:,1); |
| 25 |
%this can create isolated Canyons and a singular matrix: fld(fld==0)=NaN; |
| 26 |
|
| 27 |
mskWet=fld; mskWet(~isnan(fld))=1; |
| 28 |
mskDry=1*isnan(mskWet); mskDry(mskDry==0)=NaN; |
| 29 |
|
| 30 |
%check for domain edge points where no exchange is possible: |
| 31 |
tmp1=mskWet; tmp1(:)=1; tmp2=exch_T_N(tmp1); |
| 32 |
for iF=1:mskWet.nFaces; |
| 33 |
tmp3=mskWet{iF}; tmp4=tmp2{iF}; |
| 34 |
tmp4=tmp4(2:end-1,1:end-2)+tmp4(2:end-1,3:end)+tmp4(1:end-2,2:end-1)+tmp4(3:end,2:end-1); |
| 35 |
if ~isempty(find(isnan(tmp4)&~isnan(tmp3))); fprintf('warning: mask was modified\n'); end; |
| 36 |
tmp3(isnan(tmp4))=NaN; mskWet{iF}=tmp3; |
| 37 |
end; |
| 38 |
%why does this matters... see set-ups with open boundary conditions... |
| 39 |
|
| 40 |
%put 0 first guess if needed and switch land mask: |
| 41 |
fld(find(isnan(fld)))=0; fld=fld.*mskWet; |
| 42 |
|
| 43 |
%define mapping from global array to (no nan points) global vector |
| 44 |
tmp1=convert2array(mskWet); |
| 45 |
kk=find(~isnan(tmp1)); nn=length(kk); |
| 46 |
KK=tmp1; KK(:)=0; KK(kk)=kk; KKfaces=convert2array(KK); %global array indices |
| 47 |
LL=tmp1; LL(:)=0; LL(kk)=[1:nn]; LLfaces=convert2array(LL); %global vector indices |
| 48 |
|
| 49 |
%note: diffsmooth2D_extrap_inv.m uses a Dirichlet boundary condition |
| 50 |
% so that it uses a nan/1 mask in KK/LL // here we use a |
| 51 |
% Neumann boundary condition so we use a 0/1 mask (see below also) |
| 52 |
|
| 53 |
%form matrix problem: |
| 54 |
A=sparse([],[],[],nn,nn,nn*5); |
| 55 |
for iF=1:fld.nFaces; for ii=1:3; for jj=1:3; |
| 56 |
|
| 57 |
%1) seed points (FLDones) and neighborhood of influence (FLDkkFROM) |
| 58 |
FLDones=fld; FLDones(:)=0; |
| 59 |
FLDones{iF}(ii:3:end,jj:3:end)=1; |
| 60 |
FLDones(KKfaces==0)=0; |
| 61 |
|
| 62 |
FLDkkFROMtmp=fld; FLDkkFROMtmp(:)=0; |
| 63 |
FLDkkFROMtmp{iF}(ii:3:end,jj:3:end)=KKfaces{iF}(ii:3:end,jj:3:end); |
| 64 |
FLDkkFROMtmp(find(isnan(fld)))=0; |
| 65 |
|
| 66 |
FLDkkFROM=exch_T_N(FLDkkFROMtmp); FLDkkFROM(isnan(FLDkkFROM))=0; |
| 67 |
for iF2=1:fld.nFaces; |
| 68 |
tmp1=FLDkkFROM{iF2}; tmp2=zeros(size(tmp1)-2); |
| 69 |
for ii2=1:3; for jj2=1:3; tmp2=tmp2+tmp1(ii2:end-3+ii2,jj2:end-3+jj2); end; end; |
| 70 |
FLDkkFROM{iF2}=tmp2; |
| 71 |
end; |
| 72 |
|
| 73 |
%2) compute effect of each point on neighboring target point: |
| 74 |
[tmpU,tmpV]=calc_T_grad(FLDones,0); |
| 75 |
%unlike calc_T_grad, we work in grid point index, so we need to omit grid factors |
| 76 |
tmpU=tmpU.*mygrid.DXC; tmpV=tmpV.*mygrid.DYC; |
| 77 |
%and accordingly we use no grid scaling factor in div. |
| 78 |
[dFLDdt]=calc_UV_conv(tmpU,tmpV,{}); |
| 79 |
|
| 80 |
%note: diffsmooth2D_extrap_inv.m uses a Dirichlet boundary condition |
| 81 |
% so that it needs to apply mskFreeze to dFLDdt // here we use a |
| 82 |
% Neumann boundary condition so we do not mask dFLDdt (see above also) |
| 83 |
|
| 84 |
%3) include seed contributions in matrix: |
| 85 |
FLDkkFROM=convert2array(FLDkkFROM); |
| 86 |
%3.1) for wet points |
| 87 |
dFLDdtWet=convert2array(dFLDdt.*mskWet); |
| 88 |
tmp1=find(dFLDdtWet~=0&~isnan(dFLDdtWet)); |
| 89 |
dFLDdtWet=dFLDdtWet(tmp1); FLDkkFROMtmp=FLDkkFROM(tmp1); FLDkkTOtmp=KK(tmp1); |
| 90 |
A=A+sparse(LL(FLDkkTOtmp),LL(FLDkkFROMtmp),dFLDdtWet,nn,nn); |
| 91 |
%3.2) for dry points (-> this part reflects the neumann boundary condition) |
| 92 |
dFLDdtDry=convert2array(dFLDdt.*mskDry); |
| 93 |
tmp1=find(dFLDdtDry~=0&~isnan(dFLDdtDry)); |
| 94 |
dFLDdtDry=dFLDdtDry(tmp1); FLDkkFROMtmp=FLDkkFROM(tmp1); |
| 95 |
A=A+sparse(LL(FLDkkFROMtmp),LL(FLDkkFROMtmp),dFLDdtDry,nn,nn); |
| 96 |
|
| 97 |
end; end; end; |
| 98 |
|
| 99 |
%to check results: |
| 100 |
%figure; spy(A); |
| 101 |
|
| 102 |
%4) solve for potential: |
| 103 |
yy=convert2array(fld); yy=yy(find(KK~=0)); |
| 104 |
xx=A\yy; |
| 105 |
yyFROMxx=A*xx; |
| 106 |
|
| 107 |
%5) prepare output: |
| 108 |
fldDivPot=0*convert2array(fld); fldDivPot(find(KK~=0))=xx; |
| 109 |
fldDivPot=convert2array(fldDivPot); |
| 110 |
[fldUdiv,fldVdiv]=calc_T_grad(fldDivPot,0); |
| 111 |
%unlike calc_T_grad, we work in grid point index, so we need to omit grid factors |
| 112 |
fldUdiv=fldUdiv.*mygrid.DXC; fldVdiv=fldVdiv.*mygrid.DYC; |
| 113 |
|
| 114 |
%to check the results: |
| 115 |
%fld1=0*convert2array(fld); fld1(find(KK~=0))=xx; |
| 116 |
%fld2=0*convert2array(fld); fld2(find(KK~=0))=yyFROMxx; |
| 117 |
%fld3=convert2array(calc_UV_conv(fldU,fldV,{}));%no grid scaling factor in div. or transports |
| 118 |
|
| 119 |
warning('on','MATLAB:divideByZero'); |
| 120 |
warning('on','MATLAB:nearlySingularMatrix'); |
| 121 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph