Changeset 14769

Timestamp:

Apr 26, 2013, 3:33:39 PM (13 years ago)

Author:

bdef

Message:

NEW: split HydrologyDC in 2 sub analyses, and added penalties for the sediment layer

Location:

issm/trunk-jpl/src/c

Files:

• EnumDefinitions/EnumDefinitions.h (modified) (2 diffs)
• Makefile.am (modified) (1 diff)
• classes/objects/Elements/Tria.cpp (modified) (15 diffs)
• classes/objects/Elements/Tria.h (modified) (1 diff)
• classes/objects/Loads/Pengrid.cpp (modified) (7 diffs)
• classes/objects/Loads/Pengrid.h (modified) (2 diffs)
• modules/EnumToStringx/EnumToStringx.cpp (modified) (2 diffs)
• modules/ModelProcessorx/CreateDataSets.cpp (modified) (1 diff)
• modules/ModelProcessorx/CreateParameters.cpp (modified) (1 diff)
• modules/ModelProcessorx/DistributeNumDofs.cpp (modified) (1 diff)
• modules/ModelProcessorx/HydrologyDCEfficient (added)
• modules/ModelProcessorx/HydrologyDCEfficient/CMakeLists.txt (added)
• modules/ModelProcessorx/HydrologyDCEfficient/CreateConstraintsHydrologyDCEfficient.cpp (added)
• modules/ModelProcessorx/HydrologyDCEfficient/CreateLoadsHydrologyDCEfficient.cpp (added)
• modules/ModelProcessorx/HydrologyDCEfficient/CreateNodesHydrologyDCEfficient.cpp (added)
• modules/ModelProcessorx/HydrologyDCEfficient/CreateParametersHydrologyDCEfficient.cpp (added)
• modules/ModelProcessorx/HydrologyDCEfficient/UpdateElementsHydrologyDCEfficient.cpp (added)
• modules/ModelProcessorx/HydrologyDCInefficient (added)
• modules/ModelProcessorx/HydrologyDCInefficient/CMakeLists.txt (added)
• modules/ModelProcessorx/HydrologyDCInefficient/CreateConstraintsHydrologyDCInefficient.cpp (added)
• modules/ModelProcessorx/HydrologyDCInefficient/CreateLoadsHydrologyDCInefficient.cpp (added)
• modules/ModelProcessorx/HydrologyDCInefficient/CreateNodesHydrologyDCInefficient.cpp (added)
• modules/ModelProcessorx/HydrologyDCInefficient/CreateParametersHydrologyDCInefficient.cpp (added)
• modules/ModelProcessorx/HydrologyDCInefficient/UpdateElementsHydrologyDCInefficient.cpp (added)
• modules/ModelProcessorx/HydrologyShreve (moved) (moved from issm/trunk-jpl/src/c/modules/ModelProcessorx/Hydrology )
• modules/ModelProcessorx/HydrologyShreve/CMakeLists.txt (modified) (1 diff)
• modules/ModelProcessorx/HydrologyShreve/CreateConstraintsHydrologyShreve.cpp (moved) (moved from issm/trunk-jpl/src/c/modules/ModelProcessorx/Hydrology/CreateConstraintsHydrology.cpp ) (3 diffs)
• modules/ModelProcessorx/HydrologyShreve/CreateLoadsHydrologyShreve.cpp (moved) (moved from issm/trunk-jpl/src/c/modules/ModelProcessorx/Hydrology/CreateLoadsHydrology.cpp ) (2 diffs)
• modules/ModelProcessorx/HydrologyShreve/CreateNodesHydrologyShreve.cpp (moved) (moved from issm/trunk-jpl/src/c/modules/ModelProcessorx/Hydrology/CreateNodesHydrology.cpp ) (4 diffs)
• modules/ModelProcessorx/HydrologyShreve/CreateParametersHydrologyShreve.cpp (moved) (moved from issm/trunk-jpl/src/c/modules/ModelProcessorx/Hydrology/CreateParametersHydrology.cpp ) (3 diffs)
• modules/ModelProcessorx/HydrologyShreve/UpdateElementsHydrologyShreve.cpp (moved) (moved from issm/trunk-jpl/src/c/modules/ModelProcessorx/Hydrology/UpdateElementsHydrology.cpp ) (4 diffs)
• modules/ModelProcessorx/ModelProcessorx.h (modified) (2 diffs)
• modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp (modified) (2 diffs)
• modules/StringToEnumx/StringToEnumx.cpp (modified) (8 diffs)
• solutions/AnalysisConfiguration.cpp (modified) (1 diff)
• solutions/hydrology_core.cpp (modified) (2 diffs)
• solvers/solver_hydro_nonlinear.cpp (added)
• solvers/solver_thermal_nonlinear.cpp (modified) (2 diffs)
• solvers/solvers.h (modified) (1 diff)

issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h

@@ -103 +103 @@
         HydrologydcEplTransmitivityEnum,
         HydrologydcIsefficientlayerEnum,
+        HydrologyLayerEnum,
+        HydrologySedimentEnum,
+        HydrologyEfficientEnum,
+        HydrologySedimentKmaxEnum,
         IndependentObjectEnum,
         InversionControlParametersEnum,
@@ -274 +278 @@
         FlaimAnalysisEnum,
         FlaimSolutionEnum,
-        HydrologyAnalysisEnum,
+        HydrologyShreveAnalysisEnum,
+        HydrologyDCInefficientAnalysisEnum,
+        HydrologyDCEfficientAnalysisEnum,
         HydrologySolutionEnum,
         MeltingAnalysisEnum,

issm/trunk-jpl/src/c/Makefile.am

@@ -454 +454 @@
 #}}}
 #Hydrology sources  {{{
-hydrology_sources  = ./modules/ModelProcessorx/Hydrology/UpdateElementsHydrology.cpp\
-                                              ./modules/ModelProcessorx/Hydrology/CreateNodesHydrology.cpp\
-                                              ./modules/ModelProcessorx/Hydrology/CreateConstraintsHydrology.cpp\
-                                              ./modules/ModelProcessorx/Hydrology/CreateLoadsHydrology.cpp \
-                                                        ./modules/ModelProcessorx/Hydrology/CreateParametersHydrology.cpp \
-                                                  ./solutions/hydrology_core.cpp
+hydrology_sources  = ./modules/ModelProcessorx/HydrologyShreve/UpdateElementsHydrologyShreve.cpp\
+                                              ./modules/ModelProcessorx/HydrologyShreve/CreateNodesHydrologyShreve.cpp\
+                                              ./modules/ModelProcessorx/HydrologyShreve/CreateConstraintsHydrologyShreve.cpp\
+                                              ./modules/ModelProcessorx/HydrologyShreve/CreateLoadsHydrologyShreve.cpp \
+                                                        ./modules/ModelProcessorx/HydrologyShreve/CreateParametersHydrologyShreve.cpp \
+                                                        ./modules/ModelProcessorx/HydrologyDCInefficient/UpdateElementsHydrologyDCInefficient.cpp\
+                                                        ./modules/ModelProcessorx/HydrologyDCInefficient/CreateNodesHydrologyDCInefficient.cpp\
+                                                        ./modules/ModelProcessorx/HydrologyDCInefficient/CreateConstraintsHydrologyDCInefficient.cpp\
+                                                        ./modules/ModelProcessorx/HydrologyDCInefficient/CreateLoadsHydrologyDCInefficient.cpp \
+                                                        ./modules/ModelProcessorx/HydrologyDCInefficient/CreateParametersHydrologyDCInefficient.cpp \
+                                                        ./modules/ModelProcessorx/HydrologyDCEfficient/UpdateElementsHydrologyDCEfficient.cpp\
+                                                        ./modules/ModelProcessorx/HydrologyDCEfficient/CreateNodesHydrologyDCEfficient.cpp\
+                                                        ./modules/ModelProcessorx/HydrologyDCEfficient/CreateConstraintsHydrologyDCEfficient.cpp\
+                                                        ./modules/ModelProcessorx/HydrologyDCEfficient/CreateLoadsHydrologyDCEfficient.cpp \
+                                                        ./modules/ModelProcessorx/HydrologyDCEfficient/CreateParametersHydrologyDCEfficient.cpp \
+                                                  ./solutions/hydrology_core.cpp\
+                                                  ./solvers/solver_hydro_nonlinear.cpp
 #}}}
 #Diagnostic sources  {{{

issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp

@@ -252 +252 @@
                         break;
                 #ifdef _HAVE_HYDROLOGY_
-                case HydrologyAnalysisEnum:
-                        Ke=CreateKMatrixHydrology();
+                case HydrologyShreveAnalysisEnum:
+                        Ke=CreateKMatrixHydrologyShreve();
+                        break;
+                case HydrologyDCInefficientAnalysisEnum:
+                        Ke=CreateKMatrixHydrologyDCInefficient();
+                        break;
+                case HydrologyDCEfficientAnalysisEnum:
+                        Ke=CreateKMatrixHydrologyDCEfficient();
                         break;
                 #endif
@@ -609 +615 @@
                         break;
                 #ifdef _HAVE_HYDROLOGY_
-                case HydrologyAnalysisEnum:
-                        pe=CreatePVectorHydrology();
+                case HydrologyShreveAnalysisEnum:
+                        pe=CreatePVectorHydrologyShreve();
+                        break;
+                case HydrologyDCInefficientAnalysisEnum:
+                        pe=CreatePVectorHydrologyDCInefficient();
+                        break;
+                case HydrologyDCEfficientAnalysisEnum:
+                        pe=CreatePVectorHydrologyDCEfficient();
                         break;
                 #endif
@@ -1383 +1395 @@
         #endif
         #ifdef _HAVE_HYDROLOGY_
-        case HydrologyAnalysisEnum:
-                GetSolutionFromInputsHydrology(solution);
+        case HydrologyShreveAnalysisEnum:
+                GetSolutionFromInputsHydrologyShreve(solution);
                 break;
         #endif
@@ -1721 +1733 @@
                 #ifdef _HAVE_DIAGNOSTIC_
                 case DiagnosticHorizAnalysisEnum:
-                        InputUpdateFromSolutionDiagnosticHoriz( solution);
+                        InputUpdateFromSolutionDiagnosticHoriz(solution);
                         break;
                 case DiagnosticHutterAnalysisEnum:
-                        InputUpdateFromSolutionDiagnosticHoriz( solution);
+                        InputUpdateFromSolutionDiagnosticHoriz(solution);
                         break;
                 #endif
                 #ifdef _HAVE_CONTROL_
                 case AdjointHorizAnalysisEnum:
-                        InputUpdateFromSolutionAdjointHoriz( solution);
+                        InputUpdateFromSolutionAdjointHoriz(solution);
                         break;
                 case AdjointBalancethicknessAnalysisEnum:
-                        InputUpdateFromSolutionAdjointBalancethickness( solution);
+                        InputUpdateFromSolutionAdjointBalancethickness(solution);
                         break;
                 #endif
                 #ifdef _HAVE_HYDROLOGY_ 
-                case HydrologyAnalysisEnum:
-                        InputUpdateFromSolutionHydrology(solution);
-                        break ;
+                case HydrologyShreveAnalysisEnum:
+                        InputUpdateFromSolutionHydrologyShreve(solution);
+                        break;
+                case HydrologyDCInefficientAnalysisEnum:
+                        InputUpdateFromSolutionHydrologyDCInefficient(solution);
+                        break;
+                case HydrologyDCEfficientAnalysisEnum:
+                        InputUpdateFromSolutionHydrologyDCEfficient(solution);
+                        break;
                 #endif
                 #ifdef _HAVE_BALANCED_
@@ -5802 +5820 @@
 }
 /*}}}*/
-/*FUNCTION Tria::CreateKMatrixHydrology{{{*/
-ElementMatrix* Tria::CreateKMatrixHydrology(void){
-
-        int  hydrology_model;
-        bool isefficientlayer;
-        parameters->FindParam(&hydrology_model,HydrologyEnum);
-        
-
-        switch(hydrology_model){
-                case HydrologyshreveEnum:
-                        return CreateKMatrixHydrologyShreve();
-                case HydrologydcEnum:
-                        parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
-                        if(!isefficientlayer)
-                                return CreateKMatrixHydrologyDCsediment();
-                        else if(isefficientlayer)
-                                return CreateKMatrixHydrologyDCepl();
-                default:
-                        _error_("Approximation " << EnumToStringx(hydrology_model) << " not supported yet");
-        }
-
-}
-/*}}}*/
-/*FUNCTION Tria::CreateKMatrixHydrologyShreve(void){{{*/
+/*FUNCTION Tria::CreateKMatrixHydrologyShreve{{{*/
 ElementMatrix* Tria::CreateKMatrixHydrologyShreve(void){
         
@@ -5881 +5876 @@
                 
                 TripleMultiply( &L[0],1,numdof,1,
-                                                                                &DL_scalar,1,1,0,
-                                                                                &L[0],1,numdof,0,
-                                                                                &Ke->values[0],1);
+                                        &DL_scalar,1,1,0,
+                                        &L[0],1,numdof,0,
+                                        &Ke->values[0],1);
                 
                 GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss);
@@ -5929 +5924 @@
 }
 /*}}}*/
