Changeset 28039


Ignore:
Timestamp:
01/04/24 05:49:30 (15 months ago)
Author:
seroussi
Message:

NEW: extended CalvingMIP to 3d Penta

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

Legend:

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

    r27964 r28039  
    286286        IssmDouble  groundedice,bed,sealevel;
    287287
    288         /*Depth average B for stress calculation*/
     288        /*Depth average velocity for stress calculation*/
    289289        this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
    290290        this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
     
    625625                delete gauss;
    626626        }
     627}
     628/*}}}*/
     629void       Penta::CalvingRateCalvingMIP(){/*{{{*/
     630
     631        if(!this->IsOnBase()) return;
     632
     633        IssmDouble  calvingrate[NUMVERTICES];
     634        IssmDouble  calvingratex[NUMVERTICES];
     635        IssmDouble  calvingratey[NUMVERTICES];
     636        int                     experiment = 1;  /* exp:1 by default */
     637        int         dim, domaintype;
     638        IssmDouble      vx, vy, vel, c, wrate;
     639        IssmDouble  time, groundedice, yts;
     640
     641        /*Get problem dimension and whether there is moving front or not*/
     642        this->FindParam(&domaintype,DomainTypeEnum);
     643        this->FindParam(&time,TimeEnum);
     644        this->FindParam(&yts,ConstantsYtsEnum);
     645
     646        switch(domaintype){
     647                case Domain2DverticalEnum:   dim = 1; break;
     648                case Domain2DhorizontalEnum: dim = 2; break;
     649                case Domain3DEnum:           dim = 2; break;
     650                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     651        }
     652        if(dim==1) _error_("not implemented in 1D...");
     653
     654        /*Depth average velocity for stress calculation*/
     655        this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
     656        this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
     657
     658        /*Retrieve all inputs and parameters we will need*/
     659        Input *vx_input      = this->GetInput(VxAverageEnum);                                _assert_(vx_input);
     660        Input *vy_input      = this->GetInput(VyAverageEnum);                                _assert_(vy_input);
     661        Input *wrate_input   = this->GetInput(CalvingAblationrateEnum);               _assert_(wrate_input);
     662        Input* gr_input      = this->GetInput(MaskOceanLevelsetEnum);                                           _assert_(gr_input);
     663
     664        /* Use which experiment: use existing Enum */
     665        this->FindParam(&experiment, CalvingUseParamEnum);
     666
     667        /* Start looping on the number of vertices: */
     668        GaussPenta gauss;
     669        for(int iv=0;iv<3;iv++){
     670                gauss.GaussVertex(iv);
     671
     672                /*Get velocity components */
     673                vx_input->GetInputValue(&vx,&gauss);
     674                vy_input->GetInputValue(&vy,&gauss);
     675                vel=sqrt(vx*vx+vy*vy)+1.e-14;
     676
     677                /* no calving for grounded ice in EXP4 */
     678                gr_input->GetInputValue(&groundedice,&gauss);
     679
     680                switch (experiment) {
     681                        case 1:
     682                        case 3:
     683                                /* Exp 1 and 3: set c=v-wrate, wrate=0, so that w=0 */
     684                                wrate = 0.0;
     685                                break;
     686                        case 2:
     687                                /* Exp 2: set c=v-wrate(given)*/
     688                                wrate_input->GetInputValue(&wrate,&gauss);
     689                                break;
     690                        case 4:
     691                                /* Exp 4: set c=v-wrate(given), for the first 500 years, then c=0 for the second 500 years*/
     692                                if((groundedice<0) && (time<=500.0*yts)) {
     693                                //      wrate_input->GetInputValue(&wrate,&gauss);
     694                                        wrate = -750*sin(2.0*M_PI*time/yts/1000)/yts;  // m/a -> m/s
     695                                }
     696                                else {
     697                                        /* no calving on the grounded ice*/
     698                                        wrate = vel;
     699                                }
     700                                break;
     701                        default:
     702                                _error_("The experiment is not supported yet!");
     703                }
     704
     705                calvingrate[iv] = vel - wrate;
     706                calvingratex[iv] = vx - wrate*vx/vel;
     707                calvingratey[iv] = vy - wrate*vy/vel;
     708        }
     709        /*Add input*/
     710        this->AddBasalInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum);
     711        this->AddBasalInput(CalvingratexEnum,&calvingratex[0],P1DGEnum);
     712        this->AddBasalInput(CalvingrateyEnum,&calvingratey[0],P1DGEnum);
     713
     714        /*Extrude*/
     715        this->InputExtrude(CalvingCalvingrateEnum,-1);
     716        this->InputExtrude(CalvingratexEnum,-1);
     717        this->InputExtrude(CalvingrateyEnum,-1);
    627718}
    628719/*}}}*/
     
    28102901                case CalvingVonmisesEnum:
    28112902                case CalvingLevermannEnum:
     2903                case CalvingCalvingMIPEnum:
    28122904                        calvingratex_input=this->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
    28132905                        calvingratey_input=this->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
     
    28472939                        case CalvingVonmisesEnum:
    28482940                        case CalvingLevermannEnum:
     2941                        case CalvingCalvingMIPEnum:
    28492942                                calvingratex_input->GetInputValue(&c[0],&gauss);
    28502943                                calvingratey_input->GetInputValue(&c[1],&gauss);
  • issm/trunk-jpl/src/c/classes/Elements/Penta.h

    r27964 r28039  
    5858                void           CalvingFluxLevelset();
    5959                void           CalvingMeltingFluxLevelset();
     60                void           CalvingRateCalvingMIP();
    6061                IssmDouble     CharacteristicLength(void){_error_("not implemented yet");};
    6162                void           ComputeBasalStress(void);
Note: See TracChangeset for help on using the changeset viewer.