Changeset 26463


Ignore:
Timestamp:
09/29/21 07:15:33 (3 years ago)
Author:
Cheng Gong
Message:

CHG: use Friction.cpp to take care of the choice of basal velocity, set the pointers of the inputs as public attributes of Friction, provide functions to fetch the basal velocities that used for calculating the friciton

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

Legend:

Unmodified
Added
Removed

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

    r26329 r26463  
    409409        Input* pressure_input       = element->GetInput(PressureEnum);                    _assert_(pressure_input);
    410410        Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
    411         Input* vx_input             = element->GetInput(VxEnum);                          _assert_(vx_input);
    412         Input* vy_input             = element->GetInput(VyEnum);                          _assert_(vy_input);
    413         Input* vz_input             = element->GetInput(VzEnum);                          _assert_(vz_input);
    414411
    415412        /*Build friction element, needed later: */
    416         Friction* friction=new Friction(element,dim);
     413        Friction* friction=new Friction(element,3);
    417414
    418415        /******** MELTING RATES  ************************************//*{{{*/
    … …  
    470467                        /*basal friction*/
    471468                        friction->GetAlpha2(&alpha2,gauss);
    472                         vx_input->GetInputValue(&vx,gauss);             vy_input->GetInputValue(&vy,gauss);             vz_input->GetInputValue(&vz,gauss);
     469                        friction->GetBasalSlidingSpeeds(&vx, &vy, &vz, gauss);
    473470                        basalfriction=alpha2*(vx*vx + vy*vy + vz*vz);
    474471                        geothermalflux_input->GetInputValue(&geothermalflux,gauss);
    … …  
    947944        if(dt==0. && !converged) enthalpy_enum=EnthalpyPicardEnum; // use enthalpy from last iteration
    948945        else enthalpy_enum=EnthalpyEnum; // use enthalpy from last time step
    949         Input* vx_input             = element->GetInput(VxEnum);                          _assert_(vx_input);
    950         Input* vy_input             = element->GetInput(VyEnum);                          _assert_(vy_input);
    951         Input* vz_input             = element->GetInput(VzEnum);                          _assert_(vz_input);
    952946        Input* enthalpy_input            = element->GetInput(enthalpy_enum);                                     _assert_(enthalpy_input);
    953947        Input* pressure_input            = element->GetInput(PressureEnum);                                                      _assert_(pressure_input);
    … …  
    988982                                geothermalflux_input->GetInputValue(&geothermalflux,gauss);
    989983                                friction->GetAlpha2(&alpha2,gauss);
    990                                 vx_input->GetInputValue(&vx,gauss);
    991                                 vy_input->GetInputValue(&vy,gauss);
    992                                 vz_input->GetInputValue(&vz,gauss);
     984                                friction->GetBasalSlidingSpeeds(&vx, &vy, &vz, gauss);
    993985                                basalfriction=alpha2*(vx*vx+vy*vy+vz*vz);
    994986                                heatflux=(basalfriction+geothermalflux)/(rho_ice);

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

    r26460 r26463  
    360360        IssmDouble e_v       = element->FindParam(HydrologyEnglacialVoidRatioEnum);
    361361        Input* hr_input     = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input);
    362         Input* vx_input     = element->GetInput(VxBaseEnum);_assert_(vx_input);
    363         Input* vy_input     = element->GetInput(VyBaseEnum);_assert_(vy_input);
    364362        Input* h_input      = element->GetInput(HydrologySheetThicknessEnum);_assert_(h_input);
    365363        Input* H_input      = element->GetInput(ThicknessEnum); _assert_(H_input);
    … …  
    383381
    384382                /*Get input values at gauss points*/
    385                 vx_input->GetInputValue(&vx,gauss);
    386                 vy_input->GetInputValue(&vy,gauss);
    387383                h_input->GetInputValue(&h,gauss);
    388384                G_input->GetInputValue(&G,gauss);
    … …  
    395391
    396392                /*Get basal velocity*/
     393                friction->GetBasalSlidingSpeeds(&vx, &vy ,gauss);
    397394                ub = sqrt(vx*vx + vy*vy);
    398395

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

    r26329 r26463  
    119119        iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum);
    120120        iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyEnum);
     121        if(iomodel->domaintype==Domain2DhorizontalEnum){
     122                iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxBaseEnum);
     123                iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyBaseEnum);
     124        }
    121125        iomodel->FindConstant(&frictionlaw,"md.friction.law");
    122126
    … …  
    254258        Input* n_input              = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
    255259        Input* englacial_input      = element->GetInput(HydrologyEnglacialInputEnum);    _assert_(englacial_input);
    256         Input* vx_input             = element->GetInput(VxEnum);                         _assert_(vx_input);
    257         Input* vy_input             = element->GetInput(VyEnum);                         _assert_(vy_input);
    258260        Input* lr_input             = element->GetInput(HydrologyBumpSpacingEnum);       _assert_(lr_input);
    259261        Input* br_input             = element->GetInput(HydrologyBumpHeightEnum);        _assert_(br_input);
    … …  
    287289                lr_input->GetInputValue(&lr,gauss);
    288290                br_input->GetInputValue(&br,gauss);
    289                 vx_input->GetInputValue(&vx,gauss);
    290                 vy_input->GetInputValue(&vy,gauss);
    291291      headold_input->GetInputValue(&head_old,gauss);
    292292
    … …  
    304304                /*Compute frictional heat flux*/
    305305                friction->GetAlpha2(&alpha2,gauss);
    306                 vx_input->GetInputValue(&vx,gauss);
    307                 vy_input->GetInputValue(&vy,gauss);
     306                friction->GetBasalSlidingSpeeds(&vx, &vy, gauss);
    308307                frictionheat=alpha2*(vx*vx+vy*vy);
    309308
    … …  
    494493        Input* B_input              = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
    495494        Input* n_input              = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
    496         Input* vx_input             = element->GetInput(VxEnum);                         _assert_(vx_input);
    497         Input* vy_input             = element->GetInput(VyEnum);                         _assert_(vy_input);
    498495        Input* lr_input             = element->GetInput(HydrologyBumpSpacingEnum);       _assert_(lr_input);
    499496        Input* br_input             = element->GetInput(HydrologyBumpHeightEnum);        _assert_(br_input);
    … …  
    523520                lr_input->GetInputValue(&lr,gauss);
    524521                br_input->GetInputValue(&br,gauss);
    525                 vx_input->GetInputValue(&vx,gauss);
    526                 vy_input->GetInputValue(&vy,gauss);
    527522
    528523                /*Get ice A parameter*/
    … …  
    539534                /*Compute frictional heat flux*/
    540535                friction->GetAlpha2(&alpha2,gauss);
    541                 vx_input->GetInputValue(&vx,gauss);
    542                 vy_input->GetInputValue(&vy,gauss);
     536                friction->GetBasalSlidingSpeeds(&vx, &vy, gauss);
    543537                frictionheat=alpha2*(vx*vx+vy*vy);
    544538

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

    r26457 r26463  
    34893489        /*build friction object, used later on: */
    34903490        /*dim=4 is special for HO, which is actually 2.5D*/
    3491         Friction* friction=new Friction(element,dim==3?4:1);
     3491        Friction* friction=new Friction(element,dim==3?2.5:1);
    34923492
    34933493        /*Recover portion of element that is grounded*/

  • issm/trunk-jpl/src/c/classes/Loads/Friction.cpp

    r26457 r26463  
    1818Friction::Friction(){/*{{{*/
    1919        this->element=NULL;
    20         this->dim=0;
    2120        this->law=0;
     21        this->apply_dim = 1;
     22        this->domaintype=-1;
     23        this->vx_input=NULL;
     24        this->vy_input=NULL;
     25        this->vz_input=NULL;
    2226}
    2327/*}}}*/
    24 Friction::Friction(Element* element_in,int dim_in){/*{{{*/
    25 
     28Friction::Friction(Element* element_in){/*{{{*/
     29        /* Determine the dimension according to the domain type automatically.
     30         * There are exceptions, e.g. HO, which needs the user to specify the dimension used in Friciton.*/
    2631        this->element=element_in;
    27         this->dim=dim_in;
     32
     33        /* Load necessary parameters */
    2834        element_in->FindParam(&this->law,FrictionLawEnum);
     35        element_in->FindParam(&this->domaintype,DomainTypeEnum);
     36       
     37        /* Load VxBase and VyBase for this special case */
     38        switch(this->domaintype){
     39                case Domain2DhorizontalEnum:
     40                        this->apply_dim = 2;
     41                        this->vx_input = element_in->GetInput(VxBaseEnum);                      _assert_(this->vx_input);
     42                        this->vy_input = element_in->GetInput(VyBaseEnum);                      _assert_(this->vy_input);
     43                        this->vz_input = NULL;
     44                        break;
     45      case Domain2DverticalEnum:
     46                        this->apply_dim = 2;
     47                        this->vx_input = element_in->GetInput(VxEnum);                          _assert_(this->vx_input);
     48                        this->vy_input = element_in->GetInput(VyEnum);                          _assert_(this->vy_input);
     49                        this->vz_input = NULL;
     50                        break;
     51      case Domain3DEnum:           
     52                        this->apply_dim = 3;
     53                        this->vx_input = element_in->GetInput(VxEnum);                          _assert_(this->vx_input);
     54                        this->vy_input = element_in->GetInput(VyEnum);                          _assert_(this->vy_input);
     55                        this->vz_input = element_in->GetInput(VzEnum);                          _assert_(this->vz_input);
     56                        break;
     57      default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     58        }
     59}
     60/*}}}*/
     61Friction::Friction(Element* element_in, double dim){/*{{{*/
     62
     63        this->element=element_in;
     64        this->apply_dim = dim;
     65
     66        /* Load necessary parameters */
     67        element_in->FindParam(&this->law,FrictionLawEnum);
     68        element_in->FindParam(&this->domaintype,DomainTypeEnum);
     69       
     70        /* Load VxBase and VyBase for this special case */
     71        switch(this->domaintype){
     72                case Domain2DhorizontalEnum:
     73                        this->vx_input = element_in->GetInput(VxBaseEnum);                      _assert_(this->vx_input);
     74                        this->vy_input = element_in->GetInput(VyBaseEnum);                      _assert_(this->vy_input);
     75                        this->vz_input = NULL;
     76                        break;
     77      case Domain2DverticalEnum:
     78                        this->vx_input = element_in->GetInput(VxEnum);                          _assert_(this->vx_input);
     79                        this->vy_input = element_in->GetInput(VyEnum);                          _assert_(this->vy_input);
     80                        this->vz_input = NULL;
     81                        break;
     82      case Domain3DEnum:           
     83                        this->vx_input = element_in->GetInput(VxEnum);                          _assert_(this->vx_input);
     84                        this->vy_input = element_in->GetInput(VyEnum);                          _assert_(this->vy_input);
     85                        this->vz_input = element_in->GetInput(VzEnum);                          _assert_(this->vz_input);
     86                        break;
     87      default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     88        }
    2989}
    3090/*}}}*/
    … …  
    3696void Friction::Echo(void){/*{{{*/
    3797        _printf_("Friction:\n");
    38         _printf_("   dim: " << this->dim<< "\n");
     98        _printf_("   Domain type: " << this->domaintype<< "\n");
    3999}
    40100/*}}}*/
    … …  
    855915        /*diverse*/
    856916        IssmDouble vx,vy,vz,vmag;
     917
     918        this->vx_input->GetInputValue(&vx, gauss);
     919        this->vy_input->GetInputValue(&vy, gauss);
     920
     921        if ((this->vz_input == NULL) || (this->apply_dim<3)) vz = 0.0;
     922        else this->vz_input->GetInputValue(&vz, gauss);
    857923       
    858         switch(dim){
    859                 case 1:
    860                         element->GetInputValue(&vx,gauss,VxEnum);
    861                         vmag=sqrt(vx*vx);
    862                         break;
    863                 case 2:
    864                         element->GetInputValue(&vx,gauss,VxBaseEnum);
    865                         element->GetInputValue(&vy,gauss,VyBaseEnum);
    866                         vmag=sqrt(vx*vx+vy*vy);
    867                         break;
    868                 case 3:
    869                         element->GetInputValue(&vx,gauss,VxEnum);
    870                         element->GetInputValue(&vy,gauss,VyEnum);
    871                         element->GetInputValue(&vz,gauss,VzEnum);
    872                         vmag=sqrt(vx*vx+vy*vy+vz*vz);
    873                         break;
    874                 case 4:
    875                         /* This is for HO 3D case, since it requires the horizontal velocities */
    876                         element->GetInputValue(&vx,gauss,VxEnum);
    877                         element->GetInputValue(&vy,gauss,VyEnum);
    878                         vmag=sqrt(vx*vx+vy*vy);
    879                         break;
    880                 default:
    881                         _error_("not supported");
    882         }
    883 
     924        if (this->apply_dim<2) vy = 0.0;
     925
     926        vmag = sqrt(vx*vx+vy*vy+vz*vz);
    884927        return vmag;
    885 
    886 }/*}}}*/
     928}/*}}}*/
     929void Friction::GetBasalSlidingSpeeds(IssmDouble* pvx, Gauss* gauss){/*{{{*/
     930
     931        this->vx_input->GetInputValue(pvx, gauss);
     932        /*Checks*/
     933        _assert_(!xIsNan<IssmDouble>(*pvx));
     934        _assert_(!xIsInf<IssmDouble>(*pvx));
     935}/*}}}*/
     936void Friction::GetBasalSlidingSpeeds(IssmDouble* pvx, IssmDouble* pvy, Gauss* gauss){/*{{{*/
     937
     938        this->vx_input->GetInputValue(pvx, gauss);
     939        this->vy_input->GetInputValue(pvy, gauss);
     940        /*Checks*/
     941        _assert_(!xIsNan<IssmDouble>(*pvx));
     942        _assert_(!xIsInf<IssmDouble>(*pvx));
     943        _assert_(!xIsNan<IssmDouble>(*pvy));
     944        _assert_(!xIsInf<IssmDouble>(*pvy));
     945}/*}}}*/
     946void Friction::GetBasalSlidingSpeeds(IssmDouble* pvx, IssmDouble* pvy, IssmDouble* pvz, Gauss* gauss){/*{{{*/
     947
     948        this->vx_input->GetInputValue(pvx, gauss);
     949        this->vy_input->GetInputValue(pvy, gauss);
     950        this->vz_input->GetInputValue(pvz, gauss);
     951        /*Checks*/
     952        _assert_(!xIsNan<IssmDouble>(*pvx));
     953        _assert_(!xIsInf<IssmDouble>(*pvx));
     954        _assert_(!xIsNan<IssmDouble>(*pvy));
     955        _assert_(!xIsInf<IssmDouble>(*pvy));
     956        _assert_(!xIsNan<IssmDouble>(*pvz));
     957        _assert_(!xIsInf<IssmDouble>(*pvz));
     958}/*}}}*/

  • issm/trunk-jpl/src/c/classes/Loads/Friction.h

    r26457 r26463  
    1515        public:
    1616                Element* element;
    17                 int      dim;
    1817                int      law;
     18                int      domaintype;
     19                double  apply_dim;
     20                Input*  vx_input;
     21                Input*  vy_input;
     22                Input*  vz_input;
    1923
    2024                /*methods: */
    2125                Friction();
    22                 Friction(Element* element_in,int dim_in);
     26                Friction(Element* element_in);
     27                Friction(Element* element_in, double dim);
    2328                ~Friction();
    2429
    … …  
    4752                IssmDouble EffectivePressure(Gauss* gauss);
    4853                IssmDouble VelMag(Gauss* gauss);
     54                void GetBasalSlidingSpeeds(IssmDouble* pvx, Gauss* gauss);
     55                void GetBasalSlidingSpeeds(IssmDouble* pvx, IssmDouble* pvy, Gauss* gauss);
     56                void GetBasalSlidingSpeeds(IssmDouble* pvx, IssmDouble* pvy, IssmDouble* pvz, Gauss* gauss);
    4957};
    5058
Note: See TracChangeset for help on using the changeset viewer.