Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 22180) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 22181) @@ -68,4 +68,11 @@ iomodel->FetchDataToInput(elements,"md.calving.meltingrate",CalvingMeltingrateEnum); break; + case CalvingHabEnum: + iomodel->FetchDataToInput(elements,"md.calving.meltingrate",CalvingMeltingrateEnum); + break; + case CalvingCrevasseDepthEnum: + iomodel->FetchDataToInput(elements,"md.calving.meltingrate",CalvingMeltingrateEnum); + iomodel->FetchDataToInput(elements,"md.calving.water_height",WaterheightEnum); + break; default: _error_("Calving law "<AddObject(iomodel->CopyConstantObject("md.calving.min_thickness",CalvingMinthicknessEnum)); + break; + case CalvingHabEnum: + parameters->AddObject(iomodel->CopyConstantObject("md.calving.flotation_fraction",CalvingHabEnum)); + break; + case CalvingCrevasseDepthEnum: + parameters->AddObject(iomodel->CopyConstantObject("md.calving.critical_fraction",CalvingCrevasseDepthEnum)); break; default: @@ -235,4 +248,15 @@ meltingrate_input = basalelement->GetInput(CalvingMeltingrateEnum); _assert_(meltingrate_input); break; + case CalvingHabEnum: + lsf_slopex_input = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input); + if(dim==2) lsf_slopey_input = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input); + meltingrate_input = basalelement->GetInput(CalvingMeltingrateEnum); _assert_(meltingrate_input); + break; + case CalvingCrevasseDepthEnum: + 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); + meltingrate_input = basalelement->GetInput(CalvingMeltingrateEnum); _assert_(meltingrate_input); + break; default: _error_("Calving law "<GetInputValue(&dlsf[0],gauss); + if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss); + meltingrate_input->GetInputValue(&meltingrate,gauss); + + norm_dlsf=0.; + for(i=0;i1.e-10) + for(i=0;iGetInputValue(&dlsf[0],gauss); + if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss); + 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.; + for(i=0;i1.e-10) + for(i=0;iparameters->FindParam(&calvinglaw,CalvingLawEnum); @@ -640,4 +715,40 @@ } } + + if(calvinglaw==CalvingHabEnum){ + + /*Get the fraction of the flotation thickness at the terminus*/ + femmodel->parameters->FindParam(&flotation_fraction,CalvingHabEnum); + + /*Loop over all elements of this partition*/ + for(int i=0;ielements->Size();i++){ + Element* element = xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + + rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); + rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum); + + int numnodes = element->GetNumberOfNodes(); + Gauss* gauss = element->NewGauss(); + Input* H_input = element->GetInput(ThicknessEnum); _assert_(H_input); + Input* bed = element->GetInput(BedEnum); _assert_(bed); + + /*Potentially constrain nodes of this element*/ + for(int in=0;inGaussNode(element->GetElementType(),in); + Node* node=element->GetNode(in); + H_input->GetInputValue(&thickness,gauss); + bed->GetInputValue(&water_depth,gauss); + if(thickness<((rho_water/rho_ice)*(1+flotation_fraction)*-water_depth)){ + node->ApplyConstraint(0,+1.); + } + else { + /* no ice, set no spc */ + node->DofInFSet(0); + } + } + delete gauss; + } + } + /*Default, do nothing*/ return; Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22180) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22181) @@ -2356,4 +2356,7 @@ this->CalvingRateDev(); break; + case CalvingCrevasseDepthEnum: + this->CalvingCrevasseDepth(); + break; default: _error_("Calving law "<StrainRateparallel(); break; case StrainRateperpendicularEnum: this->StrainRateperpendicular(); break; + case SurfaceCrevasseEnum: this->CalvingCrevasseDepth(); break; + case SigmaVMEnum: this->CalvingRateDev(); break; } Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 22180) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 22181) @@ -189,4 +189,5 @@ virtual void BasalNodeIndices(int* pnumindices,int** pindices,int finiteelement){_error_("not implemented yet");}; virtual void CalvingRateDev(void){_error_("not implemented yet");}; + virtual void CalvingCrevasseDepth(void){_error_("not implemented yet");}; virtual void CalvingRateLevermann(void)=0; virtual IssmDouble CharacteristicLength(void)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 22180) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 22181) @@ -216,5 +216,6 @@ IssmDouble calvingrate[NUMVERTICES]; IssmDouble lambda1,lambda2,ex,ey,vx,vy,vel; - IssmDouble sigma_vm,sigma_max,sigma_max_floating,sigma_max_grounded; + IssmDouble sigma_vm[NUMVERTICES]; + IssmDouble sigma_max,sigma_max_floating,sigma_max_grounded; IssmDouble epse_2,groundedice,bed; @@ -260,5 +261,5 @@ /*Calculate sigma_vm*/ epse_2 = 1./2. *(lambda1*lambda1 + lambda2*lambda2); - sigma_vm = sqrt(3.) * B * pow(epse_2,1./(2.*n)); + sigma_vm[iv] = sqrt(3.) * B * pow(epse_2,1./(2.*n)); /*OLD (keep for a little bit)*/ @@ -279,6 +280,6 @@ } else{ - calvingratex[iv]=vx*sigma_vm/sigma_max; - calvingratey[iv]=vy*sigma_vm/sigma_max; + calvingratex[iv]=vx*sigma_vm[iv]/sigma_max; + calvingratey[iv]=vy*sigma_vm[iv]/sigma_max; } calvingrate[iv] =sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]); @@ -289,6 +290,100 @@ this->inputs->AddInput(new TriaInput(CalvingrateyEnum,&calvingratey[0],P1Enum)); this->inputs->AddInput(new TriaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum)); + this->inputs->AddInput(new TriaInput(SigmaVMEnum,&sigma_vm[0],P1Enum)); /*Clean up and return*/ + delete gauss; +} +/*}}}*/ +void Tria::CalvingCrevasseDepth(){/*{{{*/ + + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble calvingrate[NUMVERTICES]; + IssmDouble vx,vy,vel; + IssmDouble critical_fraction,water_height; + IssmDouble bathymetry,Ho,thickness,float_depth,groundedice; + IssmDouble surface_crevasse[NUMVERTICES], basal_crevasse[NUMVERTICES], crevasse_depth[NUMVERTICES], H_surf, H_surfbasal; + IssmDouble strainparallel, straineffective; + IssmDouble yts; + + /* Get node coordinates and dof list: */ + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + + /*Get the critical fraction of thickness surface and basal crevasses penetrate for calving onset*/ + this->parameters->FindParam(&critical_fraction,CalvingCrevasseDepthEnum); + this->parameters->FindParam(&yts,ConstantsYtsEnum); + + IssmDouble rho_ice = this->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble rho_seawater = this->GetMaterialParameter(MaterialsRhoSeawaterEnum); + IssmDouble rho_freshwater = this->GetMaterialParameter(MaterialsRhoFreshwaterEnum); + IssmDouble constant_g = this->GetMaterialParameter(ConstantsGEnum); + IssmDouble rheology_B = this->GetMaterialParameter(MaterialsRheologyBbarEnum); + IssmDouble rheology_n = this->GetMaterialParameter(MaterialsRheologyNEnum); + + Input* H_input = inputs->GetInput(ThicknessEnum); _assert_(H_input); + Input* bed = inputs->GetInput(BedEnum); _assert_(bed); + Input* surface = inputs->GetInput(SurfaceEnum); _assert_(surface); + Input* strainrateparallel = inputs->GetInput(StrainRateparallelEnum); _assert_(strainrateparallel); + Input* strainrateeffective = inputs->GetInput(StrainRateeffectiveEnum); _assert_(strainrateeffective); + Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = inputs->GetInput(VxEnum); _assert_(vy_input); + Input* gr_input = inputs->GetInput(MaskGroundediceLevelsetEnum); _assert_(gr_input); + Input* waterheight_input = inputs->GetInput(WaterheightEnum); _assert_(waterheight_input); + + /*Loop over all elements of this partition*/ + GaussTria* gauss=new GaussTria(); + for (int iv=0;ivGaussVertex(iv); + + H_input->GetInputValue(&thickness,gauss); + bed->GetInputValue(&bathymetry,gauss); + surface->GetInputValue(&float_depth,gauss); + strainrateparallel->GetInputValue(&strainparallel,gauss); + strainrateeffective->GetInputValue(&straineffective,gauss); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + gr_input->GetInputValue(&groundedice,gauss); + waterheight_input->GetInputValue(&water_height,gauss); + vel=sqrt(vx*vx+vy*vy)+1.e-14; + + Ho = thickness - (rho_seawater/rho_ice) * (-bathymetry); + if(Ho<0.) Ho=0.; + + /*Otero2010: balance between the tensile deviatoric stress and ice overburden pressure*/ + /*surface crevasse*/ + //surface_crevasse[iv] = rheology_B * strainparallel * pow(straineffective, ((1 / rheology_n)-1)) / (rho_ice * constant_g); + surface_crevasse[iv] = 2 * rheology_B * pow(max(strainparallel,0.),(1/rheology_n)) / (rho_ice * constant_g); + if (surface_crevasse[iv]<0.) { + surface_crevasse[iv]=0.; + water_height = 0.; + } + if (surface_crevasse[iv]0.) basal_crevasse[iv] = 0.; + + H_surf = surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height - critical_fraction*float_depth; + H_surfbasal = (surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height + basal_crevasse[iv])-(critical_fraction*thickness); + + crevasse_depth[iv] = max(H_surf,H_surfbasal); + + /*Assign values */ + if(crevasse_depth[iv]>=0. && bathymetry<=0.){ + calvingrate[iv] = vel+3000./yts; + } + else + calvingrate[iv]=0.; + } + + this->inputs->AddInput(new TriaInput(SurfaceCrevasseEnum,&surface_crevasse[0],P1Enum)); + this->inputs->AddInput(new TriaInput(BasalCrevasseEnum,&basal_crevasse[0],P1Enum)); + this->inputs->AddInput(new TriaInput(CrevasseDepthEnum,&crevasse_depth[0],P1Enum)); + this->inputs->AddInput(new TriaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum)); + delete gauss; } @@ -331,8 +426,10 @@ /*Calving rate proportionnal to the positive product of the strain rate along the ice flow direction and the strain rate perpendicular to the ice flow */ - calvingrate[iv]=propcoeff*strainparallel*strainperpendicular; - if(calvingrate[iv]<0. || bed>0.){ + if(strainparallel>0. && strainperpendicular>0. && bed<=0.){ + calvingrate[iv]=propcoeff*strainparallel*strainperpendicular; + } + else calvingrate[iv]=0.; - } + calvingratex[iv]=calvingrate[iv]*vx/(sqrt(vel)+1.e-14); calvingratey[iv]=calvingrate[iv]*vy/(sqrt(vel)+1.e-14); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 22180) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 22181) @@ -53,4 +53,5 @@ void AverageOntoPartition(Vector* partition_contributions,Vector* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part); void CalvingRateDev(); + void CalvingCrevasseDepth(); void CalvingRateLevermann(); IssmDouble CharacteristicLength(void); Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 22180) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 22181) @@ -1940,4 +1940,12 @@ } /*}}}*/ +void FemModel::StrainRateeffectivex(){/*{{{*/ + + for(int i=0;iSize();i++){ + Element* element=dynamic_cast(this->elements->GetObjectByOffset(i)); + element->ComputeStrainRate(); + } +} +/*}}}*/ void FemModel::StressIntensityFactorx(){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 22180) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 22181) @@ -110,4 +110,5 @@ void StrainRateparallelx(); void StrainRateperpendicularx(); + void StrainRateeffectivex(); void StressIntensityFactorx(); void TotalFloatingBmbx(IssmDouble* pFbmb); Index: /issm/trunk-jpl/src/c/modules/Calvingx/Calvingx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Calvingx/Calvingx.cpp (revision 22180) +++ /issm/trunk-jpl/src/c/modules/Calvingx/Calvingx.cpp (revision 22181) @@ -18,4 +18,11 @@ case CalvingMinthicknessEnum: break; + case CalvingHabEnum: + break; + case CalvingCrevasseDepthEnum: + femmodel->StrainRateparallelx(); + femmodel->StrainRateeffectivex(); + femmodel->ElementOperationx(&Element::CalvingCrevasseDepth); + break; case CalvingLevermannEnum: if(VerboseModule()) _printf0_(" computing Levermann's calving rate\n"); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22180) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22181) @@ -258,4 +258,10 @@ CalvingDevEnum, CalvingMinthicknessEnum, + CalvingHabEnum, + CalvingCrevasseDepthEnum, + SurfaceCrevasseEnum, + BasalCrevasseEnum, + CrevasseDepthEnum, + WaterheightEnum, DefaultCalvingEnum, CalvinglevermannCoeffEnum, @@ -269,4 +275,5 @@ CalvingStressThresholdFloatingiceEnum, CalvingMaxEnum, + SigmaVMEnum, StrainRateparallelEnum, StrainRateperpendicularEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22180) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22181) @@ -264,4 +264,10 @@ case CalvingDevEnum : return "CalvingDev"; case CalvingMinthicknessEnum : return "CalvingMinthickness"; + case CalvingHabEnum : return "CalvingHab"; + case CalvingCrevasseDepthEnum : return "CalvingCrevasseDepth"; + case SurfaceCrevasseEnum : return "SurfaceCrevasse"; + case BasalCrevasseEnum : return "BasalCrevasse"; + case CrevasseDepthEnum : return "CrevasseDepth"; + case WaterheightEnum : return "Waterheight"; case DefaultCalvingEnum : return "DefaultCalving"; case CalvinglevermannCoeffEnum : return "CalvinglevermannCoeff"; @@ -275,4 +281,5 @@ case CalvingStressThresholdFloatingiceEnum : return "CalvingStressThresholdFloatingice"; case CalvingMaxEnum : return "CalvingMax"; + case SigmaVMEnum : return "SigmaVM"; 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 22180) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 22181) @@ -270,4 +270,13 @@ else if (strcmp(name,"CalvingDev")==0) return CalvingDevEnum; else if (strcmp(name,"CalvingMinthickness")==0) return CalvingMinthicknessEnum; + else stage=3; + } + if(stage==3){ + if (strcmp(name,"CalvingHab")==0) return CalvingHabEnum; + else if (strcmp(name,"CalvingCrevasseDepth")==0) return CalvingCrevasseDepthEnum; + else if (strcmp(name,"SurfaceCrevasse")==0) return SurfaceCrevasseEnum; + else if (strcmp(name,"BasalCrevasse")==0) return BasalCrevasseEnum; + else if (strcmp(name,"CrevasseDepth")==0) return CrevasseDepthEnum; + else if (strcmp(name,"Waterheight")==0) return WaterheightEnum; else if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum; else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum; @@ -281,4 +290,5 @@ else if (strcmp(name,"CalvingStressThresholdFloatingice")==0) return CalvingStressThresholdFloatingiceEnum; else if (strcmp(name,"CalvingMax")==0) return CalvingMaxEnum; + else if (strcmp(name,"SigmaVM")==0) return SigmaVMEnum; else if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum; else if (strcmp(name,"StrainRateperpendicular")==0) return StrainRateperpendicularEnum; @@ -757,5 +767,8 @@ if (strcmp(name,"RootPath")==0) return RootPathEnum; else if (strcmp(name,"OutputFileName")==0) return OutputFileNameEnum; - else if (strcmp(name,"InputFileName")==0) return InputFileNameEnum; + else stage=7; + } + if(stage==7){ + if (strcmp(name,"InputFileName")==0) return InputFileNameEnum; else if (strcmp(name,"LockFileName")==0) return LockFileNameEnum; else if (strcmp(name,"RestartFileName")==0) return RestartFileNameEnum; @@ -919,5 +932,8 @@ else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum; else if (strcmp(name,"ElementSId")==0) return ElementSIdEnum; - else if (strcmp(name,"VectorParam")==0) return VectorParamEnum; + else stage=8; + } + if(stage==8){ + if (strcmp(name,"VectorParam")==0) return VectorParamEnum; else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum; else if (strcmp(name,"Segment")==0) return SegmentEnum; @@ -1041,5 +1057,8 @@ else if (strcmp(name,"P1bubblecondensed")==0) return P1bubblecondensedEnum; else if (strcmp(name,"P2")==0) return P2Enum; - else if (strcmp(name,"P2bubble")==0) return P2bubbleEnum; + else stage=9; + } + if(stage==9){ + if (strcmp(name,"P2bubble")==0) return P2bubbleEnum; else if (strcmp(name,"P2bubblecondensed")==0) return P2bubblecondensedEnum; else if (strcmp(name,"P2xP1")==0) return P2xP1Enum; Index: /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 22180) +++ /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 22181) @@ -151,4 +151,6 @@ case 3: return CalvingLevermannEnum; case 4: return CalvingMinthicknessEnum; + case 5: return CalvingHabEnum; + case 6: return CalvingCrevasseDepthEnum; default: _error_("Marshalled Calving law code \""<