Changeset 16590

Timestamp:

10/31/13 02:40:38 (11 years ago)

Author:

jbondzio

Message:

ADD: thermal conductivity of temperate ice "materials.temperateiceconductivity", minor bugs

Location:

issm/trunk-jpl

Files:

• src/c/classes/Elements/Penta.cpp (modified) (3 diffs)
• src/c/classes/Materials/Matpar.cpp (modified) (5 diffs)
• src/c/classes/Materials/Matpar.h (modified) (2 diffs)
• src/c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• src/c/shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• src/c/shared/Enum/StringToEnumx.cpp (modified) (4 diffs)
• src/m/classes/matice.m (modified) (4 diffs)
• src/m/classes/matice.py (modified) (4 diffs)
• src/m/enum/EnumDefinitions.py (modified) (1 diff)
• src/m/enum/MaterialsTemperateiceconductivityEnum.m (added)
• test/Archives/Archive120.nc (modified) ( previous)
• test/Archives/Archive121.nc (modified) ( previous)
• test/Archives/Archive122.nc (modified) ( previous)
• test/Archives/Archive325.nc (modified) ( previous)
• test/Archives/Archive326.nc (modified) ( previous)
• test/Archives/Archive327.nc (modified) ( previous)
• test/Archives/Archive431.nc (modified) ( previous)
• test/Archives/Archive432.nc (modified) ( previous)

issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

@@ -3815 +3815 @@
 
         /* Start  looping on the number of gaussian points: */
-        gauss=new GaussPenta(2,2);
+        gauss=new GaussPenta(3,3);
         for(int ig=gauss->begin();ig<gauss->end();ig++){
 
@@ -3828 +3828 @@
                 pressure_input->GetInputValue(&pressure, gauss);
                 kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure); _assert_(kappa>0.);
+
                 D_scalar_conduct=gauss->weight*Jdet*kappa/rho_ice;
                 if(reCast<bool,IssmDouble>(dt)) D_scalar_conduct=D_scalar_conduct*dt;
@@ -3895 +3896 @@
                 else if(stabilization==2){
                         GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
-                        tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,kappa);
+                        tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,kappa/rho_ice);
 
                         for(i=0;i<numdof;i++){

issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp

@@ -34 +34 @@
         iomodel->Constant(&this->heatcapacity,MaterialsHeatcapacityEnum);
         iomodel->Constant(&this->thermalconductivity,MaterialsThermalconductivityEnum);
+        iomodel->Constant(&this->temperateiceconductivity,MaterialsTemperateiceconductivityEnum);
         iomodel->Constant(&this->latentheat,MaterialsLatentheatEnum);
         iomodel->Constant(&this->beta,MaterialsBetaEnum);
@@ -100 +101 @@
         _printf_("   heatcapacity: " << heatcapacity << "\n");
         _printf_("   thermalconductivity: " << thermalconductivity << "\n");
+        _printf_("   temperateiceconductivity: " << temperateiceconductivity << "\n");
         _printf_("   latentheat: " << latentheat << "\n");
         _printf_("   beta: " << beta << "\n");
@@ -169 +171 @@
                         this->heatcapacity=constant;
                         break;
-                case MaterialsThermalconductivityEnum:
+                case MaterialsThermalconductivityEnum:
                         this->thermalconductivity=constant;
+                        break;
+                case MaterialsTemperateiceconductivityEnum:
+                        this->temperateiceconductivity=constant;
                         break;
                 case  MaterialsLatentheatEnum:
@@ -306 +311 @@
 IssmDouble Matpar::GetThermalConductivity(){
         return thermalconductivity;
+}
+/*}}}*/
+/*FUNCTION Matpar::GetTemperateIceConductivity {{{*/
+IssmDouble Matpar::GetTemperateIceConductivity(){
+        return temperateiceconductivity;
 }
 /*}}}*/
@@ -383 +393 @@
         */
 
-        IssmDouble eps=0.05*heatcapacity;
+        IssmDouble eps=0.1*heatcapacity;
         IssmDouble hpmp=PureIceEnthalpy(pressure);
         IssmDouble kappa_c=thermalconductivity/heatcapacity;
-        IssmDouble kappa_t=thermalconductivity/heatcapacity*pow(10.,-1);
-        return 1./(1.+exp(-(enthalpy-(hpmp))/eps))*(kappa_t-kappa_c) + kappa_c;
+        IssmDouble kappa_t=temperateiceconductivity/heatcapacity;
+        
+        if(enthalpy<=hpmp-eps)
+                return kappa_c;
+        else if(enthalpy>=hpmp+eps)
+                return kappa_t;
+        else {
+                IssmDouble xi=enthalpy-hpmp;
+                IssmDouble pi=3.141592653589793238462643;
+                return kappa_c + (kappa_t-kappa_c)*((xi+eps)/(2*eps) + sin(pi*xi/eps)/(2*pi));
+        }
+
+        //return 1./(1.+exp(-(enthalpy-(hpmp))/eps))*(kappa_t-kappa_c) + kappa_c;
 }
 /*}}}*/

