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
Line 
1%SEALEVELMODEL class definition
2%
3% Usage:
4% slm = sealevelmodel(varargin)
5%
6% where varargin is a variable list of options:
7%
8% Example:
9% slm = sealevel('icecap',md_greenland,'icecap',md_antarctica,'earth',md_earth);
10
11classdef sealevelmodel < handle
12 properties (SetAccess=public) %Model fields
13 % {{{
14 icecaps = {}; % list of land/ice models, name should be change longer term.
15 earth = 0; % model for the whole earth
16 basins = {}; % list of basins, matching icecaps, where shapefile info is held
17 cluster = 0;
18 miscellaneous = 0;
19 settings = 0;
20 private = 0;
21 mergedcaps = 0;
22 transitions = {};
23 eltransitions = {};
24 planet = '';
25 %}}}
26 end
27 methods
28 function slm = sealevelmodel(varargin) % {{{
29 slm=setdefaultparameters(slm);
30
31 if nargin==1,
32
33 options=pairoptions(varargin{:});
34
35 %recover all the icecap models:
36 slm.icecaps=getfieldvalues(options,'ice_cap',{});
37
38 %recover the earth model:
39 slm.earth = getfieldvalue(options,'earth',0);
40
41 %set planet type:
42 slm.planet=getfieldvalue(options,'planet','earth');
43
44 end
45 end
46 %}}}
47 function checkconsistency(slm,solutiontype) % {{{
48
49 %is the coupler turned on?
50 for i=1:length(slm.icecaps),
51 if slm.icecaps{i}.transient.iscoupler==0,
52 warning(sprintf('sealevelmodel.m::checkconsistency: icecap model %s should have the transient coupler option turned on!',slm.icecaps{i}.miscellaneous.name));
53 end
54 end
55
56 if slm.earth.transient.iscoupler==0,
57 warning('sealevelmodel.m::checkconsistency: earth model should have the transient coupler option turned on!');
58 end
59
60 %check that the transition vectors have the right size:
61 for i=1:length(slm.icecaps),
62 if slm.icecaps{i}.mesh.numberofvertices ~= length(slm.earth.solidearth.transitions{i}),
63 error(['sealevelmodel.m::checkconsistency: issue with size of transition vector for ice cap: ' num2str(i) ' name: ' slm.icecaps{i}.miscellaneous.name]);
64 end
65 end
66
67 %check that runfrequency is the same everywhere:
68 for i=1:length(slm.icecaps),
69 if slm.icecaps{i}.solidearth.settings.runfrequency~=slm.earth.solidearth.settings.runfrequency,
70 error(sprintf('sealevelmodel.m::checkconsistency: icecap model %s should have the same run frequency as earth!',slm.icecaps{i}.miscellaneous.name));
71 end
72 end
73
74 %make sure steric_rate is the same everywhere:
75 for i=1:length(slm.icecaps),
76 md= slm.icecaps{i};
77 if ~isempty(find(md.dsl.steric_rate - slm.earth.dsl.steric_rate(slm.earth.dsl.transitions{i}))),
78 error(sprintf('sealevelmodel.m::checkconsistency: steric rate on ice cap %s is not the same as for the earth\n',md.miscellaneous.name));
79 end
80 end
81
82 %make sure grd is the same everywhere:
83 for i=1:length(slm.icecaps),
84 md= slm.icecaps{i};
85 if md.solidearthsettings.isgrd~=slm.earth.solidearthsettings.isgrd
86 error(sprintf('sealevelmodel.m::checkconsistency: isgrd on ice cap %s is not the same as for the earth\n',md.miscellaneous.name));
87 end
88 end
89
90 %make sure that there is no solid earth external forcing on the basins:
91 for i=1:length(slm.icecaps),
92 md= slm.icecaps{i};
93 if ~isempty(md.solidearth.external),
94 error('sealevelmodel.m::checkconsistency: cannot run external forcings on an ice sheet when running a coupling earth/ice sheet model');
95 end
96
97 end
98 %make sure that we have the right grd model for computing out sealevel patterns:
99 for i=1:length(slm.icecaps),
100 md= slm.icecaps{i};
101 if md.solidearth.settings.grdmodel~=0
102 error(sprintf('sealevelmodel.m::checkconsistency: ice sheets do not run GRD module, specify solidearth.settings.grdmodel=0 on ice cap %i',i));
103 end
104 end
105
106 end
107 %}}}
108 function slm = setdefaultparameters(slm) % {{{
109
110 %initialize subclasses
111 slm.icecaps = {};
112 slm.earth = {};
113 slm.miscellaneous = miscellaneous();
114 slm.settings = issmsettings();
115 slm.private = private();
116 slm.cluster = generic();
117 slm.transitions = {};
118 slm.eltransitions = {};
119 slm.planet = 'earth';
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 % }}}
129 function self=mergeresults(self) % {{{
130 champs=fieldnames(self.icecaps{1}.results.TransientSolution);
131 for i=1:length(self.mergedcaps)/2,
132 md=self.mergedcaps{2*(i-1)+1}; trans=self.mergedcaps{2*(i-1)+2};
133 %icecaps=self.icecaps(self.range{2*(i-1)+2});
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'),
137 %vertex or element?
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
144 else
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
160 end
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 % }}}
170 function n=ncaps(self) % {{{
171 n=length(self.icecaps);
172 end % }}}
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 % }}}
189 function addbasin(self,bas) % {{{
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;
194 end % }}}
195 function intersections2d(self,varargin) % {{{
196
197 options=pairoptions(varargin{:});
198 force=getfieldvalue(options,'force',0);
199
200 %initialize, to avoid issues of having more transitions than meshes.
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
211 %for elements:
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 % }}}
222 function intersections(self,varargin) % {{{
223
224 options=pairoptions(varargin{:});
225 force=getfieldvalue(options,'force',0);
226
227 %initialize, to avoid issues of having more transitions than meshes.
228 self.transitions={};
229 self.eltransitions={};
230
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;
235
236 for i=1:length(self.icecaps),
237 mdi=self.icecaps{i};
238 mdi=TwoDToThreeD(mdi,self.planet);
239
240
241 %for elements:
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
251 end % }}}
252 function checkintersections(self) % {{{
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');
258
259 end % }}}
260 function checkbasinconsistency(self) % {{{
261 for i=1:length(self.basins),
262 self.basins{i}.checkconsistency();
263 end
264
265 end % }}}
266 function baslist=basinindx(self,varargin) % {{{
267 options=pairoptions(varargin{:});
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,
274 if strcmpi(continent{1},'all'),
275 %need to transform this into a list of continents:
276 continent={};
277 for i=1:length(self.basins),
278 continent{end+1}=self.basins{i}.continent;
279 end
280 continent=unique(continent);
281 else
282 %nothing to do: assume we have a list of continents
283 end
284 else
285 %nothing to do: assume we have a list of continents
286 end
287 else
288 if strcmpi(continent,'all'),
289 %need to transform this into a list of continents:
290 continent={};
291 for i=1:length(self.basins),
292 continent{end+1}=self.basins{i}.continent;
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,
303 if strcmpi(bas{1},'all'),
304 %need to transform this into a list of basins:
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
322
323 end
324 else
325 %we have a list of basin names:
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) % {{{
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;
368 end % }}}
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 % }}}
379 function basinsplot3d(self,varargin) % {{{
380 for i=1:length(self.basins),
381 self.basins{i}.plot3d(varargin{:});
382 end
383 end % }}}
384 function caticecaps(self,varargin) % {{{
385
386 %recover options:
387 options=pairoptions(varargin{:});
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),
394 models{i}=TwoDToThreeD(models{i},self.planet);
395 end
396
397 %Plug all models together:
398 md=models{1};
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
417 %Plug into earth:
418 self.earth=md;
419
420 %Create mesh radius:
421 self.earth.mesh.r=planetradius('earth')*ones(md.mesh.numberofvertices,1);
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);
424
425
426 end % }}}
427 function caticecaps2d(self,varargin) % {{{
428
429 %recover options:
430 options=pairoptions(varargin{:});
431 tolerance=getfieldvalue(options,'tolerance',1e-5);
432 loneedgesdetect=getfieldvalue(options,'loneedgesdetect',0);
433 models=self.icecaps;
434
435 %Plug all models together:
436 md=models{1};
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*-');
449 end
450 end %}}}
451
452 %Plug into earth:
453 self.earth=md;
454
455 end % }}}
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 % }}}
466 function maxtimestep(self) % {{{
467 for i=1:length(self.icecaps),
468 ic=self.icecaps{i};
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 % }}}
477 function transfer(self,string) % {{{
478 %Recover field size in one icecap:
479 eval(['n=size(self.icecaps{1}.' string ',1);']);
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
485 elseif n==(self.icecaps{1}.mesh.numberofvertices+1),
486 %dealing with a transient dataset.
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,:)),
493 error(['Time stamps for ' string ' field are different between icecaps ' num2str(i) ' and ' num2str(j)]);
494 end
495 end
496 end
497 eval(['capfield1= self.icecaps{1}.' string ';']);
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 %}}}
514 elseif n==self.icecaps{1}.mesh.numberofelements,
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 % }}}
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),
529 ic=self.icecaps{i};
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),
537 ic=self.icecaps{i};
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 % }}}
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 % }}}
555 end
556end
Note: See TracBrowser for help on using the repository browser.