Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 23606) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 23607) @@ -107,5 +107,5 @@ 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: @@ -671,5 +671,6 @@ /*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; @@ -758,14 +759,17 @@ 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*/ @@ -776,5 +780,9 @@ 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); } @@ -791,10 +799,13 @@ 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*/ @@ -816,6 +827,9 @@ crevassedepth_input->GetInputValue(&crevassedepth,gauss); bed_input->GetInputValue(&bed,gauss); - - if(crevassedepth>0. && bed<0. && levelset>-300. && levelset<0. && constraint_nodes[node->Sid()]==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. && levelset>-300. && levelset<0. && constraint_nodes[node->Sid()]==0.){ local_nflipped++; vec_constraint_nodes->SetValue(node->Sid(),1.0,INS_VAL); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 23606) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 23607) @@ -332,15 +332,15 @@ 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); @@ -361,4 +361,5 @@ 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*/ @@ -369,5 +370,5 @@ 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); @@ -378,4 +379,5 @@ 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; @@ -388,26 +390,29 @@ 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 */ + }