Index: /issm/trunk-jpl/src/m/classes/materials.m =================================================================== --- /issm/trunk-jpl/src/m/classes/materials.m (revision 22304) +++ /issm/trunk-jpl/src/m/classes/materials.m (revision 22305) @@ -16,16 +16,118 @@ self.nature=varargin; end - + + %check this is acceptable: + 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'')'); + end + end + %start filling in the dynamic fields: for i=1:length(self.nature), - 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 - + 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) % {{{ @@ -33,8 +135,37 @@ for i=1:length(self.nature), - 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}); + 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'')'); end end @@ -42,34 +173,43 @@ function md = checkconsistency(self,md,solution,analyses) % {{{ - 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.']) + 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); + + 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'')'); + end end @@ -81,118 +221,109 @@ WriteData(fid,prefix,'name','md.materials.nature','data',naturetointeger(self.nature),'format','IntMat','mattype',3); - 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 % }}} + 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'')'); + end + 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 + for i=1:length(self.nature), + nat=self.nature{i}; + 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 + 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 list = listoffields(self,nat) % {{{ - - switch nat - case 'ice' - 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 % }}} + + 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'')'); - end - end - end % }}} - end -end + 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 % }}}