source: issm/oecreview/Archive/23390-24306/ISSM-23606-23607.diff@ 24307

Last change on this file since 24307 was 24307, checked in by Mathieu Morlighem, 5 years ago

NEW: adding Archive/23390-24306
File size: 10.2 KB

  • ../trunk-jpl/src/c/classes/Elements/Tria.cpp

     
    331331        IssmDouble  xyz_list[NUMVERTICES][3];
    332332        IssmDouble  calvingrate[NUMVERTICES];
    333333        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;
    336335        IssmDouble  surface_crevasse[NUMVERTICES], basal_crevasse[NUMVERTICES], crevasse_depth[NUMVERTICES], H_surf, H_surfbasal;
    337         IssmDouble  strainparallel, straineffective;
     336        IssmDouble  B, strainparallel, straineffective;
    338337        IssmDouble  s_xx,s_xy,s_yy,s1,s2,stmp;
    339 
     338        int crevasse_opening_stress;
     339       
    340340        /* Get node coordinates and dof list: */
    341341        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    342342
    343         /*Get the critical fraction of thickness surface and basal crevasses penetrate for calving onset*/
    344         this->parameters->FindParam(&critical_fraction,CalvingCrevasseDepthEnum);
     343        /*retrieve the type of crevasse_opening_stress*/
     344        this->parameters->FindParam(&crevasse_opening_stress,CalvingCrevasseDepthEnum);
    345345
    346346        IssmDouble rho_ice        = this->GetMaterialParameter(MaterialsRhoIceEnum);
    347347        IssmDouble rho_seawater   = this->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    … …  
    360360        Input*   s_xx_input              = inputs->GetInput(DeviatoricStressxxEnum);     _assert_(s_xx_input);
    361361        Input*   s_xy_input              = inputs->GetInput(DeviatoricStressxyEnum);     _assert_(s_xy_input);
    362362        Input*   s_yy_input              = inputs->GetInput(DeviatoricStressyyEnum);     _assert_(s_yy_input);
     363        Input*  B_input  = inputs->GetInput(MaterialsRheologyBbarEnum);   _assert_(B_input);
    363364
    364365        /*Loop over all elements of this partition*/
    365366        GaussTria* gauss=new GaussTria();
    … …  
    368369
    369370                H_input->GetInputValue(&thickness,gauss);
    370371                bed_input->GetInputValue(&bed,gauss);
    371                 surface_input->GetInputValue(&float_depth,gauss);
     372                surface_input->GetInputValue(&surface,gauss);
    372373                strainrateparallel_input->GetInputValue(&strainparallel,gauss);
    373374                strainrateeffective_input->GetInputValue(&straineffective,gauss);
    374375                vx_input->GetInputValue(&vx,gauss);
    … …  
    377378                s_xx_input->GetInputValue(&s_xx,gauss);
    378379                s_xy_input->GetInputValue(&s_xy,gauss);
    379380                s_yy_input->GetInputValue(&s_yy,gauss);
     381                B_input->GetInputValue(&B,gauss);
    380382
    381383                vel=sqrt(vx*vx+vy*vy)+1.e-14;
    382384
    … …  
    387389                Ho = thickness - (rho_seawater/rho_ice) * (-bed);
    388390                if(Ho<0.)  Ho=0.;
    389391
    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 */
    394402                if (surface_crevasse[iv]<0.) {
    395403                        surface_crevasse[iv]=0.;
    396404                        water_height = 0.;
    397405                }
     406                if (basal_crevasse[iv]<0.) basal_crevasse[iv]=0.;
     407                if (bed>0.) basal_crevasse[iv] = 0.;
     408               
    398409                //if (surface_crevasse[iv]<water_height){
    399410                //      water_height = surface_crevasse[iv];
    400411                //}
    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       
    412417        }
    413418
    414419        this->inputs->AddInput(new TriaInput(SurfaceCrevasseEnum,&surface_crevasse[0],P1Enum));

  • ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp

     
    106106                case CalvingHabEnum:
    107107                        break;
    108108                case CalvingCrevasseDepthEnum:
    109                         parameters->AddObject(iomodel->CopyConstantObject("md.calving.critical_fraction",CalvingCrevasseDepthEnum));
     109                        parameters->AddObject(iomodel->CopyConstantObject("md.calving.crevasse_opening_stress",CalvingCrevasseDepthEnum));
    110110                        break;
    111111                case CalvingDev2Enum:
    112112                        parameters->AddObject(iomodel->CopyConstantObject("md.calving.height_above_floatation",CalvingHeightAboveFloatationEnum));
    … …  
    670670
    671671        /*Intermediaries*/
    672672        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;
    674675        IssmDouble  rho_ice,rho_water;
    675676        IssmDouble  bed,water_depth;
    676677        IssmDouble  levelset;
    … …  
    757758                /*Get the DistanceToCalvingfront*/
    758759                femmodel->elements->InputDuplicate(MaskIceLevelsetEnum,DistanceToCalvingfrontEnum);
    759760                femmodel->DistanceToFieldValue(MaskIceLevelsetEnum,0,DistanceToCalvingfrontEnum);
    760 
     761               
    761762                /*Vector of size number of nodes*/
    762763                vec_constraint_nodes=new Vector<IssmDouble>(femmodel->nodes->NumberOfNodes());
    763764
    764765                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);
    770774
    771775                        /*First, look at ice front and figure out if any of the nodes will be calved*/
    772776                        if(element->IsIcefront()){
    … …  
    775779                                        Node* node=element->GetNode(in);
    776780                                        crevassedepth_input->GetInputValue(&crevassedepth,gauss);
    777781                                        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.){
    779787                                                vec_constraint_nodes->SetValue(node->Sid(),1.0,INS_VAL);
    780788                                        }
    781789                                }
    … …  
    790798                while(nflipped){
    791799                        local_nflipped=0;
    792800                        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);
    799810
    800811                                /*Is this element connected to a node that should be calved*/
    801812                                bool isconnected = false;
    … …  
    815826                                                levelset_input->GetInputValue(&levelset,gauss);
    816827                                                crevassedepth_input->GetInputValue(&crevassedepth,gauss);
    817828                                                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);
    818832
    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.){
    820834                                                        local_nflipped++;
    821835                                                        vec_constraint_nodes->SetValue(node->Sid(),1.0,INS_VAL);
    822836                                                }
Note: See TracBrowser for help on using the repository browser.