Changeset 18659

Timestamp:

10/20/14 09:14:10 (10 years ago)

Author:

seroussi

Message:

BUG: fixed spc of basal conditions using a vector

Location:

issm/trunk-jpl/src/c/analyses

Files:

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

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

@@ -1189 +1189 @@
 }/*}}}*/
 void EnthalpyAnalysis::UpdateBasalConstraints(FemModel* femmodel){/*{{{*/
+
         /*Update basal dirichlet BCs for enthalpy: */
+        Vector<IssmDouble>* spc           = NULL;
+        IssmDouble*         serial_spc    = NULL;
+
+        spc=new Vector<IssmDouble>(femmodel->nodes->NumberOfNodes(EnthalpyAnalysisEnum));
+        /*First create a vector to figure out what elements should be constrained*/
         for(int i=0;i<femmodel->elements->Size();i++){
                 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
-                UpdateBasalConstraints(element);
-        }
+                GetBasalConstraints(spc,element);
+        }
+
+        /*Assemble and serialize*/
+        spc->Assemble();
+        serial_spc=spc->ToMPISerial();
+        delete spc;
+
+        /*Then update basal constraints nodes accordingly*/
+        for(int i=0;i<femmodel->elements->Size();i++){
+                Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+                ApplyBasalConstraints(serial_spc,element);
+        }
+
         femmodel->UpdateConstraintsx();
-}/*}}}*/
-void EnthalpyAnalysis::UpdateBasalConstraints(Element* element){/*{{{*/
+
+        /*Delete*/
+        xDelete<IssmDouble>(serial_spc);
+}/*}}}*/
+void EnthalpyAnalysis::ApplyBasalConstraints(IssmDouble* serial_spc,Element* element){/*{{{*/
+
+        /* Check if ice in element */
+        if(!element->IsIceInElement()) return;
+
+        /* Only update Constraints at the base of grounded ice*/
+        if(!(element->IsOnBase()) || element->IsFloating()) return;
+
+        /*Intermediary*/
+        bool        isdynamicbasalspc;
+        int         numindices;
+        int        *indices = NULL;
+        Node*       node = NULL;
+        IssmDouble      pressure;
+
+        /*Check wether dynamic basal boundary conditions are activated */
+        element->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
+        if(!isdynamicbasalspc) return;
+
+        /*Get parameters and inputs: */
+        Input* pressure_input            = element->GetInput(PressureEnum);                                                      _assert_(pressure_input);
+
+        /*Fetch indices of basal & surface nodes for this finite element*/
+        Penta *penta =  (Penta *) element; // TODO: add Basal-/SurfaceNodeIndices to element.h, and change this to Element*
+        penta->BasalNodeIndices(&numindices,&indices,element->GetElementType());
+
+        GaussPenta* gauss=new GaussPenta();
+        for(int i=0;i<numindices;i++){
+                gauss->GaussNode(element->GetElementType(),indices[i]);
+
+                pressure_input->GetInputValue(&pressure,gauss);
+
+                /*apply or release spc*/
+                node=element->GetNode(indices[i]);
+                if(serial_spc[node->Sid()]==1.){
+                        pressure_input->GetInputValue(&pressure, gauss);
+                        node->ApplyConstraint(0,PureIceEnthalpy(element,pressure));
+                }
+                else                    
+                        node->DofInFSet(0);
+        }
+
+        /*Free ressources:*/
+        xDelete<int>(indices);
+        delete gauss;
+}/*}}}*/
+void EnthalpyAnalysis::GetBasalConstraints(Vector<IssmDouble>* vec_spc,Element* element){/*{{{*/
 
         /* Check if ice in element */
@@ -1214 +1281 @@
         element->FindParam(&dt,TimesteppingTimeStepEnum);
         if(dt==0.){
-                UpdateBasalConstraintsSteadystate(element);
+                GetBasalConstraintsSteadystate(vec_spc,element);
         }
         else{
-                UpdateBasalConstraintsTransient(element);
-        }
-}/*}}}*/
-void EnthalpyAnalysis::UpdateBasalConstraintsTransient(Element* element){/*{{{*/
+                GetBasalConstraintsTransient(vec_spc,element);
+        }
+}/*}}}*/
+void EnthalpyAnalysis::GetBasalConstraintsTransient(Vector<IssmDouble>* vec_spc,Element* element){/*{{{*/
 
         /* Check if ice in element */
@@ -1229 +1296 @@
 
         /*Intermediary*/
-        bool        isdynamicbasalspc,setspc;
         int         numindices, numindicesup, state;
         int        *indices = NULL, *indicesup = NULL;
-        Node*       node = NULL;
         IssmDouble      enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate;
-
-        /*Check wether dynamic basal boundary conditions are activated */
-        element->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
-        if(!isdynamicbasalspc) return;
 
         /*Get parameters and inputs: */
@@ -1266 +1327 @@
                 state=GetThermalBasalCondition(element, enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate);
 
-                setspc=false;
                 switch (state) {
                         case 0:
                                 // cold, dry base: apply basal surface forcing
+                                vec_spc->SetValue(element->nodes[i]->Sid(),0.,INS_VAL);
                                 break;
                         case 1:
                                 // cold, wet base: keep at pressure melting point 
-                                setspc=true;
+                                vec_spc->SetValue(element->nodes[i]->Sid(),1.,INS_VAL);
                                 break;
                         case 2:
                                 // temperate, thin refreezing base: release spc
+                                vec_spc->SetValue(element->nodes[i]->Sid(),0.,INS_VAL);
                                 break;
                         case 3:
                                 // temperate, thin melting base: set spc
-                                setspc=true;
+                                vec_spc->SetValue(element->nodes[i]->Sid(),1.,INS_VAL);
                                 break;
                         case 4:
                                 // temperate, thick melting base: set grad H*n=0
+                                vec_spc->SetValue(element->nodes[i]->Sid(),0.,INS_VAL);
                                 break;
                         default:
@@ -1289 +1352 @@
                 }
 
-                /*apply or release spc*/
-                node=element->GetNode(indices[i]);
-                if(setspc){
-                        pressure_input->GetInputValue(&pressure, gauss);
-                        node->ApplyConstraint(0,PureIceEnthalpy(element,pressure));
-                }
-                else                    
-                        node->DofInFSet(0);
         }
 
