Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 23606) +++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 23607) @@ -331,17 +331,17 @@ IssmDouble xyz_list[NUMVERTICES][3]; IssmDouble calvingrate[NUMVERTICES]; IssmDouble vx,vy,vel; - IssmDouble critical_fraction,water_height; - IssmDouble bed,Ho,thickness,float_depth; + IssmDouble water_height, bed,Ho,thickness,surface; IssmDouble surface_crevasse[NUMVERTICES], basal_crevasse[NUMVERTICES], crevasse_depth[NUMVERTICES], H_surf, H_surfbasal; - IssmDouble strainparallel, straineffective; + IssmDouble B, strainparallel, straineffective; IssmDouble s_xx,s_xy,s_yy,s1,s2,stmp; - + int crevasse_opening_stress; + /* 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); + /*retrieve the type of crevasse_opening_stress*/ + this->parameters->FindParam(&crevasse_opening_stress,CalvingCrevasseDepthEnum); IssmDouble rho_ice = this->GetMaterialParameter(MaterialsRhoIceEnum); IssmDouble rho_seawater = this->GetMaterialParameter(MaterialsRhoSeawaterEnum); @@ -360,6 +360,7 @@ Input* s_xx_input = inputs->GetInput(DeviatoricStressxxEnum); _assert_(s_xx_input); Input* s_xy_input = inputs->GetInput(DeviatoricStressxyEnum); _assert_(s_xy_input); Input* s_yy_input = inputs->GetInput(DeviatoricStressyyEnum); _assert_(s_yy_input); + Input* B_input = inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(B_input); /*Loop over all elements of this partition*/ GaussTria* gauss=new GaussTria(); @@ -368,7 +369,7 @@ H_input->GetInputValue(&thickness,gauss); bed_input->GetInputValue(&bed,gauss); - surface_input->GetInputValue(&float_depth,gauss); + surface_input->GetInputValue(&surface,gauss); strainrateparallel_input->GetInputValue(&strainparallel,gauss); strainrateeffective_input->GetInputValue(&straineffective,gauss); vx_input->GetInputValue(&vx,gauss); @@ -377,6 +378,7 @@ s_xx_input->GetInputValue(&s_xx,gauss); s_xy_input->GetInputValue(&s_xy,gauss); s_yy_input->GetInputValue(&s_yy,gauss); + B_input->GetInputValue(&B,gauss); vel=sqrt(vx*vx+vy*vy)+1.e-14; @@ -387,28 +389,31 @@ Ho = thickness - (rho_seawater/rho_ice) * (-bed); 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] = s1 / (rho_ice*constant_g); + if(crevasse_opening_stress==0){ /*Otero2010: balance between the tensile deviatoric stress and ice overburden pressure*/ + surface_crevasse[iv] = B * strainparallel * pow(straineffective, ((1 / rheology_n)-1)) / (rho_ice * constant_g); + basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice)) * (B * strainparallel * pow(straineffective,((1/rheology_n)-1)) / (rho_ice*constant_g) - Ho); + } + else if(crevasse_opening_stress==1){ /*Benn2017,Todd2018: maximum principal stress */ + surface_crevasse[iv] = s1 / (rho_ice*constant_g); + basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice))* (s1/ (rho_ice*constant_g)-Ho); + } + + /* some constraints */ if (surface_crevasse[iv]<0.) { surface_crevasse[iv]=0.; water_height = 0.; } + if (basal_crevasse[iv]<0.) basal_crevasse[iv]=0.; + if (bed>0.) basal_crevasse[iv] = 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); + + /* add water in surface crevasse */ + surface_crevasse[iv] = surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height; /* surface crevasse + water */ + crevasse_depth[iv] = surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height + basal_crevasse[iv]; /* surface crevasse + basal crevasse + water */ + } this->inputs->AddInput(new TriaInput(SurfaceCrevasseEnum,&surface_crevasse[0],P1Enum)); Index: ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 23606) +++ ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 23607) @@ -106,7 +106,7 @@ case CalvingHabEnum: break; case CalvingCrevasseDepthEnum: - parameters->AddObject(iomodel->CopyConstantObject("md.calving.critical_fraction",CalvingCrevasseDepthEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.calving.crevasse_opening_stress",CalvingCrevasseDepthEnum)); break; case CalvingDev2Enum: parameters->AddObject(iomodel->CopyConstantObject("md.calving.height_above_floatation",CalvingHeightAboveFloatationEnum)); @@ -670,7 +670,8 @@ /*Intermediaries*/ int calvinglaw; - IssmDouble min_thickness,thickness,hab_fraction,crevassedepth; + IssmDouble min_thickness,thickness,hab_fraction; + IssmDouble crevassedepth,surface_crevasse,surface,critical_fraction; IssmDouble rho_ice,rho_water; IssmDouble bed,water_depth; IssmDouble levelset; @@ -757,16 +758,19 @@ /*Get the DistanceToCalvingfront*/ femmodel->elements->InputDuplicate(MaskIceLevelsetEnum,DistanceToCalvingfrontEnum); femmodel->DistanceToFieldValue(MaskIceLevelsetEnum,0,DistanceToCalvingfrontEnum); - + /*Vector of size number of nodes*/ vec_constraint_nodes=new Vector(femmodel->nodes->NumberOfNodes()); for(int i=0;ielements->Size();i++){ - Element* element = xDynamicCast(femmodel->elements->GetObjectByOffset(i)); - int numnodes = element->GetNumberOfNodes(); - Gauss* gauss = element->NewGauss(); - Input* crevassedepth_input = element->GetInput(CrevasseDepthEnum); _assert_(crevassedepth_input); - Input* bed_input = element->GetInput(BedEnum); _assert_(bed_input); + Element* element = xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + int numnodes = element->GetNumberOfNodes(); + Gauss* gauss = element->NewGauss(); + Input* crevassedepth_input = element->GetInput(CrevasseDepthEnum); _assert_(crevassedepth_input); + Input* bed_input = element->GetInput(BedEnum); _assert_(bed_input); + Input* surface_crevasse_input = element->GetInput(SurfaceCrevasseEnum); _assert_(surface_crevasse_input); + Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input); /*First, look at ice front and figure out if any of the nodes will be calved*/ if(element->IsIcefront()){ @@ -775,7 +779,11 @@ Node* node=element->GetNode(in); crevassedepth_input->GetInputValue(&crevassedepth,gauss); bed_input->GetInputValue(&bed,gauss); - if(crevassedepth>0. && bed<0.){ + surface_crevasse_input->GetInputValue(&surface_crevasse,gauss); + thickness_input->GetInputValue(&thickness,gauss); + surface_input->GetInputValue(&surface,gauss); + + if((surface_crevasse-surface>0. || crevassedepth-thickness>0.) && bed<0.){ vec_constraint_nodes->SetValue(node->Sid(),1.0,INS_VAL); } } @@ -790,12 +798,15 @@ while(nflipped){ local_nflipped=0; for(int i=0;ielements->Size();i++){ - Element* element = xDynamicCast(femmodel->elements->GetObjectByOffset(i)); - int numnodes = element->GetNumberOfNodes(); - Gauss* gauss = element->NewGauss(); - Input* levelset_input = element->GetInput(DistanceToCalvingfrontEnum); _assert_(levelset_input); - Input* crevassedepth_input = element->GetInput(CrevasseDepthEnum); _assert_(crevassedepth_input); - Input* bed_input = element->GetInput(BedEnum); _assert_(bed_input); + Element* element = xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + int numnodes = element->GetNumberOfNodes(); + Gauss* gauss = element->NewGauss(); + Input* levelset_input = element->GetInput(DistanceToCalvingfrontEnum); _assert_(levelset_input); + Input* crevassedepth_input = element->GetInput(CrevasseDepthEnum); _assert_(crevassedepth_input); + Input* bed_input = element->GetInput(BedEnum); _assert_(bed_input); + Input* surface_crevasse_input = element->GetInput(SurfaceCrevasseEnum); _assert_(surface_crevasse_input); + Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input); /*Is this element connected to a node that should be calved*/ bool isconnected = false; @@ -815,8 +826,11 @@ levelset_input->GetInputValue(&levelset,gauss); crevassedepth_input->GetInputValue(&crevassedepth,gauss); bed_input->GetInputValue(&bed,gauss); + surface_crevasse_input->GetInputValue(&surface_crevasse,gauss); + thickness_input->GetInputValue(&thickness,gauss); + surface_input->GetInputValue(&surface,gauss); - if(crevassedepth>0. && bed<0. && levelset>-300. && levelset<0. && constraint_nodes[node->Sid()]==0.){ + if((surface_crevasse-surface>0. || crevassedepth-thickness>0.) && bed<0. && levelset>-300. && levelset<0. && constraint_nodes[node->Sid()]==0.){ local_nflipped++; vec_constraint_nodes->SetValue(node->Sid(),1.0,INS_VAL); }