Context Navigation

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

Last change on this file since 22107 was 22107, checked in by adhikari, 8 years ago
CHG: reltol set as default criterion for RSL convergence
File size: 9.8 KB

Rev	Line
[19984]	1	%SLR class definition
	2	%
	3	% Usage:
	4	% slr=slr();
	5
	6	classdef slr
	7	properties (SetAccess=public)
	8	deltathickness = NaN;
[20154]	9	sealevel = NaN;
[19984]	10	maxiter = 0;
	11	reltol = 0;
	12	abstol = 0;
	13	love_h = 0; %provided by PREM model
	14	love_k = 0; %ideam
[21000]	15	love_l = 0; %ideam
[20349]	16	tide_love_k = 0; %ideam
	17	tide_love_h = 0; %ideam
[21345]	18	fluid_love = 0;
	19	equatorial_moi = 0;
	20	polar_moi = 0;
[21331]	21	angular_velocity = 0;
[20349]	22	rigid = 0;
	23	elastic = 0;
	24	rotation = 0;
[21295]	25	ocean_area_scaling = 0;
[21752]	26	steric_rate = 0; %rate of ocean expansion from steric effects.
[20033]	27	degacc = 0;
[20036]	28	requested_outputs = {};
[20137]	29	transitions = {};
[19984]	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
[22107]	43	self.reltol=0.01; % 1 per cent
	44	self.abstol=NaN; % default
[19984]	45
	46	%maximum of non-linear iterations.
[22107]	47	self.maxiter=5;
[19984]	48
	49	%computational flags:
	50	self.rigid=1;
	51	self.elastic=1;
[21331]	52	self.rotation=0;
[21309]	53	self.ocean_area_scaling=0;
[20033]	54
[20349]	55	%tidal love numbers:
	56	self.tide_love_h=0.6149; %degree 2
	57	self.tide_love_k=0.3055; % degree 2
[21345]	58
	59	%secular fluid love number:
	60	self.fluid_love=0.942;
[21331]	61
[21345]	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
[21331]	66	% mean rotational velocity of earth
	67	self.angular_velocity=7.2921*10^-5; % [s^-1]
[20349]	68
[20033]	69	%numerical discretization accuracy
	70	self.degacc=.01;
[21752]	71
	72	%steric:
	73	self.steric_rate=0;
	74
	75
[20036]	76	%output default:
	77	self.requested_outputs={'default'};
[20137]	78
	79	%transitions should be a cell array of vectors:
	80	self.transitions={};
[20036]	81
[19984]	82	end % }}}
	83	function md = checkconsistency(self,md,solution,analyses) % {{{
	84
[21049]	85	if ~ismember('SealevelriseAnalysis',analyses), return; end
[20322]	86	md = checkfield(md,'fieldname','slr.deltathickness','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]);
[20154]	87	md = checkfield(md,'fieldname','slr.sealevel','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
[19984]	88	md = checkfield(md,'fieldname','slr.love_h','NaN',1,'Inf',1);
	89	md = checkfield(md,'fieldname','slr.love_k','NaN',1,'Inf',1);
[21000]	90	md = checkfield(md,'fieldname','slr.love_l','NaN',1,'Inf',1);
[20349]	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);
[21345]	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);
[21331]	96	md = checkfield(md,'fieldname','slr.angular_velocity','NaN',1,'Inf',1);
[19984]	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);
[21752]	100	md = checkfield(md,'fieldname','slr.steric_rate','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
[20033]	101	md = checkfield(md,'fieldname','slr.degacc','size',[1 1],'>=',1e-10);
[20036]	102	md = checkfield(md,'fieldname','slr.requested_outputs','stringrow',1);
[19984]	103
	104	%check that love numbers are provided at the same level of accuracy:
[21000]	105	if (size(self.love_h,1)~=size(self.love_k,1) \| size(self.love_h,1)~=size(self.love_l,1)),
[19984]	106	error('slr error message: love numbers should be provided at the same level of accuracy');
	107	end
	108
[20322]	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
[19984]	117	end % }}}
[20036]	118	function list=defaultoutputs(self,md) % {{{
[20154]	119	list = {'Sealevel'};
[20036]	120	end % }}}
[19984]	121	function disp(self) % {{{
	122	disp(sprintf(' slr parameters:'));
	123
	124	fielddisplay(self,'deltathickness','thickness change (main loading of the slr solution core [m]');
[20154]	125	fielddisplay(self,'sealevel','current sea level (prior to computation) [m]');
[22107]	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');
[19984]	128	fielddisplay(self,'maxiter','maximum number of nonlinear iterations');
[21000]	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)');
[21345]	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]');
[21331]	137	fielddisplay(self,'angular_velocity','mean rotational velocity of earth [per second]');
[19984]	138	fielddisplay(self,'rigid','rigid earth graviational potential perturbation');
	139	fielddisplay(self,'elastic','elastic earth graviational potential perturbation');
[21295]	140	fielddisplay(self,'rotation','earth rotational potential perturbation');
[21309]	141	fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area [default: No correction]');
[21752]	142	fielddisplay(self,'steric_rate','rate of steric ocean expansion (in mm/yr)');
[20033]	143	fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions');
[20137]	144	fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps');
[20036]	145	fielddisplay(self,'requested_outputs','additional outputs requested');
[19984]	146
	147	end % }}}
