Index: /issm/trunk-jpl/src/m/classes/balancethickness.m =================================================================== --- /issm/trunk-jpl/src/m/classes/balancethickness.m (revision 22303) +++ /issm/trunk-jpl/src/m/classes/balancethickness.m (revision 22304) @@ -11,4 +11,6 @@ omega = NaN; + slopex = NaN; + slopey = NaN; end methods @@ -53,4 +55,6 @@ WriteData(fid,prefix,'object',self,'fieldname','stabilization','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','slopex','format','DoubleMat','mattype',1); + WriteData(fid,prefix,'object',self,'fieldname','slopey','format','DoubleMat','mattype',1); WriteData(fid,prefix,'object',self,'fieldname','omega','format','DoubleMat','mattype',1); end % }}} Index: /issm/trunk-jpl/src/m/classes/flowequation.m =================================================================== --- /issm/trunk-jpl/src/m/classes/flowequation.m (revision 22303) +++ /issm/trunk-jpl/src/m/classes/flowequation.m (revision 22304) @@ -99,5 +99,5 @@ md = checkfield(md,'fieldname','flowequation.isFS','numel',[1],'values',[0 1]); md = checkfield(md,'fieldname','flowequation.fe_SSA','values',{'P1','P1bubble','P1bubblecondensed','P2','P2bubble'}); - md = checkfield(md,'fieldname','flowequation.fe_HO' ,'values',{'P1','P1bubble','P1bubblecondensed','P1xP2','P2xP1','P2','P2bubble','P1xP3','P2xP4'}); + md = checkfield(md,'fieldname','flowequation.fe_HO' ,'values',{'P1','P1bubble','P1bubblecondensed','P1xP2','P2xP1','P2','P2bubble','P1xP3','P1xP4','P2xP4'}); md = checkfield(md,'fieldname','flowequation.fe_FS' ,'values',{'P1P1','P1P1GLS','MINIcondensed','MINI','TaylorHood','LATaylorHood','XTaylorHood','OneLayerP4z','CrouzeixRaviart','LACrouzeixRaviart'}); md = checkfield(md,'fieldname','flowequation.augmented_lagrangian_r','numel',[1],'>=',0.); Index: /issm/trunk-jpl/src/m/classes/materials.m =================================================================== --- /issm/trunk-jpl/src/m/classes/materials.m (revision 22303) +++ /issm/trunk-jpl/src/m/classes/materials.m (revision 22304) @@ -16,118 +16,16 @@ self.nature=varargin; end - - %check this is acceptable: + + %start filling in the dynamic fields: for i=1:length(self.nature), - if ~(strcmpi(self.nature{i},'litho') | strcmpi(self.nature{i},'ice')), - error('materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'')'); + list = listoffields(self,self.nature{i}); + for j=1:size(list,1), + %Add property to materials + self.addprop(list{j,1}); + %Set default value according to table + self.(list{j,1}) = list{j,2}; end end - - %start filling in the dynamic fields: - for i=1:length(self.nature), - nat=self.nature{i}; - switch nat - case 'ice' - self.addprop('rho_ice'); - self.addprop('rho_water'); - self.addprop('rho_freshwater'); - self.addprop('mu_water'); - self.addprop('heatcapacity'); - self.addprop('latentheat'); - self.addprop('thermalconductivity'); - self.addprop('temperateiceconductivity'); - self.addprop('meltingpoint'); - self.addprop('beta'); - self.addprop('mixed_layer_capacity'); - self.addprop('thermal_exchange_velocity'); - self.addprop('rheology_B'); - self.addprop('rheology_n'); - self.addprop('rheology_law'); - case 'litho' - self.addprop('numlayers'); - self.addprop('radius'); - self.addprop('viscosity'); - self.addprop('lame_lambda'); - self.addprop('lame_mu'); - self.addprop('burgers_viscosity'); - self.addprop('burgers_mu'); - self.addprop('isburgers'); - self.addprop('density'); - self.addprop('issolid'); - otherwise - error('materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'')'); - end - end - %set default parameters: - self.setdefaultparameters(); - end % }}} - function self = setdefaultparameters(self) % {{{ - - for i=1:length(self.nature), - nat=self.nature{i}; - switch nat - case 'ice' - %ice density (kg/m^3) - self.rho_ice=917.; - - %ocean water density (kg/m^3) - self.rho_water=1023.; - - %fresh water density (kg/m^3) - self.rho_freshwater=1000.; - - %water viscosity (N.s/m^2) - self.mu_water=0.001787; - - %ice heat capacity cp (J/kg/K) - self.heatcapacity=2093.; - - %ice latent heat of fusion L (J/kg) - self.latentheat=3.34*10^5; - - %ice thermal conductivity (W/m/K) - self.thermalconductivity=2.4; - - %wet ice thermal conductivity (W/m/K) - self.temperateiceconductivity=.24; - - %the melting point of ice at 1 atmosphere of pressure in K - self.meltingpoint=273.15; - - %rate of change of melting point with pressure (K/Pa) - self.beta=9.8*10^-8; - - %mixed layer (ice-water interface) heat capacity (J/kg/K) - self.mixed_layer_capacity=3974.; - - %thermal exchange velocity (ice-water interface) (m/s) - self.thermal_exchange_velocity=1.00*10^-4; - - %Rheology law: what is the temperature dependence of B with T - %available: none, paterson and arrhenius - self.rheology_law='Paterson'; - - case 'litho' - %we default to a configuration that enables running GIA solutions using giacaron and/or giaivins. - self.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) - self.radius=[1e3;6278*1e3;6378*1e3]; - - self.viscosity=[1e21;1e40]; %mantle and lithosphere viscosity (respectively) [Pa.s] - self.lame_mu=[1.45*1e11;6.7*10^10]; % (Pa) %lithosphere and mantle shear modulus (respectively) [Pa] - self.lame_lambda=self.lame_mu; % (Pa) %mantle and lithosphere lamba parameter (respectively) [Pa] - self.burgers_viscosity=[NaN;NaN]; - self.burgers_mu=[NaN;NaN]; - self.isburgers=[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. - - otherwise - error('materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'')'); - end - end end % }}} function disp(self) % {{{ @@ -135,37 +33,8 @@ for i=1:length(self.nature), - nat=self.nature{i}; - switch nat - case 'ice' - disp(sprintf(' \nIce:')); - fielddisplay(self,'rho_ice','ice density [kg/m^3]'); - fielddisplay(self,'rho_water','ocean water density [kg/m^3]'); - fielddisplay(self,'rho_freshwater','fresh water density [kg/m^3]'); - fielddisplay(self,'mu_water','water viscosity [N s/m^2]'); - fielddisplay(self,'heatcapacity','heat capacity [J/kg/K]'); - fielddisplay(self,'thermalconductivity',['ice thermal conductivity [W/m/K]']); - fielddisplay(self,'temperateiceconductivity','temperate ice thermal conductivity [W/m/K]'); - fielddisplay(self,'meltingpoint','melting point of ice at 1atm in K'); - fielddisplay(self,'latentheat','latent heat of fusion [J/kg]'); - fielddisplay(self,'beta','rate of change of melting point with pressure [K/Pa]'); - fielddisplay(self,'mixed_layer_capacity','mixed layer capacity [W/kg/K]'); - fielddisplay(self,'thermal_exchange_velocity','thermal exchange velocity [m/s]'); - fielddisplay(self,'rheology_B','flow law parameter [Pa/s^(1/n)]'); - fielddisplay(self,'rheology_n','Glen''s flow law exponent'); - fielddisplay(self,'rheology_law',['law for the temperature dependance of the rheology: ''None'', ''BuddJacka'', Cuffey'', ''CuffeyTemperate'', ''Paterson'', ''Arrhenius'' or ''LliboutryDuval''']); - case 'litho' - disp(sprintf(' \nLitho:')); - fielddisplay(self,'numlayers','number of layers (default 2)'); - fielddisplay(self,'radius','array describing the radius for each interface (numlayers+1) [m]'); - fielddisplay(self,'viscosity','array describing each layer''s viscosity (numlayers) [Pa.s]'); - fielddisplay(self,'lame_lambda','array describing the lame lambda parameter (numlayers) [Pa]'); - fielddisplay(self,'lame_mu','array describing the shear modulus for each layers (numlayers) [Pa]'); - 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,'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)'); - otherwise - error('materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'')'); + list = listoffields(self,self.nature{i}); + disp(sprintf('\n %s:',upper(self.nature{i}))); + for j=1:size(list,1), + fielddisplay(self,list{j,1},list{j,3}); end end @@ -173,43 +42,34 @@ function md = checkconsistency(self,md,solution,analyses) % {{{ - for i=1:length(self.nature), - nat=self.nature{i}; - switch nat - case 'ice' - md = checkfield(md,'fieldname','materials.rho_ice','>',0); - md = checkfield(md,'fieldname','materials.rho_water','>',0); - md = checkfield(md,'fieldname','materials.rho_freshwater','>',0); - md = checkfield(md,'fieldname','materials.mu_water','>',0); - md = checkfield(md,'fieldname','materials.rheology_B','>',0,'timeseries',1,'NaN',1,'Inf',1); - md = checkfield(md,'fieldname','materials.rheology_n','>',0,'size',[md.mesh.numberofelements 1]); - md = checkfield(md,'fieldname','materials.rheology_law','values',{'None' 'BuddJacka' 'Cuffey' 'CuffeyTemperate' 'Paterson' 'Arrhenius' 'LliboutryDuval'}); - case 'litho' - if ~ismember('LoveAnalysis',analyses), return; end - md = checkfield(md,'fieldname','materials.numlayers','NaN',1,'Inf',1,'>',0,'numel',1); - md = checkfield(md,'fieldname','materials.radius','NaN',1,'Inf',1,'size',[md.materials.numlayers+1 1],'>',0); - md = checkfield(md,'fieldname','materials.lame_mu','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>=',0); - md = checkfield(md,'fieldname','materials.lame_lambda','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>=',0); - md = checkfield(md,'fieldname','materials.issolid','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>=',0,'<',2); - 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.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); + if isprop(self,'rho_ice'), md = checkfield(md,'fieldname','materials.rho_ice','>',0); end + if isprop(self,'rho_water'), md = checkfield(md,'fieldname','materials.rho_water','>',0); end + if isprop(self,'rho_freshwater'), md = checkfield(md,'fieldname','materials.rho_freshwater','>',0); end + if isprop(self,'mu_water'), md = checkfield(md,'fieldname','materials.mu_water','>',0); end + if isprop(self,'rheology_B'), md = checkfield(md,'fieldname','materials.rheology_B','>',0,'timeseries',1,'NaN',1,'Inf',1); end + if isprop(self,'rheology_n'), md = checkfield(md,'fieldname','materials.rheology_n','>',0,'size',[md.mesh.numberofelements 1]); end + if isprop(self,'rheology_law'), md = checkfield(md,'fieldname','materials.rheology_law','values',{'None' 'BuddJacka' 'Cuffey' 'CuffeyTemperate' 'Paterson' 'Arrhenius' 'LliboutryDuval'}); end + if ~ismember('LoveAnalysis',analyses), return; end + if isprop(self,'numlayers'), md = checkfield(md,'fieldname','materials.numlayers','NaN',1,'Inf',1,'>',0,'numel',1); end + if isprop(self,'radius'), md = checkfield(md,'fieldname','materials.radius','NaN',1,'Inf',1,'size',[md.materials.numlayers+1 1],'>',0); end + if isprop(self,'lame_mu'), md = checkfield(md,'fieldname','materials.lame_mu','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>=',0); end + if isprop(self,'lame_lambda'), md = checkfield(md,'fieldname','materials.lame_lambda','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>=',0); end + if isprop(self,'issolid'), md = checkfield(md,'fieldname','materials.issolid','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>=',0,'<',2); end + if isprop(self,'density'), md = checkfield(md,'fieldname','materials.density','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>',0); end + if isprop(self,'viscosity'), md = checkfield(md,'fieldname','materials.viscosity','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>=',0); end + if isprop(self,'isburgers'), md = checkfield(md,'fieldname','materials.isburgers','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>=',0,'<=',1); end + if isprop(self,'burgers_viscosity'), md = checkfield(md,'fieldname','materials.burgers_viscosity','Inf',1,'size',[md.materials.numlayers 1],'>=',0); end + if isprop(self,'burgers_mu'), md = checkfield(md,'fieldname','materials.burgers_mu','Inf',1,'size',[md.materials.numlayers 1],'>=',0); end - 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'); - end - end - if md.materials.issolid(1)==0 | md.materials.lame_mu(1)==0 - error('First layer must be solid (issolid(1) > 0 AND lame_mu(1) > 0). Add a weak inner core if necessary.'); - end - ind=find(md.materials.issolid==0); - 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.']) - end - otherwise - error('materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'')'); + 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'); end + end + if md.materials.issolid(1)==0 | md.materials.lame_mu(1)==0 + error('First layer must be solid (issolid(1) > 0 AND lame_mu(1) > 0). Add a weak inner core if necessary.'); + end + ind=find(md.materials.issolid==0); + 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.']) end @@ -221,109 +81,118 @@ WriteData(fid,prefix,'name','md.materials.nature','data',naturetointeger(self.nature),'format','IntMat','mattype',3); - for i=1:length(self.nature), - nat=self.nature{i}; - switch nat - case 'ice' - WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_ice','format','Double'); - WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_water','format','Double'); - WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_freshwater','format','Double'); - WriteData(fid,prefix,'object',self,'class','materials','fieldname','mu_water','format','Double'); - WriteData(fid,prefix,'object',self,'class','materials','fieldname','heatcapacity','format','Double'); - WriteData(fid,prefix,'object',self,'class','materials','fieldname','latentheat','format','Double'); - WriteData(fid,prefix,'object',self,'class','materials','fieldname','thermalconductivity','format','Double'); - WriteData(fid,prefix,'object',self,'class','materials','fieldname','temperateiceconductivity','format','Double'); - WriteData(fid,prefix,'object',self,'class','materials','fieldname','meltingpoint','format','Double'); - WriteData(fid,prefix,'object',self,'class','materials','fieldname','beta','format','Double'); - WriteData(fid,prefix,'object',self,'class','materials','fieldname','mixed_layer_capacity','format','Double'); - WriteData(fid,prefix,'object',self,'class','materials','fieldname','thermal_exchange_velocity','format','Double'); - WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_B','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); - WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_n','format','DoubleMat','mattype',2); - WriteData(fid,prefix,'data',self.rheology_law,'name','md.materials.rheology_law','format','String'); - case 'litho' - WriteData(fid,prefix,'object',self,'class','materials','fieldname','numlayers','format','Integer'); - WriteData(fid,prefix,'object',self,'class','materials','fieldname','radius','format','DoubleMat','mattype',3); - WriteData(fid,prefix,'object',self,'class','materials','fieldname','lame_mu','format','DoubleMat','mattype',3); - WriteData(fid,prefix,'object',self,'class','materials','fieldname','lame_lambda','format','DoubleMat','mattype',3); - WriteData(fid,prefix,'object',self,'class','materials','fieldname','issolid','format','DoubleMat','mattype',3); - 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','burgers_viscosity','format','DoubleMat','mattype',3); - WriteData(fid,prefix,'object',self,'class','materials','fieldname','burgers_mu','format','DoubleMat','mattype',3); - otherwise - error('materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'')'); + if isprop(self,'rho_ice'), WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_ice','format','Double'); end + if isprop(self,'rho_water'), WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_water','format','Double'); end + if isprop(self,'rho_freshwater'), WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_freshwater','format','Double'); end + if isprop(self,'mu_water'), WriteData(fid,prefix,'object',self,'class','materials','fieldname','mu_water','format','Double'); end + if isprop(self,'heatcapacity'), WriteData(fid,prefix,'object',self,'class','materials','fieldname','heatcapacity','format','Double'); end + if isprop(self,'latentheat'), WriteData(fid,prefix,'object',self,'class','materials','fieldname','latentheat','format','Double'); end + if isprop(self,'thermalconductivity'), WriteData(fid,prefix,'object',self,'class','materials','fieldname','thermalconductivity','format','Double'); end + if isprop(self,'temperateiceconductivity'), WriteData(fid,prefix,'object',self,'class','materials','fieldname','temperateiceconductivity','format','Double'); end + if isprop(self,'meltingpoint'), WriteData(fid,prefix,'object',self,'class','materials','fieldname','meltingpoint','format','Double'); end + if isprop(self,'beta'), WriteData(fid,prefix,'object',self,'class','materials','fieldname','beta','format','Double'); end + if isprop(self,'mixed_layer_capacity'), WriteData(fid,prefix,'object',self,'class','materials','fieldname','mixed_layer_capacity','format','Double'); end + if isprop(self,'thermal_exchange_velocity'), WriteData(fid,prefix,'object',self,'class','materials','fieldname','thermal_exchange_velocity','format','Double'); end + if isprop(self,'rheology_B'), WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_B','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); end + if isprop(self,'rheology_n'), WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_n','format','DoubleMat','mattype',2); end + if isprop(self,'rheology_law'), Write Data(fid,prefix,'data',self.rheology_law,'nae','md.materials.rheology_law','format','String'); end + if isprop(self,'numlayers'), WriteData(fid,prefix,'object',self,'class','materials','fieldname','numlayers','format','Integer'); end + if isprop(self,'radius'), WriteData(fid,prefix,'object',self,'class','materials','fieldname','radius','format','DoubleMat','mattype',3); end + if isprop(self,'lame_mu'), WriteData(fid,prefix,'object',self,'class','materials','fieldname','lame_mu','format','DoubleMat','mattype',3); end + if isprop(self,'lame_lambda'), WriteData(fid,prefix,'object',self,'class','materials','fieldname','lame_lambda','format','DoubleMat','mattype',3); end + if isprop(self,'issolid'), WriteData(fid,prefix,'object',self,'class','materials','fieldname','issolid','format','DoubleMat','mattype',3); end + if isprop(self,'density'), WriteData(fid,prefix,'object',self,'class','materials','fieldname','density','format','DoubleMat','mattype',3); end + if isprop(self,'viscosity'), WriteData(fid,prefix,'object',self,'class','materials','fieldname','viscosity','format','DoubleMat','mattype',3); end + if isprop(self,'isburgers'), WriteData(fid,prefix,'object',self,'class','materials','fieldname','isburgers','format','DoubleMat','mattype',3); end + if isprop(self,'burgers_viscosity'), WriteData(fid,prefix,'object',self,'class','materials','fieldname','burgers_viscosity','format','DoubleMat','mattype',3); end + if isprop(self,'burgers_mu'), WriteData(fid,prefix,'object',self,'class','materials','fieldname','burgers_mu','format','DoubleMat','mattype',3); end + end % }}} + function self = extrude(self,md) % {{{ + if isprop(self,'rheology_B'), self.rheology_B=project3d(md,'vector',self.rheology_B,'type','node'); end + if isprop(self,'rheology_n'), self.rheology_n=project3d(md,'vector',self.rheology_n,'type','element'); end + end % }}} + function savemodeljs(self,fid,modelname) % {{{ + + proplist = properties(self); + + for i=1:numel(proplist), + if numel(self.(proplist{i}))==1, + writejsdouble(fid,[modelname '.materials.' proplist{i}],self.(proplist{i})); + elseif size(self.(proplist{i}),2)==1, + writejs1Darray(fid,[modelname '.materials.' proplist{i}],self.(proplist{i})); + else + error('not supported yet'); end end end % }}} - function self = extrude(self,md) % {{{ - for i=1:length(self.nature), - nat=self.nature{i}; - switch nat + function list = listoffields(self,nat) % {{{ + + switch nat case 'ice' - self.rheology_B=project3d(md,'vector',self.rheology_B,'type','node'); - self.rheology_n=project3d(md,'vector',self.rheology_n,'type','element'); - end + list={... + 'rho_ice' , 917., 'ice density [kg/m^3]'; ... + 'rho_water' ,1023., 'ocean water density [kg/m^3]'; ... + 'rho_freshwater' ,1000., 'fresh water density [kg/m^3]';... + 'mu_water' ,0.001787, 'water viscosity [N s/m^2]';... + 'heatcapacity' ,2093., 'heat capacity [J/kg/K]';... + 'latentheat' ,3.34e5, 'latent heat of fusion [J/kg]';... + 'thermalconductivity' ,2.4, 'ice thermal conductivity [W/m/K]';... + 'temperateiceconductivity' ,.24, 'temperate ice thermal conductivity [W/m/K]';... + 'meltingpoint' ,273.15, 'melting point of ice at 1 atm in K';... + 'beta' ,9.8e-8 'rate of change of melting point with pressure [K/Pa]';... + 'mixed_layer_capacity' ,3974., 'mixed layer capacity [W/kg/K]'; ... + 'thermal_exchange_velocity' ,1e-4 'thermal exchange velocity [m/s]'; ... + 'rheology_B' ,NaN, 'flow law parameter [Pa/s^(1/n)]'; ... + 'rheology_n' ,NaN, 'Glen''s flow law exponent'; ... + 'rheology_law' ,'Paterson','law for the temperature dependance of the rheology: ''None'' , ''BuddJacka'' , Cuffey'' , ''CuffeyTemperate'' , ''Paterson'' , ''Arrhenius'' or ''LliboutryDuval'''}; + case 'damageice' + list = listoffields(self,'ice'), + case 'enhancedice' + list = listoffields(self,'ice'); + list(end+1,:) = {'rheology_E',NaN,'enhancement factor'}; + case 'estarice' + list = listoffields(self,'ice'); + list(end+1,:) = {'rheology_Ec',NaN,'compressive enhancement factor'}; + list(end+1,:) = {'rheology_Es',NaN,'shear enhancement factor'}; + list(find(strcmp(list(:,1),'rheology_n')),:) = []; + case 'litho' + %Note for the radius: + % center of the earth (approximation, must not be 0), then the lab (lithosphere/asthenosphere boundary) then the surface + % (with 1d3 to avoid numerical singularities) + %Note: here we consider by default two layers: mantle and lithosphere + list={... + 'numlayers' ,2, 'number of layers (default 2)'; ... + 'radius' ,[1e3;6278*1e3;6378*1e3] 'array describing the radius for each interface (numlayers+1) [m]'; ... + 'viscosity' ,[1e21;1e40], 'array describing each layer''s viscosity (numlayers) [Pa.s]'; ... + 'lame_lambda' ,[1.45e11;6.7e10], 'array describing the lame lambda parameter (numlayers) [Pa]'; ... + 'lame_mu' ,[1.45e11;6.7e10], 'array describing the shear modulus for each layers (numlayers) [Pa]'; ... + 'burgers_viscosity' ,[NaN;NaN], 'array describing each layer''s transient viscosity, only for Burgers rheologies (numlayers) [Pa.s]',... + 'burgers_mu' ,[0;0], 'array describing each layer''s transient shear modulus, only for Burgers rheologies (numlayers) [Pa]'; ... + 'isburgers' ,0, 'array describing whether we adopt a MaxWell (0) or Burgers (1) rheology (default 0)'; ... + 'density' ,[5.51*1e3;5.50*1e3], 'array describing each layer''s density (numlayers) [kg/m^3]'; ... + 'issolid' ,[1;1], 'array describing whether the layer is solid or liquid (default 1) (numlayers)' + } + otherwise + error('nature of the material not supported yet!'); end end % }}} - function savemodeljs(self,fid,modelname) % {{{ - - for i=1:length(self.nature), - nat=self.nature{i}; - switch nat - case 'ice' - writejsdouble(fid,[modelname '.materials.rho_ice'],self.rho_ice); - writejsdouble(fid,[modelname '.materials.rho_water'],self.rho_water); - writejsdouble(fid,[modelname '.materials.rho_freshwater'],self.rho_freshwater); - writejsdouble(fid,[modelname '.materials.mu_water'],self.mu_water); - writejsdouble(fid,[modelname '.materials.heatcapacity'],self.heatcapacity); - writejsdouble(fid,[modelname '.materials.latentheat'],self.latentheat); - writejsdouble(fid,[modelname '.materials.thermalconductivity'],self.thermalconductivity); - writejsdouble(fid,[modelname '.materials.temperateiceconductivity'],self.temperateiceconductivity); - writejsdouble(fid,[modelname '.materials.meltingpoint'],self.meltingpoint); - writejsdouble(fid,[modelname '.materials.beta'],self.beta); - writejsdouble(fid,[modelname '.materials.mixed_layer_capacity'],self.mixed_layer_capacity); - writejsdouble(fid,[modelname '.materials.thermal_exchange_velocity'],self.thermal_exchange_velocity); - writejsdouble(fid,[modelname '.materials.mixed_layer_capacity'],self.mixed_layer_capacity); - writejs1Darray(fid,[modelname '.materials.rheology_B'],self.rheology_B); - writejs1Darray(fid,[modelname '.materials.rheology_n'],self.rheology_n); - writejsstring(fid,[modelname '.materials.rheology_law'],self.rheology_law); - case 'litho' - writejsdouble(fid,[modelname '.materials.numlayers'],self.numlayers); - writejsdouble(fid,[modelname '.materials.radius'],self.radius); - writejsdouble(fid,[modelname '.materials.lame_mu'],self.lame_mu); - writejsdouble(fid,[modelname '.materials.lame_lambda'],self.lame_lambda); - writejsdouble(fid,[modelname '.materials.issolid'],self.issolid); - 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.burgers_viscosity'],self.burgers_viscosity); - writejsdouble(fid,[modelname '.materials.burgers_mu'],self.burgers_mu); - otherwise - error('materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'')'); - end - end - - - - end % }}} - end -end function intnat = naturetointeger(strnat) % {{{ intnat=zeros(length(strnat),1); for i=1:length(strnat), switch strnat{i}, - case 'damageice' - intnat(i)=1; - case 'estar' - intnat(i)=2; - case 'ice' - intnat(i)=3; - case 'enhancedice' - intnat(i)=4; - case 'litho' - intnat(i)=5; - otherwise - error('materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'')'); + case 'damageice' + intnat(i)=1; + case 'estar' + intnat(i)=2; + case 'ice' + intnat(i)=3; + case 'enhancedice' + intnat(i)=4; + case 'litho' + intnat(i)=5; + otherwise + error('materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'')'); end end end % }}} + end +end