Changeset 21964

Timestamp:

08/17/17 05:27:12 (8 years ago)

Author:

bdef

Message:

NEW:adding a bypass for unconfined scheme (transitory for nondev use

Location:

issm/trunk-jpl/src/c

Files:

• analyses/HydrologyDCInefficientAnalysis.cpp (modified) (8 diffs)
• shared/Enum/EnumDefinitions.h (modified) (1 diff)
• shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• shared/Enum/StringToEnumx.cpp (modified) (8 diffs)

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

@@ -16 +16 @@
         int         sedimentlimit_flag;
         int         transfer_flag;
+        int         unconfined_flag;
         int         penalty_lock;
         int         hydro_maxiter;
@@ -32 +33 @@
         iomodel->FetchData(&sedimentlimit_flag, "md.hydrology.sedimentlimit_flag" );
         iomodel->FetchData(&transfer_flag,      "md.hydrology.transfer_flag" );
+        iomodel->FetchData(&unconfined_flag,      "md.hydrology.unconfined_flag" );
         iomodel->FetchData(&penalty_factor,     "md.hydrology.penalty_factor" );
         iomodel->FetchData(&isefficientlayer,   "md.hydrology.isefficientlayer");
@@ -41 +43 @@
         parameters->AddObject(new IntParam(HydrologydcSedimentlimitFlagEnum,sedimentlimit_flag));
         parameters->AddObject(new IntParam(HydrologydcTransferFlagEnum,transfer_flag));
+        parameters->AddObject(new IntParam(HydrologydcUnconfinedFlagEnum,unconfined_flag));
         parameters->AddObject(new IntParam(HydrologydcPenaltyLockEnum,penalty_lock));
         parameters->AddObject(new IntParam(HydrologydcMaxIterEnum,hydro_maxiter));
@@ -79 +82 @@
                 }
         }
-
         iomodel->FetchDataToInput(elements,"md.geometry.thickness",ThicknessEnum);
         iomodel->FetchDataToInput(elements,"md.geometry.base",BaseEnum);
@@ -220 +222 @@
         /* Start  looping on the number of gaussian points: */
         Gauss* gauss=basalelement->NewGauss(2);
+
         for(int ig=gauss -> begin();ig<gauss->end();ig++){
                 gauss          -> GaussPoint(ig);
@@ -317 +320 @@
         Input* base_input                 = basalelement->GetInput(BaseEnum);
         Input* water_input        = basalelement->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(water_input);
+
         if(dt!= 0.){
                 old_wh_input = basalelement->GetInput(SedimentHeadOldEnum);                  _assert_(old_wh_input);
@@ -513 +517 @@
 
 IssmDouble HydrologyDCInefficientAnalysis::SedimentStoring(Element* element,Gauss* gauss,Input* sed_head_input, Input* base_input){/*{{{*/
+        int unconf_scheme;
+        IssmDouble expfac; 
         IssmDouble sediment_storing;
         IssmDouble storing;
@@ -522 +528 @@
         IssmDouble sediment_compressibility = element->GetMaterialParameter(HydrologydcSedimentCompressibilityEnum);
         IssmDouble water_compressibility    = element->GetMaterialParameter(HydrologydcWaterCompressibilityEnum);
-        base_input->GetInputValue(&base_elev,gauss);
-        sed_head_input->GetInputValue(&prestep_head,gauss);
-        water_sheet=max(0.0,(prestep_head-(base_elev-sediment_thickness)));
-
-        storing=rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));
+        element->FindParam(&unconf_scheme,HydrologydcUnconfinedFlagEnum);
+        switch(unconf_scheme){
+        case 0:
+                sediment_storing=rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));
+                break;
+        case 1:
+                base_input->GetInputValue(&base_elev,gauss);
+                sed_head_input->GetInputValue(&prestep_head,gauss);
+                water_sheet=max(0.0,(prestep_head-(base_elev-sediment_thickness)));
+
+                storing=rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));
         
