Changeset 16034

Timestamp:

08/30/13 10:24:28 (12 years ago)

Author:

seroussi

Message:

CHG: using GetInputOnVertices in enthalpy modules

File:

• issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (modified) (8 diffs)

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

@@ -4768 +4768 @@
                                 for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
                         }
-                        else{
-                                /*do nothing (water layer acts as insulation)*/ 
-                        }
+                        else{  /*do nothing (water layer acts as insulation)*/  }
                 }
         }
@@ -5308 +5306 @@
         IssmDouble  xyz_list[NUMVERTICES][3];
         IssmDouble  xyz_list_tria[NUMVERTICES2D][3];
-        IssmDouble  heatflux, geothermalflux_value,kappa;
+        IssmDouble  heatflux,kappa;
         IssmDouble  vec_heatflux[3];
         IssmDouble  normal_base[3], d1enthalpy[3];
-        IssmDouble  basalmeltingrate[NUMVERTICES], watercolumn[NUMVERTICES], meltingrate_enthalpy;
+        IssmDouble  basalmeltingrate[NUMVERTICES], watercolumn[NUMVERTICES];
         IssmDouble  enthalpy[NUMVERTICES], enthalpyup;
-        IssmDouble  pressure, pressureup;
+        IssmDouble  pressure[NUMVERTICES], pressureup;
         IssmDouble  temperature, waterfraction;
         IssmDouble  latentheat, rho_ice;
         IssmDouble  basalfriction, alpha2;
-        IssmDouble  vx,vy,vz,dt,connectivity;
+        IssmDouble  vx[NUMVERTICES],vy[NUMVERTICES],vz[NUMVERTICES];
+        IssmDouble  geothermalflux[NUMVERTICES];
+        IssmDouble  dt,meltingrate_enthalpy;
         Friction   *friction  = NULL;
 
@@ -5330 +5330 @@
         latentheat=matpar->GetLatentHeat();
         rho_ice=this->matpar->GetRhoIce();
-        Input* watercolumn_input=inputs->GetInput(WatercolumnEnum);                    _assert_(watercolumn_input);
-        Input* basalmeltingrate_input=inputs->GetInput(BasalforcingsMeltingRateEnum);  _assert_(basalmeltingrate_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);
+        GetInputListOnVertices(&vx[0],VxEnum);
+        GetInputListOnVertices(&vy[0],VyEnum);
+        GetInputListOnVertices(&vz[0],VzEnum);
+        GetInputListOnVertices(&enthalpy[0],EnthalpyEnum);
+        GetInputListOnVertices(&pressure[0],PressureEnum);
+        GetInputListOnVertices(&watercolumn[0],WatercolumnEnum);
+        GetInputListOnVertices(&basalmeltingrate[0],BasalforcingsMeltingRateEnum);
+        GetInputListOnVertices(&geothermalflux[0],BasalforcingsGeothermalfluxEnum);
+        Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input);
 
         for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
@@ -5349 +5350 @@
         GaussPenta* gaussup=new GaussPenta();
         for(int iv=0;iv<3;iv++){
+
                 gauss->GaussVertex(iv);
                 gaussup->GaussVertex(iv+3);
-
                 checkpositivethickness=true;
-                watercolumn_input->GetInputValue(&watercolumn[iv], gauss);
-                basalmeltingrate_input->GetInputValue(&basalmeltingrate[iv],gauss);
-                enthalpy_input->GetInputValue(&enthalpy[iv], gauss);
-                pressure_input->GetInputValue(&pressure, gauss);
 
                 /*Calculate basal meltingrate after Fig.5 of A.Aschwanden 2012*/
                 meltingrate_enthalpy=0.;
-                if((watercolumn[iv]>0.) && (enthalpy[iv]<matpar->PureIceEnthalpy(pressure))){
+                if((watercolumn[iv]>0.) && (enthalpy[iv]<matpar->PureIceEnthalpy(pressure[iv]))){
                         /*ensure that no ice is at T<Tm(p), if water layer present*/
-                        enthalpy[iv]=matpar->PureIceEnthalpy(pressure); 
+                        enthalpy[iv]=matpar->PureIceEnthalpy(pressure[iv]); 
                         //meltingrate_enthalpy=0.;
                 }
-                else if(enthalpy[iv]<matpar->PureIceEnthalpy(pressure)){
+                else if(enthalpy[iv]<matpar->PureIceEnthalpy(pressure[iv])){
                         /*cold base: set q*n=q_geo*n+frictionheating as Neumann BC in Penta::CreatePVectorEnthalpySheet*/
                         meltingrate_enthalpy=0.;
@@ -5375 +5372 @@
                         /*ok, from here on all basal ice is temperate. Check for positive thickness of layer of temperate ice. */
                         istemperatelayer=false;
-                        enthalpy_input->GetInputValue(&enthalpyup, gaussup);
-                        pressure_input->GetInputValue(&pressureup, gaussup);    
-                        if(enthalpyup>=matpar->PureIceEnthalpy(pressureup)) istemperatelayer=true;
+                        if(enthalpy[iv+3]>=matpar->PureIceEnthalpy(pressure[iv+3])) istemperatelayer=true;
                         if(istemperatelayer) 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[iv],pressure);
+                                kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy[iv],pressure[iv]);
                                 for(i=0;i<3;i++) vec_heatflux[i]=-kappa*d1enthalpy[i];
                         }
@@ -5389 +5384 @@
                         heatflux=0.;
                         for(i=0;i<3;i++) heatflux+=(vec_heatflux[i])*normal_base[i];
-                        geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
 
                         /*basal friction*/
                         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));
-
-                        matpar->EnthalpyToThermal(&temperature, &waterfraction, enthalpy[iv],pressure);
-                        meltingrate_enthalpy=(basalfriction-(heatflux-geothermalflux_value))/((1-waterfraction)*latentheat*rho_ice); // m/yr water equivalent 
+                        basalfriction=alpha2*(pow(vx[iv],2.0)+pow(vy[iv],2.0)+pow(vz[iv],2.0));
+
+                        matpar->EnthalpyToThermal(&temperature, &waterfraction, enthalpy[iv],pressure[iv]);
+                        meltingrate_enthalpy=(basalfriction-(heatflux-geothermalflux[iv]))/((1-waterfraction)*latentheat*rho_ice); // m/yr water equivalent 
                 }
 
                 /*Update water column, basal meltingrate*/
-                connectivity=IssmDouble(vertices[iv]->Connectivity());
-                meltingrate_enthalpy/=connectivity; // divide meltingrate_enthalpy by connectivity to address multiple iterations over vertex
                 basalmeltingrate[iv]+=meltingrate_enthalpy;
                 this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
@@ -5431 +5420 @@
         IssmDouble enthalpy[NUMVERTICES], pressure[NUMVERTICES]; 
         IssmDouble latentheat, dt;
-        IssmDouble connectivity;
         GaussPenta* gauss;
 
-        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* 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);
 
         /*Check wether enthalpy is activated*/
@@ -5453 +5441 @@
 
                 /*drain water fraction & update enthalpy*/
-                // TODO: replace connectivity-wise draining by draining once per node per timestep
-                connectivity=IssmDouble(vertices[iv]->Connectivity());
-                waterfraction[iv]-=DrainageFunctionWaterfraction(waterfraction[iv], dt)/connectivity;
+                waterfraction[iv]-=DrainageFunctionWaterfraction(waterfraction[iv], dt);
                 matpar->ThermalToEnthalpy(&enthalpy[iv], temperature[iv], waterfraction[iv], pressure[iv]);        
         }