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

Last change on this file since 22808 was 22808, checked in by adhikari, 7 years ago

minor
File size: 9.8 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 love_l = 0; %ideam
16 tide_love_k = 0; %ideam
17 tide_love_h = 0; %ideam
18 fluid_love = 0;
19 equatorial_moi = 0;
20 polar_moi = 0;
21 angular_velocity = 0;
22 rigid = 0;
23 elastic = 0;
24 rotation = 0;
25 ocean_area_scaling = 0;
26 steric_rate = 0; %rate of ocean expansion from steric effects.
27 degacc = 0;
28 requested_outputs = {};
29 transitions = {};
30 end
31 methods
32 function self = slr(varargin) % {{{
33 switch nargin
34 case 0
35 self=setdefaultparameters(self);
36 otherwise
37 error('constructor not supported');
38 end
39 end % }}}
40 function self = setdefaultparameters(self) % {{{
41
42 %Convergence criterion: absolute, relative and residual
43 self.reltol=0.01; % 1 per cent
44 self.abstol=NaN; % default
45
46 %maximum of non-linear iterations.
47 self.maxiter=5;
48
49 %computational flags:
50 self.rigid=1;
51 self.elastic=1;
52 self.rotation=0;
53 self.ocean_area_scaling=0;
54
55 %tidal love numbers:
56 self.tide_love_h=0.6149; %degree 2
57 self.tide_love_k=0.3055; % degree 2
58
59 %secular fluid love number:
60 self.fluid_love=0.942;
61
62 %moment of inertia:
63 self.equatorial_moi=8.0077*10^37; % [kg m^2]
64 self.polar_moi =8.0345*10^37; % [kg m^2]
65
66 % mean rotational velocity of earth
67 self.angular_velocity=7.2921*10^-5; % [s^-1]
68
69 %numerical discretization accuracy
70 self.degacc=.01;
71
72 %steric:
73 self.steric_rate=0;
74
75
76 %output default:
77 self.requested_outputs={'default'};
78
79 %transitions should be a cell array of vectors:
80 self.transitions={};
81
82 end % }}}
83 function md = checkconsistency(self,md,solution,analyses) % {{{
84
85 if ~ismember('SealevelriseAnalysis',analyses), return; end
86 md = checkfield(md,'fieldname','slr.deltathickness','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]);
87 md = checkfield(md,'fieldname','slr.sealevel','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
88 md = checkfield(md,'fieldname','slr.love_h','NaN',1,'Inf',1);
89 md = checkfield(md,'fieldname','slr.love_k','NaN',1,'Inf',1);
90 md = checkfield(md,'fieldname','slr.love_l','NaN',1,'Inf',1);
91 md = checkfield(md,'fieldname','slr.tide_love_h','NaN',1,'Inf',1);
92 md = checkfield(md,'fieldname','slr.tide_love_k','NaN',1,'Inf',1);
93 md = checkfield(md,'fieldname','slr.fluid_love','NaN',1,'Inf',1);
94 md = checkfield(md,'fieldname','slr.equatorial_moi','NaN',1,'Inf',1);
95 md = checkfield(md,'fieldname','slr.polar_moi','NaN',1,'Inf',1);
96 md = checkfield(md,'fieldname','slr.angular_velocity','NaN',1,'Inf',1);
97 md = checkfield(md,'fieldname','slr.reltol','size',[1 1]);
98 md = checkfield(md,'fieldname','slr.abstol','size',[1 1]);
99 md = checkfield(md,'fieldname','slr.maxiter','size',[1 1],'>=',1);
100 md = checkfield(md,'fieldname','slr.steric_rate','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
101 md = checkfield(md,'fieldname','slr.degacc','size',[1 1],'>=',1e-10);
102 md = checkfield(md,'fieldname','slr.requested_outputs','stringrow',1);
103
104 %check that love numbers are provided at the same level of accuracy:
105 if (size(self.love_h,1)~=size(self.love_k,1) | size(self.love_h,1)~=size(self.love_l,1)),
106 error('slr error message: love numbers should be provided at the same level of accuracy');
107 end
108
109 %cross check that whereever we have an ice load, the mask is <0 on each vertex:
110 pos=find(self.deltathickness);
111 maskpos=md.mask.ice_levelset(md.mesh.elements(pos,:));
112 [els,vertices]=find(maskpos>0);
113 if length(els),
114 error('slr checkconsistency fail: there are elements with ice loads where some vertices are not on the ice!');
115 end
116
117 end % }}}
118 function list=defaultoutputs(self,md) % {{{
119 list = {'Sealevel'};
120 end % }}}
121 function disp(self) % {{{
122 disp(sprintf(' slr parameters:'));
123
124 fielddisplay(self,'deltathickness','thickness change: ice height equivalent [m]');
125 fielddisplay(self,'sealevel','current sea level (prior to computation) [m]');
126 fielddisplay(self,'reltol','sea level rise relative convergence criterion, (NaN: not applied)');
127 fielddisplay(self,'abstol','sea level rise absolute convergence criterion, (default, NaN: not applied)');
128 fielddisplay(self,'maxiter','maximum number of nonlinear iterations');
129 fielddisplay(self,'love_h','load Love number for radial displacement');
130 fielddisplay(self,'love_k','load Love number for gravitational potential perturbation');
131 fielddisplay(self,'love_l','load Love number for horizontal displacements');
132 fielddisplay(self,'tide_love_k','tidal load Love number (deg 2)');
133 fielddisplay(self,'tide_love_h','tidal load Love number (deg 2)');
134 fielddisplay(self,'fluid_love','secular fluid Love number');
135 fielddisplay(self,'equatorial_moi','mean equatorial moment of inertia [kg m^2]');
136 fielddisplay(self,'polar_moi','polar moment of inertia [kg m^2]');
137 fielddisplay(self,'angular_velocity','mean rotational velocity of earth [per second]');
138 fielddisplay(self,'rigid','rigid earth graviational potential perturbation');
139 fielddisplay(self,'elastic','elastic earth graviational potential perturbation');
140 fielddisplay(self,'rotation','earth rotational potential perturbation');
141 fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area (default: No correction)');
142 fielddisplay(self,'steric_rate','rate of steric ocean expansion [mm/yr]');
143 fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions');
144 fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps');
145 fielddisplay(self,'requested_outputs','additional outputs requested');
146
147 end % }}}
148 function marshall(self,prefix,md,fid) % {{{
149 %WriteData(fid,prefix,'object',self,'fieldname','deltathickness','format','DoubleMat','mattype',2);
150 WriteData(fid,prefix,'object',self,'fieldname','deltathickness','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
151 WriteData(fid,prefix,'object',self,'fieldname','sealevel','mattype',1,'format','DoubleMat','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
152 WriteData(fid,prefix,'object',self,'fieldname','reltol','format','Double');
153 WriteData(fid,prefix,'object',self,'fieldname','abstol','format','Double');
154 WriteData(fid,prefix,'object',self,'fieldname','maxiter','format','Integer');
155 WriteData(fid,prefix,'object',self,'fieldname','love_h','format','DoubleMat','mattype',1);
156 WriteData(fid,prefix,'object',self,'fieldname','love_k','format','DoubleMat','mattype',1);
157 WriteData(fid,prefix,'object',self,'fieldname','love_l','format','DoubleMat','mattype',1);
158 WriteData(fid,prefix,'object',self,'fieldname','tide_love_h','format','Double');
159 WriteData(fid,prefix,'object',self,'fieldname','tide_love_k','format','Double');
160 WriteData(fid,prefix,'object',self,'fieldname','fluid_love','format','Double');
161 WriteData(fid,prefix,'object',self,'fieldname','equatorial_moi','format','Double');
162 WriteData(fid,prefix,'object',self,'fieldname','polar_moi','format','Double');
163 WriteData(fid,prefix,'object',self,'fieldname','angular_velocity','format','Double');
164 WriteData(fid,prefix,'object',self,'fieldname','rigid','format','Boolean');
165 WriteData(fid,prefix,'object',self,'fieldname','elastic','format','Boolean');
166 WriteData(fid,prefix,'object',self,'fieldname','rotation','format','Boolean');
167 WriteData(fid,prefix,'object',self,'fieldname','ocean_area_scaling','format','Boolean');
168 WriteData(fid,prefix,'object',self,'fieldname','steric_rate','format','DoubleMat','mattype',1,'scale',1e-3/md.constants.yts);
169 WriteData(fid,prefix,'object',self,'fieldname','degacc','format','Double');
170 WriteData(fid,prefix,'object',self,'fieldname','transitions','format','MatArray');
171
172 %process requested outputs
173 outputs = self.requested_outputs;
174 pos = find(ismember(outputs,'default'));
175 if ~isempty(pos),
176 outputs(pos) = []; %remove 'default' from outputs
177 outputs = [outputs defaultoutputs(self,md)]; %add defaults
178 end
179 WriteData(fid,prefix,'data',outputs,'name','md.slr.requested_outputs','format','StringArray');
180
181 end % }}}
182 function savemodeljs(self,fid,modelname) % {{{
183
184 writejs1Darray(fid,[modelname '.slr.deltathickness'],self.deltathickness);
185 writejs1Darray(fid,[modelname '.slr.sealevel'],self.sealevel);
186 writejsdouble(fid,[modelname '.slr.maxiter'],self.maxiter);
187 writejsdouble(fid,[modelname '.slr.reltol'],self.reltol);
188 writejsdouble(fid,[modelname '.slr.abstol'],self.abstol);
189 writejs1Darray(fid,[modelname '.slr.love_h'],self.love_h);
190 writejs1Darray(fid,[modelname '.slr.love_k'],self.love_k);
191 writejs1Darray(fid,[modelname '.slr.love_l'],self.love_l);
192 writejsdouble(fid,[modelname '.slr.tide_love_k'],self.tide_love_k);
193 writejsdouble(fid,[modelname '.slr.tide_love_h'],self.tide_love_h);
194 writejsdouble(fid,[modelname '.slr.fluid_love'],self.fluid_love);
195 writejsdouble(fid,[modelname '.slr.equatorial_moi'],self.equatorial_moi);
196 writejsdouble(fid,[modelname '.slr.polar_moi'],self.polar_moi);
197 writejsdouble(fid,[modelname '.slr.angular_velocity'],self.angular_velocity);
198 writejsdouble(fid,[modelname '.slr.rigid'],self.rigid);
199 writejsdouble(fid,[modelname '.slr.elastic'],self.elastic);
200 writejsdouble(fid,[modelname '.slr.rotation'],self.rotation);
201 writejsdouble(fid,[modelname '.slr.ocean_area_scaling'],self.ocean_area_scaling);
202 writejs1Darray(fid,[modelname '.slr.steric_rate'],self.steric_rate);
203 writejsdouble(fid,[modelname '.slr.degacc'],self.degacc);
204 writejscellstring(fid,[modelname '.slr.requested_outputs'],self.requested_outputs);
205 writejscellarray(fid,[modelname '.slr.transitions'],self.transitions);
206 end % }}}
207 end
208end
Note: See TracBrowser for help on using the repository browser.