Changeset 25839


Ignore:
Timestamp:
12/08/20 13:10:07 (4 years ago)
Author:
Cheng Gong
Message:

CHG: move the calving laws out from LevelsetAnalysis. The advection velocity in LevelsetAnalysis is now defined by MovingFrontalVxEnum and MovingFrontalVyEnum.

Location:
issm/trunk-jpl
Files:
24 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/Makefile.am

    r25811 r25839  
    498498issm_sources += ./analyses/LevelsetAnalysis.cpp
    499499issm_sources += ./modules/Calvingx/Calvingx.cpp
     500issm_sources += ./modules/MovingFrontalVelx/MovingFrontalVelx.cpp
    500501issm_sources += ./modules/KillIcebergsx/KillIcebergsx.cpp
    501502endif

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

    r25741 r25839  
    186186
    187187        /*Intermediaries */
    188         int  stabilization,dim, domaintype, calvinglaw;
     188        int  stabilization,dim,domaintype;
    189189        int i,j,k,row, col;
    190190        IssmDouble kappa;
    191191        IssmDouble Jdet, dt, D_scalar;
    192192        IssmDouble h,hx,hy,hz;
    193         IssmDouble vel,v[3],w[3],c[3],m[3],dlsf[3];
    194         IssmDouble norm_dlsf, norm_calving, calvingrate, meltingrate, groundedice;
    195         IssmDouble migrationmax, calvinghaf, heaviside, haf_eps;
     193        IssmDouble vel,w[3];
     194        IssmDouble migrationmax;
    196195        IssmDouble* xyz_list = NULL;
    197196
    198197        /*Get problem dimension and whether there is moving front or not*/
    199198        basalelement->FindParam(&domaintype,DomainTypeEnum);
    200         basalelement->FindParam(&calvinglaw,CalvingLawEnum);
    201199        basalelement->FindParam(&stabilization,LevelsetStabilizationEnum);
    202200        switch(domaintype){
    … …  
    206204                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
    207205        }
    208 
    209         /*Calving threshold*/
    210 
    211206        /*Fetch number of nodes and dof for this finite element*/
    212207        int numnodes    = basalelement->GetNumberOfNodes();
    … …  
    225220        basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
    226221        basalelement->FindParam(&migrationmax,MigrationMaxEnum);
    227         Input* vx_input           = NULL;
    228         Input* vy_input           = NULL;
    229         Input* calvingratex_input = NULL;
    230         Input* calvingratey_input = NULL;
    231         Input* lsf_slopex_input   = NULL;
    232         Input* lsf_slopey_input   = NULL;
    233         Input* calvingrate_input  = NULL;
    234         Input* meltingrate_input  = NULL;
    235         Input* gr_input           = NULL;
     222
     223        Input* mf_vx_input        = NULL;
     224        Input* mf_vy_input        = NULL;
    236225
    237226        /*Load velocities*/
    238227        switch(domaintype){
    239228                case Domain2DverticalEnum:
    240                         vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
     229                        mf_vx_input=basalelement->GetInput(MovingFrontalVxEnum); _assert_(mf_vx_input);
    241230                        break;
    242231                case Domain2DhorizontalEnum:
    243                         vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
    244                         vy_input=basalelement->GetInput(VyEnum); _assert_(vy_input);
    245                         gr_input=basalelement->GetInput(MaskOceanLevelsetEnum); _assert_(gr_input);
     232                        mf_vx_input=basalelement->GetInput(MovingFrontalVxEnum); _assert_(mf_vx_input);
     233                        mf_vy_input=basalelement->GetInput(MovingFrontalVyEnum); _assert_(mf_vy_input);
    246234                        break;
    247235                case Domain3DEnum:
    248                         vx_input=basalelement->GetInput(VxAverageEnum); _assert_(vx_input);
    249                         vy_input=basalelement->GetInput(VyAverageEnum); _assert_(vy_input);
    250                         gr_input=basalelement->GetInput(MaskOceanLevelsetEnum); _assert_(gr_input);
     236                        mf_vx_input=basalelement->GetInput(MovingFrontalVxEnum); _assert_(mf_vx_input);
     237                        mf_vy_input=basalelement->GetInput(MovingFrontalVyEnum); _assert_(mf_vy_input);
    251238                        break;
    252239                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
    253         }
    254 
    255         /*Load calving inputs*/
    256         switch(calvinglaw){
    257                 case DefaultCalvingEnum:
    258                 case CalvingVonmisesEnum:
    259                         lsf_slopex_input  = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
    260                         if(dim==2) lsf_slopey_input  = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
    261                         calvingrate_input = basalelement->GetInput(CalvingCalvingrateEnum);     _assert_(calvingrate_input);
    262                         meltingrate_input = basalelement->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
    263                         break;
    264                 case CalvingLevermannEnum:
    265                         switch(domaintype){
    266                                 case Domain2DverticalEnum:
    267                                         calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
    268                                         break;
    269                                 case Domain2DhorizontalEnum:
    270                                         calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
    271                                         calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
    272                                         break;
    273                                 case Domain3DEnum:
    274                                         calvingratex_input=basalelement->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
    275                                         calvingratey_input=basalelement->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
    276                                         break;
    277                                 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
    278                         }
    279                         meltingrate_input = basalelement->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
    280                         break;
    281                 case CalvingMinthicknessEnum:
    282                         lsf_slopex_input  = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
    283                         if(dim==2) lsf_slopey_input  = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
    284                         meltingrate_input = basalelement->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
    285                         break;
    286                 case CalvingHabEnum:
    287                         lsf_slopex_input  = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
    288                         if(dim==2) lsf_slopey_input  = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
    289                         meltingrate_input = basalelement->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
    290                         break;
    291                 case CalvingCrevasseDepthEnum:
    292                         lsf_slopex_input  = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
    293                         if(dim==2) lsf_slopey_input  = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
    294                         meltingrate_input = basalelement->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
    295                         break;
    296                 case CalvingDev2Enum:
    297                         basalelement->FindParam(&calvinghaf,CalvingHeightAboveFloatationEnum);
    298                         lsf_slopex_input  = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
    299                         if(dim==2) lsf_slopey_input  = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
    300                         calvingrate_input = basalelement->GetInput(CalvingCalvingrateEnum);     _assert_(calvingrate_input);
    301                         meltingrate_input = basalelement->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
    302                         break;
    303                 default:
    304                         _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
    305240        }
    306241
    … …  
    324259                }
    325260
    326                 /* Advection */
    327                 vx_input->GetInputValue(&v[0],gauss);
    328                 vy_input->GetInputValue(&v[1],gauss);
    329                 gr_input->GetInputValue(&groundedice,gauss);
    330 
    331                 /*Get calving speed*/
    332                 switch(calvinglaw){
    333                         case DefaultCalvingEnum:
    334                         case CalvingVonmisesEnum:
    335                                 lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
    336                                 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
    337                                 calvingrate_input->GetInputValue(&calvingrate,gauss);
    338                                 meltingrate_input->GetInputValue(&meltingrate,gauss);
    339 
    340                                 /*Limit calving rate to c <= v + 3 km/yr */
    341                                 vel=sqrt(v[0]*v[0] + v[1]*v[1]);
    342                                 if(calvingrate>migrationmax+vel) calvingrate = vel+migrationmax;
    343                                 if(groundedice<0) meltingrate = 0.;
    344 
    345                                 norm_dlsf=0.;
    346                                 for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
    347                                 norm_dlsf=sqrt(norm_dlsf);
    348 
    349                                 if(norm_dlsf>1.e-10)
    350                                  for(i=0;i<dim;i++){
    351                                          c[i]=calvingrate*dlsf[i]/norm_dlsf; m[i]=meltingrate*dlsf[i]/norm_dlsf;
    352                                  }
    353                                 else
    354                                  for(i=0;i<dim;i++){
    355                                          c[i]=0.; m[i]=0.;
    356                                  }
    357                                 break;
    358 
    359                         case CalvingLevermannEnum:
    360                                 calvingratex_input->GetInputValue(&c[0],gauss);
    361                                 if(dim==2) calvingratey_input->GetInputValue(&c[1],gauss);
    362                                 meltingrate_input->GetInputValue(&meltingrate,gauss);
    363                                 norm_calving=0.;
    364                                 for(i=0;i<dim;i++) norm_calving+=pow(c[i],2);
    365                                 norm_calving=sqrt(norm_calving)+1.e-14;
    366                                 for(i=0;i<dim;i++) m[i]=meltingrate*c[i]/norm_calving;
    367                                 break;
    368 
    369                         case CalvingMinthicknessEnum:
    370                                 lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
    371                                 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
    372                                 meltingrate_input->GetInputValue(&meltingrate,gauss);
    373 
    374                                 norm_dlsf=0.;
    375                                 for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
    376                                 norm_dlsf=sqrt(norm_dlsf);
    377 
    378                                 if(norm_dlsf>1.e-10)
    379                                  for(i=0;i<dim;i++){
    380                                          c[i]=0.;
    381                                          m[i]=meltingrate*dlsf[i]/norm_dlsf;
    382                                  }
    383                                 else
    384                                  for(i=0;i<dim;i++){
    385                                          c[i]=0.;
    386                                          m[i]=0.;
    387                                  }
    388                                 break;
    389 
    390                         case CalvingHabEnum:
    391                                 lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
    392                                 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
    393                                 meltingrate_input->GetInputValue(&meltingrate,gauss);
    394 
    395                                 norm_dlsf=0.;
    396                                 for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
    397                                 norm_dlsf=sqrt(norm_dlsf);
    398 
    399                                 if(norm_dlsf>1.e-10)
    400                                  for(i=0;i<dim;i++){
    401                                          c[i]=0.;
    402                                          m[i]=meltingrate*dlsf[i]/norm_dlsf;
    403                                  }
    404                                 else
    405                                  for(i=0;i<dim;i++){
    406                                          c[i]=0.;
    407                                          m[i]=0.;
    408                                  }
    409                                 break;
    410 
    411                         case CalvingCrevasseDepthEnum:
    412                                 lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
    413                                 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
    414                                 meltingrate_input->GetInputValue(&meltingrate,gauss);
    415 
    416                                 if(groundedice<0) meltingrate = 0.;
    417 
    418                                 norm_dlsf=0.;
    419                                 for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
    420                                 norm_dlsf=sqrt(norm_dlsf);
    421 
    422                                 if(norm_dlsf>1.e-10)
    423                                  for(i=0;i<dim;i++){
    424                                          c[i]=0.;
    425                                          m[i]=meltingrate*dlsf[i]/norm_dlsf;
    426                                  }
    427                                 else
    428                                  for(i=0;i<dim;i++){
    429                                          c[i]=0.;
    430                                          m[i]=0.;
    431                                  }
    432                                 break;
    433 
    434                         case CalvingDev2Enum:
    435                                   {
    436                                         lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
    437                                         if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
    438                                         calvingrate_input->GetInputValue(&calvingrate,gauss);
    439                                         meltingrate_input->GetInputValue(&meltingrate,gauss);
    440                                         gr_input->GetInputValue(&groundedice,gauss);
    441 
    442                                         //idea: no retreat on ice above critical calving height "calvinghaf" . Limit using regularized Heaviside function.
    443                                         vel=sqrt(v[0]*v[0] + v[1]*v[1]);
    444                                         haf_eps=10.;
    445                                         if(groundedice-calvinghaf<=-haf_eps){
    446                                                 // ice floats freely below calvinghaf: calve freely
    447                                                 // undercutting has no effect:
    448                                                 meltingrate=0.;
    449                                         }
    450                                         else if(groundedice-calvinghaf>=haf_eps){
    451                                                 // ice is well above calvinghaf -> no calving back, i.e. limit calving rate to ice velocity
    452                                                 calvingrate=min(calvingrate,vel);
    453                                                 // ice is almost grounded: frontal undercutting has maximum effect (do nothing).
    454                                         }
    455                                         else{ // ice is close to calvinghaf: smooth transition between limitation and free calving.
    456                                                 //heaviside: 0 for floating, 1 for grounded
    457                                                 heaviside=(groundedice-calvinghaf+haf_eps)/(2.*haf_eps) + sin(PI*(groundedice-calvinghaf)/haf_eps)/(2.*PI);
    458                                                 calvingrate=heaviside*(min(calvingrate,vel)-calvingrate)+calvingrate;
    459                                                 meltingrate=heaviside*meltingrate+0.;
    460                                         }
    461 
    462                                         norm_dlsf=0.;
    463                                         for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
    464                                         norm_dlsf=sqrt(norm_dlsf);
    465 
    466                                         if(norm_dlsf>1.e-10)
    467                                          for(i=0;i<dim;i++){
    468                                                  c[i]=calvingrate*dlsf[i]/norm_dlsf;
    469                                                  m[i]=meltingrate*dlsf[i]/norm_dlsf;
    470                                          }
    471                                         else
    472                                          for(i=0;i<dim;i++){
    473                                                  c[i]=0.;
    474                                                  m[i]=0.;
    475                                          }
    476                                         break;
    477                                   }
    478 
    479                         default:
    480                                 _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
    481                 }
    482 
    483                 /*Levelset speed is ice velocity - calving rate*/
    484                 for(i=0;i<dim;i++) w[i]=v[i]-c[i]-m[i];
     261                /* Levelset speed */
     262                mf_vx_input->GetInputValue(&w[0], gauss);
     263                mf_vy_input->GetInputValue(&w[1], gauss);
     264
     265                /* Apply limiter to the migration rate */               
     266                vel = 0.;
     267                for(i=0;i<dim;i++) vel += w[i]*w[i];
     268                vel = sqrt(vel)+1e-14;
     269                /* !!NOTE: This is different from the previous version 25838 (and before). The current threshold restrict both advance and retreat velocity. */
     270                if (vel > migrationmax) {
     271                        for(i=0;i<dim;i++) w[i] = w[i]/vel*migrationmax;
     272                }
    485273
    486274                /*Compute D*/

  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r25752 r25839  
    297297                virtual IssmDouble Misfit(int modelenum,int observationenum,int weightsenum)=0;
    298298                virtual IssmDouble MisfitArea(int weightsenum)=0;
     299                virtual void       MovingFrontalVelocity(void){_error_("not implemented yet");};
    299300                virtual Gauss*     NewGauss(void)=0;
    300301                virtual Gauss*     NewGauss(int order)=0;

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

    r25723 r25839  
    27452745
    27462746        return minlength;
     2747}
     2748/*}}}*/
     2749void       Penta::MovingFrontalVelocity(void){/*{{{*/
     2750        if(!this->IsOnBase()) return;
     2751        int  dim, domaintype, calvinglaw, i;
     2752        IssmDouble v[3],w[3],c[3],m[3],dlsf[3];
     2753        IssmDouble norm_dlsf, norm_calving, calvingrate, meltingrate, groundedice;
     2754        IssmDouble migrationmax, calvinghaf, heaviside, haf_eps;
     2755        IssmDouble  xyz_list[NUMVERTICES][3];
     2756        IssmDouble  movingfrontvx[NUMVERTICES];
     2757        IssmDouble  movingfrontvy[NUMVERTICES];
     2758        IssmDouble  vel;
     2759       
     2760        /* Get node coordinates and dof list: */
     2761        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     2762
     2763        Input* vx_input           = NULL;
     2764        Input* vy_input           = NULL;
     2765        Input* calvingratex_input = NULL;
     2766        Input* calvingratey_input = NULL;
     2767        Input* lsf_slopex_input   = NULL;
     2768        Input* lsf_slopey_input   = NULL;
     2769        Input* calvingrate_input  = NULL;
     2770        Input* meltingrate_input  = NULL;
     2771        Input* gr_input           = NULL;
     2772
     2773        /*Get problem dimension and whether there is moving front or not*/
     2774        this->FindParam(&domaintype,DomainTypeEnum);
     2775        this->FindParam(&calvinglaw,CalvingLawEnum);
     2776
     2777        switch(domaintype){
     2778                case Domain2DverticalEnum:   dim = 1; break;
     2779                case Domain2DhorizontalEnum: dim = 2; break;
     2780                case Domain3DEnum:           dim = 2; break;
     2781                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     2782        }
     2783        /*Load velocities*/
     2784        switch(domaintype){
     2785                case Domain2DverticalEnum:
     2786                        vx_input=this->GetInput(VxEnum); _assert_(vx_input);
     2787                        break;
     2788                case Domain2DhorizontalEnum:
     2789                        vx_input=this->GetInput(VxEnum); _assert_(vx_input);
     2790                        vy_input=this->GetInput(VyEnum); _assert_(vy_input);
     2791                        gr_input=this->GetInput(MaskOceanLevelsetEnum); _assert_(gr_input);
     2792                        break;
     2793                case Domain3DEnum:
     2794                        vx_input=this->GetInput(VxAverageEnum); _assert_(vx_input);
     2795                        vy_input=this->GetInput(VyAverageEnum); _assert_(vy_input);
     2796                        gr_input=this->GetInput(MaskOceanLevelsetEnum); _assert_(gr_input);
     2797                        break;
     2798                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     2799        }
     2800
     2801        switch(calvinglaw){
     2802                case DefaultCalvingEnum:
     2803                case CalvingVonmisesEnum:
     2804                        lsf_slopex_input  = this->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
     2805                        if(dim==2) lsf_slopey_input  = this->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
     2806                        calvingrate_input = this->GetInput(CalvingCalvingrateEnum);     _assert_(calvingrate_input);
     2807                        meltingrate_input = this->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
     2808                        break;
     2809                case CalvingLevermannEnum:
     2810                        switch(domaintype){
     2811                                case Domain2DverticalEnum:
     2812                                        calvingratex_input=this->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
     2813                                        break;
     2814                                case Domain2DhorizontalEnum:
     2815                                        calvingratex_input=this->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
     2816                                        calvingratey_input=this->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
     2817                                        break;
     2818                                case Domain3DEnum:
     2819                                        calvingratex_input=this->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
     2820                                        calvingratey_input=this->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
     2821                                        break;
     2822                                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     2823                        }
     2824                        meltingrate_input = this->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
     2825                        break;
     2826                case CalvingMinthicknessEnum:
     2827                        lsf_slopex_input  = this->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
     2828                        if(dim==2) lsf_slopey_input  = this->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
     2829                        meltingrate_input = this->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
     2830                        break;
     2831                case CalvingHabEnum:
     2832                        lsf_slopex_input  = this->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
     2833                        if(dim==2) lsf_slopey_input  = this->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
     2834                        meltingrate_input = this->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
     2835                        break;
     2836                case CalvingCrevasseDepthEnum:
     2837                        lsf_slopex_input  = this->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
     2838                        if(dim==2) lsf_slopey_input  = this->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
     2839                        meltingrate_input = this->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
     2840                        break;
     2841                case CalvingDev2Enum:
     2842                        this->FindParam(&calvinghaf,CalvingHeightAboveFloatationEnum);
     2843                        lsf_slopex_input  = this->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
     2844                        if(dim==2) lsf_slopey_input  = this->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
     2845                        calvingrate_input = this->GetInput(CalvingCalvingrateEnum);     _assert_(calvingrate_input);
     2846                        meltingrate_input = this->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
     2847                        break;
     2848                default:
     2849                        _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
     2850        }
     2851
     2852        /* Start looping on the number of vertices: */
     2853        GaussPenta* gauss=new GaussPenta();
     2854        for(int iv=0;iv<NUMVERTICES;iv++){
     2855                gauss->GaussVertex(iv);
     2856
     2857                /* Advection */
     2858                vx_input->GetInputValue(&v[0],gauss);
     2859                vy_input->GetInputValue(&v[1],gauss);
     2860                gr_input->GetInputValue(&groundedice,gauss);
     2861
     2862                /*Get calving speed*/
     2863                switch(calvinglaw){
     2864                        case DefaultCalvingEnum:
     2865                        case CalvingVonmisesEnum:
     2866                                lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
     2867                                if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
     2868                                calvingrate_input->GetInputValue(&calvingrate,gauss);
     2869                                meltingrate_input->GetInputValue(&meltingrate,gauss);
     2870
     2871                                if(groundedice<0) meltingrate = 0.;
     2872
     2873                                norm_dlsf=0.;
     2874                                for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
     2875                                norm_dlsf=sqrt(norm_dlsf);
     2876
     2877                                if(norm_dlsf>1.e-10)
     2878                                 for(i=0;i<dim;i++){
     2879                                         c[i]=calvingrate*dlsf[i]/norm_dlsf; m[i]=meltingrate*dlsf[i]/norm_dlsf;
     2880                                 }
     2881                                else
     2882                                 for(i=0;i<dim;i++){
     2883                                         c[i]=0.; m[i]=0.;
     2884                                 }
     2885                                break;
     2886
     2887                        case CalvingLevermannEnum:
     2888                                calvingratex_input->GetInputValue(&c[0],gauss);
     2889                                if(dim==2) calvingratey_input->GetInputValue(&c[1],gauss);
     2890                                meltingrate_input->GetInputValue(&meltingrate,gauss);
     2891                                norm_calving=0.;
     2892                                for(i=0;i<dim;i++) norm_calving+=pow(c[i],2);
     2893                                norm_calving=sqrt(norm_calving)+1.e-14;
     2894                                for(i=0;i<dim;i++) m[i]=meltingrate*c[i]/norm_calving;
     2895                                break;
     2896
     2897                        case CalvingMinthicknessEnum:
     2898                                lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
     2899                                if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
     2900                                meltingrate_input->GetInputValue(&meltingrate,gauss);
     2901
     2902                                norm_dlsf=0.;
     2903                                for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
     2904                                norm_dlsf=sqrt(norm_dlsf);
     2905
     2906                                if(norm_dlsf>1.e-10)
     2907                                 for(i=0;i<dim;i++){
     2908                                         c[i]=0.;
     2909                                         m[i]=meltingrate*dlsf[i]/norm_dlsf;
     2910                                 }
     2911                                else
     2912                                 for(i=0;i<dim;i++){
     2913                                         c[i]=0.;
     2914                                         m[i]=0.;
     2915                                 }
     2916                                break;
     2917
     2918                        case CalvingHabEnum:
     2919                                lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
     2920                                if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
     2921                                meltingrate_input->GetInputValue(&meltingrate,gauss);
     2922
     2923                                norm_dlsf=0.;
     2924                                for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
     2925                                norm_dlsf=sqrt(norm_dlsf);
     2926
     2927                                if(norm_dlsf>1.e-10)
     2928                                 for(i=0;i<dim;i++){
     2929                                         c[i]=0.;
     2930                                         m[i]=meltingrate*dlsf[i]/norm_dlsf;
     2931                                 }
     2932                                else
     2933                                 for(i=0;i<dim;i++){
     2934                                         c[i]=0.;
     2935                                         m[i]=0.;
     2936                                 }
     2937                                break;
     2938
     2939                        case CalvingCrevasseDepthEnum:
     2940                                lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
     2941                                if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
     2942                                meltingrate_input->GetInputValue(&meltingrate,gauss);
     2943
     2944                                if(groundedice<0) meltingrate = 0.;
     2945
     2946                                norm_dlsf=0.;
     2947                                for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
     2948                                norm_dlsf=sqrt(norm_dlsf);
     2949
     2950                                if(norm_dlsf>1.e-10)
     2951                                 for(i=0;i<dim;i++){
     2952                                         c[i]=0.;
     2953                                         m[i]=meltingrate*dlsf[i]/norm_dlsf;
     2954                                 }
     2955                                else
     2956                                 for(i=0;i<dim;i++){
     2957                                         c[i]=0.;
     2958                                         m[i]=0.;
     2959                                 }
     2960                                break;
     2961
     2962                        case CalvingDev2Enum:
     2963                                  {
     2964                                        lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
     2965                                        if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
     2966                                        calvingrate_input->GetInputValue(&calvingrate,gauss);
     2967                                        meltingrate_input->GetInputValue(&meltingrate,gauss);
     2968                                        gr_input->GetInputValue(&groundedice,gauss);
     2969
     2970                                        //idea: no retreat on ice above critical calving height "calvinghaf" . Limit using regularized Heaviside function.
     2971                                        vel=sqrt(v[0]*v[0] + v[1]*v[1]);
     2972                                        haf_eps=10.;
     2973                                        if(groundedice-calvinghaf<=-haf_eps){
     2974                                                // ice floats freely below calvinghaf: calve freely
     2975                                                // undercutting has no effect:
     2976                                                meltingrate=0.;
     2977                                        }
     2978                                        else if(groundedice-calvinghaf>=haf_eps){
     2979                                                // ice is well above calvinghaf -> no calving back, i.e. limit calving rate to ice velocity
     2980                                                calvingrate=min(calvingrate,vel);
     2981                                                // ice is almost grounded: frontal undercutting has maximum effect (do nothing).
     2982                                        }
     2983                                        else{ // ice is close to calvinghaf: smooth transition between limitation and free calving.
     2984                                                //heaviside: 0 for floating, 1 for grounded
     2985                                                heaviside=(groundedice-calvinghaf+haf_eps)/(2.*haf_eps) + sin(PI*(groundedice-calvinghaf)/haf_eps)/(2.*PI);
     2986                                                calvingrate=heaviside*(min(calvingrate,vel)-calvingrate)+calvingrate;
     2987                                                meltingrate=heaviside*meltingrate+0.;
     2988                                        }
     2989
     2990                                        norm_dlsf=0.;
     2991                                        for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
     2992                                        norm_dlsf=sqrt(norm_dlsf);
     2993
     2994                                        if(norm_dlsf>1.e-10)
     2995                                         for(i=0;i<dim;i++){
     2996                                                 c[i]=calvingrate*dlsf[i]/norm_dlsf;
     2997                                                 m[i]=meltingrate*dlsf[i]/norm_dlsf;
     2998                                         }
     2999                                        else
     3000                                         for(i=0;i<dim;i++){
     3001                                                 c[i]=0.;
     3002                                                 m[i]=0.;
     3003                                         }
     3004                                        break;
     3005                                  }
     3006
     3007                        default:
     3008                                _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
     3009                }
     3010
     3011                for(i=0;i<dim;i++) w[i]=v[i]-c[i]-m[i];
     3012
     3013                movingfrontvx[iv] = w[0];
     3014                movingfrontvy[iv] = w[1];               
     3015        }
     3016
     3017        /*Add input*/
     3018        this->AddInput(MovingFrontalVxEnum,&movingfrontvx[0],P1DGEnum);
     3019        this->AddInput(MovingFrontalVyEnum,&movingfrontvy[0],P1DGEnum);
     3020
     3021        this->InputExtrude(MovingFrontalVxEnum,-1);
     3022        this->InputExtrude(MovingFrontalVyEnum,-1);
     3023        /*Clean up and return*/
     3024        delete gauss;
    27473025}
    27483026/*}}}*/

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

    r25752 r25839  
    129129                IssmDouble     Misfit(int modelenum,int observationenum,int weightsenum){_error_("not implemented yet");};
    130130                IssmDouble     MisfitArea(int weightsenum){_error_("not implemented yet");};
     131                void           MovingFrontalVelocity(void);
    131132                Gauss*         NewGauss(void);
    132133                Gauss*         NewGauss(int order);

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

    r25774 r25839  
    34423442        delete gauss;
    34433443        return Jelem;
     3444}
     3445/*}}}*/
     3446void       Tria::MovingFrontalVelocity(void){/*{{{*/
     3447
     3448        int  dim, domaintype, calvinglaw, i;
     3449        IssmDouble v[3],w[3],c[3],m[3],dlsf[3];
     3450        IssmDouble norm_dlsf, norm_calving, calvingrate, meltingrate, groundedice;
     3451        IssmDouble migrationmax, calvinghaf, heaviside, haf_eps;
     3452        IssmDouble  xyz_list[NUMVERTICES][3];
     3453        IssmDouble  movingfrontvx[NUMVERTICES];
     3454        IssmDouble  movingfrontvy[NUMVERTICES];
     3455        IssmDouble  vel;
     3456       
     3457        /* Get node coordinates and dof list: */
     3458        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     3459
     3460        Input* vx_input           = NULL;
     3461        Input* vy_input           = NULL;
     3462        Input* calvingratex_input = NULL;
     3463        Input* calvingratey_input = NULL;
     3464        Input* lsf_slopex_input   = NULL;
     3465        Input* lsf_slopey_input   = NULL;
     3466        Input* calvingrate_input  = NULL;
     3467        Input* meltingrate_input  = NULL;
     3468        Input* gr_input           = NULL;
     3469
     3470        /*Get problem dimension and whether there is moving front or not*/
     3471        this->FindParam(&domaintype,DomainTypeEnum);
     3472        this->FindParam(&calvinglaw,CalvingLawEnum);
     3473
     3474        switch(domaintype){
     3475                case Domain2DverticalEnum:   dim = 1; break;
     3476                case Domain2DhorizontalEnum: dim = 2; break;
     3477                case Domain3DEnum:           dim = 2; break;
     3478                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     3479        }
     3480        /*Load velocities*/
     3481        switch(domaintype){
     3482                case Domain2DverticalEnum:
     3483                        vx_input=this->GetInput(VxEnum); _assert_(vx_input);
     3484                        break;
     3485                case Domain2DhorizontalEnum:
     3486                        vx_input=this->GetInput(VxEnum); _assert_(vx_input);
     3487                        vy_input=this->GetInput(VyEnum); _assert_(vy_input);
     3488                        gr_input=this->GetInput(MaskOceanLevelsetEnum); _assert_(gr_input);
     3489                        break;
     3490                case Domain3DEnum:
     3491                        vx_input=this->GetInput(VxAverageEnum); _assert_(vx_input);
     3492                        vy_input=this->GetInput(VyAverageEnum); _assert_(vy_input);
     3493                        gr_input=this->GetInput(MaskOceanLevelsetEnum); _assert_(gr_input);
     3494                        break;
     3495                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     3496        }
     3497
     3498        switch(calvinglaw){
     3499                case DefaultCalvingEnum:
     3500                case CalvingVonmisesEnum:
     3501                        lsf_slopex_input  = this->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
     3502                        if(dim==2) lsf_slopey_input  = this->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
     3503                        calvingrate_input = this->GetInput(CalvingCalvingrateEnum);     _assert_(calvingrate_input);
     3504                        meltingrate_input = this->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
     3505                        break;
     3506                case CalvingLevermannEnum:
     3507                        switch(domaintype){
     3508                                case Domain2DverticalEnum:
     3509                                        calvingratex_input=this->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
     3510                                        break;
     3511                                case Domain2DhorizontalEnum:
     3512                                        calvingratex_input=this->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
     3513                                        calvingratey_input=this->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
     3514                                        break;
     3515                                case Domain3DEnum:
     3516                                        calvingratex_input=this->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
     3517                                        calvingratey_input=this->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
     3518                                        break;
     3519                                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     3520                        }
     3521                        meltingrate_input = this->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
     3522                        break;
     3523                case CalvingMinthicknessEnum:
     3524                        lsf_slopex_input  = this->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
     3525                        if(dim==2) lsf_slopey_input  = this->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
     3526                        meltingrate_input = this->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
     3527                        break;
     3528                case CalvingHabEnum:
     3529                        lsf_slopex_input  = this->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
     3530                        if(dim==2) lsf_slopey_input  = this->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
     3531                        meltingrate_input = this->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
     3532                        break;
     3533                case CalvingCrevasseDepthEnum:
     3534                        lsf_slopex_input  = this->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
     3535                        if(dim==2) lsf_slopey_input  = this->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
     3536                        meltingrate_input = this->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
     3537                        break;
     3538                case CalvingDev2Enum:
     3539                        this->FindParam(&calvinghaf,CalvingHeightAboveFloatationEnum);
     3540                        lsf_slopex_input  = this->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
     3541                        if(dim==2) lsf_slopey_input  = this->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
     3542                        calvingrate_input = this->GetInput(CalvingCalvingrateEnum);     _assert_(calvingrate_input);
     3543                        meltingrate_input = this->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
     3544                        break;
     3545                default:
     3546                        _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
     3547        }
     3548
     3549        /* Start looping on the number of vertices: */
     3550        GaussTria* gauss=new GaussTria();
     3551        for(int iv=0;iv<NUMVERTICES;iv++){
     3552                gauss->GaussVertex(iv);
     3553
     3554                /* Advection */
     3555                vx_input->GetInputValue(&v[0],gauss);
     3556                vy_input->GetInputValue(&v[1],gauss);
     3557                gr_input->GetInputValue(&groundedice,gauss);
     3558
     3559                /*Get calving speed*/
     3560                switch(calvinglaw){
     3561                        case DefaultCalvingEnum:
     3562                        case CalvingVonmisesEnum:
     3563                                lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
     3564                                if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
     3565                                calvingrate_input->GetInputValue(&calvingrate,gauss);
     3566                                meltingrate_input->GetInputValue(&meltingrate,gauss);
     3567                                if(groundedice<0) meltingrate = 0.;
     3568
     3569                                norm_dlsf=0.;
     3570                                for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
     3571                                norm_dlsf=sqrt(norm_dlsf);
     3572
     3573                                if(norm_dlsf>1.e-10)
     3574                                 for(i=0;i<dim;i++){
     3575                                         c[i]=calvingrate*dlsf[i]/norm_dlsf; m[i]=meltingrate*dlsf[i]/norm_dlsf;
     3576                                 }
     3577                                else
     3578                                 for(i=0;i<dim;i++){
     3579                                         c[i]=0.; m[i]=0.;
     3580                                 }
     3581                                break;
     3582
     3583                        case CalvingLevermannEnum:
     3584                                calvingratex_input->GetInputValue(&c[0],gauss);
     3585                                if(dim==2) calvingratey_input->GetInputValue(&c[1],gauss);
     3586                                meltingrate_input->GetInputValue(&meltingrate,gauss);
     3587                                norm_calving=0.;
     3588                                for(i=0;i<dim;i++) norm_calving+=pow(c[i],2);
     3589                                norm_calving=sqrt(norm_calving)+1.e-14;
     3590                                for(i=0;i<dim;i++) m[i]=meltingrate*c[i]/norm_calving;
     3591                                break;
     3592
     3593                        case CalvingMinthicknessEnum:
     3594                                lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
     3595                                if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
     3596                                meltingrate_input->GetInputValue(&meltingrate,gauss);
     3597
     3598                                norm_dlsf=0.;
     3599                                for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
     3600                                norm_dlsf=sqrt(norm_dlsf);
     3601
     3602                                if(norm_dlsf>1.e-10)
     3603                                 for(i=0;i<dim;i++){
     3604                                         c[i]=0.;
     3605                                         m[i]=meltingrate*dlsf[i]/norm_dlsf;
     3606                                 }
     3607                                else
     3608                                 for(i=0;i<dim;i++){
     3609                                         c[i]=0.;
     3610                                         m[i]=0.;
     3611                                 }
     3612                                break;
     3613
     3614                        case CalvingHabEnum:
     3615                                lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
     3616                                if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
     3617                                meltingrate_input->GetInputValue(&meltingrate,gauss);
     3618
     3619                                norm_dlsf=0.;
     3620                                for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
     3621                                norm_dlsf=sqrt(norm_dlsf);
     3622
     3623                                if(norm_dlsf>1.e-10)
     3624                                 for(i=0;i<dim;i++){
     3625                                         c[i]=0.;
     3626                                         m[i]=meltingrate*dlsf[i]/norm_dlsf;
     3627                                 }
     3628                                else
     3629                                 for(i=0;i<dim;i++){
     3630                                         c[i]=0.;
     3631                                         m[i]=0.;
     3632                                 }
     3633                                break;
     3634
     3635                        case CalvingCrevasseDepthEnum:
     3636                                lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
     3637                                if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
     3638                                meltingrate_input->GetInputValue(&meltingrate,gauss);
     3639
     3640                                if(groundedice<0) meltingrate = 0.;
     3641
     3642                                norm_dlsf=0.;
     3643                                for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
     3644                                norm_dlsf=sqrt(norm_dlsf);
     3645
     3646                                if(norm_dlsf>1.e-10)
     3647                                 for(i=0;i<dim;i++){
     3648                                         c[i]=0.;
     3649                                         m[i]=meltingrate*dlsf[i]/norm_dlsf;
     3650                                 }
     3651                                else
     3652                                 for(i=0;i<dim;i++){
     3653                                         c[i]=0.;
     3654                                         m[i]=0.;
     3655                                 }
     3656                                break;
     3657
     3658                        case CalvingDev2Enum:
     3659                                  {
     3660                                        lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
     3661                                        if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
     3662                                        calvingrate_input->GetInputValue(&calvingrate,gauss);
     3663                                        meltingrate_input->GetInputValue(&meltingrate,gauss);
     3664                                        gr_input->GetInputValue(&groundedice,gauss);
     3665
     3666                                        //idea: no retreat on ice above critical calving height "calvinghaf" . Limit using regularized Heaviside function.
     3667                                        vel=sqrt(v[0]*v[0] + v[1]*v[1]);
     3668                                        haf_eps=10.;
     3669                                        if(groundedice-calvinghaf<=-haf_eps){
     3670                                                // ice floats freely below calvinghaf: calve freely
     3671                                                // undercutting has no effect:
     3672                                                meltingrate=0.;
     3673                                        }
     3674                                        else if(groundedice-calvinghaf>=haf_eps){
     3675                                                // ice is well above calvinghaf -> no calving back, i.e. limit calving rate to ice velocity
     3676                                                calvingrate=min(calvingrate,vel);
     3677                                                // ice is almost grounded: frontal undercutting has maximum effect (do nothing).
     3678                                        }
     3679                                        else{ // ice is close to calvinghaf: smooth transition between limitation and free calving.
     3680                                                //heaviside: 0 for floating, 1 for grounded
     3681                                                heaviside=(groundedice-calvinghaf+haf_eps)/(2.*haf_eps) + sin(PI*(groundedice-calvinghaf)/haf_eps)/(2.*PI);
     3682                                                calvingrate=heaviside*(min(calvingrate,vel)-calvingrate)+calvingrate;
     3683                                                meltingrate=heaviside*meltingrate+0.;
     3684                                        }
     3685
     3686                                        norm_dlsf=0.;
     3687                                        for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
     3688                                        norm_dlsf=sqrt(norm_dlsf);
     3689
     3690                                        if(norm_dlsf>1.e-10)
     3691                                         for(i=0;i<dim;i++){
     3692                                                 c[i]=calvingrate*dlsf[i]/norm_dlsf;
     3693                                                 m[i]=meltingrate*dlsf[i]/norm_dlsf;
     3694                                         }
     3695                                        else
     3696                                         for(i=0;i<dim;i++){
     3697                                                 c[i]=0.;
     3698                                                 m[i]=0.;
     3699                                         }
     3700                                        break;
     3701                                  }
     3702
     3703                        default:
     3704                                _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
     3705                }
     3706                for(i=0;i<dim;i++) w[i]=v[i]-c[i]-m[i];
     3707                               
     3708                movingfrontvx[iv] = w[0];
     3709                movingfrontvy[iv] = w[1];               
     3710        }
     3711
     3712        /*Add input*/
     3713        this->AddInput(MovingFrontalVxEnum,&movingfrontvx[0],P1DGEnum);
     3714        this->AddInput(MovingFrontalVyEnum,&movingfrontvy[0],P1DGEnum);
     3715
     3716        /*Clean up and return*/
     3717        delete gauss;
    34443718}
    34453719/*}}}*/

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

    r25752 r25839  
    202202                void           JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
    203203                IssmDouble     MinEdgeLength(IssmDouble* xyz_list){_error_("not implemented yet");};
     204                void           MovingFrontalVelocity(void);
    204205                Gauss*         NewGauss(void);
    205206                Gauss*         NewGauss(int order);

  • issm/trunk-jpl/src/c/cores/movingfront_core.cpp

    r25680 r25839  
    111111        }
    112112
     113        /* Calculate the frontal velocity for levelset function */
     114        MovingFrontalVelx(femmodel);
    113115
    114116        /* solve level set equation */

  • issm/trunk-jpl/src/c/modules/modules.h

    r25627 r25839  
    7272#include "./SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.h"
    7373#include "./ModelProcessorx/ModelProcessorx.h"
     74#include "./MovingFrontalVelx/MovingFrontalVelx.h"
    7475#include "./ParseToolkitsOptionsx/ParseToolkitsOptionsx.h"
    7576#include "./NodalValuex/NodalValuex.h"

  • issm/trunk-jpl/src/c/shared/Enum/Enum.vim

    r25807 r25839  
    689689syn keyword cConstant MeshVertexonsurfaceEnum
    690690syn keyword cConstant MisfitEnum
     691syn keyword cConstant MovingFrontalVxEnum
     692syn keyword cConstant MovingFrontalVyEnum
    691693syn keyword cConstant NeumannfluxEnum
    692694syn keyword cConstant NewDamageEnum

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r25807 r25839  
    685685        MeshVertexonsurfaceEnum,
    686686        MisfitEnum,
     687        MovingFrontalVxEnum,
     688        MovingFrontalVyEnum,
    687689        NeumannfluxEnum,
    688690        NewDamageEnum,

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r25807 r25839  
    691691                case MeshVertexonsurfaceEnum : return "MeshVertexonsurface";
    692692                case MisfitEnum : return "Misfit";
     693                case MovingFrontalVxEnum : return "MovingFrontalVx";
     694                case MovingFrontalVyEnum : return "MovingFrontalVy";
    693695                case NeumannfluxEnum : return "Neumannflux";
    694696                case NewDamageEnum : return "NewDamage";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r25807 r25839  
    706706              else if (strcmp(name,"MeshVertexonsurface")==0) return MeshVertexonsurfaceEnum;
    707707              else if (strcmp(name,"Misfit")==0) return MisfitEnum;
     708              else if (strcmp(name,"MovingFrontalVx")==0) return MovingFrontalVxEnum;
     709              else if (strcmp(name,"MovingFrontalVy")==0) return MovingFrontalVyEnum;
    708710              else if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum;
    709711              else if (strcmp(name,"NewDamage")==0) return NewDamageEnum;
    … …  
    750752              else if (strcmp(name,"SealevelriseG")==0) return SealevelriseGEnum;
    751753              else if (strcmp(name,"SealevelriseGU")==0) return SealevelriseGUEnum;
    752               else if (strcmp(name,"SealevelriseGE")==0) return SealevelriseGEEnum;
    753               else if (strcmp(name,"SealevelriseGN")==0) return SealevelriseGNEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum;
     757              if (strcmp(name,"SealevelriseGE")==0) return SealevelriseGEEnum;
     758              else if (strcmp(name,"SealevelriseGN")==0) return SealevelriseGNEnum;
     759              else if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum;
    758760              else if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum;
    759761              else if (strcmp(name,"SedimentHeadSubstep")==0) return SedimentHeadSubstepEnum;
    … …  
    873875              else if (strcmp(name,"StressMaxPrincipal")==0) return StressMaxPrincipalEnum;
    874876              else if (strcmp(name,"StressTensorxx")==0) return StressTensorxxEnum;
    875               else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum;
    876               else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum;
     880              if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum;
     881              else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum;
     882              else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum;
    881883              else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum;
    882884              else if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum;
    … …  
    996998              else if (strcmp(name,"Outputdefinition69")==0) return Outputdefinition69Enum;
    997999              else if (strcmp(name,"Outputdefinition6")==0) return Outputdefinition6Enum;
    998               else if (strcmp(name,"Outputdefinition70")==0) return Outputdefinition70Enum;
    999               else if (strcmp(name,"Outputdefinition71")==0) return Outputdefinition71Enum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"Outputdefinition72")==0) return Outputdefinition72Enum;
     1003              if (strcmp(name,"Outputdefinition70")==0) return Outputdefinition70Enum;
     1004              else if (strcmp(name,"Outputdefinition71")==0) return Outputdefinition71Enum;
     1005              else if (strcmp(name,"Outputdefinition72")==0) return Outputdefinition72Enum;
    10041006              else if (strcmp(name,"Outputdefinition73")==0) return Outputdefinition73Enum;
    10051007              else if (strcmp(name,"Outputdefinition74")==0) return Outputdefinition74Enum;
    … …  
    11191121              else if (strcmp(name,"EsaTransitions")==0) return EsaTransitionsEnum;
    11201122              else if (strcmp(name,"ExternalResult")==0) return ExternalResultEnum;
    1121               else if (strcmp(name,"ExtrapolationAnalysis")==0) return ExtrapolationAnalysisEnum;
    1122               else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum;
     1126              if (strcmp(name,"ExtrapolationAnalysis")==0) return ExtrapolationAnalysisEnum;
     1127              else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum;
     1128              else if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum;
    11271129              else if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum;
    11281130              else if (strcmp(name,"FSSolver")==0) return FSSolverEnum;
    … …  
    12421244              else if (strcmp(name,"MaxDivergence")==0) return MaxDivergenceEnum;
    12431245              else if (strcmp(name,"MaxVel")==0) return MaxVelEnum;
    1244               else if (strcmp(name,"MaxVx")==0) return MaxVxEnum;
    1245               else if (strcmp(name,"MaxVy")==0) return MaxVyEnum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"MaxVz")==0) return MaxVzEnum;
     1249              if (strcmp(name,"MaxVx")==0) return MaxVxEnum;
     1250              else if (strcmp(name,"MaxVy")==0) return MaxVyEnum;
     1251              else if (strcmp(name,"MaxVz")==0) return MaxVzEnum;
    12501252              else if (strcmp(name,"Melange")==0) return MelangeEnum;
    12511253              else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum;
    … …  
    13651367              else if (strcmp(name,"SubelementMelt1")==0) return SubelementMelt1Enum;
    13661368              else if (strcmp(name,"SubelementMelt2")==0) return SubelementMelt2Enum;
    1367               else if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum;
    1368               else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum;
    13691369         else stage=12;
    13701370   }
    13711371   if(stage==12){
    1372               if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum;
     1372              if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum;
     1373              else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum;
     1374              else if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum;
    13731375              else if (strcmp(name,"Tetra")==0) return TetraEnum;
    13741376              else if (strcmp(name,"TetraInput")==0) return TetraInputEnum;

  • issm/trunk-jpl/test/NightlyRun/test540.m

    r23829 r25839  
    1616pos = find(md.mesh.vertexonboundary);
    1717md.levelset.spclevelset(pos) = md.mask.ice_levelset(pos);
     18md.levelset.migration_max = 1e10;
    1819
    1920%Force MUMPS sequential analysis

  • issm/trunk-jpl/test/NightlyRun/test541.m

    r23688 r25839  
    1616pos = find(md.mesh.vertexonboundary);
    1717md.levelset.spclevelset(pos) = md.mask.ice_levelset(pos);
     18md.levelset.migration_max = 1e10;
    1819
    1920%Force MUMPS sequential analysis

  • issm/trunk-jpl/test/NightlyRun/test804.m

    r23652 r25839  
    1717md.calving.calvingrate=1000.*ones(md.mesh.numberofvertices,1);
    1818md.frontalforcings.meltingrate=zeros(md.mesh.numberofvertices,1);
     19md.levelset.migration_max = 1e10;
    1920
    2021md=solve(md,'Transient');

  • issm/trunk-jpl/test/NightlyRun/test805.m

    r23652 r25839  
    2626md.groundingline.melt_interpolation='SubelementMelt1';
    2727md.levelset.stabilization=2;
     28md.levelset.migration_max = 1e10;
    2829
    2930md=solve(md,'Transient');

  • issm/trunk-jpl/test/NightlyRun/test806.m

    r24043 r25839  
    3232md.frontalforcings.meltingrate=zeros(md.mesh.numberofvertices,1);
    3333md.levelset.spclevelset=NaN(md.mesh.numberofvertices,1);
     34md.levelset.migration_max = 1e8;
    3435
    3536md.transient.requested_outputs={'default','StrainRateparallel','StrainRateperpendicular','Calvingratex','Calvingratey','CalvingCalvingrate'};

  • issm/trunk-jpl/test/NightlyRun/test808.m

    r24043 r25839  
    3333md.levelset.spclevelset=NaN(md.mesh.numberofvertices,1);
    3434md.levelset.reinit_frequency=1;
     35md.levelset.migration_max = 1e10;
    3536
    3637md=solve(md,'Transient');

  • issm/trunk-jpl/test/NightlyRun/test809.m

    r23652 r25839  
    2323md.levelset.spclevelset=NaN(md.mesh.numberofvertices,1);
    2424md.levelset.reinit_frequency=1;
     25md.levelset.migration_max = 1e10;
    2526
    2627md=solve(md,'Transient');
Note: See TracChangeset for help on using the changeset viewer.