issm/trunk-jpl/src/c/classes/Materials/Matpar.h

@@ -22 +22 @@
                 IssmDouble  heatcapacity;
                 IssmDouble  thermalconductivity;
+                IssmDouble  temperateiceconductivity;
                 IssmDouble  latentheat;
                 IssmDouble  beta;
@@ -109 +110 @@
                 IssmDouble GetHeatCapacity();
                 IssmDouble GetThermalConductivity();
+                IssmDouble GetTemperateIceConductivity();
                 IssmDouble GetLatentHeat();
                 IssmDouble GetBeta();

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

@@ -178 +178 @@
         MaterialsThermalExchangeVelocityEnum,
         MaterialsThermalconductivityEnum,
+  MaterialsTemperateiceconductivityEnum,
         MaterialsLithosphereShearModulusEnum,
         MaterialsLithosphereDensityEnum,

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

@@ -186 +186 @@
                 case MaterialsThermalExchangeVelocityEnum : return "MaterialsThermalExchangeVelocity";
                 case MaterialsThermalconductivityEnum : return "MaterialsThermalconductivity";
+                case MaterialsTemperateiceconductivityEnum : return "MaterialsTemperateiceconductivity";
                 case MaterialsLithosphereShearModulusEnum : return "MaterialsLithosphereShearModulus";
                 case MaterialsLithosphereDensityEnum : return "MaterialsLithosphereDensity";

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

@@ -189 +189 @@
               else if (strcmp(name,"MaterialsThermalExchangeVelocity")==0) return MaterialsThermalExchangeVelocityEnum;
               else if (strcmp(name,"MaterialsThermalconductivity")==0) return MaterialsThermalconductivityEnum;
+              else if (strcmp(name,"MaterialsTemperateiceconductivity")==0) return MaterialsTemperateiceconductivityEnum;
               else if (strcmp(name,"MaterialsLithosphereShearModulus")==0) return MaterialsLithosphereShearModulusEnum;
               else if (strcmp(name,"MaterialsLithosphereDensity")==0) return MaterialsLithosphereDensityEnum;
@@ -259 +260 @@
               else if (strcmp(name,"SurfaceforcingsMassBalance")==0) return SurfaceforcingsMassBalanceEnum;
               else if (strcmp(name,"SurfaceforcingsIspdd")==0) return SurfaceforcingsIspddEnum;
