![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
source:
issm/oecreview/Archive/23390-24306/ISSM-23606-23607.diff
| Last change on this file was 24307, checked in by , 5 years ago | |
|---|---|
| File size: 10.2 KB | |
-
../trunk-jpl/src/c/classes/Elements/Tria.cpp
331 331 IssmDouble xyz_list[NUMVERTICES][3]; 332 332 IssmDouble calvingrate[NUMVERTICES]; 333 333 IssmDouble vx,vy,vel; 334 IssmDouble critical_fraction,water_height; 335 IssmDouble bed,Ho,thickness,float_depth; 334 IssmDouble water_height, bed,Ho,thickness,surface; 336 335 IssmDouble surface_crevasse[NUMVERTICES], basal_crevasse[NUMVERTICES], crevasse_depth[NUMVERTICES], H_surf, H_surfbasal; 337 IssmDouble strainparallel, straineffective;336 IssmDouble B, strainparallel, straineffective; 338 337 IssmDouble s_xx,s_xy,s_yy,s1,s2,stmp; 339 338 int crevasse_opening_stress; 339 340 340 /* Get node coordinates and dof list: */ 341 341 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); 342 342 343 /* Get the critical fraction of thickness surface and basal crevasses penetrate for calving onset*/344 this->parameters->FindParam(&cr itical_fraction,CalvingCrevasseDepthEnum);343 /*retrieve the type of crevasse_opening_stress*/ 344 this->parameters->FindParam(&crevasse_opening_stress,CalvingCrevasseDepthEnum); 345 345 346 346 IssmDouble rho_ice = this->GetMaterialParameter(MaterialsRhoIceEnum); 347 347 IssmDouble rho_seawater = this->GetMaterialParameter(MaterialsRhoSeawaterEnum); … … 360 360 Input* s_xx_input = inputs->GetInput(DeviatoricStressxxEnum); _assert_(s_xx_input); 361 361 Input* s_xy_input = inputs->GetInput(DeviatoricStressxyEnum); _assert_(s_xy_input); 362 362 Input* s_yy_input = inputs->GetInput(DeviatoricStressyyEnum); _assert_(s_yy_input); 363 Input* B_input = inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(B_input); 363 364 364 365 /*Loop over all elements of this partition*/ 365 366 GaussTria* gauss=new GaussTria(); … … 368 369 369 370 H_input->GetInputValue(&thickness,gauss); 370 371 bed_input->GetInputValue(&bed,gauss); 371 surface_input->GetInputValue(& float_depth,gauss);372 surface_input->GetInputValue(&surface,gauss); 372 373 strainrateparallel_input->GetInputValue(&strainparallel,gauss); 373 374 strainrateeffective_input->GetInputValue(&straineffective,gauss); 374 375 vx_input->GetInputValue(&vx,gauss); … … 377 378 s_xx_input->GetInputValue(&s_xx,gauss); 378 379 s_xy_input->GetInputValue(&s_xy,gauss); 379 380 s_yy_input->GetInputValue(&s_yy,gauss); 381 B_input->GetInputValue(&B,gauss); 380 382 381 383 vel=sqrt(vx*vx+vy*vy)+1.e-14; 382 384 … … 387 389 Ho = thickness - (rho_seawater/rho_ice) * (-bed); 388 390 if(Ho<0.) Ho=0.; 389 391 390 /*Otero2010: balance between the tensile deviatoric stress and ice overburden pressure*/ 391 /*surface crevasse*/ 392 //surface_crevasse[iv] = rheology_B * strainparallel * pow(straineffective, ((1 / rheology_n)-1)) / (rho_ice * constant_g); 393 surface_crevasse[iv] = s1 / (rho_ice*constant_g); 392 if(crevasse_opening_stress==0){ /*Otero2010: balance between the tensile deviatoric stress and ice overburden pressure*/ 393 surface_crevasse[iv] = B * strainparallel * pow(straineffective, ((1 / rheology_n)-1)) / (rho_ice * constant_g); 394 basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice)) * (B * strainparallel * pow(straineffective,((1/rheology_n)-1)) / (rho_ice*constant_g) - Ho); 395 } 396 else if(crevasse_opening_stress==1){ /*Benn2017,Todd2018: maximum principal stress */ 397 surface_crevasse[iv] = s1 / (rho_ice*constant_g); 398 basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice))* (s1/ (rho_ice*constant_g)-Ho); 399 } 400 401 /* some constraints */ 394 402 if (surface_crevasse[iv]<0.) { 395 403 surface_crevasse[iv]=0.; 396 404 water_height = 0.; 397 405 } 406 if (basal_crevasse[iv]<0.) basal_crevasse[iv]=0.; 407 if (bed>0.) basal_crevasse[iv] = 0.; 408 398 409 //if (surface_crevasse[iv]<water_height){ 399 410 // water_height = surface_crevasse[iv]; 400 411 //} 401 402 /*basal crevasse*/ 403 //basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice)) * (rheology_B * strainparallel * pow(straineffective,((1/rheology_n)-1)) / (rho_ice*constant_g) - Ho); 404 basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice))* (s1/ (rho_ice*constant_g)-Ho); 405 if (basal_crevasse[iv]<0.) basal_crevasse[iv]=0.; 406 if (bed>0.) basal_crevasse[iv] = 0.; 407 408 H_surf = surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height - critical_fraction*float_depth; 409 H_surfbasal = (surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height + basal_crevasse[iv])-(critical_fraction*thickness); 410 411 crevasse_depth[iv] = max(H_surf,H_surfbasal); 412 413 /* add water in surface crevasse */ 414 surface_crevasse[iv] = surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height; /* surface crevasse + water */ 415 crevasse_depth[iv] = surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height + basal_crevasse[iv]; /* surface crevasse + basal crevasse + water */ 416 412 417 } 413 418 414 419 this->inputs->AddInput(new TriaInput(SurfaceCrevasseEnum,&surface_crevasse[0],P1Enum)); -
../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
106 106 case CalvingHabEnum: 107 107 break; 108 108 case CalvingCrevasseDepthEnum: 109 parameters->AddObject(iomodel->CopyConstantObject("md.calving.cr itical_fraction",CalvingCrevasseDepthEnum));109 parameters->AddObject(iomodel->CopyConstantObject("md.calving.crevasse_opening_stress",CalvingCrevasseDepthEnum)); 110 110 break; 111 111 case CalvingDev2Enum: 112 112 parameters->AddObject(iomodel->CopyConstantObject("md.calving.height_above_floatation",CalvingHeightAboveFloatationEnum)); … … 670 670 671 671 /*Intermediaries*/ 672 672 int calvinglaw; 673 IssmDouble min_thickness,thickness,hab_fraction,crevassedepth; 673 IssmDouble min_thickness,thickness,hab_fraction; 674 IssmDouble crevassedepth,surface_crevasse,surface,critical_fraction; 674 675 IssmDouble rho_ice,rho_water; 675 676 IssmDouble bed,water_depth; 676 677 IssmDouble levelset; … … 757 758 /*Get the DistanceToCalvingfront*/ 758 759 femmodel->elements->InputDuplicate(MaskIceLevelsetEnum,DistanceToCalvingfrontEnum); 759 760 femmodel->DistanceToFieldValue(MaskIceLevelsetEnum,0,DistanceToCalvingfrontEnum); 760 761 761 762 /*Vector of size number of nodes*/ 762 763 vec_constraint_nodes=new Vector<IssmDouble>(femmodel->nodes->NumberOfNodes()); 763 764 764 765 for(int i=0;i<femmodel->elements->Size();i++){ 765 Element* element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i)); 766 int numnodes = element->GetNumberOfNodes(); 767 Gauss* gauss = element->NewGauss(); 768 Input* crevassedepth_input = element->GetInput(CrevasseDepthEnum); _assert_(crevassedepth_input); 769 Input* bed_input = element->GetInput(BedEnum); _assert_(bed_input); 766 Element* element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i)); 767 int numnodes = element->GetNumberOfNodes(); 768 Gauss* gauss = element->NewGauss(); 769 Input* crevassedepth_input = element->GetInput(CrevasseDepthEnum); _assert_(crevassedepth_input); 770 Input* bed_input = element->GetInput(BedEnum); _assert_(bed_input); 771 Input* surface_crevasse_input = element->GetInput(SurfaceCrevasseEnum); _assert_(surface_crevasse_input); 772 Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); 773 Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input); 770 774 771 775 /*First, look at ice front and figure out if any of the nodes will be calved*/ 772 776 if(element->IsIcefront()){ … … 775 779 Node* node=element->GetNode(in); 776 780 crevassedepth_input->GetInputValue(&crevassedepth,gauss); 777 781 bed_input->GetInputValue(&bed,gauss); 778 if(crevassedepth>0. && bed<0.){ 782 surface_crevasse_input->GetInputValue(&surface_crevasse,gauss); 783 thickness_input->GetInputValue(&thickness,gauss); 784 surface_input->GetInputValue(&surface,gauss); 785 786 if((surface_crevasse-surface>0. || crevassedepth-thickness>0.) && bed<0.){ 779 787 vec_constraint_nodes->SetValue(node->Sid(),1.0,INS_VAL); 780 788 } 781 789 } … … 790 798 while(nflipped){ 791 799 local_nflipped=0; 792 800 for(int i=0;i<femmodel->elements->Size();i++){ 793 Element* element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i)); 794 int numnodes = element->GetNumberOfNodes(); 795 Gauss* gauss = element->NewGauss(); 796 Input* levelset_input = element->GetInput(DistanceToCalvingfrontEnum); _assert_(levelset_input); 797 Input* crevassedepth_input = element->GetInput(CrevasseDepthEnum); _assert_(crevassedepth_input); 798 Input* bed_input = element->GetInput(BedEnum); _assert_(bed_input); 801 Element* element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i)); 802 int numnodes = element->GetNumberOfNodes(); 803 Gauss* gauss = element->NewGauss(); 804 Input* levelset_input = element->GetInput(DistanceToCalvingfrontEnum); _assert_(levelset_input); 805 Input* crevassedepth_input = element->GetInput(CrevasseDepthEnum); _assert_(crevassedepth_input); 806 Input* bed_input = element->GetInput(BedEnum); _assert_(bed_input); 807 Input* surface_crevasse_input = element->GetInput(SurfaceCrevasseEnum); _assert_(surface_crevasse_input); 808 Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); 809 Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input); 799 810 800 811 /*Is this element connected to a node that should be calved*/ 801 812 bool isconnected = false; … … 815 826 levelset_input->GetInputValue(&levelset,gauss); 816 827 crevassedepth_input->GetInputValue(&crevassedepth,gauss); 817 828 bed_input->GetInputValue(&bed,gauss); 829 surface_crevasse_input->GetInputValue(&surface_crevasse,gauss); 830 thickness_input->GetInputValue(&thickness,gauss); 831 surface_input->GetInputValue(&surface,gauss); 818 832 819 if( crevassedepth>0.&& bed<0. && levelset>-300. && levelset<0. && constraint_nodes[node->Sid()]==0.){833 if((surface_crevasse-surface>0. || crevassedepth-thickness>0.) && bed<0. && levelset>-300. && levelset<0. && constraint_nodes[node->Sid()]==0.){ 820 834 local_nflipped++; 821 835 vec_constraint_nodes->SetValue(node->Sid(),1.0,INS_VAL); 822 836 }
Note:
See TracBrowser
for help on using the repository browser.
