Changeset 15895

Timestamp:

08/23/13 08:36:26 (12 years ago)

Author:

jbondzio

Message:

CHG: updated thermal values now fed back into model.

File:

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

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

@@ -5245 +5245 @@
         IssmDouble  enthalpy, enthalpyup;
         IssmDouble  pressure, pressureup;
-        IssmDouble  meltrate, watercolumn;
         IssmDouble  temperature, waterfraction;
         IssmDouble  latentheat;
@@ -5282 +5281 @@
         SurfaceNodeIndices(&numindicesup,&indicesup,this->element_type);
         _assert_(numindices==numindicesup); 
+    
+        IssmDouble  meltrate[numindices], watercolumn[numindices];
 
         /*Ok, get meltrates now from basal conditions*/
@@ -5287 +5288 @@
         GaussPenta* gaussup=new GaussPenta();
         for(ig=0;ig<numindices;ig++){
-                gauss->GaussNode(this->element_type,indices[ig]);
+                gauss->GaussVertex(this->element_type,indices[ig]);
                 gaussup->GaussNode(this->element_type,indicesup[ig]);
 
-                watercolumn_input->GetInputValue(&watercolumn, gauss);
+        // TODO: make sure that no node is computed twice (insert mask)
+
+                watercolumn_input->GetInputValue(&watercolumn[indices[ig]], gauss);
                 enthalpy_input->GetInputValue(&enthalpy, gauss);
                 pressure_input->GetInputValue(&pressure, gauss);
 
                 /*Calculate basal meltrate*/
-                if((watercolumn>0.) && (enthalpy<matpar->PureIceEnthalpy(pressure))){
+                if((watercolumn[indices[ig]]>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.;
+                        meltrate[indices[ig]]=0.;   //TODO: set zero meltrate and watercolumn in model
+                        watercolumn[indices[ig]]=0.;
                         return;
                 }
@@ -5330 +5333 @@
                 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; 
+                meltrate[indices[ig]]=(basalfriction-(heatflux-geothermalflux_value))/(1-waterfraction)/latentheat; 
 
                 /*Update water column*/
                 this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
                 if(reCast<bool,IssmDouble>(dt))
-                 watercolumn+=dt*meltrate;
+            watercolumn[indices[ig]]+=dt*meltrate[indices[ig]];
                 else
-                 watercolumn=meltrate;
-                // TODO: feed meltrate & watercolumn back to model
+            watercolumn[indices[ig]]=meltrate[indices[ig]];
         }  
+    /*update meltrate and watercolumn*/
+    this->inputs->AddInput(new PentaInput(WatercolumnEnum, watercolumn, P1Enum));
+    this->inputs->AddInput(new PentaInput(BasalMeltrateEnum, meltrate, P1Enum));
 
         /*Clean up and return*/
@@ -5354 +5359 @@
 
     /*Intermediaries*/
-    const int numdof=NDOF1*NUMVERTICES;
     int ig;
     bool isenthalpy;
-    IssmDouble waterfraction, temperature;
-    IssmDouble enthalpy, pressure; 
+    IssmDouble waterfraction_array[NUMVERTICES], temperature[NUMVERTICES];
+    IssmDouble enthalpy[NUMVERTICES], pressure[NUMVERTICES]; 
     IssmDouble latentheat, dt;
-
+    GaussPenta* gauss;
+    
     Input* watercolumn_input=inputs->GetInput(WatercolumnEnum);       _assert_(watercolumn_input);
         Input* enthalpy_input=inputs->GetInput(EnthalpyEnum);             _assert_(enthalpy_input);
@@ -5368 +5373 @@
         parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
         if(!isenthalpy) return;       
+    
     this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-    
-    GaussPenta* gauss=new GaussPenta(2,3);
-    for(ig=gauss->begin();ig<gauss->end();ig++){ 
-        gauss->GaussPoint(ig);
-        
-        enthalpy_input->GetInputValue(&enthalpy, gauss);
-        pressure_input->GetInputValue(&pressure, gauss);
-        matpar->EnthalpyToThermal(&temperature, &waterfraction, enthalpy,pressure); 
+    latentheat=matpar->GetLatentHeat();
+    gauss=new GaussPenta(2,3);
+    for(ig=0;ig<NUMVERTICES;ig++){ 
+        gauss->GaussVertex(ig);
+        /*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]); 
     
         /*drain water fraction*/
-        waterfraction-=dt*DrainageFunctionWaterfraction(waterfraction);
-        if(waterfraction<0) waterfraction=0.;
+        waterfraction[ig]-=dt*DrainageFunctionWaterfraction(waterfraction[ig]);
+        if(waterfraction[ig]<0) waterfraction[ig]=0.;
         /*update enthalpy*/
-        latentheat=matpar->GetLatentHeat();
-        matpar->ThermalToEnthalpy(&enthalpy, temperature, waterfraction, pressure);
-        //TODO feed result back into model
+        matpar->ThermalToEnthalpy(&enthalpy[ig], temperature[ig], waterfraction[ig], pressure[ig]);        
     }
+    /*feed updated results back into model*/
+    this->inputs->AddInput(new PentaInput(EnthalpyEnum,enthalpy,P1Enum));
+    this->inputs->AddInput(new PentaInput(WaterfractionEnum,waterfraction,P1Enum));
+    this->inputs->AddInput(new PentaInput(TemperatureEnum,temperature,P1Enum));    
 }
 /*}}}*/