@@ -1305 +1360 @@
         delete gaussup;
 }/*}}}*/
-void EnthalpyAnalysis::UpdateBasalConstraintsSteadystate(Element* element){/*{{{*/
+void EnthalpyAnalysis::GetBasalConstraintsSteadystate(Vector<IssmDouble>* vec_spc,Element* element){/*{{{*/
 
         /* Check if ice in element */
@@ -1314 +1369 @@
 
         /*Intermediary*/
-        bool        isdynamicbasalspc,setspc;
         int         numindices, numindicesup, state;
         int        *indices = NULL, *indicesup = NULL;
-        Node*       node = NULL;
         IssmDouble      enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate;
-
-        /*Check wether dynamic basal boundary conditions are activated */
-        element->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
-        if(!isdynamicbasalspc) return;
 
         /*Get parameters and inputs: */
@@ -1349 +1398 @@
 
                 state=GetThermalBasalCondition(element, enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate);
-                setspc=false;
                 switch (state) {
                         case 0:
                                 // cold, dry base: apply basal surface forcing
+                                vec_spc->SetValue(element->nodes[i]->Sid(),0.,INS_VAL);
                                 break;
                         case 1:
                                 // cold, wet base: keep at pressure melting point 
-                                setspc=true;
+                                vec_spc->SetValue(element->nodes[i]->Sid(),1.,INS_VAL);
                                 break;
                         case 2:
                                 // temperate, thin refreezing base: release spc
+                                vec_spc->SetValue(element->nodes[i]->Sid(),0.,INS_VAL);
                                 break;
                         case 3:
                                 // temperate, thin melting base: set spc
-                                setspc=true;
+                                vec_spc->SetValue(element->nodes[i]->Sid(),1.,INS_VAL);
                                 break;
                         case 4:
                                 // temperate, thick melting base: s
-                                setspc=true;
+                                vec_spc->SetValue(element->nodes[i]->Sid(),1.,INS_VAL);
                                 break;
                         default:
                                 _printf0_("     unknown thermal basal state found!");
                 }
-
-                /*apply or release spc*/
-                node=element->GetNode(indices[i]);
-                if(setspc){
-                        pressure_input->GetInputValue(&pressure, gauss);
-                        node->ApplyConstraint(0,PureIceEnthalpy(element,pressure));
-                }
-                else                    
-                        node->DofInFSet(0);
         }
 

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

@@ -48 +48 @@
                 static void DrainWaterfraction(Element* element, IssmDouble* pdrainrate_element);
                 static void UpdateBasalConstraints(FemModel* femmodel);
-                static void UpdateBasalConstraints(Element* element);
-                static void UpdateBasalConstraintsTransient(Element* element);
-                static void UpdateBasalConstraintsSteadystate(Element* element);
+                static void ApplyBasalConstraints(IssmDouble* serial_spc,Element* element);
+                static void GetBasalConstraints(Vector<IssmDouble>* vec_spc,Element* element);
+                static void GetBasalConstraintsTransient(Vector<IssmDouble>* vec_spc,Element* element);
+                static void GetBasalConstraintsSteadystate(Vector<IssmDouble>* vec_spc,Element* element);
                 static int GetThermalBasalCondition(Element* element, IssmDouble enthalpy, IssmDouble enthalpy_up, IssmDouble pressure, IssmDouble pressure_up, IssmDouble watercolumn, IssmDouble meltingrate);
                 static IssmDouble GetWetIceConductivity(Element* element, IssmDouble enthalpy, IssmDouble pressure);