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