Context Navigation

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

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

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

Download in other formats: