Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 24139) +++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 24140) @@ -275,4 +275,6 @@ parameters->AddObject(iomodel->CopyConstantObject("md.thermal.isenthalpy",ThermalIsenthalpyEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.thermal.isdynamicbasalspc",ThermalIsdynamicbasalspcEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.thermal.isdrainicecolumn",ThermalIsdrainicecolumnEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.thermal.watercolumn_upperlimit",ThermalWatercolumnUpperlimitEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.friction.law",FrictionLawEnum)); @@ -424,4 +426,6 @@ element->GetInputListOnNodes(basalmeltingrates,BasalforcingsGroundediceMeltingRateEnum); + IssmDouble watercolumnupperlimit = element->FindParam(ThermalWatercolumnUpperlimitEnum); + Gauss* gauss=element->NewGauss(); for(is=0;iswatercolumnupperlimit) watercolumns[nodedown]=watercolumnupperlimit; } else{ @@ -1168,5 +1173,9 @@ for(k=0;kThermalToEnthalpy(&enthalpies[k], temperatures[k], waterfractions[k], pressures[k]); } @@ -1187,10 +1196,8 @@ IssmDouble thermalconductivity = element->FindParam(MaterialsThermalconductivityEnum); - if(enthalpy < PureIceEnthalpy(element,pressure)){ + if(enthalpy < PureIceEnthalpy(element,pressure)) return thermalconductivity/heatcapacity; - } - else{ + else return temperateiceconductivity/heatcapacity; - } }/*}}}*/ IssmDouble EnthalpyAnalysis::EnthalpyDiffusionParameterVolume(Element* element,int enthalpy_enum){/*{{{*/ @@ -1678,10 +1685,11 @@ /*Intermediaries*/ bool computebasalmeltingrates=true; - bool drainicecolumn=true; + bool isdrainicecolumn; IssmDouble dt; femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - - if(drainicecolumn && (dt>0.)) DrainWaterfraction(femmodel); + femmodel->parameters->FindParam(&isdrainicecolumn,ThermalIsdrainicecolumnEnum); + + if(isdrainicecolumn) DrainWaterfraction(femmodel); if(computebasalmeltingrates) ComputeBasalMeltingrate(femmodel); Index: /issm/trunk-jpl/src/c/shared/Elements/DrainageFunctionWaterfraction.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/DrainageFunctionWaterfraction.cpp (revision 24139) +++ /issm/trunk-jpl/src/c/shared/Elements/DrainageFunctionWaterfraction.cpp (revision 24140) @@ -31,8 +31,8 @@ /*drain only up to w0*/ if(dt==0.){ - if((waterfraction>w0) && (waterfraction-Dret*ytsw0) return waterfraction-w0; else - return Dret*yts; + return Dret; } else{ Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24139) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24140) @@ -389,4 +389,5 @@ StressbalanceRiftPenaltyThresholdEnum, StressbalanceShelfDampeningEnum, + ThermalIsdrainicecolumnEnum, ThermalIsdynamicbasalspcEnum, ThermalIsenthalpyEnum, @@ -399,4 +400,5 @@ ThermalRequestedOutputsEnum, ThermalStabilizationEnum, + ThermalWatercolumnUpperlimitEnum, TimeEnum, TimesteppingCflCoefficientEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24139) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24140) @@ -397,4 +397,5 @@ case StressbalanceRiftPenaltyThresholdEnum : return "StressbalanceRiftPenaltyThreshold"; case StressbalanceShelfDampeningEnum : return "StressbalanceShelfDampening"; + case ThermalIsdrainicecolumnEnum : return "ThermalIsdrainicecolumn"; case ThermalIsdynamicbasalspcEnum : return "ThermalIsdynamicbasalspc"; case ThermalIsenthalpyEnum : return "ThermalIsenthalpy"; @@ -407,4 +408,5 @@ case ThermalRequestedOutputsEnum : return "ThermalRequestedOutputs"; case ThermalStabilizationEnum : return "ThermalStabilization"; + case ThermalWatercolumnUpperlimitEnum : return "ThermalWatercolumnUpperlimit"; case TimeEnum : return "Time"; case TimesteppingCflCoefficientEnum : return "TimesteppingCflCoefficient"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24139) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24140) @@ -406,4 +406,5 @@ else if (strcmp(name,"StressbalanceRiftPenaltyThreshold")==0) return StressbalanceRiftPenaltyThresholdEnum; else if (strcmp(name,"StressbalanceShelfDampening")==0) return StressbalanceShelfDampeningEnum; + else if (strcmp(name,"ThermalIsdrainicecolumn")==0) return ThermalIsdrainicecolumnEnum; else if (strcmp(name,"ThermalIsdynamicbasalspc")==0) return ThermalIsdynamicbasalspcEnum; else if (strcmp(name,"ThermalIsenthalpy")==0) return ThermalIsenthalpyEnum; @@ -416,4 +417,5 @@ else if (strcmp(name,"ThermalRequestedOutputs")==0) return ThermalRequestedOutputsEnum; else if (strcmp(name,"ThermalStabilization")==0) return ThermalStabilizationEnum; + else if (strcmp(name,"ThermalWatercolumnUpperlimit")==0) return ThermalWatercolumnUpperlimitEnum; else if (strcmp(name,"Time")==0) return TimeEnum; else if (strcmp(name,"TimesteppingCflCoefficient")==0) return TimesteppingCflCoefficientEnum; @@ -504,10 +506,10 @@ else if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum; else if (strcmp(name,"DeviatoricStressxx")==0) return DeviatoricStressxxEnum; - else if (strcmp(name,"DeviatoricStressxy")==0) return DeviatoricStressxyEnum; - else if (strcmp(name,"DeviatoricStressxz")==0) return DeviatoricStressxzEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum; + if (strcmp(name,"DeviatoricStressxy")==0) return DeviatoricStressxyEnum; + else if (strcmp(name,"DeviatoricStressxz")==0) return DeviatoricStressxzEnum; + else if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum; else if (strcmp(name,"DeviatoricStressyz")==0) return DeviatoricStressyzEnum; else if (strcmp(name,"DeviatoricStresszz")==0) return DeviatoricStresszzEnum; @@ -627,10 +629,10 @@ else if (strcmp(name,"MaterialsRheologyN")==0) return MaterialsRheologyNEnum; else if (strcmp(name,"MeshScaleFactor")==0) return MeshScaleFactorEnum; - else if (strcmp(name,"MeshVertexonbase")==0) return MeshVertexonbaseEnum; - else if (strcmp(name,"MeshVertexonboundary")==0) return MeshVertexonboundaryEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"MeshVertexonsurface")==0) return MeshVertexonsurfaceEnum; + if (strcmp(name,"MeshVertexonbase")==0) return MeshVertexonbaseEnum; + else if (strcmp(name,"MeshVertexonboundary")==0) return MeshVertexonboundaryEnum; + else if (strcmp(name,"MeshVertexonsurface")==0) return MeshVertexonsurfaceEnum; else if (strcmp(name,"Misfit")==0) return MisfitEnum; else if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum; @@ -750,10 +752,10 @@ else if (strcmp(name,"SmbTini")==0) return SmbTiniEnum; else if (strcmp(name,"SmbTmean")==0) return SmbTmeanEnum; - else if (strcmp(name,"SmbTz")==0) return SmbTzEnum; - else if (strcmp(name,"SmbV")==0) return SmbVEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"SmbVmean")==0) return SmbVmeanEnum; + if (strcmp(name,"SmbTz")==0) return SmbTzEnum; + else if (strcmp(name,"SmbV")==0) return SmbVEnum; + else if (strcmp(name,"SmbVmean")==0) return SmbVmeanEnum; else if (strcmp(name,"SmbVz")==0) return SmbVzEnum; else if (strcmp(name,"SmbW")==0) return SmbWEnum; @@ -873,10 +875,10 @@ else if (strcmp(name,"Outputdefinition52")==0) return Outputdefinition52Enum; else if (strcmp(name,"Outputdefinition53")==0) return Outputdefinition53Enum; - else if (strcmp(name,"Outputdefinition54")==0) return Outputdefinition54Enum; - else if (strcmp(name,"Outputdefinition55")==0) return Outputdefinition55Enum; else stage=8; } if(stage==8){ - if (strcmp(name,"Outputdefinition56")==0) return Outputdefinition56Enum; + if (strcmp(name,"Outputdefinition54")==0) return Outputdefinition54Enum; + else if (strcmp(name,"Outputdefinition55")==0) return Outputdefinition55Enum; + else if (strcmp(name,"Outputdefinition56")==0) return Outputdefinition56Enum; else if (strcmp(name,"Outputdefinition57")==0) return Outputdefinition57Enum; else if (strcmp(name,"Outputdefinition58")==0) return Outputdefinition58Enum; @@ -996,10 +998,10 @@ else if (strcmp(name,"Domain3Dsurface")==0) return Domain3DsurfaceEnum; else if (strcmp(name,"DoubleArrayInput")==0) return DoubleArrayInputEnum; - else if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum; - else if (strcmp(name,"DoubleInput")==0) return DoubleInputEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"DoubleMatArrayParam")==0) return DoubleMatArrayParamEnum; + if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum; + else if (strcmp(name,"DoubleInput")==0) return DoubleInputEnum; + else if (strcmp(name,"DoubleMatArrayParam")==0) return DoubleMatArrayParamEnum; else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum; else if (strcmp(name,"DoubleMatParam")==0) return DoubleMatParamEnum; @@ -1119,10 +1121,10 @@ else if (strcmp(name,"Matdamageice")==0) return MatdamageiceEnum; else if (strcmp(name,"Matenhancedice")==0) return MatenhancediceEnum; - else if (strcmp(name,"Materials")==0) return MaterialsEnum; - else if (strcmp(name,"Matestar")==0) return MatestarEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"Matice")==0) return MaticeEnum; + if (strcmp(name,"Materials")==0) return MaterialsEnum; + else if (strcmp(name,"Matestar")==0) return MatestarEnum; + else if (strcmp(name,"Matice")==0) return MaticeEnum; else if (strcmp(name,"Matlitho")==0) return MatlithoEnum; else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum; @@ -1242,10 +1244,10 @@ else if (strcmp(name,"StressbalanceSolution")==0) return StressbalanceSolutionEnum; else if (strcmp(name,"StressbalanceVerticalAnalysis")==0) return StressbalanceVerticalAnalysisEnum; - else if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum; - else if (strcmp(name,"StringExternalResult")==0) return StringExternalResultEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"StringParam")==0) return StringParamEnum; + if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum; + else if (strcmp(name,"StringExternalResult")==0) return StringExternalResultEnum; + else if (strcmp(name,"StringParam")==0) return StringParamEnum; else if (strcmp(name,"SubelementFriction1")==0) return SubelementFriction1Enum; else if (strcmp(name,"SubelementFriction2")==0) return SubelementFriction2Enum; Index: /issm/trunk-jpl/src/m/classes/thermal.js =================================================================== --- /issm/trunk-jpl/src/m/classes/thermal.js (revision 24139) +++ /issm/trunk-jpl/src/m/classes/thermal.js (revision 24140) @@ -28,4 +28,10 @@ //will basal boundary conditions be set dynamically this.isdynamicbasalspc=0; + + //wether waterfraction drainage is enabled + this.isdrainicecolumn=1; + + //set an upper limit for local stored watercolumn + this.watercolumn_upperlimit=1000; //Linear elements by default @@ -41,5 +47,5 @@ fielddisplay(this,'spctemperature','temperature constraints (NaN means no constraint) [K]'); - fielddisplay(this,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG, 3: anisotropic SUPG'); + fielddisplay(this,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG'); fielddisplay(this,'reltol','relative tolerance convergence criterion for enthalpy'); fielddisplay(this,'maxiter','maximum number of non linear iterations'); @@ -49,4 +55,6 @@ fielddisplay(this,'isenthalpy','use an enthalpy formulation to include temperate ice (default is 0)'); fielddisplay(this,'isdynamicbasalspc','enable dynamic setting of basal forcing. required for enthalpy formulation (default is 0)'); + fielddisplay(this,'isdrainicecolumn','wether waterfraction drainage is enabled for enthalpy formulation (default is 1)'); + fielddisplay(this,'watercolumn_upperlimit','upper limit of basal watercolumn for enthalpy formulation (default is 1000m)'); fielddisplay(this,'fe','Finite Element type: "P1" (default), "P1xP2"'); fielddisplay(this,'requested_outputs','additional outputs requested'); @@ -56,9 +64,9 @@ return "thermal"; }// }}} - this.extrude = function(md) {//{{{ - this.spctemperature=project3d(md,'vector',this.spctemperature,'type','node','layer',md.mesh.numberoflayers,'padding',NaN); - if (md.initialization.temperature.length===md.mesh.numberofvertices) { - this.spctemperature = NewArrayFill(md.mesh.numberofvertices, NaN); - var pos=ArrayFindNot(md.mesh.vertexonsurface, 0); + this.extrude = function(md) {//{{{ + this.spctemperature=project3d(md,'vector',this.spctemperature,'type','node','layer',md.mesh.numberoflayers,'padding',NaN); + if (md.initialization.temperature.length===md.mesh.numberofvertices) { + this.spctemperature = NewArrayFill(md.mesh.numberofvertices, NaN); + var pos=ArrayFindNot(md.mesh.vertexonsurface, 0); // impose observed temperature on surface for (var i=0,posIndex=0,count=0;i=',0); + for(var i=0;i=',0); %Make sure the spc are less than melting point @@ -99,5 +109,5 @@ md = checkfield(md,'fieldname','thermal.reltol','>',0.,'message','reltol must be larger than zero'); end - end + end md = checkfield(md,'fieldname','thermal.requested_outputs','stringrow',1); @@ -114,5 +124,7 @@ fielddisplay(self,'penalty_factor','scaling exponent (default is 3)'); fielddisplay(self,'isenthalpy','use an enthalpy formulation to include temperate ice (default is 0)'); - fielddisplay(self,'isdynamicbasalspc',['enable dynamic setting of basal forcing. required for enthalpy formulation (default is 0)']); + fielddisplay(self,'isdynamicbasalspc','enable dynamic setting of basal forcing. required for enthalpy formulation (default is 0)'); + fielddisplay(self,'isdrainicecolumn','wether waterfraction drainage is enabled for enthalpy formulation (default is 1)'); + fielddisplay(self,'watercolumn_upperlimit','upper limit of basal watercolumn for enthalpy formulation (default is 1000m)'); fielddisplay(self,'fe','Finite Element type: ''P1'' (default), ''P1xP2'''); fielddisplay(self,'requested_outputs','additional outputs requested'); @@ -128,4 +140,6 @@ WriteData(fid,prefix,'object',self,'fieldname','penalty_factor','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','isenthalpy','format','Boolean'); + WriteData(fid,prefix,'object',self,'fieldname','isdrainicecolumn','format','Boolean'); + WriteData(fid,prefix,'object',self,'fieldname','watercolumn_upperlimit','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','fe','format','String'); WriteData(fid,prefix,'object',self,'fieldname','isdynamicbasalspc','format','Boolean'); @@ -150,4 +164,6 @@ writejsdouble(fid,[modelname '.thermal.penalty_factor'],self.penalty_factor); writejsdouble(fid,[modelname '.thermal.isenthalpy'],self.isenthalpy); + writejsdouble(fid,[modelname '.thermal.isdrainicecolumn'],self.isdrainicecolumn); + writejsdouble(fid,[modelname '.thermal.watercolumn_upperlimit'],self.watercolumn_upperlimit); writejsdouble(fid,[modelname '.thermal.isdynamicbasalspc'],self.isdynamicbasalspc); writejscellstring(fid,[modelname '.thermal.requested_outputs'],self.requested_outputs); Index: /issm/trunk-jpl/src/m/classes/thermal.py =================================================================== --- /issm/trunk-jpl/src/m/classes/thermal.py (revision 24139) +++ /issm/trunk-jpl/src/m/classes/thermal.py (revision 24140) @@ -23,4 +23,6 @@ self.isenthalpy = 0 self.isdynamicbasalspc = 0 + isdrainicecolumn = 0; + watercolumn_upperlimit = 0; self.fe = 'P1' self.requested_outputs = [] @@ -33,5 +35,5 @@ string=' Thermal solution parameters:' string="%s\n%s"%(string,fielddisplay(self,'spctemperature','temperature constraints (NaN means no constraint) [K]')) - string="%s\n%s"%(string,fielddisplay(self,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG, 3: anisotropic SUPG')) + string="%s\n%s"%(string,fielddisplay(self,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG')) string="%s\n%s"%(string,fielddisplay(self,'maxiter','maximum number of non linear iterations')) string="%s\n%s"%(string,fielddisplay(self,'reltol','relative tolerance criterion')) @@ -40,4 +42,6 @@ string="%s\n%s"%(string,fielddisplay(self,'isenthalpy','use an enthalpy formulation to include temperate ice (default is 0)')) string="%s\n%s"%(string,fielddisplay(self,'isdynamicbasalspc','enable dynamic setting of basal forcing. required for enthalpy formulation (default is 0)')) + string="%s\n%s"%(string,fielddisplay(self,'isdrainicecolumn','wether waterfraction drainage is enabled for enthalpy formulation (default is 1)') + string="%s\n%s"%(string,fielddisplay(self,'watercolumn_upperlimit','upper limit of basal watercolumn for enthalpy formulation (default is 1000m)') string="%s\n%s"%(string,fielddisplay(self,'requested_outputs','additional outputs requested')) return string @@ -82,4 +86,10 @@ self.isdynamicbasalspc=0 + #wether waterfraction drainage is enabled + self.isdrainicecolumn=1 + + #set an upper limit for local stored watercolumn + self.watercolumn_upperlimit=1000 + #Finite element interpolation self.fe='P1' @@ -96,9 +106,12 @@ return md - md = checkfield(md,'fieldname','thermal.stabilization','numel',[1],'values',[0,1,2,3]) + md = checkfield(md,'fieldname','thermal.stabilization','numel',[1],'values',[0,1,2]) md = checkfield(md,'fieldname','thermal.spctemperature','Inf',1,'timeseries',1) md = checkfield(md,'fieldname','thermal.requested_outputs','stringrow',1) if 'EnthalpyAnalysis' in analyses and md.thermal.isenthalpy and md.mesh.dimension()==3: + md = checkfield(md,'fieldname','thermal.isdrainicecolumn','numel',[1],'values',[0,1]) + md = checkfield(md,'fieldname','thermal.watercolumn_upperlimit','>=',0) + TEMP = md.thermal.spctemperature[:-1].flatten(-1) pos=np.where(~np.isnan(TEMP)) @@ -114,9 +127,9 @@ md = checkfield(md,'fieldname','thermal.spctemperature','field',md.thermal.spctemperature.flatten(-1)[pos],'<=',control[pos],'message',"spctemperature should be below the adjusted melting point") md = checkfield(md,'fieldname','thermal.isenthalpy','numel',[1],'values',[0,1]) - md = checkfield(md,'fieldname','thermal.isdynamicbasalspc','numel',[1],'values',[0,1]); + md = checkfield(md,'fieldname','thermal.isdynamicbasalspc','numel',[1],'values',[0,1]) if(md.thermal.isenthalpy): if np.isnan(md.stressbalance.reltol): md.checkmessage("for a steadystate computation, thermal.reltol (relative convergence criterion) must be defined!") - md = checkfield(md,'fieldname','thermal.reltol','>',0.,'message',"reltol must be larger than zero"); + md = checkfield(md,'fieldname','thermal.reltol','>',0.,'message',"reltol must be larger than zero") @@ -127,11 +140,13 @@ WriteData(fid,prefix,'object',self,'fieldname','penalty_threshold','format','Integer') WriteData(fid,prefix,'object',self,'fieldname','stabilization','format','Integer') - WriteData(fid,prefix,'object',self,'fieldname','reltol','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','reltol','format','Double') WriteData(fid,prefix,'object',self,'fieldname','maxiter','format','Integer') WriteData(fid,prefix,'object',self,'fieldname','penalty_lock','format','Integer') WriteData(fid,prefix,'object',self,'fieldname','penalty_factor','format','Double') WriteData(fid,prefix,'object',self,'fieldname','isenthalpy','format','Boolean') - WriteData(fid,prefix,'object',self,'fieldname','fe','format','String'); - WriteData(fid,prefix,'object',self,'fieldname','isdynamicbasalspc','format','Boolean'); + WriteData(fid,prefix,'object',self,'fieldname','isdrainicecolumn','format','Boolean') + WriteData(fid,prefix,'object',self,'fieldname','watercolumn_upperlimit','format','Double') + WriteData(fid,prefix,'object',self,'fieldname','fe','format','String') + WriteData(fid,prefix,'object',self,'fieldname','isdynamicbasalspc','format','Boolean') #process requested outputs