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');
warning('off','MATLAB:nearlySingularMatrix');

fld=calc_UV_conv(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_conv(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_conv(fldU,fldV,{}));%no grid scaling factor in div. or transports

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

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	gforget	1.5	warning('off','MATLAB:nearlySingularMatrix');
22	gforget	1.2
23	gforget	1.4	fld=calc_UV_conv(fldU,fldV,{});%no grid scaling factor in div. or transports
24	gforget	1.1	fld=fld.*mygrid.mskC(:,:,1);
25	gforget	1.2	fld(fld==0)=NaN;
26
27	gforget	1.1	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	gforget	1.3	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	gforget	1.1
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	gforget	1.4	[dFLDdt]=calc_UV_conv(tmpU,tmpV,{});
79	gforget	1.1
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	gforget	1.3	fldDivPot=convert2array(fldDivPot);
110	gforget	1.1	[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	gforget	1.4	%fld3=convert2array(calc_UV_conv(fldU,fldV,{}));%no grid scaling factor in div. or transports
118	gforget	1.1
119	gforget	1.2	warning('on','MATLAB:divideByZero');
120	gforget	1.5	warning('on','MATLAB:nearlySingularMatrix');
121	gforget	1.1