Changeset 23697

Timestamp:

02/06/19 03:09:00 (6 years ago)

Author:

rueckamp

Message:

CHG: added different computations for the Effective Conductivity.

Location:

issm/trunk-jpl/src

Files:

• c/analyses/EnthalpyAnalysis.cpp (modified) (4 diffs)
• c/modules/ModelProcessorx/CreateParameters.cpp (modified) (2 diffs)
• c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• c/shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• c/shared/Enum/StringToEnumx.cpp (modified) (9 diffs)
• m/classes/matdamageice.m (modified) (5 diffs)
• m/classes/matdamageice.py (modified) (5 diffs)
• m/classes/matenhancedice.m (modified) (6 diffs)
• m/classes/matenhancedice.py (modified) (5 diffs)
• m/classes/matestar.m (modified) (6 diffs)
• m/classes/matestar.py (modified) (5 diffs)
• m/classes/matice.m (modified) (6 diffs)
• m/classes/matice.py (modified) (5 diffs)

issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp

@@ -1172 +1172 @@
         /*Get pressures and enthalpies on vertices*/
         int         numvertices = element->GetNumberOfVertices();
+        int         effectiveconductivity_averaging;
         IssmDouble* pressures   = xNew<IssmDouble>(numvertices);
         IssmDouble* enthalpies  = xNew<IssmDouble>(numvertices);
@@ -1178 +1179 @@
         element->GetInputListOnVertices(pressures,PressureEnum);
         element->GetInputListOnVertices(enthalpies,enthalpy_enum);
