Changeset 24140


Ignore:
Timestamp:
09/11/19 09:30:58 (6 years ago)
Author:
rueckamp
Message:

CHG: add drainiceolumn and watercolumn_upperlimit as user input; improve drainage for steadystate case

Location:
issm/trunk-jpl
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp

    r24136 r24140  
    275275        parameters->AddObject(iomodel->CopyConstantObject("md.thermal.isenthalpy",ThermalIsenthalpyEnum));
    276276        parameters->AddObject(iomodel->CopyConstantObject("md.thermal.isdynamicbasalspc",ThermalIsdynamicbasalspcEnum));
     277        parameters->AddObject(iomodel->CopyConstantObject("md.thermal.isdrainicecolumn",ThermalIsdrainicecolumnEnum));
     278        parameters->AddObject(iomodel->CopyConstantObject("md.thermal.watercolumn_upperlimit",ThermalWatercolumnUpperlimitEnum));
    277279        parameters->AddObject(iomodel->CopyConstantObject("md.friction.law",FrictionLawEnum));
    278280
    … …  
    424426        element->GetInputListOnNodes(basalmeltingrates,BasalforcingsGroundediceMeltingRateEnum);
    425427
     428        IssmDouble watercolumnupperlimit = element->FindParam(ThermalWatercolumnUpperlimitEnum);
     429       
    426430        Gauss* gauss=element->NewGauss();
    427431        for(is=0;is<numsegments;is++){
    … …  
    484488                                watercolumns[nodedown]+=dt*meltingrate_enthalpy[is];
    485489                        }
     490                        if(watercolumns[nodedown]>watercolumnupperlimit) watercolumns[nodedown]=watercolumnupperlimit;
    486491                }
    487492                else{
    … …  
    11681173
    11691174                for(k=0;k<numnodes;k++){
    1170                         waterfractions[k]-=dt*drainage[k];
     1175                        if(dt==0.)
     1176                                waterfractions[k]-=drainage[k];
     1177                        else
     1178                                waterfractions[k]-=dt*drainage[k];
     1179                       
    11711180                        element->ThermalToEnthalpy(&enthalpies[k], temperatures[k], waterfractions[k], pressures[k]);
    11721181                }
    … …  
    11871196        IssmDouble thermalconductivity      = element->FindParam(MaterialsThermalconductivityEnum);
    11881197
    1189         if(enthalpy < PureIceEnthalpy(element,pressure)){
     1198        if(enthalpy < PureIceEnthalpy(element,pressure))
    11901199                return thermalconductivity/heatcapacity;
    1191         }
    1192         else{
     1200        else
    11931201                return temperateiceconductivity/heatcapacity;
    1194         }
    11951202}/*}}}*/
    11961203IssmDouble     EnthalpyAnalysis::EnthalpyDiffusionParameterVolume(Element* element,int enthalpy_enum){/*{{{*/
    … …  
    16781685        /*Intermediaries*/
    16791686        bool computebasalmeltingrates=true;
    1680         bool drainicecolumn=true;
     1687        bool isdrainicecolumn;
    16811688        IssmDouble dt;
    16821689
    16831690        femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    1684 
    1685         if(drainicecolumn && (dt>0.))   DrainWaterfraction(femmodel);
     1691        femmodel->parameters->FindParam(&isdrainicecolumn,ThermalIsdrainicecolumnEnum);
     1692
     1693        if(isdrainicecolumn)    DrainWaterfraction(femmodel);
    16861694        if(computebasalmeltingrates)    ComputeBasalMeltingrate(femmodel);
    16871695

  • issm/trunk-jpl/src/c/shared/Elements/DrainageFunctionWaterfraction.cpp

    r23066 r24140  
    3131        /*drain only up to w0*/
    3232        if(dt==0.){
    33                 if((waterfraction>w0) && (waterfraction-Dret*yts<w0))
     33                if(waterfraction>w0)
    3434                        return waterfraction-w0;
    3535                else
    36                         return Dret*yts;
     36                        return Dret;
    3737        }
    3838        else{

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

    r24139 r24140  
    389389        StressbalanceRiftPenaltyThresholdEnum,
    390390        StressbalanceShelfDampeningEnum,
     391        ThermalIsdrainicecolumnEnum,
    391392        ThermalIsdynamicbasalspcEnum,
    392393        ThermalIsenthalpyEnum,
    … …  
    399400        ThermalRequestedOutputsEnum,
    400401        ThermalStabilizationEnum,
     402        ThermalWatercolumnUpperlimitEnum,
    401403        TimeEnum,
    402404        TimesteppingCflCoefficientEnum,

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

    r24139 r24140  
    397397                case StressbalanceRiftPenaltyThresholdEnum : return "StressbalanceRiftPenaltyThreshold";
    398398                case StressbalanceShelfDampeningEnum : return "StressbalanceShelfDampening";
     399                case ThermalIsdrainicecolumnEnum : return "ThermalIsdrainicecolumn";
    399400                case ThermalIsdynamicbasalspcEnum : return "ThermalIsdynamicbasalspc";
    400401                case ThermalIsenthalpyEnum : return "ThermalIsenthalpy";
    … …  
    407408                case ThermalRequestedOutputsEnum : return "ThermalRequestedOutputs";
    408409                case ThermalStabilizationEnum : return "ThermalStabilization";
     410                case ThermalWatercolumnUpperlimitEnum : return "ThermalWatercolumnUpperlimit";
    409411                case TimeEnum : return "Time";
    410412                case TimesteppingCflCoefficientEnum : return "TimesteppingCflCoefficient";

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

    r24139 r24140  
    406406              else if (strcmp(name,"StressbalanceRiftPenaltyThreshold")==0) return StressbalanceRiftPenaltyThresholdEnum;
    407407              else if (strcmp(name,"StressbalanceShelfDampening")==0) return StressbalanceShelfDampeningEnum;
     408              else if (strcmp(name,"ThermalIsdrainicecolumn")==0) return ThermalIsdrainicecolumnEnum;
    408409              else if (strcmp(name,"ThermalIsdynamicbasalspc")==0) return ThermalIsdynamicbasalspcEnum;
    409410              else if (strcmp(name,"ThermalIsenthalpy")==0) return ThermalIsenthalpyEnum;
    … …  
    416417              else if (strcmp(name,"ThermalRequestedOutputs")==0) return ThermalRequestedOutputsEnum;
    417418              else if (strcmp(name,"ThermalStabilization")==0) return ThermalStabilizationEnum;
     419              else if (strcmp(name,"ThermalWatercolumnUpperlimit")==0) return ThermalWatercolumnUpperlimitEnum;
    418420              else if (strcmp(name,"Time")==0) return TimeEnum;
    419421              else if (strcmp(name,"TimesteppingCflCoefficient")==0) return TimesteppingCflCoefficientEnum;
    … …  
    504506              else if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum;
    505507              else if (strcmp(name,"DeviatoricStressxx")==0) return DeviatoricStressxxEnum;
    506               else if (strcmp(name,"DeviatoricStressxy")==0) return DeviatoricStressxyEnum;
    507               else if (strcmp(name,"DeviatoricStressxz")==0) return DeviatoricStressxzEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum;
     511              if (strcmp(name,"DeviatoricStressxy")==0) return DeviatoricStressxyEnum;
     512              else if (strcmp(name,"DeviatoricStressxz")==0) return DeviatoricStressxzEnum;
     513              else if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum;
    512514              else if (strcmp(name,"DeviatoricStressyz")==0) return DeviatoricStressyzEnum;
    513515              else if (strcmp(name,"DeviatoricStresszz")==0) return DeviatoricStresszzEnum;
    … …  
    627629              else if (strcmp(name,"MaterialsRheologyN")==0) return MaterialsRheologyNEnum;
    628630              else if (strcmp(name,"MeshScaleFactor")==0) return MeshScaleFactorEnum;
    629               else if (strcmp(name,"MeshVertexonbase")==0) return MeshVertexonbaseEnum;
    630               else if (strcmp(name,"MeshVertexonboundary")==0) return MeshVertexonboundaryEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"MeshVertexonsurface")==0) return MeshVertexonsurfaceEnum;
     634              if (strcmp(name,"MeshVertexonbase")==0) return MeshVertexonbaseEnum;
     635              else if (strcmp(name,"MeshVertexonboundary")==0) return MeshVertexonboundaryEnum;
     636              else if (strcmp(name,"MeshVertexonsurface")==0) return MeshVertexonsurfaceEnum;
    635637              else if (strcmp(name,"Misfit")==0) return MisfitEnum;
    636638              else if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum;
    … …  
    750752              else if (strcmp(name,"SmbTini")==0) return SmbTiniEnum;
    751753              else if (strcmp(name,"SmbTmean")==0) return SmbTmeanEnum;
    752               else if (strcmp(name,"SmbTz")==0) return SmbTzEnum;
    753               else if (strcmp(name,"SmbV")==0) return SmbVEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"SmbVmean")==0) return SmbVmeanEnum;
     757              if (strcmp(name,"SmbTz")==0) return SmbTzEnum;
     758              else if (strcmp(name,"SmbV")==0) return SmbVEnum;
     759              else if (strcmp(name,"SmbVmean")==0) return SmbVmeanEnum;
    758760              else if (strcmp(name,"SmbVz")==0) return SmbVzEnum;
    759761              else if (strcmp(name,"SmbW")==0) return SmbWEnum;
    … …  
    873875              else if (strcmp(name,"Outputdefinition52")==0) return Outputdefinition52Enum;
    874876              else if (strcmp(name,"Outputdefinition53")==0) return Outputdefinition53Enum;
    875               else if (strcmp(name,"Outputdefinition54")==0) return Outputdefinition54Enum;
    876               else if (strcmp(name,"Outputdefinition55")==0) return Outputdefinition55Enum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"Outputdefinition56")==0) return Outputdefinition56Enum;
     880              if (strcmp(name,"Outputdefinition54")==0) return Outputdefinition54Enum;
     881              else if (strcmp(name,"Outputdefinition55")==0) return Outputdefinition55Enum;
     882              else if (strcmp(name,"Outputdefinition56")==0) return Outputdefinition56Enum;
    881883              else if (strcmp(name,"Outputdefinition57")==0) return Outputdefinition57Enum;
    882884              else if (strcmp(name,"Outputdefinition58")==0) return Outputdefinition58Enum;
    … …  
    996998              else if (strcmp(name,"Domain3Dsurface")==0) return Domain3DsurfaceEnum;
    997999              else if (strcmp(name,"DoubleArrayInput")==0) return DoubleArrayInputEnum;
    998               else if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum;
    999               else if (strcmp(name,"DoubleInput")==0) return DoubleInputEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"DoubleMatArrayParam")==0) return DoubleMatArrayParamEnum;
     1003              if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum;
     1004              else if (strcmp(name,"DoubleInput")==0) return DoubleInputEnum;
     1005              else if (strcmp(name,"DoubleMatArrayParam")==0) return DoubleMatArrayParamEnum;
    10041006              else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
    10051007              else if (strcmp(name,"DoubleMatParam")==0) return DoubleMatParamEnum;
    … …  
    11191121              else if (strcmp(name,"Matdamageice")==0) return MatdamageiceEnum;
    11201122              else if (strcmp(name,"Matenhancedice")==0) return MatenhancediceEnum;
    1121               else if (strcmp(name,"Materials")==0) return MaterialsEnum;
    1122               else if (strcmp(name,"Matestar")==0) return MatestarEnum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"Matice")==0) return MaticeEnum;
     1126              if (strcmp(name,"Materials")==0) return MaterialsEnum;
     1127              else if (strcmp(name,"Matestar")==0) return MatestarEnum;
     1128              else if (strcmp(name,"Matice")==0) return MaticeEnum;
    11271129              else if (strcmp(name,"Matlitho")==0) return MatlithoEnum;
    11281130              else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum;
    … …  
    12421244              else if (strcmp(name,"StressbalanceSolution")==0) return StressbalanceSolutionEnum;
    12431245              else if (strcmp(name,"StressbalanceVerticalAnalysis")==0) return StressbalanceVerticalAnalysisEnum;
    1244               else if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum;
    1245               else if (strcmp(name,"StringExternalResult")==0) return StringExternalResultEnum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"StringParam")==0) return StringParamEnum;
     1249              if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum;
     1250              else if (strcmp(name,"StringExternalResult")==0) return StringExternalResultEnum;
     1251              else if (strcmp(name,"StringParam")==0) return StringParamEnum;
    12501252              else if (strcmp(name,"SubelementFriction1")==0) return SubelementFriction1Enum;
    12511253              else if (strcmp(name,"SubelementFriction2")==0) return SubelementFriction2Enum;

  • issm/trunk-jpl/src/m/classes/thermal.js

    r24136 r24140  
    2828                //will basal boundary conditions be set dynamically
    2929                this.isdynamicbasalspc=0;
     30
     31                //wether waterfraction drainage is enabled
     32                this.isdrainicecolumn=1;
     33
     34                //set an upper limit for local stored watercolumn
     35                this.watercolumn_upperlimit=1000;
    3036               
    3137                //Linear elements by default
    … …  
    4147
    4248                fielddisplay(this,'spctemperature','temperature constraints (NaN means no constraint) [K]');
    43                 fielddisplay(this,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG, 3: anisotropic SUPG');
     49                fielddisplay(this,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG');
    4450                fielddisplay(this,'reltol','relative tolerance convergence criterion for enthalpy');
    4551                fielddisplay(this,'maxiter','maximum number of non linear iterations');
    … …  
    4955                fielddisplay(this,'isenthalpy','use an enthalpy formulation to include temperate ice (default is 0)');
    5056                fielddisplay(this,'isdynamicbasalspc','enable dynamic setting of basal forcing. required for enthalpy formulation (default is 0)');
     57                fielddisplay(this,'isdrainicecolumn','wether waterfraction drainage is enabled for enthalpy formulation (default is 1)');
     58                fielddisplay(this,'watercolumn_upperlimit','upper limit of basal watercolumn for enthalpy formulation (default is 1000m)');
    5159                fielddisplay(this,'fe','Finite Element type: "P1" (default), "P1xP2"');
    5260                fielddisplay(this,'requested_outputs','additional outputs requested');
    … …  
    5664                return "thermal";
    5765        }// }}}
    58     this.extrude = function(md) {//{{{
    59         this.spctemperature=project3d(md,'vector',this.spctemperature,'type','node','layer',md.mesh.numberoflayers,'padding',NaN);
    60         if (md.initialization.temperature.length===md.mesh.numberofvertices) {
    61             this.spctemperature = NewArrayFill(md.mesh.numberofvertices, NaN);
    62             var pos=ArrayFindNot(md.mesh.vertexonsurface, 0);
     66        this.extrude = function(md) {//{{{
     67                this.spctemperature=project3d(md,'vector',this.spctemperature,'type','node','layer',md.mesh.numberoflayers,'padding',NaN);
     68                if (md.initialization.temperature.length===md.mesh.numberofvertices) {
     69                        this.spctemperature = NewArrayFill(md.mesh.numberofvertices, NaN);
     70                        var pos=ArrayFindNot(md.mesh.vertexonsurface, 0);
    6371                        // impose observed temperature on surface
    6472                        for (var i=0,posIndex=0,count=0;i<md.initialization.temperature.length;i++){
    … …  
    8189//                              }
    8290//                      }
    83         }
     91                }
    8492
    85         return this;
    86     }//}}}
     93                return this;
     94        }//}}}
    8795        this.checkconsistency = function(md,solution,analyses){ // {{{
    8896
    … …  
    9098                if(!ArrayAnyEqual(ArrayIsMember('ThermalAnalysis',analyses),1) & !ArrayAnyEqual(ArrayIsMember('EnthalpyAnalysis',analyses),1)  | (solution == 'TransientSolution' & md.trans.isthermal==0)) return;
    9199
    92                 checkfield(md,'fieldname','thermal.stabilization','numel',[1],'values',[0 ,1, 2, 3]);
     100                checkfield(md,'fieldname','thermal.stabilization','numel',[1],'values',[0 ,1, 2]);
    93101                checkfield(md,'fieldname','thermal.spctemperature','Inf',1,'timeseries',1);
    94102                checkfield(md,'fieldname','thermal.fe','values',['P1','P1xP2','P1xP3']);
    95103                if(ArrayAnyEqual(ArrayIsMember('EnthalpyAnalysis',analyses),1) & md.thermal.isenthalpy & md.mesh.dimension() == 3){
    96                        
     104                        checkfield(md,'fieldname','thermal.isdrainicecolumn','numel',[1],'values',[0, 1]);
     105                        checkfield(md,'fieldname','thermal.watercolumn_upperlimit','>=',0);
     106
    97107                        for(var i=0;i<md.mesh.numberofvertices;i++){
    98108                                for(var j=0;j<md.thermal.spctemperature[0].length;j++){
    … …  
    127137                        WriteData(fid,prefix,'object',this,'fieldname','penalty_factor','format','Double');
    128138                        WriteData(fid,prefix,'object',this,'fieldname','isenthalpy','format','Boolean');
     139                        WriteData(fid,prefix,'object',this,'fieldname','isdrainicecolumn','format','Boolean');
     140                        WriteData(fid,prefix,'object',this,'fieldname','watercolumn_upperlimit','format','Double');
    129141                        WriteData(fid,prefix,'object',this,'fieldname','fe','format','String');
    130142                        WriteData(fid,prefix,'object',this,'fieldname','isdynamicbasalspc','format','Boolean');
    … …  
    155167        this.penalty_threshold = 0;
    156168        this.stabilization     = 0;
    157         this.reltol                        = 0;
     169        this.reltol        = 0;
    158170        this.maxiter           = 0;
    159171        this.penalty_lock      = 0;
    … …  
    161173        this.isenthalpy        = 0;
    162174        this.isdynamicbasalspc = 0;
     175        this.isdrainicecolumn  = 0;
     176        this.watercolumn_upperlimit=0;
    163177        this.fe                = 'P1';
    164178        this.requested_outputs = [];

  • issm/trunk-jpl/src/m/classes/thermal.m

    r24136 r24140  
    1515                isenthalpy        = 0;
    1616                isdynamicbasalspc = 0;
     17                isdrainicecolumn   = 0;
     18                watercolumn_upperlimit= 0;
    1719                fe                = 'P1';
    1820                requested_outputs = {};
    … …  
    6668                        %will basal boundary conditions be set dynamically
    6769                        self.isdynamicbasalspc=0;
     70               
     71                        %wether waterfraction drainage is enabled
     72                        self.isdrainicecolumn=1;
     73
     74                        %set an upper limit for local stored watercolumn
     75                        self.watercolumn_upperlimit=1000;
    6876
    6977                        %Linear elements by default
    … …  
    8290                        md = checkfield(md,'fieldname','thermal.fe','values',{'P1','P1xP2','P1xP3'});
    8391                        if (ismember('EnthalpyAnalysis',analyses) & md.thermal.isenthalpy & dimension(md.mesh)==3),
     92                                md = checkfield(md,'fieldname','thermal.isdrainicecolumn','numel',[1],'values',[0 1]);
     93                                md = checkfield(md,'fieldname','thermal.watercolumn_upperlimit','>=',0);
    8494
    8595                                %Make sure the spc are less than melting point
    … …  
    99109                                        md = checkfield(md,'fieldname','thermal.reltol','>',0.,'message','reltol must be larger than zero');
    100110                                end
    101             end
     111                        end
    102112
    103113                 md = checkfield(md,'fieldname','thermal.requested_outputs','stringrow',1);
    … …  
    114124                        fielddisplay(self,'penalty_factor','scaling exponent (default is 3)');
    115125                        fielddisplay(self,'isenthalpy','use an enthalpy formulation to include temperate ice (default is 0)');
    116                         fielddisplay(self,'isdynamicbasalspc',['enable dynamic setting of basal forcing. required for enthalpy formulation (default is 0)']);
     126                        fielddisplay(self,'isdynamicbasalspc','enable dynamic setting of basal forcing. required for enthalpy formulation (default is 0)');
     127                        fielddisplay(self,'isdrainicecolumn','wether waterfraction drainage is enabled for enthalpy formulation (default is 1)');
     128                        fielddisplay(self,'watercolumn_upperlimit','upper limit of basal watercolumn for enthalpy formulation (default is 1000m)');
    117129                        fielddisplay(self,'fe','Finite Element type: ''P1'' (default), ''P1xP2''');
    118130                        fielddisplay(self,'requested_outputs','additional outputs requested');
    … …  
    128140                        WriteData(fid,prefix,'object',self,'fieldname','penalty_factor','format','Double');
    129141                        WriteData(fid,prefix,'object',self,'fieldname','isenthalpy','format','Boolean');
     142                        WriteData(fid,prefix,'object',self,'fieldname','isdrainicecolumn','format','Boolean');
     143                        WriteData(fid,prefix,'object',self,'fieldname','watercolumn_upperlimit','format','Double');
    130144                        WriteData(fid,prefix,'object',self,'fieldname','fe','format','String');
    131145                        WriteData(fid,prefix,'object',self,'fieldname','isdynamicbasalspc','format','Boolean');
    … …  
    150164                        writejsdouble(fid,[modelname '.thermal.penalty_factor'],self.penalty_factor);
    151165                        writejsdouble(fid,[modelname '.thermal.isenthalpy'],self.isenthalpy);
     166                        writejsdouble(fid,[modelname '.thermal.isdrainicecolumn'],self.isdrainicecolumn);
     167                        writejsdouble(fid,[modelname '.thermal.watercolumn_upperlimit'],self.watercolumn_upperlimit);
    152168                        writejsdouble(fid,[modelname '.thermal.isdynamicbasalspc'],self.isdynamicbasalspc);
    153169                        writejscellstring(fid,[modelname '.thermal.requested_outputs'],self.requested_outputs);

  • issm/trunk-jpl/src/m/classes/thermal.py

    r24136 r24140  
    2323                self.isenthalpy        = 0
    2424                self.isdynamicbasalspc = 0
     25                isdrainicecolumn       = 0;
     26                watercolumn_upperlimit = 0;
    2527                self.fe                = 'P1'
    2628                self.requested_outputs = []
    … …  
    3335                string='   Thermal solution parameters:'
    3436                string="%s\n%s"%(string,fielddisplay(self,'spctemperature','temperature constraints (NaN means no constraint) [K]'))
    35                 string="%s\n%s"%(string,fielddisplay(self,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG, 3: anisotropic SUPG'))
     37                string="%s\n%s"%(string,fielddisplay(self,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG'))
    3638                string="%s\n%s"%(string,fielddisplay(self,'maxiter','maximum number of non linear iterations'))
    3739                string="%s\n%s"%(string,fielddisplay(self,'reltol','relative tolerance criterion'))
    … …  
    4042                string="%s\n%s"%(string,fielddisplay(self,'isenthalpy','use an enthalpy formulation to include temperate ice (default is 0)'))
    4143                string="%s\n%s"%(string,fielddisplay(self,'isdynamicbasalspc','enable dynamic setting of basal forcing. required for enthalpy formulation (default is 0)'))
     44                string="%s\n%s"%(string,fielddisplay(self,'isdrainicecolumn','wether waterfraction drainage is enabled for enthalpy formulation (default is 1)')
     45                string="%s\n%s"%(string,fielddisplay(self,'watercolumn_upperlimit','upper limit of basal watercolumn for enthalpy formulation (default is 1000m)')
    4246                string="%s\n%s"%(string,fielddisplay(self,'requested_outputs','additional outputs requested'))
    4347                return string
    … …  
    8286                self.isdynamicbasalspc=0
    8387
     88                #wether waterfraction drainage is enabled
     89                self.isdrainicecolumn=1
     90
     91                #set an upper limit for local stored watercolumn
     92                self.watercolumn_upperlimit=1000
     93
    8494                #Finite element interpolation
    8595                self.fe='P1'
    … …  
    96106                        return md
    97107
    98                 md = checkfield(md,'fieldname','thermal.stabilization','numel',[1],'values',[0,1,2,3])
     108                md = checkfield(md,'fieldname','thermal.stabilization','numel',[1],'values',[0,1,2])
    99109                md = checkfield(md,'fieldname','thermal.spctemperature','Inf',1,'timeseries',1)
    100110                md = checkfield(md,'fieldname','thermal.requested_outputs','stringrow',1)
    101111
    102112                if 'EnthalpyAnalysis' in analyses and md.thermal.isenthalpy and md.mesh.dimension()==3:
     113                        md = checkfield(md,'fieldname','thermal.isdrainicecolumn','numel',[1],'values',[0,1])
     114                        md = checkfield(md,'fieldname','thermal.watercolumn_upperlimit','>=',0)
     115
    103116                        TEMP = md.thermal.spctemperature[:-1].flatten(-1)
    104117                        pos=np.where(~np.isnan(TEMP))
    … …  
    114127                        md = checkfield(md,'fieldname','thermal.spctemperature','field',md.thermal.spctemperature.flatten(-1)[pos],'<=',control[pos],'message',"spctemperature should be below the adjusted melting point")
    115128                        md = checkfield(md,'fieldname','thermal.isenthalpy','numel',[1],'values',[0,1])
    116                         md = checkfield(md,'fieldname','thermal.isdynamicbasalspc','numel',[1],'values',[0,1]);
     129                        md = checkfield(md,'fieldname','thermal.isdynamicbasalspc','numel',[1],'values',[0,1])
    117130                        if(md.thermal.isenthalpy):
    118131                                if np.isnan(md.stressbalance.reltol):
    119132                                        md.checkmessage("for a steadystate computation, thermal.reltol (relative convergence criterion) must be defined!")
    120                                 md = checkfield(md,'fieldname','thermal.reltol','>',0.,'message',"reltol must be larger than zero");
     133                                md = checkfield(md,'fieldname','thermal.reltol','>',0.,'message',"reltol must be larger than zero")
    121134
    122135
    … …  
    127140                WriteData(fid,prefix,'object',self,'fieldname','penalty_threshold','format','Integer')
    128141                WriteData(fid,prefix,'object',self,'fieldname','stabilization','format','Integer')
    129                 WriteData(fid,prefix,'object',self,'fieldname','reltol','format','Double');
     142                WriteData(fid,prefix,'object',self,'fieldname','reltol','format','Double')
    130143                WriteData(fid,prefix,'object',self,'fieldname','maxiter','format','Integer')
    131144                WriteData(fid,prefix,'object',self,'fieldname','penalty_lock','format','Integer')
    132145                WriteData(fid,prefix,'object',self,'fieldname','penalty_factor','format','Double')
    133146                WriteData(fid,prefix,'object',self,'fieldname','isenthalpy','format','Boolean')
    134                 WriteData(fid,prefix,'object',self,'fieldname','fe','format','String');
    135                 WriteData(fid,prefix,'object',self,'fieldname','isdynamicbasalspc','format','Boolean');
     147                WriteData(fid,prefix,'object',self,'fieldname','isdrainicecolumn','format','Boolean')
     148                WriteData(fid,prefix,'object',self,'fieldname','watercolumn_upperlimit','format','Double')
     149                WriteData(fid,prefix,'object',self,'fieldname','fe','format','String')
     150                WriteData(fid,prefix,'object',self,'fieldname','isdynamicbasalspc','format','Boolean')
    136151
    137152                #process requested outputs
Note: See TracChangeset for help on using the changeset viewer.