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

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

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