Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 17360) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 17361) @@ -410,6 +410,8 @@ break; case 1: - - element->GetInputValue(&epl_thickness,gauss,HydrologydcEplThicknessEnum); + /* _assert_(input) */ + /* get input */ + + element->GetInputValue(&epl_thickness,gauss,HydrologydcEplThicknessEnum); element->GetInputValue(&sed_head,gauss,SedimentHeadEnum); element->GetInputValue(&epl_head,gauss,EplHeadEnum); @@ -481,4 +483,89 @@ return transfer; }/*}}}*/ + +void HydrologyDCEfficientAnalysis::ComputeEPLThickness(FemModel* femmodel){/*{{{*/ + + bool active_element; + int meshtype; + IssmDouble dt,A,B; + IssmDouble EPLgrad2; + IssmDouble EPL_N; + + femmodel->parameters->FindParam(&meshtype,MeshTypeEnum); + + for(int j=0;jelements->Size();j++){ + + Element* element=(Element*)femmodel->elements->GetObjectByOffset(j); + + switch(meshtype){ + case Mesh2DhorizontalEnum: + if(!element->IsOnBed()) return; + B = element->GetMaterialParameter(MaterialsRheologyBbarEnum); + break; + case Mesh3DEnum: + B = element->GetMaterialParameter(MaterialsRheologyBEnum); + break; + default: + _error_("not Implemented Yet"); + } + + int numnodes = element->GetNumberOfNodes(); + IssmDouble thickness[numnodes]; + IssmDouble eplhead[numnodes]; + IssmDouble epl_slopeX[numnodes],epl_slopeY[numnodes]; + IssmDouble old_thickness[numnodes]; + IssmDouble ice_thickness[numnodes],bed[numnodes]; + + Input* active_element_input=element->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); + active_element_input->GetInputValue(&active_element); + element->FindParam(&dt,TimesteppingTimeStepEnum); + + /*For now, assuming just one way to compute EPL thickness*/ + IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum); + IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum); + IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble n = element->GetMaterialParameter(MaterialsRheologyNEnum); + IssmDouble latentheat = element->GetMaterialParameter(MaterialsLatentheatEnum); + IssmDouble epl_conductivity = element->GetMaterialParameter(HydrologydcEplConductivityEnum); + IssmDouble init_thick = element->GetMaterialParameter(HydrologydcEplInitialThicknessEnum); + + A=pow(B,-n); + + element->GetInputListOnVertices(&eplhead[0],EplHeadEnum); + element->GetInputListOnVertices(&epl_slopeX[0],EplHeadSlopeXEnum); + element->GetInputListOnVertices(&epl_slopeY[0],EplHeadSlopeYEnum); + element->GetInputListOnVertices(&old_thickness[0],HydrologydcEplThicknessOldEnum); + element->GetInputListOnVertices(&ice_thickness[0],ThicknessEnum); + element->GetInputListOnVertices(&bed[0],BedEnum); + + if(!active_element){ + + /*Keeping thickness to initial value if EPL is not active*/ + for(int i=0;i10.0*init_thick){ + thickness[i] = 10.0*init_thick; + } + } + } + element->AddInput(HydrologydcEplThicknessEnum,thickness,P1Enum); + } +} +/*}}}*/ void HydrologyDCEfficientAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h (revision 17360) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h (revision 17361) @@ -37,4 +37,5 @@ IssmDouble GetHydrologyKMatrixTransfer(Element* element, Gauss* gauss); IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss); + void ComputeEPLThickness(FemModel* femmodel); }; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 17360) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 17361) @@ -286,5 +286,5 @@ virtual void HydrologyEPLGetMask(Vector* mask)=0; virtual void HydrologyEPLGetActive(Vector* active)=0; - virtual void ComputeEPLThickness(void)=0; + // virtual void ComputeEPLThickness(void)=0; virtual void MigrateGroundingLine(IssmDouble* sheet_ungrounding)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 17360) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 17361) @@ -4220,80 +4220,80 @@ } /*}}}*/ -/*FUNCTION Penta::ComputeEPLThickness{{{*/ -void Penta::ComputeEPLThickness(void){ - - int i; - const int numdof = NDOF1 *NUMVERTICES; - const int numdof2d = NDOF1 *NUMVERTICES2D; - bool isefficientlayer; - IssmDouble n,A,dt,init_thick; - IssmDouble rho_water,rho_ice; - IssmDouble gravity,latentheat,EPLgrad; - IssmDouble EPL_N,epl_conductivity; - IssmDouble activeEpl[numdof],thickness[numdof]; - IssmDouble eplhead[numdof], old_thickness[numdof]; - IssmDouble epl_slopeX[numdof],epl_slopeY[numdof]; - IssmDouble ice_thickness[numdof],bed[numdof]; - Penta *penta = NULL; - /*If not on bed, return*/ - if (!IsOnBed())return; - - /*Get the flag to know if the efficient layer is present*/ - this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - - if(isefficientlayer){ - /*For now, assuming just one way to compute EPL thickness*/ - rho_water = matpar->GetRhoWater(); - rho_ice = matpar->GetRhoIce(); - gravity = matpar->GetG(); - latentheat = matpar->GetLatentHeat(); - epl_conductivity = matpar->GetEplConductivity(); - init_thick = matpar->GetEplInitialThickness(); - n = material->GetN(); - A = material->GetA(); +/* /\*FUNCTION Penta::ComputeEPLThickness{{{*\/ */ +/* void Penta::ComputeEPLThickness(void){ */ + +/* int i; */ +/* const int numdof = NDOF1 *NUMVERTICES; */ +/* const int numdof2d = NDOF1 *NUMVERTICES2D; */ +/* bool isefficientlayer; */ +/* IssmDouble n,A,dt,init_thick; */ +/* IssmDouble rho_water,rho_ice; */ +/* IssmDouble gravity,latentheat,EPLgrad; */ +/* IssmDouble EPL_N,epl_conductivity; */ +/* IssmDouble activeEpl[numdof],thickness[numdof]; */ +/* IssmDouble eplhead[numdof], old_thickness[numdof]; */ +/* IssmDouble epl_slopeX[numdof],epl_slopeY[numdof]; */ +/* IssmDouble ice_thickness[numdof],bed[numdof]; */ +/* Penta *penta = NULL; */ +/* /\*If not on bed, return*\/ */ +/* if (!IsOnBed())return; */ + +/* /\*Get the flag to know if the efficient layer is present*\/ */ +/* this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); */ +/* this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); */ + +/* if(isefficientlayer){ */ +/* /\*For now, assuming just one way to compute EPL thickness*\/ */ +/* rho_water = matpar->GetRhoWater(); */ +/* rho_ice = matpar->GetRhoIce(); */ +/* gravity = matpar->GetG(); */ +/* latentheat = matpar->GetLatentHeat(); */ +/* epl_conductivity = matpar->GetEplConductivity(); */ +/* init_thick = matpar->GetEplInitialThickness(); */ +/* n = material->GetN(); */ +/* A = material->GetA(); */ - GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveNodeEnum); - GetInputListOnVertices(&eplhead[0],EplHeadEnum); - GetInputListOnVertices(&epl_slopeX[0],EplHeadSlopeXEnum); - GetInputListOnVertices(&epl_slopeY[0],EplHeadSlopeYEnum); - GetInputListOnVertices(&old_thickness[0],HydrologydcEplThicknessOldEnum); - GetInputListOnVertices(&ice_thickness[0],ThicknessEnum); - GetInputListOnVertices(&bed[0],BedEnum); +/* GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveNodeEnum); */ +/* GetInputListOnVertices(&eplhead[0],EplHeadEnum); */ +/* GetInputListOnVertices(&epl_slopeX[0],EplHeadSlopeXEnum); */ +/* GetInputListOnVertices(&epl_slopeY[0],EplHeadSlopeYEnum); */ +/* GetInputListOnVertices(&old_thickness[0],HydrologydcEplThicknessOldEnum); */ +/* GetInputListOnVertices(&ice_thickness[0],ThicknessEnum); */ +/* GetInputListOnVertices(&bed[0],BedEnum); */ - for(int i=0;iinputs->AddInput(new PentaInput(HydrologydcEplThicknessEnum,thickness,P1Enum)); +/* /\*And proceed to the real thing*\/ */ +/* thickness[i] = old_thickness[i]*(1+((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*pow(EPLgrad,2.0)-2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n))); */ +/* thickness[i+numdof2d]=thickness[i]; */ +/* } */ +/* } */ +/* penta=this; */ +/* for(;;){ */ + +/* /\*Add input to the element: *\/ */ +/* penta->inputs->AddInput(new PentaInput(HydrologydcEplThicknessEnum,thickness,P1Enum)); */ - /*Stop if we have reached the surface*/ - if (penta->IsOnSurface()) break; +/* /\*Stop if we have reached the surface*\/ */ +/* if (penta->IsOnSurface()) break; */ - /* get upper Penta*/ - penta=penta->GetUpperPenta(); _assert_(penta->Id()!=this->id); - } - } -} -/*}}}*/ +/* /\* get upper Penta*\/ */ +/* penta=penta->GetUpperPenta(); _assert_(penta->Id()!=this->id); */ +/* } */ +/* } */ +/* } */ +/* /\*}}}*\/ */ /*FUNCTION Penta::MigrateGroundingLine{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 17360) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 17361) @@ -244,5 +244,5 @@ void HydrologyEPLGetActive(Vector* active_vec); void HydrologyEPLGetMask(Vector* vec_mask); - void ComputeEPLThickness(void); + // void ComputeEPLThickness(void); void UpdateConstraintsExtrudeFromBase(void); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 17360) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 17361) @@ -153,5 +153,6 @@ void HydrologyEPLGetActive(Vector* active_vec){_error_("not implemented yet");}; void HydrologyEPLGetMask(Vector* vec_mask){_error_("not implemented yet");}; - void ComputeEPLThickness(void){_error_("not implemented yet");}; + // void ComputeEPLThickness(void){_error_("not implemented yet");}; + void GetSolutionFromInputsOneDof(Vector* solution,int enum_type){_error_("not implemented yet");}; void GetVectorFromInputs(Vector* vector, int name_enum){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17360) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17361) @@ -4546,75 +4546,75 @@ } /*}}}*/ -/*FUNCTION Tria::ComputeEPLThickness{{{*/ -void Tria::ComputeEPLThickness(void){ - - const int numdof = NDOF1 *NUMVERTICES; - bool isefficientlayer; - bool active_element; - IssmDouble n,A,dt,init_thick; - IssmDouble rho_water,rho_ice; - IssmDouble gravity,latentheat,EPLgrad2; - IssmDouble EPL_N,epl_conductivity; - IssmDouble thickness[numdof]; - IssmDouble eplhead[numdof]; - IssmDouble epl_slopeX[numdof],epl_slopeY[numdof]; - IssmDouble old_thickness[numdof]; - IssmDouble ice_thickness[numdof],bed[numdof]; - - Input* active_element_input=NULL; - - /*Get the flag to know if the efficient layer is present*/ - this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - - if(isefficientlayer){ - - /*For now, assuming just one way to compute EPL thickness*/ - rho_water = matpar->GetRhoWater(); - rho_ice = matpar->GetRhoIce(); - gravity = matpar->GetG(); - latentheat = matpar->GetLatentHeat(); - epl_conductivity = matpar->GetEplConductivity(); - init_thick = matpar->GetEplInitialThickness(); - n = material->GetN(); - A = material->GetAbar(); +//* *FUNCTION Tria::ComputeEPLThickness{{{*\/ */ +/* void Tria::ComputeEPLThickness(void){ */ + +/* const int numdof = NDOF1 *NUMVERTICES; */ +/* bool isefficientlayer; */ +/* bool active_element; */ +/* IssmDouble n,A,dt,init_thick; */ +/* IssmDouble rho_water,rho_ice; */ +/* IssmDouble gravity,latentheat,EPLgrad2; */ +/* IssmDouble EPL_N,epl_conductivity; */ +/* IssmDouble thickness[numdof]; */ +/* IssmDouble eplhead[numdof]; */ +/* IssmDouble epl_slopeX[numdof],epl_slopeY[numdof]; */ +/* IssmDouble old_thickness[numdof]; */ +/* IssmDouble ice_thickness[numdof],bed[numdof]; */ + +/* Input* active_element_input=NULL; */ + +/* /\*Get the flag to know if the efficient layer is present*\/ */ +/* this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); */ +/* this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); */ + +/* if(isefficientlayer){ */ + +/* /\*For now, assuming just one way to compute EPL thickness*\/ */ +/* rho_water = matpar->GetRhoWater(); */ +/* rho_ice = matpar->GetRhoIce(); */ +/* gravity = matpar->GetG(); */ +/* latentheat = matpar->GetLatentHeat(); */ +/* epl_conductivity = matpar->GetEplConductivity(); */ +/* init_thick = matpar->GetEplInitialThickness(); */ +/* n = material->GetN(); */ +/* A = material->GetAbar(); */ - active_element_input=inputs->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); - active_element_input->GetInputValue(&active_element); +/* active_element_input=inputs->GetInput(HydrologydcMaskEplactiveEltEnum); _assert_(active_element_input); */ +/* active_element_input->GetInputValue(&active_element); */ - GetInputListOnVertices(&eplhead[0],EplHeadEnum); - GetInputListOnVertices(&epl_slopeX[0],EplHeadSlopeXEnum); - GetInputListOnVertices(&epl_slopeY[0],EplHeadSlopeYEnum); - GetInputListOnVertices(&old_thickness[0],HydrologydcEplThicknessOldEnum); - GetInputListOnVertices(&ice_thickness[0],ThicknessEnum); - GetInputListOnVertices(&bed[0],BedEnum); +/* GetInputListOnVertices(&eplhead[0],EplHeadEnum); */ +/* GetInputListOnVertices(&epl_slopeX[0],EplHeadSlopeXEnum); */ +/* GetInputListOnVertices(&epl_slopeY[0],EplHeadSlopeYEnum); */ +/* GetInputListOnVertices(&old_thickness[0],HydrologydcEplThicknessOldEnum); */ +/* GetInputListOnVertices(&ice_thickness[0],ThicknessEnum); */ +/* GetInputListOnVertices(&bed[0],BedEnum); */ - if(!active_element){ - /*Keeping thickness to initial value if EPL is not active*/ - for(int i=0;i10.0*init_thick){ - thickness[i] = 10.0*init_thick; - } - } - } - this->inputs->AddInput(new TriaInput(HydrologydcEplThicknessEnum,thickness,P1Enum)); - } -} +/* /\*Take care of otherthikening*\/ */ +/* if(thickness[i]>10.0*init_thick){ */ +/* thickness[i] = 10.0*init_thick; */ +/* } */ +/* } */ +/* } */ +/* this->inputs->AddInput(new TriaInput(HydrologydcEplThicknessEnum,thickness,P1Enum)); */ +/* } */ +/* } */ /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 17360) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 17361) @@ -253,8 +253,9 @@ void CreateHydrologyWaterVelocityInput(void); + void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode); void HydrologyEPLGetActive(Vector* active_vec); void HydrologyEPLGetMask(Vector* vec_mask); - void ComputeEPLThickness(void); + // void ComputeEPLThickness(void); /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 17360) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 17361) @@ -1407,11 +1407,11 @@ } /*}}}*/ -void FemModel::HydrologyEPLThicknessx(void){ /*{{{*/ - - for (int i=0;iSize();i++){ - Element* element=dynamic_cast(elements->GetObjectByOffset(i)); - element->ComputeEPLThickness(); - } -} +/* void FemModel::HydrologyEPLThicknessx(void){ /\*{{{*\/ */ + +/* for (int i=0;iSize();i++){ */ +/* Element* element=dynamic_cast(elements->GetObjectByOffset(i)); */ +/* element->ComputeEPLThickness(); */ +/* } */ +/* } */ /*}}}*/ void FemModel::UpdateConstraintsL2ProjectionEPLx(void){ /*{{{*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 17360) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 17361) @@ -96,5 +96,5 @@ void UpdateConstraintsExtrudeFromTopx(); void HydrologyEPLupdateDomainx(void); - void HydrologyEPLThicknessx(void); + // void HydrologyEPLThicknessx(void); void UpdateConstraintsL2ProjectionEPLx(void); }; Index: /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 17360) +++ /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 17361) @@ -639,5 +639,5 @@ if(zigzag_counter>penalty_lock){ unstable=0; - active=1; + active=0; } } Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 17360) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 17361) @@ -33,4 +33,7 @@ Vector* df=NULL; + HydrologyDCInefficientAnalysis* inefanalysis = NULL; + HydrologyDCEfficientAnalysis* effanalysis = NULL; + bool sedconverged,eplconverged,hydroconverged; bool isefficientlayer; @@ -60,11 +63,10 @@ if(isefficientlayer) { + inefanalysis = new HydrologyDCInefficientAnalysis(); + effanalysis = new HydrologyDCEfficientAnalysis(); femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum); GetSolutionFromInputsx(&ug_epl,femmodel); - /*Initialize the element mask*/ - HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis(); - analysis->ElementizeEplMask(femmodel); - delete analysis; + inefanalysis->ElementizeEplMask(femmodel); } /*The real computation starts here, outermost loop is on the two layer system*/ @@ -115,4 +117,5 @@ if(num_unstable_constraints==0) sedconverged = true; if (sedcount>=hydro_maxiter){ + //sedconverged = true; _error_(" maximum number of Sediment iterations (" << hydro_maxiter << ") exceeded"); } @@ -126,7 +129,5 @@ /*Update EPL mask*/ if(isefficientlayer){ - HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis(); - analysis->ElementizeEplMask(femmodel); - delete analysis; + inefanalysis->ElementizeEplMask(femmodel); } sedconverged=false; @@ -180,12 +181,10 @@ solutionsequence_linear(femmodel); femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum); - femmodel->HydrologyEPLThicknessx(); - + + effanalysis->ComputeEPLThickness(femmodel); + //updating mask after the computation of the epl thickness (Allow to close too thin EPL) femmodel->HydrologyEPLupdateDomainx(); - - HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis(); - analysis->ElementizeEplMask(femmodel); - delete analysis; + inefanalysis->ElementizeEplMask(femmodel); /* }}} */ @@ -281,4 +280,5 @@ delete uf_epl; delete uf_epl_sub_iter; - + delete inefanalysis; + delete effanalysis; }