-        //Heavyside approximation (1/(1+exp(-2kx))) with k=25 centering at thickness minus 1%
-        sediment_storing=(sediment_porosity*exp(-20.*(water_sheet-0.99*sediment_thickness))+storing)/(1+exp(-20.*(water_sheet-0.99*sediment_thickness)));
-
+                //      using logistic function for heavyside approximation
+                expfac=10.;
+                sediment_storing=sediment_porosity+(storing-sediment_porosity)/(1+exp(-2*expfac*(water_sheet-0.99*sediment_thickness)));
+                break;
+        default:
+                _error_("UnconfinedFlag is 0 or 1");
+        }
         return sediment_storing;
 }/*}}}*/
 
 IssmDouble HydrologyDCInefficientAnalysis::SedimentTransmitivity(Element* element,Gauss* gauss,Input* sed_head_input, Input* base_input,Input* SedTrans_input){/*{{{*/
+        int unconf_scheme;
+        IssmDouble ratio,expfac;
         IssmDouble sediment_transmitivity;
         IssmDouble FullLayer_transmitivity;
         IssmDouble base_elev,prestep_head,water_sheet;
         IssmDouble sediment_thickness       = element->GetMaterialParameter(HydrologydcSedimentThicknessEnum);
-        base_input->GetInputValue(&base_elev,gauss);
-        sed_head_input->GetInputValue(&prestep_head,gauss);
+
+        element->FindParam(&unconf_scheme,HydrologydcUnconfinedFlagEnum);
         SedTrans_input->GetInputValue(&FullLayer_transmitivity,gauss);
-        water_sheet=max(0.0,(prestep_head-(base_elev-sediment_thickness)));
-
-        if (water_sheet<=sediment_thickness){
-                sediment_transmitivity=FullLayer_transmitivity*water_sheet/sediment_thickness;
-        }
-        else{
+        switch(unconf_scheme){
+        case 0:
                 sediment_transmitivity=FullLayer_transmitivity;
+                break;
+        case 1:
+                base_input->GetInputValue(&base_elev,gauss);
+                sed_head_input->GetInputValue(&prestep_head,gauss);
+                water_sheet=max(0.0,(prestep_head-(base_elev-sediment_thickness)));
+
+                if (water_sheet<=sediment_thickness){
+                        sediment_transmitivity=FullLayer_transmitivity*water_sheet/sediment_thickness;
+                }
+                else{
+                        sediment_transmitivity=FullLayer_transmitivity;
+                }
+                break;
+        default:
+                _error_("UnconfinedFlag is 0 or 1");
         }
         return sediment_transmitivity;

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

@@ -162 +162 @@
         HydrologydcSedimentlimitEnum,
         HydrologydcTransferFlagEnum,
+        HydrologydcUnconfinedFlagEnum,
         HydrologydcLeakageFactorEnum,
         HydrologydcPenaltyFactorEnum,

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

@@ -168 +168 @@
                 case HydrologydcSedimentlimitEnum : return "HydrologydcSedimentlimit";
                 case HydrologydcTransferFlagEnum : return "HydrologydcTransferFlag";
+                case HydrologydcUnconfinedFlagEnum : return "HydrologydcUnconfinedFlag";
                 case HydrologydcLeakageFactorEnum : return "HydrologydcLeakageFactor";
                 case HydrologydcPenaltyFactorEnum : return "HydrologydcPenaltyFactor";

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

@@ -171 +171 @@
               else if (strcmp(name,"HydrologydcSedimentlimit")==0) return HydrologydcSedimentlimitEnum;
               else if (strcmp(name,"HydrologydcTransferFlag")==0) return HydrologydcTransferFlagEnum;
+              else if (strcmp(name,"HydrologydcUnconfinedFlag")==0) return HydrologydcUnconfinedFlagEnum;
               else if (strcmp(name,"HydrologydcLeakageFactor")==0) return HydrologydcLeakageFactorEnum;
               else if (strcmp(name,"HydrologydcPenaltyFactor")==0) return HydrologydcPenaltyFactorEnum;
@@ -259 +260 @@
               else if (strcmp(name,"Damage")==0) return DamageEnum;
               else if (strcmp(name,"NewDamage")==0) return NewDamageEnum;
-              else if (strcmp(name,"StressIntensityFactor")==0) return StressIntensityFactorEnum;
          else stage=3;
    }
    if(stage==3){
-              if (strcmp(name,"CalvingLaw")==0) return CalvingLawEnum;
+              if (strcmp(name,"StressIntensityFactor")==0) return StressIntensityFactorEnum;
+              else if (strcmp(name,"CalvingLaw")==0) return CalvingLawEnum;
               else if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum;
               else if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum;
@@ -382 +383 @@
               else if (strcmp(name,"BalancethicknessDiffusionCoefficient")==0) return BalancethicknessDiffusionCoefficientEnum;
               else if (strcmp(name,"BalancethicknessOmega")==0) return BalancethicknessOmegaEnum;
-              else if (strcmp(name,"BalancethicknessD0")==0) return BalancethicknessD0Enum;
          else stage=4;
    }
    if(stage==4){
-              if (strcmp(name,"Smb")==0) return SmbEnum;
+              if (strcmp(name,"BalancethicknessD0")==0) return BalancethicknessD0Enum;
+              else if (strcmp(name,"Smb")==0) return SmbEnum;
               else if (strcmp(name,"SmbAnalysis")==0) return SmbAnalysisEnum;
               else if (strcmp(name,"SmbSolution")==0) return SmbSolutionEnum;
@@ -505 +506 @@
               else if (strcmp(name,"Internal")==0) return InternalEnum;
               else if (strcmp(name,"MassFlux")==0) return MassFluxEnum;
-              else if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum;
          else stage=5;
    }
    if(stage==5){
-              if (strcmp(name,"Misfit")==0) return MisfitEnum;
+              if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum;
+              else if (strcmp(name,"Misfit")==0) return MisfitEnum;
               else if (strcmp(name,"Pressure")==0) return PressureEnum;
               else if (strcmp(name,"PressurePicard")==0) return PressurePicardEnum;
@@ -628 +629 @@
               else if (strcmp(name,"Outputdefinition20")==0) return Outputdefinition20Enum;
               else if (strcmp(name,"Outputdefinition21")==0) return Outputdefinition21Enum;
-              else if (strcmp(name,"Outputdefinition22")==0) return Outputdefinition22Enum;
          else stage=6;
    }
    if(stage==6){
-              if (strcmp(name,"Outputdefinition23")==0) return Outputdefinition23Enum;
+              if (strcmp(name,"Outputdefinition22")==0) return Outputdefinition22Enum;
+              else if (strcmp(name,"Outputdefinition23")==0) return Outputdefinition23Enum;
               else if (strcmp(name,"Outputdefinition24")==0) return Outputdefinition24Enum;
               else if (strcmp(name,"Outputdefinition25")==0) return Outputdefinition25Enum;
@@ -751 +752 @@
               else if (strcmp(name,"OutputFilePointer")==0) return OutputFilePointerEnum;
               else if (strcmp(name,"ToolkitsFileName")==0) return ToolkitsFileNameEnum;
-              else if (strcmp(name,"RootPath")==0) return RootPathEnum;
          else stage=7;
    }
    if(stage==7){
-              if (strcmp(name,"OutputFileName")==0) return OutputFileNameEnum;
+              if (strcmp(name,"RootPath")==0) return RootPathEnum;
+              else if (strcmp(name,"OutputFileName")==0) return OutputFileNameEnum;
               else if (strcmp(name,"InputFileName")==0) return InputFileNameEnum;
               else if (strcmp(name,"LockFileName")==0) return LockFileNameEnum;
@@ -874 +875 @@
               else if (strcmp(name,"ElementHook")==0) return ElementHookEnum;
               else if (strcmp(name,"Hook")==0) return HookEnum;
-              else if (strcmp(name,"ExternalResult")==0) return ExternalResultEnum;
          else stage=8;
    }
    if(stage==8){
-              if (strcmp(name,"FileParam")==0) return FileParamEnum;
+              if (strcmp(name,"ExternalResult")==0) return ExternalResultEnum;
+              else if (strcmp(name,"FileParam")==0) return FileParamEnum;
               else if (strcmp(name,"Input")==0) return InputEnum;
               else if (strcmp(name,"IntInput")==0) return IntInputEnum;
@@ -997 +998 @@
               else if (strcmp(name,"MinVel")==0) return MinVelEnum;
               else if (strcmp(name,"MaxVel")==0) return MaxVelEnum;
-              else if (strcmp(name,"MinVx")==0) return MinVxEnum;
          else stage=9;
    }
    if(stage==9){
-              if (strcmp(name,"MaxVx")==0) return MaxVxEnum;
+              if (strcmp(name,"MinVx")==0) return MinVxEnum;
+              else if (strcmp(name,"MaxVx")==0) return MaxVxEnum;
               else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum;
               else if (strcmp(name,"MinVy")==0) return MinVyEnum;