Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26232) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26233) @@ -250,13 +250,23 @@ LockFileNameEnum, LoveAllowLayerDeletionEnum, + LoveCoreMantleBoundaryEnum, + LoveEarthMassEnum, LoveForcingTypeEnum, LoveFrequenciesEnum, + LoveIsTemporalEnum, LoveG0Enum, + LoveGravitationalConstantEnum, + LoveInnerCoreBoundaryEnum, + LoveIntStepsPerLayerEnum, LoveKernelsEnum, LoveMu0Enum, LoveNfreqEnum, + LoveNTemporalIterationsEnum, + LoveNYiEquationsEnum, LoveR0Enum, LoveShNmaxEnum, LoveShNminEnum, + LoveStartingLayerEnum, + LoveUnderflowTolEnum, MassFluxSegmentsEnum, MassFluxSegmentsPresentEnum, @@ -779,4 +789,13 @@ SealevelGrotm2Enum, SealevelGrotm3Enum, + SealevelGUrotm1Enum, + SealevelGUrotm2Enum, + SealevelGUrotm3Enum, + SealevelGNrotm1Enum, + SealevelGNrotm2Enum, + SealevelGNrotm3Enum, + SealevelGErotm1Enum, + SealevelGErotm2Enum, + SealevelGErotm3Enum, SealevelRSLBarystaticEnum, SealevelRSLRateEnum, Index: /issm/trunk-jpl/src/m/classes/fourierlove.m =================================================================== --- /issm/trunk-jpl/src/m/classes/fourierlove.m (revision 26232) +++ /issm/trunk-jpl/src/m/classes/fourierlove.m (revision 26233) @@ -6,14 +6,22 @@ classdef fourierlove properties (SetAccess=public) - nfreq = 0; - frequencies = 0; - sh_nmax = 0; - sh_nmin = 0; - g0 = 0; - r0 = 0; - mu0 = 0; - allow_layer_deletion = 0; - love_kernels = 0; - forcing_type = 0; + nfreq = 0; + frequencies = 0; + sh_nmax = 0; + sh_nmin = 0; + g0 = 0; + r0 = 0; + mu0 = 0; + Gravitational_Constant = 0; + allow_layer_deletion = 0; + underflow_tol = 0; + integration_steps_per_layer= 0; + istemporal = 0; + n_temporal_iterations = 0; + time = 0; + love_kernels = 0; + forcing_type = 0; + inner_core_boundary = 0; + core_mantle_boundary = 0; end methods (Static) @@ -44,7 +52,15 @@ self.r0=6371*1e3; %m; self.mu0=10^11; % Pa - self.allow_layer_deletion=1; + self.Gravitational_Constant=6.67259e-11; % m^3 kg^-1 s^-2 + self.allow_layer_deletion=1; + self.underflow_tol=1e-16; %threshold of deep to surface love number ratio to trigger the deletion of layer + self.integration_steps_per_layer=100; + self.istemporal=0; + self.n_temporal_iterations=8; + self.time=[0]; %s self.love_kernels=0; - self.forcing_type = 11; + self.forcing_type = 11; % surface loading + self.inner_core_boundary=1; + self.core_mantle_boundary=2; end % }}} function disp(self) % {{{ @@ -54,9 +70,17 @@ fielddisplay(self,'sh_nmin','minimum spherical harmonic degree (default 1)'); fielddisplay(self,'g0','adimensioning constant for gravity (default 10) [m/s^2]'); - fielddisplay(self,'r0','adimensioning constant for radius (default 6378*10^3) [m]'); + fielddisplay(self,'r0','adimensioning constant for radius (default 6371*10^3) [m]'); fielddisplay(self,'mu0','adimensioning constant for stress (default 10^11) [Pa]'); + fielddisplay(self,'Gravitational_Constant','Newtonian constant of gravitation (default 6.67259e-11 [m^3 kg^-1 s^-2])'); fielddisplay(self,'allow_layer_deletion','allow for migration of the integration boundary with increasing spherical harmonics degree (default 1)'); + fielddisplay(self,'underflow_tol','threshold of deep to surface love number ratio to trigger the deletion of layers (default 2.2204460492503131E-016)'); + fielddisplay(self,'integration_steps_per_layer','number of radial steps to propagate the yi system from the bottom to the top of each layer (default 100)'); + fielddisplay(self,'istemporal',{'1 for time-dependent love numbers, 0 for frequency-dependent or elastic love numbers (default 0)', 'If 1: use fourierlove function build_frequencies_from_time to meet consistency'}); + fielddisplay(self,'n_temporal_iterations','max number of iterations in the inverse Laplace transform. Also the number of spectral samples per time step requested (default 8)'); + fielddisplay(self,'time','time vector for deformation if istemporal (default 0) [s]'); fielddisplay(self,'love_kernels','compute love numbers at depth? (default 0)'); - fielddisplay(self,'forcing_type',{'integer indicating the nature and depth of the forcing for the Love number calculation (default 11) :','1: Inner core boundary -- Volumic Potential','2: Inner core boundary -- Pressure','3: Inner core boundary -- Loading','4: Inner core boundary -- Tangential traction','5: Core mantle boundary -- Volumic Potential','6: Core mantle boundary -- Pressure','7: Core mantle boundary -- Loading','8: Core mantle boundary -- Tangential traction','9: Surface -- Volumic Potential','10: Surface -- Pressure','11: Surface -- Loading','12: Surface -- Tangential traction '}); + fielddisplay(self,'forcing_type',{'integer indicating the nature and depth of the forcing for the Love number calculation (default 11) :','1: Inner core boundary -- Volumic Potential','2: Inner core boundary -- Pressure','3: Inner core boundary -- Loading','4: Inner core boundary -- Tangential traction','5: Core mantle boundary -- Volumic Potential','6: Core mantle boundary -- Pressure','7: Core mantle boundary -- Loading','8: Core mantle boundary -- Tangential traction','9: Surface -- Volumic Potential','10: Surface -- Pressure','11: Surface -- Loading','12: Surface -- Tangential traction '}); + fielddisplay(self,'inner_core_boundary','interface index in materials.radius locating forcing. Only used for forcing_type 1--4 (default 1)'); + fielddisplay(self,'core_mantle_boundary','interface index in materials.radius locating forcing. Only used for forcing_type 5--8 (default 2)'); end % }}} @@ -72,8 +96,17 @@ md = checkfield(md,'fieldname','love.r0','NaN',1,'Inf',1,'numel',1,'>',0); md = checkfield(md,'fieldname','love.mu0','NaN',1,'Inf',1,'numel',1,'>',0); + md = checkfield(md,'fieldname','love.Gravitational_Constant','NaN',1,'Inf',1,'numel',1,'>',0); md = checkfield(md,'fieldname','love.allow_layer_deletion','values',[0 1]); + md = checkfield(md,'fieldname','love.underflow_tol','NaN',1,'Inf',1,'numel',1,'>',0); + md = checkfield(md,'fieldname','love.integration_steps_per_layer','NaN',1,'Inf',1,'numel',1,'>',0); md = checkfield(md,'fieldname','love.love_kernels','values',[0 1]); md = checkfield(md,'fieldname','love.forcing_type','NaN',1,'Inf',1,'numel',1,'>',0, '<=', 12); - if md.love.sh_nmin<=1 & md.love.forcing_type==9 + + md = checkfield(md,'fieldname','love.istemporal','values',[0 1]); + if md.love.istemporal==1 + md = checkfield(md,'fieldname','love.n_temporal_iterations','NaN',1,'Inf',1,'numel',1,'>',0); + md = checkfield(md,'fieldname','love.time','NaN',1,'Inf',1,'numel',md.love.nfreq/2/md.love.n_temporal_iterations); + end + if md.love.sh_nmin<=1 & (md.love.forcing_type==9 || md.love.forcing_type==5 || md.love.forcing_type==1) error('Degree 1 not supported for Volumetric Potential forcing. Use sh_min>=2 for this kind of calculation.') end @@ -82,4 +115,11 @@ if ~isa(md.materials,'materials') | ~sum(strcmpi(md.materials.nature,'litho')) error('Need a ''litho'' material to run a Fourier Love number analysis'); + end + + mat=find(strcmpi(md.materials.nature,'litho')); + if md.love.forcing_type<=4 + md = checkfield(md,'fieldname','love.inner_core_boundary','NaN',1,'Inf',1,'numel',1,'>',0, '<=', md.materials(mat).numlayers); + elseif md.love.forcing_type<=8 + md = checkfield(md,'fieldname','love.core_mantle_boundary','NaN',1,'Inf',1,'numel',1,'>',0, '<=', md.materials(mat).numlayers); end @@ -94,7 +134,15 @@ WriteData(fid,prefix,'object',self,'fieldname','r0','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','mu0','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','Gravitational_Constant','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','allow_layer_deletion','format','Boolean'); + WriteData(fid,prefix,'object',self,'fieldname','underflow_tol','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','integration_steps_per_layer','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','istemporal','format','Boolean'); + WriteData(fid,prefix,'object',self,'fieldname','n_temporal_iterations','format','Integer'); + %note: no need to marshall the time vector, we have frequencies WriteData(fid,prefix,'object',self,'fieldname','love_kernels','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','forcing_type','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','inner_core_boundary','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','core_mantle_boundary','format','Integer'); end % }}} @@ -102,4 +150,17 @@ error('not implemented yet!'); end % }}} + function self=build_frequencies_from_time(self) % {{{ + if ~self.istemporal + error('cannot build frequencies for temporal love numbers if love.istemporal==0') + end + disp('Temporal love numbers: Overriding md.love.nfreq and md.love.frequencies'); + self.nfreq=length(self.time)*2*self.n_temporal_iterations; + self.frequencies=zeros(self.nfreq,1); + for i=1:length(self.time) + for j=1:2*self.n_temporal_iterations + self.frequencies((i-1)*2*self.n_temporal_iterations +j) = j*log(2)/self.time(i)/2/pi; + end + end + end % }}} end end Index: /issm/trunk-jpl/src/m/classes/lovenumbers.m =================================================================== --- /issm/trunk-jpl/src/m/classes/lovenumbers.m (revision 26232) +++ /issm/trunk-jpl/src/m/classes/lovenumbers.m (revision 26233) @@ -57,4 +57,5 @@ md = checkfield(md,'fieldname','solidearth.lovenumbers.th','NaN',1,'Inf',1); md = checkfield(md,'fieldname','solidearth.lovenumbers.tk','NaN',1,'Inf',1); + md = checkfield(md,'fieldname','solidearth.lovenumbers.tl','NaN',1,'Inf',1); md = checkfield(md,'fieldname','solidearth.lovenumbers.tk2secular','NaN',1,'Inf',1); @@ -78,4 +79,5 @@ fielddisplay(self,'th','tidal load Love number (deg 2)'); fielddisplay(self,'tk','tidal load Love number (deg 2)'); + fielddisplay(self,'tl','tidal load Love number (deg 2)'); fielddisplay(self,'tk2secular','secular fluid Love number'); Index: /issm/trunk-jpl/src/m/classes/materials.m =================================================================== --- /issm/trunk-jpl/src/m/classes/materials.m (revision 26232) +++ /issm/trunk-jpl/src/m/classes/materials.m (revision 26233) @@ -53,5 +53,9 @@ self.addprop('burgers_viscosity'); self.addprop('burgers_mu'); - self.addprop('isburgers'); + self.addprop('ebm_alpha'); + self.addprop('ebm_delta'); + self.addprop('ebm_taul'); + self.addprop('ebm_tauh'); + self.addprop('rheologymodel'); self.addprop('density'); self.addprop('issolid'); @@ -136,5 +140,10 @@ self.burgers_viscosity=[NaN;NaN]; self.burgers_mu=[NaN;NaN]; - self.isburgers=[0;0]; + + self.ebm_alpha=[NaN;NaN]; + self.ebm_delta=[NaN;NaN]; + self.ebm_taul=[NaN;NaN]; + self.ebm_tauh=[NaN;NaN]; + self.rheologymodel=[0;0]; self.density=[5.51*1e3;5.50*1e3]; % (Pa) %mantle and lithosphere density [kg/m^3] self.issolid=[1;1]; % is layer solid or liquid. @@ -191,5 +200,12 @@ fielddisplay(self,'burgers_viscosity','array describing each layer''s transient viscosity, only for Burgers rheologies (numlayers) [Pa.s]'); fielddisplay(self,'burgers_mu','array describing each layer''s transient shear modulus, only for Burgers rheologies (numlayers) [Pa]'); - fielddisplay(self,'isburgers','array describing whether we adopt a MaxWell (0) or Burgers (1) rheology (default 0)'); + + fielddisplay(self,'ebm_alpha','array describing each layer''s exponent parameter controlling the shape of shear modulus curve between taul and tauh, only for EBM rheology (numlayers)'); + fielddisplay(self,'ebm_delta','array describing each layer''s amplitude of the transient relaxation (ratio between elastic rigity to pre-maxwell relaxation rigity), only for EBM rheology (numlayers)'); + fielddisplay(self,'ebm_taul','array describing each layer''s starting period for transient relaxation, only for EBM rheology (numlayers) [s]'); + fielddisplay(self,'ebm_tauh','array describing each layer''s array describing each layer''s end period for transient relaxation, only for Burgers rheology (numlayers) [s]'); + + + fielddisplay(self,'rheologymodel','array describing whether we adopt a MaxWell (0), Burgers (1) or EBM (2) rheology (default 0)'); fielddisplay(self,'density','array describing each layer''s density (numlayers) [kg/m^3]'); fielddisplay(self,'issolid','array describing whether the layer is solid or liquid (default 1) (numlayers)'); @@ -228,11 +244,18 @@ md = checkfield(md,'fieldname','materials.density','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>',0); md = checkfield(md,'fieldname','materials.viscosity','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>=',0); - md = checkfield(md,'fieldname','materials.isburgers','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>=',0,'<=',1); + md = checkfield(md,'fieldname','materials.rheologymodel','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>=',0,'<=',2); md = checkfield(md,'fieldname','materials.burgers_viscosity','Inf',1,'size',[md.materials.numlayers 1],'>=',0); md = checkfield(md,'fieldname','materials.burgers_mu','Inf',1,'size',[md.materials.numlayers 1],'>=',0); + md = checkfield(md,'fieldname','materials.ebm_alpha','Inf',1,'size',[md.materials.numlayers 1],'>=',0); + md = checkfield(md,'fieldname','materials.ebm_delta','Inf',1,'size',[md.materials.numlayers 1],'>=',0); + md = checkfield(md,'fieldname','materials.ebm_taul','Inf',1,'size',[md.materials.numlayers 1],'>=',0); + md = checkfield(md,'fieldname','materials.ebm_tauh','Inf',1,'size',[md.materials.numlayers 1],'>=',0); for i=1:md.materials.numlayers, - if md.materials.isburgers(i) & (isnan(md.materials.burgers_viscosity(i) | isnan(md.materials.burgers_mu(i)))), - error('materials checkconsistency error message: Litho burgers_viscosity or burgers_mu has NaN values, inconsistent with isburgers choice'); + if md.materials.rheologymodel(i)==1 & (isnan(md.materials.burgers_viscosity(i) | isnan(md.materials.burgers_mu(i)))), + error('materials checkconsistency error message: Litho burgers_viscosity or burgers_mu has NaN values, inconsistent with rheologymodel choice'); + end + if md.materials.rheologymodel(i)==2 & (isnan(md.materials.ebm_alpha(i)) | isnan(md.materials.ebm_delta(i)) | isnan(md.materials.ebm_taul(i)) | isnan(md.materials.ebm_tauh(i)) ), + error('materials checkconsistency error message: Litho ebm_alpha, ebm_delta, ebm_taul or ebm_tauh has NaN values, inconsistent with rheologymodel choice'); end end @@ -290,7 +313,11 @@ WriteData(fid,prefix,'object',self,'class','materials','fieldname','density','format','DoubleMat','mattype',3); WriteData(fid,prefix,'object',self,'class','materials','fieldname','viscosity','format','DoubleMat','mattype',3); - WriteData(fid,prefix,'object',self,'class','materials','fieldname','isburgers','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheologymodel','format','DoubleMat','mattype',3); WriteData(fid,prefix,'object',self,'class','materials','fieldname','burgers_viscosity','format','DoubleMat','mattype',3); WriteData(fid,prefix,'object',self,'class','materials','fieldname','burgers_mu','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',self,'class','materials','fieldname','ebm_alpha','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',self,'class','materials','fieldname','ebm_delta','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',self,'class','materials','fieldname','ebm_taul','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',self,'class','materials','fieldname','ebm_tauh','format','DoubleMat','mattype',3); %compute earth density compatible with our layer density distribution: earth_density=0; @@ -351,7 +378,11 @@ writejsdouble(fid,[modelname '.materials.density'],self.density); writejsdouble(fid,[modelname '.materials.viscosity'],self.viscosity); - writejsdouble(fid,[modelname '.materials.isburgers'],self.isburgers); + writejsdouble(fid,[modelname '.materials.rheologymodel'],self.rheologymodel); writejsdouble(fid,[modelname '.materials.burgers_viscosity'],self.burgers_viscosity); writejsdouble(fid,[modelname '.materials.burgers_mu'],self.burgers_mu); + writejsdouble(fid,[modelname '.materials.ebm_alpha'],self.ebm_alpha); + writejsdouble(fid,[modelname '.materials.ebm_delta'],self.ebm_delta); + writejsdouble(fid,[modelname '.materials.ebm_taul'],self.ebm_taul); + writejsdouble(fid,[modelname '.materials.ebm_tauh'],self.ebm_tauh); case 'hydro'