[20690]	148	function marshall(self,prefix,md,fid) % {{{
[21939]	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);
[20902]	151	WriteData(fid,prefix,'object',self,'fieldname','sealevel','mattype',1,'format','DoubleMat','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
[20690]	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);
[21000]	157	WriteData(fid,prefix,'object',self,'fieldname','love_l','format','DoubleMat','mattype',1);
[20891]	158	WriteData(fid,prefix,'object',self,'fieldname','tide_love_h','format','Double');
[20690]	159	WriteData(fid,prefix,'object',self,'fieldname','tide_love_k','format','Double');
[21345]	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');
[21331]	163	WriteData(fid,prefix,'object',self,'fieldname','angular_velocity','format','Double');
[20690]	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');
[21295]	167	WriteData(fid,prefix,'object',self,'fieldname','ocean_area_scaling','format','Boolean');
[21752]	168	WriteData(fid,prefix,'object',self,'fieldname','steric_rate','format','DoubleMat','mattype',1,'scale',1e-3/md.constants.yts);
[20690]	169	WriteData(fid,prefix,'object',self,'fieldname','degacc','format','Double');
	170	WriteData(fid,prefix,'object',self,'fieldname','transitions','format','MatArray');
[20036]	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
[20690]	179	WriteData(fid,prefix,'data',outputs,'name','md.slr.requested_outputs','format','StringArray');
[20036]	180
[19984]	181	end % }}}
	182	function savemodeljs(self,fid,modelname) % {{{
	183
[20181]	184	writejs1Darray(fid,[modelname '.slr.deltathickness'],self.deltathickness);
	185	writejs1Darray(fid,[modelname '.slr.sealevel'],self.sealevel);
[20262]	186	writejsdouble(fid,[modelname '.slr.maxiter'],self.maxiter);
[19984]	187	writejsdouble(fid,[modelname '.slr.reltol'],self.reltol);
	188	writejsdouble(fid,[modelname '.slr.abstol'],self.abstol);
[20181]	189	writejs1Darray(fid,[modelname '.slr.love_h'],self.love_h);
	190	writejs1Darray(fid,[modelname '.slr.love_k'],self.love_k);
[21000]	191	writejs1Darray(fid,[modelname '.slr.love_l'],self.love_l);
[20349]	192	writejsdouble(fid,[modelname '.slr.tide_love_k'],self.tide_love_k);
	193	writejsdouble(fid,[modelname '.slr.tide_love_h'],self.tide_love_h);
[21345]	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);
[21331]	197	writejsdouble(fid,[modelname '.slr.angular_velocity'],self.angular_velocity);
[19984]	198	writejsdouble(fid,[modelname '.slr.rigid'],self.rigid);
[21295]	199	writejsdouble(fid,[modelname '.slr.elastic'],self.elastic);
[20349]	200	writejsdouble(fid,[modelname '.slr.rotation'],self.rotation);
[21295]	201	writejsdouble(fid,[modelname '.slr.ocean_area_scaling'],self.ocean_area_scaling);
[21752]	202	writejs1Darray(fid,[modelname '.slr.steric_rate'],self.steric_rate);
[20033]	203	writejsdouble(fid,[modelname '.slr.degacc'],self.degacc);
[20036]	204	writejscellstring(fid,[modelname '.slr.requested_outputs'],self.requested_outputs);
[20262]	205	writejscellarray(fid,[modelname '.slr.transitions'],self.transitions);
[19984]	206	end % }}}
	207	end
	208	end

Note: See TracBrowser for help on using the repository browser.

Download in other formats: