matlab_class/gcmfaces_calc/gcmfaces_timestep.m

function [FLD]=gcmfaces_timestep(myOp,fld,fldRelax);
%object:    time step a 2D tracer diffusion/relaxation (could later be
%               extended with 3D case, advection and surface forcing)
%inputs:    myOp specifies the equation with the following fields
%               dt   time step
%               nbt  number of time steps
%               eps  termination criterium
%               Kux, Kvy diffusion coefficients (along grid axes; m2/s)
%               tau  relaxation time scale
%               Kuy, Kvx diffusion coefficients (only for rotated diff)
%               doStep   is a field of 0 (tracer stays unchanged) 
%                        and 1 (tracer evolves).
%           fld is the initial tracer field
%           fldRelax is the field to relax to (if tau<Inf)
%output:    FLD is the final tracer field
%
%note : case of 3D tracer is not coded yet
%note : advection (m/s) and forcing (tracer/area/s) are not coded yet


gcmfaces_global;

if isempty(who('fldRelax'));   fldRelax=[]; end;
if isempty(who('fldForcing')); fldForcing=[]; end;

%rotated diffusion:
if isfield(myOp,'Kuy')&isfield(myOp,'Kvx');
    if ~isempty(myOp.Kuy)&~isempty(myOp.Kvx);
        doExtraDiag=1;
        if myenv.verbose;
            gcmfaces_msg(['* extra-diagonal diffusion included.']);
        end;
    else;
        doExtraDiag=0;
    end;
else;
    doExtraDiag=0;
end;

%relaxation term:
if isfield(myOp,'tau');
    doRelax=1;
    if isempty(fldRelax);
        error('relaxation field must be specified');
    end;
    if ~isa(myOp.tau,'double')&~isa(myOp.tau,'gcmfaces');
        error('mispecified relaxation time scale (must be doule or gcmfaces)');
    end;
    
    if myenv.verbose;
        gcmfaces_msg(['* relaxation term included.']);
    end;
    test1=1*(myOp.tau<myOp.dt);
    if isa(myOp.tau,'gcmfaces'); test1=convert2vector(test1); end;
    test1=~isempty(find(test1==1));
    if test1;
        error('tau cannot exceed 1');
    end;
else;
    doRelax=0;
end;

%mask of frozen/evolving points:
if isfield(myOp,'doStep');
    doStep=myOp.doStep;
    test1=sum(doStep~=1&doStep~=0);
    if test1;
        error('myOp.doStep must be a field of 0/1');
    end;
else;
    doStep=mygrid.XC; doStep(:)=1;
end;

%forcing term:
if ~isempty(fldForcing);
    doForcing=1;
    if myenv.verbose;
        gcmfaces_msg(['* forcing term included.']);
    end;
else;
    doForcing=0;
end;

%
if isfield(myOp,'eps');
    nbt=Inf;
    eps=myOp.eps;
    if myenv.verbose;
        gcmfaces_msg(['* specified termination : ' num2str(eps) ' .']);
    end;
else;
    nbt=myOp.nbt;
    eps=0;
    if myenv.verbose;
        gcmfaces_msg(['* specified number of time steps : ' num2str(nbt) ' .']);
    end;
end;

%initialization
FLD=fld; it=0; doStop=0; nrm=[]; nrmRef=Inf;

%time-stepping loop:
while ~doStop;
    
    it=it+1;
    
    [dTdxAtU,dTdyAtV]=calc_T_grad(FLD,0);
    tmpU=dTdxAtU.*myOp.Kux;
    tmpV=dTdyAtV.*myOp.Kvy;
    
    if doExtraDiag;
        dTdyAtU=dTdxAtU; dTdxAtV=dTdyAtV;
        [dTdxAtU,dTdyAtV]=exch_UV_N(dTdxAtU,dTdyAtV);
        for iF=1:FLD.nFaces;
            msk=dTdxAtU{iF}(2:end-1,2:end-1); msk(~isnan(msk))=1;
            tmp1=dTdyAtV{iF}; tmp1(isnan(tmp1))=0;
            dTdyAtU{iF}=0.25*msk.*(tmp1(1:end-2,2:end-1)+tmp1(1:end-2,3:end)+...
                tmp1(2:end-1,2:end-1)+tmp1(2:end-1,3:end));
            
            msk=dTdyAtV{iF}(2:end-1,2:end-1); msk(~isnan(msk))=1;
            tmp1=dTdxAtU{iF}; tmp1(isnan(tmp1))=0;
            dTdxAtV{iF}=0.25*msk.*(tmp1(2:end-1,1:end-2)+tmp1(3:end,1:end-2)+...
                tmp1(2:end-1,2:end-1)+tmp1(3:end,2:end-1));
        end;
        dTdxAtU=cut_T_N(dTdxAtU);
        dTdyAtV=cut_T_N(dTdyAtV);
        
        tmpU=tmpU+dTdyAtU.*myOp.Kuy;
        tmpV=tmpV+dTdxAtV.*myOp.Kvx;
    end;
    
    [fldDIV]=calc_UV_div(tmpU,tmpV,{'dh'});
    dFLDdt=-myOp.dt*fldDIV./mygrid.RAC;
    
    if doRelax;
        dFLDdt=dFLDdt+myOp.dt*(fldRelax-FLD)./myOp.tau;
    end;
    
    %mask of frozen/evolving points:
    dFLDdt=dFLDdt.*doStep;
    
    %now step forward
    FLD=FLD+dFLDdt;
    
    %test for termination criteria
    nrm=[nrm;sqrt(nanmean(dFLDdt.^2))];
    if it==1; nrmRef=nrm(it); end;
    doStop=((nrm(it)<nrmRef*eps)|(it==nbt));
    
    %monitor convergence
    if mod(it,50)==0|it==1|doStop;
        tmp1=sprintf('it=%04i  100*nrm/nrmRef=%4.4g.',it,100*nrm(it)/nrmRef);
        if myenv.verbose;
            gcmfaces_msg(['* convergence monitor : ' tmp1]);
        end;        
    end;
    
    if 0;%monitor convergence
        if mod(it,1000)==0|doStop;
            figure; plot(log10(nrm)); ylabel('log10(increment)');
            eval(['FLD_' num2str(it) '=FLD;']);
        end;
    end;
    
    if 0;%test of positivity
        test1=convert2vector(FLD<0); test1=~isempty(find(test1(:)==1));
        if test1;
            fprintf(['it=' num2str(it) ' min(FLD)=' num2str(min(FLD)) '\n']);
        end;
    end;
    
end;

