Index: /issm/trunk-jpl/src/m/classes/flowequation.js =================================================================== --- /issm/trunk-jpl/src/m/classes/flowequation.js (revision 26304) +++ /issm/trunk-jpl/src/m/classes/flowequation.js (revision 26305) @@ -142,5 +142,5 @@ this.isFS = 0; this.isNitscheBC = 0; - this.FSNitscheGamma = 0; + this.FSNitscheGamma = 1e6; this.fe_SSA = ''; this.fe_HO = ''; Index: /issm/trunk-jpl/src/m/classes/fourierlove.js =================================================================== --- /issm/trunk-jpl/src/m/classes/fourierlove.js (revision 26304) +++ /issm/trunk-jpl/src/m/classes/fourierlove.js (revision 26305) @@ -18,5 +18,5 @@ this.g0=9.81; // m/s^2; this.r0=6371*1e3; //m; - this.mu0=10^11; // Pa + this.mu0=1e11; // Pa this.Gravitational_Constant=6.67259e-11; // m^3 kg^-1 s^-2 this.allow_layer_deletion=1; @@ -86,5 +86,5 @@ /* if ~isa(md.materials,'materials') | ~sum(strcmpi(md.materials.nature,'litho')) - error('Need a ''litho'' material to run a Fourier Love number analysis'); + error('Need a "litho" material to run a Fourier Love number analysis'); end Index: /issm/trunk-jpl/src/m/classes/friction.js =================================================================== --- /issm/trunk-jpl/src/m/classes/friction.js (revision 26304) +++ /issm/trunk-jpl/src/m/classes/friction.js (revision 26305) @@ -67,5 +67,5 @@ WriteData(fid,prefix,'name','md.friction.law','data',1,'format','Integer'); let mattype,tsl; - if ((size(this.coefficient,1)==md.mesh.numberofvertices | size(this.coefficient,1)==md.mesh.numberofvertices+1)) { + if ((size(this.coefficient,0)==md.mesh.numberofvertices | size(this.coefficient,1)==md.mesh.numberofvertices+1)) { mattype=1; tsl = md.mesh.numberofvertices; @@ -81,4 +81,5 @@ switch (this.coupling) { case 0: + break; case 1: WriteData(fid,prefix,'class','friction','object',this,'fieldname','effective_pressure','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); Index: /issm/trunk-jpl/src/m/classes/frontalforcings.js =================================================================== --- /issm/trunk-jpl/src/m/classes/frontalforcings.js (revision 26305) +++ /issm/trunk-jpl/src/m/classes/frontalforcings.js (revision 26305) @@ -0,0 +1,41 @@ +//FRONTAL FORCINGS Class definition +// +// Usage: +// frontalforcings=frontalforcings(); + +function frontalforcings(){ + //methods + this.classname = function(){ // {{{ + return "frontalforcings"; + } // }}} + this.extrude = function(md) {//{{{ + this.meltingrate=project3d(md,'vector',this.meltingrate,'type','node'); + return this; + }//}}} + this.setdefaultparameters = function(){// {{{ + this.meltingrate = NaN; + } // }}} + this.checkconsistency = function(md,solution,analyses) { //{{{ + //Early return + if (!solution=='TransientSolution' || md.transient.ismovingfront == 0) return; + + md = checkfield(md,'fieldname','frontalforcings.meltingrate','NaN',1,'Inf',1,'timeseries',1,'>=',0); + + } // }}} + this.disp = function(){ // {{{ + console.log(sprintf(' Frontalforcings parameters:')); + fielddisplay(this,'meltingrate','melting rate at given location [m/a]'); + } // }}} + this.marshall=function(md,prefix,fid) { //{{{ + let yts=md.constants.yts; + WriteData(fid,prefix,'name','md.frontalforcings.parameterization','data',1,'format','Integer'); + WriteData(fid,prefix,'object',this,'fieldname','meltingrate','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1./yts); + }//}}} + this.fix=function() { //{{{ + }//}}} + //properties + // {{{ + this.meltingrate = NaN; + this.setdefaultparameters(); + // }}} +} Index: /issm/trunk-jpl/src/m/classes/geometry.js =================================================================== --- /issm/trunk-jpl/src/m/classes/geometry.js (revision 26304) +++ /issm/trunk-jpl/src/m/classes/geometry.js (revision 26305) @@ -48,5 +48,5 @@ } // }}} this.marshall=function(md,prefix,fid) { //{{{ - let length_thickness=size(this.thickness,1); + let length_thickness=size(this.thickness,0); if (length_thickness==md.mesh.numberofvertices || length_thickness==md.mesh.numberofvertices+1) { WriteData(fid,prefix,'object',this,'fieldname','thickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); Index: /issm/trunk-jpl/src/m/classes/groundingline.js =================================================================== --- /issm/trunk-jpl/src/m/classes/groundingline.js (revision 26304) +++ /issm/trunk-jpl/src/m/classes/groundingline.js (revision 26305) @@ -47,6 +47,6 @@ this.marshall=function(md,prefix,fid) { //{{{ WriteData(fid,prefix,'data',this.migration,'name','md.groundingline.migration','format','String'); - WriteData(fid,prefix,'data',this.migration,'name','md.groundingline.friction_interpolation','format','String'); - WriteData(fid,prefix,'data',this.migration,'name','md.groundingline.melt_interpolation','format','String'); + WriteData(fid,prefix,'data',this.friction_interpolation,'name','md.groundingline.friction_interpolation','format','String'); + WriteData(fid,prefix,'data',this.melt_interpolation,'name','md.groundingline.melt_interpolation','format','String'); }//}}} this.fix=function() { //{{{ Index: /issm/trunk-jpl/src/m/classes/levelset.js =================================================================== --- /issm/trunk-jpl/src/m/classes/levelset.js (revision 26304) +++ /issm/trunk-jpl/src/m/classes/levelset.js (revision 26305) @@ -8,6 +8,6 @@ this.setdefaultparameters = function(){// {{{ - //stabilization = 2 by default - this.stabilization = 2; + //stabilization = 1 by default + this.stabilization = 1; this.reinit_frequency = 5; this.kill_icebergs = 1; Index: /issm/trunk-jpl/src/m/classes/materials.js =================================================================== --- /issm/trunk-jpl/src/m/classes/materials.js (revision 26305) +++ /issm/trunk-jpl/src/m/classes/materials.js (revision 26305) @@ -0,0 +1,389 @@ +//MATERIALS class definition +// +// Usage: +// materials=materials(); + +class materials { + //methods + constructor() {//{{{ + this.nature = []; + + let nargs = arguments.length; + if (nargs == 0) { + this.nature=['ice']; + } else { + this.nature=arguments; + } + + //check this is acceptable: + for (let i=0; i < length(this.nature); i++){ + if (!(strcmpi(this.nature[i],'litho') || strcmpi(this.nature[i],'ice') || strcmpi(this.nature[i],'hydro'))) { + error('materials constructor error message: nature of the material not supported yet! ("ice" or "litho" or "hydro")'); + } + } + + //start filling in the dynamic fields: + for (let i=0; i < length(this.nature); i++) { + let nat=this.nature[i]; + switch (nat) { + case 'ice': + this.rho_ice = 0; + this.rho_water = 0; + this.rho_freshwater = 0; + this.mu_water = 0; + this.heatcapacity = 0; + this.latentheat = 0; + this.thermalconductivity = 0; + this.temperateiceconductivity = 0; + this.effectiveconductivity_averaging = 0; + this.meltingpoint = 0; + this.beta = 0; + this.mixed_layer_capacity = 0; + this.thermal_exchange_velocity = 0; + this.rheology_B = 0; + this.rheology_n = 0; + this.rheology_law = 0; + break; + case 'litho': + this.numlayers = 0; + this.radius = 0; + this.viscosity = 0; + this.lame_lambda = 0; + this.lame_mu = 0; + this.burgers_viscosity = 0; + this.burgers_mu = 0; + this.ebm_alpha = 0; + this.ebm_delta = 0; + this.ebm_taul = 0; + this.ebm_tauh = 0; + this.rheologymodel = 0; + this.density = 0; + this.issolid = 0; + break; + case 'hydro': + this.rho_ice = 0; + this.rho_water = 0; + this.rho_freshwater = 0; + break; + default: + error('materials constructor error message: nature of the material not supported yet! ("ice" or "litho" or "hydro")'); + } + } + } //}}} + setdefaultparameters() {//{{{ + for (let i=0; i < length(this.nature); i++) { + let nat=this.nature[i]; + switch (nat) { + case 'ice': + //ice density (kg/m^3) + this.rho_ice=917.; + + //ocean water density (kg/m^3) + this.rho_water=1023.; + + //fresh water density (kg/m^3) + this.rho_freshwater=1000.; + + //water viscosity (N.s/m^2) + this.mu_water=0.001787; + + //ice heat capacity cp (J/kg/K) + this.heatcapacity=2093.; + + //ice latent heat of fusion L (J/kg) + this.latentheat=3.34*1e5; + + //ice thermal conductivity (W/m/K) + this.thermalconductivity=2.4; + + //wet ice thermal conductivity (W/m/K) + this.temperateiceconductivity=.24; + + //computation of effective conductivity + this.effectiveconductivity_averaging=1; + + //the melting point of ice at 1 atmosphere of pressure in K + this.meltingpoint=273.15; + + //rate of change of melting point with pressure (K/Pa) + this.beta=9.8*1e-8; + + //mixed layer (ice-water interface) heat capacity (J/kg/K) + this.mixed_layer_capacity=3974.; + + //thermal exchange velocity (ice-water interface) (m/s) + this.thermal_exchange_velocity=1.00*1e-4; + + //Rheology law: what is the temperature dependence of B with T + //available: none, paterson and arrhenius + this.rheology_law='Paterson'; + + //Rheology fields default: + this.rheology_B = 1 * 1e8; + this.rheology_n = 3; + break; + + case 'litho': + //we default to a configuration that enables running GIA solutions using giacaron and/or giaivins. + this.numlayers=2; + + //center of the earth (approximation, must not be 0), then the lab (lithosphere/asthenosphere boundary) then the surface + //(with 1d3 to avoid numerical singularities) + this.radius=[1e3,6278*1e3,6378*1e3]; + + this.viscosity=[1e21,1e40]; //mantle and lithosphere viscosity (respectively) [Pa.s] + this.lame_mu=[1.45*1e11,6.7*1e10]; // (Pa) //lithosphere and mantle shear modulus (respectively) [Pa] + this.lame_lambda=this.lame_mu; // (Pa) //mantle and lithosphere lamba parameter (respectively) [Pa] + this.burgers_viscosity=[NaN,NaN]; + this.burgers_mu=[NaN,NaN]; + + this.ebm_alpha=[NaN,NaN]; + this.ebm_delta=[NaN,NaN]; + this.ebm_taul=[NaN,NaN]; + this.ebm_tauh=[NaN,NaN]; + this.rheologymodel=[0,0]; + this.density=[5.51*1e3,5.50*1e3]; // (Pa) //mantle and lithosphere density [kg/m^3] + this.issolid=[1,1]; // is layer solid or liquid. + break; + + case 'hydro': + //ice density (kg/m^3) + this.rho_ice=917.; + + //ocean water density (kg/m^3) + this.rho_water=1023.; + + //fresh water density (kg/m^3) + this.rho_freshwater=1000.; + break; + + default: + error('materials constructor error message: nature of the material not supported yet! ("ice" or "litho" or "hydro")'); + } + + // average density of the Earth (kg/m^3) + this.earth_density=5512; + + } + } //}}} + disp() {//{{{ + console.log(sprintf(' Materials:')); + + for (let i=0; i < length(this.nature); i++) { + let nat=this.nature[i]; + switch (nat) { + case 'ice': + console.log(sprintf(' \nIce:')); + fielddisplay(this,'rho_ice','ice density [kg/m^3]'); + fielddisplay(this,'rho_water','ocean water density [kg/m^3]'); + fielddisplay(this,'rho_freshwater','fresh water density [kg/m^3]'); + fielddisplay(this,'mu_water','water viscosity [N s/m^2]'); + fielddisplay(this,'heatcapacity','heat capacity [J/kg/K]'); + fielddisplay(this,'thermalconductivity','ice thermal conductivity [W/m/K]'); + fielddisplay(this,'temperateiceconductivity','temperate ice thermal conductivity [W/m/K]'); + fielddisplay(this,'meltingpoint','melting point of ice at 1atm in K'); + fielddisplay(this,'latentheat','latent heat of fusion [J/kg]'); + fielddisplay(this,'beta','rate of change of melting point with pressure [K/Pa]'); + fielddisplay(this,'mixed_layer_capacity','mixed layer capacity [W/kg/K]'); + fielddisplay(this,'thermal_exchange_velocity','thermal exchange velocity [m/s]'); + fielddisplay(this,'rheology_B','flow law parameter [Pa s^(1/n)]'); + fielddisplay(this,'rheology_n','Glen\'s flow law exponent'); + fielddisplay(this,'rheology_law','law for the temperature dependance of the rheology: "None", "BuddJacka", Cuffey", "CuffeyTemperate", "Paterson", "Arrhenius", "LliboutryDuval", "NyeCO2", or "NyeH2O"'); + break; + case 'litho': + console.log(sprintf(' \nLitho:')); + fielddisplay(this,'numlayers','number of layers (default: 2)'); + fielddisplay(this,'radius','array describing the radius for each interface (numlayers+1) [m]'); + fielddisplay(this,'viscosity','array describing each layer\'s viscosity (numlayers) [Pa.s]'); + fielddisplay(this,'lame_lambda','array describing the lame lambda parameter (numlayers) [Pa]'); + fielddisplay(this,'lame_mu','array describing the shear modulus for each layers (numlayers) [Pa]'); + fielddisplay(this,'burgers_viscosity','array describing each layer\'s transient viscosity, only for Burgers rheologies (numlayers) [Pa.s]'); + fielddisplay(this,'burgers_mu','array describing each layer\'s transient shear modulus, only for Burgers rheologies (numlayers) [Pa]'); + + fielddisplay(this,'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(this,'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(this,'ebm_taul','array describing each layer\'s starting period for transient relaxation, only for EBM rheology (numlayers) [s]'); + fielddisplay(this,'ebm_tauh','array describing each layer\'s array describing each layer\'s end period for transient relaxation, only for Burgers rheology (numlayers) [s]'); + + + fielddisplay(this,'rheologymodel','array describing whether we adopt a Maxwell (0), Burgers (1) or EBM (2) rheology (default: 0)'); + fielddisplay(this,'density','array describing each layer\'s density (numlayers) [kg/m^3]'); + fielddisplay(this,'issolid','array describing whether the layer is solid or liquid (default 1) (numlayers)'); + break; + case 'hydro': + console.log(sprintf(' \nHydro:')); + fielddisplay(this,'rho_ice','ice density [kg/m^3]'); + fielddisplay(this,'rho_water','ocean water density [kg/m^3]'); + fielddisplay(this,'earth_density','mantle density [kg/m^3]'); + fielddisplay(this,'rho_freshwater','fresh water density [kg/m^3]'); + break; + default: + error('materials constructor error message: nature of the material not supported yet! ("ice" or "litho" or "hydro")'); + } + } + + } // }}} + extrude() {//{{{ + for (let i=0; i < length(this.nature); i++) { + let nat=this.nature[i]; + switch (nat) { + case 'ice': + this.rheology_B=project3d(md,'vector',this.rheology_B,'type','node'); + this.rheology_n=project3d(md,'vector',this.rheology_n,'type','element'); + break; + } + } + return this; + }//}}} + checkconsistency(md,solution,analyses) { // {{{ + + for (let i=0; i < length(this.nature); i++) { + let nat=this.nature[i]; + switch (nat) { + case 'ice': + checkfield(md,'fieldname','materials.rho_ice','>',0); + checkfield(md,'fieldname','materials.rho_water','>',0); + checkfield(md,'fieldname','materials.rho_freshwater','>',0); + checkfield(md,'fieldname','materials.mu_water','>',0); + checkfield(md,'fieldname','materials.rheology_B','>',0,'timeseries',1,'NaN',1,'Inf',1); + checkfield(md,'fieldname','materials.rheology_n','>',0,'size',[md.mesh.numberofelements, 1]); + checkfield(md,'fieldname','materials.rheology_law','values',['None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval', 'NyeCO2', 'NyeH2O']); + break; + case 'litho': + if (!ismember('LoveAnalysis',analyses)) return; end + checkfield(md,'fieldname','materials.numlayers','NaN',1,'Inf',1,'>',0,'numel',1); + checkfield(md,'fieldname','materials.radius','NaN',1,'Inf',1,'size',[md.materials.numlayers+1, 1],'>',0); + checkfield(md,'fieldname','materials.lame_mu','NaN',1,'Inf',1,'size',[md.materials.numlayers, 1],'>=',0); + checkfield(md,'fieldname','materials.lame_lambda','NaN',1,'Inf',1,'size',[md.materials.numlayers, 1],'>=',0); + checkfield(md,'fieldname','materials.issolid','NaN',1,'Inf',1,'size',[md.materials.numlayers, 1],'>=',0,'<',2); + checkfield(md,'fieldname','materials.density','NaN',1,'Inf',1,'size',[md.materials.numlayers, 1],'>',0); + checkfield(md,'fieldname','materials.viscosity','NaN',1,'Inf',1,'size',[md.materials.numlayers, 1],'>=',0); + checkfield(md,'fieldname','materials.rheologymodel','NaN',1,'Inf',1,'size',[md.materials.numlayers, 1],'>=',0,'<=',2); + checkfield(md,'fieldname','materials.burgers_viscosity','Inf',1,'size',[md.materials.numlayers, 1],'>=',0); + checkfield(md,'fieldname','materials.burgers_mu','Inf',1,'size',[md.materials.numlayers, 1],'>=',0); + checkfield(md,'fieldname','materials.ebm_alpha','Inf',1,'size',[md.materials.numlayers, 1],'>=',0); + checkfield(md,'fieldname','materials.ebm_delta','Inf',1,'size',[md.materials.numlayers, 1],'>=',0); + checkfield(md,'fieldname','materials.ebm_taul','Inf',1,'size',[md.materials.numlayers, 1],'>=',0); + checkfield(md,'fieldname','materials.ebm_tauh','Inf',1,'size',[md.materials.numlayers, 1],'>=',0); + + for (let i=0; i < md.materials.numlayers; i++) { + 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'); + } + 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'); + } + } + if (md.materials.issolid[0]==0 || md.materials.lame_mu[0]==0) { + error('First layer must be solid (issolid(1) > 0 AND lame_mu(1) > 0). Add a weak inner core if necessary.'); + } + ind=find(ArrayEqual(md.materials.issolid,0)); //ind=find(md.materials.issolid==0); + console.log('materials consistency warning: not checking for fluid layers (not implemented)'); + //if sum(ismember(diff(ind),1)>=1) %if there are at least two consecutive indices that contain issolid = 0 + // error(['Fluid layers detected at layers #', num2str(ind'), ', but having 2 or more adjacent fluid layers is not supported yet. Consider merging them.']) + //} + break; + case 'hydro': + checkfield(md,'fieldname','materials.rho_ice','>',0); + checkfield(md,'fieldname','materials.rho_water','>',0); + checkfield(md,'fieldname','materials.earth_density','>',0,'numel',1); + checkfield(md,'fieldname','materials.rho_freshwater','>',0); + break; + default: + error('materials constructor error message: nature of the material not supported yet! ("ice" or "litho" or "hydro")'); + } + } + + } // }}} + marshall(md,prefix,fid) { //{{{ + //1: MatdamageiceEnum 2: MatestarEnum 3: MaticeEnum 4: MatenhancediceEnum 5: MaterialsEnum + WriteData(fid,prefix,'name','md.materials.nature','data',naturetointeger(this.nature),'format','IntMat','mattype',3); + WriteData(fid,prefix,'name','md.materials.type','data',5,'format','Integer'); //DANGER: this can evolve if you have classes. + for (let i=0; i < length(this.nature); i++) { + let nat=this.nature[i]; + switch (nat) { + case 'ice': + WriteData(fid,prefix,'object',this,'class','materials','fieldname','rho_ice','format','Double'); + WriteData(fid,prefix,'object',this,'class','materials','fieldname','rho_water','format','Double'); + WriteData(fid,prefix,'object',this,'class','materials','fieldname','rho_freshwater','format','Double'); + WriteData(fid,prefix,'object',this,'class','materials','fieldname','mu_water','format','Double'); + WriteData(fid,prefix,'object',this,'class','materials','fieldname','heatcapacity','format','Double'); + WriteData(fid,prefix,'object',this,'class','materials','fieldname','latentheat','format','Double'); + WriteData(fid,prefix,'object',this,'class','materials','fieldname','thermalconductivity','format','Double'); + WriteData(fid,prefix,'object',this,'class','materials','fieldname','temperateiceconductivity','format','Double'); + WriteData(fid,prefix,'object',this,'class','materials','fieldname','effectiveconductivity_averaging','format','Integer'); + WriteData(fid,prefix,'object',this,'class','materials','fieldname','meltingpoint','format','Double'); + WriteData(fid,prefix,'object',this,'class','materials','fieldname','beta','format','Double'); + WriteData(fid,prefix,'object',this,'class','materials','fieldname','mixed_layer_capacity','format','Double'); + WriteData(fid,prefix,'object',this,'class','materials','fieldname','thermal_exchange_velocity','format','Double'); + WriteData(fid,prefix,'object',this,'class','materials','fieldname','rheology_B','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',this,'class','materials','fieldname','rheology_n','format','DoubleMat','mattype',2); + WriteData(fid,prefix,'data',this.rheology_law,'name','md.materials.rheology_law','format','String'); + break; + case 'litho': + WriteData(fid,prefix,'object',this,'class','materials','fieldname','numlayers','format','Integer'); + WriteData(fid,prefix,'object',this,'class','materials','fieldname','radius','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',this,'class','materials','fieldname','lame_mu','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',this,'class','materials','fieldname','lame_lambda','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',this,'class','materials','fieldname','issolid','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',this,'class','materials','fieldname','density','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',this,'class','materials','fieldname','viscosity','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',this,'class','materials','fieldname','rheologymodel','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',this,'class','materials','fieldname','burgers_viscosity','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',this,'class','materials','fieldname','burgers_mu','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',this,'class','materials','fieldname','ebm_alpha','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',this,'class','materials','fieldname','ebm_delta','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',this,'class','materials','fieldname','ebm_taul','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',this,'class','materials','fieldname','ebm_tauh','format','DoubleMat','mattype',3); + //compute earth density compatible with our layer density distribution: + let earth_density=0; + for (let i=0; i < this.numlayers; i++){ + earth_density=earth_density + (this.radius[i+1]**3-this.radius[i]**3)*this.density[i]; + } + earth_density=earth_density/this.radius(this.numlayers+1)**3; + console.log(earth_density); + this.earth_density=earth_density; + break; + case 'hydro': + WriteData(fid,prefix,'object',this,'class','materials','fieldname','rho_ice','format','Double'); + WriteData(fid,prefix,'object',this,'class','materials','fieldname','rho_water','format','Double'); + WriteData(fid,prefix,'object',this,'class','materials','fieldname','rho_freshwater','format','Double'); + break; + default: + error('materials constructor error message: nature of the material not supported yet! ("ice" or "litho" or "hydro")'); + } + } + WriteData(fid,prefix,'data',this.earth_density,'name','md.materials.earth_density','format','Double'); + }//}}} +} +function naturetointeger(strnat) {// {{{ + let intnat=zeros(length(strnat),1); + for (let i=0; i < length(strnat); i++){ + switch (strnat[i]) { + case 'damageice': + intnat[i]=1; + break; + case 'estar': + intnat[i]=2; + break; + case 'ice': + intnat[i]=3; + break; + case 'enhancedice': + intnat[i]=4; + break; + //case 'materials': //this case will never happen, kept to ensure equivalent of codes between IoCodeToMaterialsEnum and IoCodeToNatureEnum + // intnat[i]=5; + // break; + case 'litho': + intnat[i]=6; + break; + case 'hydro': + intnat[i]=7; + break; + default: + error('materials constructor error message: nature of the material not supported yet! ("ice" or "litho" or "hydro")'); + } + } + return intnat; +}// }}} Index: /issm/trunk-jpl/src/m/classes/mesh3dsurface.js =================================================================== --- /issm/trunk-jpl/src/m/classes/mesh3dsurface.js (revision 26304) +++ /issm/trunk-jpl/src/m/classes/mesh3dsurface.js (revision 26305) @@ -87,4 +87,5 @@ WriteData(fid,prefix,'object',this,'fieldname','long','format','DoubleMat','mattype',1); WriteData(fid,prefix,'object',this,'fieldname','r','format','DoubleMat','mattype',1); + //WriteData(fid,prefix,'object',this,'fieldname','z','format','DoubleMat','mattype',1); WriteData(fid,prefix,'object',this,'fieldname','elements','format','DoubleMat','mattype',2); WriteData(fid,prefix,'object',this,'fieldname','numberofelements','format','Integer'); @@ -92,5 +93,5 @@ WriteData(fid,prefix,'object',this,'fieldname','average_vertex_connectivity','format','Integer'); WriteData(fid,prefix,'object',this,'fieldname','vertexonboundary','format','DoubleMat','mattype',1); - WriteData(fid,prefix,'object',this,'fieldname','segments','format','DoubleMat','mattype',3); + //WriteData(fid,prefix,'object',this,'fieldname','segments','format','DoubleMat','mattype',3); }//}}} this.fix=function() { //{{{ Index: /issm/trunk-jpl/src/m/classes/model.js =================================================================== --- /issm/trunk-jpl/src/m/classes/model.js (revision 26304) +++ /issm/trunk-jpl/src/m/classes/model.js (revision 26305) @@ -86,4 +86,5 @@ this.calving = new calving(); this.gia = new giaivins(); + this.frontforcings = new frontforcings(); this.love = new fourierlove(); this.esa = new esa(); @@ -759,4 +760,5 @@ this.levelset = 0; this.calving = 0; + this.frontforcings = 0; this.love = 0; this.gia = 0; Index: /issm/trunk-jpl/src/m/classes/sampling.js =================================================================== --- /issm/trunk-jpl/src/m/classes/sampling.js (revision 26305) +++ /issm/trunk-jpl/src/m/classes/sampling.js (revision 26305) @@ -0,0 +1,99 @@ +//SAMPLING class definition +// +// Usage: +// sampling=new sampling(); + +function sampling (){ + //methods + this.defaultoutputs = function(){// {{{ + return []; + }// }}} + this.setdefaultparameters = function(){// {{{ + + //Scaling coefficient + this.tau=1; + + //Apply Robin boundary conditions + this.robin=0; + + //Temporal correlation factor + this.phi=0; + + //Exponent in fraction SPDE (default=2, biLaplacian covariance + //operator) + this.alpha=2; // Default + + //Seed for pseudorandom number generator (default -1 for random seed) + this.seed=-1; + + //default output + this.requested_outputs=['default']; + }// }}} + this.disp= function(){// {{{ + console.log(sprintf(' Sampling parameters:')); + + console.log(sprintf('\n %s','Parameters of PDE operator (kappa^2 I-Laplacian)^(alpha/2)(tau):')); + fielddisplay(this,'kappa','coefficient of the identity operator'); + fielddisplay(this,'tau','scaling coefficient of the solution (default 1.0)'); + fielddisplay(this,'alpha','exponent in PDE operator, (default 2.0, BiLaplacian covariance operator)'); + + console.log(sprintf('\n %s','Parameters of Robin boundary conditions nabla () \cdot normvec + beta ():')); + fielddisplay(this,'robin','Apply Robin boundary conditions (1 if applied and 0 for homogenous Neumann boundary conditions) (default 0)'); + fielddisplay(this,'beta','Coefficient in Robin boundary conditions (to be defined for robin = 1)'); + + console.log(sprintf('\n %s','Parameters for first-order autoregressive process (X_t = phi X_{t-1} + noise) (if transient):')); + fielddisplay(this,'phi','Temporal correlation factor (|phi|<1 for stationary process, phi = 1 for random walk process) (default 0)'); + + console.log(sprintf('\n %s','Other parameters of stochastic sampler:')); + fielddisplay(this,'seed','Seed for pseudorandom number generator (given seed if >=0 and random seed if <0) (default -1)'); + fielddisplay(this,'requested_outputs','additional outputs requested (not implemented yet)'); + }// }}} + this.classname= function(){// {{{ + return "sampling"; + }// }}} + this.checkconsistency = function(md,solution,analyses) { //{{{ + if (!ismember('SamplingAnalysis',analyses)) return; + + checkfield(md,'fieldname','sampling.kappa','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices,1],'>',0); + checkfield(md,'fieldname','sampling.tau','NaN',1,'Inf',1,'numel',1,'>',0); + checkfield(md,'fieldname','sampling.robin','numel',1,'values',[0,1]); + if (md.sampling.robin) { + checkfield(md,'fieldname','sampling.beta','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices,1],'>',0); + } + checkfield(md,'fieldname','sampling.phi','NaN',1,'Inf',1,'numel',1,'>=',0); + checkfield(md,'fieldname','sampling.alpha','NaN',1,'Inf',1,'numel',1,'>',0); + checkfield(md,'fieldname','sampling.seed','NaN',1,'Inf',1,'numel',1); + checkfield(md,'fieldname','sampling.requested_outputs','stringrow',1); + } // }}} + this.marshall=function(md,prefix,fid) { //{{{ + WriteData(fid,prefix,'object',this,'fieldname','kappa','format','DoubleMat','mattype',1); + WriteData(fid,prefix,'object',this,'fieldname','tau','format','Double'); + WriteData(fid,prefix,'object',this,'fieldname','beta','format','DoubleMat','mattype',1); + WriteData(fid,prefix,'object',this,'fieldname','phi','format','Double'); + WriteData(fid,prefix,'object',this,'fieldname','alpha','format','Integer'); + WriteData(fid,prefix,'object',this,'fieldname','robin','format','Boolean'); + WriteData(fid,prefix,'object',this,'fieldname','seed','format','Integer'); + + //process requested outputs + outputs = this.requested_outputs; + pos = find(ismember(outputs,'default')); + if (!isempty(pos)) { + ArrayIndex(outputs,pos,[]); //remove 'default' from outputs + outputs = ArrayConcat(outputs, this.defaultoutputs()); //add defaults + } + WriteData(fid,prefix,'data',outputs,'name','md.sampling.requested_outputs','format','StringArray'); + }//}}} + //properties + // {{{ + this.kappa = NaN; + this.tau = 0; + this.beta = NaN; + this.phi = 0; + this.alpha = 0; + this.robin = 0; + this.seed = 0; + this.requested_outputs = []; + + this.setdefaultparameters(); + //}}} +} Index: sm/trunk-jpl/src/m/classes/trans.js =================================================================== --- /issm/trunk-jpl/src/m/classes/trans.js (revision 26304) +++ (revision ) @@ -1,132 +1,0 @@ -//TRANS class definition -// -// Usage: -// trans =new trans(); - -function trans (){ - //methods - this.setdefaultparameters = function(){// {{{ - - //full analysis: Stressbalance, Masstransport and Thermal but no groundingline migration for now - this.issmb = 1; - this.ismasstransport = 1; - this.isstressbalance = 1; - this.isthermal = 0; - this.isgroundingline = 0; - this.isgia = 0; - this.isesa = 0; - this.isdamageevolution = 0; - this.ismovingfront = 0; - this.ishydrology = 0; - this.isslr = 0; - this.isoceancoupling = 0; - this.iscoupler = 0; - this.amr_frequency = 0; - - //default output - this.requested_outputs=['default']; - - }// }}} - this.disp= function(){// {{{ - - console.log(sprintf(' transient solution parameters:')); - - fielddisplay(this,'issmb','indicates whether a surface mass balance solution is used in the transient'); - fielddisplay(this,'ismasstransport','indicates whether a masstransport solution is used in the transient'); - fielddisplay(this,'isstressbalance','indicates whether a stressbalance solution is used in the transient'); - fielddisplay(this,'isthermal','indicates whether a thermal solution is used in the transient'); - fielddisplay(this,'isgroundingline','indicates whether a groundingline migration is used in the transient'); - fielddisplay(this,'isgia','indicates whether a postglacial rebound model is used in the transient'); - fielddisplay(this,'isesa','indicates whether an elastic adjustment model is used in the transient'); - fielddisplay(this,'isdamageevolution','indicates whether damage evolution is used in the transient'); - fielddisplay(this,'ismovingfront','indicates whether a moving front capability is used in the transient'); - fielddisplay(this,'ishydrology','indicates whether an hydrology model is used'); - fielddisplay(this,'isslr','indicates whether a sea-level rise model is used'); - fielddisplay(this,'isoceancoupling','indicates whether a coupling with an ocean model is used in the transient'); - fielddisplay(this,'iscoupler','indicates whether different models are being run with need for coupling'); - fielddisplay(this,'amr_frequency','frequency at which mesh is refined in simulations with multiple time_steps'); - fielddisplay(this,'requested_outputs','list of additional outputs requested'); - - - }// }}} - this.classname= function(){// {{{ - return "transient"; - }// }}} - this.checkconsistency = function(md,solution,analyses) { // {{{ - - //Early return - if (solution!='TransientSolution') return; - - checkfield(md,'fieldname','trans.issmb','numel',[1],'values',[0, 1]); - checkfield(md,'fieldname','trans.ismasstransport','numel',[1],'values',[0, 1]); - checkfield(md,'fieldname','trans.isstressbalance','numel',[1],'values',[0, 1]); - checkfield(md,'fieldname','trans.isthermal','numel',[1],'values',[0, 1]); - checkfield(md,'fieldname','trans.isgroundingline','numel',[1],'values',[0, 1]); - checkfield(md,'fieldname','trans.isgia','numel',[1],'values',[0, 1]); - checkfield(md,'fieldname','trans.isesa','numel',[1],'values',[0, 1]); - checkfield(md,'fieldname','trans.isdamageevolution','numel',[1],'values',[0, 1]); - checkfield(md,'fieldname','trans.ismovingfront','numel',[1],'values',[0, 1]); - checkfield(md,'fieldname','trans.ishydrology','numel',[1],'values',[0, 1]); - checkfield(md,'fieldname','trans.isslr','numel',[1],'values',[0, 1]); - checkfield(md,'fieldname','trans.isoceancoupling','numel',[1],'values',[0, 1]); - checkfield(md,'fieldname','trans.iscoupler','numel',[1],'values',[0, 1]); - checkfield(md,'fieldname','trans.amr_frequency','numel',[1],'>=',0,'NaN',1,'Inf',1); - checkfield(md,'fieldname','trans.requested_outputs','stringrow',1); - } // }}} - this.marshall=function(md,prefix,fid) { //{{{ - - prefix='md.transient'; - WriteData(fid,prefix,'object',this,'fieldname','issmb','format','Boolean'); - WriteData(fid,prefix,'object',this,'fieldname','ismasstransport','format','Boolean'); - WriteData(fid,prefix,'object',this,'fieldname','isoceantransport','format','Boolean'); - WriteData(fid,prefix,'object',this,'fieldname','isstressbalance','format','Boolean'); - WriteData(fid,prefix,'object',this,'fieldname','isthermal','format','Boolean'); - WriteData(fid,prefix,'object',this,'fieldname','isgroundingline','format','Boolean'); - WriteData(fid,prefix,'object',this,'fieldname','isesa','format','Boolean'); - WriteData(fid,prefix,'object',this,'fieldname','isdamageevolution','format','Boolean'); - WriteData(fid,prefix,'object',this,'fieldname','ishydrology','format','Boolean'); - WriteData(fid,prefix,'object',this,'fieldname','ismovingfront','format','Boolean'); - WriteData(fid,prefix,'object',this,'fieldname','issampling','format','Boolean'); - WriteData(fid,prefix,'object',this,'fieldname','isslc','format','Boolean'); - WriteData(fid,prefix,'object',this,'fieldname','isoceancoupling','format','Boolean'); - WriteData(fid,prefix,'object',this,'fieldname','amr_frequency','format','Integer'); - - //process requested outputs - var outputs = this.requested_outputs; - for (var i=0;i=',0,'NaN',1,'Inf',1); + checkfield(md,'fieldname','transient.requested_outputs','stringrow',1); + } // }}} + this.marshall=function(md,prefix,fid) { //{{{ + + prefix='md.transient'; + WriteData(fid,prefix,'object',this,'fieldname','issmb','format','Boolean'); + WriteData(fid,prefix,'object',this,'fieldname','ismasstransport','format','Boolean'); + WriteData(fid,prefix,'object',this,'fieldname','isoceantransport','format','Boolean'); + WriteData(fid,prefix,'object',this,'fieldname','isstressbalance','format','Boolean'); + WriteData(fid,prefix,'object',this,'fieldname','isthermal','format','Boolean'); + WriteData(fid,prefix,'object',this,'fieldname','isgroundingline','format','Boolean'); + WriteData(fid,prefix,'object',this,'fieldname','isesa','format','Boolean'); + WriteData(fid,prefix,'object',this,'fieldname','isdamageevolution','format','Boolean'); + WriteData(fid,prefix,'object',this,'fieldname','ishydrology','format','Boolean'); + WriteData(fid,prefix,'object',this,'fieldname','ismovingfront','format','Boolean'); + WriteData(fid,prefix,'object',this,'fieldname','issampling','format','Boolean'); + WriteData(fid,prefix,'object',this,'fieldname','isslc','format','Boolean'); + WriteData(fid,prefix,'object',this,'fieldname','isoceancoupling','format','Boolean'); + WriteData(fid,prefix,'object',this,'fieldname','amr_frequency','format','Integer'); + + //process requested outputs + var outputs = this.requested_outputs; + for (var i=0;i