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

%SLR class definition
%
%   Usage:
%      slr=slr();

classdef slr
        properties (SetAccess=public) 
                deltathickness = NaN;
                sealevel       = NaN; 
                maxiter        = 0;
                reltol         = 0;
                abstol         = 0;
                love_h         = 0; %provided by PREM model
                love_k         = 0; %ideam
                rigid         = 0;
                elastic         = 0;
                eustatic         = 0;
                degacc         = 0;
                requested_outputs      = {};
                transitions    = {};
        end
        methods
                function self = slr(varargin) % {{{
                        switch nargin
                                case 0
                                        self=setdefaultparameters(self);
                                otherwise
                                        error('constructor not supported');
                        end
                end % }}}
                function self = setdefaultparameters(self) % {{{
                
                %Convergence criterion: absolute, relative and residual
                self.reltol=NaN; %default
                self.abstol=0.001; %1 mm of sea level rise

                %maximum of non-linear iterations.
                self.maxiter=10;

                %computational flags: 
                self.rigid=1;
                self.elastic=1;
                self.eustatic=1;

                %numerical discretization accuracy
                self.degacc=.01;
                
                %output default:
                self.requested_outputs={'default'};

                %transitions should be a cell array of vectors: 
                self.transitions={};
                
                end % }}}
                function md = checkconsistency(self,md,solution,analyses) % {{{

                        if ~ismember(SealevelriseAnalysisEnum(),analyses), return; end
                        md = checkfield(md,'fieldname','slr.deltathickness','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]);
                        %md = checkfield(md,'fieldname','slr.deltathickness','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
                        md = checkfield(md,'fieldname','slr.sealevel','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
                        md = checkfield(md,'fieldname','slr.love_h','NaN',1,'Inf',1);
                        md = checkfield(md,'fieldname','slr.love_k','NaN',1,'Inf',1);
                        md = checkfield(md,'fieldname','slr.reltol','size',[1 1]);
                        md = checkfield(md,'fieldname','slr.abstol','size',[1 1]);
                        md = checkfield(md,'fieldname','slr.maxiter','size',[1 1],'>=',1);
                        md = checkfield(md,'fieldname','slr.degacc','size',[1 1],'>=',1e-10);
                        md = checkfield(md,'fieldname','slr.requested_outputs','stringrow',1);

                        %check that love numbers are provided at the same level of accuracy: 
                        if (size(self.love_h,1) ~= size(self.love_k,1)),
                                error('slr error message: love numbers should be provided at the same level of accuracy');
                        end

                        %cross check that whereever we have an ice load, the mask is <0 on each vertex: 
                        pos=find(self.deltathickness);
                        maskpos=md.mask.ice_levelset(md.mesh.elements(pos,:)); 
                        [els,vertices]=find(maskpos>0);
                        if length(els),
                                error('slr checkconsistency fail: there are elements with ice loads where some vertices are not on the ice!');
                        end

                end % }}}
                function list=defaultoutputs(self,md) % {{{
                        list = {'Sealevel'};
                end % }}}
                function disp(self) % {{{
                        disp(sprintf('   slr parameters:'));

                        fielddisplay(self,'deltathickness','thickness change (main loading of the slr solution core [m]');
                        fielddisplay(self,'sealevel','current sea level (prior to computation) [m]');
                        fielddisplay(self,'reltol','sea level rise relative convergence criterion, (default, NaN: not applied)');
                        fielddisplay(self,'abstol','sea level rise absolute convergence criterion, NaN: not applied');
                        fielddisplay(self,'maxiter','maximum number of nonlinear iterations');
                        fielddisplay(self,'love_h','love load number for radial displacement');
                        fielddisplay(self,'love_k','love load number for gravitational potential perturbation');
                        fielddisplay(self,'rigid','rigid earth graviational potential perturbation');
                        fielddisplay(self,'elastic','elastic earth graviational potential perturbation');
                        fielddisplay(self,'eustatic','eustatic sea level rise');
                        fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions');
                        fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps');
                        fielddisplay(self,'requested_outputs','additional outputs requested');

                end % }}}
                function marshall(self,md,fid) % {{{
                        WriteData(fid,'object',self,'class','sealevelrise','fieldname','deltathickness','format','DoubleMat','mattype',2);
                        %WriteData(fid,'object',self,'class','sealevelrise','fieldname','deltathickness','format','DoubleMat','mattype',1);
                        WriteData(fid,'data',self.sealevel,'mattype',1,'format','DoubleMat','enum',SealevelEnum(),'timeserieslength',md.mesh.numberofvertices+1);
                        WriteData(fid,'object',self,'class','sealevelrise','fieldname','reltol','format','Double');
                        WriteData(fid,'object',self,'class','sealevelrise','fieldname','abstol','format','Double');
                        WriteData(fid,'object',self,'class','sealevelrise','fieldname','maxiter','format','Integer');
                        WriteData(fid,'object',self,'class','sealevelrise','fieldname','love_h','format','DoubleMat','mattype',1);
                        WriteData(fid,'object',self,'class','sealevelrise','fieldname','love_k','format','DoubleMat','mattype',1);
                        WriteData(fid,'object',self,'class','sealevelrise','fieldname','rigid','format','Boolean');
                        WriteData(fid,'object',self,'class','sealevelrise','fieldname','elastic','format','Boolean');
                        WriteData(fid,'object',self,'class','sealevelrise','fieldname','eustatic','format','Boolean');
                        WriteData(fid,'object',self,'class','sealevelrise','fieldname','degacc','format','Double');
                        WriteData(fid,'object',self,'class','sealevelrise','fieldname','transitions','format','MatArray');
                        
                        %process requested outputs
                        outputs = self.requested_outputs;
                        pos  = find(ismember(outputs,'default'));
                        if ~isempty(pos),
                                outputs(pos) = [];                         %remove 'default' from outputs
                                outputs      = [outputs defaultoutputs(self,md)]; %add defaults
                        end
                        WriteData(fid,'data',outputs,'enum',SealevelriseRequestedOutputsEnum,'format','StringArray');

                end % }}}
                function savemodeljs(self,fid,modelname) % {{{
                
                        writejs1Darray(fid,[modelname '.slr.deltathickness'],self.deltathickness);
                        writejs1Darray(fid,[modelname '.slr.sealevel'],self.sealevel);
                        writejsdouble(fid,[modelname '.slr.maxiter'],self.maxiter);
                        writejsdouble(fid,[modelname '.slr.reltol'],self.reltol);
                        writejsdouble(fid,[modelname '.slr.abstol'],self.abstol);
                        writejs1Darray(fid,[modelname '.slr.love_h'],self.love_h);
                        writejs1Darray(fid,[modelname '.slr.love_k'],self.love_k);
                        writejsdouble(fid,[modelname '.slr.rigid'],self.rigid);
                        writejsdouble(fid,[modelname '.slr.eustatic'],self.eustatic);
                        writejsdouble(fid,[modelname '.slr.elastic'],self.elastic);
                        writejsdouble(fid,[modelname '.slr.degacc'],self.degacc);
                        writejscellstring(fid,[modelname '.slr.requested_outputs'],self.requested_outputs);
                        writejscellarray(fid,[modelname '.slr.transitions'],self.transitions);
                end % }}}
        end
end