1	gforget	1.1	function [FLD]=gcmfaces_timestep(myOp,fld,fldRelax);
2			%object: time step a 2D tracer diffusion/relaxation (could later be
3			% extended with 3D case, advection and surface forcing)
4			%inputs: myOp specifies the equation with the following fields
5			% dt time step
6			% nbt number of time steps
7			% eps termination criterium
8			% Kux, Kvy diffusion coefficients (along grid axes; m2/s)
9			% tau relaxation time scale
10			% Kuy, Kvx diffusion coefficients (only for rotated diff)
11	gforget	1.2	% doStep is a field of 0 (tracer stays unchanged)
12			% and 1 (tracer evolves).
13	gforget	1.1	% fld is the initial tracer field
14			% fldRelax is the field to relax to (if tau<Inf)
15			%output: FLD is the final tracer field
16			%
17			%note : case of 3D tracer is not coded yet
18			%note : advection (m/s) and forcing (tracer/area/s) are not coded yet
19
20
21			gcmfaces_global;
22
23	gforget	1.2	if isempty(who('fldRelax')); fldRelax=[]; end;
24	gforget	1.1	if isempty(who('fldForcing')); fldForcing=[]; end;
25
26			%rotated diffusion:
27			if isfield(myOp,'Kuy')&isfield(myOp,'Kvx');
28			if ~isempty(myOp.Kuy)&~isempty(myOp.Kvx);
29			doExtraDiag=1;
30			if myenv.verbose;
31			gcmfaces_msg(['* extra-diagonal diffusion included.']);
32			end;
33			else;
34			doExtraDiag=0;
35			end;
36			else;
37			doExtraDiag=0;
38			end;
39
40			%relaxation term:
41			if isfield(myOp,'tau');
42			doRelax=1;
43			if isempty(fldRelax);
44			error('relaxation field must be specified');
45			end;
46			if ~isa(myOp.tau,'double')&~isa(myOp.tau,'gcmfaces');
47			error('mispecified relaxation time scale (must be doule or gcmfaces)');
48			end;
49
50			if myenv.verbose;
51			gcmfaces_msg(['* relaxation term included.']);
52			end;
53			test1=1*(myOp.tau<myOp.dt);
54			if isa(myOp.tau,'gcmfaces'); test1=convert2vector(test1); end;
55			test1=~isempty(find(test1==1));
56			if test1;
57			error('tau cannot exceed 1');
58			end;
59			else;
60			doRelax=0;
61			end;
62
63	gforget	1.2	%mask of frozen/evolving points:
64			if isfield(myOp,'doStep');
65			doStep=myOp.doStep;
66			test1=sum(doStep~=1&doStep~=0);
67			if test1;
68			error('myOp.doStep must be a field of 0/1');
69			end;
70			else;
71			doStep=mygrid.XC; doStep(:)=1;
72			end;
73
74	gforget	1.1	%forcing term:
75			if ~isempty(fldForcing);
76			doForcing=1;
77			if myenv.verbose;
78			gcmfaces_msg(['* forcing term included.']);
79			end;
80			else;
81			doForcing=0;
82			end;
83
84			%
85			if isfield(myOp,'eps');
86			nbt=Inf;
87			eps=myOp.eps;
88			if myenv.verbose;
89			gcmfaces_msg(['* specified termination : ' num2str(eps) ' .']);
90			end;
91			else;
92			nbt=myOp.nbt;
93			eps=0;
94			if myenv.verbose;
95			gcmfaces_msg(['* specified number of time steps : ' num2str(nbt) ' .']);
96			end;
97			end;
98
99			%initialization
100			FLD=fld; it=0; doStop=0; nrm=[]; nrmRef=Inf;
101
102			%time-stepping loop:
103			while ~doStop;
104
105			it=it+1;
106
107			[dTdxAtU,dTdyAtV]=calc_T_grad(FLD,0);
108			tmpU=dTdxAtU.*myOp.Kux;
109			tmpV=dTdyAtV.*myOp.Kvy;
110
111			if doExtraDiag;
112			dTdyAtU=dTdxAtU; dTdxAtV=dTdyAtV;
113			[dTdxAtU,dTdyAtV]=exch_UV_N(dTdxAtU,dTdyAtV);
114			for iF=1:FLD.nFaces;
115			msk=dTdxAtU{iF}(2:end-1,2:end-1); msk(~isnan(msk))=1;
116			tmp1=dTdyAtV{iF}; tmp1(isnan(tmp1))=0;
117			dTdyAtU{iF}=0.25msk.(tmp1(1:end-2,2:end-1)+tmp1(1:end-2,3:end)+...
118			tmp1(2:end-1,2:end-1)+tmp1(2:end-1,3:end));
119
120			msk=dTdyAtV{iF}(2:end-1,2:end-1); msk(~isnan(msk))=1;
121			tmp1=dTdxAtU{iF}; tmp1(isnan(tmp1))=0;
122			dTdxAtV{iF}=0.25msk.(tmp1(2:end-1,1:end-2)+tmp1(3:end,1:end-2)+...
123			tmp1(2:end-1,2:end-1)+tmp1(3:end,2:end-1));
124			end;
125			dTdxAtU=cut_T_N(dTdxAtU);
126			dTdyAtV=cut_T_N(dTdyAtV);
127
128			tmpU=tmpU+dTdyAtU.*myOp.Kuy;
129			tmpV=tmpV+dTdxAtV.*myOp.Kvx;
130			end;
131
132			[fldDIV]=calc_UV_div(tmpU,tmpV,{'dh'});
133			dFLDdt=-myOp.dt*fldDIV./mygrid.RAC;
134
135			if doRelax;
136			dFLDdt=dFLDdt+myOp.dt*(fldRelax-FLD)./myOp.tau;
137			end;
138
139	gforget	1.2	%mask of frozen/evolving points:
140			dFLDdt=dFLDdt.*doStep;
141
142			%now step forward
143	gforget	1.1	FLD=FLD+dFLDdt;
144
145			%test for termination criteria
146			nrm=[nrm;sqrt(nanmean(dFLDdt.^2))];
147			if it==1; nrmRef=nrm(it); end;
148			doStop=((nrm(it)<nrmRef*eps)\|(it==nbt));
149
150			%monitor convergence
151			if mod(it,50)==0\|it==1\|doStop;
152			tmp1=sprintf('it=%04i 100nrm/nrmRef=%4.4g.',it,100nrm(it)/nrmRef);
153			if myenv.verbose;
154			gcmfaces_msg(['* convergence monitor : ' tmp1]);
155			end;
156			end;
157
158			if 0;%monitor convergence
159			if mod(it,1000)==0\|doStop;
160			figure; plot(log10(nrm)); ylabel('log10(increment)');
161			eval(['FLD_' num2str(it) '=FLD;']);
162			end;
163			end;
164
165			if 0;%test of positivity
166			test1=convert2vector(FLD<0); test1=~isempty(find(test1(:)==1));
167			if test1;
168			fprintf(['it=' num2str(it) ' min(FLD)=' num2str(min(FLD)) '\n']);
169			end;
170			end;
171
172			end;
173