| 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 |
% doStep is a field of 0 (tracer stays unchanged) |
| 12 |
% and 1 (tracer evolves). |
| 13 |
% 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 |
if isempty(who('fldRelax')); fldRelax=[]; end; |
| 24 |
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 |
%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 |
%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.25*msk.*(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.25*msk.*(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_conv(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 |
%mask of frozen/evolving points: |
| 140 |
dFLDdt=dFLDdt.*doStep; |
| 141 |
|
| 142 |
%now step forward |
| 143 |
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 100*nrm/nrmRef=%4.4g.',it,100*nrm(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 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph