Changeset 15893

Timestamp:

08/23/13 07:26:55 (12 years ago)

Author:

jbondzio

Message:

ADD: Penta::DrainWaterfraction implements drainage of excess water fraction in enthalpy module.

File:

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

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

@@ -5236 +5236 @@
         /* Intermediaries */
         bool        isenthalpy;
-        int         i,j,analysis_type,numindices,numindicesup;
+        int         i,j,ig,analysis_type,numindices,numindicesup;
         IssmDouble  xyz_list[NUMVERTICES][3];
         IssmDouble  xyz_list_tria[NUMVERTICES2D][3];
@@ -5283 +5283 @@
         _assert_(numindices==numindicesup); 
 
-        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
+        /*Ok, get meltrates now from basal conditions*/
         GaussPenta* gauss=new GaussPenta();
         GaussPenta* gaussup=new GaussPenta();
-        for(i=0;i<numindices;i++){
-                gauss->GaussNode(this->element_type,indices[i]);
-                gaussup->GaussNode(this->element_type,indicesup[i]);
+        for(ig=0;ig<numindices;ig++){
+                gauss->GaussNode(this->element_type,indices[ig]);
+                gaussup->GaussNode(this->element_type,indicesup[ig]);
 
                 watercolumn_input->GetInputValue(&watercolumn, gauss);
@@ -5308 +5308 @@
                 pressure_input->GetInputValue(&pressureup, gaussup);    
                 if(enthalpyup>=matpar->PureIceEnthalpy(pressureup)){
-                 for(i=0;i<3;i++) vec_heatflux[i]=0.;
+            for(i=0;i<3;i++) vec_heatflux[i]=0.;
                 }
                 else{
@@ -5346 +5346 @@
         delete gauss;
         delete gaussup;
+}
+/*}}}*/
+/*FUNCTION Penta::DrainWaterfraction{{{*/
+void Penta::DrainWaterfraction(void){
+    
+// TODO: create vector to mark whether node has been drained already i.o. to not drain nodes multiple times
+
+    /*Intermediaries*/
+    const int numdof=NDOF1*NUMVERTICES;
+    int ig;
+    bool isenthalpy;
+    IssmDouble waterfraction, temperature;
+    IssmDouble enthalpy, pressure; 
+    IssmDouble latentheat, dt;
+
+    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*/
+        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); 
+    
+        /*drain water fraction*/
+        waterfraction-=dt*DrainageFunctionWaterfraction(waterfraction);
+        if(waterfraction<0) waterfraction=0.;
+        /*update enthalpy*/
+        latentheat=matpar->GetLatentHeat();
+        matpar->ThermalToEnthalpy(&enthalpy, temperature, waterfraction, pressure);
+        //TODO feed result back into model
+    }
 }
 /*}}}*/