matlab_class/gcmfaces_smooth/diffsmooth2D_div_inv.m

function [fldUdiv,fldVdiv,fldDivPot]=diffsmooth2D_div_inv(fldU,fldV);

%object: compute the divergent part of a tranport vector field
%       by solving a Laplace equation with Neumann B.C.
%
%input: fldU,fldV       transport vectors (masked with NaNs)
%output:fldUdiv,fldVdiv divergent part 
%       fldDivPot       Potential for the divergent flow (optional)     
%
%note 1: here we do everything in terms of transports, so we omit all
%       grid factors (i.e. with dx=dy=1), and the resulting potential 
%       has usits of transports.
%note 2: this routine was derived from diffsmooth2D_extrap_inv.m but
%       it differs in several aspects. It : 1) omits grid factors;
%       2) uses of a neumann, rather than dirichlet, boundary condition;
%       3) gets the mask directly from fld.

global mygrid;

warning('off','MATLAB:divideByZero');

fld=calc_UV_div(fldU,fldV,{});%no grid scaling factor in div. or transports
fld=fld.*mygrid.mskC(:,:,1);
fld(fld==0)=NaN;
 
mskWet=fld; mskWet(~isnan(fld))=1;
mskDry=1*isnan(mskWet); mskDry(mskDry==0)=NaN;

%check for domain edge points where no exchange is possible:
tmp1=mskWet; tmp1(:)=1; tmp2=exch_T_N(tmp1);
for iF=1:mskWet.nFaces;
tmp3=mskWet{iF}; tmp4=tmp2{iF}; 
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);
if ~isempty(find(isnan(tmp4)&~isnan(tmp3))); fprintf('warning: mask was modified\n'); end;
tmp3(isnan(tmp4))=NaN; mskWet{iF}=tmp3; 
end;
%why does this matters... see set-ups with open boundary conditions...

%put 0 first guess if needed and switch land mask:
fld(find(isnan(fld)))=0; fld=fld.*mskWet;

%define mapping from global array to (no nan points) global vector
tmp1=convert2array(mskWet); 
kk=find(~isnan(tmp1)); nn=length(kk);
KK=tmp1; KK(:)=0; KK(kk)=kk; KKfaces=convert2array(KK); %global array indices
LL=tmp1; LL(:)=0; LL(kk)=[1:nn]; LLfaces=convert2array(LL); %global vector indices

%note:  diffsmooth2D_extrap_inv.m uses a Dirichlet boundary condition
%       so that it uses a nan/1 mask in KK/LL // here we use a 
%       Neumann boundary condition so we use a 0/1 mask (see below also)

%form matrix problem:
A=sparse([],[],[],nn,nn,nn*5);
for iF=1:fld.nFaces; for ii=1:3; for jj=1:3;

%1) seed points (FLDones) and neighborhood of influence (FLDkkFROM)
    FLDones=fld; FLDones(:)=0; 
    FLDones{iF}(ii:3:end,jj:3:end)=1; 
    FLDones(KKfaces==0)=0;

    FLDkkFROMtmp=fld; FLDkkFROMtmp(:)=0;
    FLDkkFROMtmp{iF}(ii:3:end,jj:3:end)=KKfaces{iF}(ii:3:end,jj:3:end);
    FLDkkFROMtmp(find(isnan(fld)))=0;

    FLDkkFROM=exch_T_N(FLDkkFROMtmp); FLDkkFROM(isnan(FLDkkFROM))=0; 
    for iF2=1:fld.nFaces; 
       tmp1=FLDkkFROM{iF2}; tmp2=zeros(size(tmp1)-2);
       for ii2=1:3; for jj2=1:3; tmp2=tmp2+tmp1(ii2:end-3+ii2,jj2:end-3+jj2); end; end;
       FLDkkFROM{iF2}=tmp2; 
    end;

%2) compute effect of each point on neighboring target point:
    [tmpU,tmpV]=calc_T_grad(FLDones,0);
%unlike calc_T_grad, we work in grid point index, so we need to omit grid factors
    tmpU=tmpU.*mygrid.DXC; tmpV=tmpV.*mygrid.DYC;
%and accordingly we use no grid scaling factor in div.
    [dFLDdt]=calc_UV_div(tmpU,tmpV,{});

