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

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

CHG: fixed issue in sealevelmodel that resulted in wrong concatenation of icecaps and wrong intersections. Also fixed an issue in the solidearthsolution where horizontals were being computed without need for it.
File size: 17.7 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 % }}}
369 function basinsplot3d(self,varargin) % {{{
[25136]370 for i=1:length(self.basins),
371 self.basins{i}.plot3d(varargin{:});
372 end
[24156]373 end % }}}
374 function caticecaps(self,varargin) % {{{
375
[26358]376 %recover options:
377 options=pairoptions(varargin{:});
[24156]378 tolerance=getfieldvalue(options,'tolerance',.65);
379 loneedgesdetect=getfieldvalue(options,'loneedgesdetect',0);
380
381 %make 3D model:
382 models=self.icecaps;
383 for i=1:length(models),
[24779]384 models{i}=TwoDToThreeD(models{i},self.planet);
[24156]385 end
386
387 %Plug all models together:
[26358]388 md=models{1};
[24156]389 for i=2:length(models),
390 md=modelmerge3d(md,models{i},'tolerance',tolerance);
391 md.private.bamg.landmask=[md.private.bamg.landmask;models{i}.private.bamg.landmask];
392 end
393
394 %Look for lone edges if asked for it: {{{
395 if loneedgesdetect,
396 edges=loneedges(md);
397 plotmodel(md,'data',md.mask.land_levelset);
398 hold on;
399 for i=1:length(edges),
400 ind1=edges(i,1);
401 ind2=edges(i,2);
402 %plot([md.mesh.x(ind1),md.mesh.x(ind2)],[md.mesh.y(ind1),md.mesh.y(ind2)],'r*-');
403 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*-');
404 end
405 end %}}}
406
[26358]407 %Plug into earth:
[24156]408 self.earth=md;
409
[24886]410 %Create mesh radius:
411 self.earth.mesh.r=planetradius('earth')*ones(md.mesh.numberofvertices,1);
[27058]412 self.earth.mesh.lat = asind(self.earth.mesh.z./self.earth.mesh.r);
413 self.earth.mesh.long = atan2d(self.earth.mesh.y,self.earth.mesh.x);
[24886]414
[27058]415
[24156]416 end % }}}
[25620]417 function caticecaps2d(self,varargin) % {{{
418
[26358]419 %recover options:
420 options=pairoptions(varargin{:});
[25620]421 tolerance=getfieldvalue(options,'tolerance',1e-5);
422 loneedgesdetect=getfieldvalue(options,'loneedgesdetect',0);
423 models=self.icecaps;
424
425 %Plug all models together:
[26358]426 md=models{1};
[25620]427 for i=2:length(models),
428 md=modelmerge2d(md,models{i},'tolerance',tolerance);
429 end
430
431 %Look for lone edges if asked for it: {{{
432 if loneedgesdetect,
433 edges=loneedges(md);
434 hold on;
435 for i=1:length(edges),
436 ind1=edges(i,1);
437 ind2=edges(i,2);
438 plot([md.mesh.x(ind1),md.mesh.x(ind2)],[md.mesh.y(ind1),md.mesh.y(ind2)],'g*-');
[22955]439 end
[25620]440 end %}}}
441
[26358]442 %Plug into earth:
[25620]443 self.earth=md;
444
[22955]445 end % }}}
[25688]446 function viscousiterations(self) % {{{
447 for i=1:length(self.icecaps),
448 ic=self.icecaps{i};
449 mvi=ic.results.TransientSolution(1).StressbalanceConvergenceNumSteps;
450 for j=2:length(ic.results.TransientSolution)-1,
451 mvi=max(mvi,ic.results.TransientSolution(j).StressbalanceConvergenceNumSteps);
452 end
453 disp(sprintf('%i, %s: %i',i,self.icecaps{i}.miscellaneous.name,mvi));
454 end
455 end % }}}
[22955]456 function maxtimestep(self) % {{{
457 for i=1:length(self.icecaps),
[25688]458 ic=self.icecaps{i};
[22955]459 mvi=length(ic.results.TransientSolution);
460 timei=ic.results.TransientSolution(end).time;
461 disp(sprintf('%i, %s: %i/%g',i,self.icecaps{i}.miscellaneous.name,mvi,timei));
462 end
463 mvi=length(self.earth.results.TransientSolution);
464 timei=self.earth.results.TransientSolution(end).time;
465 disp(sprintf('Earth: %i/%g',mvi,timei));
466 end % }}}
[24165]467 function transfer(self,string) % {{{
[26358]468 %Recover field size in one icecap:
[24779]469 eval(['n=size(self.icecaps{1}.' string ',1);']);
[24165]470 if n==self.icecaps{1}.mesh.numberofvertices,
471 eval(['self.earth.' string '=zeros(self.earth.mesh.numberofvertices,1);']);
472 for i=1:length(self.icecaps),
473 eval(['self.earth.' string '(self.transitions{' num2str(i) '})=self.icecaps{' num2str(i) '}.' string ';']);
474 end
[24779]475 elseif n==(self.icecaps{1}.mesh.numberofvertices+1),
[26358]476 %dealing with a transient dataset.
[24779]477 %check that all timetags are similar between all icecaps: %{{{
478 for i=1:length(self.icecaps),
479 eval(['capfieldi= self.icecaps{' num2str(i) '}.' string ';']);
480 for j=(i+1):length(self.icecaps),
481 eval(['capfieldj= self.icecaps{' num2str(j) '}.' string ';']);
482 if ~isequal(capfieldi(end,:),capfieldj(end,:)),
[25499]483 error(['Time stamps for ' string ' field are different between icecaps ' num2str(i) ' and ' num2str(j)]);
[24779]484 end
485 end
486 end
[26358]487 eval(['capfield1= self.icecaps{1}.' string ';']);
[24779]488 times=capfield1(end,:);
489 nsteps=length(times);
490 %}}}
491 %initialize: %{{{
492 eval(['field=zeros(self.earth.mesh.numberofvertices+1,' num2str(nsteps) ');']);
493 eval(['field(end,:)=times;']); %transfer the times only, not the values
494 %}}}
495 %transfer all time fields: {{{
496 for i=1:length(self.icecaps),
497 eval(['capfieldi= self.icecaps{' num2str(i) '}.' string ';']);
498 for j=1:nsteps,
499 eval(['field(self.transitions{' num2str(i) '},' num2str(j) ')=capfieldi(1:end-1,' num2str(j) ');']); %transfer only the values, not the time.
500 end
501 end
502 eval(['self.earth.' string '=field;']); %do not forget to plug the field variable into its final location
503 %}}}
[25455]504 elseif n==self.icecaps{1}.mesh.numberofelements,
[24165]505 eval(['self.earth.' string '=zeros(self.earth.mesh.numberofelements,1);']);
506 for i=1:length(self.icecaps),
507 eval(['self.earth.' string '(self.eltransitions{' num2str(i) '})=self.icecaps{' num2str(i) '}.' string ';']);
508 end
509 else
510 error('not supported yet');
511 end
512 end % }}}
[22955]513 function self=homogeneize(self,noearth) % {{{
514 if nargin==1,
515 noearth=0;
516 end
517 mintimestep=Inf;
518 for i=1:length(self.icecaps),
[26358]519 ic=self.icecaps{i};
[22955]520 mintimestep=min(mintimestep, length(ic.results.TransientSolution));
521 end
522 if ~noearth,
523 mintimestep=min(mintimestep, length(self.earth.results.TransientSolution));
524 end
525
526 for i=1:length(self.icecaps),
[26358]527 ic=self.icecaps{i};
[22955]528 ic.results.TransientSolution=ic.results.TransientSolution(1:mintimestep);
529 self.icecaps{i}=ic;
530 end
531 ic=self.earth;
532 if ~noearth,
533 ic.results.TransientSolution=ic.results.TransientSolution(1:mintimestep);
534 end
535 self.earth=ic;
536 end % }}}
[24886]537 function self=initializemodels(self) % {{{
538
539 for i=1:length(self.basins),
540 md=model();
541 md.miscellaneous.name=self.basins{i}.name;
542 self.addicecap(md);
543 end
544 end % }}}
[20110]545 end
546end
Note: See TracBrowser for help on using the repository browser.