Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 22198) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 22199) @@ -75,4 +75,8 @@ iomodel->FetchDataToInput(elements,"md.calving.water_height",WaterheightEnum); break; + case CalvingDev2Enum: + iomodel->FetchDataToInput(elements,"md.calving.stress_threshold_groundedice",CalvingStressThresholdGroundediceEnum); + iomodel->FetchDataToInput(elements,"md.calving.stress_threshold_floatingice",CalvingStressThresholdFloatingiceEnum); + break; default: _error_("Calving law "<AddObject(iomodel->CopyConstantObject("md.calving.critical_fraction",CalvingCrevasseDepthEnum)); break; + case CalvingDev2Enum: + parameters->AddObject(iomodel->CopyConstantObject("md.calving.height_above_floatation",CalvingHeightAboveFloatationEnum)); + break; default: _error_("Calving law "<GetInput(CalvingMeltingrateEnum); _assert_(meltingrate_input); break; + case CalvingDev2Enum: + basalelement->FindParam(&calvinghaf,CalvingHeightAboveFloatationEnum); + lsf_slopex_input = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input); + if(dim==2) lsf_slopey_input = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input); + calvingrate_input = basalelement->GetInput(CalvingCalvingrateEnum); _assert_(calvingrate_input); + gr_input=basalelement->GetInput(MaskGroundediceLevelsetEnum); _assert_(gr_input); + break; default: _error_("Calving law "<GetInputValue(&c[0],gauss); @@ -394,4 +408,43 @@ } break; + + case CalvingDev2Enum: + { + lsf_slopex_input->GetInputValue(&dlsf[0],gauss); + if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss); + calvingrate_input->GetInputValue(&calvingrate,gauss); + gr_input->GetInputValue(&groundedice,gauss); + + //idea: no retreat on ice above critical calving height "calvinghaf" . Limit using regularized Heaviside function. + vel=sqrt(v[0]*v[0] + v[1]*v[1]); + haf_eps=10.; + if(groundedice-calvinghaf<=-haf_eps){ + // ice floats freely below calvinghaf: calve freely + } + else if(groundedice-calvinghaf>=haf_eps){ + // ice is well above calvinghaf -> no calving back, i.e. limit calving rate to ice velocity + calvingrate=min(calvingrate,vel); + } + else{ // ice is close to calvinghaf: smooth transition between limitation and free calving. + heaviside=(groundedice-calvinghaf+haf_eps)/(2.*haf_eps) + sin(PI*(groundedice-calvinghaf)/haf_eps)/(2.*PI); + calvingrate=heaviside*(min(calvingrate,vel)-calvingrate)+calvingrate; + } + + norm_dlsf=0.; + for(i=0;i1.e-10) + for(i=0;iCalvingRateVonmises(); break; Index: /issm/trunk-jpl/src/c/modules/Calvingx/Calvingx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Calvingx/Calvingx.cpp (revision 22198) +++ /issm/trunk-jpl/src/c/modules/Calvingx/Calvingx.cpp (revision 22199) @@ -32,4 +32,5 @@ break; case CalvingVonmisesEnum: + case CalvingDev2Enum: femmodel->CalvingRateVonmisesx(); femmodel->ElementOperationx(&Element::CalvingRateVonmises); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22198) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22199) @@ -256,4 +256,5 @@ CalvingCalvingrateEnum, CalvingMeltingrateEnum, + CalvingHeightAboveFloatationEnum, CalvingLevermannEnum, CalvingVonmisesEnum, @@ -261,4 +262,5 @@ CalvingHabEnum, CalvingCrevasseDepthEnum, + CalvingDev2Enum, SurfaceCrevasseEnum, BasalCrevasseEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22198) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22199) @@ -262,4 +262,5 @@ case CalvingCalvingrateEnum : return "CalvingCalvingrate"; case CalvingMeltingrateEnum : return "CalvingMeltingrate"; + case CalvingHeightAboveFloatationEnum : return "CalvingHeightAboveFloatation"; case CalvingLevermannEnum : return "CalvingLevermann"; case CalvingVonmisesEnum : return "CalvingVonmises"; @@ -267,4 +268,5 @@ case CalvingHabEnum : return "CalvingHab"; case CalvingCrevasseDepthEnum : return "CalvingCrevasseDepth"; + case CalvingDev2Enum : return "CalvingDev2"; case SurfaceCrevasseEnum : return "SurfaceCrevasse"; case BasalCrevasseEnum : return "BasalCrevasse"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 22198) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 22199) @@ -268,4 +268,5 @@ else if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum; else if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum; + else if (strcmp(name,"CalvingHeightAboveFloatation")==0) return CalvingHeightAboveFloatationEnum; else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum; else if (strcmp(name,"CalvingVonmises")==0) return CalvingVonmisesEnum; @@ -273,4 +274,5 @@ else if (strcmp(name,"CalvingHab")==0) return CalvingHabEnum; else if (strcmp(name,"CalvingCrevasseDepth")==0) return CalvingCrevasseDepthEnum; + else if (strcmp(name,"CalvingDev2")==0) return CalvingDev2Enum; else if (strcmp(name,"SurfaceCrevasse")==0) return SurfaceCrevasseEnum; else if (strcmp(name,"BasalCrevasse")==0) return BasalCrevasseEnum; @@ -381,10 +383,10 @@ else if (strcmp(name,"TransientIsthermal")==0) return TransientIsthermalEnum; else if (strcmp(name,"TransientIsgia")==0) return TransientIsgiaEnum; - else if (strcmp(name,"TransientIsesa")==0) return TransientIsesaEnum; - else if (strcmp(name,"TransientIsdamageevolution")==0) return TransientIsdamageevolutionEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"TransientIshydrology")==0) return TransientIshydrologyEnum; + if (strcmp(name,"TransientIsesa")==0) return TransientIsesaEnum; + else if (strcmp(name,"TransientIsdamageevolution")==0) return TransientIsdamageevolutionEnum; + else if (strcmp(name,"TransientIshydrology")==0) return TransientIshydrologyEnum; else if (strcmp(name,"TransientIsmovingfront")==0) return TransientIsmovingfrontEnum; else if (strcmp(name,"TransientIsslr")==0) return TransientIsslrEnum; @@ -504,10 +506,10 @@ else if (strcmp(name,"SmbBMin")==0) return SmbBMinEnum; else if (strcmp(name,"Adjointp")==0) return AdjointpEnum; - else if (strcmp(name,"Adjointx")==0) return AdjointxEnum; - else if (strcmp(name,"Adjointy")==0) return AdjointyEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"Adjointz")==0) return AdjointzEnum; + if (strcmp(name,"Adjointx")==0) return AdjointxEnum; + else if (strcmp(name,"Adjointy")==0) return AdjointyEnum; + else if (strcmp(name,"Adjointz")==0) return AdjointzEnum; else if (strcmp(name,"BalancethicknessMisfit")==0) return BalancethicknessMisfitEnum; else if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum; @@ -627,10 +629,10 @@ else if (strcmp(name,"Outputdefinition10")==0) return Outputdefinition10Enum; else if (strcmp(name,"Outputdefinition11")==0) return Outputdefinition11Enum; - else if (strcmp(name,"Outputdefinition12")==0) return Outputdefinition12Enum; - else if (strcmp(name,"Outputdefinition13")==0) return Outputdefinition13Enum; else stage=6; } if(stage==6){ - if (strcmp(name,"Outputdefinition14")==0) return Outputdefinition14Enum; + if (strcmp(name,"Outputdefinition12")==0) return Outputdefinition12Enum; + else if (strcmp(name,"Outputdefinition13")==0) return Outputdefinition13Enum; + else if (strcmp(name,"Outputdefinition14")==0) return Outputdefinition14Enum; else if (strcmp(name,"Outputdefinition15")==0) return Outputdefinition15Enum; else if (strcmp(name,"Outputdefinition16")==0) return Outputdefinition16Enum; @@ -750,10 +752,10 @@ else if (strcmp(name,"Gradient")==0) return GradientEnum; else if (strcmp(name,"GroundinglineMigration")==0) return GroundinglineMigrationEnum; - else if (strcmp(name,"Gset")==0) return GsetEnum; - else if (strcmp(name,"Index")==0) return IndexEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"Indexed")==0) return IndexedEnum; + if (strcmp(name,"Gset")==0) return GsetEnum; + else if (strcmp(name,"Index")==0) return IndexEnum; + else if (strcmp(name,"Indexed")==0) return IndexedEnum; else if (strcmp(name,"Intersect")==0) return IntersectEnum; else if (strcmp(name,"Nodal")==0) return NodalEnum; @@ -873,10 +875,10 @@ else if (strcmp(name,"AmrIceFrontResolution")==0) return AmrIceFrontResolutionEnum; else if (strcmp(name,"AmrIceFrontDistance")==0) return AmrIceFrontDistanceEnum; - else if (strcmp(name,"AmrThicknessErrorResolution")==0) return AmrThicknessErrorResolutionEnum; - else if (strcmp(name,"AmrThicknessErrorThreshold")==0) return AmrThicknessErrorThresholdEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"AmrDeviatoricErrorResolution")==0) return AmrDeviatoricErrorResolutionEnum; + if (strcmp(name,"AmrThicknessErrorResolution")==0) return AmrThicknessErrorResolutionEnum; + else if (strcmp(name,"AmrThicknessErrorThreshold")==0) return AmrThicknessErrorThresholdEnum; + else if (strcmp(name,"AmrDeviatoricErrorResolution")==0) return AmrDeviatoricErrorResolutionEnum; else if (strcmp(name,"AmrDeviatoricErrorThreshold")==0) return AmrDeviatoricErrorThresholdEnum; else if (strcmp(name,"DeviatoricStressErrorEstimator")==0) return DeviatoricStressErrorEstimatorEnum; @@ -996,10 +998,10 @@ else if (strcmp(name,"GiaIvinsAnalysis")==0) return GiaIvinsAnalysisEnum; else if (strcmp(name,"EsaSolution")==0) return EsaSolutionEnum; - else if (strcmp(name,"EsaAnalysis")==0) return EsaAnalysisEnum; - else if (strcmp(name,"LoveSolution")==0) return LoveSolutionEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"LoveAnalysis")==0) return LoveAnalysisEnum; + if (strcmp(name,"EsaAnalysis")==0) return EsaAnalysisEnum; + else if (strcmp(name,"LoveSolution")==0) return LoveSolutionEnum; + else if (strcmp(name,"LoveAnalysis")==0) return LoveAnalysisEnum; else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum; else if (strcmp(name,"ExtrapolationAnalysis")==0) return ExtrapolationAnalysisEnum; Index: /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 22198) +++ /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 22199) @@ -153,4 +153,5 @@ case 5: return CalvingHabEnum; case 6: return CalvingCrevasseDepthEnum; + case 7: return CalvingDev2Enum; default: _error_("Marshalled Calving law code \""<',0,'nan',1,'Inf',1); + md = checkfield(md,'fieldname','calving.stress_threshold_floatingice','>',0,'nan',1,'Inf',1); + md = checkfield(md,'fieldname','calving.height_above_floatation','<=',0,'singletimeseries',1); + + end % }}} + function disp(self) % {{{ + disp(sprintf(' Calving VonMises parameters:')); + 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,'height_above_floatation','height above floatation that the ice tongue must reach before it can calve [m]'); + + end % }}} + function marshall(self,prefix,md,fid) % {{{ + yts=md.constants.yts; + WriteData(fid,prefix,'name','md.calving.law','data',7,'format','Integer'); + 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','height_above_floatation', 'format','Double'); + end % }}} + end +end