%note:  diffsmooth2D_extrap_inv.m uses a Dirichlet boundary condition
%       so that it needs to apply mskFreeze to dFLDdt // here we use a 
%       Neumann boundary condition so we do not mask dFLDdt (see above also)

%3) include seed contributions in matrix:
    FLDkkFROM=convert2array(FLDkkFROM);
%3.1) for wet points
    dFLDdtWet=convert2array(dFLDdt.*mskWet);
    tmp1=find(dFLDdtWet~=0&~isnan(dFLDdtWet)); 
    dFLDdtWet=dFLDdtWet(tmp1); FLDkkFROMtmp=FLDkkFROM(tmp1); FLDkkTOtmp=KK(tmp1);
    A=A+sparse(LL(FLDkkTOtmp),LL(FLDkkFROMtmp),dFLDdtWet,nn,nn);
%3.2) for dry points (-> this part reflects the neumann boundary condition)
    dFLDdtDry=convert2array(dFLDdt.*mskDry);
    tmp1=find(dFLDdtDry~=0&~isnan(dFLDdtDry));
    dFLDdtDry=dFLDdtDry(tmp1); FLDkkFROMtmp=FLDkkFROM(tmp1);
    A=A+sparse(LL(FLDkkFROMtmp),LL(FLDkkFROMtmp),dFLDdtDry,nn,nn);

end; end; end;

%to check results: 
%figure; spy(A);

%4) solve for potential:
yy=convert2array(fld); yy=yy(find(KK~=0));
xx=A\yy;
yyFROMxx=A*xx; 

%5) prepare output:
fldDivPot=0*convert2array(fld); fldDivPot(find(KK~=0))=xx;
fldDivPot=convert2array(fldDivPot);
[fldUdiv,fldVdiv]=calc_T_grad(fldDivPot,0);
%unlike calc_T_grad, we work in grid point index, so we need to omit grid factors
fldUdiv=fldUdiv.*mygrid.DXC; fldVdiv=fldVdiv.*mygrid.DYC;

%to check the results:
%fld1=0*convert2array(fld); fld1(find(KK~=0))=xx;
%fld2=0*convert2array(fld); fld2(find(KK~=0))=yyFROMxx;
%fld3=convert2array(calc_UV_div(fldU,fldV,{}));%no grid scaling factor in div. or transports

warning('on','MATLAB:divideByZero');

