Context Navigation

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

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

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

Download in other formats: