Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 22487) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 22488) @@ -264,5 +264,4 @@ 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; @@ -386,10 +385,6 @@ lsf_slopex_input->GetInputValue(&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.; @@ -400,5 +395,5 @@ if(norm_dlsf>1.e-10) for(i=0;iparameters->FindParam(&calvinglaw,CalvingLawEnum); @@ -793,5 +788,5 @@ Gauss* gauss = element->NewGauss(); Input* H_input = element->GetInput(ThicknessEnum); _assert_(H_input); - Input* bed = element->GetInput(BedEnum); _assert_(bed); + Input* bed_input = element->GetInput(BedEnum); _assert_(bed_input); Input* hab_fraction_input = element->GetInput(CalvingHabFractionEnum); _assert_(hab_fraction_input); Input* ls_input = element->GetInput(DistanceToCalvingfrontEnum); _assert_(ls_input); @@ -802,9 +797,9 @@ Node* node=element->GetNode(in); H_input->GetInputValue(&thickness,gauss); - bed->GetInputValue(&water_depth,gauss); + bed_input->GetInputValue(&water_depth,gauss); ls_input->GetInputValue(&levelset,gauss); hab_fraction_input->GetInputValue(&hab_fraction,gauss); - if(thickness<((rho_water/rho_ice)*(1+hab_fraction)*-water_depth) && levelset>-300){ + if(thickness<((rho_water/rho_ice)*(1+hab_fraction)*-water_depth) && levelset>-300. && levelset<0.){ node->ApplyConstraint(0,+1.); } @@ -818,4 +813,115 @@ } + if(calvinglaw==CalvingCrevasseDepthEnum){ + + int nflipped,local_nflipped; + Vector* vec_constraint_nodes = NULL; + IssmDouble* constraint_nodes = NULL; + + /*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(LevelsetAnalysisEnum)); + + 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); + + /*First, look at ice front and figure out if any of the nodes will be calved*/ + if(element->IsIcefront()){ + for(int in=0;inGaussNode(element->GetElementType(),in); + Node* node=element->GetNode(in); + crevassedepth_input->GetInputValue(&crevassedepth,gauss); + bed_input->GetInputValue(&bed,gauss); + if(crevassedepth>0. && bed<0.){ + vec_constraint_nodes->SetValue(node->Sid(),1.0,INS_VAL); + } + } + } + delete gauss; + } + /*Assemble vector and serialize: */ + vec_constraint_nodes->Assemble(); + constraint_nodes=vec_constraint_nodes->ToMPISerial(); + + nflipped=1; + 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); + + /*Is this element connected to a node that should be calved*/ + bool isconnected = false; + for(int in=0;inGetNode(in); + if(constraint_nodes[node->Sid()]==1.){ + isconnected = true; + break; + } + } + + /*Check status if connected*/ + if(isconnected){ + for(int in=0;inGaussNode(element->GetElementType(),in); + Node* node=element->GetNode(in); + levelset_input->GetInputValue(&levelset,gauss); + crevassedepth_input->GetInputValue(&crevassedepth,gauss); + bed_input->GetInputValue(&bed,gauss); + + if(crevassedepth>0. && bed<0. && levelset>-300. && levelset<0. && constraint_nodes[node->Sid()]==0.){ + local_nflipped++; + vec_constraint_nodes->SetValue(node->Sid(),1.0,INS_VAL); + } + } + } + } + + /*Count how many new nodes were found*/ + ISSM_MPI_Allreduce(&local_nflipped,&nflipped,1,ISSM_MPI_INT,ISSM_MPI_SUM,IssmComm::GetComm()); + //_printf0_("Found "<Assemble(); + xDelete(constraint_nodes); + constraint_nodes=vec_constraint_nodes->ToMPISerial(); + } + /*Free ressources:*/ + delete vec_constraint_nodes; + + /*Contrain the nodes that will be calved*/ + for(int i=0;ielements->Size();i++){ + Element* element = xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + int numnodes = element->GetNumberOfNodes(); + Gauss* gauss = element->NewGauss(); + /*Potentially constrain nodes of this element*/ + for(int in=0;inGaussNode(element->GetElementType(),in); + Node* node=element->GetNode(in); + + if(constraint_nodes[node->Sid()]){ + node->ApplyConstraint(0,+1.); + } + else { + /* no ice, set no spc */ + node->DofInFSet(0); + } + } + delete gauss; + } + xDelete(constraint_nodes); + } + /*Default, do nothing*/ return; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 22487) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 22488) @@ -302,9 +302,9 @@ IssmDouble calvingrate[NUMVERTICES]; IssmDouble vx,vy,vel; - IssmDouble rheology_B,critical_fraction,water_height; - IssmDouble bathymetry,Ho,thickness,float_depth,groundedice; + IssmDouble critical_fraction,water_height; + IssmDouble bed,Ho,thickness,float_depth; IssmDouble surface_crevasse[NUMVERTICES], basal_crevasse[NUMVERTICES], crevasse_depth[NUMVERTICES], H_surf, H_surfbasal; IssmDouble strainparallel, straineffective; - IssmDouble yts; + IssmDouble s_xx,s_xy,s_yy,s1,s2,stmp; /* Get node coordinates and dof list: */ @@ -313,5 +313,4 @@ /*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); @@ -321,15 +320,16 @@ IssmDouble rheology_n = this->GetMaterialParameter(MaterialsRheologyNEnum); - Input* B_input = inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(B_input); - 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); - + Input* H_input = inputs->GetInput(ThicknessEnum); _assert_(H_input); + Input* bed_input = inputs->GetInput(BedEnum); _assert_(bed_input); + Input* surface_input = inputs->GetInput(SurfaceEnum); _assert_(surface_input); + Input* strainrateparallel_input = inputs->GetInput(StrainRateparallelEnum); _assert_(strainrateparallel_input); + Input* strainrateeffective_input = inputs->GetInput(StrainRateeffectiveEnum); _assert_(strainrateeffective_input); + Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = inputs->GetInput(VxEnum); _assert_(vy_input); + Input* waterheight_input = inputs->GetInput(WaterheightEnum); _assert_(waterheight_input); + 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); + /*Loop over all elements of this partition*/ GaussTria* gauss=new GaussTria(); @@ -337,36 +337,42 @@ gauss->GaussVertex(iv); - B_input->GetInputValue(&rheology_B,gauss); H_input->GetInputValue(&thickness,gauss); - bed->GetInputValue(&bathymetry,gauss); - surface->GetInputValue(&float_depth,gauss); - strainrateparallel->GetInputValue(&strainparallel,gauss); - strainrateeffective->GetInputValue(&straineffective,gauss); + bed_input->GetInputValue(&bed,gauss); + surface_input->GetInputValue(&float_depth,gauss); + strainrateparallel_input->GetInputValue(&strainparallel,gauss); + strainrateeffective_input->GetInputValue(&straineffective,gauss); vx_input->GetInputValue(&vx,gauss); vy_input->GetInputValue(&vy,gauss); - gr_input->GetInputValue(&groundedice,gauss); waterheight_input->GetInputValue(&water_height,gauss); + s_xx_input->GetInputValue(&s_xx,gauss); + s_xy_input->GetInputValue(&s_xy,gauss); + s_yy_input->GetInputValue(&s_yy,gauss); + vel=sqrt(vx*vx+vy*vy)+1.e-14; - Ho = thickness - (rho_seawater/rho_ice) * (-bathymetry); + s1=(s_xx+s_yy)/2.+sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2)); + s2=(s_xx+s_yy)/2.-sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2)); + if(fabs(s2)>fabs(s1)){stmp=s2; s2=s1; s1=stmp;} + + 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] = 2 * rheology_B * pow(max(strainparallel,0.),(1/rheology_n)) / (rho_ice * constant_g); + //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 (surface_crevasse[iv]<0.) { surface_crevasse[iv]=0.; water_height = 0.; } - if (surface_crevasse[iv]0.) basal_crevasse[iv] = 0.; + if (bed>0.) basal_crevasse[iv] = 0.; H_surf = surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height - critical_fraction*float_depth; @@ -374,17 +380,9 @@ 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; Index: /issm/trunk-jpl/src/c/modules/Calvingx/Calvingx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Calvingx/Calvingx.cpp (revision 22487) +++ /issm/trunk-jpl/src/c/modules/Calvingx/Calvingx.cpp (revision 22488) @@ -23,4 +23,5 @@ femmodel->StrainRateparallelx(); femmodel->StrainRateeffectivex(); + femmodel->DeviatoricStressx(); femmodel->ElementOperationx(&Element::CalvingCrevasseDepth); break;