Index: ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 16599) +++ ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 16600) @@ -13,6 +13,7 @@ int hydrology_model; int sedimentlimit_flag; int transfer_flag; + int penalty_lock; bool isefficientlayer; IssmDouble sedimentlimit; IssmDouble penalty_factor; @@ -30,6 +31,7 @@ iomodel->FetchData(&transfer_flag,HydrologydcTransferFlagEnum); iomodel->FetchData(&penalty_factor,HydrologydcPenaltyFactorEnum); iomodel->FetchData(&rel_tol,HydrologydcRelTolEnum); + iomodel->FetchData(&penalty_lock,HydrologydcPenaltyLockEnum); if(sedimentlimit_flag==1){ iomodel->FetchData(&sedimentlimit,HydrologydcSedimentlimitEnum); @@ -47,6 +49,7 @@ parameters->AddObject(new IntParam(HydrologydcSedimentlimitFlagEnum,sedimentlimit_flag)); parameters->AddObject(new IntParam(HydrologydcTransferFlagEnum,transfer_flag)); parameters->AddObject(new DoubleParam(HydrologydcRelTolEnum,rel_tol)); + parameters->AddObject(new IntParam(HydrologydcPenaltyLockEnum,penalty_lock)); }/*}}}*/ void HydrologyDCInefficientAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16599) +++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16600) @@ -120,6 +120,7 @@ HydrologydcTransferFlagEnum, HydrologydcLeakageFactorEnum, HydrologydcPenaltyFactorEnum, + HydrologydcPenaltyLockEnum, HydrologyLayerEnum, HydrologySedimentEnum, HydrologyEfficientEnum, Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16599) +++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16600) @@ -128,6 +128,7 @@ case HydrologydcTransferFlagEnum : return "HydrologydcTransferFlag"; case HydrologydcLeakageFactorEnum : return "HydrologydcLeakageFactor"; case HydrologydcPenaltyFactorEnum : return "HydrologydcPenaltyFactor"; + case HydrologydcPenaltyLockEnum : return "HydrologydcPenaltyLock"; case HydrologyLayerEnum : return "HydrologyLayer"; case HydrologySedimentEnum : return "HydrologySediment"; case HydrologyEfficientEnum : return "HydrologyEfficient"; Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16599) +++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16600) @@ -128,6 +128,7 @@ else if (strcmp(name,"HydrologydcTransferFlag")==0) return HydrologydcTransferFlagEnum; else if (strcmp(name,"HydrologydcLeakageFactor")==0) return HydrologydcLeakageFactorEnum; else if (strcmp(name,"HydrologydcPenaltyFactor")==0) return HydrologydcPenaltyFactorEnum; + else if (strcmp(name,"HydrologydcPenaltyLock")==0) return HydrologydcPenaltyLockEnum; else if (strcmp(name,"HydrologyLayer")==0) return HydrologyLayerEnum; else if (strcmp(name,"HydrologySediment")==0) return HydrologySedimentEnum; else if (strcmp(name,"HydrologyEfficient")==0) return HydrologyEfficientEnum; @@ -135,11 +136,11 @@ else if (strcmp(name,"WaterTransfer")==0) return WaterTransferEnum; else if (strcmp(name,"IndependentObject")==0) return IndependentObjectEnum; else if (strcmp(name,"InversionControlParameters")==0) return InversionControlParametersEnum; - else if (strcmp(name,"InversionCostFunctionThreshold")==0) return InversionCostFunctionThresholdEnum; else stage=2; } if(stage==2){ - if (strcmp(name,"InversionCostFunctionsCoefficients")==0) return InversionCostFunctionsCoefficientsEnum; + if (strcmp(name,"InversionCostFunctionThreshold")==0) return InversionCostFunctionThresholdEnum; + else if (strcmp(name,"InversionCostFunctionsCoefficients")==0) return InversionCostFunctionsCoefficientsEnum; else if (strcmp(name,"InversionCostFunctions")==0) return InversionCostFunctionsEnum; else if (strcmp(name,"InversionGradientScaling")==0) return InversionGradientScalingEnum; else if (strcmp(name,"InversionIscontrol")==0) return InversionIscontrolEnum; @@ -258,11 +259,11 @@ else if (strcmp(name,"SteadystateReltol")==0) return SteadystateReltolEnum; else if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum; else if (strcmp(name,"Surface")==0) return SurfaceEnum; - else if (strcmp(name,"SurfaceforcingsPrecipitation")==0) return SurfaceforcingsPrecipitationEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"SurfaceforcingsMassBalance")==0) return SurfaceforcingsMassBalanceEnum; + if (strcmp(name,"SurfaceforcingsPrecipitation")==0) return SurfaceforcingsPrecipitationEnum; + else if (strcmp(name,"SurfaceforcingsMassBalance")==0) return SurfaceforcingsMassBalanceEnum; else if (strcmp(name,"SurfaceforcingsIspdd")==0) return SurfaceforcingsIspddEnum; else if (strcmp(name,"SurfaceforcingsDesfac")==0) return SurfaceforcingsDesfacEnum; else if (strcmp(name,"SurfaceforcingsS0p")==0) return SurfaceforcingsS0pEnum; @@ -381,11 +382,11 @@ else if (strcmp(name,"Element")==0) return ElementEnum; else if (strcmp(name,"ExternalResult")==0) return ExternalResultEnum; else if (strcmp(name,"FileParam")==0) return FileParamEnum; - else if (strcmp(name,"Input")==0) return InputEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"IntInput")==0) return IntInputEnum; + if (strcmp(name,"Input")==0) return InputEnum; + else if (strcmp(name,"IntInput")==0) return IntInputEnum; else if (strcmp(name,"InputToExtrude")==0) return InputToExtrudeEnum; else if (strcmp(name,"InputToL2Project")==0) return InputToL2ProjectEnum; else if (strcmp(name,"IntParam")==0) return IntParamEnum; @@ -504,11 +505,11 @@ else if (strcmp(name,"QmuTemperature")==0) return QmuTemperatureEnum; else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum; else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum; - else if (strcmp(name,"StressTensor")==0) return StressTensorEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"StressTensorxx")==0) return StressTensorxxEnum; + if (strcmp(name,"StressTensor")==0) return StressTensorEnum; + else if (strcmp(name,"StressTensorxx")==0) return StressTensorxxEnum; else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum; else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum; else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum; Index: ../trunk-jpl/src/c/classes/Loads/Pengrid.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 16599) +++ ../trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 16600) @@ -639,8 +639,9 @@ void Pengrid::ConstraintActivateHydrologyDCInefficient(int* punstable){ // The penalty is stable if it doesn't change during two consecutive iterations. - int unstable = 0; + int unstable=0; int new_active; + int penalty_lock; IssmDouble pressure; IssmDouble h; IssmDouble h_max; @@ -655,16 +656,28 @@ /*Get sediment water head h*/ element->GetInputValue(&h,node,SedimentHeadEnum); element->GetHydrologyDCInefficientHmax(&h_max,node); + parameters->FindParam(&penalty_lock,HydrologydcPenaltyLockEnum); + if (h>h_max) new_active=1; else new_active=0; - if(this->active==new_active) - unstable=0; - else - unstable=1; + if(this->active==new_active){ + unstable=0; + } + else{ + unstable=1; + if(penalty_lock)zigzag_counter++; + } + /*If penalty keeps zigzagging more than penalty_lock times: */ + if(penalty_lock){ + if(zigzag_counter>penalty_lock){ + unstable=0; + active=1; + } + } /*Set penalty flag*/ this->active=new_active; Index: ../trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 16599) +++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 16600) @@ -1411,4 +1411,12 @@ delete transferg; } /*}}}*/ +void FemModel::HydrologyEPLThicknessx(void){ /*{{{*/ + + for (int i=0;iSize();i++){ + Element* element=dynamic_cast(elements->GetObjectByOffset(i)); + element->ComputeEPLThickness(); + } +} +/*}}}*/ #endif Index: ../trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Element.h (revision 16599) +++ ../trunk-jpl/src/c/classes/Elements/Element.h (revision 16600) @@ -135,6 +135,7 @@ virtual void GetHydrologyTransfer(Vector* transfer)=0; virtual void HydrologyEPLGetMask(Vector* mask)=0; virtual void HydrologyEPLGetActive(Vector* active)=0; + virtual void ComputeEPLThickness(void)=0; #endif #ifdef _HAVE_GROUNDINGLINE_ Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16599) +++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16600) @@ -7223,6 +7223,12 @@ } } /*}}}*/ +/*FUNCTION Tria::ComputeEPLThickness{{{*/ +void Tria::ComputeEPLThickness(void){ +} +/*}}}*/ + + #endif #ifdef _HAVE_MASSTRANSPORT_ Index: ../trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 16599) +++ ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 16600) @@ -302,6 +302,7 @@ void GetHydrologyTransfer(Vector* transfer); void HydrologyEPLGetActive(Vector* active_vec); void HydrologyEPLGetMask(Vector* vec_mask); + void ComputeEPLThickness(void); bool AllActive(void); bool AnyActive(void); #endif Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16599) +++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16600) @@ -10844,6 +10844,17 @@ } /*}}}*/ +/*FUNCTION Penta::ComputeEPLThickness{{{*/ +void Penta::ComputeEPLThickness(void){ + + if (!IsOnBed())return; + + Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. + tria->ComputeEPLThickness(); + delete tria->material; delete tria; + +} +/*}}}*/ /*FUNCTION Penta::InputUpdateFromSolutionHydrologyDCInefficient{{{*/ void Penta::InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution){ Index: ../trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 16599) +++ ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 16600) @@ -322,6 +322,7 @@ void HydrologyEPLGetActive(Vector* active_vec); void HydrologyEPLGetMask(Vector* vec_mask); void InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution); + void ComputeEPLThickness(void); #endif void UpdateConstraintsExtrudeFromBase(void){_error_("not implemented yet");}; Index: ../trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 16599) +++ ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 16600) @@ -100,6 +100,7 @@ void GetHydrologyTransfer(Vector* transfer){_error_("not implemented yet");}; 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");}; #endif void GetSolutionFromInputs(Vector* solution){_error_("not implemented yet");}; void GetVectorFromInputs(Vector* vector, int name_enum){_error_("not implemented yet");}; Index: ../trunk-jpl/src/c/classes/FemModel.h =================================================================== --- ../trunk-jpl/src/c/classes/FemModel.h (revision 16599) +++ ../trunk-jpl/src/c/classes/FemModel.h (revision 16600) @@ -101,6 +101,7 @@ #ifdef _HAVE_HYDROLOGY_ void HydrologyTransferx(void); void HydrologyEPLupdateDomainx(void); + void HydrologyEPLThicknessx(void); #endif }; Index: ../trunk-jpl/src/m/classes/hydrologydc.m =================================================================== --- ../trunk-jpl/src/m/classes/hydrologydc.m (revision 16599) +++ ../trunk-jpl/src/m/classes/hydrologydc.m (revision 16600) @@ -8,6 +8,7 @@ water_compressibility = 0; isefficientlayer = 0; penalty_factor = 0; + penalty_lock = 0; rel_tol = 0; sedimentlimit_flag = 0; sedimentlimit = 0; @@ -106,6 +107,7 @@ fielddisplay(obj,'water_compressibility','compressibility of water [Pa^-1]'); fielddisplay(obj,'isefficientlayer','do we use an efficient drainage system [1: true; 0: false]'); fielddisplay(obj,'penalty_factor','exponent of the value used in the penalisation method [dimensionless]'); + fielddisplay(obj,'penalty_lock','stabilize unstable constraints that keep zigzagging after n iteration (default is 0, no stabilization)'); fielddisplay(obj,'rel_tol','tolerance of the nonlinear iteration for the transfer between layers [dimensionless]'); fielddisplay(obj,'sedimentlimit_flag','what kind of upper limit is applied for the inefficient layer'); disp(sprintf('%55s 0: no limit',' ')); @@ -146,6 +148,7 @@ WriteData(fid,'object',obj,'fieldname','water_compressibility','format','Double'); WriteData(fid,'object',obj,'fieldname','isefficientlayer','format','Boolean'); WriteData(fid,'object',obj,'fieldname','penalty_factor','format','Double'); + WriteData(fid,'object',obj,'fieldname','penalty_lock','format','Integer'); WriteData(fid,'object',obj,'fieldname','rel_tol','format','Double'); WriteData(fid,'object',obj,'fieldname','sedimentlimit_flag','format','Integer'); WriteData(fid,'object',obj,'fieldname','transfer_flag','format','Integer'); Index: ../trunk-jpl/src/m/enum/HydrologydcPenaltyLockEnum.m =================================================================== --- ../trunk-jpl/src/m/enum/HydrologydcPenaltyLockEnum.m (revision 0) +++ ../trunk-jpl/src/m/enum/HydrologydcPenaltyLockEnum.m (revision 16600) @@ -0,0 +1,11 @@ +function macro=HydrologydcPenaltyLockEnum() +%HYDROLOGYDCPENALTYLOCKENUM - Enum of HydrologydcPenaltyLock +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=HydrologydcPenaltyLockEnum() + +macro=StringToEnum('HydrologydcPenaltyLock'); Index: ../trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- ../trunk-jpl/src/m/enum/EnumDefinitions.py (revision 16599) +++ ../trunk-jpl/src/m/enum/EnumDefinitions.py (revision 16600) @@ -120,6 +120,7 @@ def HydrologydcTransferFlagEnum(): return StringToEnum("HydrologydcTransferFlag")[0] def HydrologydcLeakageFactorEnum(): return StringToEnum("HydrologydcLeakageFactor")[0] def HydrologydcPenaltyFactorEnum(): return StringToEnum("HydrologydcPenaltyFactor")[0] +def HydrologydcPenaltyLockEnum(): return StringToEnum("HydrologydcPenaltyLock")[0] def HydrologyLayerEnum(): return StringToEnum("HydrologyLayer")[0] def HydrologySedimentEnum(): return StringToEnum("HydrologySediment")[0] def HydrologyEfficientEnum(): return StringToEnum("HydrologyEfficient")[0]