Context Navigation

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

Last change on this file since 24729 was 24729, checked in by Eric.Larour, 5 years ago
CHG: added a check on the type of mask that should be initialized for an slr solution.
File size: 13.6 KB

Rev	Line
[19984]	1	%SLR class definition
	2	%
	3	% Usage:
	4	% slr=slr();
	5
	6	classdef slr
	7	properties (SetAccess=public)
[22970]	8	deltathickness = NaN;
	9	sealevel = NaN;
	10	spcthickness = NaN;
	11	maxiter = 0;
	12	reltol = 0;
	13	abstol = 0;
	14	love_h = 0; %provided by PREM model
	15	love_k = 0; %ideam
	16	love_l = 0; %ideam
	17	tide_love_k = 0; %ideam
	18	tide_love_h = 0; %ideam
	19	fluid_love = 0;
	20	equatorial_moi = 0;
	21	polar_moi = 0;
	22	angular_velocity = 0;
	23	rigid = 0;
	24	elastic = 0;
	25	rotation = 0;
	26	ocean_area_scaling = 0;
[24259]	27	hydro_rate = 0; %rate of steric expansion from hydrological effects.
[22970]	28	geodetic_run_frequency = 1; %how many time steps we skip before we run the geodetic part of the solver during transient
	29	geodetic = 0; %compute geodetic SLR? (in addition to steric?)
	30	degacc = 0;
	31	loop_increment = 0;
	32	horiz = 0;
	33	Ngia = NaN;
	34	Ugia = NaN;
