source: issm/branches/trunk-larour-SLPS2022/src/m/classes/sealevelmodel.m@ 27124

Last change on this file since 27124 was 27124, checked in by Eric.Larour, 3 years ago

CHG: new diagnostics routine.
File size: 18.1 KB
RevLine 
[20110]1%SEALEVELMODEL class definition
2%
3% Usage:
4% slm = sealevelmodel(varargin)
5%
[26358]6% where varargin is a variable list of options:
[20110]7%
[26358]8% Example:
[20110]9% slm = sealevel('icecap',md_greenland,'icecap',md_antarctica,'earth',md_earth);
10
[24156]11classdef sealevelmodel < handle
[20110]12 properties (SetAccess=public) %Model fields
13 % {{{
[24728]14 icecaps = {}; % list of land/ice models, name should be change longer term.
[20110]15 earth = 0; % model for the whole earth
[26307]16 basins = {}; % list of basins, matching icecaps, where shapefile info is held
[20110]17 cluster = 0;
18 miscellaneous = 0;
19 settings = 0;
20 private = 0;
[24156]21 mergedcaps = 0;
22 transitions = {};
[25455]23 eltransitions = {};
[24779]24 planet = '';
[20110]25 %}}}
26 end
27 methods
28 function slm = sealevelmodel(varargin) % {{{
[25688]29 slm=setdefaultparameters(slm);
[20110]30
[25688]31 if nargin==1,
[20110]32
[25688]33 options=pairoptions(varargin{:});
34
[26358]35 %recover all the icecap models:
36 slm.icecaps=getfieldvalues(options,'ice_cap',{});
[20110]37
38 %recover the earth model:
[24779]39 slm.earth = getfieldvalue(options,'earth',0);
40
[26358]41 %set planet type:
[24779]42 slm.planet=getfieldvalue(options,'planet','earth');
43
[20110]44 end
45 end
46 %}}}
47 function checkconsistency(slm,solutiontype) % {{{
48
[26358]49 %is the coupler turned on?
[20110]50 for i=1:length(slm.icecaps),
51 if slm.icecaps{i}.transient.iscoupler==0,
[26352]52 warning(sprintf('sealevelmodel.m::checkconsistency: icecap model %s should have the transient coupler option turned on!',slm.icecaps{i}.miscellaneous.name));
[20110]53 end
54 end
55
56 if slm.earth.transient.iscoupler==0,
[26352]57 warning('sealevelmodel.m::checkconsistency: earth model should have the transient coupler option turned on!');
[20110]58 end
59
[26358]60 %check that the transition vectors have the right size:
[20131]61 for i=1:length(slm.icecaps),
[25956]62 if slm.icecaps{i}.mesh.numberofvertices ~= length(slm.earth.solidearth.transitions{i}),
[26352]63 error(['sealevelmodel.m::checkconsistency: issue with size of transition vector for ice cap: ' num2str(i) ' name: ' slm.icecaps{i}.miscellaneous.name]);
[20131]64 end
65 end
[22955]66
[26358]67 %check that runfrequency is the same everywhere:
[22955]68 for i=1:length(slm.icecaps),
[25956]69 if slm.icecaps{i}.solidearth.settings.runfrequency~=slm.earth.solidearth.settings.runfrequency,
[26352]70 error(sprintf('sealevelmodel.m::checkconsistency: icecap model %s should have the same run frequency as earth!',slm.icecaps{i}.miscellaneous.name));
[22955]71 end
72 end
[20110]73
[26358]74 %make sure steric_rate is the same everywhere:
[22955]75 for i=1:length(slm.icecaps),
[26358]76 md= slm.icecaps{i};
[25956]77 if ~isempty(find(md.dsl.steric_rate - slm.earth.dsl.steric_rate(slm.earth.dsl.transitions{i}))),
[26352]78 error(sprintf('sealevelmodel.m::checkconsistency: steric rate on ice cap %s is not the same as for the earth\n',md.miscellaneous.name));
[22955]79 end
80 end
[20131]81
[26047]82 %make sure grd is the same everywhere:
83 for i=1:length(slm.icecaps),
[26358]84 md= slm.icecaps{i};
[26080]85 if md.solidearthsettings.isgrd~=slm.earth.solidearthsettings.isgrd
[26352]86 error(sprintf('sealevelmodel.m::checkconsistency: isgrd on ice cap %s is not the same as for the earth\n',md.miscellaneous.name));
[26047]87 end
88 end
[26307]89
[26358]90 %make sure that there is no solid earth external forcing on the basins:
[26047]91 for i=1:length(slm.icecaps),
[26358]92 md= slm.icecaps{i};
[26047]93 if ~isempty(md.solidearth.external),
[26352]94 error('sealevelmodel.m::checkconsistency: cannot run external forcings on an ice sheet when running a coupling earth/ice sheet model');
[26047]95 end
96
97 end
[26358]98 %make sure that we have the right grd model for computing out sealevel patterns:
[26047]99 for i=1:length(slm.icecaps),
[26358]100 md= slm.icecaps{i};
[26047]101 if md.solidearth.settings.grdmodel~=0
[26352]102 error(sprintf('sealevelmodel.m::checkconsistency: ice sheets do not run GRD module, specify solidearth.settings.grdmodel=0 on ice cap %i',i));
[26047]103 end
104 end
105
[20110]106 end
107 %}}}
108 function slm = setdefaultparameters(slm) % {{{
109
110 %initialize subclasses
111 slm.icecaps = {};
112 slm.earth = {};
113 slm.miscellaneous = miscellaneous();
[24165]114 slm.settings = issmsettings();
[20110]115 slm.private = private();
116 slm.cluster = generic();
[24156]117 slm.transitions = {};
[24779]118 slm.eltransitions = {};
119 slm.planet = 'earth';
[20110]120 end
121 %}}}
122 function disp(self) % {{{
123 disp(sprintf('%19s: %-22s -- %s','icecaps' ,['[' num2str(length(self.icecaps)) 'x1 ' class(self.icecaps) ']'],'ice caps'));
124 disp(sprintf('%19s: %-22s -- %s','earth' ,['[1x1 ' class(self.earth) ']'],'earth'));
125 disp(sprintf('%19s: %-22s -- %s','settings' ,['[1x1 ' class(self.settings) ']'],'settings properties'));
126 disp(sprintf('%19s: %-22s -- %s','cluster' ,['[1x1 ' class(self.cluster) ']'],'cluster parameters (number of cpus...)'));
127 disp(sprintf('%19s: %-22s -- %s','miscellaneous' ,['[1x1 ' class(self.miscellaneous) ']'],'miscellaneous fields'));
128 end % }}}
[22955]129 function self=mergeresults(self) % {{{
[25502]130 champs=fieldnames(self.icecaps{1}.results.TransientSolution);
[22955]131 for i=1:length(self.mergedcaps)/2,
132 md=self.mergedcaps{2*(i-1)+1}; trans=self.mergedcaps{2*(i-1)+2};
[24886]133 %icecaps=self.icecaps(self.range{2*(i-1)+2});
[22955]134 for j=1:length(self.icecaps{1}.results.TransientSolution),
135 for k=1:length(champs),
136 if strcmpi(class(icecaps{1}.results.TransientSolution(j).(champs{k})),'double'),
[26358]137 %vertex or element?
[22955]138 if length(icecaps{1}.results.TransientSolution(j).(champs{k}))==icecaps{1}.mesh.numberofvertices,
139 md.results.TransientSolution(j).(champs{k})=zeros(md.mesh.numberofvertices,1);
140 for l=1:length(trans),
141 resultcap=icecaps{l}.results.TransientSolution(j).(champs{k});
142 md.results.TransientSolution(j).(champs{k})(trans{l})=resultcap;
143 end
[26358]144 else
[22955]145 if strcmpi(champs{k},'IceVolume') | strcmpi(champs{k},'IceVolumeAboveFloatation') ,
146 md.results.TransientSolution(j).(champs{k})=0;
147 for l=1:length(trans),
148 resultcap=icecaps{l}.results.TransientSolution(j).(champs{k});
149 md.results.TransientSolution(j).(champs{k})= md.results.TransientSolution(j).(champs{k})+resultcap;
150 end
151 elseif strcmpi(champs{k},'time'),
152 md.results.TransientSolution(j).(champs{k})= icecaps{1}.results.TransientSolution(j).(champs{k});
153 else
154 continue;
155 end
156 end
157 else
158 continue;
159 end
[24156]160 end
[22955]161 end
162 self.mergedcaps{2*(i-1)+1}=md;
163 end
164 end % }}}
165 function listcaps(self) % {{{
166 for i=1:length(self.icecaps),
167 disp(sprintf('%i: %s',i,self.icecaps{i}.miscellaneous.name));
168 end
169 end % }}}
[25677]170 function n=ncaps(self) % {{{
171 n=length(self.icecaps);
[25675]172 end % }}}
[24886]173 function list=continents(self) % {{{
174 list={};
175 for i=1:length(self.basins),
176 list{end+1}=self.basins{i}.continent;
177 end
178 list=unique(list);
179 end % }}}
180 function list=basinsfromcontinent(self,continent) % {{{
181 list={};
182 for i=1:length(self.icecaps),
183 if strcmpi(self.basins{i}.continent,continent),
184 list{end+1}=self.basins{i}.name;
185 end
186 end
187 list=unique(list);
188 end % }}}
[24156]189 function addbasin(self,bas) % {{{
[25065]190 if ~strcmpi(class(bas),'basin')
191 error('addbasin method only takes a ''basin'' class object as input');
192 end;
193 self.basins{end+1}=bas;
[24156]194 end % }}}
[25620]195 function intersections2d(self,varargin) % {{{
196
[26358]197 options=pairoptions(varargin{:});
[25620]198 force=getfieldvalue(options,'force',0);
199
[26358]200 %initialize, to avoid issues of having more transitions than meshes.
[25620]201 self.transitions={};
202 self.eltransitions={};
203
204 %for elements:
205 xe=self.earth.mesh.x(self.earth.mesh.elements)*[1;1;1]/3;
206 ye=self.earth.mesh.y(self.earth.mesh.elements)*[1;1;1]/3;
207
208 for i=1:length(self.icecaps),
209 mdi=self.icecaps{i};
210
[26358]211 %for elements:
[25620]212 xei=mdi.mesh.x(mdi.mesh.elements)*[1;1;1]/3;
213 yei=mdi.mesh.y(mdi.mesh.elements)*[1;1;1]/3;
214
215 disp(sprintf('Computing vertex intersections for basin %s',self.basins{i}.name));
216
217 self.transitions{end+1}=meshintersect2d(self.earth.mesh.x,self.earth.mesh.y,mdi.mesh.x,mdi.mesh.y,'force',force);
218
219 self.eltransitions{end+1}=meshintersect2d(xe,ye,xei,yei,'force',force);
220 end
221 end % }}}
[24156]222 function intersections(self,varargin) % {{{
223
[26358]224 options=pairoptions(varargin{:});
[25455]225 force=getfieldvalue(options,'force',0);
226
[26358]227 %initialize, to avoid issues of having more transitions than meshes.
[25499]228 self.transitions={};
229 self.eltransitions={};
[24156]230
[25455]231 %for elements:
232 xe=self.earth.mesh.x(self.earth.mesh.elements)*[1;1;1]/3;
233 ye=self.earth.mesh.y(self.earth.mesh.elements)*[1;1;1]/3;
234 ze=self.earth.mesh.z(self.earth.mesh.elements)*[1;1;1]/3;
[25499]235
[25455]236 for i=1:length(self.icecaps),
237 mdi=self.icecaps{i};
238 mdi=TwoDToThreeD(mdi,self.planet);
[27058]239
240
[26358]241 %for elements:
[25455]242 xei=mdi.mesh.x(mdi.mesh.elements)*[1;1;1]/3;
243 yei=mdi.mesh.y(mdi.mesh.elements)*[1;1;1]/3;
244 zei=mdi.mesh.z(mdi.mesh.elements)*[1;1;1]/3;
245
246 disp(sprintf('Computing vertex intersections for basin %s',self.basins{i}.name));
247
248 self.transitions{end+1}=meshintersect3d(self.earth.mesh.x,self.earth.mesh.y,self.earth.mesh.z,mdi.mesh.x,mdi.mesh.y,mdi.mesh.z,'force',force);
249 self.eltransitions{end+1}=meshintersect3d(xe,ye,ze,xei,yei,zei,'force',force);
250 end
[24156]251 end % }}}
252 function checkintersections(self) % {{{
[25455]253 flags=zeros(self.earth.mesh.numberofvertices,1);
254 for i=1:length(self.basins),
255 flags(self.transitions{i})=i;
256 end
257 plotmodel(self.earth,'data',flags,'coastline','on');
[24156]258
259 end % }}}
[24779]260 function checkbasinconsistency(self) % {{{
261 for i=1:length(self.basins),
262 self.basins{i}.checkconsistency();
263 end
264
265 end % }}}
[24156]266 function baslist=basinindx(self,varargin) % {{{
[26358]267 options=pairoptions(varargin{:});
[24156]268 continent=getfieldvalue(options,'continent','all');
269 bas=getfieldvalue(options,'basin','all');
270
271 %expand continent list: {{{
272 if iscell(continent),
273 if length(continent)==1,
[25125]274 if strcmpi(continent{1},'all'),
[26358]275 %need to transform this into a list of continents:
276 continent={};
[25125]277 for i=1:length(self.basins),
278 continent{end+1}=self.basins{i}.continent;
279 end
280 continent=unique(continent);
[25455]281 else
282 %nothing to do: assume we have a list of continents
[25125]283 end
[24156]284 else
[25455]285 %nothing to do: assume we have a list of continents
[24156]286 end
287 else
288 if strcmpi(continent,'all'),
[26358]289 %need to transform this into a list of continents:
290 continent={};
[24156]291 for i=1:length(self.basins),
[25125]292 continent{end+1}=self.basins{i}.continent;
[24156]293 end
294 continent=unique(continent);
295 else
296 continent={continent};
297 end
298 end
299 %}}}
300 %expand basins list using the continent list above and the extra bas discriminator: %{{{
301 if iscell(bas),
302 if length(bas)==1,
[25125]303 if strcmpi(bas{1},'all'),
[26358]304 %need to transform this into a list of basins:
[25125]305 baslist=[];
306 for i=1:length(self.basins),
307 if self.basins{i}.iscontinentany(continent{:}),
308 baslist(end+1)=i;
309 end
310 end
311 baslist=unique(baslist);
312 else
313 bas=bas{1};
314 baslist=[];
315 for i=1:length(self.basins),
316 if self.basins{i}.iscontinentany(continent{:}),
317 if self.basins{i}.isnameany(bas),
318 baslist(end+1)=i;
319 end
320 end
321 end
[24156]322
[25125]323 end
[24156]324 else
[26358]325 %we have a list of basin names:
[24156]326 baslist=[];
327 for i=1:length(bas),
328 basname=bas{i};
329 for j=1:length(self.basins),
330 if self.basins{j}.iscontinentany(continent{:}),
331 if self.basins{j}.isnameany(basname),
332 baslist(end+1)=j;
333 end
334 end
335 end
336 baslist=unique(baslist);
337 end
338 end
339 else
340 if strcmpi(bas,'all'),
341 baslist=[];
342 for i=1:length(self.basins),
343 if self.basins{i}.iscontinentany(continent{:}),
344 baslist(end+1)=i;
345 end
346 end
347 baslist=unique(baslist);
348 else
349 baslist=[];
350 for i=1:length(self.basins),
351 if self.basins{i}.iscontinentany(continent{:}),
352 if self.basins{i}.isnameany(bas),
353 baslist(end+1)=i;
354 end
355 end
356 end
357 baslist=unique(baslist);
358 end
359 end
360 %}}}
361
362 end % }}}
363 function addicecap(self,md) % {{{
[25136]364 if ~strcmpi(class(md),'model')
365 error('addicecap method only takes a ''model'' class object as input');
366 end
367 self.icecaps{end+1}=md;
[24156]368 end % }}}
[27124]369 function icecapsinfo(self,varargin) % {{{
370 totalnel=0;
371 disp('Sealevel model basins mesh info:');
372 for i=1:length(self.icecaps),
373 md=self.icecaps{i};
374 disp(sprintf(' %i(%s): %i',i,md.miscellaneous.name,md.mesh.numberofelements));
375 totalnel=totalnel+md.mesh.numberofelements;
376 end
377 disp(sprintf('Total number of elements:%i',totalnel));
378 end % }}}
[24156]379 function basinsplot3d(self,varargin) % {{{
[25136]380 for i=1:length(self.basins),
381 self.basins{i}.plot3d(varargin{:});
382 end
[24156]383 end % }}}
384 function caticecaps(self,varargin) % {{{
385
[26358]386 %recover options:
387 options=pairoptions(varargin{:});
[24156]388 tolerance=getfieldvalue(options,'tolerance',.65);
389 loneedgesdetect=getfieldvalue(options,'loneedgesdetect',0);
390
391 %make 3D model:
392 models=self.icecaps;
393 for i=1:length(models),
[24779]394 models{i}=TwoDToThreeD(models{i},self.planet);
[24156]395 end
396
397 %Plug all models together:
[26358]398 md=models{1};
[24156]399 for i=2:length(models),
400 md=modelmerge3d(md,models{i},'tolerance',tolerance);
401 md.private.bamg.landmask=[md.private.bamg.landmask;models{i}.private.bamg.landmask];
402 end
403
404 %Look for lone edges if asked for it: {{{
405 if loneedgesdetect,
406 edges=loneedges(md);
407 plotmodel(md,'data',md.mask.land_levelset);
408 hold on;
409 for i=1:length(edges),
410 ind1=edges(i,1);
411 ind2=edges(i,2);
412 %plot([md.mesh.x(ind1),md.mesh.x(ind2)],[md.mesh.y(ind1),md.mesh.y(ind2)],'r*-');
413 plot3([md.mesh.x(ind1),md.mesh.x(ind2)],[md.mesh.y(ind1),md.mesh.y(ind2)],[md.mesh.z(ind1),md.mesh.z(ind2)],'g*-');
414 end
415 end %}}}
416
[26358]417 %Plug into earth:
[24156]418 self.earth=md;
419
[24886]420 %Create mesh radius:
421 self.earth.mesh.r=planetradius('earth')*ones(md.mesh.numberofvertices,1);
[27058]422 self.earth.mesh.lat = asind(self.earth.mesh.z./self.earth.mesh.r);
423 self.earth.mesh.long = atan2d(self.earth.mesh.y,self.earth.mesh.x);
[24886]424
[27058]425
[24156]426 end % }}}
[25620]427 function caticecaps2d(self,varargin) % {{{
428
[26358]429 %recover options:
430 options=pairoptions(varargin{:});
[25620]431 tolerance=getfieldvalue(options,'tolerance',1e-5);
432 loneedgesdetect=getfieldvalue(options,'loneedgesdetect',0);
433 models=self.icecaps;
434
435 %Plug all models together:
[26358]436 md=models{1};
[25620]437 for i=2:length(models),
438 md=modelmerge2d(md,models{i},'tolerance',tolerance);
439 end
440
441 %Look for lone edges if asked for it: {{{
442 if loneedgesdetect,
443 edges=loneedges(md);
444 hold on;
445 for i=1:length(edges),
446 ind1=edges(i,1);
447 ind2=edges(i,2);
448 plot([md.mesh.x(ind1),md.mesh.x(ind2)],[md.mesh.y(ind1),md.mesh.y(ind2)],'g*-');
[22955]449 end
[25620]450 end %}}}
451
[26358]452 %Plug into earth:
[25620]453 self.earth=md;
454
[22955]455 end % }}}
[25688]456 function viscousiterations(self) % {{{
457 for i=1:length(self.icecaps),
458 ic=self.icecaps{i};
459 mvi=ic.results.TransientSolution(1).StressbalanceConvergenceNumSteps;
460 for j=2:length(ic.results.TransientSolution)-1,
461 mvi=max(mvi,ic.results.TransientSolution(j).StressbalanceConvergenceNumSteps);
462 end
463 disp(sprintf('%i, %s: %i',i,self.icecaps{i}.miscellaneous.name,mvi));
464 end
465 end % }}}
[22955]466 function maxtimestep(self) % {{{
467 for i=1:length(self.icecaps),
[25688]468 ic=self.icecaps{i};
[22955]469 mvi=length(ic.results.TransientSolution);
470 timei=ic.results.TransientSolution(end).time;
471 disp(sprintf('%i, %s: %i/%g',i,self.icecaps{i}.miscellaneous.name,mvi,timei));
472 end
473 mvi=length(self.earth.results.TransientSolution);
474 timei=self.earth.results.TransientSolution(end).time;
475 disp(sprintf('Earth: %i/%g',mvi,timei));
476 end % }}}
[24165]477 function transfer(self,string) % {{{
[26358]478 %Recover field size in one icecap:
[24779]479 eval(['n=size(self.icecaps{1}.' string ',1);']);
[24165]480 if n==self.icecaps{1}.mesh.numberofvertices,
481 eval(['self.earth.' string '=zeros(self.earth.mesh.numberofvertices,1);']);
482 for i=1:length(self.icecaps),
483 eval(['self.earth.' string '(self.transitions{' num2str(i) '})=self.icecaps{' num2str(i) '}.' string ';']);
484 end
[24779]485 elseif n==(self.icecaps{1}.mesh.numberofvertices+1),
[26358]486 %dealing with a transient dataset.
[24779]487 %check that all timetags are similar between all icecaps: %{{{
488 for i=1:length(self.icecaps),
489 eval(['capfieldi= self.icecaps{' num2str(i) '}.' string ';']);
490 for j=(i+1):length(self.icecaps),
491 eval(['capfieldj= self.icecaps{' num2str(j) '}.' string ';']);
492 if ~isequal(capfieldi(end,:),capfieldj(end,:)),
[25499]493 error(['Time stamps for ' string ' field are different between icecaps ' num2str(i) ' and ' num2str(j)]);
[24779]494 end
495 end
496 end
[26358]497 eval(['capfield1= self.icecaps{1}.' string ';']);
[24779]498 times=capfield1(end,:);
499 nsteps=length(times);
500 %}}}
501 %initialize: %{{{
502 eval(['field=zeros(self.earth.mesh.numberofvertices+1,' num2str(nsteps) ');']);
503 eval(['field(end,:)=times;']); %transfer the times only, not the values
504 %}}}
505 %transfer all time fields: {{{
506 for i=1:length(self.icecaps),
507 eval(['capfieldi= self.icecaps{' num2str(i) '}.' string ';']);
508 for j=1:nsteps,
509 eval(['field(self.transitions{' num2str(i) '},' num2str(j) ')=capfieldi(1:end-1,' num2str(j) ');']); %transfer only the values, not the time.
510 end
511 end
512 eval(['self.earth.' string '=field;']); %do not forget to plug the field variable into its final location
513 %}}}
[25455]514 elseif n==self.icecaps{1}.mesh.numberofelements,
[24165]515 eval(['self.earth.' string '=zeros(self.earth.mesh.numberofelements,1);']);
516 for i=1:length(self.icecaps),
517 eval(['self.earth.' string '(self.eltransitions{' num2str(i) '})=self.icecaps{' num2str(i) '}.' string ';']);
518 end
519 else
520 error('not supported yet');
521 end
522 end % }}}
[22955]523 function self=homogeneize(self,noearth) % {{{
524 if nargin==1,
525 noearth=0;
526 end
527 mintimestep=Inf;
528 for i=1:length(self.icecaps),
[26358]529 ic=self.icecaps{i};
[22955]530 mintimestep=min(mintimestep, length(ic.results.TransientSolution));
531 end
532 if ~noearth,
533 mintimestep=min(mintimestep, length(self.earth.results.TransientSolution));
534 end
535
536 for i=1:length(self.icecaps),
[26358]537 ic=self.icecaps{i};
[22955]538 ic.results.TransientSolution=ic.results.TransientSolution(1:mintimestep);
539 self.icecaps{i}=ic;
540 end
541 ic=self.earth;
542 if ~noearth,
543 ic.results.TransientSolution=ic.results.TransientSolution(1:mintimestep);
544 end
545 self.earth=ic;
546 end % }}}
[24886]547 function self=initializemodels(self) % {{{
548
549 for i=1:length(self.basins),
550 md=model();
551 md.miscellaneous.name=self.basins{i}.name;
552 self.addicecap(md);
553 end
554 end % }}}
[20110]555 end
556end
Note: See TracBrowser for help on using the repository browser.