Changeset 26876


Ignore:
Timestamp:
02/10/22 11:46:50 (3 years ago)
Author:
Mathieu Morlighem
Message:

CHG: simplifying calving

Location:
issm/trunk-jpl/src/c/classes/Elements
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r26875 r26876  
    28732873void       Penta::MovingFrontalVelocity(void){/*{{{*/
    28742874        if(!this->IsOnBase()) return;
    2875         int  dim, domaintype, calvinglaw, i;
     2875        int        domaintype, calvinglaw, i;
    28762876        IssmDouble v[3],w[3],c[3],m[3],dlsf[3];
    28772877        IssmDouble norm_dlsf, norm_calving, calvingrate, meltingrate, groundedice;
    28782878        IssmDouble migrationmax, calvinghaf, heaviside, haf_eps;
    2879         IssmDouble  xyz_list[NUMVERTICES][3];
    2880         IssmDouble  movingfrontvx[NUMVERTICES];
    2881         IssmDouble  movingfrontvy[NUMVERTICES];
    2882         IssmDouble  vel;
     2879        IssmDouble xyz_list[NUMVERTICES][3];
     2880        IssmDouble movingfrontvx[NUMVERTICES];
     2881        IssmDouble movingfrontvy[NUMVERTICES];
     2882        IssmDouble vel;
     2883        int        dim=2;
    28832884
    28842885        /* Get node coordinates and dof list: */
    28852886        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    28862887
    2887         Input* vx_input           = NULL;
    2888         Input* vy_input           = NULL;
    28892888        Input* calvingratex_input = NULL;
    28902889        Input* calvingratey_input = NULL;
    2891         Input* lsf_slopex_input   = NULL;
    2892         Input* lsf_slopey_input   = NULL;
    28932890        Input* calvingrate_input  = NULL;
    28942891        Input* meltingrate_input  = NULL;
    2895         Input* gr_input           = NULL;
     2892
     2893        /*Load levelset function gradients*/
     2894        Input *lsf_slopex_input = this->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
     2895        Input *lsf_slopey_input = this->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
     2896        Input *vx_input         = this->GetInput(VxEnum);                     _assert_(vx_input);
     2897        Input *vy_input         = this->GetInput(VyEnum);                     _assert_(vy_input);
     2898        Input *gr_input         = this->GetInput(MaskOceanLevelsetEnum);      _assert_(gr_input);
    28962899
    28972900        /*Get problem dimension and whether there is moving front or not*/
    2898         this->FindParam(&domaintype,DomainTypeEnum);
    28992901        this->FindParam(&calvinglaw,CalvingLawEnum);
    2900 
    2901         switch(domaintype){
    2902                 case Domain2DverticalEnum:   dim = 1; break;
    2903                 case Domain2DhorizontalEnum: dim = 2; break;
    2904                 case Domain3DEnum:           dim = 2; break;
    2905                 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
    2906         }
    2907         /*Load velocities*/
    2908         switch(domaintype){
    2909                 case Domain2DverticalEnum:
    2910                         vx_input=this->GetInput(VxEnum); _assert_(vx_input);
    2911                         break;
    2912                 case Domain2DhorizontalEnum:
    2913                         vx_input=this->GetInput(VxEnum); _assert_(vx_input);
    2914                         vy_input=this->GetInput(VyEnum); _assert_(vy_input);
    2915                         gr_input=this->GetInput(MaskOceanLevelsetEnum); _assert_(gr_input);
    2916                         break;
    2917                 case Domain3DEnum:
    2918                         vx_input=this->GetInput(VxAverageEnum); _assert_(vx_input);
    2919                         vy_input=this->GetInput(VyAverageEnum); _assert_(vy_input);
    2920                         gr_input=this->GetInput(MaskOceanLevelsetEnum); _assert_(gr_input);
    2921                         break;
    2922                 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
    2923         }
    2924 
    29252902        switch(calvinglaw){
    29262903                case DefaultCalvingEnum:
    29272904                case CalvingVonmisesEnum:
    2928                         lsf_slopex_input  = this->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
    2929                         if(dim==2) lsf_slopey_input  = this->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
    29302905                        calvingrate_input = this->GetInput(CalvingCalvingrateEnum);     _assert_(calvingrate_input);
    29312906                        meltingrate_input = this->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
    29322907                        break;
    29332908                case CalvingLevermannEnum:
    2934                         switch(domaintype){
    2935                                 case Domain2DverticalEnum:
    2936                                         calvingratex_input=this->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
    2937                                         break;
    2938                                 case Domain2DhorizontalEnum:
    2939                                         calvingratex_input=this->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
    2940                                         calvingratey_input=this->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
    2941                                         break;
    2942                                 case Domain3DEnum:
    2943                                         calvingratex_input=this->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
    2944                                         calvingratey_input=this->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
    2945                                         break;
    2946                                 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
    2947                         }
     2909                        calvingratex_input=this->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
     2910                        calvingratey_input=this->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
    29482911                        meltingrate_input = this->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
    29492912                        break;
    … …  
    29512914                case CalvingHabEnum:
    29522915                case CalvingCrevasseDepthEnum:
    2953                         lsf_slopex_input  = this->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
    2954                         if(dim==2) lsf_slopey_input  = this->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
    29552916                        meltingrate_input = this->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
    29562917                        break;
    29572918                case CalvingDev2Enum:
    29582919                        this->FindParam(&calvinghaf,CalvingHeightAboveFloatationEnum);
    2959                         lsf_slopex_input  = this->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
    2960                         if(dim==2) lsf_slopey_input  = this->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
    29612920                        calvingrate_input = this->GetInput(CalvingCalvingrateEnum);     _assert_(calvingrate_input);
    29622921                        meltingrate_input = this->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
    … …  
    29752934                vy_input->GetInputValue(&v[1],&gauss);
    29762935                gr_input->GetInputValue(&groundedice,&gauss);
     2936      lsf_slopex_input->GetInputValue(&dlsf[0],&gauss);
     2937      lsf_slopey_input->GetInputValue(&dlsf[1],&gauss);
     2938      norm_dlsf=sqrt(dlsf[0]*dlsf[0] + dlsf[1]*dlsf[1]);
    29772939
    29782940                /*Get calving speed*/
    … …  
    29802942                        case DefaultCalvingEnum:
    29812943                        case CalvingVonmisesEnum:
    2982                                 lsf_slopex_input->GetInputValue(&dlsf[0],&gauss);
    2983                                 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss);
    29842944                                calvingrate_input->GetInputValue(&calvingrate,&gauss);
    29852945                                meltingrate_input->GetInputValue(&meltingrate,&gauss);
    29862946
    29872947                                if(groundedice<0) meltingrate = 0.;
    2988 
    2989                                 norm_dlsf=0.;
    2990                                 for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
    2991                                 norm_dlsf=sqrt(norm_dlsf);
    29922948
    29932949                                if(norm_dlsf>1.e-10)
    … …  
    30032959                        case CalvingLevermannEnum:
    30042960                                calvingratex_input->GetInputValue(&c[0],&gauss);
    3005                                 if(dim==2) calvingratey_input->GetInputValue(&c[1],&gauss);
     2961                                calvingratey_input->GetInputValue(&c[1],&gauss);
    30062962                                meltingrate_input->GetInputValue(&meltingrate,&gauss);
    30072963                                norm_calving=0.;
    … …  
    30142970                        case CalvingHabEnum:
    30152971                        case CalvingCrevasseDepthEnum:
    3016                                 lsf_slopex_input->GetInputValue(&dlsf[0],&gauss);
    3017                                 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss);
    30182972                                meltingrate_input->GetInputValue(&meltingrate,&gauss);
    30192973                                if(groundedice<0) meltingrate = 0.;
    3020 
    3021                                 norm_dlsf=0.;
    3022                                 for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
    3023                                 norm_dlsf=sqrt(norm_dlsf);
    30242974
    30252975                                if(norm_dlsf>1.e-10)
    … …  
    30372987                        case CalvingDev2Enum:
    30382988                                  {
    3039                                         lsf_slopex_input->GetInputValue(&dlsf[0],&gauss);
    3040                                         if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss);
    30412989                                        calvingrate_input->GetInputValue(&calvingrate,&gauss);
    30422990                                        meltingrate_input->GetInputValue(&meltingrate,&gauss);
    … …  
    30623010                                                meltingrate=heaviside*meltingrate+0.;
    30633011                                        }
    3064 
    3065                                         norm_dlsf=0.;
    3066                                         for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
    3067                                         norm_dlsf=sqrt(norm_dlsf);
    30683012
    30693013                                        if(norm_dlsf>1.e-10)

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

    r26875 r26876  
    43204320        Input* calvingratex_input = NULL;
    43214321        Input* calvingratey_input = NULL;
    4322         Input* lsf_slopex_input   = NULL;
    43234322        Input* lsf_slopey_input   = NULL;
    43244323        Input* calvingrate_input  = NULL;
    43254324        Input* meltingrate_input  = NULL;
    43264325        Input* gr_input           = NULL;
     4326
    43274327
    43284328        /*Get problem dimension and whether there is moving front or not*/
    … …  
    43364336                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
    43374337        }
     4338
     4339        /*Load levelset function gradients*/
     4340        Input* lsf_slopex_input  = this->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
     4341        if(dim==2){lsf_slopey_input  = this->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);}
     4342
    43384343        /*Load velocities*/
    43394344        switch(domaintype){
    … …  
    43574362                case DefaultCalvingEnum:
    43584363                case CalvingVonmisesEnum:
    4359                         lsf_slopex_input  = this->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
    4360                         if(dim==2) lsf_slopey_input  = this->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
    43614364                        calvingrate_input = this->GetInput(CalvingCalvingrateEnum);     _assert_(calvingrate_input);
    43624365                        meltingrate_input = this->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
    … …  
    43824385                case CalvingHabEnum:
    43834386                case CalvingCrevasseDepthEnum:
    4384                         lsf_slopex_input  = this->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
    4385                         if(dim==2) lsf_slopey_input  = this->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
    4386                         meltingrate_input = this->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
    43874387                        break;
    43884388                case CalvingDev2Enum:
    43894389                        this->FindParam(&calvinghaf,CalvingHeightAboveFloatationEnum);
    4390                         lsf_slopex_input  = this->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
    4391                         if(dim==2) lsf_slopey_input  = this->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
    43924390                        calvingrate_input = this->GetInput(CalvingCalvingrateEnum);     _assert_(calvingrate_input);
    43934391                        meltingrate_input = this->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
    … …  
    43984396                        break;
    43994397                case CalvingParameterizationEnum:
    4400                         lsf_slopex_input  = this->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
    4401                         if(dim==2) lsf_slopey_input  = this->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
    44024398                        calvingrate_input = this->GetInput(CalvingCalvingrateEnum);     _assert_(calvingrate_input);
    44034399                        meltingrate_input = this->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
    … …  
    44164412                vy_input->GetInputValue(&v[1],&gauss);
    44174413                gr_input->GetInputValue(&groundedice,&gauss);
     4414                lsf_slopex_input->GetInputValue(&dlsf[0],&gauss);
     4415                if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss);
     4416
     4417                norm_dlsf=0.;
     4418                for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
     4419                norm_dlsf=sqrt(norm_dlsf);
    44184420
    44194421                /*Get calving speed*/
    … …  
    44214423                        case DefaultCalvingEnum:
    44224424                        case CalvingVonmisesEnum:
    4423                                 lsf_slopex_input->GetInputValue(&dlsf[0],&gauss);
    4424                                 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss);
    44254425                                calvingrate_input->GetInputValue(&calvingrate,&gauss);
    44264426                                meltingrate_input->GetInputValue(&meltingrate,&gauss);
    44274427                                if(groundedice<0) meltingrate = 0.;
    4428 
    4429                                 norm_dlsf=0.;
    4430                                 for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
    4431                                 norm_dlsf=sqrt(norm_dlsf);
    44324428
    44334429                                if(norm_dlsf>1.e-10)
    … …  
    44544450                        case CalvingHabEnum:
    44554451                        case CalvingCrevasseDepthEnum:
    4456                                 lsf_slopex_input->GetInputValue(&dlsf[0],&gauss);
    4457                                 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss);
    44584452                                meltingrate_input->GetInputValue(&meltingrate,&gauss);
    4459 
    4460                                 norm_dlsf=0.;
    4461                                 for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
    4462                                 norm_dlsf=sqrt(norm_dlsf);
    44634453
    44644454                                if(norm_dlsf>1.e-10)
    … …  
    44764466                        case CalvingDev2Enum:
    44774467                                  {
    4478                                         lsf_slopex_input->GetInputValue(&dlsf[0],&gauss);
    4479                                         if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss);
    44804468                                        calvingrate_input->GetInputValue(&calvingrate,&gauss);
    44814469                                        meltingrate_input->GetInputValue(&meltingrate,&gauss);
    … …  
    45024490                                        }
    45034491
    4504                                         norm_dlsf=0.;
    4505                                         for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
    4506                                         norm_dlsf=sqrt(norm_dlsf);
    4507 
    45084492                                        if(norm_dlsf>1.e-10)
    45094493                                         for(i=0;i<dim;i++){
    … …  
    45244508                                break;
    45254509                        case CalvingParameterizationEnum:
    4526             lsf_slopex_input->GetInputValue(&dlsf[0],&gauss);
    4527             if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss);
    45284510            calvingrate_input->GetInputValue(&calvingrate,&gauss);
    45294511            meltingrate_input->GetInputValue(&meltingrate,&gauss);
    45304512                                meltingrate = 0.;
    4531 
    4532             norm_dlsf=0.;
    4533             for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
    4534             norm_dlsf=sqrt(norm_dlsf);
    45354513
    45364514            if(norm_dlsf>1.e-10)
Note: See TracChangeset for help on using the changeset viewer.