source: issm/trunk-jpl/src/m/classes/slr.m@ 20322

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

CHG: slr: check against the mask to ensure conservation of mass.

model: additional checks.
File size: 6.6 KB
Line 
1%SLR class definition
2%
3% Usage:
4% slr=slr();
5
6classdef slr
7 properties (SetAccess=public)
8 deltathickness = NaN;
9 sealevel = NaN;
10 maxiter = 0;
11 reltol = 0;
12 abstol = 0;
13 love_h = 0; %provided by PREM model
14 love_k = 0; %ideam
15 rigid = 0;
16 elastic = 0;
17 eustatic = 0;
18 degacc = 0;
19 requested_outputs = {};
20 transitions = {};
21 end
22 methods
23 function self = slr(varargin) % {{{
24 switch nargin
25 case 0
26 self=setdefaultparameters(self);
27 otherwise
28 error('constructor not supported');
29 end
30 end % }}}
31 function self = setdefaultparameters(self) % {{{
32
33 %Convergence criterion: absolute, relative and residual
34 self.reltol=NaN; %default
35 self.abstol=0.001; %1 mm of sea level rise
36
37 %maximum of non-linear iterations.
38 self.maxiter=10;
39
40 %computational flags:
41 self.rigid=1;
42 self.elastic=1;
43 self.eustatic=1;
44
45 %numerical discretization accuracy
46 self.degacc=.01;
47
48 %output default:
49 self.requested_outputs={'default'};
50
51 %transitions should be a cell array of vectors:
52 self.transitions={};
53
54 end % }}}
55 function md = checkconsistency(self,md,solution,analyses) % {{{
56
57 if ~ismember(SealevelriseAnalysisEnum(),analyses), return; end
58 md = checkfield(md,'fieldname','slr.deltathickness','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]);
59 %md = checkfield(md,'fieldname','slr.deltathickness','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
60 md = checkfield(md,'fieldname','slr.sealevel','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
61 md = checkfield(md,'fieldname','slr.love_h','NaN',1,'Inf',1);
62 md = checkfield(md,'fieldname','slr.love_k','NaN',1,'Inf',1);
63 md = checkfield(md,'fieldname','slr.reltol','size',[1 1]);
64 md = checkfield(md,'fieldname','slr.abstol','size',[1 1]);
65 md = checkfield(md,'fieldname','slr.maxiter','size',[1 1],'>=',1);
66 md = checkfield(md,'fieldname','slr.degacc','size',[1 1],'>=',1e-10);
67 md = checkfield(md,'fieldname','slr.requested_outputs','stringrow',1);
68
69 %check that love numbers are provided at the same level of accuracy:
70 if (size(self.love_h,1) ~= size(self.love_k,1)),
71 error('slr error message: love numbers should be provided at the same level of accuracy');
72 end
73
74 %cross check that whereever we have an ice load, the mask is <0 on each vertex:
75 pos=find(self.deltathickness);
76 maskpos=md.mask.ice_levelset(md.mesh.elements(pos,:));
77 [els,vertices]=find(maskpos>0);
78 if length(els),
79 error('slr checkconsistency fail: there are elements with ice loads where some vertices are not on the ice!');
80 end
81
82 end % }}}
83 function list=defaultoutputs(self,md) % {{{
84 list = {'Sealevel'};
85 end % }}}
86 function disp(self) % {{{
87 disp(sprintf(' slr parameters:'));
88
89 fielddisplay(self,'deltathickness','thickness change (main loading of the slr solution core [m]');
90 fielddisplay(self,'sealevel','current sea level (prior to computation) [m]');
91 fielddisplay(self,'reltol','sea level rise relative convergence criterion, (default, NaN: not applied)');
92 fielddisplay(self,'abstol','sea level rise absolute convergence criterion, NaN: not applied');
93 fielddisplay(self,'maxiter','maximum number of nonlinear iterations');
94 fielddisplay(self,'love_h','love load number for radial displacement');
95 fielddisplay(self,'love_k','love load number for gravitational potential perturbation');
96 fielddisplay(self,'rigid','rigid earth graviational potential perturbation');
97 fielddisplay(self,'elastic','elastic earth graviational potential perturbation');
98 fielddisplay(self,'eustatic','eustatic sea level rise');
99 fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions');
100 fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps');
101 fielddisplay(self,'requested_outputs','additional outputs requested');
102
103 end % }}}
104 function marshall(self,md,fid) % {{{
105 WriteData(fid,'object',self,'class','sealevelrise','fieldname','deltathickness','format','DoubleMat','mattype',2);
106 %WriteData(fid,'object',self,'class','sealevelrise','fieldname','deltathickness','format','DoubleMat','mattype',1);
107 WriteData(fid,'data',self.sealevel,'mattype',1,'format','DoubleMat','enum',SealevelEnum(),'timeserieslength',md.mesh.numberofvertices+1);
108 WriteData(fid,'object',self,'class','sealevelrise','fieldname','reltol','format','Double');
109 WriteData(fid,'object',self,'class','sealevelrise','fieldname','abstol','format','Double');
110 WriteData(fid,'object',self,'class','sealevelrise','fieldname','maxiter','format','Integer');
111 WriteData(fid,'object',self,'class','sealevelrise','fieldname','love_h','format','DoubleMat','mattype',1);
112 WriteData(fid,'object',self,'class','sealevelrise','fieldname','love_k','format','DoubleMat','mattype',1);
113 WriteData(fid,'object',self,'class','sealevelrise','fieldname','rigid','format','Boolean');
114 WriteData(fid,'object',self,'class','sealevelrise','fieldname','elastic','format','Boolean');
115 WriteData(fid,'object',self,'class','sealevelrise','fieldname','eustatic','format','Boolean');
116 WriteData(fid,'object',self,'class','sealevelrise','fieldname','degacc','format','Double');
117 WriteData(fid,'object',self,'class','sealevelrise','fieldname','transitions','format','MatArray');
118
119 %process requested outputs
120 outputs = self.requested_outputs;
121 pos = find(ismember(outputs,'default'));
122 if ~isempty(pos),
123 outputs(pos) = []; %remove 'default' from outputs
124 outputs = [outputs defaultoutputs(self,md)]; %add defaults
125 end
126 WriteData(fid,'data',outputs,'enum',SealevelriseRequestedOutputsEnum,'format','StringArray');
127
128 end % }}}
129 function savemodeljs(self,fid,modelname) % {{{
130
131 writejs1Darray(fid,[modelname '.slr.deltathickness'],self.deltathickness);
132 writejs1Darray(fid,[modelname '.slr.sealevel'],self.sealevel);
133 writejsdouble(fid,[modelname '.slr.maxiter'],self.maxiter);
134 writejsdouble(fid,[modelname '.slr.reltol'],self.reltol);
135 writejsdouble(fid,[modelname '.slr.abstol'],self.abstol);
136 writejs1Darray(fid,[modelname '.slr.love_h'],self.love_h);
137 writejs1Darray(fid,[modelname '.slr.love_k'],self.love_k);
138 writejsdouble(fid,[modelname '.slr.rigid'],self.rigid);
139 writejsdouble(fid,[modelname '.slr.eustatic'],self.eustatic);
140 writejsdouble(fid,[modelname '.slr.elastic'],self.elastic);
141 writejsdouble(fid,[modelname '.slr.degacc'],self.degacc);
142 writejscellstring(fid,[modelname '.slr.requested_outputs'],self.requested_outputs);
143 writejscellarray(fid,[modelname '.slr.transitions'],self.transitions);
144 end % }}}
145 end
146end
Note: See TracBrowser for help on using the repository browser.