Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 21625) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 21626) @@ -61,5 +61,4 @@ break; case CalvingDevEnum: - iomodel->FetchDataToInput(elements,"md.calving.coeff",CalvingdevCoeffEnum); iomodel->FetchDataToInput(elements,"md.calving.meltingrate",CalvingMeltingrateEnum); break; @@ -75,4 +74,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.levelset.stabilization",LevelsetStabilizationEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.levelset.reinit_frequency",LevelsetReinitFrequencyEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.levelset.calving_max",CalvingMaxEnum)); int calvinglaw; iomodel->FindConstant(&calvinglaw,"md.calving.law"); @@ -80,5 +80,8 @@ case DefaultCalvingEnum: case CalvingLevermannEnum: + break; case CalvingDevEnum: + parameters->AddObject(iomodel->CopyConstantObject("md.calving.stress_threshold_groundedice",CalvingStressThresholdGroundediceEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.calving.stress_threshold_floatingice",CalvingStressThresholdFloatingiceEnum)); break; case CalvingMinthicknessEnum: @@ -138,5 +141,6 @@ IssmDouble h,hx,hy,hz; IssmDouble vel; - IssmDouble norm_dlsf, norm_calving, calvingrate, meltingrate; + IssmDouble norm_dlsf, norm_calving, calvingrate, meltingrate, groundedice; + IssmDouble calvingmax; IssmDouble* xyz_list = NULL; @@ -152,4 +156,6 @@ } + /*Calving threshold*/ + /*Fetch number of nodes and dof for this finite element*/ int numnodes = basalelement->GetNumberOfNodes(); @@ -170,4 +176,5 @@ basalelement->GetVerticesCoordinates(&xyz_list); basalelement->FindParam(&dt,TimesteppingTimeStepEnum); + basalelement->FindParam(&calvingmax,CalvingMaxEnum); Input* vx_input = NULL; Input* vy_input = NULL; @@ -178,4 +185,5 @@ Input* calvingrate_input = NULL; Input* meltingrate_input = NULL; + Input* gr_input = NULL; /*Load velocities*/ @@ -187,8 +195,10 @@ vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input); vy_input=basalelement->GetInput(VyEnum); _assert_(vy_input); + gr_input=basalelement->GetInput(MaskGroundediceLevelsetEnum); _assert_(gr_input); break; case Domain3DEnum: vx_input=basalelement->GetInput(VxAverageEnum); _assert_(vx_input); vy_input=basalelement->GetInput(VyAverageEnum); _assert_(vy_input); + gr_input=basalelement->GetInput(MaskGroundediceLevelsetEnum); _assert_(gr_input); break; default: _error_("mesh "<GetInputValue(&v[0],gauss); vy_input->GetInputValue(&v[1],gauss); + gr_input->GetInputValue(&groundedice,gauss); /*Get calving speed*/ @@ -262,4 +273,9 @@ calvingrate_input->GetInputValue(&calvingrate,gauss); meltingrate_input->GetInputValue(&meltingrate,gauss); + + /*Limit calving rate to c <= v + 3 km/yr */ + vel=sqrt(v[0]*v[0] + v[1]*v[1]); + if(calvingrate>calvingmax+vel) calvingrate = vel+calvingmax; + if(groundedice<0) meltingrate = 0; norm_dlsf=0.; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 21625) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 21626) @@ -186,5 +186,6 @@ IssmDouble calvingrate[NUMVERTICES]; IssmDouble lambda1,lambda2,ex,ey,vx,vy,vel; - IssmDouble sigma_vm,sigma_max,epse_2,groundedice; + IssmDouble sigma_vm,sigma_max,sigma_max_floating,sigma_max_grounded; + IssmDouble epse_2,groundedice; /* Get node coordinates and dof list: */ @@ -202,4 +203,6 @@ IssmDouble B = this->GetMaterialParameter(MaterialsRheologyBbarEnum); IssmDouble n = this->GetMaterialParameter(MaterialsRheologyNEnum); + this->parameters->FindParam(&sigma_max_floating,CalvingStressThresholdFloatingiceEnum); + this->parameters->FindParam(&sigma_max_grounded,CalvingStressThresholdGroundediceEnum); /* Start looping on the number of vertices: */ @@ -229,7 +232,10 @@ epse_2 = 1./2. *(lambda1*lambda1 + lambda2*lambda2); sigma_vm = sqrt(3.) * B * pow(epse_2,1./(2.*n)); - sigma_max = 1000.e+3; - - if(groundedice<0) sigma_max=200.e+3; + + /*Tensile stress threshold*/ + if(groundedice<0) + sigma_max = sigma_max_floating; + else + sigma_max = sigma_max_grounded; /*Assign values*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 21625) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 21626) @@ -216,5 +216,6 @@ IssmDouble calvingrate[NUMVERTICES]; IssmDouble lambda1,lambda2,ex,ey,vx,vy,vel; - IssmDouble sigma_vm,sigma_max,epse_2,groundedice; + IssmDouble sigma_vm,sigma_max,sigma_max_floating,sigma_max_grounded; + IssmDouble epse_2,groundedice; /* Get node coordinates and dof list: */ @@ -227,4 +228,6 @@ IssmDouble B = this->GetMaterialParameter(MaterialsRheologyBbarEnum); IssmDouble n = this->GetMaterialParameter(MaterialsRheologyNEnum); + this->parameters->FindParam(&sigma_max_floating,CalvingStressThresholdFloatingiceEnum); + this->parameters->FindParam(&sigma_max_grounded,CalvingStressThresholdGroundediceEnum); /* Start looping on the number of vertices: */ @@ -254,17 +257,15 @@ epse_2 = 1./2. *(lambda1*lambda1 + lambda2*lambda2); sigma_vm = sqrt(3.) * B * pow(epse_2,1./(2.*n)); - //sigma_max = 125.e+3; - sigma_max = 350.e+3; - sigma_max = 450.e+3; - sigma_max = 800.e+3; //too much - //sigma_max = 700.e+3; - //sigma_max = 670.e+3; - //sigma_max = 550.e+3; - sigma_max = 750.e+3; //too high - sigma_max = 850.e+3; //too low - sigma_max = 800.e+3; //IUGG previous test - sigma_max = 1000.e+3; //850 seems small - - if(groundedice<0) sigma_max=200.e+3; + + /*OLD (keep for a little bit)*/ + //sigma_max = 800.e+3; //IUGG previous test + //sigma_max = 1000.e+3; //GRL + //if(groundedice<0) sigma_max=150.e+3; + + /*Tensile stress threshold*/ + if(groundedice<0) + sigma_max = sigma_max_floating; + else + sigma_max = sigma_max_grounded; /*Assign values*/ @@ -294,5 +295,4 @@ IssmDouble calvingratey[NUMVERTICES]; IssmDouble calvingrate[NUMVERTICES]; - /* Get node coordinates and dof list: */ Index: /issm/trunk-jpl/src/c/cores/levelsetfunctionslope_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/levelsetfunctionslope_core.cpp (revision 21625) +++ /issm/trunk-jpl/src/c/cores/levelsetfunctionslope_core.cpp (revision 21626) @@ -33,5 +33,5 @@ } if(domaintype==Domain2DverticalEnum){ - femmodel->parameters->SetParam(LevelsetfunctionSlopeXEnum,InputToExtrudeEnum); + femmodel->parameters->SetParam(LevelsetfunctionSlopeXEnum,InputToExtrudeEnum); extrudefrombase_core(femmodel); } Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21625) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21626) @@ -259,4 +259,7 @@ CalvingratexAverageEnum, CalvingrateyAverageEnum, + CalvingStressThresholdGroundediceEnum, + CalvingStressThresholdFloatingiceEnum, + CalvingMaxEnum, StrainRateparallelEnum, StrainRateperpendicularEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21625) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21626) @@ -265,4 +265,7 @@ case CalvingratexAverageEnum : return "CalvingratexAverage"; case CalvingrateyAverageEnum : return "CalvingrateyAverage"; + case CalvingStressThresholdGroundediceEnum : return "CalvingStressThresholdGroundedice"; + case CalvingStressThresholdFloatingiceEnum : return "CalvingStressThresholdFloatingice"; + case CalvingMaxEnum : return "CalvingMax"; case StrainRateparallelEnum : return "StrainRateparallel"; case StrainRateperpendicularEnum : return "StrainRateperpendicular"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21625) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21626) @@ -271,4 +271,7 @@ else if (strcmp(name,"CalvingratexAverage")==0) return CalvingratexAverageEnum; else if (strcmp(name,"CalvingrateyAverage")==0) return CalvingrateyAverageEnum; + else if (strcmp(name,"CalvingStressThresholdGroundedice")==0) return CalvingStressThresholdGroundediceEnum; + else if (strcmp(name,"CalvingStressThresholdFloatingice")==0) return CalvingStressThresholdFloatingiceEnum; + else if (strcmp(name,"CalvingMax")==0) return CalvingMaxEnum; else if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum; else if (strcmp(name,"StrainRateperpendicular")==0) return StrainRateperpendicularEnum; @@ -380,11 +383,11 @@ else if (strcmp(name,"SmbRequestedOutputs")==0) return SmbRequestedOutputsEnum; else if (strcmp(name,"SmbIsInitialized")==0) return SmbIsInitializedEnum; - else if (strcmp(name,"SmbDzini")==0) return SmbDziniEnum; - else if (strcmp(name,"SmbDini")==0) return SmbDiniEnum; - else if (strcmp(name,"SmbReini")==0) return SmbReiniEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"SmbGdnini")==0) return SmbGdniniEnum; + if (strcmp(name,"SmbDzini")==0) return SmbDziniEnum; + else if (strcmp(name,"SmbDini")==0) return SmbDiniEnum; + else if (strcmp(name,"SmbReini")==0) return SmbReiniEnum; + else if (strcmp(name,"SmbGdnini")==0) return SmbGdniniEnum; else if (strcmp(name,"SmbGspini")==0) return SmbGspiniEnum; else if (strcmp(name,"SmbECini")==0) return SmbECiniEnum; @@ -503,11 +506,11 @@ else if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum; else if (strcmp(name,"SurfaceArea")==0) return SurfaceAreaEnum; - else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum; - else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum; - else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum; + if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum; + else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum; + else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum; + else if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum; else if (strcmp(name,"SurfaceSlopeX")==0) return SurfaceSlopeXEnum; else if (strcmp(name,"SurfaceSlopeY")==0) return SurfaceSlopeYEnum; @@ -626,11 +629,11 @@ else if (strcmp(name,"Outputdefinition31")==0) return Outputdefinition31Enum; else if (strcmp(name,"Outputdefinition32")==0) return Outputdefinition32Enum; - else if (strcmp(name,"Outputdefinition33")==0) return Outputdefinition33Enum; - else if (strcmp(name,"Outputdefinition34")==0) return Outputdefinition34Enum; - else if (strcmp(name,"Outputdefinition35")==0) return Outputdefinition35Enum; else stage=6; } if(stage==6){ - if (strcmp(name,"Outputdefinition36")==0) return Outputdefinition36Enum; + if (strcmp(name,"Outputdefinition33")==0) return Outputdefinition33Enum; + else if (strcmp(name,"Outputdefinition34")==0) return Outputdefinition34Enum; + else if (strcmp(name,"Outputdefinition35")==0) return Outputdefinition35Enum; + else if (strcmp(name,"Outputdefinition36")==0) return Outputdefinition36Enum; else if (strcmp(name,"Outputdefinition37")==0) return Outputdefinition37Enum; else if (strcmp(name,"Outputdefinition38")==0) return Outputdefinition38Enum; @@ -749,11 +752,11 @@ else if (strcmp(name,"Scaled")==0) return ScaledEnum; else if (strcmp(name,"Separate")==0) return SeparateEnum; - else if (strcmp(name,"Sset")==0) return SsetEnum; - else if (strcmp(name,"Dense")==0) return DenseEnum; - else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum; + if (strcmp(name,"Sset")==0) return SsetEnum; + else if (strcmp(name,"Dense")==0) return DenseEnum; + else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum; + else if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum; else if (strcmp(name,"Seq")==0) return SeqEnum; else if (strcmp(name,"Mpi")==0) return MpiEnum; @@ -872,11 +875,11 @@ else if (strcmp(name,"Penta")==0) return PentaEnum; else if (strcmp(name,"PentaInput")==0) return PentaInputEnum; - else if (strcmp(name,"Vertex")==0) return VertexEnum; - else if (strcmp(name,"VertexPId")==0) return VertexPIdEnum; - else if (strcmp(name,"VertexSId")==0) return VertexSIdEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"Option")==0) return OptionEnum; + if (strcmp(name,"Vertex")==0) return VertexEnum; + else if (strcmp(name,"VertexPId")==0) return VertexPIdEnum; + else if (strcmp(name,"VertexSId")==0) return VertexSIdEnum; + else if (strcmp(name,"Option")==0) return OptionEnum; else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum; else if (strcmp(name,"OptionCell")==0) return OptionCellEnum; @@ -995,11 +998,11 @@ else if (strcmp(name,"LACrouzeixRaviart")==0) return LACrouzeixRaviartEnum; else if (strcmp(name,"Closed")==0) return ClosedEnum; - else if (strcmp(name,"Free")==0) return FreeEnum; - else if (strcmp(name,"Open")==0) return OpenEnum; - else if (strcmp(name,"Air")==0) return AirEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"Ice")==0) return IceEnum; + if (strcmp(name,"Free")==0) return FreeEnum; + else if (strcmp(name,"Open")==0) return OpenEnum; + else if (strcmp(name,"Air")==0) return AirEnum; + else if (strcmp(name,"Ice")==0) return IceEnum; else if (strcmp(name,"Melange")==0) return MelangeEnum; else if (strcmp(name,"Water")==0) return WaterEnum; Index: /issm/trunk-jpl/src/m/classes/calvingdev.m =================================================================== --- /issm/trunk-jpl/src/m/classes/calvingdev.m (revision 21625) +++ /issm/trunk-jpl/src/m/classes/calvingdev.m (revision 21626) @@ -6,5 +6,6 @@ classdef calvingdev properties (SetAccess=public) - coeff = NaN; + stress_threshold_groundedice = 0.; + stress_threshold_floatingice = 0.; meltingrate = NaN; end @@ -29,5 +30,4 @@ end % }}} function self = extrude(self,md) % {{{ - self.coeff=project3d(md,'vector',self.coeff,'type','node'); self.meltingrate=project3d(md,'vector',self.meltingrate,'type','node'); end % }}} @@ -35,5 +35,9 @@ %Proportionality coefficient in Pi model - self.coeff=2e13; + %self.coeff=2e13; + + %Default sigma max + self.stress_threshold_groundedice = 1e6; + self.stress_threshold_floatingice = 150e3; end % }}} function md = checkconsistency(self,md,solution,analyses) % {{{ @@ -41,10 +45,12 @@ if (~strcmp(solution,'TransientSolution') | md.transient.ismovingfront==0), return; end - md = checkfield(md,'fieldname','calving.coeff','>',0,'size',[md.mesh.numberofvertices 1]); + md = checkfield(md,'fieldname','calving.stress_threshold_groundedice','>',0,'numel',1,'nan',1,'Inf',1); + md = checkfield(md,'fieldname','calving.stress_threshold_floatingice','>',0,'numel',1,'nan',1,'Inf',1); md = checkfield(md,'fieldname','calving.meltingrate','NaN',1,'Inf',1,'timeseries',1,'>=',0); end % }}} function disp(self) % {{{ disp(sprintf(' Calving Pi parameters:')); - fielddisplay(self,'coeff','proportionality coefficient in Pi model'); + fielddisplay(self,'stress_threshold_groundedice','sigma_max applied to grounded ice only [Pa]'); + fielddisplay(self,'stress_threshold_floatingice','sigma_max applied to floating ice only [Pa]'); fielddisplay(self,'meltingrate','melting rate at given location [m/a]'); @@ -53,6 +59,7 @@ yts=md.constants.yts; WriteData(fid,prefix,'name','md.calving.law','data',2,'format','Integer'); - WriteData(fid,prefix,'object',self,'fieldname','coeff','format','DoubleMat','mattype',1); - WriteData(fid,prefix,'object',self,'fieldname','meltingrate','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1,'scale',1./yts); + WriteData(fid,prefix,'object',self,'fieldname','stress_threshold_groundedice','format','DoubleMat','mattype',1); + WriteData(fid,prefix,'object',self,'fieldname','stress_threshold_floatingice','format','DoubleMat','mattype',1); + WriteData(fid,prefix,'object',self,'fieldname','meltingrate','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'scale',1./yts); end % }}} end Index: /issm/trunk-jpl/src/m/classes/calvinglevermann.m =================================================================== --- /issm/trunk-jpl/src/m/classes/calvinglevermann.m (revision 21625) +++ /issm/trunk-jpl/src/m/classes/calvinglevermann.m (revision 21626) @@ -42,5 +42,5 @@ md = checkfield(md,'fieldname','calving.coeff','>',0,'size',[md.mesh.numberofvertices 1]); - md = checkfield(md,'fieldname','calving.meltingrate','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>=',0); + md = checkfield(md,'fieldname','calving.meltingrate','NaN',1,'Inf',1,'timeseries',1,'>=',0); end % }}} function disp(self) % {{{ @@ -54,5 +54,5 @@ WriteData(fid,prefix,'name','md.calving.law','data',3,'format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','coeff','format','DoubleMat','mattype',1); - WriteData(fid,prefix,'object',self,'fieldname','meltingrate','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1./yts); + WriteData(fid,prefix,'object',self,'fieldname','meltingrate','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'scale',1./yts); end % }}} end Index: /issm/trunk-jpl/src/m/classes/calvingminthickness.m =================================================================== --- /issm/trunk-jpl/src/m/classes/calvingminthickness.m (revision 21625) +++ /issm/trunk-jpl/src/m/classes/calvingminthickness.m (revision 21626) @@ -53,5 +53,5 @@ WriteData(fid,prefix,'name','md.calving.law','data',4,'format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','min_thickness','format','Double'); - WriteData(fid,prefix,'object',self,'fieldname','meltingrate','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1,'scale',1./yts); + WriteData(fid,prefix,'object',self,'fieldname','meltingrate','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'scale',1./yts); end % }}} end Index: /issm/trunk-jpl/src/m/classes/levelset.m =================================================================== --- /issm/trunk-jpl/src/m/classes/levelset.m (revision 21625) +++ /issm/trunk-jpl/src/m/classes/levelset.m (revision 21626) @@ -9,4 +9,5 @@ spclevelset = NaN; reinit_frequency = 5; + calving_max = 0.; end methods @@ -36,6 +37,7 @@ %stabilization = 2 by default - self.stabilization = 2; + self.stabilization = 2; self.reinit_frequency = 5; + self.calving_max = 3000.; end % }}} @@ -46,4 +48,5 @@ md = checkfield(md,'fieldname','levelset.spclevelset','Inf',1,'timeseries',1); md = checkfield(md,'fieldname','levelset.stabilization','values',[0 1 2]); + md = checkfield(md,'fieldname','levelset.calving_max','numel',1,'NaN',1,'Inf',1,'>',0); end % }}} function disp(self) % {{{ @@ -52,9 +55,14 @@ fielddisplay(self,'spclevelset','Levelset constraints (NaN means no constraint)'); fielddisplay(self,'reinit_frequency','Amount of time steps after which the levelset function in re-initialized'); + fielddisplay(self,'calving_max','maximum allowed calving rate (m/a)'); end % }}} function marshall(self,prefix,md,fid) % {{{ + + yts=md.constants.yts; + WriteData(fid,prefix,'object',self,'fieldname','stabilization','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','spclevelset','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); WriteData(fid,prefix,'object',self,'fieldname','reinit_frequency','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','calving_max','format','Double','scale',1./yts); end % }}} function savemodeljs(self,fid,modelname) % {{{ @@ -63,4 +71,5 @@ writejs1Darray(fid,[modelname '.levelset.spclevelset'],self.spclevelset); writejs1Darray(fid,[modelname '.levelset.reinit_frequency'],self.reinit_frequency); + writejsdouble(fid,[modelname '.levelset.calving_max'],self.calving_max); end % }}} Index: /issm/trunk-jpl/src/m/classes/levelset.py =================================================================== --- /issm/trunk-jpl/src/m/classes/levelset.py (revision 21625) +++ /issm/trunk-jpl/src/m/classes/levelset.py (revision 21626) @@ -14,7 +14,8 @@ def __init__(self): # {{{ - self.stabilization = 0 - self.spclevelset = float('NaN') + self.stabilization = 0 + self.spclevelset = float('NaN') self.reinit_frequency = 0 + self.calving_max = 0. #set defaults @@ -27,4 +28,5 @@ string="%s\n%s"%(string,fielddisplay(self,'spclevelset','levelset constraints (NaN means no constraint)')) string="%s\n%s"%(string,fielddisplay(self,'reinit_frequency','Amount of time steps after which the levelset function in re-initialized')) + string="%s\n%s"%(string,fielddisplay(self,'calving_max','maximum allowed calving rate (m/a)')) return string @@ -39,4 +41,5 @@ self.stabilization = 2 self.reinit_frequency = 5 + self.calving_max = 3000 return self @@ -50,4 +53,5 @@ md = checkfield(md,'fieldname','levelset.spclevelset','Inf',1,'timeseries',1) md = checkfield(md,'fieldname','levelset.stabilization','values',[0,1,2]); + md = checkfield(md,'fieldname','levelset.calving_max','numel',1,'NaN',1,'Inf',1,'>',0); return md @@ -55,6 +59,9 @@ def marshall(self,prefix,md,fid): # {{{ + yts=md.constants.yts; + WriteData(fid,prefix,'object',self,'fieldname','stabilization','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','spclevelset','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); WriteData(fid,prefix,'object',self,'fieldname','reinit_frequency','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','calving_max','format','Double','scale',1./yts); # }}}