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

Last change on this file since 21309 was 21309, checked in by adhikari, 8 years ago

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