Changeset 15896

Timestamp:

08/23/13 08:57:45 (12 years ago)

Author:

jbondzio

Message:

BUG: edited bugs in Penta::ComputeBasalMeltrate and Penta::DrainWaterfraction

File:

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

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

@@ -5236 +5236 @@
         /* Intermediaries */
         bool        isenthalpy;
-        int         i,j,ig,analysis_type,numindices,numindicesup;
+        int         i,j,analysis_type;
         IssmDouble  xyz_list[NUMVERTICES][3];
         IssmDouble  xyz_list_tria[NUMVERTICES2D][3];
@@ -5243 +5243 @@
         IssmDouble  normal_base[3], d1enthalpy[3];
         IssmDouble  kappa;
+    IssmDouble  meltrate[3], watercolumn[3];
         IssmDouble  enthalpy, enthalpyup;
         IssmDouble  pressure, pressureup;
@@ -5250 +5251 @@
         IssmDouble  vx,vy,vz;
         IssmDouble  dt;
-        int        *indices   = NULL;
-        int        *indicesup = NULL;
         Friction   *friction  = NULL;
 
@@ -5277 +5276 @@
         friction=new Friction("3d",inputs,matpar,analysis_type);
 
-        /*Fetch indices of basal and surface nodes for this finite element*/
-        BasalNodeIndices(&numindices,&indices,this->element_type);
-        SurfaceNodeIndices(&numindicesup,&indicesup,this->element_type);
-        _assert_(numindices==numindicesup); 
-    
-        IssmDouble  meltrate[numindices], watercolumn[numindices];
-
         /*Ok, get meltrates now from basal conditions*/
         GaussPenta* gauss=new GaussPenta();
         GaussPenta* gaussup=new GaussPenta();
-        for(ig=0;ig<numindices;ig++){
-                gauss->GaussVertex(this->element_type,indices[ig]);
-                gaussup->GaussNode(this->element_type,indicesup[ig]);
+        for(int iv=0;iv<3;iv++){
+                gauss->GaussVertex(iv);
+                gaussup->GaussNode(iv+3);
 
         // TODO: make sure that no node is computed twice (insert mask)
 
-                watercolumn_input->GetInputValue(&watercolumn[indices[ig]], gauss);
+                watercolumn_input->GetInputValue(&watercolumn[iv], gauss);
                 enthalpy_input->GetInputValue(&enthalpy, gauss);
                 pressure_input->GetInputValue(&pressure, gauss);
 
                 /*Calculate basal meltrate*/
-                if((watercolumn[indices[ig]]>0.) && (enthalpy<matpar->PureIceEnthalpy(pressure))){
+                if((watercolumn[iv]>0.) && (enthalpy<matpar->PureIceEnthalpy(pressure))){
                         enthalpy=matpar->PureIceEnthalpy(pressure);
                 }
                 else if(enthalpy<matpar->PureIceEnthalpy(pressure)){
-                        meltrate[indices[ig]]=0.;   //TODO: set zero meltrate and watercolumn in model
-                        watercolumn[indices[ig]]=0.;
-                        return;
+                        meltrate[iv]=0.;   
+                        watercolumn[iv]=0.;
+                        continue;
                 }
 
@@ -5333 +5325 @@
                 vz_input->GetInputValue(&vz,gauss);
                 basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0)+pow(vz,2.0));
-                meltrate[indices[ig]]=(basalfriction-(heatflux-geothermalflux_value))/(1-waterfraction)/latentheat; 
+                meltrate[iv]=(basalfriction-(heatflux-geothermalflux_value))/(1-waterfraction)/latentheat; 
 
                 /*Update water column*/
                 this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
                 if(reCast<bool,IssmDouble>(dt))
-            watercolumn[indices[ig]]+=dt*meltrate[indices[ig]];
+            watercolumn[iv]+=dt*meltrate[iv];
                 else
-            watercolumn[indices[ig]]=meltrate[indices[ig]];
+            watercolumn[iv]=meltrate[iv];
         }  
     /*update meltrate and watercolumn*/
@@ -5347 +5339 @@
 
         /*Clean up and return*/
-        xDelete<int>(indices);
-        xDelete<int>(indicesup);
         delete gauss;
         delete gaussup;
@@ -5359 +5349 @@
 
     /*Intermediaries*/
-    int ig;
     bool isenthalpy;
     IssmDouble waterfraction_array[NUMVERTICES], temperature[NUMVERTICES];
@@ -5376 +5365 @@
     this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     latentheat=matpar->GetLatentHeat();
-    gauss=new GaussPenta(2,3);
-    for(ig=0;ig<NUMVERTICES;ig++){ 
-        gauss->GaussVertex(ig);
+    gauss=new GaussPenta();
+    for(int iv=0;iv<NUMVERTICES;iv++){ 
+        gauss->GaussVertex(iv);
         /*TODO: Check whether Vertex has been drained already*/
-        enthalpy_input->GetInputValue(&enthalpy[ig], gauss);
-        pressure_input->GetInputValue(&pressure[ig], gauss);
-        matpar->EnthalpyToThermal(&temperature[ig], &waterfraction[ig], enthalpy[ig],pressure[ig]); 
+        enthalpy_input->GetInputValue(&enthalpy[iv], gauss);
+        pressure_input->GetInputValue(&pressure[iv], gauss);
+        matpar->EnthalpyToThermal(&temperature[iv], &waterfraction[iv], enthalpy[iv],pressure[iv]); 
     
         /*drain water fraction*/
-        waterfraction[ig]-=dt*DrainageFunctionWaterfraction(waterfraction[ig]);
-        if(waterfraction[ig]<0) waterfraction[ig]=0.;
+        waterfraction[iv]-=dt*DrainageFunctionWaterfraction(waterfraction[iv]);
+        if(waterfraction[iv]<0.) waterfraction[iv]=0.;
         /*update enthalpy*/
-        matpar->ThermalToEnthalpy(&enthalpy[ig], temperature[ig], waterfraction[ig], pressure[ig]);        
+        matpar->ThermalToEnthalpy(&enthalpy[iv], temperature[iv], waterfraction[iv], pressure[iv]);        
     }
     /*feed updated results back into model*/