-/*FUNCTION Tria::CreateKMatrixHydrologyDCsediment{{{*/
-ElementMatrix* Tria::CreateKMatrixHydrologyDCsediment(void){
+/*FUNCTION Tria::CreateKMatrixHydrologyDCInefficient{{{*/
+ElementMatrix* Tria::CreateKMatrixHydrologyDCInefficient(void){
 
         /*constants: */
@@ -5990 +5985 @@
 }
 /*}}}*/
-/*FUNCTION Tria::CreateKMatrixHydrologyDCepl{{{*/
-ElementMatrix* Tria::CreateKMatrixHydrologyDCepl(void){
+/*FUNCTION Tria::CreateKMatrixHydrologyDCEfficient{{{*/
+ElementMatrix* Tria::CreateKMatrixHydrologyDCEfficient(void){
 
         /*constants: */
@@ -6051 +6046 @@
 }
 /*}}}*/
-/*FUNCTION Tria::CreatePVectorHydrology{{{*/
-ElementVector* Tria::CreatePVectorHydrology(void){
-
-        int  hydrology_model;
-        bool isefficientlayer;
-        parameters->FindParam(&hydrology_model,HydrologyEnum);
-
-        switch(hydrology_model){
-                case HydrologyshreveEnum:
-                        return CreatePVectorHydrologyShreve();
-                case HydrologydcEnum:
-                        parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
-                        if(!isefficientlayer)
-                                return CreatePVectorHydrologyDCsediment();
-                        else if(isefficientlayer)
-                                return CreatePVectorHydrologyDCepl();
-                default:
-                        _error_("Approximation " << EnumToStringx(hydrology_model) << " not supported yet");
-        }
-
-}
-/*}}}*/
 /*FUNCTION Tria::CreatePVectorHydrologyShreve {{{*/
 ElementVector* Tria::CreatePVectorHydrologyShreve(void){
@@ -6122 +6095 @@
 }
 /*}}}*/
-/*FUNCTION Tria::CreatePVectorHydrologyDCsediment {{{*/
-ElementVector* Tria::CreatePVectorHydrologyDCsediment(void){
+/*FUNCTION Tria::CreatePVectorHydrologyDCInefficient {{{*/
+ElementVector* Tria::CreatePVectorHydrologyDCInefficient(void){
 
         /*Constants*/
@@ -6179 +6152 @@
 }
 /*}}}*/
-/*FUNCTION Tria::CreatePVectorHydrologyDCepl {{{*/
-ElementVector* Tria::CreatePVectorHydrologyDCepl(void){
+/*FUNCTION Tria::CreatePVectorHydrologyDCEfficient {{{*/
+ElementVector* Tria::CreatePVectorHydrologyDCEfficient(void){
 
         /*Constants*/
@@ -6236 +6209 @@
 }
 /*}}}*/
-/*FUNCTION Tria::GetSolutionFromInputsHydrology{{{*/
-void  Tria::GetSolutionFromInputsHydrology(Vector<IssmDouble>* solution){
+/*FUNCTION Tria::GetSolutionFromInputsHydrologyShreve{{{*/
+void  Tria::GetSolutionFromInputsHydrologyShreve(Vector<IssmDouble>* solution){
 
         const int    numdof=NDOF1*NUMVERTICES;
 
-        int i;
-        int*         doflist=NULL;
-        IssmDouble       watercolumn;
-        IssmDouble       values[numdof];
-        GaussTria*   gauss=NULL;
+        int         i;
+        int        *doflist = NULL;
+        IssmDouble  watercolumn;
+        IssmDouble  values[numdof];
+        GaussTria  *gauss   = NULL;
 
         /*Get dof list: */
@@ -6270 +6243 @@
         delete gauss;
         xDelete<int>(doflist);
-}
-/*}}}*/
-/*FUNCTION Tria::InputUpdateFromSolutionHydrology{{{*/
-void  Tria::InputUpdateFromSolutionHydrology(IssmDouble* solution){
-
-        int hydrology_model;
-  parameters->FindParam(&hydrology_model,HydrologyEnum);
-        
-  switch(hydrology_model){
-        case HydrologyshreveEnum:
-                return InputUpdateFromSolutionHydrologyShreve(solution);
-        case HydrologydcEnum:
-                return InputUpdateFromSolutionHydrologyDC(solution);
-        default:
-                _error_("Approximation " << EnumToStringx(hydrology_model) << " not supported yet");
-  }
 }
 /*}}}*/
@@ -6313 +6270 @@
 }
 /*}}}*/
-/*FUNCTION Tria::InputUpdateFromSolutionHydrologyDC{{{*/
-void  Tria::InputUpdateFromSolutionHydrologyDC(IssmDouble* solution){
+/*FUNCTION Tria::InputUpdateFromSolutionHydrologyDCInefficient{{{*/
+void  Tria::InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution){
 
         /*Intermediaries*/
@@ -6332 +6289 @@
         /*Add input to the element: */
         this->inputs->AddInput(new TriaP1Input(SedimentHeadEnum,values));
+
+        /*Free ressources:*/
+        xDelete<int>(doflist);
+}
+/*}}}*/
+/*FUNCTION Tria::InputUpdateFromSolutionHydrologyDCEfficient{{{*/
+void  Tria::InputUpdateFromSolutionHydrologyDCEfficient(IssmDouble* solution){
+
+        /*Intermediaries*/
+        const int   numdof         = NDOF1 *NUMVERTICES;
+        int        *doflist        = NULL;
+        IssmDouble  values[numdof];
+
+        /*Get dof list: */
+        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+
+        /*Use the dof list to index into the solution vector: */
+        for(int i=0;i<numdof;i++){
+                values[i]=solution[doflist[i]];
+                if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
+        }
+
+        /*Add input to the element: */
+        this->inputs->AddInput(new TriaP1Input(EplHeadEnum,values));
 
         /*Free ressources:*/

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

@@ -239 +239 @@
 
                 #ifdef _HAVE_HYDROLOGY_
-                ElementMatrix* CreateKMatrixHydrology(void);
                 ElementMatrix* CreateKMatrixHydrologyShreve(void);
-                ElementMatrix* CreateKMatrixHydrologyDCsediment(void);
-                ElementMatrix* CreateKMatrixHydrologyDCepl(void);
-                ElementVector* CreatePVectorHydrology(void);
+                ElementMatrix* CreateKMatrixHydrologyDCInefficient(void);
+                ElementMatrix* CreateKMatrixHydrologyDCEfficient(void);
                 ElementVector* CreatePVectorHydrologyShreve(void);
-                ElementVector* CreatePVectorHydrologyDCsediment(void);
-                ElementVector* CreatePVectorHydrologyDCepl(void);
+                ElementVector* CreatePVectorHydrologyDCInefficient(void);
+                ElementVector* CreatePVectorHydrologyDCEfficient(void);
+                void      GetSolutionFromInputsHydrologyShreve(Vector<IssmDouble>* solution);
                 void    CreateHydrologyWaterVelocityInput(void);
-                void      GetSolutionFromInputsHydrology(Vector<IssmDouble>* solution);
                 void      InputUpdateFromSolutionHydrology(IssmDouble* solution);
                 void      InputUpdateFromSolutionHydrologyShreve(IssmDouble* solution);
-                void      InputUpdateFromSolutionHydrologyDC(IssmDouble* solution);
+                void    InputUpdateFromSolutionHydrologyDC(IssmDouble* solution);
+                void      InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution);
+                void      InputUpdateFromSolutionHydrologyDCEfficient(IssmDouble* solution);
                 #endif
                 #ifdef _HAVE_BALANCED_

issm/trunk-jpl/src/c/classes/objects/Loads/Pengrid.cpp

@@ -246 +246 @@
                         break;
                 #endif
+                #ifdef _HAVE_HYDROLOGY_
+                case HydrologyDCInefficientAnalysisEnum:
+                        Ke=PenaltyCreateKMatrixHydrologyDCInefficient(kmax);
+                        break;
+                #endif
                 default:
                         _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
@@ -274 +279 @@
                         break;
                 case DiagnosticHorizAnalysisEnum: case AdjointHorizAnalysisEnum:
+                        break;
+                #endif
+                #ifdef _HAVE_HYDROLOGY_
+                case HydrologyDCInefficientAnalysisEnum:
+                        pe=PenaltyCreatePVectorHydrologyDCInefficient(kmax);
                         break;
                 #endif
@@ -429 +439 @@
                 return;
         }
+        #ifdef _HAVE_THERMAL_
         else if (analysis_type==ThermalAnalysisEnum){
                 ConstraintActivateThermal(punstable);
@@ -436 +447 @@
                 return;
         }
+        #endif
+        #ifdef _HAVE_HYDROLOGY_
+        else if (analysis_type==HydrologyDCInefficientAnalysisEnum){
+                ConstraintActivateHydrologyDCInefficient(punstable);
+                return;
+        }
+        #endif
         else{
                 _error_("analysis: " << EnumToStringx(analysis_type) << " not supported yet");
@@ -442 +460 @@
 }
 /*}}}*/
+#ifdef _HAVE_DIAGNOSTIC_
+/*FUNCTION Pengrid::PenaltyCreateKMatrixDiagnosticStokes {{{*/
+ElementMatrix* Pengrid::PenaltyCreateKMatrixDiagnosticStokes(IssmDouble kmax){
+
+        const int numdof = NUMVERTICES *NDOF4;
+        IssmDouble    slope[2];
+        IssmDouble    penalty_offset;
+        int       approximation;
+
+        Penta* penta=(Penta*)element;
+
+        /*Initialize Element vector and return if necessary*/
+        penta->inputs->GetInputValue(&approximation,ApproximationEnum);
+        if(approximation!=StokesApproximationEnum &&  approximation!=PattynStokesApproximationEnum) return NULL;
+        ElementMatrix* Ke=new ElementMatrix(&node,1,this->parameters,StokesApproximationEnum);
+
+        /*Retrieve all inputs and parameters*/
+        parameters->FindParam(&penalty_offset,DiagnosticPenaltyFactorEnum);
+        penta->GetInputValue(&slope[0],node,BedSlopeXEnum);
+        penta->GetInputValue(&slope[1],node,BedSlopeYEnum);
+
+        /*Create elementary matrix: add penalty to constrain wb (wb=ub*db/dx+vb*db/dy)*/
+        Ke->values[2*NDOF4+0]=-slope[0]*kmax*pow((IssmDouble)10.0,penalty_offset);
+        Ke->values[2*NDOF4+1]=-slope[1]*kmax*pow((IssmDouble)10.0,penalty_offset);
+        Ke->values[2*NDOF4+2]= kmax*pow((IssmDouble)10,penalty_offset);
+
+        /*Transform Coordinate System*/
+        TransformStiffnessMatrixCoord(Ke,&node,NUMVERTICES,XYZPEnum);
+
+        /*Clean up and return*/
+        return Ke;
+}
+/*}}}*/
+#endif
+#ifdef _HAVE_THERMAL_
 /*FUNCTION Pengrid::ConstraintActivateThermal {{{*/
 void  Pengrid::ConstraintActivateThermal(int* punstable){
@@ -510 +563 @@
 }
 /*}}}*/
-#ifdef _HAVE_DIAGNOSTIC_
-/*FUNCTION Pengrid::PenaltyCreateKMatrixDiagnosticStokes {{{*/
-ElementMatrix* Pengrid::PenaltyCreateKMatrixDiagnosticStokes(IssmDouble kmax){
-
-        const int numdof = NUMVERTICES *NDOF4;
-        IssmDouble    slope[2];
-        IssmDouble    penalty_offset;
-        int       approximation;
-
-        Penta* penta=(Penta*)element;
-
-        /*Initialize Element vector and return if necessary*/
-        penta->inputs->GetInputValue(&approximation,ApproximationEnum);
-        if(approximation!=StokesApproximationEnum &&  approximation!=PattynStokesApproximationEnum) return NULL;
-        ElementMatrix* Ke=new ElementMatrix(&node,1,this->parameters,StokesApproximationEnum);
-
-        /*Retrieve all inputs and parameters*/
-        parameters->FindParam(&penalty_offset,DiagnosticPenaltyFactorEnum);
-        penta->GetInputValue(&slope[0],node,BedSlopeXEnum);
-        penta->GetInputValue(&slope[1],node,BedSlopeYEnum);
-
-        /*Create elementary matrix: add penalty to constrain wb (wb=ub*db/dx+vb*db/dy)*/
-        Ke->values[2*NDOF4+0]=-slope[0]*kmax*pow((IssmDouble)10.0,penalty_offset);
-        Ke->values[2*NDOF4+1]=-slope[1]*kmax*pow((IssmDouble)10.0,penalty_offset);
-        Ke->values[2*NDOF4+2]= kmax*pow((IssmDouble)10,penalty_offset);
-
-        /*Transform Coordinate System*/
-        TransformStiffnessMatrixCoord(Ke,&node,NUMVERTICES,XYZPEnum);
-
-        /*Clean up and return*/
-        return Ke;
-}
-/*}}}*/
-#endif
-#ifdef _HAVE_THERMAL_
 /*FUNCTION Pengrid::PenaltyCreateKMatrixMelting {{{*/
 ElementMatrix* Pengrid::PenaltyCreateKMatrixMelting(IssmDouble kmax){
@@ -665 +683 @@
 /*}}}*/
 #endif
+#ifdef _HAVE_HYDROLOGY_
+/*FUNCTION Pengrid::ConstraintActivateHydrologyDCInefficient{{{*/
+void  Pengrid::ConstraintActivateHydrologyDCInefficient(int* punstable){
+
+        //   The penalty is stable if it doesn't change during two consecutive iterations.   
+        int        found           = 0;
+        const int  numnodes        = 1;
+        IssmDouble pressure;
+        IssmDouble h;
+        IssmDouble h_max;
+        int        new_active;
+        int        unstable        = 0;
+        int        reset_penalties = 0;
+        int        penalty_lock;
+
+        /*check that pengrid is not a clone (penalty to be added only once)*/
+        if(node->IsClone()){
+                unstable=0;
+                *punstable=unstable;
+                return;
+        }
+
+        /*Get sediment water head h*/
+        element->GetInputValue(&h,node,SedimentHeadEnum);
+        h_max=10000;
+
+        if (h>h_max)
+         new_active=1;
+        else
+         new_active=0;
+
+        if(this->active==new_active)
+         unstable=0;
+        else
+         unstable=1;
+
+        /*Set penalty flag*/
+        this->active=new_active;
+
+        /*Assign output pointers:*/
+        *punstable=unstable;
+}
+/*}}}*/
+/*FUNCTION Pengrid::PenaltyCreateKMatrixHydrologyDCInefficient {{{*/
+ElementMatrix* Pengrid::PenaltyCreateKMatrixHydrologyDCInefficient(IssmDouble kmax){
+
+        const int numdof=NUMVERTICES*NDOF1;
+        IssmDouble    penalty_factor = 3.;
+
+        /*Initialize Element matrix and return if necessary*/
+        if(!this->active) return NULL;
+        ElementMatrix* Ke=new ElementMatrix(&node,NUMVERTICES,this->parameters);
+
+        Ke->values[0]=kmax*pow(10.,penalty_factor);
+
+        /*Clean up and return*/
+        return Ke;
+}
+/*}}}*/
+/*FUNCTION Pengrid::PenaltyCreatePVectorHydrologyDCInefficient {{{*/
+ElementVector* Pengrid::PenaltyCreatePVectorHydrologyDCInefficient(IssmDouble kmax){
+
+        const int  numdof=NUMVERTICES*NDOF1;
+        IssmDouble h_max;
+        IssmDouble penalty_factor=3.;
+
+        /*Initialize Element matrix and return if necessary*/
+        if(!this->active) return NULL;
+        ElementVector* pe=new ElementVector(&node,1,this->parameters);
+
+        /*Get h_max*/
+        h_max = 10000.;
+
+        pe->values[0]=kmax*pow(10.,penalty_factor)*h_max;
+
+        /*Clean up and return*/
+        return pe;
+}
+/*}}}*/
+#endif
 /*FUNCTION Pengrid::ResetConstraint {{{*/
 void  Pengrid::ResetConstraint(void){

issm/trunk-jpl/src/c/classes/objects/Loads/Pengrid.h

@@ -7 +7 @@
 /*Headers:*/
 /*{{{*/
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
 #include "./Load.h"
 class Hook;
@@ -90 +95 @@
                 ElementVector* PenaltyCreatePVectorThermal(IssmDouble kmax);
                 ElementVector* PenaltyCreatePVectorMelting(IssmDouble kmax);
+                void           ConstraintActivateThermal(int* punstable);
+                #endif
+                #ifdef _HAVE_HYDROLOGY_
+                ElementMatrix* PenaltyCreateKMatrixHydrologyDCInefficient(IssmDouble kmax);
+                ElementVector* PenaltyCreatePVectorHydrologyDCInefficient(IssmDouble kmax);
+                void           ConstraintActivateHydrologyDCInefficient(int* punstable);
                 #endif
                 void  ConstraintActivate(int* punstable);
-                void  ConstraintActivateThermal(int* punstable);
                 void  UpdateInputs(IssmDouble* solution);
                 void  ResetConstraint(void);

issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp

@@ -108 +108 @@
                 case HydrologydcEplTransmitivityEnum : return "HydrologydcEplTransmitivity";
                 case HydrologydcIsefficientlayerEnum : return "HydrologydcIsefficientlayer";
+                case HydrologyLayerEnum : return "HydrologyLayer";
+                case HydrologySedimentEnum : return "HydrologySediment";
+                case HydrologyEfficientEnum : return "HydrologyEfficient";
+                case HydrologySedimentKmaxEnum : return "HydrologySedimentKmax";
                 case IndependentObjectEnum : return "IndependentObject";
                 case InversionControlParametersEnum : return "InversionControlParameters";
@@ -277 +281 @@
                 case FlaimAnalysisEnum : return "FlaimAnalysis";
                 case FlaimSolutionEnum : return "FlaimSolution";
-                case HydrologyAnalysisEnum : return "HydrologyAnalysis";
+                case HydrologyShreveAnalysisEnum : return "HydrologyShreveAnalysis";
+                case HydrologyDCInefficientAnalysisEnum : return "HydrologyDCInefficientAnalysis";
+                case HydrologyDCEfficientAnalysisEnum : return "HydrologyDCEfficientAnalysis";
                 case HydrologySolutionEnum : return "HydrologySolution";
                 case MeltingAnalysisEnum : return "MeltingAnalysis";

issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateDataSets.cpp

@@ -56 +56 @@
 
                 #ifdef _HAVE_HYDROLOGY_
-                case HydrologyAnalysisEnum:
-                        CreateNodesHydrology(pnodes, iomodel);
-                        CreateConstraintsHydrology(pconstraints,iomodel);
-                        CreateLoadsHydrology(ploads,iomodel);
-                        UpdateElementsHydrology(elements,iomodel,analysis_counter,analysis_type);
+                case HydrologyShreveAnalysisEnum:
+                        CreateNodesHydrologyShreve(pnodes, iomodel);
+                        CreateConstraintsHydrologyShreve(pconstraints,iomodel);
+                        CreateLoadsHydrologyShreve(ploads,iomodel);
+                        UpdateElementsHydrologyShreve(elements,iomodel,analysis_counter,analysis_type);
+                        break;
+                case HydrologyDCInefficientAnalysisEnum:
+                        CreateNodesHydrologyDCInefficient(pnodes, iomodel);
+                        CreateConstraintsHydrologyDCInefficient(pconstraints,iomodel);
+                        CreateLoadsHydrologyDCInefficient(ploads,iomodel);
+                        UpdateElementsHydrologyDCInefficient(elements,iomodel,analysis_counter,analysis_type);
+                        break;
+                case HydrologyDCEfficientAnalysisEnum:
+                        CreateNodesHydrologyDCEfficient(pnodes, iomodel);
+                        CreateConstraintsHydrologyDCEfficient(pconstraints,iomodel);
+                        CreateLoadsHydrologyDCEfficient(ploads,iomodel);
+                        UpdateElementsHydrologyDCEfficient(elements,iomodel,analysis_counter,analysis_type);
                         break;
                 #endif

issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp

@@ -226 +226 @@
         /*Solution specific parameters (FIXME: extend to other params)*/
         #ifdef _HAVE_HYDROLOGY_
-        CreateParametersHydrology(&parameters,iomodel,solution_type,analysis_type);
+        CreateParametersHydrologyShreve(&parameters,iomodel,solution_type,analysis_type);
+        CreateParametersHydrologyDCInefficient(&parameters,iomodel,solution_type,analysis_type);
+        CreateParametersHydrologyDCEfficient(&parameters,iomodel,solution_type,analysis_type);
         #endif
 

issm/trunk-jpl/src/c/modules/ModelProcessorx/DistributeNumDofs.cpp

@@ -80 +80 @@
                 numdofs=1;
         }
-        else if (analysis_type==HydrologyAnalysisEnum){
+        else if (analysis_type==HydrologyDCInefficientAnalysisEnum){
+                numdofs=1;
+        }
+        else if (analysis_type==HydrologyDCEfficientAnalysisEnum){
+                numdofs=1;
+        }
+        else if (analysis_type==HydrologyShreveAnalysisEnum){
                 numdofs=1;
         }

issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/CMakeLists.txt

@@ -2 +2 @@
 # }}}
 # HYDROLOGY_SOURCES {{{
-set(HYDROLOGY_SOURCES $ENV{ISSM_DIR}/src/c/modules/ModelProcessorx/Hydrology/CreateConstraintsHydrology.cpp
-                            $ENV{ISSM_DIR}/src/c/modules/ModelProcessorx/Hydrology/CreateLoadsHydrology.cpp
-                            $ENV{ISSM_DIR}/src/c/modules/ModelProcessorx/Hydrology/CreateNodesHydrology.cpp
-                         $ENV{ISSM_DIR}/src/c/modules/ModelProcessorx/Hydrology/UpdateElementsHydrology.cpp PARENT_SCOPE)
+set(HYDROLOGY_SOURCES $ENV{ISSM_DIR}/src/c/modules/ModelProcessorx/HydrologyShreve/CreateConstraintsHydrologyShreve.cpp
+                            $ENV{ISSM_DIR}/src/c/modules/ModelProcessorx/HydrologyShreve/CreateLoadsHydrologyShreve.cpp
+                            $ENV{ISSM_DIR}/src/c/modules/ModelProcessorx/HydrologyShreve/CreateNodesHydrologyShreve.cpp
+                         $ENV{ISSM_DIR}/src/c/modules/ModelProcessorx/HydrologyShreve/UpdateElementsHydrologyShreve.cpp PARENT_SCOPE)
 # }}}

issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/CreateConstraintsHydrologyShreve.cpp

@@ -1 +1 @@
 /*
- * CreateConstraintsHydrology.c:
+ * CreateConstraintsHydrologyShreve.c:
  */
 
@@ -12 +12 @@
 #include "../ModelProcessorx.h"
 
-void    CreateConstraintsHydrology(Constraints** pconstraints, IoModel* iomodel){
+void    CreateConstraintsHydrologyShreve(Constraints** pconstraints, IoModel* iomodel){
 
         /*Recover pointer: */
@@ -25 +25 @@
         if(!constraints) constraints = new Constraints();
 
-        if(hydrology_model==HydrologyshreveEnum){
-                IoModelToConstraintsx(constraints,iomodel,HydrologyshreveSpcwatercolumnEnum,HydrologyAnalysisEnum);
+        if(hydrology_model!=HydrologyshreveEnum){
+                *pconstraints=constraints;
+                return;
         }
-        else if(hydrology_model==HydrologydcEnum){
-                iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum);
-                IoModelToConstraintsx(constraints,iomodel,HydrologydcSpcsedimentHeadEnum,HydrologyAnalysisEnum);
-                
-                if(isefficientlayer){
-                        IoModelToConstraintsx(constraints,iomodel,HydrologydcSpceplHeadEnum,HydrologyAnalysisEnum);
-                }
-        }
-        else{
-                _error_("Hydrology model "<<EnumToStringx(hydrology_model)<<" not supported yet");
-        }
+
+        IoModelToConstraintsx(constraints,iomodel,HydrologyshreveSpcwatercolumnEnum,HydrologyShreveAnalysisEnum);
 
         /*Assign output pointer: */

issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/CreateLoadsHydrologyShreve.cpp

@@ -1 @@
-/*! \file CreateLoadsHydrology.c:
+/*! \file CreateLoadsHydrologyShreve.c:
  */
 
@@ -11 +11 @@
 #include "../ModelProcessorx.h"
 
-void    CreateLoadsHydrology(Loads** ploads, IoModel* iomodel){
+void    CreateLoadsHydrologyShreve(Loads** ploads, IoModel* iomodel){
+
+        /*Intermediary*/
+        int      numberofvertices;
+        Pengrid *pengrid = NULL;
 
         /*Recover pointer: */

issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/CreateNodesHydrologyShreve.cpp

@@ -1 +1 @@
 /*
- * CreateNodesHydrology.c:
+ * CreateNodesHydrologyShreve.c:
  */
 
@@ -13 +13 @@
 #include "../ModelProcessorx.h"
 
-void    CreateNodesHydrology(Nodes** pnodes, IoModel* iomodel){
+void    CreateNodesHydrologyShreve(Nodes** pnodes, IoModel* iomodel){
 
         /*Intermediary*/
-        int i;
+        int  i;
+        int  hydrology_model;
         bool continuous_galerkin=true;
-        int    numberofvertices;
+        int  numberofvertices;
 
         /*Fetch parameters: */
         iomodel->Constant(&numberofvertices,MeshNumberofverticesEnum);
+        iomodel->Constant(&hydrology_model,HydrologyEnum);
 
         /*Recover pointer: */
@@ -28 +30 @@
         /*Create nodes if they do not exist yet*/
         if(!nodes) nodes = new Nodes();
+
+        /*Now, do we really want Shreve?*/
+        if(hydrology_model!=HydrologyshreveEnum){
+                *pnodes=nodes;
+                return;
+        }
 
         /*Continuous Galerkin partition of nodes: */
@@ -38 +46 @@
                 if(iomodel->my_vertices[i]){
                         /*Add node to nodes dataset: */
-                        nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i+1,i,iomodel,HydrologyAnalysisEnum));
-
+                        nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i+1,i,iomodel,HydrologyShreveAnalysisEnum));
                 }
         }

issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/CreateParametersHydrologyShreve.cpp

@@ -1 @@
-/*!\file: CreateParametersHydrology.cpp
+/*!\file: CreateParametersHydrologyShreve.cpp
  * \brief driver for creating parameters dataset, for control analysis.
  */ 
@@ -12 +12 @@
 #include "../ModelProcessorx.h"
 
-void CreateParametersHydrology(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type){
+void CreateParametersHydrologyShreve(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type){
 
         Parameters *parameters = NULL;
@@ -24 +24 @@
         iomodel->Constant(&hydrology_model,HydrologyEnum);
 
+        /*Now, do we really want Shreve?*/
+        if(hydrology_model!=HydrologyshreveEnum){
+                *pparameters=parameters;
+                return;
+        }
+
         parameters->AddObject(new IntParam(HydrologyEnum,hydrology_model));
-        if(hydrology_model==HydrologyshreveEnum){
-                parameters->AddObject(iomodel->CopyConstantObject(HydrologyshreveStabilizationEnum));
-        }
-        else if(hydrology_model==HydrologydcEnum){
-                iomodel->FetchData(&isefficientlayer,HydrologydcIsefficientlayerEnum);
-                        parameters->AddObject(new BoolParam(HydrologydcIsefficientlayerEnum,isefficientlayer));
-        }
-        else{
-                _error_("Hydrology model "<<EnumToStringx(hydrology_model)<<" not supported yet");
-        }
+        parameters->AddObject(iomodel->CopyConstantObject(HydrologyshreveStabilizationEnum));
 
         /*Assign output pointer: */

issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/UpdateElementsHydrologyShreve.cpp

@@ -1 +1 @@
 /*
- * UpdateElementsHydrology:
+ * UpdateElementsHydrologyShreve:
  */
 
@@ -14 +14 @@
 #include "../ModelProcessorx.h"
 
-void    UpdateElementsHydrology(Elements* elements, IoModel* iomodel,int analysis_counter,int analysis_type){
+void    UpdateElementsHydrologyShreve(Elements* elements, IoModel* iomodel,int analysis_counter,int analysis_type){
 
         int    hydrology_model;
@@ -23 +23 @@
         iomodel->Constant(&numberofelements,MeshNumberofelementsEnum);
         iomodel->FetchData(1,MeshElementsEnum);
+
+        /*Now, do we really want Shreve?*/
+        if(hydrology_model!=HydrologyshreveEnum) return;
 
         /*Update elements: */
@@ -42 +45 @@
         iomodel->FetchDataToInput(elements,MaskElementonwaterEnum);
         iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
-        if(hydrology_model==HydrologyshreveEnum){
-                iomodel->FetchDataToInput(elements,WatercolumnEnum);
-                elements->InputDuplicate(WatercolumnEnum,WaterColumnOldEnum);
-        }
-        else if(hydrology_model==HydrologydcEnum){
-                iomodel->FetchDataToInput(elements,SedimentHeadEnum);
-        }
-        else{
-                _error_("Hydrology model "<<EnumToStringx(hydrology_model)<<" not supported yet");
-        }
+        iomodel->FetchDataToInput(elements,WatercolumnEnum);
+
+        elements->InputDuplicate(WatercolumnEnum,WaterColumnOldEnum);
 
         /*Free data: */

issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h

@@ -24 +24 @@
 void  CreateParametersControl(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type);
 void  CreateParametersDakota(Parameters** pparameters,IoModel* iomodel,char* rootpath,int solution_type,int analysis_type);
-void  CreateParametersHydrology(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type);
+void  CreateParametersHydrologyShreve(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type);
+void  CreateParametersHydrologyDCInefficient(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type);
+void  CreateParametersHydrologyDCEfficient(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type);
 void  UpdateElementsAndMaterialsControl(Elements* elements,Materials* materials, IoModel* iomodel);
 
@@ -79 +81 @@
 void    UpdateElementsEnthalpy(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
 
-/*hydrology:*/
-void    CreateNodesHydrology(Nodes** pnodes,IoModel* iomodel);
-void    CreateConstraintsHydrology(Constraints** pconstraints,IoModel* iomodel);
-void  CreateLoadsHydrology(Loads** ploads, IoModel* iomodel);
-void    UpdateElementsHydrology(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
+/*hydrology Shreve:*/
+void    CreateNodesHydrologyShreve(Nodes** pnodes,IoModel* iomodel);
+void    CreateConstraintsHydrologyShreve(Constraints** pconstraints,IoModel* iomodel);
+void  CreateLoadsHydrologyShreve(Loads** ploads, IoModel* iomodel);
+void    UpdateElementsHydrologyShreve(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
+
+/*hydrology DC:*/
+void    CreateNodesHydrologyDCInefficient(Nodes** pnodes,IoModel* iomodel);
+void    CreateConstraintsHydrologyDCInefficient(Constraints** pconstraints,IoModel* iomodel);
+void  CreateLoadsHydrologyDCInefficient(Loads** ploads, IoModel* iomodel);
+void    UpdateElementsHydrologyDCInefficient(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
+void    CreateNodesHydrologyDCEfficient(Nodes** pnodes,IoModel* iomodel);
+void    CreateConstraintsHydrologyDCEfficient(Constraints** pconstraints,IoModel* iomodel);
+void  CreateLoadsHydrologyDCEfficient(Loads** ploads, IoModel* iomodel);
+void    UpdateElementsHydrologyDCEfficient(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
 
 /*melting:*/

issm/trunk-jpl/src/c/modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp

@@ -14 +14 @@
 
         /*Intermediary*/
-        int i;
-        int    dim;
-        int    numberofvertices;
-        Pengrid*    pengrid  = NULL;
+        int      dim;
+        int      numberofvertices;
+        Pengrid *pengrid          = NULL;
 
         /*Recover pointer: */
@@ -36 +35 @@
         CreateSingleNodeToElementConnectivity(iomodel);
 
-        for (i=0;i<numberofvertices;i++){
+        for(int i=0;i<numberofvertices;i++){
 
                 /*keep only this partition's nodes:*/

issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp

@@ -109 +109 @@
               else if (strcmp(name,"HydrologydcEplTransmitivity")==0) return HydrologydcEplTransmitivityEnum;
               else if (strcmp(name,"HydrologydcIsefficientlayer")==0) return HydrologydcIsefficientlayerEnum;
+              else if (strcmp(name,"HydrologyLayer")==0) return HydrologyLayerEnum;
+              else if (strcmp(name,"HydrologySediment")==0) return HydrologySedimentEnum;
+              else if (strcmp(name,"HydrologyEfficient")==0) return HydrologyEfficientEnum;
+              else if (strcmp(name,"HydrologySedimentKmax")==0) return HydrologySedimentKmaxEnum;
               else if (strcmp(name,"IndependentObject")==0) return IndependentObjectEnum;
               else if (strcmp(name,"InversionControlParameters")==0) return InversionControlParametersEnum;
@@ -134 +138 @@
               else if (strcmp(name,"MaskElementongroundedice")==0) return MaskElementongroundediceEnum;
               else if (strcmp(name,"MaskElementonwater")==0) return MaskElementonwaterEnum;
-              else if (strcmp(name,"MaskVertexonfloatingice")==0) return MaskVertexonfloatingiceEnum;
+         else stage=2;
+   }
+   if(stage==2){
+              if (strcmp(name,"MaskVertexonfloatingice")==0) return MaskVertexonfloatingiceEnum;
               else if (strcmp(name,"MaskVertexongroundedice")==0) return MaskVertexongroundediceEnum;
               else if (strcmp(name,"MaskVertexonwater")==0) return MaskVertexonwaterEnum;
               else if (strcmp(name,"MaterialsBeta")==0) return MaterialsBetaEnum;
-         else stage=2;
-   }
-   if(stage==2){
-              if (strcmp(name,"MaterialsHeatcapacity")==0) return MaterialsHeatcapacityEnum;
+              else if (strcmp(name,"MaterialsHeatcapacity")==0) return MaterialsHeatcapacityEnum;
               else if (strcmp(name,"MaterialsLatentheat")==0) return MaterialsLatentheatEnum;
               else if (strcmp(name,"MaterialsMeltingpoint")==0) return MaterialsMeltingpointEnum;
@@ -257 +261 @@
               else if (strcmp(name,"TransientRequestedOutputs")==0) return TransientRequestedOutputsEnum;
               else if (strcmp(name,"SolutionType")==0) return SolutionTypeEnum;
-              else if (strcmp(name,"AnalysisType")==0) return AnalysisTypeEnum;
+         else stage=3;
+   }
+   if(stage==3){
+              if (strcmp(name,"AnalysisType")==0) return AnalysisTypeEnum;
               else if (strcmp(name,"ConfigurationType")==0) return ConfigurationTypeEnum;
               else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
               else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum;
-         else stage=3;
-   }
-   if(stage==3){
-              if (strcmp(name,"AdjointSolution")==0) return AdjointSolutionEnum;
+              else if (strcmp(name,"AdjointSolution")==0) return AdjointSolutionEnum;
               else if (strcmp(name,"AnalysisCounter")==0) return AnalysisCounterEnum;
               else if (strcmp(name,"NoneAnalysis")==0) return NoneAnalysisEnum;
@@ -284 +288 @@
               else if (strcmp(name,"FlaimAnalysis")==0) return FlaimAnalysisEnum;
               else if (strcmp(name,"FlaimSolution")==0) return FlaimSolutionEnum;
-              else if (strcmp(name,"HydrologyAnalysis")==0) return HydrologyAnalysisEnum;
+              else if (strcmp(name,"HydrologyShreveAnalysis")==0) return HydrologyShreveAnalysisEnum;
+              else if (strcmp(name,"HydrologyDCInefficientAnalysis")==0) return HydrologyDCInefficientAnalysisEnum;
+              else if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum;
               else if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum;
               else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum;
@@ -378 +384 @@
               else if (strcmp(name,"Water")==0) return WaterEnum;
               else if (strcmp(name,"Closed")==0) return ClosedEnum;
-              else if (strcmp(name,"Free")==0) return FreeEnum;
+         else stage=4;
+   }
+   if(stage==4){
+              if (strcmp(name,"Free")==0) return FreeEnum;
               else if (strcmp(name,"Open")==0) return OpenEnum;
               else if (strcmp(name,"Adjointp")==0) return AdjointpEnum;
@@ -384 +393 @@
               else if (strcmp(name,"Adjointy")==0) return AdjointyEnum;
               else if (strcmp(name,"Adjointz")==0) return AdjointzEnum;
-         else stage=4;
-   }
-   if(stage==4){
-              if (strcmp(name,"BalancethicknessMisfit")==0) return BalancethicknessMisfitEnum;
+              else if (strcmp(name,"BalancethicknessMisfit")==0) return BalancethicknessMisfitEnum;
               else if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum;
               else if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum;
@@ -501 +507 @@
               else if (strcmp(name,"MinVz")==0) return MinVzEnum;
               else if (strcmp(name,"MaxVz")==0) return MaxVzEnum;
-              else if (strcmp(name,"MaxAbsVz")==0) return MaxAbsVzEnum;
+         else stage=5;
+   }
+   if(stage==5){
+              if (strcmp(name,"MaxAbsVz")==0) return MaxAbsVzEnum;
               else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum;
               else if (strcmp(name,"TotalSmb")==0) return TotalSmbEnum;
@@ -507 +516 @@
               else if (strcmp(name,"Absolute")==0) return AbsoluteEnum;
               else if (strcmp(name,"Incremental")==0) return IncrementalEnum;
-         else stage=5;
-   }
-   if(stage==5){
-              if (strcmp(name,"AgressiveMigration")==0) return AgressiveMigrationEnum;
+              else if (strcmp(name,"AgressiveMigration")==0) return AgressiveMigrationEnum;
               else if (strcmp(name,"None")==0) return NoneEnum;
               else if (strcmp(name,"SoftMigration")==0) return SoftMigrationEnum;

issm/trunk-jpl/src/c/solutions/AnalysisConfiguration.cpp

@@ -64 +64 @@
 
                 case HydrologySolutionEnum:
-                        numanalyses=3;
+                        numanalyses=5;
                         analyses=xNew<int>(numanalyses);
-                        analyses[0]=HydrologyAnalysisEnum;
-                        analyses[1]=SurfaceSlopeAnalysisEnum;
-                        analyses[2]=BedSlopeAnalysisEnum;
+                        analyses[0]=HydrologyShreveAnalysisEnum;
+                        analyses[1]=HydrologyDCInefficientAnalysisEnum;
+                        analyses[2]=HydrologyDCEfficientAnalysisEnum;
+                        analyses[3]=SurfaceSlopeAnalysisEnum;
+                        analyses[4]=BedSlopeAnalysisEnum;
                         break;
 

issm/trunk-jpl/src/c/solutions/hydrology_core.cpp

@@ -61 +61 @@
                 if (hydrology_model==HydrologyshreveEnum){
                         if(VerboseSolution()) _pprintLine_("   computing water column");
-                        femmodel->SetCurrentConfiguration(HydrologyAnalysisEnum);
+                        femmodel->SetCurrentConfiguration(HydrologyShreveAnalysisEnum);
                         solver_nonlinear(femmodel,modify_loads);
 
@@ -81 +81 @@
                 else if (hydrology_model==HydrologydcEnum){
                         if(VerboseSolution()) _pprintLine_("   computing water head");
-                        femmodel->SetCurrentConfiguration(HydrologyAnalysisEnum);
-                        solver_linear(femmodel);
+                        solver_hydro_nonlinear(femmodel);
                         if(save_results && ((i+1)%output_frequency==0 || (i+1)==nsteps)){
                                 if(VerboseSolution()) _pprintLine_("   saving results ");

issm/trunk-jpl/src/c/solvers/solver_thermal_nonlinear.cpp

@@ -3 +3 @@
  */ 
 
+#include "./solvers.h"
 #include "../toolkits/toolkits.h"
 #include "../classes/objects/objects.h"
@@ -32 +33 @@
 
         /*parameters:*/
-        int kflag,pflag;
         int  configuration_type;
 
         /*Recover parameters: */
-        kflag=1; pflag=1;
         femmodel->parameters->FindParam(&thermal_penalty_threshold,ThermalPenaltyThresholdEnum);
         femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);

issm/trunk-jpl/src/c/solvers/solvers.h

@@ -13 +13 @@
 
 void solver_thermal_nonlinear(FemModel* femmodel);
+void solver_hydro_nonlinear(FemModel* femmodel);
 void solver_nonlinear(FemModel* femmodel,bool conserve_loads);
 void solver_newton(FemModel* femmodel);