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

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

CHG: new dsl class to take care of the dynamic sea level in our sealevelrise_core. We made a new field
to the model class to fit the dsl fields. The fields come from the slr steric rates class essentially.
We also added GetStericRate and GetDynamicRate to the sea level core, which splits what used to be only
the steric rate field.

File size: 13.3 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), return; end
101 md = checkfield(md,'fieldname','slr.deltathickness','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]);
102 md = checkfield(md,'fieldname','slr.sealevel','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
103 md = checkfield(md,'fieldname','slr.spcthickness','Inf',1,'timeseries',1);
104 md = checkfield(md,'fieldname','slr.love_h','NaN',1,'Inf',1);
105 md = checkfield(md,'fieldname','slr.love_k','NaN',1,'Inf',1);
106 md = checkfield(md,'fieldname','slr.love_l','NaN',1,'Inf',1);
107 md = checkfield(md,'fieldname','slr.tide_love_h','NaN',1,'Inf',1);
108 md = checkfield(md,'fieldname','slr.tide_love_k','NaN',1,'Inf',1);
109 md = checkfield(md,'fieldname','slr.fluid_love','NaN',1,'Inf',1);
110 md = checkfield(md,'fieldname','slr.equatorial_moi','NaN',1,'Inf',1);
111 md = checkfield(md,'fieldname','slr.polar_moi','NaN',1,'Inf',1);
112 md = checkfield(md,'fieldname','slr.angular_velocity','NaN',1,'Inf',1);
113 md = checkfield(md,'fieldname','slr.reltol','size',[1 1]);
114 md = checkfield(md,'fieldname','slr.abstol','size',[1 1]);
115 md = checkfield(md,'fieldname','slr.maxiter','size',[1 1],'>=',1);
116 md = checkfield(md,'fieldname','slr.geodetic_run_frequency','size',[1 1],'>=',1);
117 md = checkfield(md,'fieldname','slr.hydro_rate','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
118 md = checkfield(md,'fieldname','slr.degacc','size',[1 1],'>=',1e-10);
119 md = checkfield(md,'fieldname','slr.requested_outputs','stringrow',1);
120 md = checkfield(md,'fieldname','slr.loop_increment','NaN',1,'Inf',1,'>=',1);
121 md = checkfield(md,'fieldname','slr.horiz','NaN',1,'Inf',1,'values',[0 1]);
122 md = checkfield(md,'fieldname','slr.Ngia','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
123 md = checkfield(md,'fieldname','slr.Ugia','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
124
125 %check that love numbers are provided at the same level of accuracy:
126 if (size(self.love_h,1)~=size(self.love_k,1) | size(self.love_h,1)~=size(self.love_l,1)),
127 error('slr error message: love numbers should be provided at the same level of accuracy');
128 end
129
130 %cross check that whereever we have an ice load, the mask is <0 on each vertex:
131 pos=find(self.deltathickness);
132 maskpos=md.mask.ice_levelset(md.mesh.elements(pos,:));
133 [els,vertices]=find(maskpos>0);
134 if length(els),
135 warning('slr checkconsistency fail: there are elements with ice loads where some vertices are not on the ice!');
136 end
137
138 %check that if geodetic is requested, we are a mesh3dsurface model (planet), or if we are not,
139 %a coupler to a planet model is provided.
140 if self.geodetic,
141 if md.transient.iscoupler,
142 %we are good;
143 else
144 if strcmpi(class(md.mesh),'mesh3dsurface'),
145 %we are good
146 else
147 error('model is requesting geodetic computations without being a mesh3dsurface, or being coupled to one!');
148 end
149 end
150 end
151
152
153 end % }}}
154 function list=defaultoutputs(self,md) % {{{
155 list = {'Sealevel'};
156 end % }}}
157 function disp(self) % {{{
158 disp(sprintf(' slr parameters:'));
159
160 fielddisplay(self,'deltathickness','thickness change: ice height equivalent [m]');
161 fielddisplay(self,'sealevel','current sea level (prior to computation) [m]');
162 fielddisplay(self,'spcthickness','thickness constraints (NaN means no constraint) [m]');
163 fielddisplay(self,'reltol','sea level rise relative convergence criterion, (default, NaN: not applied)');
164 fielddisplay(self,'abstol','sea level rise absolute convergence criterion, NaN: not applied');
165 fielddisplay(self,'maxiter','maximum number of nonlinear iterations');
166 fielddisplay(self,'love_h','load Love number for radial displacement');
167 fielddisplay(self,'love_k','load Love number for gravitational potential perturbation');
168 fielddisplay(self,'love_l','load Love number for horizontal displacements');
169 fielddisplay(self,'tide_love_k','tidal load Love number (deg 2)');
170 fielddisplay(self,'tide_love_h','tidal load Love number (deg 2)');
171 fielddisplay(self,'fluid_love','secular fluid Love number');
172 fielddisplay(self,'equatorial_moi','mean equatorial moment of inertia [kg m^2]');
173 fielddisplay(self,'polar_moi','polar moment of inertia [kg m^2]');
174 fielddisplay(self,'angular_velocity','mean rotational velocity of earth [per second]');
175 fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area [default: No correction]');
176 fielddisplay(self,'hydro_rate','rate of hydrological expansion (in mm/yr)');
177 fielddisplay(self,'Ngia','rate of viscous (GIA) geoid expansion (in mm/yr)');
178 fielddisplay(self,'Ugia','rate of viscous (GIA) bedrock uplift (in mm/yr)');
179 fielddisplay(self,'loop_increment','vector assembly (in the convolution) framentation');
180 fielddisplay(self,'geodetic','compute geodetic SLR? ( in addition to steric?) default 0');
181 fielddisplay(self,'geodetic_run_frequency','how many time steps we skip before we run SLR solver during transient (default: 1)');
182 fielddisplay(self,'rigid','rigid earth graviational potential perturbation');
183 fielddisplay(self,'elastic','elastic earth graviational potential perturbation');
184 fielddisplay(self,'rotation','earth rotational potential perturbation');
185 fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions');
186 fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps');
187 fielddisplay(self,'requested_outputs','additional outputs requested');
188
189 end % }}}
190 function marshall(self,prefix,md,fid) % {{{
191 %WriteData(fid,prefix,'object',self,'fieldname','deltathickness','format','DoubleMat','mattype',2);
192 WriteData(fid,prefix,'object',self,'fieldname','deltathickness','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
193 %WriteData(fid,prefix,'object',self,'fieldname','deltathickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofelements+1);
194 WriteData(fid,prefix,'object',self,'fieldname','sealevel','mattype',1,'format','DoubleMat','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
195 WriteData(fid,prefix,'object',self,'fieldname','spcthickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
196 WriteData(fid,prefix,'object',self,'fieldname','reltol','format','Double');
197 WriteData(fid,prefix,'object',self,'fieldname','abstol','format','Double');
198 WriteData(fid,prefix,'object',self,'fieldname','maxiter','format','Integer');
199 WriteData(fid,prefix,'object',self,'fieldname','love_h','format','DoubleMat','mattype',1);
200 WriteData(fid,prefix,'object',self,'fieldname','love_k','format','DoubleMat','mattype',1);
201 WriteData(fid,prefix,'object',self,'fieldname','love_l','format','DoubleMat','mattype',1);
202 WriteData(fid,prefix,'object',self,'fieldname','tide_love_h','format','Double');
203 WriteData(fid,prefix,'object',self,'fieldname','tide_love_k','format','Double');
204 WriteData(fid,prefix,'object',self,'fieldname','fluid_love','format','Double');
205 WriteData(fid,prefix,'object',self,'fieldname','equatorial_moi','format','Double');
206 WriteData(fid,prefix,'object',self,'fieldname','polar_moi','format','Double');
207 WriteData(fid,prefix,'object',self,'fieldname','angular_velocity','format','Double');
208 WriteData(fid,prefix,'object',self,'fieldname','rigid','format','Boolean');
209 WriteData(fid,prefix,'object',self,'fieldname','elastic','format','Boolean');
210 WriteData(fid,prefix,'object',self,'fieldname','rotation','format','Boolean');
211 WriteData(fid,prefix,'object',self,'fieldname','ocean_area_scaling','format','Boolean');
212 WriteData(fid,prefix,'object',self,'fieldname','geodetic_run_frequency','format','Integer');
213 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);
214 WriteData(fid,prefix,'object',self,'fieldname','Ngia','format','DoubleMat','mattype',1,'scale',1e-3/md.constants.yts);
215 WriteData(fid,prefix,'object',self,'fieldname','Ugia','format','DoubleMat','mattype',1,'scale',1e-3/md.constants.yts);
216 WriteData(fid,prefix,'object',self,'fieldname','degacc','format','Double');
217 WriteData(fid,prefix,'object',self,'fieldname','transitions','format','MatArray');
218 WriteData(fid,prefix,'object',self,'fieldname','loop_increment','format','Integer');
219 WriteData(fid,prefix,'object',self,'fieldname','horiz','format','Integer');
220 WriteData(fid,prefix,'object',self,'fieldname','geodetic','format','Integer');
221
222 %process requested outputs
223 outputs = self.requested_outputs;
224 pos = find(ismember(outputs,'default'));
225 if ~isempty(pos),
226 outputs(pos) = []; %remove 'default' from outputs
227 outputs = [outputs defaultoutputs(self,md)]; %add defaults
228 end
229 WriteData(fid,prefix,'data',outputs,'name','md.slr.requested_outputs','format','StringArray');
230
231 end % }}}
232 function savemodeljs(self,fid,modelname) % {{{
233
234 writejs1Darray(fid,[modelname '.slr.deltathickness'],self.deltathickness);
235 writejs1Darray(fid,[modelname '.slr.sealevel'],self.sealevel);
236 writejs1Darray(fid,[modelname '.slr.spcthickness'],self.spcthickness);
237 writejsdouble(fid,[modelname '.slr.maxiter'],self.maxiter);
238 writejsdouble(fid,[modelname '.slr.reltol'],self.reltol);
239 writejsdouble(fid,[modelname '.slr.abstol'],self.abstol);
240 writejs1Darray(fid,[modelname '.slr.love_h'],self.love_h);
241 writejs1Darray(fid,[modelname '.slr.love_k'],self.love_k);
242 writejs1Darray(fid,[modelname '.slr.love_l'],self.love_l);
243 writejsdouble(fid,[modelname '.slr.tide_love_k'],self.tide_love_k);
244 writejsdouble(fid,[modelname '.slr.tide_love_h'],self.tide_love_h);
245 writejsdouble(fid,[modelname '.slr.fluid_love'],self.fluid_love);
246 writejsdouble(fid,[modelname '.slr.equatorial_moi'],self.equatorial_moi);
247 writejsdouble(fid,[modelname '.slr.polar_moi'],self.polar_moi);
248 writejsdouble(fid,[modelname '.slr.angular_velocity'],self.angular_velocity);
249 writejsdouble(fid,[modelname '.slr.rigid'],self.rigid);
250 writejsdouble(fid,[modelname '.slr.elastic'],self.elastic);
251 writejsdouble(fid,[modelname '.slr.rotation'],self.rotation);
252 writejsdouble(fid,[modelname '.slr.ocean_area_scaling'],self.ocean_area_scaling);
253 writejsdouble(fid,[modelname '.slr.geodetic_run_frequency'],self.geodetic_run_frequency);
254 writejs1Darray(fid,[modelname '.slr.hydro_rate'],self.hydro_rate);
255 writejsdouble(fid,[modelname '.slr.degacc'],self.degacc);
256 writejscellstring(fid,[modelname '.slr.requested_outputs'],self.requested_outputs);
257 writejscellarray(fid,[modelname '.slr.transitions'],self.transitions);
258 end % }}}
259 function self = extrude(self,md) % {{{
260 self.sealevel=project3d(md,'vector',self.sealevel,'type','node');
261 end % }}}
262 end
263end
Note: See TracBrowser for help on using the repository browser.