Changeset 15870

Timestamp:

08/22/13 08:20:35 (12 years ago)

Author:

Mathieu Morlighem

Message:

CHG: fixed some stuff from Johannes

Location:

issm/trunk-jpl/src/c/classes

Files:

• Elements/Penta.cpp (modified) (1 diff)
• Materials/Matpar.cpp (modified) (1 diff)

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

@@ -5058 +5058 @@
 /*FUNCTION Penta::ComputeBasalMeltrate{{{*/
 void Penta::ComputeBasalMeltrate(void){
-  /*Calculate the basal melt rates of the enthalpy model after Aschwanden 2012*/
-
-  /*Parameters*/
-  int dim;
-  parameters->FindParam(&dim, MeshDimensionEnum);
-
-  /* Intermediaries */
-  bool isenthalpy;
-  int i, analysis_type;
-  IssmDouble xyz_list[NUMVERTICES][3];
-  IssmDouble heatflux, geothermalflux_value;
-  IssmDouble vec_heatflux[dim];
-  IssmDouble normal_base[dim], d1enthalpy[dim];
-  IssmDouble kappa;
-  IssmDouble enthalpy, enthalpyup;
-  IssmDouble pressure, pressureup;
-  IssmDouble meltrate, watercolumn;
-  IssmDouble temperature, waterfraction;
-  IssmDouble latentheat;
-  IssmDouble basalfriction, alpha2;
-  IssmDouble vx,vy,vz;
-  IssmDouble dt;
-  Friction*  friction=NULL;
-  GaussPenta* gauss=NULL;
-  GaussPenta* gaussup=NULL;
-
-  /*check that dimension is 3*/
-  if(dim!=3) return;
-  /* Only compute melt rates at the base of grounded ice*/
-  if(!IsOnBed() || IsFloating()) return;
-  /*Check wether enthalpy is activated*/
-  parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
-  if(!isenthalpy) return;
-
-  /*Fetch parameters: */
-  latentheat=matpar->GetLatentHeat();
-
-  Input* watercolumn_input=inputs->GetInput(WatercolumnEnum);       _assert_(watercolumn_input);
-  Input* enthalpy_input=inputs->GetInput(EnthalpyEnum);             _assert_(enthalpy_input);
-  Input* pressure_input=inputs->GetInput(PressureEnum);             _assert_(pressure_input);
-  Input* vx_input=inputs->GetInput(VxEnum);                         _assert_(vx_input);
-  Input* vy_input=inputs->GetInput(VyEnum);                         _assert_(vy_input);
-  Input* vz_input=inputs->GetInput(VzEnum);                         _assert_(vz_input);
-  Input* geothermalflux_input=inputs->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
-
-  /*Build friction element, needed later: */
-  parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-  friction=new Friction("3d",inputs,matpar,analysis_type);
-
-  /* Start  looping on the number of gaussian points: */
-  gauss=new GaussPenta(0,1,2,2);
-  gaussup=new GaussPenta(3,4,5,2);
-  for(int ig=gauss->begin();ig<gauss->end();ig++){
-    gauss->GaussPoint(ig);
-    gaussup->GaussPoint(ig);
-    
-    watercolumn_input->GetInputValue(&watercolumn, gauss);
-    enthalpy_input->GetInputValue(&enthalpy, gauss);
-    pressure_input->GetInputValue(&pressure, gauss);
-
-    /*Calculate basal meltrate*/
-    if((watercolumn>0) && (enthalpy<matpar->PureIceEnthalpy(pressure)))
-      enthalpy=matpar->PureIceEnthalpy(pressure);
-    else if(enthalpy<matpar->PureIceEnthalpy(pressure)){
-      meltrate=0.;   //TODO: set zero meltrate and watercolumn in model
-      watercolumn=0.;
-      return;
-    }
-
-    /*ok, from here on all basal ice is temperate. Check for positive thickness of layer of temperate ice) */
-    enthalpy_input->GetInputValue(&enthalpyup, gaussup);
-    pressure_input->GetInputValue(&pressureup, gaussup);    
-    if(enthalpyup>=matpar->PureIceEnthalpy(pressureup))
-      for(i=0;i<dim;i++)
-        vec_heatflux[i]=0.;
-    else{
-      enthalpy_input->GetInputDerivativeValue(&d1enthalpy[0],&xyz_list[0][0],gauss);
-      kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure);
-      for(i=0;i<dim;i++)
-        vec_heatflux[i]=-1.*kappa*d1enthalpy[i];
-    }
-    /*Get normal at base*/
-    normal_base[0]=0.; normal_base[1]=0.; normal_base[2]=1.; // TODO: get normal from model geometry
-    
-    heatflux=0.;
-    for(i=0;i<dim;i++)
-      heatflux+=(vec_heatflux[i])*normal_base[i];
-    geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
-    
-    matpar->EnthalpyToThermal(&temperature, &waterfraction, enthalpy,pressure);
-    
-    friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
-    vx_input->GetInputValue(&vx,gauss);
-    vy_input->GetInputValue(&vy,gauss);
-    vz_input->GetInputValue(&vz,gauss);
-    basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0)+pow(vz,2.0));
-    meltrate=(basalfriction-(heatflux-geothermalflux_value))/(1-waterfraction)/latentheat; 
-    
-    /*Update water column*/
-    this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-    if(reCast<bool,IssmDouble>(dt))
-      watercolumn+=dt*meltrate;
-    else
-      watercolumn=meltrate;
-    // TODO: feed meltrate & watercolumn back to model
-  }  
+        /*Calculate the basal melt rates of the enthalpy model after Aschwanden 2012*/
+
+        /* Intermediaries */
+        bool        isenthalpy;
+        int         i,j,analysis_type;
+        IssmDouble  xyz_list[NUMVERTICES][3];
+        IssmDouble  xyz_list_tria[NUMVERTICES2D][3];
+        IssmDouble  heatflux, geothermalflux_value;
+        IssmDouble  vec_heatflux[3];
+        IssmDouble  normal_base[3], d1enthalpy[3];
+        IssmDouble  kappa;
+        IssmDouble  enthalpy, enthalpyup;
+        IssmDouble  pressure, pressureup;
+        IssmDouble  meltrate, watercolumn;
+        IssmDouble  temperature, waterfraction;
+        IssmDouble  latentheat;
+        IssmDouble  basalfriction, alpha2;
+        IssmDouble  vx,vy,vz;
+        IssmDouble  dt;
+        Friction   *friction                 = NULL;
+        GaussPenta *gauss                    = NULL;
+        GaussPenta *gaussup                  = NULL;
+
+        /* Only compute melt rates at the base of grounded ice*/
+        if(!IsOnBed() || IsFloating()) return;
+
+        /*Check wether enthalpy is activated*/
+        parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
+        if(!isenthalpy) return;
+
+        /*Fetch parameters and inputs */
+        latentheat=matpar->GetLatentHeat();
+        Input* watercolumn_input=inputs->GetInput(WatercolumnEnum);       _assert_(watercolumn_input);
+        Input* enthalpy_input=inputs->GetInput(EnthalpyEnum);             _assert_(enthalpy_input);
+        Input* pressure_input=inputs->GetInput(PressureEnum);             _assert_(pressure_input);
+        Input* vx_input=inputs->GetInput(VxEnum);                         _assert_(vx_input);
+        Input* vy_input=inputs->GetInput(VyEnum);                         _assert_(vy_input);
+        Input* vz_input=inputs->GetInput(VzEnum);                         _assert_(vz_input);
+        Input* geothermalflux_input=inputs->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
+
+        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
+
+        /*Build friction element, needed later: */
+        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+        friction=new Friction("3d",inputs,matpar,analysis_type);
+
+        /* Start  looping on the number of gaussian points: */
+        gauss=new GaussPenta(0,1,2,2);
+        gaussup=new GaussPenta(3,4,5,2);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+                gaussup->GaussPoint(ig);
+
+                watercolumn_input->GetInputValue(&watercolumn, gauss);
+                enthalpy_input->GetInputValue(&enthalpy, gauss);
+                pressure_input->GetInputValue(&pressure, gauss);
+
+                /*Calculate basal meltrate*/
+                if((watercolumn>0.) && (enthalpy<matpar->PureIceEnthalpy(pressure))){
+                        enthalpy=matpar->PureIceEnthalpy(pressure);
+                }
+                else if(enthalpy<matpar->PureIceEnthalpy(pressure)){
+                        meltrate=0.;   //TODO: set zero meltrate and watercolumn in model
+                        watercolumn=0.;
+                        return;
+                }
+
+                /*ok, from here on all basal ice is temperate. Check for positive thickness of layer of temperate ice) */
+                enthalpy_input->GetInputValue(&enthalpyup, gaussup);
+                pressure_input->GetInputValue(&pressureup, gaussup);    
+                if(enthalpyup>=matpar->PureIceEnthalpy(pressureup)){
+                 for(i=0;i<3;i++) vec_heatflux[i]=0.;
+                }
+                else{
+                        enthalpy_input->GetInputDerivativeValue(&d1enthalpy[0],&xyz_list[0][0],gauss);
+                        kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure);
+                        for(i=0;i<3;i++) vec_heatflux[i]=-kappa*d1enthalpy[i];
+                }
+
+                /*Get normal vector to the bed */
+                BedNormal(&normal_base[0],xyz_list_tria);
+
+                heatflux=0.;
+                for(i=0;i<3;i++) heatflux+=(vec_heatflux[i])*normal_base[i];
+                geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
+
+                matpar->EnthalpyToThermal(&temperature, &waterfraction, enthalpy,pressure);
+
+                friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
+                vx_input->GetInputValue(&vx,gauss);
+                vy_input->GetInputValue(&vy,gauss);
+                vz_input->GetInputValue(&vz,gauss);
+                basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0)+pow(vz,2.0));
+                meltrate=(basalfriction-(heatflux-geothermalflux_value))/(1-waterfraction)/latentheat; 
+
+                /*Update water column*/
+                this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+                if(reCast<bool,IssmDouble>(dt))
+                 watercolumn+=dt*meltrate;
+                else
+                 watercolumn=meltrate;
+                // TODO: feed meltrate & watercolumn back to model
+        }  
 }
 /*}}}*/

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

@@ -398 +398 @@
 IssmDouble Matpar::GetEnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){
         if(enthalpy<PureIceEnthalpy(pressure)){
-          return thermalconductivity/(rho_ice*heatcapacity);
+                return thermalconductivity/(rho_ice*heatcapacity);
         }
         else{
-          return 1.045*pow(10,-4.)/rho_ice; // K0=1.045*1e-4 from Aschwanden 2012. TODO: fetch K0 from model
+                return 1.045*1.e-4/rho_ice; // K0=1.045*1e-4 from Aschwanden 2012. TODO: fetch K0 from model
         }
 }