Changeset 15853

Timestamp:

08/21/13 11:27:27 (12 years ago)

Author:

Mathieu Morlighem

Message:

NEW: preparing Dynamic Basal Dirichlet conditions for enthalpy model (not working for now, preparing code for Johannes)

Location:

issm/trunk-jpl/src/c

Files:

• analyses/thermal_core.cpp (modified) (1 diff)
• classes/Elements/Element.h (modified) (1 diff)
• classes/Elements/Penta.cpp (modified) (1 diff)
• classes/Elements/Penta.h (modified) (1 diff)
• classes/Elements/Tria.h (modified) (1 diff)
• classes/FemModel.cpp (modified) (2 diffs)
• classes/FemModel.h (modified) (1 diff)

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

@@ -27 +27 @@
         }
 
+        if(VerboseSolution()) _printf0_("   setting basal Dirichlet boundary conditions\n");
+        femmodel->UpdateThermalBasalConstraintsx();
+
         if(VerboseSolution()) _printf0_("   computing temperatures\n");
         femmodel->SetCurrentConfiguration(ThermalAnalysisEnum);

issm/trunk-jpl/src/c/classes/Elements/Element.h ¶

@@ -127 +127 @@
                 #endif
 
+                #ifdef _HAVE_THERMAL_
+                virtual void UpdateThermalBasalConstraints(void)=0;
+                #endif
+
                 #ifdef _HAVE_HYDROLOGY_
                 virtual void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode)=0;

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

@@ -5003 +5003 @@
         /*Free ressources:*/
         xDelete<int>(doflist);
+}
+/*}}}*/
+/*FUNCTION Penta::UpdateThermalBasalConstraints{{{*/
+void  Penta::UpdateThermalBasalConstraints(void){
+
+        /*Intermediary*/
+        bool        isenthalpy,isdynamicbasalspc,istemperatelayer;
+        int         numindices;
+        IssmDouble  h_pmp,pressure;
+        int        *indices = NULL;
+
+        /* Only update Constraints at the base of grounded ice*/
+        if(!IsOnBed() || !IsFloating()) return;
+
+        /*Check wether dynamic basal boudary conditions are activated -> TODO: Johannes :) */
+        parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
+        if(!isenthalpy) return;
+        //parameters->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
+        isdynamicbasalspc = true;
+        if(!isdynamicbasalspc) return;
+
+
+        /*Fetch indices of basal nodes for this finite element*/
+        BasalNodeIndices(&numindices,&indices,this->VelocityInterpolation());
+
+        /*Get parameters and inputs: */
+        Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input);
+
+        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
+        GaussPenta* gauss=new GaussPenta();
+        for(int i=0;i<numindices;i++){
+                gauss->GaussNode(this->element_type,indices[i]);
+
+                /*Check wether there is a temperate layer at the base or not -> TODO: Johannes:) */
+                istemperatelayer = false;
+
+                /*Add Dirichlet constraint to this node if there is a positive thickness of temperate ice*/
+                if(istemperatelayer){
+
+                        /*Calculate enthalpy at pressure melting point */
+                        pressure_input->GetInputValue(&pressure,gauss);
+                        h_pmp=matpar->PureIceEnthalpy(pressure);
+
+
+                        /*Apply Dirichlet condition (dof = 0 here, since there is only one degree of freedom per node)*/
+                        nodes[indices[i]]->ApplyConstraint(0,h_pmp);
+                }
+        }
+
+        /*Free ressources:*/
+        xDelete<int>(indices);
 }
 /*}}}*/

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

@@ -349 +349 @@
                 ElementVector* CreatePVectorThermalSheet(void);
                 void             GetSolutionFromInputsThermal(Vector<IssmDouble>* solutiong);
-                void           InputUpdateFromSolutionThermal( IssmDouble* solutiong);
-                void           InputUpdateFromSolutionEnthalpy( IssmDouble* solutiong);
+                void           InputUpdateFromSolutionThermal(IssmDouble* solutiong);
+                void           InputUpdateFromSolutionEnthalpy(IssmDouble* solutiong);
+                void           UpdateThermalBasalConstraints(void);
                 #endif
 

issm/trunk-jpl/src/c/classes/Elements/Tria.h ¶

@@ -246 +246 @@
                 #endif
 
+                #ifdef _HAVE_THERMAL_
+                void UpdateThermalBasalConstraints(void){_error_("not implemented yet");};
+                #endif
+
                 #ifdef _HAVE_HYDROLOGY_
                 ElementMatrix* CreateKMatrixHydrologyShreve(void);

issm/trunk-jpl/src/c/classes/FemModel.cpp ¶

@@ -1557 +1557 @@
 /*}}}*/
 #endif
-
 #ifdef _HAVE_HYDROLOGY_
 void FemModel::HydrologyEPLupdateDomainx(void){ /*{{{*/
@@ -1640 +1639 @@
 /*}}}*/
 #endif
+#ifdef _HAVE_THERMAL_
+void FemModel::UpdateThermalBasalConstraintsx(void){ /*{{{*/
+
+        for(int i=0;i<elements->Size();i++){
+                Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
+                element->UpdateThermalBasalConstraints();
+        }
+
+}
+/*}}}*/
+#endif

issm/trunk-jpl/src/c/classes/FemModel.h ¶

@@ -99 +99 @@
                 void UpdateConstraintsx(void);
                 int  UpdateVertexPositionsx(void);
-                void ParEplMask(void);          
+                #ifdef _HAVE_HYDROLOGY_
                 void HydrologyTransferx(void);
                 void HydrologyEPLupdateDomainx(void);
+                #endif
+                #ifdef _HAVE_THERMAL_
+                void UpdateThermalBasalConstraintsx(void);
+                #endif
 };