+        element->FindParam(&effectiveconductivity_averaging,MaterialsEffectiveconductivityAveragingEnum);
+
         for(iv=0;iv<numvertices;iv++){
                 PIE[iv]   = PureIceEnthalpy(element,pressures[iv]);
@@ -1195 +1198 @@
         }
         else{
-                /* return harmonic mean of thermal conductivities, weighted by fraction of cold/temperate ice,
-                        cf Patankar 1980, pp44 */
                 kappa_c = EnthalpyDiffusionParameter(element,PureIceEnthalpy(element,0.)-1.,0.);
                 kappa_t = EnthalpyDiffusionParameter(element,PureIceEnthalpy(element,0.)+1.,0.);
+
                 Hc=0.; Ht=0.;
                 for(iv=0; iv<numvertices;iv++){
@@ -1210 +1212 @@
                 _assert_(lambda>=0.);
                 _assert_(lambda<=1.);
-                kappa  = kappa_c*kappa_t/(lambda*kappa_t+(1.-lambda)*kappa_c); // ==(lambda/kappa_c + (1.-lambda)/kappa_t)^-1
+
+                if(effectiveconductivity_averaging==0){
+                        /* return arithmetic mean (volume average) of thermal conductivities, weighted by fraction of cold/temperate ice */
+                        kappa=kappa_c*lambda+(1.-lambda)*kappa_t;
+                }
+                else if(effectiveconductivity_averaging==1){
+                        /* return harmonic mean (reciprocal avarage) of thermal conductivities, weighted by fraction of cold/temperate ice, cf Patankar 1980, pp44 */
+                        kappa=kappa_c*kappa_t/(lambda*kappa_t+(1.-lambda)*kappa_c); 
+                }
+                else if(effectiveconductivity_averaging==2){
+                        /* return geometric mean (power law) of thermal conductivities, weighted by fraction of cold/temperate ice */
+                        kappa=pow(kappa_c,lambda)*pow(kappa_t,1.-lambda); 
+                }
+                else{
+                        _error_("effectiveconductivity_averaging not supported yet");
+                }
         }       
 

issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp

@@ -288 +288 @@
                         parameters->AddObject(iomodel->CopyConstantObject("md.materials.thermalconductivity",MaterialsThermalconductivityEnum));
                         parameters->AddObject(iomodel->CopyConstantObject("md.materials.temperateiceconductivity",MaterialsTemperateiceconductivityEnum));
+                        parameters->AddObject(iomodel->CopyConstantObject("md.materials.effectiveconductivity_averaging",MaterialsEffectiveconductivityAveragingEnum));
                         parameters->AddObject(iomodel->CopyConstantObject("md.materials.latentheat",MaterialsLatentheatEnum));
                         parameters->AddObject(iomodel->CopyConstantObject("md.materials.beta",MaterialsBetaEnum));
@@ -325 +326 @@
                                                 parameters->AddObject(iomodel->CopyConstantObject("md.materials.thermalconductivity",MaterialsThermalconductivityEnum));
                                                 parameters->AddObject(iomodel->CopyConstantObject("md.materials.temperateiceconductivity",MaterialsTemperateiceconductivityEnum));
+                                                parameters->AddObject(iomodel->CopyConstantObject("md.materials.effectiveconductivity_averaging",MaterialsEffectiveconductivityAveragingEnum));
                                                 parameters->AddObject(iomodel->CopyConstantObject("md.materials.latentheat",MaterialsLatentheatEnum));
                                                 parameters->AddObject(iomodel->CopyConstantObject("md.materials.beta",MaterialsBetaEnum));

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

@@ -240 +240 @@
         MaterialsRhoSeawaterEnum,
         MaterialsTemperateiceconductivityEnum,
+        MaterialsEffectiveconductivityAveragingEnum,
         MaterialsThermalconductivityEnum,
         MaterialsThermalExchangeVelocityEnum,

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

@@ -248 +248 @@
                 case MaterialsRhoSeawaterEnum : return "MaterialsRhoSeawater";
                 case MaterialsTemperateiceconductivityEnum : return "MaterialsTemperateiceconductivity";
+                case MaterialsEffectiveconductivityAveragingEnum : return "MaterialsEffectiveconductivityAveraging";
                 case MaterialsThermalconductivityEnum : return "MaterialsThermalconductivity";
                 case MaterialsThermalExchangeVelocityEnum : return "MaterialsThermalExchangeVelocity";

issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

@@ -251 +251 @@
               else if (strcmp(name,"MaterialsRhoSeawater")==0) return MaterialsRhoSeawaterEnum;
               else if (strcmp(name,"MaterialsTemperateiceconductivity")==0) return MaterialsTemperateiceconductivityEnum;
+              else if (strcmp(name,"MaterialsEffectiveconductivityAveraging")==0) return MaterialsEffectiveconductivityAveragingEnum;
               else if (strcmp(name,"MaterialsThermalconductivity")==0) return MaterialsThermalconductivityEnum;
               else if (strcmp(name,"MaterialsThermalExchangeVelocity")==0) return MaterialsThermalExchangeVelocityEnum;
@@ -259 +260 @@
               else if (strcmp(name,"MeshNumberofvertices")==0) return MeshNumberofverticesEnum;
               else if (strcmp(name,"ModelId")==0) return ModelIdEnum;
-              else if (strcmp(name,"Nodes")==0) return NodesEnum;
          else stage=3;
    }
    if(stage==3){
-              if (strcmp(name,"NumModels")==0) return NumModelsEnum;
+              if (strcmp(name,"Nodes")==0) return NodesEnum;
+              else if (strcmp(name,"NumModels")==0) return NumModelsEnum;
               else if (strcmp(name,"OceanGridNx")==0) return OceanGridNxEnum;
               else if (strcmp(name,"OceanGridNy")==0) return OceanGridNyEnum;
@@ -382 +383 @@
               else if (strcmp(name,"StressbalanceMaxiter")==0) return StressbalanceMaxiterEnum;
               else if (strcmp(name,"StressbalanceNumRequestedOutputs")==0) return StressbalanceNumRequestedOutputsEnum;
-              else if (strcmp(name,"StressbalancePenaltyFactor")==0) return StressbalancePenaltyFactorEnum;
          else stage=4;
    }
    if(stage==4){
-              if (strcmp(name,"StressbalanceReltol")==0) return StressbalanceReltolEnum;
+              if (strcmp(name,"StressbalancePenaltyFactor")==0) return StressbalancePenaltyFactorEnum;
+              else if (strcmp(name,"StressbalanceReltol")==0) return StressbalanceReltolEnum;
               else if (strcmp(name,"StressbalanceRequestedOutputs")==0) return StressbalanceRequestedOutputsEnum;
               else if (strcmp(name,"StressbalanceRestol")==0) return StressbalanceRestolEnum;
@@ -505 +506 @@
               else if (strcmp(name,"FrictionAs")==0) return FrictionAsEnum;
               else if (strcmp(name,"FrictionC")==0) return FrictionCEnum;
-              else if (strcmp(name,"FrictionCoefficientcoulomb")==0) return FrictionCoefficientcoulombEnum;
          else stage=5;
    }
    if(stage==5){
-              if (strcmp(name,"FrictionCoefficient")==0) return FrictionCoefficientEnum;
+              if (strcmp(name,"FrictionCoefficientcoulomb")==0) return FrictionCoefficientcoulombEnum;
+              else if (strcmp(name,"FrictionCoefficient")==0) return FrictionCoefficientEnum;
               else if (strcmp(name,"FrictionEffectivePressure")==0) return FrictionEffectivePressureEnum;
               else if (strcmp(name,"FrictionM")==0) return FrictionMEnum;
@@ -628 +629 @@
               else if (strcmp(name,"SmbD")==0) return SmbDEnum;
               else if (strcmp(name,"SmbDini")==0) return SmbDiniEnum;
-              else if (strcmp(name,"SmbDlwrf")==0) return SmbDlwrfEnum;
          else stage=6;
    }
    if(stage==6){
-              if (strcmp(name,"SmbDswrf")==0) return SmbDswrfEnum;
+              if (strcmp(name,"SmbDlwrf")==0) return SmbDlwrfEnum;
+              else if (strcmp(name,"SmbDswrf")==0) return SmbDswrfEnum;
               else if (strcmp(name,"SmbDz")==0) return SmbDzEnum;
               else if (strcmp(name,"SmbDzini")==0) return SmbDziniEnum;
@@ -751 +752 @@
               else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
               else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum;
-              else if (strcmp(name,"Adjointp")==0) return AdjointpEnum;
          else stage=7;
    }
    if(stage==7){
-              if (strcmp(name,"AggressiveMigration")==0) return AggressiveMigrationEnum;
+              if (strcmp(name,"Adjointp")==0) return AdjointpEnum;
+              else if (strcmp(name,"AggressiveMigration")==0) return AggressiveMigrationEnum;
               else if (strcmp(name,"AmrBamg")==0) return AmrBamgEnum;
               else if (strcmp(name,"AmrNeopz")==0) return AmrNeopzEnum;
@@ -874 +875 @@
               else if (strcmp(name,"GroundedArea")==0) return GroundedAreaEnum;
               else if (strcmp(name,"GroundedAreaScaled")==0) return GroundedAreaScaledEnum;
-              else if (strcmp(name,"GroundingOnly")==0) return GroundingOnlyEnum;
          else stage=8;
    }
    if(stage==8){
-              if (strcmp(name,"Gset")==0) return GsetEnum;
+              if (strcmp(name,"GroundingOnly")==0) return GroundingOnlyEnum;
+              else if (strcmp(name,"Gset")==0) return GsetEnum;
               else if (strcmp(name,"Gsl")==0) return GslEnum;
               else if (strcmp(name,"HOApproximation")==0) return HOApproximationEnum;
@@ -997 +998 @@
               else if (strcmp(name,"Outputdefinition13")==0) return Outputdefinition13Enum;
               else if (strcmp(name,"Outputdefinition14")==0) return Outputdefinition14Enum;
-              else if (strcmp(name,"Outputdefinition15")==0) return Outputdefinition15Enum;
          else stage=9;
    }
    if(stage==9){
-              if (strcmp(name,"Outputdefinition16")==0) return Outputdefinition16Enum;
+              if (strcmp(name,"Outputdefinition15")==0) return Outputdefinition15Enum;
+              else if (strcmp(name,"Outputdefinition16")==0) return Outputdefinition16Enum;
               else if (strcmp(name,"Outputdefinition17")==0) return Outputdefinition17Enum;
               else if (strcmp(name,"Outputdefinition18")==0) return Outputdefinition18Enum;
@@ -1120 +1121 @@
               else if (strcmp(name,"Regular")==0) return RegularEnum;
               else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
-              else if (strcmp(name,"Scaled")==0) return ScaledEnum;
          else stage=10;
    }
    if(stage==10){
-              if (strcmp(name,"SealevelAbsolute")==0) return SealevelAbsoluteEnum;
+              if (strcmp(name,"Scaled")==0) return ScaledEnum;
+              else if (strcmp(name,"SealevelAbsolute")==0) return SealevelAbsoluteEnum;
               else if (strcmp(name,"SealevelEmotion")==0) return SealevelEmotionEnum;
               else if (strcmp(name,"SealevelInertiaTensorXZ")==0) return SealevelInertiaTensorXZEnum;

issm/trunk-jpl/src/m/classes/matdamageice.m

@@ -14 +14 @@
                 thermalconductivity        = 0.;
                 temperateiceconductivity   = 0.;
+                effectiveconductivity_averaging = 0.;
                 meltingpoint               = 0.;
                 beta                       = 0.;
@@ -81 +82 @@
                         self.temperateiceconductivity=.24;
 
+                        %computation of effective conductivity
+                        self.effectiveconductivity_averaging=2;
+            
                         %the melting point of ice at 1 atmosphere of pressure in K
                         self.meltingpoint=273.15;
@@ -115 +119 @@
                         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'});
-
+                        md = checkfield(md,'fieldname','materials.effectiveconductivity_averaging','numel',[1],'values',[0 1 2]);
+            
                         if ismember('GiaAnalysis',analyses),
                                 md = checkfield(md,'fieldname','materials.lithosphere_shear_modulus','>',0,'numel',1);
@@ -138 +143 @@
                         fielddisplay(self,'thermalconductivity',['ice thermal conductivity [W/m/K]']);
                         fielddisplay(self,'temperateiceconductivity','temperate ice thermal conductivity [W/m/K]');
+                        fielddisplay(self,'effectiveconductivity_averaging','computation of effective conductivity: (0) arithmetic mean, (1) harmonic mean, (2) geometric mean (default)');
                         fielddisplay(self,'meltingpoint','melting point of ice at 1atm in K');
                         fielddisplay(self,'latentheat','latent heat of fusion [J/kg]');
@@ -162 +168 @@
                         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','effectiveconductivity_averaging','format','Integer');
                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','meltingpoint','format','Double');
                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','beta','format','Double');

issm/trunk-jpl/src/m/classes/matdamageice.py

@@ -21 +21 @@
                 self.thermalconductivity       = 0.
                 self.temperateiceconductivity  = 0.
+                self.effectiveconductivity_averaging = 0.
                 self.meltingpoint              = 0.
                 self.beta                      = 0.
@@ -51 +52 @@
                 string="%s\n%s"%(string,fielddisplay(self,"thermalconductivity","ice thermal conductivity [W/m/K]"))
                 string="%s\n%s"%(string,fielddisplay(self,"temperateiceconductivity","temperate ice thermal conductivity [W/m/K]"))
+                string="%s\n%s"%(string,fielddisplay(self,"effectiveconductivity_averaging","computation of effective conductivity: (0) arithmetic mean, (1) harmonic mean, (2) geometric mean (default)")
                 string="%s\n%s"%(string,fielddisplay(self,"meltingpoint","melting point of ice at 1atm in K"))
                 string="%s\n%s"%(string,fielddisplay(self,"latentheat","latent heat of fusion [J/m^3]"))
@@ -98 +100 @@
                 self.temperateiceconductivity=0.24
 
+                #computation of effective conductivity
+                self.effectiveconductivity_averaging=2
+        
                 #the melting point of ice at 1 atmosphere of pressure in K
                 self.meltingpoint=273.15
@@ -139 +144 @@
                 md = checkfield(md,'fieldname','materials.mantle_density','>',0,'numel',[1]);
                 md = checkfield(md,'fieldname','materials.earth_density','>',0,'numel',[1]);
+                md = checkfield(md,'fieldname','materials.effectiveconductivity_averaging','numel',[1],'values',[0 1 2])
 
                 return md
@@ -152 +158 @@
                 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','effectiveconductivity_averaging','format','Integer');
                 WriteData(fid,prefix,'object',self,'class','materials','fieldname','meltingpoint','format','Double')
                 WriteData(fid,prefix,'object',self,'class','materials','fieldname','beta','format','Double')

issm/trunk-jpl/src/m/classes/matenhancedice.m

@@ -14 +14 @@
                 thermalconductivity        = 0.;
                 temperateiceconductivity   = 0.;
+                effectiveconductivity_averaging = 0.;
                 meltingpoint               = 0.;
                 beta                       = 0.;
@@ -83 +84 @@
                         self.temperateiceconductivity=.24;
 
+                        %computation of effective conductivity
+                        self.effectiveconductivity_averaging=2;
+            
                         %the melting point of ice at 1 atmosphere of pressure in K
                         self.meltingpoint=273.15;
@@ -118 +122 @@
                         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'});
-
+                        md = checkfield(md,'fieldname','materials.effectiveconductivity_averaging','numel',[1],'values',[0 1 2]);
+            
                         if ismember('GiaAnalysis',analyses),
                                 md = checkfield(md,'fieldname','materials.lithosphere_shear_modulus','>',0,'numel',1);
@@ -140 +145 @@
                         fielddisplay(self,'thermalconductivity',['ice thermal conductivity [W/m/K]']);
                         fielddisplay(self,'temperateiceconductivity','temperate ice thermal conductivity [W/m/K]');
+                        fielddisplay(self,'effectiveconductivity_averaging','computation of effective conductivity: (0) arithmetic mean, (1) harmonic mean, (2) geometric mean (default)');
                         fielddisplay(self,'meltingpoint','melting point of ice at 1atm in K');
                         fielddisplay(self,'latentheat','latent heat of fusion [J/kg]');
@@ -165 +171 @@
                         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','effectiveconductivity_averaging','format','Integer');
                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','meltingpoint','format','Double');
                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','beta','format','Double');
@@ -190 +197 @@
                         writejsdouble(fid,[modelname '.materials.thermalconductivity'],self.thermalconductivity);
                         writejsdouble(fid,[modelname '.materials.temperateiceconductivity'],self.temperateiceconductivity);
+                        writejsdouble(fid,[modelname '.materials.effectiveconductivity_averaging'],self.effectiveconductivity_averaging);
                         writejsdouble(fid,[modelname '.materials.meltingpoint'],self.meltingpoint);
                         writejsdouble(fid,[modelname '.materials.beta'],self.beta);

issm/trunk-jpl/src/m/classes/matenhancedice.py

@@ -21 +21 @@
                 self.thermalconductivity       = 0.
                 self.temperateiceconductivity  = 0.
+                self.effectiveconductivity_averaging = 0.
                 self.meltingpoint              = 0.
                 self.beta                      = 0.
@@ -51 +52 @@
                 string="%s\n%s"%(string,fielddisplay(self,"thermalconductivity","ice thermal conductivity [W/m/K]"))
                 string="%s\n%s"%(string,fielddisplay(self,"temperateiceconductivity","temperate ice thermal conductivity [W/m/K]"))
+        string="%s\n%s"%(string,fielddisplay(self,"effectiveconductivity_averaging","computation of effective conductivity: (0) arithmetic mean, (1) harmonic mean, (2) geometric mean (default)")
                 string="%s\n%s"%(string,fielddisplay(self,"meltingpoint","melting point of ice at 1atm in K"))
                 string="%s\n%s"%(string,fielddisplay(self,"latentheat","latent heat of fusion [J/m^3]"))
@@ -98 +100 @@
                 #temperate ice thermal conductivity (W/m/K)
                 self.temperateiceconductivity=0.24
+        
+                #computation of effective conductivity
+                self.effectiveconductivity_averaging=2
 
                 #the melting point of ice at 1 atmosphere of pressure in K
@@ -135 +140 @@
                 md = checkfield(md,'fieldname','materials.rheology_n','>',0,'size',[md.mesh.numberofelements])
                 md = checkfield(md,'fieldname','materials.rheology_law','values',['None','BuddJacka','Cuffey','CuffeyTemperate','Paterson','Arrhenius','LliboutryDuval'])
-
+                md = checkfield(md,'fieldname','materials.effectiveconductivity_averaging','numel',[1],'values',[0 1 2])
+        
                 if 'GiaAnalysis' in analyses:
                         md = checkfield(md,'fieldname','materials.lithosphere_shear_modulus','>',0,'numel',1)
@@ -155 +161 @@
                 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','effectiveconductivity_averaging','format','Integer');
                 WriteData(fid,prefix,'object',self,'class','materials','fieldname','meltingpoint','format','Double')
                 WriteData(fid,prefix,'object',self,'class','materials','fieldname','beta','format','Double')

issm/trunk-jpl/src/m/classes/matestar.m

@@ -14 +14 @@
                 thermalconductivity        = 0.;
                 temperateiceconductivity   = 0.;
+                effectiveconductivity_averaging = 0.;
                 meltingpoint               = 0.;
                 beta                       = 0.;
@@ -90 +91 @@
                         self.temperateiceconductivity=.24;
 
+                        %computation of effective conductivity
+                        self.effectiveconductivity_averaging=2;
+            
                         %the melting point of ice at 1 atmosphere of pressure in K
                         self.meltingpoint=273.15;
@@ -125 +129 @@
                         md = checkfield(md,'fieldname','materials.rheology_Es','>',0,'size',[md.mesh.numberofvertices 1],'NaN',1,'Inf',1);
                         md = checkfield(md,'fieldname','materials.rheology_law','values',{'None' 'BuddJacka' 'Cuffey' 'CuffeyTemperate' 'Paterson' 'Arrhenius' 'LliboutryDuval'});
-
+                        md = checkfield(md,'fieldname','materials.effectiveconductivity_averaging','numel',[1],'values',[0 1 2]);
+            
                         if ismember('GiaAnalysis',analyses),
                                 md = checkfield(md,'fieldname','materials.lithosphere_shear_modulus','>',0,'numel',1);
@@ -147 +152 @@
                         fielddisplay(self,'thermalconductivity',['ice thermal conductivity [W/m/K]']);
                         fielddisplay(self,'temperateiceconductivity','temperate ice thermal conductivity [W/m/K]');
+                        fielddisplay(self,'effectiveconductivity_averaging','computation of effective conductivity: (0) arithmetic mean, (1) harmonic mean, (2) geometric mean (default)');
                         fielddisplay(self,'meltingpoint','melting point of ice at 1atm in K');
                         fielddisplay(self,'latentheat','latent heat of fusion [J/kg]');
@@ -172 +178 @@
                         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','effectiveconductivity_averaging','format','Integer');
                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','meltingpoint','format','Double');
                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','beta','format','Double');
@@ -197 +204 @@
                         writejsdouble(fid,[modelname '.materials.thermalconductivity'],self.thermalconductivity);
                         writejsdouble(fid,[modelname '.materials.temperateiceconductivity'],self.temperateiceconductivity);
+                        writejsdouble(fid,[modelname '.materials.effectiveconductivity_averaging'],self.effectiveconductivity_averaging);
                         writejsdouble(fid,[modelname '.materials.meltingpoint'],self.meltingpoint);
                         writejsdouble(fid,[modelname '.materials.beta'],self.beta);

issm/trunk-jpl/src/m/classes/matestar.py

@@ -23 +23 @@
                 thermalconductivity        = 0.
                 temperateiceconductivity   = 0.
+                self.effectiveconductivity_averaging = 0.
                 meltingpoint               = 0.
                 beta                       = 0.
@@ -54 +55 @@
                 string="%s\n%s"%(string,fielddisplay(self,'thermalconductivity',['ice thermal conductivity [W/m/K]']))
                 string="%s\n%s"%(string,fielddisplay(self,'temperateiceconductivity','temperate ice thermal conductivity [W/m/K]'))
+                string="%s\n%s"%(string,fielddisplay(self,"effectiveconductivity_averaging","computation of effective conductivity: (0) arithmetic mean, (1) harmonic mean, (2) geometric mean (default)")
                 string="%s\n%s"%(string,fielddisplay(self,'meltingpoint','melting point of ice at 1atm in K'))
                 string="%s\n%s"%(string,fielddisplay(self,'latentheat','latent heat of fusion [J/kg]'))
@@ -102 +104 @@
                 self.temperateiceconductivity=.24
 
+                #computation of effective conductivity
+                self.effectiveconductivity_averaging=2
+        
                 #the melting point of ice at 1 atmosphere of pressure in K
                 self.meltingpoint=273.15
@@ -138 +143 @@
                 md = checkfield(md,'fieldname','materials.rheology_Es','>',0,'size',[md.mesh.numberofvertices],'NaN',1,'Inf',1)
                 md = checkfield(md,'fieldname','materials.rheology_law','values',['None','BuddJacka', 'Cuffey','CuffeyTemperate','Paterson','Arrhenius','LliboutryDuval'])
+                md = checkfield(md,'fieldname','materials.effectiveconductivity_averaging','numel',[1],'values',[0 1 2])
 
                 if 'GiaAnalysis' in analyses:
@@ -159 +165 @@
                 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','effectiveconductivity_averaging','format','Integer');
                 WriteData(fid,prefix,'object',self,'class','materials','fieldname','meltingpoint','format','Double')
                 WriteData(fid,prefix,'object',self,'class','materials','fieldname','beta','format','Double')

issm/trunk-jpl/src/m/classes/matice.m

@@ -14 +14 @@
                 thermalconductivity        = 0.;
                 temperateiceconductivity   = 0.;
+                effectiveconductivity_averaging = 0.;
                 meltingpoint               = 0.;
                 beta                       = 0.;
@@ -80 +81 @@
                         %wet ice thermal conductivity (W/m/K)
                         self.temperateiceconductivity=.24;
+                        
+                        %computation of effective conductivity
+                        self.effectiveconductivity_averaging=2;
 
                         %the melting point of ice at 1 atmosphere of pressure in K
@@ -115 +119 @@
                         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'});
+                        md = checkfield(md,'fieldname','materials.effectiveconductivity_averaging','numel',[1],'values',[0 1 2]);
 
                         if ismember('GiaAnalysis',analyses),
@@ -137 +142 @@
                         fielddisplay(self,'thermalconductivity',['ice thermal conductivity [W/m/K]']);
                         fielddisplay(self,'temperateiceconductivity','temperate ice thermal conductivity [W/m/K]');
+                        fielddisplay(self,'effectiveconductivity_averaging','computation of effective conductivity: (0) arithmetic mean, (1) harmonic mean, (2) geometric mean (default)');
                         fielddisplay(self,'meltingpoint','melting point of ice at 1atm in K');
                         fielddisplay(self,'latentheat','latent heat of fusion [J/kg]');
@@ -161 +167 @@
                         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','effectiveconductivity_averaging','format','Integer');
                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','meltingpoint','format','Double');
                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','beta','format','Double');
@@ -185 +192 @@
                         writejsdouble(fid,[modelname '.materials.thermalconductivity'],self.thermalconductivity);
                         writejsdouble(fid,[modelname '.materials.temperateiceconductivity'],self.temperateiceconductivity);
+                        writejsdouble(fid,[modelname '.materials.effectiveconductivity_averaging'],self.effectiveconductivity_averaging);
                         writejsdouble(fid,[modelname '.materials.meltingpoint'],self.meltingpoint);
                         writejsdouble(fid,[modelname '.materials.beta'],self.beta);

issm/trunk-jpl/src/m/classes/matice.py

@@ -21 +21 @@
                 self.thermalconductivity       = 0.
                 self.temperateiceconductivity  = 0.
+                self.effectiveconductivity_averaging = 0.
                 self.meltingpoint              = 0.
                 self.beta                      = 0.
@@ -52 +53 @@
                 string="%s\n%s"%(string,fielddisplay(self,"thermalconductivity","ice thermal conductivity [W/m/K]"))
                 string="%s\n%s"%(string,fielddisplay(self,"temperateiceconductivity","temperate ice thermal conductivity [W/m/K]"))
+                string="%s\n%s"%(string,fielddisplay(self,"effectiveconductivity_averaging","computation of effective conductivity: (0) arithmetic mean, (1) harmonic mean, (2) geometric mean (default)")
                 string="%s\n%s"%(string,fielddisplay(self,"meltingpoint","melting point of ice at 1atm in K"))
                 string="%s\n%s"%(string,fielddisplay(self,"latentheat","latent heat of fusion [J/m^3]"))
@@ -99 +101 @@
                 self.temperateiceconductivity=0.24
 
+                #computation of effective conductivity
+                self.effectiveconductivity_averaging=2
+
                 #the melting point of ice at 1 atmosphere of pressure in K
                 self.meltingpoint=273.15
@@ -140 +145 @@
                 md = checkfield(md,'fieldname','materials.mantle_density','>',0,'numel',[1]);
                 md = checkfield(md,'fieldname','materials.earth_density','>',0,'numel',[1]);
+                md = checkfield(md,'fieldname','materials.effectiveconductivity_averaging','numel',[1],'values',[0 1 2])
 
                 return md
@@ -153 +159 @@
                 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','effectiveconductivity_averaging','format','Integer');       
                 WriteData(fid,prefix,'object',self,'class','materials','fieldname','meltingpoint','format','Double')
                 WriteData(fid,prefix,'object',self,'class','materials','fieldname','beta','format','Double')