Index: /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h (revision 14768) +++ /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h (revision 14769) @@ -103,4 +103,8 @@ HydrologydcEplTransmitivityEnum, HydrologydcIsefficientlayerEnum, + HydrologyLayerEnum, + HydrologySedimentEnum, + HydrologyEfficientEnum, + HydrologySedimentKmaxEnum, IndependentObjectEnum, InversionControlParametersEnum, @@ -274,5 +278,7 @@ FlaimAnalysisEnum, FlaimSolutionEnum, - HydrologyAnalysisEnum, + HydrologyShreveAnalysisEnum, + HydrologyDCInefficientAnalysisEnum, + HydrologyDCEfficientAnalysisEnum, HydrologySolutionEnum, MeltingAnalysisEnum, Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 14768) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 14769) @@ -454,10 +454,21 @@ #}}} #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 {{{ Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp (revision 14768) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp (revision 14769) @@ -252,6 +252,12 @@ 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,6 +615,12 @@ 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,6 +1395,6 @@ #endif #ifdef _HAVE_HYDROLOGY_ - case HydrologyAnalysisEnum: - GetSolutionFromInputsHydrology(solution); + case HydrologyShreveAnalysisEnum: + GetSolutionFromInputsHydrologyShreve(solution); break; #endif @@ -1721,22 +1733,28 @@ #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,28 +5820,5 @@ } /*}}}*/ -/*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,7 +5876,7 @@ 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,6 +5924,6 @@ } /*}}}*/ -/*FUNCTION Tria::CreateKMatrixHydrologyDCsediment{{{*/ -ElementMatrix* Tria::CreateKMatrixHydrologyDCsediment(void){ +/*FUNCTION Tria::CreateKMatrixHydrologyDCInefficient{{{*/ +ElementMatrix* Tria::CreateKMatrixHydrologyDCInefficient(void){ /*constants: */ @@ -5990,6 +5985,6 @@ } /*}}}*/ -/*FUNCTION Tria::CreateKMatrixHydrologyDCepl{{{*/ -ElementMatrix* Tria::CreateKMatrixHydrologyDCepl(void){ +/*FUNCTION Tria::CreateKMatrixHydrologyDCEfficient{{{*/ +ElementMatrix* Tria::CreateKMatrixHydrologyDCEfficient(void){ /*constants: */ @@ -6051,26 +6046,4 @@ } /*}}}*/ -/*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,6 +6095,6 @@ } /*}}}*/ -/*FUNCTION Tria::CreatePVectorHydrologyDCsediment {{{*/ -ElementVector* Tria::CreatePVectorHydrologyDCsediment(void){ +/*FUNCTION Tria::CreatePVectorHydrologyDCInefficient {{{*/ +ElementVector* Tria::CreatePVectorHydrologyDCInefficient(void){ /*Constants*/ @@ -6179,6 +6152,6 @@ } /*}}}*/ -/*FUNCTION Tria::CreatePVectorHydrologyDCepl {{{*/ -ElementVector* Tria::CreatePVectorHydrologyDCepl(void){ +/*FUNCTION Tria::CreatePVectorHydrologyDCEfficient {{{*/ +ElementVector* Tria::CreatePVectorHydrologyDCEfficient(void){ /*Constants*/ @@ -6236,14 +6209,14 @@ } /*}}}*/ -/*FUNCTION Tria::GetSolutionFromInputsHydrology{{{*/ -void Tria::GetSolutionFromInputsHydrology(Vector* solution){ +/*FUNCTION Tria::GetSolutionFromInputsHydrologyShreve{{{*/ +void Tria::GetSolutionFromInputsHydrologyShreve(Vector* 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,20 +6243,4 @@ delete gauss; xDelete(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,6 +6270,6 @@ } /*}}}*/ -/*FUNCTION Tria::InputUpdateFromSolutionHydrologyDC{{{*/ -void Tria::InputUpdateFromSolutionHydrologyDC(IssmDouble* solution){ +/*FUNCTION Tria::InputUpdateFromSolutionHydrologyDCInefficient{{{*/ +void Tria::InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution){ /*Intermediaries*/ @@ -6332,4 +6289,28 @@ /*Add input to the element: */ this->inputs->AddInput(new TriaP1Input(SedimentHeadEnum,values)); + + /*Free ressources:*/ + xDelete(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(values[i])) _error_("NaN found in solution vector"); + } + + /*Add input to the element: */ + this->inputs->AddInput(new TriaP1Input(EplHeadEnum,values)); /*Free ressources:*/ Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h (revision 14768) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h (revision 14769) @@ -239,17 +239,17 @@ #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* solution); void CreateHydrologyWaterVelocityInput(void); - void GetSolutionFromInputsHydrology(Vector* 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_ Index: /issm/trunk-jpl/src/c/classes/objects/Loads/Pengrid.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Loads/Pengrid.cpp (revision 14768) +++ /issm/trunk-jpl/src/c/classes/objects/Loads/Pengrid.cpp (revision 14769) @@ -246,4 +246,9 @@ break; #endif + #ifdef _HAVE_HYDROLOGY_ + case HydrologyDCInefficientAnalysisEnum: + Ke=PenaltyCreateKMatrixHydrologyDCInefficient(kmax); + break; + #endif default: _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet"); @@ -274,4 +279,9 @@ break; case DiagnosticHorizAnalysisEnum: case AdjointHorizAnalysisEnum: + break; + #endif + #ifdef _HAVE_HYDROLOGY_ + case HydrologyDCInefficientAnalysisEnum: + pe=PenaltyCreatePVectorHydrologyDCInefficient(kmax); break; #endif @@ -429,4 +439,5 @@ return; } + #ifdef _HAVE_THERMAL_ else if (analysis_type==ThermalAnalysisEnum){ ConstraintActivateThermal(punstable); @@ -436,4 +447,11 @@ return; } + #endif + #ifdef _HAVE_HYDROLOGY_ + else if (analysis_type==HydrologyDCInefficientAnalysisEnum){ + ConstraintActivateHydrologyDCInefficient(punstable); + return; + } + #endif else{ _error_("analysis: " << EnumToStringx(analysis_type) << " not supported yet"); @@ -442,4 +460,39 @@ } /*}}}*/ +#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,39 +563,4 @@ } /*}}}*/ -#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,4 +683,84 @@ /*}}}*/ #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){ Index: /issm/trunk-jpl/src/c/classes/objects/Loads/Pengrid.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Loads/Pengrid.h (revision 14768) +++ /issm/trunk-jpl/src/c/classes/objects/Loads/Pengrid.h (revision 14769) @@ -7,4 +7,9 @@ /*Headers:*/ /*{{{*/ +#ifdef HAVE_CONFIG_H +#include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif #include "./Load.h" class Hook; @@ -90,7 +95,12 @@ 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); Index: /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 14768) +++ /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 14769) @@ -108,4 +108,8 @@ 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,5 +281,7 @@ 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"; Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateDataSets.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateDataSets.cpp (revision 14768) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateDataSets.cpp (revision 14769) @@ -56,9 +56,21 @@ #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 Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 14768) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 14769) @@ -226,5 +226,7 @@ /*Solution specific parameters (FIXME: extend to other params)*/ #ifdef _HAVE_HYDROLOGY_ - CreateParametersHydrology(¶meters,iomodel,solution_type,analysis_type); + CreateParametersHydrologyShreve(¶meters,iomodel,solution_type,analysis_type); + CreateParametersHydrologyDCInefficient(¶meters,iomodel,solution_type,analysis_type); + CreateParametersHydrologyDCEfficient(¶meters,iomodel,solution_type,analysis_type); #endif Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/DistributeNumDofs.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/DistributeNumDofs.cpp (revision 14768) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/DistributeNumDofs.cpp (revision 14769) @@ -80,5 +80,11 @@ 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; } Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCEfficient/CMakeLists.txt =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCEfficient/CMakeLists.txt (revision 14769) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCEfficient/CMakeLists.txt (revision 14769) @@ -0,0 +1,8 @@ +# Subdirectories {{{ +# }}} +# HYDROLOGY_SOURCES {{{ +set(HYDROLOGY_SOURCES $ENV{ISSM_DIR}/src/c/modules/ModelProcessorx/HydrologyDCEfficient/CreateConstraintsHydrologyDCEfficient.cpp + $ENV{ISSM_DIR}/src/c/modules/ModelProcessorx/HydrologyDCEfficient/CreateLoadsHydrologyDCEfficient.cpp + $ENV{ISSM_DIR}/src/c/modules/ModelProcessorx/HydrologyDCEfficient/CreateNodesHydrologyDCEfficient.cpp + $ENV{ISSM_DIR}/src/c/modules/ModelProcessorx/HydrologyDCEfficient/UpdateElementsHydrologyDCEfficient.cpp PARENT_SCOPE) +# }}} Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCEfficient/CreateConstraintsHydrologyDCEfficient.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCEfficient/CreateConstraintsHydrologyDCEfficient.cpp (revision 14769) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCEfficient/CreateConstraintsHydrologyDCEfficient.cpp (revision 14769) @@ -0,0 +1,42 @@ +/* + * CreateConstraintsHydrologyDCEfficient.c: + */ + +#include "../../../Container/Container.h" +#include "../../../modules/modules.h" +#include "../../../toolkits/toolkits.h" +#include "../../../EnumDefinitions/EnumDefinitions.h" +#include "../../../io/io.h" +#include "../../../classes/objects/objects.h" +#include "../../../shared/shared.h" +#include "../ModelProcessorx.h" + +void CreateConstraintsHydrologyDCEfficient(Constraints** pconstraints, IoModel* iomodel){ + + /*Recover pointer: */ + bool isefficientlayer; + int hydrology_model; + Constraints* constraints=*pconstraints; + + /*Create constraints if they do not exist yet*/ + if(!constraints) constraints = new Constraints(); + + /*Do we really want DC?*/ + iomodel->Constant(&hydrology_model,HydrologyEnum); + if(hydrology_model!=HydrologydcEnum){ + *pconstraints=constraints; + return; + } + + /*Do we want an efficient layer*/ + iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum); + if(!isefficientlayer){ + *pconstraints=constraints; + return; + } + + IoModelToConstraintsx(constraints,iomodel,HydrologydcSpceplHeadEnum,HydrologyDCEfficientAnalysisEnum); + + /*Assign output pointer: */ + *pconstraints=constraints; +} Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCEfficient/CreateLoadsHydrologyDCEfficient.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCEfficient/CreateLoadsHydrologyDCEfficient.cpp (revision 14769) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCEfficient/CreateLoadsHydrologyDCEfficient.cpp (revision 14769) @@ -0,0 +1,48 @@ +/*! \file CreateLoadsHydrologyDCEfficient.c: + */ + +#include "../../../Container/Container.h" +#include "../../../toolkits/toolkits.h" +#include "../../../io/io.h" +#include "../../../EnumDefinitions/EnumDefinitions.h" +#include "../../../classes/objects/objects.h" +#include "../../../shared/shared.h" +#include "../../../include/include.h" +#include "../ModelProcessorx.h" + +void CreateLoadsHydrologyDCEfficient(Loads** ploads, IoModel* iomodel){ + + /*Intermediary*/ + bool isefficientlayer; + int hydrology_model; + int numberofvertices; + Pengrid *pengrid = NULL; + + /*Recover pointer: */ + Loads* loads=*ploads; + + /*Fetch parameters: */ + iomodel->Constant(&numberofvertices,MeshNumberofverticesEnum); + + /*Create loads if they do not exist yet*/ + if(!loads) loads = new Loads(); + + /*Do we really want DC?*/ + iomodel->Constant(&hydrology_model,HydrologyEnum); + if(hydrology_model!=HydrologydcEnum){ + *ploads=loads; + return; + } + + /*Do we want an efficient layer*/ + iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum); + if(!isefficientlayer){ + *ploads=loads; + return; + } + + /*Nothing for now*/ + + /*Assign output pointer: */ + *ploads=loads; +} Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCEfficient/CreateNodesHydrologyDCEfficient.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCEfficient/CreateNodesHydrologyDCEfficient.cpp (revision 14769) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCEfficient/CreateNodesHydrologyDCEfficient.cpp (revision 14769) @@ -0,0 +1,64 @@ +/* + * CreateNodesHydrologyDCEfficient.c: + */ + +#include "../../../Container/Container.h" +#include "../../../toolkits/toolkits.h" +#include "../../../io/io.h" +#include "../../../EnumDefinitions/EnumDefinitions.h" +#include "../../../classes/objects/objects.h" +#include "../../../include/include.h" +#include "../../../shared/shared.h" +#include "../../MeshPartitionx/MeshPartitionx.h" +#include "../ModelProcessorx.h" + +void CreateNodesHydrologyDCEfficient(Nodes** pnodes, IoModel* iomodel){ + + /*Intermediary*/ + int i; + bool isefficientlayer; + bool continuous_galerkin=true; + int hydrology_model; + int numberofvertices; + + /*Fetch parameters: */ + iomodel->Constant(&numberofvertices,MeshNumberofverticesEnum); + + /*Recover pointer: */ + Nodes* nodes=*pnodes; + + /*Create nodes if they do not exist yet*/ + if(!nodes) nodes = new Nodes(); + + /*Now, do we really want DC?*/ + iomodel->Constant(&hydrology_model,HydrologyEnum); + if(hydrology_model!=HydrologydcEnum){ + *pnodes=nodes; + return; + } + + /*Do we want an efficient layer*/ + iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum); + if(!isefficientlayer){ + *pnodes=nodes; + return; + } + + /*Continuous Galerkin partition of nodes: */ + NodesPartitioning(&iomodel->my_nodes,iomodel->my_elements,iomodel->my_vertices,iomodel,continuous_galerkin); + + /*Create nodes and vertices: */ + iomodel->FetchData(6,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,FlowequationVertexEquationEnum,MaskVertexonwaterEnum); + for (i=0;imy_vertices[i]){ + /*Add node to nodes dataset: */ + nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i+1,i,iomodel,HydrologyDCEfficientAnalysisEnum)); + + } + } + iomodel->DeleteData(6,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,FlowequationVertexEquationEnum,MaskVertexonwaterEnum); + + /*Assign output pointer: */ + *pnodes=nodes; +} Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCEfficient/CreateParametersHydrologyDCEfficient.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCEfficient/CreateParametersHydrologyDCEfficient.cpp (revision 14769) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCEfficient/CreateParametersHydrologyDCEfficient.cpp (revision 14769) @@ -0,0 +1,44 @@ +/*!\file: CreateParametersHydrologyDCEfficient.cpp + * \brief driver for creating parameters dataset, for control analysis. + */ + +#include "../../../Container/Container.h" +#include "../../../toolkits/toolkits.h" +#include "../../../io/io.h" +#include "../../../EnumDefinitions/EnumDefinitions.h" +#include "../../../classes/objects/objects.h" +#include "../../../shared/shared.h" +#include "../../../include/include.h" +#include "../ModelProcessorx.h" + +void CreateParametersHydrologyDCEfficient(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type){ + + Parameters *parameters = NULL; + int hydrology_model; + bool isefficientlayer; + + /*Get parameters: */ + parameters=*pparameters; + + /*retrieve some parameters: */ + iomodel->Constant(&hydrology_model,HydrologyEnum); + + /*Now, do we really want DC?*/ + if(hydrology_model!=HydrologydcEnum){ + *pparameters=parameters; + return; + } + + /*Do we want an efficient layer*/ + iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum); + if(!isefficientlayer){ + *pparameters=parameters; + return; + } + + + /*Nothing for now*/ + + /*Assign output pointer: */ + *pparameters=parameters; +} Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCEfficient/UpdateElementsHydrologyDCEfficient.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCEfficient/UpdateElementsHydrologyDCEfficient.cpp (revision 14769) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCEfficient/UpdateElementsHydrologyDCEfficient.cpp (revision 14769) @@ -0,0 +1,56 @@ +/* + * UpdateElementsHydrologyDCEfficient: + */ + +#include "../../../Container/Container.h" +#include "../../../toolkits/toolkits.h" +#include "../../../io/io.h" +#include "../../../modules/modules.h" +#include "../../../EnumDefinitions/EnumDefinitions.h" +#include "../../../classes/objects/objects.h" +#include "../../../shared/shared.h" +#include "../../MeshPartitionx/MeshPartitionx.h" +#include "../../../include/include.h" +#include "../ModelProcessorx.h" + +void UpdateElementsHydrologyDCEfficient(Elements* elements, IoModel* iomodel,int analysis_counter,int analysis_type){ + + bool isefficientlayer; + int hydrology_model; + int numberofelements; + + /*Fetch data needed: */ + iomodel->Constant(&numberofelements,MeshNumberofelementsEnum); + iomodel->FetchData(1,MeshElementsEnum); + + /*Now, do we really want DC?*/ + iomodel->Constant(&hydrology_model,HydrologyEnum); + if(hydrology_model!=HydrologydcEnum) return; + + /*Do we want an efficient layer*/ + iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum); + if(!isefficientlayer) return; + + /*Update elements: */ + int counter=0; + for(int i=0;imy_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type); //we need i to index into elements. + counter++; + } + } + + iomodel->FetchDataToInput(elements,ThicknessEnum); + iomodel->FetchDataToInput(elements,SurfaceEnum); + iomodel->FetchDataToInput(elements,BedEnum); + iomodel->FetchDataToInput(elements,MaskElementonfloatingiceEnum); + iomodel->FetchDataToInput(elements,MeshElementonbedEnum); + iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum); + iomodel->FetchDataToInput(elements,MaskElementonwaterEnum); + iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum); + iomodel->FetchDataToInput(elements,SedimentHeadEnum); + + /*Free data: */ + iomodel->DeleteData(1,MeshElementsEnum); +} Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CMakeLists.txt =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CMakeLists.txt (revision 14769) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CMakeLists.txt (revision 14769) @@ -0,0 +1,8 @@ +# Subdirectories {{{ +# }}} +# HYDROLOGY_SOURCES {{{ +set(HYDROLOGY_SOURCES $ENV{ISSM_DIR}/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CreateConstraintsHydrologyDCInefficient.cpp + $ENV{ISSM_DIR}/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CreateLoadsHydrologyDCInefficient.cpp + $ENV{ISSM_DIR}/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CreateNodesHydrologyDCInefficient.cpp + $ENV{ISSM_DIR}/src/c/modules/ModelProcessorx/HydrologyDCInefficient/UpdateElementsHydrologyDCInefficient.cpp PARENT_SCOPE) +# }}} Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CreateConstraintsHydrologyDCInefficient.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CreateConstraintsHydrologyDCInefficient.cpp (revision 14769) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CreateConstraintsHydrologyDCInefficient.cpp (revision 14769) @@ -0,0 +1,35 @@ +/* + * CreateConstraintsHydrologyDCInefficient.c: + */ + +#include "../../../Container/Container.h" +#include "../../../modules/modules.h" +#include "../../../toolkits/toolkits.h" +#include "../../../EnumDefinitions/EnumDefinitions.h" +#include "../../../io/io.h" +#include "../../../classes/objects/objects.h" +#include "../../../shared/shared.h" +#include "../ModelProcessorx.h" + +void CreateConstraintsHydrologyDCInefficient(Constraints** pconstraints, IoModel* iomodel){ + + /*Recover pointer: */ + int hydrology_model; + Constraints* constraints=*pconstraints; + + /*retrieve some parameters: */ + iomodel->Constant(&hydrology_model,HydrologyEnum); + + /*Create constraints if they do not exist yet*/ + if(!constraints) constraints = new Constraints(); + + if(hydrology_model!=HydrologydcEnum){ + *pconstraints=constraints; + return; + } + + IoModelToConstraintsx(constraints,iomodel,HydrologydcSpcsedimentHeadEnum,HydrologyDCInefficientAnalysisEnum); + + /*Assign output pointer: */ + *pconstraints=constraints; +} Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CreateLoadsHydrologyDCInefficient.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CreateLoadsHydrologyDCInefficient.cpp (revision 14769) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CreateLoadsHydrologyDCInefficient.cpp (revision 14769) @@ -0,0 +1,49 @@ +/*! \file CreateLoadsHydrologyDCInefficient.c: + */ + +#include "../../../Container/Container.h" +#include "../../../toolkits/toolkits.h" +#include "../../../io/io.h" +#include "../../../EnumDefinitions/EnumDefinitions.h" +#include "../../../classes/objects/objects.h" +#include "../../../shared/shared.h" +#include "../../../include/include.h" +#include "../ModelProcessorx.h" + +void CreateLoadsHydrologyDCInefficient(Loads** ploads, IoModel* iomodel){ + + /*Intermediary*/ + int hydrology_model; + int numberofvertices; + Pengrid *pengrid = NULL; + + /*Recover pointer: */ + Loads* loads=*ploads; + + /*Fetch parameters: */ + iomodel->Constant(&hydrology_model,HydrologyEnum); + iomodel->Constant(&numberofvertices,MeshNumberofverticesEnum); + + /*Create loads if they do not exist yet*/ + if(!loads) loads = new Loads(); + + if(hydrology_model!=HydrologydcEnum){ + *ploads=loads; + return; + } + + //create penalties for nodes: no node can have a temperature over the melting point + iomodel->FetchData(1,MeshElementsEnum); + CreateSingleNodeToElementConnectivity(iomodel); + + for(int i=0;imy_vertices[i]==1)){ + loads->AddObject(new Pengrid(iomodel->loadcounter+i+1,i,iomodel,HydrologyDCInefficientAnalysisEnum)); + } + } + iomodel->DeleteData(1,MeshElementsEnum); + + /*Assign output pointer: */ + *ploads=loads; +} Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CreateNodesHydrologyDCInefficient.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CreateNodesHydrologyDCInefficient.cpp (revision 14769) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CreateNodesHydrologyDCInefficient.cpp (revision 14769) @@ -0,0 +1,56 @@ +/* + * CreateNodesHydrologyDCInefficient.c: + */ + +#include "../../../Container/Container.h" +#include "../../../toolkits/toolkits.h" +#include "../../../io/io.h" +#include "../../../EnumDefinitions/EnumDefinitions.h" +#include "../../../classes/objects/objects.h" +#include "../../../include/include.h" +#include "../../../shared/shared.h" +#include "../../MeshPartitionx/MeshPartitionx.h" +#include "../ModelProcessorx.h" + +void CreateNodesHydrologyDCInefficient(Nodes** pnodes, IoModel* iomodel){ + + /*Intermediary*/ + int i; + bool continuous_galerkin=true; + int hydrology_model; + int numberofvertices; + + /*Fetch parameters: */ + iomodel->Constant(&numberofvertices,MeshNumberofverticesEnum); + iomodel->Constant(&hydrology_model,HydrologyEnum); + + /*Recover pointer: */ + Nodes* nodes=*pnodes; + + /*Create nodes if they do not exist yet*/ + if(!nodes) nodes = new Nodes(); + + /*Now, do we really want DC?*/ + if(hydrology_model!=HydrologydcEnum){ + *pnodes=nodes; + return; + } + + /*Continuous Galerkin partition of nodes: */ + NodesPartitioning(&iomodel->my_nodes,iomodel->my_elements,iomodel->my_vertices,iomodel,continuous_galerkin); + + /*Create nodes and vertices: */ + iomodel->FetchData(6,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,FlowequationVertexEquationEnum,MaskVertexonwaterEnum); + for (i=0;imy_vertices[i]){ + /*Add node to nodes dataset: */ + nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i+1,i,iomodel,HydrologyDCInefficientAnalysisEnum)); + + } + } + iomodel->DeleteData(6,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,FlowequationVertexEquationEnum,MaskVertexonwaterEnum); + + /*Assign output pointer: */ + *pnodes=nodes; +} Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CreateParametersHydrologyDCInefficient.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CreateParametersHydrologyDCInefficient.cpp (revision 14769) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CreateParametersHydrologyDCInefficient.cpp (revision 14769) @@ -0,0 +1,38 @@ +/*!\file: CreateParametersHydrologyDCInefficient.cpp + * \brief driver for creating parameters dataset, for control analysis. + */ + +#include "../../../Container/Container.h" +#include "../../../toolkits/toolkits.h" +#include "../../../io/io.h" +#include "../../../EnumDefinitions/EnumDefinitions.h" +#include "../../../classes/objects/objects.h" +#include "../../../shared/shared.h" +#include "../../../include/include.h" +#include "../ModelProcessorx.h" + +void CreateParametersHydrologyDCInefficient(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type){ + + Parameters *parameters = NULL; + int hydrology_model; + bool isefficientlayer; + + /*Get parameters: */ + parameters=*pparameters; + + /*retrieve some parameters: */ + iomodel->Constant(&hydrology_model,HydrologyEnum); + + /*Now, do we really want DC?*/ + if(hydrology_model!=HydrologydcEnum){ + *pparameters=parameters; + return; + } + + iomodel->FetchData(&isefficientlayer,HydrologydcIsefficientlayerEnum); + parameters->AddObject(new IntParam(HydrologyEnum,hydrology_model)); + parameters->AddObject(new BoolParam(HydrologydcIsefficientlayerEnum,isefficientlayer)); + + /*Assign output pointer: */ + *pparameters=parameters; +} Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/UpdateElementsHydrologyDCInefficient.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/UpdateElementsHydrologyDCInefficient.cpp (revision 14769) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/UpdateElementsHydrologyDCInefficient.cpp (revision 14769) @@ -0,0 +1,51 @@ +/* + * UpdateElementsHydrologyDCInefficient: + */ + +#include "../../../Container/Container.h" +#include "../../../toolkits/toolkits.h" +#include "../../../io/io.h" +#include "../../../modules/modules.h" +#include "../../../EnumDefinitions/EnumDefinitions.h" +#include "../../../classes/objects/objects.h" +#include "../../../shared/shared.h" +#include "../../MeshPartitionx/MeshPartitionx.h" +#include "../../../include/include.h" +#include "../ModelProcessorx.h" + +void UpdateElementsHydrologyDCInefficient(Elements* elements, IoModel* iomodel,int analysis_counter,int analysis_type){ + + int hydrology_model; + int numberofelements; + + /*Fetch data needed: */ + iomodel->Constant(&hydrology_model,HydrologyEnum); + iomodel->Constant(&numberofelements,MeshNumberofelementsEnum); + iomodel->FetchData(1,MeshElementsEnum); + + /*Now, do we really want DC?*/ + if(hydrology_model!=HydrologydcEnum) return; + + /*Update elements: */ + int counter=0; + for(int i=0;imy_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type); //we need i to index into elements. + counter++; + } + } + + iomodel->FetchDataToInput(elements,ThicknessEnum); + iomodel->FetchDataToInput(elements,SurfaceEnum); + iomodel->FetchDataToInput(elements,BedEnum); + iomodel->FetchDataToInput(elements,MaskElementonfloatingiceEnum); + iomodel->FetchDataToInput(elements,MeshElementonbedEnum); + iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum); + iomodel->FetchDataToInput(elements,MaskElementonwaterEnum); + iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum); + iomodel->FetchDataToInput(elements,SedimentHeadEnum); + + /*Free data: */ + iomodel->DeleteData(1,MeshElementsEnum); +} Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/CMakeLists.txt =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/CMakeLists.txt (revision 14769) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/CMakeLists.txt (revision 14769) @@ -0,0 +1,8 @@ +# Subdirectories {{{ +# }}} +# HYDROLOGY_SOURCES {{{ +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) +# }}} Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/CreateConstraintsHydrologyShreve.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/CreateConstraintsHydrologyShreve.cpp (revision 14769) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/CreateConstraintsHydrologyShreve.cpp (revision 14769) @@ -0,0 +1,36 @@ +/* + * CreateConstraintsHydrologyShreve.c: + */ + +#include "../../../Container/Container.h" +#include "../../../modules/modules.h" +#include "../../../toolkits/toolkits.h" +#include "../../../EnumDefinitions/EnumDefinitions.h" +#include "../../../io/io.h" +#include "../../../classes/objects/objects.h" +#include "../../../shared/shared.h" +#include "../ModelProcessorx.h" + +void CreateConstraintsHydrologyShreve(Constraints** pconstraints, IoModel* iomodel){ + + /*Recover pointer: */ + int hydrology_model; + bool isefficientlayer; + Constraints* constraints=*pconstraints; + + /*retrieve some parameters: */ + iomodel->Constant(&hydrology_model,HydrologyEnum); + + /*Create constraints if they do not exist yet*/ + if(!constraints) constraints = new Constraints(); + + if(hydrology_model!=HydrologyshreveEnum){ + *pconstraints=constraints; + return; + } + + IoModelToConstraintsx(constraints,iomodel,HydrologyshreveSpcwatercolumnEnum,HydrologyShreveAnalysisEnum); + + /*Assign output pointer: */ + *pconstraints=constraints; +} Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/CreateLoadsHydrologyShreve.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/CreateLoadsHydrologyShreve.cpp (revision 14769) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/CreateLoadsHydrologyShreve.cpp (revision 14769) @@ -0,0 +1,27 @@ +/*! \file CreateLoadsHydrologyShreve.c: + */ + +#include "../../../Container/Container.h" +#include "../../../toolkits/toolkits.h" +#include "../../../io/io.h" +#include "../../../EnumDefinitions/EnumDefinitions.h" +#include "../../../classes/objects/objects.h" +#include "../../../shared/shared.h" +#include "../../../include/include.h" +#include "../ModelProcessorx.h" + +void CreateLoadsHydrologyShreve(Loads** ploads, IoModel* iomodel){ + + /*Intermediary*/ + int numberofvertices; + Pengrid *pengrid = NULL; + + /*Recover pointer: */ + Loads* loads=*ploads; + + /*Create loads if they do not exist yet*/ + if(!loads) loads = new Loads(); + + /*Assign output pointer: */ + *ploads=loads; +} Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/CreateNodesHydrologyShreve.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/CreateNodesHydrologyShreve.cpp (revision 14769) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/CreateNodesHydrologyShreve.cpp (revision 14769) @@ -0,0 +1,55 @@ +/* + * CreateNodesHydrologyShreve.c: + */ + +#include "../../../Container/Container.h" +#include "../../../toolkits/toolkits.h" +#include "../../../io/io.h" +#include "../../../EnumDefinitions/EnumDefinitions.h" +#include "../../../classes/objects/objects.h" +#include "../../../include/include.h" +#include "../../../shared/shared.h" +#include "../../MeshPartitionx/MeshPartitionx.h" +#include "../ModelProcessorx.h" + +void CreateNodesHydrologyShreve(Nodes** pnodes, IoModel* iomodel){ + + /*Intermediary*/ + int i; + int hydrology_model; + bool continuous_galerkin=true; + int numberofvertices; + + /*Fetch parameters: */ + iomodel->Constant(&numberofvertices,MeshNumberofverticesEnum); + iomodel->Constant(&hydrology_model,HydrologyEnum); + + /*Recover pointer: */ + Nodes* nodes=*pnodes; + + /*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: */ + NodesPartitioning(&iomodel->my_nodes,iomodel->my_elements,iomodel->my_vertices,iomodel,continuous_galerkin); + + /*Create nodes and vertices: */ + iomodel->FetchData(6,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,FlowequationVertexEquationEnum,MaskVertexonwaterEnum); + for (i=0;imy_vertices[i]){ + /*Add node to nodes dataset: */ + nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i+1,i,iomodel,HydrologyShreveAnalysisEnum)); + } + } + iomodel->DeleteData(6,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,FlowequationVertexEquationEnum,MaskVertexonwaterEnum); + + /*Assign output pointer: */ + *pnodes=nodes; +} Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/CreateParametersHydrologyShreve.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/CreateParametersHydrologyShreve.cpp (revision 14769) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/CreateParametersHydrologyShreve.cpp (revision 14769) @@ -0,0 +1,37 @@ +/*!\file: CreateParametersHydrologyShreve.cpp + * \brief driver for creating parameters dataset, for control analysis. + */ + +#include "../../../Container/Container.h" +#include "../../../toolkits/toolkits.h" +#include "../../../io/io.h" +#include "../../../EnumDefinitions/EnumDefinitions.h" +#include "../../../classes/objects/objects.h" +#include "../../../shared/shared.h" +#include "../../../include/include.h" +#include "../ModelProcessorx.h" + +void CreateParametersHydrologyShreve(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type){ + + Parameters *parameters = NULL; + int hydrology_model; + bool isefficientlayer; + + /*Get parameters: */ + parameters=*pparameters; + + /*retrieve some parameters: */ + 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)); + parameters->AddObject(iomodel->CopyConstantObject(HydrologyshreveStabilizationEnum)); + + /*Assign output pointer: */ + *pparameters=parameters; +} Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/UpdateElementsHydrologyShreve.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/UpdateElementsHydrologyShreve.cpp (revision 14769) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/UpdateElementsHydrologyShreve.cpp (revision 14769) @@ -0,0 +1,53 @@ +/* + * UpdateElementsHydrologyShreve: + */ + +#include "../../../Container/Container.h" +#include "../../../toolkits/toolkits.h" +#include "../../../io/io.h" +#include "../../../modules/modules.h" +#include "../../../EnumDefinitions/EnumDefinitions.h" +#include "../../../classes/objects/objects.h" +#include "../../../shared/shared.h" +#include "../../MeshPartitionx/MeshPartitionx.h" +#include "../../../include/include.h" +#include "../ModelProcessorx.h" + +void UpdateElementsHydrologyShreve(Elements* elements, IoModel* iomodel,int analysis_counter,int analysis_type){ + + int hydrology_model; + int numberofelements; + + /*Fetch data needed: */ + iomodel->Constant(&hydrology_model,HydrologyEnum); + iomodel->Constant(&numberofelements,MeshNumberofelementsEnum); + iomodel->FetchData(1,MeshElementsEnum); + + /*Now, do we really want Shreve?*/ + if(hydrology_model!=HydrologyshreveEnum) return; + + /*Update elements: */ + int counter=0; + for(int i=0;imy_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type); //we need i to index into elements. + counter++; + } + } + + iomodel->FetchDataToInput(elements,ThicknessEnum); + iomodel->FetchDataToInput(elements,SurfaceEnum); + iomodel->FetchDataToInput(elements,BedEnum); + iomodel->FetchDataToInput(elements,MaskElementonfloatingiceEnum); + iomodel->FetchDataToInput(elements,MeshElementonbedEnum); + iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum); + iomodel->FetchDataToInput(elements,MaskElementonwaterEnum); + iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum); + iomodel->FetchDataToInput(elements,WatercolumnEnum); + + elements->InputDuplicate(WatercolumnEnum,WaterColumnOldEnum); + + /*Free data: */ + iomodel->DeleteData(1,MeshElementsEnum); +} Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 14768) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 14769) @@ -24,5 +24,7 @@ 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,9 +81,19 @@ 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:*/ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp (revision 14768) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp (revision 14769) @@ -14,8 +14,7 @@ /*Intermediary*/ - int i; - int dim; - int numberofvertices; - Pengrid* pengrid = NULL; + int dim; + int numberofvertices; + Pengrid *pengrid = NULL; /*Recover pointer: */ @@ -36,5 +35,5 @@ CreateSingleNodeToElementConnectivity(iomodel); - for (i=0;i(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; Index: /issm/trunk-jpl/src/c/solutions/hydrology_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/hydrology_core.cpp (revision 14768) +++ /issm/trunk-jpl/src/c/solutions/hydrology_core.cpp (revision 14769) @@ -61,5 +61,5 @@ if (hydrology_model==HydrologyshreveEnum){ if(VerboseSolution()) _pprintLine_(" computing water column"); - femmodel->SetCurrentConfiguration(HydrologyAnalysisEnum); + femmodel->SetCurrentConfiguration(HydrologyShreveAnalysisEnum); solver_nonlinear(femmodel,modify_loads); @@ -81,6 +81,5 @@ 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 "); Index: /issm/trunk-jpl/src/c/solvers/solver_hydro_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solvers/solver_hydro_nonlinear.cpp (revision 14769) +++ /issm/trunk-jpl/src/c/solvers/solver_hydro_nonlinear.cpp (revision 14769) @@ -0,0 +1,124 @@ +/*!\file: solver_hydro_nonlinear.cpp + * \brief: core of the hydro solution + */ + +#include "./solvers.h" +#include "../toolkits/toolkits.h" +#include "../classes/objects/objects.h" +#include "../io/io.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../modules/modules.h" + +void solver_hydro_nonlinear(FemModel* femmodel){ + + /*solution : */ + Vector* ug=NULL; + Vector* uf=NULL; + Vector* uf_old=NULL; + Vector* ys=NULL; + IssmDouble sediment_kmax; + + /*intermediary: */ + Matrix* Kff=NULL; + Matrix* Kfs=NULL; + Vector* pf=NULL; + Vector* df=NULL; + + bool converged,isefficientlayer; + int constraints_converged; + int num_unstable_constraints; + int count; + int hydro_maxiter; + + /*Recover parameters: */ + femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); + hydro_maxiter=100; + count=1; + converged=false; + + for(;;){ + + /*First layer*/ + femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum); + InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,true,ResetPenaltiesEnum); + InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,false,ConvergedEnum); + femmodel->UpdateConstraintsx(); + femmodel->parameters->SetParam(HydrologySedimentEnum,HydrologyLayerEnum); + for(;;){ + femmodel->SystemMatricesx(&Kff, &Kfs, &pf,&df, &sediment_kmax); + CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCInefficientAnalysisEnum); + Reduceloadx(pf,Kfs,ys); xdelete(&Kfs); xdelete(&uf); + Solverx(&uf, Kff, pf,uf_old, df, femmodel->parameters); + xdelete(&uf_old); uf_old=uf->Duplicate(); + xdelete(&Kff);xdelete(&pf);xdelete(&ug); xdelete(&df); + Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters); xdelete(&ys); + InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug); + + ConstraintsStatex(&constraints_converged,&num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + + if (!converged){ + if(VerboseConvergence()) _pprintLine_(" #unstable constraints = " << num_unstable_constraints); + if(num_unstable_constraints==0) converged = true; + if (count>=hydro_maxiter){ + _error_(" maximum number of iterations (" << hydro_maxiter << ") exceeded"); + } + } + count++; + + if(converged){ + InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,converged,ConvergedEnum); + InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sediment_kmax,HydrologySedimentKmaxEnum); + InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug); + break; + } + } + + /*Second layer*/ + if(!isefficientlayer) break; + femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum); + InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,true,ResetPenaltiesEnum); + InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,false,ConvergedEnum); + femmodel->UpdateConstraintsx(); + femmodel->parameters->SetParam(HydrologyEfficientEnum,HydrologyLayerEnum); + converged = false; + for(;;){ + femmodel->SystemMatricesx(&Kff, &Kfs, &pf,&df,NULL); + CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCEfficientAnalysisEnum); + Reduceloadx(pf,Kfs,ys); xdelete(&Kfs); xdelete(&uf); + Solverx(&uf, Kff, pf,uf_old, df, femmodel->parameters); + xdelete(&uf_old); uf_old=uf->Duplicate(); + xdelete(&Kff);xdelete(&pf);xdelete(&ug); xdelete(&df); + Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters); xdelete(&ys); + InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug); + + ConstraintsStatex(&constraints_converged,&num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + + if (!converged){ + if(VerboseConvergence()) _pprintLine_(" #unstable constraints = " << num_unstable_constraints); + if(num_unstable_constraints==0) converged = true; + if (count>=hydro_maxiter){ + _error_(" maximum number of iterations (" << hydro_maxiter << ") exceeded"); + } + } + count++; + + if(converged){ + InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,converged,ConvergedEnum); + InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sediment_kmax,MeltingOffsetEnum); + InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug); + break; + } + } + + /*System convergence check*/ + break; + } + + InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); + + /*Free ressources: */ + xdelete(&ug); + xdelete(&uf); + xdelete(&uf_old); + xdelete(&ys); +} Index: /issm/trunk-jpl/src/c/solvers/solver_thermal_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solvers/solver_thermal_nonlinear.cpp (revision 14768) +++ /issm/trunk-jpl/src/c/solvers/solver_thermal_nonlinear.cpp (revision 14769) @@ -3,4 +3,5 @@ */ +#include "./solvers.h" #include "../toolkits/toolkits.h" #include "../classes/objects/objects.h" @@ -32,9 +33,7 @@ /*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); Index: /issm/trunk-jpl/src/c/solvers/solvers.h =================================================================== --- /issm/trunk-jpl/src/c/solvers/solvers.h (revision 14768) +++ /issm/trunk-jpl/src/c/solvers/solvers.h (revision 14769) @@ -13,4 +13,5 @@ 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);