Changeset 22847


Ignore:
Timestamp:
06/18/18 07:28:06 (7 years ago)
Author:
Mathieu Morlighem
Message:

CHG: added EffectivePressure function to simplify switch

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

Legend:

Unmodified
Added
Removed

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

    r22841 r22847  
    6060
    6161        /*diverse: */
    62         int         CoupledFlag;
    6362        IssmDouble  q_exp;
    6463        IssmDouble  C_param;
    6564        IssmDouble  As;
    66         IssmDouble  Neff;
    6765        IssmDouble  n;
    6866        IssmDouble  alpha;
    … …  
    7068        IssmDouble  vx,vy,vz,vmag;
    7169        IssmDouble  alpha_complement;
    72         IssmDouble  base,sealevel,thickness;
    7370
    7471        /*Recover parameters: */
    7572        element->GetInputValue(&q_exp,FrictionQEnum);
    7673        element->GetInputValue(&C_param,FrictionCEnum);
    77 
    7874        element->GetInputValue(&As,gauss,FrictionAsEnum);
    7975        element->GetInputValue(&n,gauss,MaterialsRheologyNEnum);
    80         element->GetInputValue(&thickness, gauss,ThicknessEnum);
    81         element->GetInputValue(&base, gauss,BaseEnum);
    82         element->GetInputValue(&sealevel, gauss,SealevelEnum);
    83 
    84         IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    85         IssmDouble rho_ice     = element->GetMaterialParameter(MaterialsRhoIceEnum);
    86         IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);
    87         element->parameters->FindParam(&CoupledFlag,FrictionCouplingEnum);
    88 
    89         switch(CoupledFlag){
    90                 //This should be removed at some point and the Enum uniformalized
    91                 case 0:
    92                         Neff=gravity*(rho_ice*thickness+rho_water*(base-sealevel));
    93                         break;
    94                 case 1:
    95                         element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum);
    96                         break;
    97                 case 2:
    98                         element->GetInputValue(&Neff,gauss,EffectivePressureTimeAverageEnum);
    99                         break;
    100                 default:
    101                         _error_("not supported");
    102         }
    103 
    104         if(Neff<0)Neff=0;
     76
     77        /*Get effective pressure*/
     78        IssmDouble Neff = EffectivePressure(gauss);
    10579
    10680        //We need the velocity magnitude to evaluate the basal stress:
    … …  
    180154
    181155        /*diverse: */
    182         int         CoupledFlag;
    183156        IssmDouble  r,s;
    184157        IssmDouble  vx,vy,vz,vmag;
    185158        IssmDouble  drag_p,drag_q;
    186         IssmDouble  Neff;
    187159        IssmDouble  drag_coefficient;
    188         IssmDouble  base,thickness,sealevel;
    189160        IssmDouble  alpha_complement;
    190161
    … …  
    192163        element->GetInputValue(&drag_p,FrictionPEnum);
    193164        element->GetInputValue(&drag_q,FrictionQEnum);
    194         element->GetInputValue(&thickness, gauss,ThicknessEnum);
    195         element->GetInputValue(&base, gauss,BaseEnum);
    196         element->GetInputValue(&sealevel, gauss,SealevelEnum);
    197165        element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
    198         IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    199         IssmDouble rho_ice     = element->GetMaterialParameter(MaterialsRhoIceEnum);
    200         IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);
    201         element->parameters->FindParam(&CoupledFlag,FrictionCouplingEnum);
    202         //compute r and q coefficients: */
     166
     167        /*compute r and q coefficients: */
    203168        r=drag_q/drag_p;
    204169        s=1./drag_p;
    205170
    206         //From base and thickness, compute effective pressure when drag is viscous, or get Neff from forcing or coupled to hydrologymodel:
    207         switch(CoupledFlag){
    208                 case 0:
    209                         Neff=gravity*(rho_ice*thickness+rho_water*(base-sealevel));
    210                         break;
    211                 case 1:
    212                         element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum);
    213                         break;
    214                 case 2:
    215                         element->GetInputValue(&Neff,gauss,EffectivePressureTimeAverageEnum);
    216                         break;
    217                 default:
    218                         _error_("not supported");
    219         }
    220         if(Neff<0)Neff=0;
     171        /*Get effective pressure*/
     172        IssmDouble Neff = EffectivePressure(gauss);
    221173
    222174        //We need the velocity magnitude to evaluate the basal stress:
    … …  
    290242
    291243        /*diverse: */
    292         int         CoupledFlag;
    293244        IssmDouble  r,s;
    294245        IssmDouble  drag_p, drag_q;
    295         IssmDouble  Neff;
    296246        IssmDouble  thickness,base,floatation_thickness,sealevel;
    297247        IssmDouble  vx,vy,vz,vmag;
    … …  
    307257        element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
    308258        element->GetInputValue(&drag_coefficient_coulomb, gauss,FrictionCoefficientcoulombEnum);
    309         IssmDouble rho_water        = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    310         IssmDouble rho_ice          = element->GetMaterialParameter(MaterialsRhoIceEnum);
    311         IssmDouble gravity          = element->GetMaterialParameter(ConstantsGEnum);
    312         element->parameters->FindParam(&CoupledFlag,FrictionCouplingEnum);
     259        IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
     260        IssmDouble rho_ice   = element->GetMaterialParameter(MaterialsRhoIceEnum);
     261        IssmDouble gravity   = element->GetMaterialParameter(ConstantsGEnum);
    313262
    314263        //compute r and q coefficients: */
    … …  
    316265        s=1./drag_p;
    317266
    318         //From base and thickness, compute effective pressure when drag is viscous:
    319         switch(CoupledFlag){
    320                 case 0:
    321                         Neff=gravity*(rho_ice*thickness+rho_water*(base-sealevel));
    322                         break;
    323                 case 1:
    324                         element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum);
    325                         break;
    326                 case 2:
    327                         element->GetInputValue(&Neff,gauss,EffectivePressureTimeAverageEnum);
    328                         break;
    329                 default:
    330                         _error_("not supported");
    331         }
    332 
    333         if(Neff<0)Neff=0;
    334 
     267
     268        /*Get effective pressure*/
     269        IssmDouble Neff = EffectivePressure(gauss);
     270
     271        /*Compute velocity magnitude*/
    335272        switch(dim){
    336273                case 1:
    … …  
    379316
    380317        /*diverse: */
    381         int         CoupledFlag;
    382318        IssmDouble  q_exp;
    383319        IssmDouble  C_param;
    384320        IssmDouble  As;
    385 
    386         IssmDouble  Neff;
    387321        IssmDouble  n;
    388322
    … …  
    392326        IssmDouble  vx,vy,vz,vmag;
    393327        IssmDouble  alpha2;
    394         IssmDouble  base,thickness,sealevel;
    395328
    396329        /*Recover parameters: */
    … …  
    399332        element->GetInputValue(&As,gauss,FrictionAsEnum);
    400333        element->GetInputValue(&n,gauss,MaterialsRheologyNEnum);
    401         element->GetInputValue(&thickness, gauss,ThicknessEnum);
    402         element->GetInputValue(&base, gauss,BaseEnum);
    403         element->GetInputValue(&sealevel, gauss,SealevelEnum);
    404         IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    405         IssmDouble rho_ice     = element->GetMaterialParameter(MaterialsRhoIceEnum);
    406         IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);
    407         element->parameters->FindParam(&CoupledFlag,FrictionCouplingEnum);
    408 
    409         switch(CoupledFlag){
    410                 case 0:
    411                         Neff=gravity*(rho_ice*thickness+rho_water*(base-sealevel));
    412                         break;
    413                 case 1:
    414                         element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum);
    415                         break;
    416                 case 2:
    417                         element->GetInputValue(&Neff,gauss,EffectivePressureTimeAverageEnum);
    418                         break;
    419                 default:
    420                         _error_("not supported");
    421         }
    422 
    423         if(Neff<0)Neff=0;
    424 
     334
     335        /*Get effective pressure*/
     336        IssmDouble Neff = EffectivePressure(gauss);
     337
     338        /*Compute velocity magnitude*/
    425339        switch(dim){
    426340                case 1:
    … …  
    575489
    576490        /*diverse: */
    577         int         CoupledFlag;
    578491        IssmDouble  r,s;
    579492        IssmDouble  drag_p, drag_q;
    580         IssmDouble  Neff;
    581         IssmDouble  thickness,base,sealevel;
    582493        IssmDouble  vx,vy,vz,vmag;
    583494        IssmDouble  drag_coefficient;
    … …  
    587498        element->GetInputValue(&drag_p,FrictionPEnum);
    588499        element->GetInputValue(&drag_q,FrictionQEnum);
    589         element->GetInputValue(&thickness, gauss,ThicknessEnum);
    590         element->GetInputValue(&base, gauss,BaseEnum);
    591         element->GetInputValue(&sealevel, gauss,SealevelEnum);
    592500        element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
    593         IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    594         IssmDouble rho_ice     = element->GetMaterialParameter(MaterialsRhoIceEnum);
    595         IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);
    596         element->parameters->FindParam(&CoupledFlag,FrictionCouplingEnum);
    597         //compute r and q coefficients: */
     501
     502        /*compute r and q coefficients: */
    598503        r=drag_q/drag_p;
    599504        s=1./drag_p;
    600505
    601         //From base and thickness, compute effective pressure when drag is viscous, or get Neff from forcing:
    602         switch(CoupledFlag){
    603                 case 0:
    604                         Neff=gravity*(rho_ice*thickness+rho_water*(base-sealevel));
    605                         break;
    606                 case 1:
    607                         element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum);
    608                         break;
    609                 case 2:
    610                         element->GetInputValue(&Neff,gauss,EffectivePressureTimeAverageEnum);
    611                         break;
    612                 default:
    613                         _error_("not supported");
    614         }
    615         if(Neff<0)Neff=0;
    616 
     506        /*Get effective pressure*/
     507        IssmDouble Neff = EffectivePressure(gauss);
     508
     509        /*Compute velocity magnitude*/
    617510        switch(dim){
    618511                case 1:
    … …  
    777670        *palpha2=alpha2;
    778671}/*}}}*/
     672IssmDouble Friction::EffectivePressure(Gauss* gauss){/*{{{*/
     673        /*Get effective pressure as a function of
     674         *  - coupling
     675         *  - type
     676         * Neff=rho_ice*g*thickness+rho_ice*g*base
     677         */
     678
     679
     680        /*diverse: */
     681        int         coupled_flag;
     682        IssmDouble  Neff;
     683
     684        /*Recover parameters: */
     685        element->parameters->FindParam(&coupled_flag,FrictionCouplingEnum);
     686
     687        /*From base and thickness, compute effective pressure when drag is viscous, or get Neff from forcing:*/
     688        switch(coupled_flag){
     689                case 0:{
     690                   IssmDouble  thickness,base,sealevel;
     691                        element->GetInputValue(&thickness, gauss,ThicknessEnum);
     692                        element->GetInputValue(&base, gauss,BaseEnum);
     693                        element->GetInputValue(&sealevel, gauss,SealevelEnum);
     694                        IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
     695                        IssmDouble rho_ice   = element->GetMaterialParameter(MaterialsRhoIceEnum);
     696                        IssmDouble gravity   = element->GetMaterialParameter(ConstantsGEnum);
     697                        Neff=gravity*(rho_ice*thickness+rho_water*(base-sealevel));
     698                                 }
     699                        break;
     700                case 1:
     701                        element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum);
     702                        break;
     703                case 2:
     704                        element->GetInputValue(&Neff,gauss,EffectivePressureTimeAverageEnum);
     705                        break;
     706                default:
     707                        _error_("not supported");
     708        }
     709
     710        /*Make sure Neff is positive*/
     711        if(Neff<0.) Neff=0.;
     712
     713        /*Return effective pressure*/
     714        return Neff;
     715
     716}/*}}}*/

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

    r21551 r22847  
    4343                void  GetAlpha2Weertman(IssmDouble* palpha2,Gauss* gauss);
    4444                void  GetAlpha2WeertmanTemp(IssmDouble* palpha2,Gauss* gauss);
     45
     46                IssmDouble EffectivePressure(Gauss* gauss);
    4547};
    4648
Note: See TracChangeset for help on using the changeset viewer.