Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 21963) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 21964) @@ -16,4 +16,5 @@ int sedimentlimit_flag; int transfer_flag; + int unconfined_flag; int penalty_lock; int hydro_maxiter; @@ -32,4 +33,5 @@ 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,4 +43,5 @@ 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,5 +82,4 @@ } } - iomodel->FetchDataToInput(elements,"md.geometry.thickness",ThicknessEnum); iomodel->FetchDataToInput(elements,"md.geometry.base",BaseEnum); @@ -220,4 +222,5 @@ /* Start looping on the number of gaussian points: */ Gauss* gauss=basalelement->NewGauss(2); + for(int ig=gauss -> begin();igend();ig++){ gauss -> GaussPoint(ig); @@ -317,4 +320,5 @@ 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,4 +517,6 @@ 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,31 +528,54 @@ 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; Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21963) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21964) @@ -162,4 +162,5 @@ HydrologydcSedimentlimitEnum, HydrologydcTransferFlagEnum, + HydrologydcUnconfinedFlagEnum, HydrologydcLeakageFactorEnum, HydrologydcPenaltyFactorEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21963) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21964) @@ -168,4 +168,5 @@ case HydrologydcSedimentlimitEnum : return "HydrologydcSedimentlimit"; case HydrologydcTransferFlagEnum : return "HydrologydcTransferFlag"; + case HydrologydcUnconfinedFlagEnum : return "HydrologydcUnconfinedFlag"; case HydrologydcLeakageFactorEnum : return "HydrologydcLeakageFactor"; case HydrologydcPenaltyFactorEnum : return "HydrologydcPenaltyFactor"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21963) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21964) @@ -171,4 +171,5 @@ 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,9 +260,9 @@ 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,9 +383,9 @@ 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,9 +506,9 @@ 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,9 +629,9 @@ 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,9 +752,9 @@ 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,9 +875,9 @@ 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,9 +998,9 @@ 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;