Changeset 26984


Ignore:
Timestamp:
05/04/22 05:52:26 (3 years ago)
Author:
Mathieu Morlighem
Message:

CHG: moving more calving laws away from MovingFront

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

Legend:

Unmodified
Added
Removed

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

    r26983 r26984  
    273273        this->ComputeSigmaVM();
    274274
    275         IssmDouble  calvingratex[NUMVERTICES];
    276         IssmDouble  calvingratey[NUMVERTICES];
    277275        IssmDouble  calvingrate[NUMVERTICES];
    278276        IssmDouble  vx,vy;
    … …  
    317315                /*Assign values*/
    318316                if(bed>sealevel){
    319                         calvingratex[iv]=0.;
    320                         calvingratey[iv]=0.;
     317                        calvingrate[iv] = 0.;
    321318                }
    322319                else{
    323                         calvingratex[iv]=vx*sigma_vm/sigma_max;
    324                         calvingratey[iv]=vy*sigma_vm/sigma_max;
    325                 }
    326                 calvingrate[iv] =sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]);
     320                        calvingrate[iv] = sqrt(vx*vx+vy*vy)*sigma_vm/sigma_max;
     321                }
    327322        }
    328323
    329324        /*Add input*/
    330         this->AddBasalInput(CalvingratexEnum,&calvingratex[0],P1DGEnum);
    331         this->AddBasalInput(CalvingrateyEnum,&calvingratey[0],P1DGEnum);
    332325        this->AddBasalInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum);
    333 
     326        this->CalvingFromRate();
     327
     328        /*Extrude*/
     329        this->InputExtrude(CalvingCalvingrateEnum,-1);
    334330        this->InputExtrude(CalvingratexEnum,-1);
    335331        this->InputExtrude(CalvingrateyEnum,-1);
    336         this->InputExtrude(CalvingCalvingrateEnum,-1);
    337332}
    338333/*}}}*/
    … …  
    28712866}
    28722867/*}}}*/
    2873 void       Penta::MovingFrontalVelocity(void){/*{{{*/
     2868void          Penta::MovingFrontalVelocity(void){/*{{{*/
     2869
    28742870        if(!this->IsOnBase()) return;
    28752871        int        domaintype, calvinglaw, i;
    … …  
    29022898                case DefaultCalvingEnum:
    29032899                case CalvingVonmisesEnum:
    2904                         calvingrate_input = this->GetInput(CalvingCalvingrateEnum);     _assert_(calvingrate_input);
    2905                         meltingrate_input = this->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
    2906                         break;
    29072900                case CalvingLevermannEnum:
    29082901                        calvingratex_input=this->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
    … …  
    29412934                        case DefaultCalvingEnum:
    29422935                        case CalvingVonmisesEnum:
    2943                                 calvingrate_input->GetInputValue(&calvingrate,&gauss);
    2944                                 meltingrate_input->GetInputValue(&meltingrate,&gauss);
    2945 
    2946                                 if(groundedice<0) meltingrate = 0.;
    2947 
    2948                                 if(norm_dlsf>1.e-10)
    2949                                  for(i=0;i<dim;i++){
    2950                                          c[i]=calvingrate*dlsf[i]/norm_dlsf; m[i]=meltingrate*dlsf[i]/norm_dlsf;
    2951                                  }
    2952                                 else
    2953                                  for(i=0;i<dim;i++){
    2954                                          c[i]=0.; m[i]=0.;
    2955                                  }
    2956                                 break;
    2957 
    29582936                        case CalvingLevermannEnum:
    29592937                                calvingratex_input->GetInputValue(&c[0],&gauss);
    29602938                                calvingratey_input->GetInputValue(&c[1],&gauss);
    29612939                                meltingrate_input->GetInputValue(&meltingrate,&gauss);
    2962                                 norm_calving=0.;
    2963                                 for(i=0;i<dim;i++) norm_calving+=pow(c[i],2);
    2964                                 norm_calving=sqrt(norm_calving)+1.e-14;
    2965                                 for(i=0;i<dim;i++) m[i]=meltingrate*c[i]/norm_calving;
     2940                                if(groundedice<0) meltingrate = 0.;
     2941                                m[0]=meltingrate*dlsf[0]/norm_dlsf;
     2942                                m[1]=meltingrate*dlsf[1]/norm_dlsf;
    29662943                                break;
    29672944
    … …  
    30363013        this->AddInput(MovingFrontalVxEnum,&movingfrontvx[0],P1DGEnum);
    30373014        this->AddInput(MovingFrontalVyEnum,&movingfrontvy[0],P1DGEnum);
    3038 
    30393015        this->InputExtrude(MovingFrontalVxEnum,-1);
    30403016        this->InputExtrude(MovingFrontalVyEnum,-1);

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

    r26983 r26984  
    317317
    318318        /*Now compute calving rate*/
    319         IssmDouble  calvingratex[NUMVERTICES];
    320         IssmDouble  calvingratey[NUMVERTICES];
    321319        IssmDouble  calvingrate[NUMVERTICES];
    322320        IssmDouble  sigma_vm,vx,vy;
    … …  
    325323
    326324        /*Retrieve all inputs and parameters we will need*/
    327         Input* vx_input = this->GetInput(VxEnum); _assert_(vx_input);
    328         Input* vy_input = this->GetInput(VyEnum); _assert_(vy_input);
    329         Input* gr_input = this->GetInput(MaskOceanLevelsetEnum); _assert_(gr_input);
    330         Input* bs_input = this->GetInput(BaseEnum);                    _assert_(bs_input);
    331         Input* smax_fl_input = this->GetInput(CalvingStressThresholdFloatingiceEnum); _assert_(smax_fl_input);
    332         Input* smax_gr_input = this->GetInput(CalvingStressThresholdGroundediceEnum); _assert_(smax_gr_input);
    333         Input* sl_input  = this->GetInput(SealevelEnum); _assert_(sl_input);
    334         Input* sigma_vm_input  = this->GetInput(SigmaVMEnum); _assert_(sigma_vm_input);
     325        Input* vx_input       = this->GetInput(VxEnum); _assert_(vx_input);
     326        Input* vy_input       = this->GetInput(VyEnum); _assert_(vy_input);
     327        Input* gr_input       = this->GetInput(MaskOceanLevelsetEnum); _assert_(gr_input);
     328        Input* bs_input       = this->GetInput(BaseEnum);                    _assert_(bs_input);
     329        Input* smax_fl_input  = this->GetInput(CalvingStressThresholdFloatingiceEnum); _assert_(smax_fl_input);
     330        Input* smax_gr_input  = this->GetInput(CalvingStressThresholdGroundediceEnum); _assert_(smax_gr_input);
     331        Input* sl_input       = this->GetInput(SealevelEnum); _assert_(sl_input);
     332        Input* sigma_vm_input = this->GetInput(SigmaVMEnum); _assert_(sigma_vm_input);
    335333
    336334        /* Start looping on the number of vertices: */
    … …  
    354352                else
    355353                 sigma_max = sigma_max_grounded;
     354
    356355                /*Assign values*/
    357356                if(bed>sealevel){
    358                         calvingratex[iv]=0.;
    359                         calvingratey[iv]=0.;
     357                        calvingrate[iv] = 0.;
    360358                }
    361359                else{
    362                         calvingratex[iv]=vx*sigma_vm/sigma_max;
    363                         calvingratey[iv]=vy*sigma_vm/sigma_max;
    364                 }
    365                 calvingrate[iv] =sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]);
     360                        calvingrate[iv] = sqrt(vx*vx+vy*vy)*sigma_vm/sigma_max;
     361                }
    366362        }
    367363
    368364        /*Add input*/
    369         this->AddInput(CalvingratexEnum,&calvingratex[0],P1DGEnum);
    370         this->AddInput(CalvingrateyEnum,&calvingratey[0],P1DGEnum);
    371365        this->AddInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum);
     366   this->CalvingFromRate();
    372367}
    373368/*}}}*/
    … …  
    42894284        IssmDouble norm_dlsf, norm_calving, calvingrate, meltingrate, groundedice;
    42904285        IssmDouble migrationmax, calvinghaf, heaviside, haf_eps;
    4291         IssmDouble  xyz_list[NUMVERTICES][3];
    4292         IssmDouble  movingfrontvx[NUMVERTICES];
    4293         IssmDouble  movingfrontvy[NUMVERTICES];
    4294         IssmDouble  vel;
     4286        IssmDouble xyz_list[NUMVERTICES][3];
     4287        IssmDouble movingfrontvx[NUMVERTICES];
     4288        IssmDouble movingfrontvy[NUMVERTICES];
     4289        IssmDouble vel;
    42954290
    42964291        /* Get node coordinates and dof list: */
    … …  
    43264321        switch(calvinglaw){
    43274322                case DefaultCalvingEnum:
     4323                case CalvingVonmisesEnum:
    43284324                        calvingratex_input=this->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
    43294325                        calvingratey_input=this->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
    4330                         meltingrate_input = this->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
    4331                         break;
    4332                 case CalvingVonmisesEnum:
    4333                         calvingrate_input = this->GetInput(CalvingCalvingrateEnum);     _assert_(calvingrate_input);
    43344326                        meltingrate_input = this->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
    43354327                        break;
    … …  
    43894381                switch(calvinglaw){
    43904382                        case DefaultCalvingEnum:
     4383                        case CalvingVonmisesEnum:
    43914384                                calvingratex_input->GetInputValue(&c[0],&gauss);
    43924385                                calvingratey_input->GetInputValue(&c[1],&gauss);
    … …  
    44014394                                        }
    44024395                                }
    4403                                 break;
    4404                         case CalvingVonmisesEnum:
    4405                                 calvingrate_input->GetInputValue(&calvingrate,&gauss);
    4406                                 meltingrate_input->GetInputValue(&meltingrate,&gauss);
    4407                                 if(groundedice<0) meltingrate = 0.;
    4408 
    4409                                 if(norm_dlsf>1.e-10)
    4410                                  for(i=0;i<dim;i++){
    4411                                          c[i]=calvingrate*dlsf[i]/norm_dlsf; m[i]=meltingrate*dlsf[i]/norm_dlsf;
    4412                                  }
    4413                                 else
    4414                                  for(i=0;i<dim;i++){
    4415                                          c[i]=0.; m[i]=0.;
    4416                                  }
    44174396                                break;
    44184397
Note: See TracChangeset for help on using the changeset viewer.