source: issm/branches/trunk-larour-NatGeoScience2016/src/m/classes/slr.m@ 21326

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

CHG: new steric rate formualtion.

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