Changeset 17166

Timestamp:

01/24/14 11:41:17 (11 years ago)

Author:

jbondzio

Message:

CHG: separated computation of basal melting rates and drainage of excess water in ice column. minor cleanups

Location:

issm/trunk-jpl/src/c/analyses

Files:

• EnthalpyAnalysis.cpp (modified) (7 diffs)
• EnthalpyAnalysis.h (modified) (1 diff)

issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp

@@ -821 +821 @@
 
 /*Modules*/
-/*{{{*/
-void EnthalpyAnalysis::PostProcessing(FemModel* femmodel){
+void EnthalpyAnalysis::PostProcessing(FemModel* femmodel){/*{{{*/
+
         /*Intermediaries*/
-        int solution_type;
-
-        /*Compute basal melting rates: */
-        for(int i=0;i<femmodel->elements->Size();i++){
-                Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
-                ComputeBasalMeltingrate(element);
-        }
-
-        /*Update basal dirichlet BCs for enthalpy in transient runs: */
-        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
-        if(solution_type==TransientSolutionEnum){
-                for(int i=0;i<femmodel->elements->Size();i++){
+        int solution_type, i;
+        bool computebasalmeltingrates=true;
+        bool isdrainage=true;
+        bool updatebasalconstraints=true;
+
+        if(isdrainage){
+                /*Drain excess water fraction in ice column: */
+                for(i=0;i<femmodel->elements->Size();i++){
                         Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
-                        UpdateBasalConstraints(element);
-                }
-        }
-}/*}}}*/
-
+                        DrainWaterfractionIcecolumn(element);
+                }
+        }
+
+        if(computebasalmeltingrates){
+                /*Compute basal melting rates: */
+                for(i=0;i<femmodel->elements->Size();i++){
+                        Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+                        ComputeBasalMeltingrate(element);
+                }
+        }
+
+        if(updatebasalconstraints){
+                /*Update basal dirichlet BCs for enthalpy in transient runs: */
+                femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
+                if(solution_type==TransientSolutionEnum){
+                        for(i=0;i<femmodel->elements->Size();i++){
+                                Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+                                UpdateBasalConstraints(element);
+                        }
+                }
+        }
+}/*}}}*/
 void EnthalpyAnalysis::ComputeBasalMeltingrate(Element* element){/*{{{*/
         /*Calculate the basal melt rates of the enthalpy model after Aschwanden 2012*/
@@ -946 +960 @@
                 }
         }
-
-        /******** DRAINAGE *****************************************/
-        IssmDouble* drainrate_column  = xNew<IssmDouble>(numsegments);
-        IssmDouble* drainrate_element = xNew<IssmDouble>(numsegments);
-        for(is=0;is<numsegments;is++)   drainrate_column[is]=0.;
-        Element* elementi = element;
-        for(;;){
-                for(is=0;is<numsegments;is++)   drainrate_element[is]=0.;
-                DrainWaterfraction(elementi,drainrate_element); // TODO: make sure every vertex is only drained once
-                for(is=0;is<numsegments;is++)   drainrate_column[is]+=drainrate_element[is];
-
-                if(elementi->IsOnSurface()) break;
-                elementi=elementi->GetUpperElement();                   
-        }
-        /* add drained water to melting rate*/
-        for(is=0;is<numsegments;is++) meltingrate_enthalpy[is]+=drainrate_column[is];
-
         /******** UPDATE MELTINGRATES AND WATERCOLUMN **************/
         element->FindParam(&dt,TimesteppingTimeStepEnum);
@@ -970 +967 @@
                 if(dt!=0.){
                         if(watercolumn[vertexdown]+meltingrate_enthalpy[is]*dt<0.){     // prevent too much freeze on                   
-                                melting_overshoot=watercolumn[vertexdown]+meltingrate_enthalpy[is]*dt;
-                                lambda=melting_overshoot/(meltingrate_enthalpy[is]*dt); _assert_(lambda>0); _assert_(lambda<1);
-                                basalmeltingrate[vertexdown]=(1.-lambda)*meltingrate_enthalpy[is];
+                                lambda = -watercolumn[vertexdown]/(dt*meltingrate_enthalpy[is]); _assert_(lambda>=0.); _assert_(lambda<1.);
                                 watercolumn[vertexdown]=0.;
-                                yts=365.*24.*60.*60.;
-                                enthalpy[vertexdown]+=dt/yts*lambda*heating[is];
+                                basalmeltingrate[vertexdown]=lambda*meltingrate_enthalpy[is]; // restrict freeze on only to size of watercolumn
+                                enthalpy[vertexdown]+=(1.-lambda)*meltingrate_enthalpy[is]*dt*latentheat; // use rest of energy to cool down base
                         }
                         else{
@@ -1004 +999 @@
         xDelete<IssmDouble>(meltingrate_enthalpy);
         xDelete<IssmDouble>(heating);
+        xDelete<IssmDouble>(xyz_list);
+}/*}}}*/
+void EnthalpyAnalysis::DrainWaterfractionIcecolumn(Element* element){/*{{{*/
+
+        /* Only drain waterfraction of ice column from element at base*/
+        if(!element->IsOnBed()) return; //FIXME: allow freeze on for floating elements
+
+        /* Intermediaries*/
+        int is, numvertices, numsegments;
+        int *pairindices   = NULL;
+
+        numvertices=element->GetNumberOfVertices();
+        element->VerticalSegmentIndices(&pairindices,&numsegments);
+
+        IssmDouble* watercolumn = xNew<IssmDouble>(numvertices);
+        IssmDouble* drainrate_column  = xNew<IssmDouble>(numsegments);
+        IssmDouble* drainrate_element = xNew<IssmDouble>(numsegments);
+
+        element->GetInputListOnVertices(watercolumn,WatercolumnEnum);
+
+        for(is=0;is<numsegments;is++)   drainrate_column[is]=0.;
+        Element* elementi = element;
+        for(;;){
+                for(is=0;is<numsegments;is++)   drainrate_element[is]=0.;
+                DrainWaterfraction(elementi,drainrate_element); // TODO: make sure every vertex is only drained once
+                for(is=0;is<numsegments;is++)   drainrate_column[is]+=drainrate_element[is];
+
+                if(elementi->IsOnSurface()) break;
+                elementi=elementi->GetUpperElement();                   
+        }
+        /* add drained water to water column*/
+        for(is=0;is<numsegments;is++) watercolumn[is]+=drainrate_column[is];
+        /* Feed updated water column back into model */
+        element->AddInput(WatercolumnEnum,watercolumn,P1Enum);
+
+        delete pairindices;
         xDelete<IssmDouble>(drainrate_column);
         xDelete<IssmDouble>(drainrate_element);
-        xDelete<IssmDouble>(xyz_list);
-}/*}}}*/
-
+        xDelete<IssmDouble>(watercolumn);
+}/*}}}*/
 void EnthalpyAnalysis::DrainWaterfraction(Element* element, IssmDouble* pdrainrate_element){/*{{{*/
 
@@ -1056 +1086 @@
                 height_element=fabs(xyz_list[vertexup*3+2]-xyz_list[vertexdown*3+2]);
                 pdrainrate_element[is]=(deltawaterfractions[vertexdown]+deltawaterfractions[vertexup])/2.*height_element; // return water equivalent of drainage
+                _assert_(pdrainrate_element[is]>=0.);
         }
 
@@ -1066 +1097 @@
         xDelete<IssmDouble>(deltawaterfractions);
 }/*}}}*/
-
 void EnthalpyAnalysis::UpdateBasalConstraints(Element* element){/*{{{*/
 
@@ -1120 +1150 @@
                         setspc = false;
 
-
-
                 node=element->GetNode(indices[i]);
                 if (setspc) 
-                        node->ApplyConstraint(1,h_pmp); /*apply spc*/ //nodes[indices[i]]->ApplyConstraint(1,h_pmp);
+                        node->ApplyConstraint(1,h_pmp); /*apply spc*/ 
                 else                    
-                        node->DofInFSet(0); /*remove spc*/ //nodes[indices[i]]->DofInFSet(0);
+                        node->DofInFSet(0); /*remove spc*/ 
         }
 

issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h

@@ -40 +40 @@
                 static void PostProcessing(FemModel* femmodel);
                 static void ComputeBasalMeltingrate(Element* element);
+                static void DrainWaterfractionIcecolumn(Element* element);
                 static void DrainWaterfraction(Element* element, IssmDouble* pdrainrate_element);
                 static void UpdateBasalConstraints(Element* element);