Changeset 16600

Timestamp:

11/01/13 10:09:33 (11 years ago)

Author:

bdef

Message:

Adding a ZigZag lock for hydrology and framing the variable thickness of the EPL

Location:

issm/trunk-jpl/src

Files:

• c/analyses/HydrologyDCInefficientAnalysis.cpp (modified) (3 diffs)
• c/classes/Elements/Element.h (modified) (1 diff)
• c/classes/Elements/Penta.cpp (modified) (1 diff)
• c/classes/Elements/Penta.h (modified) (1 diff)
• c/classes/Elements/Seg.h (modified) (1 diff)
• c/classes/Elements/Tria.cpp (modified) (1 diff)
• c/classes/Elements/Tria.h (modified) (1 diff)
• c/classes/FemModel.cpp (modified) (1 diff)
• c/classes/FemModel.h (modified) (1 diff)
• c/classes/Loads/Pengrid.cpp (modified) (3 diffs)
• c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• c/shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• c/shared/Enum/StringToEnumx.cpp (modified) (5 diffs)
• m/classes/hydrologydc.m (modified) (3 diffs)
• m/enum/EnumDefinitions.py (modified) (1 diff)
• m/enum/HydrologydcPenaltyLockEnum.m (added)

issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp

@@ -14 +14 @@
         int         sedimentlimit_flag;
         int         transfer_flag;
+        int         penalty_lock;
         bool        isefficientlayer;
         IssmDouble  sedimentlimit;
@@ -31 +32 @@
         iomodel->FetchData(&penalty_factor,HydrologydcPenaltyFactorEnum);
         iomodel->FetchData(&rel_tol,HydrologydcRelTolEnum);
+        iomodel->FetchData(&penalty_lock,HydrologydcPenaltyLockEnum);
 
         if(sedimentlimit_flag==1){
@@ -48 +50 @@
         parameters->AddObject(new IntParam(HydrologydcTransferFlagEnum,transfer_flag));
         parameters->AddObject(new DoubleParam(HydrologydcRelTolEnum,rel_tol));
+        parameters->AddObject(new IntParam(HydrologydcPenaltyLockEnum,penalty_lock));
 
 }/*}}}*/

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

@@ -136 +136 @@
                 virtual void HydrologyEPLGetMask(Vector<IssmDouble>* mask)=0;
                 virtual void HydrologyEPLGetActive(Vector<IssmDouble>* active)=0;
+                virtual void ComputeEPLThickness(void)=0;
                 #endif
 

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

@@ -10845 +10845 @@
 }
 /*}}}*/
+/*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){

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

@@ -323 +323 @@
                 void    HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask);
                 void    InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution);
+                void    ComputeEPLThickness(void);
                 #endif
 

issm/trunk-jpl/src/c/classes/Elements/Seg.h

@@ -101 +101 @@
                 void    HydrologyEPLGetActive(Vector<IssmDouble>* active_vec){_error_("not implemented yet");};
                 void    HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask){_error_("not implemented yet");};
+                void    ComputeEPLThickness(void){_error_("not implemented yet");};
                 #endif
                 void        GetSolutionFromInputs(Vector<IssmDouble>* solution){_error_("not implemented yet");};

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

@@ -7224 +7224 @@
 }
 /*}}}*/
+/*FUNCTION Tria::ComputeEPLThickness{{{*/
+void  Tria::ComputeEPLThickness(void){
+}
+/*}}}*/
+
+
 #endif
 

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

@@ -303 +303 @@
                 void    HydrologyEPLGetActive(Vector<IssmDouble>* active_vec);
                 void    HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask);
+                void    ComputeEPLThickness(void);
                 bool    AllActive(void);
                 bool    AnyActive(void);

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

@@ -1412 +1412 @@
 }
 /*}}}*/
+void FemModel::HydrologyEPLThicknessx(void){ /*{{{*/
+
+        for (int i=0;i<elements->Size();i++){
+                Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
+                element->ComputeEPLThickness();
+        }
+}
+/*}}}*/
 #endif

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

@@ -102 +102 @@
                 void HydrologyTransferx(void);
                 void HydrologyEPLupdateDomainx(void);
+                void HydrologyEPLThicknessx(void);
                 #endif
 };

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

@@ -640 +640 @@
 
         //   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;
@@ -656 +657 @@
         element->GetInputValue(&h,node,SedimentHeadEnum);
         element->GetHydrologyDCInefficientHmax(&h_max,node);
+        parameters->FindParam(&penalty_lock,HydrologydcPenaltyLockEnum);
+
         if (h>h_max)
          new_active=1;
@@ -661 +664 @@
          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;

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

@@ -121 +121 @@
         HydrologydcLeakageFactorEnum,
         HydrologydcPenaltyFactorEnum,
+        HydrologydcPenaltyLockEnum,
         HydrologyLayerEnum,
         HydrologySedimentEnum,

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

@@ -129 +129 @@
                 case HydrologydcLeakageFactorEnum : return "HydrologydcLeakageFactor";
                 case HydrologydcPenaltyFactorEnum : return "HydrologydcPenaltyFactor";
+                case HydrologydcPenaltyLockEnum : return "HydrologydcPenaltyLock";
                 case HydrologyLayerEnum : return "HydrologyLayer";
                 case HydrologySedimentEnum : return "HydrologySediment";

issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

@@ -129 +129 @@
               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;
@@ -136 +137 @@
               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;
@@ -259 +260 @@
               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;
@@ -382 +383 @@
               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;
@@ -505 +506 @@
               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;

issm/trunk-jpl/src/m/classes/hydrologydc.m

@@ -9 +9 @@
                 isefficientlayer         = 0;
                 penalty_factor           = 0;
+                penalty_lock             = 0;
                 rel_tol                  = 0;
                 sedimentlimit_flag       = 0;
@@ -107 +108 @@
                         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');
@@ -147 +149 @@
                         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');

issm/trunk-jpl/src/m/enum/EnumDefinitions.py

@@ -121 +121 @@
 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]