Changeset 22488


Ignore:
Timestamp:
03/01/18 14:48:25 (7 years ago)
Author:
youngmc3
Message:

CHG: updated crevasse-depth calving law

Location:
issm/trunk-jpl/src/c
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r22474 r22488  
    264264                        lsf_slopex_input  = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
    265265                        if(dim==2) lsf_slopey_input  = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
    266                         calvingrate_input = basalelement->GetInput(CalvingCalvingrateEnum);     _assert_(calvingrate_input);
    267266                        meltingrate_input = basalelement->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
    268267                        break;
    … …  
    386385                                lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
    387386                                if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
    388                                 calvingrate_input->GetInputValue(&calvingrate,gauss);
    389387                                meltingrate_input->GetInputValue(&meltingrate,gauss);
    390388
    391                                 /*Limit calving rate to c <= v + 3 km/yr */
    392                                 vel=sqrt(v[0]*v[0] + v[1]*v[1]);
    393                                 if(calvingrate>calvingmax+vel) calvingrate = vel+calvingmax;
    394389                                if(groundedice<0) meltingrate = 0.;
    395390                               
    … …  
    400395                                if(norm_dlsf>1.e-10)
    401396                                 for(i=0;i<dim;i++){
    402                                          c[i]=calvingrate*dlsf[i]/norm_dlsf;
     397                                         c[i]=0.;
    403398                                         m[i]=meltingrate*dlsf[i]/norm_dlsf;
    404399                                 }
    … …  
    740735        /*Intermediaries*/
    741736        int         calvinglaw;
    742         IssmDouble  min_thickness,thickness,hab_fraction;
    743         IssmDouble  rho_ice,rho_water,constant_g,rheology_B,rheology_n;
     737        IssmDouble  min_thickness,thickness,hab_fraction,crevassedepth;
     738        IssmDouble  rho_ice,rho_water;
    744739        IssmDouble  bed,water_depth;
    745740        IssmDouble  levelset;
    746 
     741       
    747742        femmodel->parameters->FindParam(&calvinglaw,CalvingLawEnum);
    748743
    … …  
    793788                        Gauss*   gauss              = element->NewGauss();
    794789                        Input*   H_input            = element->GetInput(ThicknessEnum); _assert_(H_input);
    795                         Input*   bed                = element->GetInput(BedEnum); _assert_(bed);
     790                        Input*   bed_input          = element->GetInput(BedEnum); _assert_(bed_input);
    796791                        Input*   hab_fraction_input = element->GetInput(CalvingHabFractionEnum); _assert_(hab_fraction_input);
    797792                        Input*   ls_input           = element->GetInput(DistanceToCalvingfrontEnum); _assert_(ls_input);
    … …  
    802797                                Node* node=element->GetNode(in);
    803798                                H_input->GetInputValue(&thickness,gauss);
    804                                 bed->GetInputValue(&water_depth,gauss);
     799                                bed_input->GetInputValue(&water_depth,gauss);
    805800                                ls_input->GetInputValue(&levelset,gauss);
    806801                                hab_fraction_input->GetInputValue(&hab_fraction,gauss);
    807802
    808                                 if(thickness<((rho_water/rho_ice)*(1+hab_fraction)*-water_depth) && levelset>-300){
     803                                if(thickness<((rho_water/rho_ice)*(1+hab_fraction)*-water_depth) && levelset>-300. && levelset<0.){
    809804                                        node->ApplyConstraint(0,+1.);
    810805                                }
    … …  
    818813        }
    819814       
     815        if(calvinglaw==CalvingCrevasseDepthEnum){
     816       
     817                int                 nflipped,local_nflipped;
     818                Vector<IssmDouble>* vec_constraint_nodes = NULL;
     819                IssmDouble* constraint_nodes = NULL;
     820               
     821                /*Get the DistanceToCalvingfront*/
     822                femmodel->elements->InputDuplicate(MaskIceLevelsetEnum,DistanceToCalvingfrontEnum);
     823                femmodel->DistanceToFieldValue(MaskIceLevelsetEnum,0,DistanceToCalvingfrontEnum);
     824
     825                /*Vector of size number of nodes*/
     826                vec_constraint_nodes=new Vector<IssmDouble>(femmodel->nodes->NumberOfNodes(LevelsetAnalysisEnum));
     827
     828                for(int i=0;i<femmodel->elements->Size();i++){
     829                        Element* element  = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
     830                        int      numnodes = element->GetNumberOfNodes();
     831                        Gauss*   gauss    = element->NewGauss();
     832                        Input*   crevassedepth_input = element->GetInput(CrevasseDepthEnum); _assert_(crevassedepth_input);
     833                        Input*   bed_input = element->GetInput(BedEnum); _assert_(bed_input);
     834
     835                        /*First, look at ice front and figure out if any of the nodes will be calved*/
     836                        if(element->IsIcefront()){
     837                                for(int in=0;in<numnodes;in++){
     838                                        gauss->GaussNode(element->GetElementType(),in);
     839                                        Node* node=element->GetNode(in);
     840                                        crevassedepth_input->GetInputValue(&crevassedepth,gauss);
     841                                        bed_input->GetInputValue(&bed,gauss);
     842                                        if(crevassedepth>0. && bed<0.){
     843                                                vec_constraint_nodes->SetValue(node->Sid(),1.0,INS_VAL);
     844                                        }
     845                                }
     846                        }
     847                        delete gauss;
     848                }
     849                /*Assemble vector and serialize: */
     850                vec_constraint_nodes->Assemble();
     851                constraint_nodes=vec_constraint_nodes->ToMPISerial();
     852
     853                nflipped=1;
     854                while(nflipped){
     855                        local_nflipped=0;
     856                        for(int i=0;i<femmodel->elements->Size();i++){
     857                                Element* element  = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
     858                                int      numnodes = element->GetNumberOfNodes();
     859                                Gauss*   gauss    = element->NewGauss();
     860                                Input*   levelset_input  = element->GetInput(DistanceToCalvingfrontEnum); _assert_(levelset_input);
     861                                Input*   crevassedepth_input = element->GetInput(CrevasseDepthEnum); _assert_(crevassedepth_input);
     862                                Input*   bed_input = element->GetInput(BedEnum); _assert_(bed_input);
     863
     864                                /*Is this element connected to a node that should be calved*/
     865                                bool isconnected = false;
     866                                for(int in=0;in<numnodes;in++){
     867                                        Node* node=element->GetNode(in);
     868                                        if(constraint_nodes[node->Sid()]==1.){
     869                                                isconnected = true;
     870                                                break;
     871                                        }
     872                                }
     873
     874                                /*Check status if connected*/
     875                                if(isconnected){
     876                                        for(int in=0;in<numnodes;in++){
     877                                                gauss->GaussNode(element->GetElementType(),in);
     878                                                Node* node=element->GetNode(in);
     879                                                levelset_input->GetInputValue(&levelset,gauss);
     880                                                crevassedepth_input->GetInputValue(&crevassedepth,gauss);
     881                                                bed_input->GetInputValue(&bed,gauss);
     882
     883                                                if(crevassedepth>0. && bed<0. && levelset>-300. && levelset<0. && constraint_nodes[node->Sid()]==0.){
     884                                                        local_nflipped++;
     885                                                        vec_constraint_nodes->SetValue(node->Sid(),1.0,INS_VAL);
     886                                                }
     887                                        }
     888                                }
     889                        }
     890
     891                        /*Count how many new nodes were found*/
     892                        ISSM_MPI_Allreduce(&local_nflipped,&nflipped,1,ISSM_MPI_INT,ISSM_MPI_SUM,IssmComm::GetComm());
     893                        //_printf0_("Found "<<nflipped<<" to flip\n");
     894
     895                        /*Assemble and serialize flag vector*/
     896                        vec_constraint_nodes->Assemble();
     897                        xDelete<IssmDouble>(constraint_nodes);
     898                        constraint_nodes=vec_constraint_nodes->ToMPISerial();
     899                }
     900                /*Free ressources:*/
     901                delete vec_constraint_nodes;
     902
     903                /*Contrain the nodes that will be calved*/
     904                for(int i=0;i<femmodel->elements->Size();i++){
     905                        Element* element  = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
     906                        int      numnodes = element->GetNumberOfNodes();
     907                        Gauss*   gauss    = element->NewGauss();
     908                        /*Potentially constrain nodes of this element*/
     909                        for(int in=0;in<numnodes;in++){
     910                                gauss->GaussNode(element->GetElementType(),in);
     911                                Node* node=element->GetNode(in);
     912                               
     913                                if(constraint_nodes[node->Sid()]){
     914                                        node->ApplyConstraint(0,+1.);
     915                                }
     916                                else {
     917                                        /* no ice, set no spc */
     918                                        node->DofInFSet(0);
     919                                }
     920                        }
     921                        delete gauss;
     922                }
     923                xDelete<IssmDouble>(constraint_nodes);
     924        }
     925
    820926        /*Default, do nothing*/
    821927        return;

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

    r22473 r22488  
    302302        IssmDouble  calvingrate[NUMVERTICES];
    303303        IssmDouble  vx,vy,vel;
    304         IssmDouble  rheology_B,critical_fraction,water_height;
    305         IssmDouble  bathymetry,Ho,thickness,float_depth,groundedice;
     304        IssmDouble  critical_fraction,water_height;
     305        IssmDouble  bed,Ho,thickness,float_depth;
    306306        IssmDouble  surface_crevasse[NUMVERTICES], basal_crevasse[NUMVERTICES], crevasse_depth[NUMVERTICES], H_surf, H_surfbasal;
    307307        IssmDouble  strainparallel, straineffective;
    308         IssmDouble  yts;
     308        IssmDouble  s_xx,s_xy,s_yy,s1,s2,stmp;
    309309
    310310        /* Get node coordinates and dof list: */
    … …  
    313313        /*Get the critical fraction of thickness surface and basal crevasses penetrate for calving onset*/
    314314        this->parameters->FindParam(&critical_fraction,CalvingCrevasseDepthEnum);
    315         this->parameters->FindParam(&yts,ConstantsYtsEnum);
    316315               
    317316        IssmDouble rho_ice        = this->GetMaterialParameter(MaterialsRhoIceEnum);
    … …  
    321320        IssmDouble rheology_n     = this->GetMaterialParameter(MaterialsRheologyNEnum);
    322321
    323         Input*   B_input       = inputs->GetInput(MaterialsRheologyBbarEnum);   _assert_(B_input);
    324         Input*   H_input       = inputs->GetInput(ThicknessEnum); _assert_(H_input);
    325         Input*   bed           = inputs->GetInput(BedEnum); _assert_(bed);
    326         Input*   surface       = inputs->GetInput(SurfaceEnum); _assert_(surface);
    327         Input*  strainrateparallel  = inputs->GetInput(StrainRateparallelEnum);  _assert_(strainrateparallel);
    328         Input*  strainrateeffective = inputs->GetInput(StrainRateeffectiveEnum); _assert_(strainrateeffective);
    329         Input*  vx_input = inputs->GetInput(VxEnum); _assert_(vx_input);
    330         Input*  vy_input = inputs->GetInput(VxEnum); _assert_(vy_input);
    331         Input*  gr_input = inputs->GetInput(MaskGroundediceLevelsetEnum); _assert_(gr_input);
    332         Input*   waterheight_input = inputs->GetInput(WaterheightEnum); _assert_(waterheight_input);
    333 
     322        Input*   H_input                 = inputs->GetInput(ThicknessEnum); _assert_(H_input);
     323        Input*   bed_input               = inputs->GetInput(BedEnum); _assert_(bed_input);
     324        Input*   surface_input           = inputs->GetInput(SurfaceEnum); _assert_(surface_input);
     325        Input*  strainrateparallel_input  = inputs->GetInput(StrainRateparallelEnum);  _assert_(strainrateparallel_input);
     326        Input*  strainrateeffective_input = inputs->GetInput(StrainRateeffectiveEnum); _assert_(strainrateeffective_input);
     327        Input*  vx_input                  = inputs->GetInput(VxEnum); _assert_(vx_input);
     328        Input*  vy_input                  = inputs->GetInput(VxEnum); _assert_(vy_input);
     329        Input*   waterheight_input       = inputs->GetInput(WaterheightEnum); _assert_(waterheight_input);
     330        Input*   s_xx_input              = inputs->GetInput(DeviatoricStressxxEnum);     _assert_(s_xx_input);
     331        Input*   s_xy_input              = inputs->GetInput(DeviatoricStressxyEnum);     _assert_(s_xy_input);
     332        Input*   s_yy_input              = inputs->GetInput(DeviatoricStressyyEnum);     _assert_(s_yy_input);
     333       
    334334        /*Loop over all elements of this partition*/
    335335        GaussTria* gauss=new GaussTria();
    … …  
    337337                gauss->GaussVertex(iv);
    338338       
    339                 B_input->GetInputValue(&rheology_B,gauss);
    340339                H_input->GetInputValue(&thickness,gauss);
    341                 bed->GetInputValue(&bathymetry,gauss);
    342                 surface->GetInputValue(&float_depth,gauss);
    343                 strainrateparallel->GetInputValue(&strainparallel,gauss);
    344                 strainrateeffective->GetInputValue(&straineffective,gauss);
     340                bed_input->GetInputValue(&bed,gauss);
     341                surface_input->GetInputValue(&float_depth,gauss);
     342                strainrateparallel_input->GetInputValue(&strainparallel,gauss);
     343                strainrateeffective_input->GetInputValue(&straineffective,gauss);
    345344                vx_input->GetInputValue(&vx,gauss);
    346345                vy_input->GetInputValue(&vy,gauss);
    347                 gr_input->GetInputValue(&groundedice,gauss);
    348346                waterheight_input->GetInputValue(&water_height,gauss);
     347                s_xx_input->GetInputValue(&s_xx,gauss);
     348                s_xy_input->GetInputValue(&s_xy,gauss);
     349                s_yy_input->GetInputValue(&s_yy,gauss);
     350               
    349351                vel=sqrt(vx*vx+vy*vy)+1.e-14;
    350352
    351                 Ho = thickness - (rho_seawater/rho_ice) * (-bathymetry);
     353                s1=(s_xx+s_yy)/2.+sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2));
     354                s2=(s_xx+s_yy)/2.-sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2));
     355                if(fabs(s2)>fabs(s1)){stmp=s2; s2=s1; s1=stmp;}
     356               
     357                Ho = thickness - (rho_seawater/rho_ice) * (-bed);
    352358                if(Ho<0.)  Ho=0.;
    353359
    354360                /*Otero2010: balance between the tensile deviatoric stress and ice overburden pressure*/
    355361                /*surface crevasse*/
    356                 surface_crevasse[iv] = rheology_B * strainparallel * pow(straineffective, ((1 / rheology_n)-1)) / (rho_ice * constant_g);
    357                 //surface_crevasse[iv] = 2 * rheology_B * pow(max(strainparallel,0.),(1/rheology_n)) / (rho_ice * constant_g);
     362                //surface_crevasse[iv] = rheology_B * strainparallel * pow(straineffective, ((1 / rheology_n)-1)) / (rho_ice * constant_g);
     363                surface_crevasse[iv] = s1 / (rho_ice*constant_g);
    358364                if (surface_crevasse[iv]<0.) {
    359365                        surface_crevasse[iv]=0.;
    360366                        water_height = 0.;
    361367                }
    362                 if (surface_crevasse[iv]<water_height){
    363                         water_height = surface_crevasse[iv];
    364                 }
     368                //if (surface_crevasse[iv]<water_height){
     369                //      water_height = surface_crevasse[iv];
     370                //}
    365371               
    366372                /*basal crevasse*/
    367                 basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice)) * (rheology_B * strainparallel * pow(straineffective,((1/rheology_n)-1)) / (rho_ice*constant_g) - Ho);
    368                 //basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice)) * (2 * rheology_B * pow(max(strainparallel,0.),(1/rheology_n)) / (rho_ice*constant_g) - Ho);
     373                //basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice)) * (rheology_B * strainparallel * pow(straineffective,((1/rheology_n)-1)) / (rho_ice*constant_g) - Ho);
     374                basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice))* (s1/ (rho_ice*constant_g)-Ho);
    369375                if (basal_crevasse[iv]<0.) basal_crevasse[iv]=0.;
    370                 if (bathymetry>0.) basal_crevasse[iv] = 0.;
     376                if (bed>0.) basal_crevasse[iv] = 0.;
    371377       
    372378                H_surf = surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height - critical_fraction*float_depth;
    … …  
    374380               
    375381                crevasse_depth[iv] = max(H_surf,H_surfbasal);
     382        }
    376383       
    377                 /*Assign values */
    378                 if(crevasse_depth[iv]>=0. && bathymetry<=0.){
    379                  calvingrate[iv] = vel+3000./yts;
    380                 }       
    381                 else
    382                  calvingrate[iv]=0.;
    383         }
    384 
    385384        this->inputs->AddInput(new TriaInput(SurfaceCrevasseEnum,&surface_crevasse[0],P1Enum));
    386385        this->inputs->AddInput(new TriaInput(BasalCrevasseEnum,&basal_crevasse[0],P1Enum));
    387386        this->inputs->AddInput(new TriaInput(CrevasseDepthEnum,&crevasse_depth[0],P1Enum));
    388         this->inputs->AddInput(new TriaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum));
    389387
    390388        delete gauss;

  • issm/trunk-jpl/src/c/modules/Calvingx/Calvingx.cpp

    r22199 r22488  
    2323                        femmodel->StrainRateparallelx();
    2424                        femmodel->StrainRateeffectivex();
     25                        femmodel->DeviatoricStressx();
    2526                        femmodel->ElementOperationx(&Element::CalvingCrevasseDepth);
    2627                        break;
Note: See TracChangeset for help on using the changeset viewer.