Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/classes/materials.m =================================================================== --- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/classes/materials.m (revision 13130) +++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/classes/materials.m (revision 13131) @@ -1,119 +0,0 @@ -%MATERIALS class definition -% -% Usage: -% materials=materials(); - -classdef materials - properties (SetAccess=public) - rho_ice = 0.; - rho_water = 0.; - rho_freshwater = 0.; - mu_water = 0.; - heatcapacity = 0.; - latentheat = 0.; - thermalconductivity = 0.; - meltingpoint = 0.; - beta = 0.; - mixed_layer_capacity = 0.; - thermal_exchange_velocity = 0.; - rheology_B = NaN; - rheology_n = NaN; - rheology_Z = NaN; - rheology_law = ''; - end - methods - function obj = materials(varargin) % {{{ - switch nargin - case 0 - obj=setdefaultparameters(obj); - otherwise - error('constructor not supported'); - end - end % }}} - function obj = setdefaultparameters(obj) % {{{ - - %ice density (kg/m^3) - obj.rho_ice=917.; - - %ocean water density (kg/m^3) - obj.rho_water=1023.; - - %fresh water density (kg/m^3) - obj.rho_freshwater=1000.; - - %water viscosity (N.s/m^2) - obj.mu_water=0.001787; - - %ice heat capacity cp (J/kg/K) - obj.heatcapacity=2093.; - - %ice latent heat of fusion L (J/kg) - obj.latentheat=3.34*10^5; - - %ice thermal conductivity (W/m/K) - obj.thermalconductivity=2.4; - - %the melting point of ice at 1 atmosphere of pressure in K - obj.meltingpoint=273.15; - - %rate of change of melting point with pressure (K/Pa) - obj.beta=9.8*10^-8; - - %mixed layer (ice-water interface) heat capacity (J/kg/K) - obj.mixed_layer_capacity=3974.; - - %thermal exchange velocity (ice-water interface) (m/s) - obj.thermal_exchange_velocity=1.00*10^-4; - - %Rheology law: what is the temperature dependence of B with T - %available: none, paterson and arrhenius - obj.rheology_law='Paterson'; - end % }}} - function md = checkconsistency(obj,md,solution,analyses) % {{{ - md = checkfield(md,'materials.rho_ice','>',0); - md = checkfield(md,'materials.rho_water','>',0); - md = checkfield(md,'materials.rho_freshwater','>',0); - md = checkfield(md,'materials.mu_water','>',0); - md = checkfield(md,'materials.rheology_B','>',0,'size',[md.mesh.numberofvertices 1]); - md = checkfield(md,'materials.rheology_n','>',0,'size',[md.mesh.numberofelements 1]); - md = checkfield(md,'materials.rheology_Z','>',0,'size',[md.mesh.numberofvertices 1]); - md = checkfield(md,'materials.rheology_law','values',{'None' 'Paterson' 'Arrhenius'}); - end % }}} - function disp(obj) % {{{ - disp(sprintf(' Materials:\n')); - - fielddisplay(obj,'rho_ice','ice density [kg/m^3]'); - fielddisplay(obj,'rho_water','ocean water density [kg/m^3]'); - fielddisplay(obj,'rho_freshwater','fresh water density [kg/m^3]'); - fielddisplay(obj,'mu_water','water viscosity [N s/m^2]'); - fielddisplay(obj,'heatcapacity','heat capacity [J/kg/K]'); - fielddisplay(obj,'thermalconductivity','ice thermal conductivity [W/m/K]'); - fielddisplay(obj,'meltingpoint','melting point of ice at 1atm in K'); - fielddisplay(obj,'latentheat','latent heat of fusion [J/m^3]'); - fielddisplay(obj,'beta','rate of change of melting point with pressure [K/Pa]'); - fielddisplay(obj,'mixed_layer_capacity','mixed layer capacity [W/kg/K]'); - fielddisplay(obj,'thermal_exchange_velocity','thermal exchange velocity [m/s]'); - fielddisplay(obj,'rheology_B','flow law parameter [Pa/s^(1/n)]'); - fielddisplay(obj,'rheology_n','Glen''s flow law exponent'); - fielddisplay(obj,'rheology_Z','rheology multiplier'); - fielddisplay(obj,'rheology_law','law for the temperature dependance of the rheology: ''None'', ''Paterson'' or ''Arrhenius'''); - end % }}} - function marshall(obj,fid) % {{{ - WriteData(fid,'object',obj,'fieldname','rho_ice','format','Double'); - WriteData(fid,'object',obj,'fieldname','rho_water','format','Double'); - WriteData(fid,'object',obj,'fieldname','rho_freshwater','format','Double'); - WriteData(fid,'object',obj,'fieldname','mu_water','format','Double'); - WriteData(fid,'object',obj,'fieldname','heatcapacity','format','Double'); - WriteData(fid,'object',obj,'fieldname','latentheat','format','Double'); - WriteData(fid,'object',obj,'fieldname','thermalconductivity','format','Double'); - WriteData(fid,'object',obj,'fieldname','meltingpoint','format','Double'); - WriteData(fid,'object',obj,'fieldname','beta','format','Double'); - WriteData(fid,'object',obj,'fieldname','mixed_layer_capacity','format','Double'); - WriteData(fid,'object',obj,'fieldname','thermal_exchange_velocity','format','Double'); - WriteData(fid,'object',obj,'fieldname','rheology_B','format','DoubleMat','mattype',1); - WriteData(fid,'object',obj,'fieldname','rheology_n','format','DoubleMat','mattype',2); - WriteData(fid,'object',obj,'fieldname','rheology_Z','format','DoubleMat','mattype',1); - WriteData(fid,'data',StringToEnum(obj.rheology_law),'enum',MaterialsRheologyLawEnum(),'format','Integer'); - end % }}} - end -end Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/classes/initialization.m =================================================================== --- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/classes/initialization.m (revision 13130) +++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/classes/initialization.m (revision 13131) @@ -81,8 +81,6 @@ WriteData(fid,'data',obj.vz,'format','DoubleMat','mattype',1,'enum',VzEnum); WriteData(fid,'data',obj.pressure,'format','DoubleMat','mattype',1,'enum',PressureEnum); WriteData(fid,'data',obj.temperature,'format','DoubleMat','mattype',1,'enum',TemperatureEnum); - WriteData(fid,'data',obj.surfacetemp,'format','DoubleMat','mattype',1,'enum',TemperatureSurfaceEnum); - WriteData(fid,'data',obj.basaltemp,'format','DoubleMat','mattype',1,'enum',TemperatureBasalEnum); WriteData(fid,'data',obj.watercolumn,'format','DoubleMat','mattype',1,'enum',WatercolumnEnum); WriteData(fid,'data',obj.waterfraction,'format','DoubleMat','mattype',1,'enum',WaterfractionEnum); end % }}} Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/classes/matice.m =================================================================== --- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/classes/matice.m (revision 0) +++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/classes/matice.m (revision 13131) @@ -0,0 +1,116 @@ +%MATICE class definition +% +% Usage: +% matice=matice(); + +classdef matice + properties (SetAccess=public) + rho_ice = 0.; + rho_water = 0.; + rho_freshwater = 0.; + mu_water = 0.; + heatcapacity = 0.; + latentheat = 0.; + thermalconductivity = 0.; + meltingpoint = 0.; + beta = 0.; + mixed_layer_capacity = 0.; + thermal_exchange_velocity = 0.; + rheology_B = NaN; + rheology_n = NaN; + rheology_law = ''; + end + methods + function obj = matice(varargin) % {{{ + switch nargin + case 0 + obj=setdefaultparameters(obj); + otherwise + error('constructor not supported'); + end + end % }}} + function obj = setdefaultparameters(obj) % {{{ + + %ice density (kg/m^3) + obj.rho_ice=917.; + + %ocean water density (kg/m^3) + obj.rho_water=1023.; + + %fresh water density (kg/m^3) + obj.rho_freshwater=1000.; + + %water viscosity (N.s/m^2) + obj.mu_water=0.001787; + + %ice heat capacity cp (J/kg/K) + obj.heatcapacity=2093.; + + %ice latent heat of fusion L (J/kg) + obj.latentheat=3.34*10^5; + + %ice thermal conductivity (W/m/K) + obj.thermalconductivity=2.4; + + %the melting point of ice at 1 atmosphere of pressure in K + obj.meltingpoint=273.15; + + %rate of change of melting point with pressure (K/Pa) + obj.beta=9.8*10^-8; + + %mixed layer (ice-water interface) heat capacity (J/kg/K) + obj.mixed_layer_capacity=3974.; + + %thermal exchange velocity (ice-water interface) (m/s) + obj.thermal_exchange_velocity=1.00*10^-4; + + %Rheology law: what is the temperature dependence of B with T + %available: none, paterson and arrhenius + obj.rheology_law='Paterson'; + end % }}} + function md = checkconsistency(obj,md,solution,analyses) % {{{ + md = checkfield(md,'materials.rho_ice','>',0); + md = checkfield(md,'materials.rho_water','>',0); + md = checkfield(md,'materials.rho_freshwater','>',0); + md = checkfield(md,'materials.mu_water','>',0); + md = checkfield(md,'materials.rheology_B','>',0,'size',[md.mesh.numberofvertices 1]); + md = checkfield(md,'materials.rheology_n','>',0,'size',[md.mesh.numberofelements 1]); + md = checkfield(md,'materials.rheology_law','values',{'None' 'Paterson' 'Arrhenius'}); + end % }}} + function disp(obj) % {{{ + disp(sprintf(' Materials:\n')); + + fielddisplay(obj,'rho_ice','ice density [kg/m^3]'); + fielddisplay(obj,'rho_water','ocean water density [kg/m^3]'); + fielddisplay(obj,'rho_freshwater','fresh water density [kg/m^3]'); + fielddisplay(obj,'mu_water','water viscosity [N s/m^2]'); + fielddisplay(obj,'heatcapacity','heat capacity [J/kg/K]'); + fielddisplay(obj,'thermalconductivity','ice thermal conductivity [W/m/K]'); + fielddisplay(obj,'meltingpoint','melting point of ice at 1atm in K'); + fielddisplay(obj,'latentheat','latent heat of fusion [J/m^3]'); + fielddisplay(obj,'beta','rate of change of melting point with pressure [K/Pa]'); + fielddisplay(obj,'mixed_layer_capacity','mixed layer capacity [W/kg/K]'); + fielddisplay(obj,'thermal_exchange_velocity','thermal exchange velocity [m/s]'); + fielddisplay(obj,'rheology_B','flow law parameter [Pa/s^(1/n)]'); + fielddisplay(obj,'rheology_n','Glen''s flow law exponent'); + fielddisplay(obj,'rheology_law','law for the temperature dependance of the rheology: ''None'', ''Paterson'' or ''Arrhenius'''); + end % }}} + function marshall(obj,fid) % {{{ + WriteData(fid,'enum',MaterialsEnum(),'data',MaticeEnum(),'format','Integer'); + WriteData(fid,'object',obj,'class','materials','fieldname','rho_ice','format','Double'); + WriteData(fid,'object',obj,'class','materials','fieldname','rho_water','format','Double'); + WriteData(fid,'object',obj,'class','materials','fieldname','rho_freshwater','format','Double'); + WriteData(fid,'object',obj,'class','materials','fieldname','mu_water','format','Double'); + WriteData(fid,'object',obj,'class','materials','fieldname','heatcapacity','format','Double'); + WriteData(fid,'object',obj,'class','materials','fieldname','latentheat','format','Double'); + WriteData(fid,'object',obj,'class','materials','fieldname','thermalconductivity','format','Double'); + WriteData(fid,'object',obj,'class','materials','fieldname','meltingpoint','format','Double'); + WriteData(fid,'object',obj,'class','materials','fieldname','beta','format','Double'); + WriteData(fid,'object',obj,'class','materials','fieldname','mixed_layer_capacity','format','Double'); + WriteData(fid,'object',obj,'class','materials','fieldname','thermal_exchange_velocity','format','Double'); + WriteData(fid,'object',obj,'class','materials','fieldname','rheology_B','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'class','materials','fieldname','rheology_n','format','DoubleMat','mattype',2); + WriteData(fid,'data',StringToEnum(obj.rheology_law),'enum',MaterialsRheologyLawEnum(),'format','Integer'); + end % }}} + end +end Index: /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/classes/model/model.m =================================================================== --- /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/classes/model/model.m (revision 13130) +++ /u/astrid-r1b/morlighe/issmuci/trunk-jpl/../trunk-jpl/src/m/classes/model/model.m (revision 13131) @@ -167,7 +167,9 @@ %materials md.materials.rheology_B=DepthAverage(md,md.materials.rheology_B); md.materials.rheology_n=project2d(md,md.materials.rheology_n,1); - md.materials.rheology_Z=DepthAverage(md,md.materials.rheology_Z); + if isa(md.materials,'matdamageice') + md.materials.rheology_Z=DepthAverage(md,md.materials.rheology_Z); + end %special for thermal modeling: md.basalforcings.melting_rate=project2d(md,md.basalforcings.melting_rate,1); @@ -716,7 +718,9 @@ %materials md.materials.rheology_B=project3d(md,'vector',md.materials.rheology_B,'type','node'); md.materials.rheology_n=project3d(md,'vector',md.materials.rheology_n,'type','element'); - md.materials.rheology_Z=project3d(md,'vector',md.materials.rheology_Z,'type','node'); + if isa(md.materials,'matdamageice') + md.materials.rheology_Z=project3d(md,'vector',md.materials.rheology_Z,'type','node'); + end %parameters md.geometry.surface=project3d(md,'vector',md.geometry.surface,'type','node'); @@ -830,10 +834,8 @@ if isfield(structmd,'penalty_offset'), md.prognostic.penalty_factor=structmd.penalty_offset; end if isfield(structmd,'B'), md.materials.rheology_B=structmd.B; end if isfield(structmd,'n'), md.materials.rheology_n=structmd.n; end - if isfield(structmd,'Z'), md.materials.rheology_Z=structmd.Z; end if isfield(structmd,'rheology_B'), md.materials.rheology_B=structmd.rheology_B; end if isfield(structmd,'rheology_n'), md.materials.rheology_n=structmd.rheology_n; end - if isfield(structmd,'rheology_Z'), md.materials.rheology_Z=structmd.rheology_Z; end if isfield(structmd,'elementoniceshelf'), md.mask.elementonfloatingice=structmd.elementoniceshelf; end if isfield(structmd,'elementonicesheet'), md.mask.elementongroundedice=structmd.elementonicesheet; end if isfield(structmd,'elementonwater'), md.mask.elementonwater=structmd.elementonwater; end @@ -1043,7 +1045,7 @@ md.rifts = rifts(); md.timestepping = timestepping(); md.groundingline = groundingline(); - md.materials = materials(); + md.materials = matice(); md.flowequation = flowequation(); md.debug = debug(); md.verbose = verbose('solution',true,'qmu',true,'control',true);