-              else if (strcmp(name,"SurfaceforcingsDesfac")==0) return SurfaceforcingsDesfacEnum;
          else stage=3;
    }
    if(stage==3){
-              if (strcmp(name,"SurfaceforcingsS0p")==0) return SurfaceforcingsS0pEnum;
+              if (strcmp(name,"SurfaceforcingsDesfac")==0) return SurfaceforcingsDesfacEnum;
+              else if (strcmp(name,"SurfaceforcingsS0p")==0) return SurfaceforcingsS0pEnum;
               else if (strcmp(name,"SurfaceforcingsIssmbgradients")==0) return SurfaceforcingsIssmbgradientsEnum;
               else if (strcmp(name,"SurfaceforcingsMonthlytemperatures")==0) return SurfaceforcingsMonthlytemperaturesEnum;
@@ -382 +383 @@
               else if (strcmp(name,"IntInput")==0) return IntInputEnum;
               else if (strcmp(name,"InputToExtrude")==0) return InputToExtrudeEnum;
-              else if (strcmp(name,"InputToL2Project")==0) return InputToL2ProjectEnum;
          else stage=4;
    }
    if(stage==4){
-              if (strcmp(name,"IntParam")==0) return IntParamEnum;
+              if (strcmp(name,"InputToL2Project")==0) return InputToL2ProjectEnum;
+              else if (strcmp(name,"IntParam")==0) return IntParamEnum;
               else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum;
               else if (strcmp(name,"TransientParam")==0) return TransientParamEnum;
@@ -505 +506 @@
               else if (strcmp(name,"StressTensorxx")==0) return StressTensorxxEnum;
               else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum;
-              else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum;
          else stage=5;
    }
    if(stage==5){
-              if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum;
+              if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum;
+              else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum;
               else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum;
               else if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum;

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

@@ -13 +13 @@
                 latentheat                 = 0.;
                 thermalconductivity        = 0.;
+                temperateiceconductivity         = 0.;
                 meltingpoint               = 0.;
                 beta                       = 0.;
@@ -69 +70 @@
                         %ice thermal conductivity (W/m/K)
                         obj.thermalconductivity=2.4;
+                        
+                        %wet ice thermal conductivity (W/m/K)
+                        obj.temperateiceconductivity=.24;
 
                         %the melting point of ice at 1 atmosphere of pressure in K
@@ -118 +122 @@
                         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,'thermalconductivity',['ice thermal conductivity [W/m/K]']);
+                        fielddisplay(obj,'temperateiceconductivity','temperate 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]');
@@ -141 +146 @@
                         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','temperateiceconductivity','format','Double');
                         WriteData(fid,'object',obj,'class','materials','fieldname','meltingpoint','format','Double');
                         WriteData(fid,'object',obj,'class','materials','fieldname','beta','format','Double');

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

@@ -21 +21 @@
                 self.latentheat                = 0.
                 self.thermalconductivity       = 0.
+                self.temperateiceconductivity  = 0.
                 self.meltingpoint              = 0.
                 self.beta                      = 0.
@@ -46 +47 @@
                 string="%s\n%s"%(string,fielddisplay(self,"heatcapacity","heat capacity [J/kg/K]"))
                 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,"meltingpoint","melting point of ice at 1atm in K"))
                 string="%s\n%s"%(string,fielddisplay(self,"latentheat","latent heat of fusion [J/m^3]"))
@@ -82 +84 @@
                 #ice thermal conductivity (W/m/K)
                 self.thermalconductivity=2.4
+
+                #temperate ice thermal conductivity (W/m/K)
+                self.temperateiceconductivity=0.24
 
                 #the melting point of ice at 1 atmosphere of pressure in K
@@ -130 +135 @@
                 WriteData(fid,'object',self,'class','materials','fieldname','latentheat','format','Double')
                 WriteData(fid,'object',self,'class','materials','fieldname','thermalconductivity','format','Double')
+                WriteData(fid,'object',self,'class','materials','fieldname','temperateiceconductivity','format','Double')
                 WriteData(fid,'object',self,'class','materials','fieldname','meltingpoint','format','Double')
                 WriteData(fid,'object',self,'class','materials','fieldname','beta','format','Double')

issm/trunk-jpl/src/m/enum/EnumDefinitions.py

@@ -178 +178 @@
 def MaterialsThermalExchangeVelocityEnum(): return StringToEnum("MaterialsThermalExchangeVelocity")[0]
 def MaterialsThermalconductivityEnum(): return StringToEnum("MaterialsThermalconductivity")[0]
+def MaterialsTemperateiceconductivityEnum(): return StringToEnum("MaterialsTemperateiceconductivity")[0]
 def MaterialsLithosphereShearModulusEnum(): return StringToEnum("MaterialsLithosphereShearModulus")[0]
 def MaterialsLithosphereDensityEnum(): return StringToEnum("MaterialsLithosphereDensity")[0]