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

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

CHG: added a check on the type of mask that should be initialized for an slr solution.
File size: 13.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 spcthickness = NaN;
11 maxiter = 0;
12 reltol = 0;
13 abstol = 0;
14 love_h = 0; %provided by PREM model
15 love_k = 0; %ideam
16 love_l = 0; %ideam
17 tide_love_k = 0; %ideam
18 tide_love_h = 0; %ideam
19 fluid_love = 0;
20 equatorial_moi = 0;
21 polar_moi = 0;
22 angular_velocity = 0;
23 rigid = 0;
24 elastic = 0;
25 rotation = 0;
26 ocean_area_scaling = 0;
27 hydro_rate = 0; %rate of steric expansion from hydrological effects.
28 geodetic_run_frequency = 1; %how many time steps we skip before we run the geodetic part of the solver during transient
29 geodetic = 0; %compute geodetic SLR? (in addition to steric?)
30 degacc = 0;
31 loop_increment = 0;
32 horiz = 0;
33 Ngia = NaN;
34 Ugia = NaN;
35
36 requested_outputs = {};
37 transitions = {};
38 end
39 methods
40 function self = slr(varargin) % {{{
41 switch nargin
42 case 0
43 self=setdefaultparameters(self);
44 otherwise
45 error('constructor not supported');
46 end
47 end % }}}
48 function self = setdefaultparameters(self) % {{{
49
50 %Convergence criterion: absolute, relative and residual
51 self.reltol=0.01; % 1 per cent
52 self.abstol=NaN; % default
53
54 %maximum of non-linear iterations.
55 self.maxiter=5;
56 self.loop_increment=200;
57
58 %computational flags:
59 self.geodetic=0;
60 self.rigid=1;
61 self.elastic=1;
62 self.ocean_area_scaling=0;
63 self.rotation=1;
64
65 %tidal love numbers:
66 self.tide_love_h=0.6149; %degree 2
67 self.tide_love_k=0.3055; % degree 2
68
69 %secular fluid love number:
70 self.fluid_love=0.942;
71
72 %moment of inertia:
73 self.equatorial_moi=8.0077*10^37; % [kg m^2]
74 self.polar_moi =8.0345*10^37; % [kg m^2]
75
76 % mean rotational velocity of earth
77 self.angular_velocity=7.2921*10^-5; % [s^-1]
78
79 %numerical discretization accuracy
80 self.degacc=.01;
81
82 %hydro
83 self.hydro_rate=0;
84
85 %how many time steps we skip before we run SLR solver during transient
86 self.geodetic_run_frequency=1;
87
88 %output default:
89 self.requested_outputs={'default'};
90
91 %transitions should be a cell array of vectors:
92 self.transitions={};
93
94 %horizontal displacement? (not by default)
95 self.horiz=0;
96
97 end % }}}
98 function md = checkconsistency(self,md,solution,analyses) % {{{
99
100 if ~ismember('SealevelriseAnalysis',analyses) | (strcmp(solution,'TransientSolution') & md.transient.isslr==0),
101 return;
102 end
103
104 %check that the mask is a maskpsl class type, otherwise, we won't be able to run very far
105 %with the slr solution:
106 if ~strcmpi(class(md.mask),'maskpsl'),
107 error('slr error message: model mask should be of the ''maskpsl'' class');
108 end
109
110 md = checkfield(md,'fieldname','slr.deltathickness','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]);
111 md = checkfield(md,'fieldname','slr.sealevel','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
112 md = checkfield(md,'fieldname','slr.spcthickness','Inf',1,'timeseries',1);
113 md = checkfield(md,'fieldname','slr.love_h','NaN',1,'Inf',1);
114 md = checkfield(md,'fieldname','slr.love_k','NaN',1,'Inf',1);
115 md = checkfield(md,'fieldname','slr.love_l','NaN',1,'Inf',1);
116 md = checkfield(md,'fieldname','slr.tide_love_h','NaN',1,'Inf',1);
117 md = checkfield(md,'fieldname','slr.tide_love_k','NaN',1,'Inf',1);
118 md = checkfield(md,'fieldname','slr.fluid_love','NaN',1,'Inf',1);
119 md = checkfield(md,'fieldname','slr.equatorial_moi','NaN',1,'Inf',1);
120 md = checkfield(md,'fieldname','slr.polar_moi','NaN',1,'Inf',1);
121 md = checkfield(md,'fieldname','slr.angular_velocity','NaN',1,'Inf',1);
122 md = checkfield(md,'fieldname','slr.reltol','size',[1 1]);
123 md = checkfield(md,'fieldname','slr.abstol','size',[1 1]);
124 md = checkfield(md,'fieldname','slr.maxiter','size',[1 1],'>=',1);
125 md = checkfield(md,'fieldname','slr.geodetic_run_frequency','size',[1 1],'>=',1);
126 md = checkfield(md,'fieldname','slr.hydro_rate','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
127 md = checkfield(md,'fieldname','slr.degacc','size',[1 1],'>=',1e-10);
128 md = checkfield(md,'fieldname','slr.requested_outputs','stringrow',1);
129 md = checkfield(md,'fieldname','slr.loop_increment','NaN',1,'Inf',1,'>=',1);
130 md = checkfield(md,'fieldname','slr.horiz','NaN',1,'Inf',1,'values',[0 1]);
131 md = checkfield(md,'fieldname','slr.Ngia','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
132 md = checkfield(md,'fieldname','slr.Ugia','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
133
134 %check that love numbers are provided at the same level of accuracy:
135 if (size(self.love_h,1)~=size(self.love_k,1) | size(self.love_h,1)~=size(self.love_l,1)),
136 error('slr error message: love numbers should be provided at the same level of accuracy');
137 end
138
139 %cross check that whereever we have an ice load, the mask is <0 on each vertex:
140 pos=find(self.deltathickness);
141 maskpos=md.mask.ice_levelset(md.mesh.elements(pos,:));
142 [els,vertices]=find(maskpos>0);
143 if length(els),
144 warning('slr checkconsistency fail: there are elements with ice loads where some vertices are not on the ice!');
145 end
146
147 %check that if geodetic is requested, we are a mesh3dsurface model (planet), or if we are not,
148 %a coupler to a planet model is provided.
149 if self.geodetic,
150 if md.transient.iscoupler,
151 %we are good;
152 else
153 if strcmpi(class(md.mesh),'mesh3dsurface'),
154 %we are good
155 else
156 error('model is requesting geodetic computations without being a mesh3dsurface, or being coupled to one!');
157 end
158 end
159 end
160
161 end % }}}
162 function list=defaultoutputs(self,md) % {{{
163 list = {'Sealevel'};
164 end % }}}
165 function disp(self) % {{{
166 disp(sprintf(' slr parameters:'));
167
168 fielddisplay(self,'deltathickness','thickness change: ice height equivalent [m]');
169 fielddisplay(self,'sealevel','current sea level (prior to computation) [m]');
170 fielddisplay(self,'spcthickness','thickness constraints (NaN means no constraint) [m]');
171 fielddisplay(self,'reltol','sea level rise relative convergence criterion, (default, NaN: not applied)');
172 fielddisplay(self,'abstol','sea level rise absolute convergence criterion, NaN: not applied');
173 fielddisplay(self,'maxiter','maximum number of nonlinear iterations');
174 fielddisplay(self,'love_h','load Love number for radial displacement');
175 fielddisplay(self,'love_k','load Love number for gravitational potential perturbation');
176 fielddisplay(self,'love_l','load Love number for horizontal displacements');
177 fielddisplay(self,'tide_love_k','tidal load Love number (deg 2)');
178 fielddisplay(self,'tide_love_h','tidal load Love number (deg 2)');
179 fielddisplay(self,'fluid_love','secular fluid Love number');
180 fielddisplay(self,'equatorial_moi','mean equatorial moment of inertia [kg m^2]');
181 fielddisplay(self,'polar_moi','polar moment of inertia [kg m^2]');
182 fielddisplay(self,'angular_velocity','mean rotational velocity of earth [per second]');
183 fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area [default: No correction]');
184 fielddisplay(self,'hydro_rate','rate of hydrological expansion (in mm/yr)');
185 fielddisplay(self,'Ngia','rate of viscous (GIA) geoid expansion (in mm/yr)');
186 fielddisplay(self,'Ugia','rate of viscous (GIA) bedrock uplift (in mm/yr)');
187 fielddisplay(self,'loop_increment','vector assembly (in the convolution) framentation');
188 fielddisplay(self,'geodetic','compute geodetic SLR? ( in addition to steric?) default 0');
189 fielddisplay(self,'geodetic_run_frequency','how many time steps we skip before we run SLR solver during transient (default: 1)');
190 fielddisplay(self,'rigid','rigid earth graviational potential perturbation');
191 fielddisplay(self,'elastic','elastic earth graviational potential perturbation');
192 fielddisplay(self,'rotation','earth rotational potential perturbation');
193 fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions');
194 fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps');
195 fielddisplay(self,'requested_outputs','additional outputs requested');
196
197 end % }}}
198 function marshall(self,prefix,md,fid) % {{{
199 %WriteData(fid,prefix,'object',self,'fieldname','deltathickness','format','DoubleMat','mattype',2);
200 WriteData(fid,prefix,'object',self,'fieldname','deltathickness','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
201 %WriteData(fid,prefix,'object',self,'fieldname','deltathickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofelements+1);
202 WriteData(fid,prefix,'object',self,'fieldname','sealevel','mattype',1,'format','DoubleMat','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
203 WriteData(fid,prefix,'object',self,'fieldname','spcthickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
204 WriteData(fid,prefix,'object',self,'fieldname','reltol','format','Double');
205 WriteData(fid,prefix,'object',self,'fieldname','abstol','format','Double');
206 WriteData(fid,prefix,'object',self,'fieldname','maxiter','format','Integer');
207 WriteData(fid,prefix,'object',self,'fieldname','love_h','format','DoubleMat','mattype',1);
208 WriteData(fid,prefix,'object',self,'fieldname','love_k','format','DoubleMat','mattype',1);
209 WriteData(fid,prefix,'object',self,'fieldname','love_l','format','DoubleMat','mattype',1);
210 WriteData(fid,prefix,'object',self,'fieldname','tide_love_h','format','Double');
211 WriteData(fid,prefix,'object',self,'fieldname','tide_love_k','format','Double');
212 WriteData(fid,prefix,'object',self,'fieldname','fluid_love','format','Double');
213 WriteData(fid,prefix,'object',self,'fieldname','equatorial_moi','format','Double');
214 WriteData(fid,prefix,'object',self,'fieldname','polar_moi','format','Double');
215 WriteData(fid,prefix,'object',self,'fieldname','angular_velocity','format','Double');
216 WriteData(fid,prefix,'object',self,'fieldname','rigid','format','Boolean');
217 WriteData(fid,prefix,'object',self,'fieldname','elastic','format','Boolean');
218 WriteData(fid,prefix,'object',self,'fieldname','rotation','format','Boolean');
219 WriteData(fid,prefix,'object',self,'fieldname','ocean_area_scaling','format','Boolean');
220 WriteData(fid,prefix,'object',self,'fieldname','geodetic_run_frequency','format','Integer');
221 WriteData(fid,prefix,'object',self,'fieldname','hydro_rate','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1e-3/md.constants.yts);
222 WriteData(fid,prefix,'object',self,'fieldname','Ngia','format','DoubleMat','mattype',1,'scale',1e-3/md.constants.yts);
223 WriteData(fid,prefix,'object',self,'fieldname','Ugia','format','DoubleMat','mattype',1,'scale',1e-3/md.constants.yts);
224 WriteData(fid,prefix,'object',self,'fieldname','degacc','format','Double');
225 WriteData(fid,prefix,'object',self,'fieldname','transitions','format','MatArray');
226 WriteData(fid,prefix,'object',self,'fieldname','loop_increment','format','Integer');
227 WriteData(fid,prefix,'object',self,'fieldname','horiz','format','Integer');
228 WriteData(fid,prefix,'object',self,'fieldname','geodetic','format','Integer');
229
230 %process requested outputs
231 outputs = self.requested_outputs;
232 pos = find(ismember(outputs,'default'));
233 if ~isempty(pos),
234 outputs(pos) = []; %remove 'default' from outputs
235 outputs = [outputs defaultoutputs(self,md)]; %add defaults
236 end
237 WriteData(fid,prefix,'data',outputs,'name','md.slr.requested_outputs','format','StringArray');
238
239 end % }}}
240 function savemodeljs(self,fid,modelname) % {{{
241
242 writejs1Darray(fid,[modelname '.slr.deltathickness'],self.deltathickness);
243 writejs1Darray(fid,[modelname '.slr.sealevel'],self.sealevel);
244 writejs1Darray(fid,[modelname '.slr.spcthickness'],self.spcthickness);
245 writejsdouble(fid,[modelname '.slr.maxiter'],self.maxiter);
246 writejsdouble(fid,[modelname '.slr.reltol'],self.reltol);
247 writejsdouble(fid,[modelname '.slr.abstol'],self.abstol);
248 writejs1Darray(fid,[modelname '.slr.love_h'],self.love_h);
249 writejs1Darray(fid,[modelname '.slr.love_k'],self.love_k);
250 writejs1Darray(fid,[modelname '.slr.love_l'],self.love_l);
251 writejsdouble(fid,[modelname '.slr.tide_love_k'],self.tide_love_k);
252 writejsdouble(fid,[modelname '.slr.tide_love_h'],self.tide_love_h);
253 writejsdouble(fid,[modelname '.slr.fluid_love'],self.fluid_love);
254 writejsdouble(fid,[modelname '.slr.equatorial_moi'],self.equatorial_moi);
255 writejsdouble(fid,[modelname '.slr.polar_moi'],self.polar_moi);
256 writejsdouble(fid,[modelname '.slr.angular_velocity'],self.angular_velocity);
257 writejsdouble(fid,[modelname '.slr.rigid'],self.rigid);
258 writejsdouble(fid,[modelname '.slr.elastic'],self.elastic);
259 writejsdouble(fid,[modelname '.slr.rotation'],self.rotation);
260 writejsdouble(fid,[modelname '.slr.ocean_area_scaling'],self.ocean_area_scaling);
261 writejsdouble(fid,[modelname '.slr.geodetic_run_frequency'],self.geodetic_run_frequency);
262 writejs1Darray(fid,[modelname '.slr.hydro_rate'],self.hydro_rate);
263 writejsdouble(fid,[modelname '.slr.degacc'],self.degacc);
264 writejscellstring(fid,[modelname '.slr.requested_outputs'],self.requested_outputs);
265 writejscellarray(fid,[modelname '.slr.transitions'],self.transitions);
266 end % }}}
267 function self = extrude(self,md) % {{{
268 self.sealevel=project3d(md,'vector',self.sealevel,'type','node');
269 end % }}}
270 end
271end
Note: See TracBrowser for help on using the repository browser.