1	gforget	1.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	gforget	1.2	warning('off','MATLAB:divideByZero');
21
22	gforget	1.1	fld=calc_UV_div(fldU,fldV,{});%no grid scaling factor in div. or transports
23			fld=fld.*mygrid.mskC(:,:,1);
24	gforget	1.2	fld(fld==0)=NaN;
25
26	gforget	1.1	mskWet=fld; mskWet(~isnan(fld))=1;
27			mskDry=1*isnan(mskWet); mskDry(mskDry==0)=NaN;
28
29			%check for domain edge points where no exchange is possible:
30			tmp1=mskWet; tmp1(:)=1; tmp2=exch_T_N(tmp1);
31			for iF=1:mskWet.nFaces;
32			tmp3=mskWet{iF}; tmp4=tmp2{iF};
33			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);
34			if ~isempty(find(isnan(tmp4)&~isnan(tmp3))); fprintf('warning: mask was modified\n'); end;
35			tmp3(isnan(tmp4))=NaN; mskWet{iF}=tmp3;
36			end;
37			%why does this matters... see set-ups with open boundary conditions...
38
39			%put 0 first guess if needed and switch land mask:
40			fld(find(isnan(fld)))=0; fld=fld.*mskWet;
41
42			%define mapping from global array to (no nan points) global vector
43			tmp1=convert2array(mskWet);
44			kk=find(~isnan(tmp1)); nn=length(kk);
45	gforget	1.3	KK=tmp1; KK(:)=0; KK(kk)=kk; KKfaces=convert2array(KK); %global array indices
46			LL=tmp1; LL(:)=0; LL(kk)=[1:nn]; LLfaces=convert2array(LL); %global vector indices
47	gforget	1.1
48			%note: diffsmooth2D_extrap_inv.m uses a Dirichlet boundary condition
49			% so that it uses a nan/1 mask in KK/LL // here we use a
50			% Neumann boundary condition so we use a 0/1 mask (see below also)
51
52			%form matrix problem:
53			A=sparse([],[],[],nn,nn,nn*5);
54			for iF=1:fld.nFaces; for ii=1:3; for jj=1:3;
55
56			%1) seed points (FLDones) and neighborhood of influence (FLDkkFROM)
57			FLDones=fld; FLDones(:)=0;
58			FLDones{iF}(ii:3:end,jj:3:end)=1;
59			FLDones(KKfaces==0)=0;
60
61			FLDkkFROMtmp=fld; FLDkkFROMtmp(:)=0;
62			FLDkkFROMtmp{iF}(ii:3:end,jj:3:end)=KKfaces{iF}(ii:3:end,jj:3:end);
63			FLDkkFROMtmp(find(isnan(fld)))=0;
64
65			FLDkkFROM=exch_T_N(FLDkkFROMtmp); FLDkkFROM(isnan(FLDkkFROM))=0;
66			for iF2=1:fld.nFaces;
67			tmp1=FLDkkFROM{iF2}; tmp2=zeros(size(tmp1)-2);
68			for ii2=1:3; for jj2=1:3; tmp2=tmp2+tmp1(ii2:end-3+ii2,jj2:end-3+jj2); end; end;
69			FLDkkFROM{iF2}=tmp2;
70			end;
71
72			%2) compute effect of each point on neighboring target point:
73			[tmpU,tmpV]=calc_T_grad(FLDones,0);
74			%unlike calc_T_grad, we work in grid point index, so we need to omit grid factors
75			tmpU=tmpU.mygrid.DXC; tmpV=tmpV.mygrid.DYC;
76			%and accordingly we use no grid scaling factor in div.
77			[dFLDdt]=calc_UV_div(tmpU,tmpV,{});
78
79			%note: diffsmooth2D_extrap_inv.m uses a Dirichlet boundary condition
80			% so that it needs to apply mskFreeze to dFLDdt // here we use a
81			% Neumann boundary condition so we do not mask dFLDdt (see above also)
82
83			%3) include seed contributions in matrix:
84			FLDkkFROM=convert2array(FLDkkFROM);
85			%3.1) for wet points
86			dFLDdtWet=convert2array(dFLDdt.*mskWet);
87			tmp1=find(dFLDdtWet~=0&~isnan(dFLDdtWet));
88			dFLDdtWet=dFLDdtWet(tmp1); FLDkkFROMtmp=FLDkkFROM(tmp1); FLDkkTOtmp=KK(tmp1);
89			A=A+sparse(LL(FLDkkTOtmp),LL(FLDkkFROMtmp),dFLDdtWet,nn,nn);
90			%3.2) for dry points (-> this part reflects the neumann boundary condition)
91			dFLDdtDry=convert2array(dFLDdt.*mskDry);
92			tmp1=find(dFLDdtDry~=0&~isnan(dFLDdtDry));
93			dFLDdtDry=dFLDdtDry(tmp1); FLDkkFROMtmp=FLDkkFROM(tmp1);
94			A=A+sparse(LL(FLDkkFROMtmp),LL(FLDkkFROMtmp),dFLDdtDry,nn,nn);
95
96			end; end; end;
97
98			%to check results:
99			%figure; spy(A);
100
101			%4) solve for potential:
102			yy=convert2array(fld); yy=yy(find(KK~=0));
103			xx=A\yy;
104			yyFROMxx=A*xx;
105
106			%5) prepare output:
107			fldDivPot=0*convert2array(fld); fldDivPot(find(KK~=0))=xx;
108	gforget	1.3	fldDivPot=convert2array(fldDivPot);
109	gforget	1.1	[fldUdiv,fldVdiv]=calc_T_grad(fldDivPot,0);
110			%unlike calc_T_grad, we work in grid point index, so we need to omit grid factors
111			fldUdiv=fldUdiv.mygrid.DXC; fldVdiv=fldVdiv.mygrid.DYC;
112
113			%to check the results:
114			%fld1=0*convert2array(fld); fld1(find(KK~=0))=xx;
115			%fld2=0*convert2array(fld); fld2(find(KK~=0))=yyFROMxx;
116			%fld3=convert2array(calc_UV_div(fldU,fldV,{}));%no grid scaling factor in div. or transports
117
118	gforget	1.2	warning('on','MATLAB:divideByZero');
119	gforget	1.1