[24469]	35
[20036]	36	requested_outputs = {};
[22970]	37	transitions = {};
[19984]	38	end
	39	methods
	40	function self = slr(varargin) % {{{
	41	switch nargin
	42	case 0
	43	self=setdefaultparameters(self);
	44	otherwise
	45	error('constructor not supported');
	46	end
	47	end % }}}
	48	function self = setdefaultparameters(self) % {{{
	49
	50	%Convergence criterion: absolute, relative and residual
[22107]	51	self.reltol=0.01; % 1 per cent
	52	self.abstol=NaN; % default
[19984]	53
	54	%maximum of non-linear iterations.
[22107]	55	self.maxiter=5;
[22955]	56	self.loop_increment=200;
[19984]	57
	58	%computational flags:
[22955]	59	self.geodetic=0;
[19984]	60	self.rigid=1;
	61	self.elastic=1;
[21309]	62	self.ocean_area_scaling=0;
[22955]	63	self.rotation=1;
[20033]	64
[20349]	65	%tidal love numbers:
	66	self.tide_love_h=0.6149; %degree 2
	67	self.tide_love_k=0.3055; % degree 2
[21345]	68
	69	%secular fluid love number:
	70	self.fluid_love=0.942;
[21331]	71
[21345]	72	%moment of inertia:
	73	self.equatorial_moi=8.0077*10^37; % [kg m^2]
	74	self.polar_moi =8.0345*10^37; % [kg m^2]
	75
[21331]	76	% mean rotational velocity of earth
	77	self.angular_velocity=7.2921*10^-5; % [s^-1]
[20349]	78
[20033]	79	%numerical discretization accuracy
	80	self.degacc=.01;
[21752]	81
[24469]	82	%hydro
[24259]	83	self.hydro_rate=0;
[21752]	84
[22955]	85	%how many time steps we skip before we run SLR solver during transient
	86	self.geodetic_run_frequency=1;
[21752]	87
[20036]	88	%output default:
	89	self.requested_outputs={'default'};
[20137]	90
	91	%transitions should be a cell array of vectors:
	92	self.transitions={};
[22955]	93
	94	%horizontal displacement? (not by default)
	95	self.horiz=0;
[20036]	96
[19984]	97	end % }}}
	98	function md = checkconsistency(self,md,solution,analyses) % {{{
	99
[24519]	100	if ~ismember('SealevelriseAnalysis',analyses) \| (strcmp(solution,'TransientSolution') & md.transient.isslr==0),
	101	return;
	102	end
[24729]	103
	104	%check that the mask is a maskpsl class type, otherwise, we won't be able to run very far
	105	%with the slr solution:
	106	if ~strcmpi(class(md.mask),'maskpsl'),
	107	error('slr error message: model mask should be of the ''maskpsl'' class');
	108	end
	109
[20322]	110	md = checkfield(md,'fieldname','slr.deltathickness','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]);
[20154]	111	md = checkfield(md,'fieldname','slr.sealevel','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
[22955]	112	md = checkfield(md,'fieldname','slr.spcthickness','Inf',1,'timeseries',1);
[19984]	113	md = checkfield(md,'fieldname','slr.love_h','NaN',1,'Inf',1);
	114	md = checkfield(md,'fieldname','slr.love_k','NaN',1,'Inf',1);
[21000]	115	md = checkfield(md,'fieldname','slr.love_l','NaN',1,'Inf',1);
[20349]	116	md = checkfield(md,'fieldname','slr.tide_love_h','NaN',1,'Inf',1);
	117	md = checkfield(md,'fieldname','slr.tide_love_k','NaN',1,'Inf',1);
[21345]	118	md = checkfield(md,'fieldname','slr.fluid_love','NaN',1,'Inf',1);
	119	md = checkfield(md,'fieldname','slr.equatorial_moi','NaN',1,'Inf',1);
	120	md = checkfield(md,'fieldname','slr.polar_moi','NaN',1,'Inf',1);
[21331]	121	md = checkfield(md,'fieldname','slr.angular_velocity','NaN',1,'Inf',1);
[19984]	122	md = checkfield(md,'fieldname','slr.reltol','size',[1 1]);
	123	md = checkfield(md,'fieldname','slr.abstol','size',[1 1]);
	124	md = checkfield(md,'fieldname','slr.maxiter','size',[1 1],'>=',1);
[22955]	125	md = checkfield(md,'fieldname','slr.geodetic_run_frequency','size',[1 1],'>=',1);
[24259]	126	md = checkfield(md,'fieldname','slr.hydro_rate','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
[20033]	127	md = checkfield(md,'fieldname','slr.degacc','size',[1 1],'>=',1e-10);
[20036]	128	md = checkfield(md,'fieldname','slr.requested_outputs','stringrow',1);
[22955]	129	md = checkfield(md,'fieldname','slr.loop_increment','NaN',1,'Inf',1,'>=',1);
	130	md = checkfield(md,'fieldname','slr.horiz','NaN',1,'Inf',1,'values',[0 1]);
	131	md = checkfield(md,'fieldname','slr.Ngia','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
	132	md = checkfield(md,'fieldname','slr.Ugia','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
[19984]	133
	134	%check that love numbers are provided at the same level of accuracy:
[21000]	135	if (size(self.love_h,1)~=size(self.love_k,1) \| size(self.love_h,1)~=size(self.love_l,1)),
[19984]	136	error('slr error message: love numbers should be provided at the same level of accuracy');
	137	end
	138
[20322]	139	%cross check that whereever we have an ice load, the mask is <0 on each vertex:
	140	pos=find(self.deltathickness);
	141	maskpos=md.mask.ice_levelset(md.mesh.elements(pos,:));
	142	[els,vertices]=find(maskpos>0);
	143	if length(els),
[22955]	144	warning('slr checkconsistency fail: there are elements with ice loads where some vertices are not on the ice!');
[20322]	145	end
	146
[22955]	147	%check that if geodetic is requested, we are a mesh3dsurface model (planet), or if we are not,
	148	%a coupler to a planet model is provided.
	149	if self.geodetic,
	150	if md.transient.iscoupler,
	151	%we are good;
	152	else
	153	if strcmpi(class(md.mesh),'mesh3dsurface'),
	154	%we are good
	155	else
	156	error('model is requesting geodetic computations without being a mesh3dsurface, or being coupled to one!');
	157	end
	158	end
	159	end
	160
[19984]	161	end % }}}
[20036]	162	function list=defaultoutputs(self,md) % {{{
[20154]	163	list = {'Sealevel'};
[20036]	164	end % }}}
[19984]	165	function disp(self) % {{{
	166	disp(sprintf(' slr parameters:'));
	167
[22808]	168	fielddisplay(self,'deltathickness','thickness change: ice height equivalent [m]');
[20154]	169	fielddisplay(self,'sealevel','current sea level (prior to computation) [m]');
[22955]	170	fielddisplay(self,'spcthickness','thickness constraints (NaN means no constraint) [m]');
	171	fielddisplay(self,'reltol','sea level rise relative convergence criterion, (default, NaN: not applied)');
	172	fielddisplay(self,'abstol','sea level rise absolute convergence criterion, NaN: not applied');
[19984]	173	fielddisplay(self,'maxiter','maximum number of nonlinear iterations');
[21000]	174	fielddisplay(self,'love_h','load Love number for radial displacement');
	175	fielddisplay(self,'love_k','load Love number for gravitational potential perturbation');
	176	fielddisplay(self,'love_l','load Love number for horizontal displacements');
	177	fielddisplay(self,'tide_love_k','tidal load Love number (deg 2)');
	178	fielddisplay(self,'tide_love_h','tidal load Love number (deg 2)');
[21345]	179	fielddisplay(self,'fluid_love','secular fluid Love number');
	180	fielddisplay(self,'equatorial_moi','mean equatorial moment of inertia [kg m^2]');
	181	fielddisplay(self,'polar_moi','polar moment of inertia [kg m^2]');
[21331]	182	fielddisplay(self,'angular_velocity','mean rotational velocity of earth [per second]');
[22955]	183	fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area [default: No correction]');
[24259]	184	fielddisplay(self,'hydro_rate','rate of hydrological expansion (in mm/yr)');
[22955]	185	fielddisplay(self,'Ngia','rate of viscous (GIA) geoid expansion (in mm/yr)');
	186	fielddisplay(self,'Ugia','rate of viscous (GIA) bedrock uplift (in mm/yr)');
	187	fielddisplay(self,'loop_increment','vector assembly (in the convolution) framentation');
	188	fielddisplay(self,'geodetic','compute geodetic SLR? ( in addition to steric?) default 0');
	189	fielddisplay(self,'geodetic_run_frequency','how many time steps we skip before we run SLR solver during transient (default: 1)');
[19984]	190	fielddisplay(self,'rigid','rigid earth graviational potential perturbation');
	191	fielddisplay(self,'elastic','elastic earth graviational potential perturbation');
[21295]	192	fielddisplay(self,'rotation','earth rotational potential perturbation');
[20033]	193	fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions');
[20137]	194	fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps');
[20036]	195	fielddisplay(self,'requested_outputs','additional outputs requested');
[19984]	196
	197	end % }}}
[20690]	198	function marshall(self,prefix,md,fid) % {{{
[21939]	199	%WriteData(fid,prefix,'object',self,'fieldname','deltathickness','format','DoubleMat','mattype',2);
	200	WriteData(fid,prefix,'object',self,'fieldname','deltathickness','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
[22955]	201	%WriteData(fid,prefix,'object',self,'fieldname','deltathickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofelements+1);
[20902]	202	WriteData(fid,prefix,'object',self,'fieldname','sealevel','mattype',1,'format','DoubleMat','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
[22955]	203	WriteData(fid,prefix,'object',self,'fieldname','spcthickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
[20690]	204	WriteData(fid,prefix,'object',self,'fieldname','reltol','format','Double');
	205	WriteData(fid,prefix,'object',self,'fieldname','abstol','format','Double');
	206	WriteData(fid,prefix,'object',self,'fieldname','maxiter','format','Integer');
	207	WriteData(fid,prefix,'object',self,'fieldname','love_h','format','DoubleMat','mattype',1);
	208	WriteData(fid,prefix,'object',self,'fieldname','love_k','format','DoubleMat','mattype',1);
[21000]	209	WriteData(fid,prefix,'object',self,'fieldname','love_l','format','DoubleMat','mattype',1);
[20891]	210	WriteData(fid,prefix,'object',self,'fieldname','tide_love_h','format','Double');
[20690]	211	WriteData(fid,prefix,'object',self,'fieldname','tide_love_k','format','Double');
[21345]	212	WriteData(fid,prefix,'object',self,'fieldname','fluid_love','format','Double');
	213	WriteData(fid,prefix,'object',self,'fieldname','equatorial_moi','format','Double');
	214	WriteData(fid,prefix,'object',self,'fieldname','polar_moi','format','Double');
[21331]	215	WriteData(fid,prefix,'object',self,'fieldname','angular_velocity','format','Double');
[20690]	216	WriteData(fid,prefix,'object',self,'fieldname','rigid','format','Boolean');
	217	WriteData(fid,prefix,'object',self,'fieldname','elastic','format','Boolean');
	218	WriteData(fid,prefix,'object',self,'fieldname','rotation','format','Boolean');
[21295]	219	WriteData(fid,prefix,'object',self,'fieldname','ocean_area_scaling','format','Boolean');
[22955]	220	WriteData(fid,prefix,'object',self,'fieldname','geodetic_run_frequency','format','Integer');
[24259]	221	WriteData(fid,prefix,'object',self,'fieldname','hydro_rate','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1e-3/md.constants.yts);
[22955]	222	WriteData(fid,prefix,'object',self,'fieldname','Ngia','format','DoubleMat','mattype',1,'scale',1e-3/md.constants.yts);
	223	WriteData(fid,prefix,'object',self,'fieldname','Ugia','format','DoubleMat','mattype',1,'scale',1e-3/md.constants.yts);
[20690]	224	WriteData(fid,prefix,'object',self,'fieldname','degacc','format','Double');
	225	WriteData(fid,prefix,'object',self,'fieldname','transitions','format','MatArray');
[22955]	226	WriteData(fid,prefix,'object',self,'fieldname','loop_increment','format','Integer');
	227	WriteData(fid,prefix,'object',self,'fieldname','horiz','format','Integer');
	228	WriteData(fid,prefix,'object',self,'fieldname','geodetic','format','Integer');
[20036]	229
	230	%process requested outputs
	231	outputs = self.requested_outputs;
	232	pos = find(ismember(outputs,'default'));
	233	if ~isempty(pos),
	234	outputs(pos) = []; %remove 'default' from outputs
	235	outputs = [outputs defaultoutputs(self,md)]; %add defaults
	236	end
[20690]	237	WriteData(fid,prefix,'data',outputs,'name','md.slr.requested_outputs','format','StringArray');
[20036]	238
[19984]	239	end % }}}
	240	function savemodeljs(self,fid,modelname) % {{{
	241
[20181]	242	writejs1Darray(fid,[modelname '.slr.deltathickness'],self.deltathickness);
	243	writejs1Darray(fid,[modelname '.slr.sealevel'],self.sealevel);
[22955]	244	writejs1Darray(fid,[modelname '.slr.spcthickness'],self.spcthickness);
[20262]	245	writejsdouble(fid,[modelname '.slr.maxiter'],self.maxiter);
[19984]	246	writejsdouble(fid,[modelname '.slr.reltol'],self.reltol);
	247	writejsdouble(fid,[modelname '.slr.abstol'],self.abstol);
[20181]	248	writejs1Darray(fid,[modelname '.slr.love_h'],self.love_h);
	249	writejs1Darray(fid,[modelname '.slr.love_k'],self.love_k);
[21000]	250	writejs1Darray(fid,[modelname '.slr.love_l'],self.love_l);
[20349]	251	writejsdouble(fid,[modelname '.slr.tide_love_k'],self.tide_love_k);
	252	writejsdouble(fid,[modelname '.slr.tide_love_h'],self.tide_love_h);
[21345]	253	writejsdouble(fid,[modelname '.slr.fluid_love'],self.fluid_love);
	254	writejsdouble(fid,[modelname '.slr.equatorial_moi'],self.equatorial_moi);
	255	writejsdouble(fid,[modelname '.slr.polar_moi'],self.polar_moi);
[21331]	256	writejsdouble(fid,[modelname '.slr.angular_velocity'],self.angular_velocity);
[19984]	257	writejsdouble(fid,[modelname '.slr.rigid'],self.rigid);
[21295]	258	writejsdouble(fid,[modelname '.slr.elastic'],self.elastic);
[20349]	259	writejsdouble(fid,[modelname '.slr.rotation'],self.rotation);
[21295]	260	writejsdouble(fid,[modelname '.slr.ocean_area_scaling'],self.ocean_area_scaling);
[22955]	261	writejsdouble(fid,[modelname '.slr.geodetic_run_frequency'],self.geodetic_run_frequency);
[24259]	262	writejs1Darray(fid,[modelname '.slr.hydro_rate'],self.hydro_rate);
[20033]	263	writejsdouble(fid,[modelname '.slr.degacc'],self.degacc);
[20036]	264	writejscellstring(fid,[modelname '.slr.requested_outputs'],self.requested_outputs);
[20262]	265	writejscellarray(fid,[modelname '.slr.transitions'],self.transitions);
[19984]	266	end % }}}
[22955]	267	function self = extrude(self,md) % {{{
	268	self.sealevel=project3d(md,'vector',self.sealevel,'type','node');
	269	end % }}}
[19984]	270	end
	271	end

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

Download in other formats: