Changeset 23839


Ignore:
Timestamp:
04/11/19 11:07:41 (6 years ago)
Author:
Mathieu Morlighem
Message:

CHG: added schoof law and centralized calculation of basal vel

Location:
issm/trunk-jpl/src
Files:
1 added
6 edited

Legend:

Unmodified
Added
Removed

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

    r23795 r23839  
    870870                        iomodel->FetchDataToInput(elements,"md.initialization.watercolumn",WatercolumnEnum,0.);
    871871                        break;
     872                case 11:
     873                        iomodel->FetchDataToInput(elements,"md.friction.C",FrictionCEnum);
     874                        iomodel->FetchDataToInput(elements,"md.friction.Cmax",FrictionCmaxEnum);
    872875                default:
    873876                        _error_("friction law "<< frictionlaw <<" not supported");
    … …  
    931934        int frictionlaw;
    932935        iomodel->FindConstant(&frictionlaw,"md.friction.law");
    933         if(frictionlaw==4 || frictionlaw==6) parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
    934         if(frictionlaw==3 || frictionlaw==4 || frictionlaw==1 || frictionlaw==7) parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
    935         if(frictionlaw==5) parameters->AddObject(iomodel->CopyConstantObject("md.friction.f",FrictionFEnum));
    936         if(frictionlaw==9) parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
    937         if(frictionlaw==10){
    938                 parameters->AddObject(iomodel->CopyConstantObject("md.friction.pseudoplasticity_exponent",FrictionPseudoplasticityExponentEnum));
    939                 parameters->AddObject(iomodel->CopyConstantObject("md.friction.threshold_speed",FrictionThresholdSpeedEnum));
    940                 parameters->AddObject(iomodel->CopyConstantObject("md.friction.delta",FrictionDeltaEnum));
    941                 parameters->AddObject(iomodel->CopyConstantObject("md.friction.void_ratio",FrictionVoidRatioEnum));
     936        switch(frictionlaw){
     937                case 1:
     938                        parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
     939                        break;
     940                case 2:
     941                        break;
     942                case 3:
     943                        parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
     944                        break;
     945                case 4:
     946                        parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
     947                        parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
     948                        break;
     949                case 5:
     950                        parameters->AddObject(iomodel->CopyConstantObject("md.friction.f",FrictionFEnum));
     951                        break;
     952                case 6:
     953                        parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
     954                        break;
     955                case 7:
     956                        parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
     957                        break;
     958                case 8:
     959                        break;
     960                case 9:
     961                        parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
     962                        break;
     963                case 10:
     964                        parameters->AddObject(iomodel->CopyConstantObject("md.friction.pseudoplasticity_exponent",FrictionPseudoplasticityExponentEnum));
     965                        parameters->AddObject(iomodel->CopyConstantObject("md.friction.threshold_speed",FrictionThresholdSpeedEnum));
     966                        parameters->AddObject(iomodel->CopyConstantObject("md.friction.delta",FrictionDeltaEnum));
     967                        parameters->AddObject(iomodel->CopyConstantObject("md.friction.void_ratio",FrictionVoidRatioEnum));
     968                        break;
     969                case 11:
     970                        parameters->AddObject(iomodel->CopyConstantObject("md.friction.m",FrictionMEnum));
     971                        break;
     972                default: _error_("Friction law "<<frictionlaw<<" not implemented yet");
    942973        }
    943974

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

    r23644 r23839  
    6666        IssmDouble  alpha;
    6767        IssmDouble  Chi,Gamma;
    68         IssmDouble  vx,vy,vz,vmag;
    6968        IssmDouble  alpha_complement;
    7069
    … …  
    7574        element->GetInputValue(&n,gauss,MaterialsRheologyNEnum);
    7675
    77         /*Get effective pressure*/
     76        /*Get effective pressure and velocity magnitude*/
    7877        IssmDouble Neff = EffectivePressure(gauss);
    79 
    80         //We need the velocity magnitude to evaluate the basal stress:
    81         switch(dim){
    82                 case 1:
    83                         element->GetInputValue(&vx,gauss,VxEnum);
    84                         vmag=sqrt(vx*vx);
    85                         break;
    86                 case 2:
    87                         element->GetInputValue(&vx,gauss,VxEnum);
    88                         element->GetInputValue(&vy,gauss,VyEnum);
    89                         vmag=sqrt(vx*vx+vy*vy);
    90                         break;
    91                 case 3:
    92                         element->GetInputValue(&vx,gauss,VxEnum);
    93                         element->GetInputValue(&vy,gauss,VyEnum);
    94                         element->GetInputValue(&vz,gauss,VzEnum);
    95                         vmag=sqrt(vx*vx+vy*vy+vz*vz);
    96                         break;
    97                 default:
    98                         _error_("not supported");
    99         }
    100         //      vmag=100./(3600.*24.*365.);
     78        IssmDouble vmag = VelMag(gauss);
    10179
    10280        if (q_exp==1){
    … …  
    155133        /*diverse: */
    156134        IssmDouble  r,s;
    157         IssmDouble  vx,vy,vz,vmag;
    158135        IssmDouble  drag_p,drag_q;
    159136        IssmDouble  drag_coefficient;
    … …  
    171148        /*Get effective pressure*/
    172149        IssmDouble Neff = EffectivePressure(gauss);
    173 
    174         //We need the velocity magnitude to evaluate the basal stress:
    175         switch(dim){
    176                 case 1:
    177                         element->GetInputValue(&vx,gauss,VxEnum);
    178                         vmag=sqrt(vx*vx);
    179                         break;
    180                 case 2:
    181                         element->GetInputValue(&vx,gauss,VxEnum);
    182                         element->GetInputValue(&vy,gauss,VyEnum);
    183                         vmag=sqrt(vx*vx+vy*vy);
    184                         break;
    185                 case 3:
    186                         element->GetInputValue(&vx,gauss,VxEnum);
    187                         element->GetInputValue(&vy,gauss,VyEnum);
    188                         element->GetInputValue(&vz,gauss,VzEnum);
    189                         vmag=sqrt(vx*vx+vy*vy+vz*vz);
    190                         break;
    191                 default:
    192                         _error_("not supported");
    193         }
     150        IssmDouble vmag = VelMag(gauss);
    194151
    195152        /*Check to prevent dividing by zero if vmag==0*/
    … …  
    237194                        GetAlpha2PISM(palpha2,gauss);
    238195                        break;
     196                case 11:
     197                        GetAlpha2Schoof(palpha2,gauss);
     198                        break;
    239199          default:
    240200                        _error_("Friction law "<< this->law <<" not supported");
    … …  
    251211        IssmDouble  drag_p, drag_q;
    252212        IssmDouble  thickness,base,floatation_thickness,sealevel;
    253         IssmDouble  vx,vy,vz,vmag;
    254213        IssmDouble  drag_coefficient,drag_coefficient_coulomb;
    255214        IssmDouble  alpha2,alpha2_coulomb;
    … …  
    273232        /*Get effective pressure*/
    274233        IssmDouble Neff = EffectivePressure(gauss);
    275 
    276         /*Compute velocity magnitude*/
    277         switch(dim){
    278                 case 1:
    279                         element->GetInputValue(&vx,gauss,VxEnum);
    280                         vmag=sqrt(vx*vx);
    281                         break;
    282                 case 2:
    283                         element->GetInputValue(&vx,gauss,VxEnum);
    284                         element->GetInputValue(&vy,gauss,VyEnum);
    285                         vmag=sqrt(vx*vx+vy*vy);
    286                         break;
    287                 case 3:
    288                         element->GetInputValue(&vx,gauss,VxEnum);
    289                         element->GetInputValue(&vy,gauss,VyEnum);
    290                         element->GetInputValue(&vz,gauss,VzEnum);
    291                         vmag=sqrt(vx*vx+vy*vy+vz*vz);
    292                         break;
    293                 default:
    294                         _error_("not supported");
    295         }
     234        IssmDouble vmag = VelMag(gauss);
    296235
    297236        /*Check to prevent dividing by zero if vmag==0*/
    … …  
    325264        IssmDouble  As;
    326265        IssmDouble  n;
    327 
    328266        IssmDouble  alpha;
    329267        IssmDouble  Chi,Gamma;
    330 
    331         IssmDouble  vx,vy,vz,vmag;
    332268        IssmDouble  alpha2;
    333269
    … …  
    340276        /*Get effective pressure*/
    341277        IssmDouble Neff = EffectivePressure(gauss);
    342 
    343         /*Compute velocity magnitude*/
    344         switch(dim){
    345                 case 1:
    346                         element->GetInputValue(&vx,gauss,VxEnum);
    347                         vmag=sqrt(vx*vx);
    348                         break;
    349                 case 2:
    350                         element->GetInputValue(&vx,gauss,VxEnum);
    351                         element->GetInputValue(&vy,gauss,VyEnum);
    352                         vmag=sqrt(vx*vx+vy*vy);
    353                         break;
    354                 case 3:
    355                         element->GetInputValue(&vx,gauss,VxEnum);
    356                         element->GetInputValue(&vy,gauss,VyEnum);
    357                         element->GetInputValue(&vz,gauss,VzEnum);
    358                         vmag=sqrt(vx*vx+vy*vy+vz*vz);
    359                         break;
    360                 default:
    361                         _error_("not supported");
    362         }
    363 
    364         //      vmag=100./(3600.*24.*365.);
     278        IssmDouble vmag = VelMag(gauss);
     279
    365280        //compute alpha and Chi coefficients: */
    366281        if (q_exp==1){
    … …  
    496411        IssmDouble  r,s;
    497412        IssmDouble  drag_p, drag_q;
    498         IssmDouble  vx,vy,vz,vmag;
    499413        IssmDouble  drag_coefficient;
    500414        IssmDouble  alpha2;
    … …  
    509423        s=1./drag_p;
    510424
    511         /*Get effective pressure*/
     425        /*Get effective pressure and basal velocity*/
    512426        IssmDouble Neff = EffectivePressure(gauss);
    513 
    514         /*Compute velocity magnitude*/
    515         switch(dim){
    516                 case 1:
    517                         element->GetInputValue(&vx,gauss,VxEnum);
    518                         vmag=sqrt(vx*vx);
    519                         break;
    520                 case 2:
    521                         element->GetInputValue(&vx,gauss,VxEnum);
    522                         element->GetInputValue(&vy,gauss,VyEnum);
    523                         vmag=sqrt(vx*vx+vy*vy);
    524                         break;
    525                 case 3:
    526                         element->GetInputValue(&vx,gauss,VxEnum);
    527                         element->GetInputValue(&vy,gauss,VyEnum);
    528                         element->GetInputValue(&vz,gauss,VzEnum);
    529                         vmag=sqrt(vx*vx+vy*vy+vz*vz);
    530                         break;
    531                 default:
    532                         _error_("not supported");
    533         }
     427        IssmDouble vmag = VelMag(gauss);
    534428
    535429        /*Check to prevent dividing by zero if vmag==0*/
    … …  
    551445        IssmDouble  Neff,F;
    552446        IssmDouble  thickness,base,sealevel;
    553         IssmDouble  vx,vy,vz,vmag;
    554447        IssmDouble  drag_coefficient,water_layer;
    555448        IssmDouble  alpha2;
    … …  
    579472        if(Neff<0) Neff=0;
    580473
    581         switch(dim){
    582                 case 1:
    583                         element->GetInputValue(&vx,gauss,VxEnum);
    584                         vmag=sqrt(vx*vx);
    585                         break;
    586                 case 2:
    587                         element->GetInputValue(&vx,gauss,VxEnum);
    588                         element->GetInputValue(&vy,gauss,VyEnum);
    589                         vmag=sqrt(vx*vx+vy*vy);
    590                         break;
    591                 case 3:
    592                         element->GetInputValue(&vx,gauss,VxEnum);
    593                         element->GetInputValue(&vy,gauss,VyEnum);
    594                         element->GetInputValue(&vz,gauss,VzEnum);
    595                         vmag=sqrt(vx*vx+vy*vy+vz*vz);
    596                         break;
    597                 default:
    598                         _error_("not supported");
    599         }
     474        IssmDouble vmag = VelMag(gauss);
    600475
    601476        /*Check to prevent dividing by zero if vmag==0*/
    … …  
    613488        /*diverse: */
    614489        IssmDouble  C,m;
    615         IssmDouble  vx,vy,vz,vmag;
    616490        IssmDouble  alpha2;
    617491
    … …  
    620494        element->GetInputValue(&m,FrictionMEnum);
    621495
    622         switch(dim){
    623                 case 1:
    624                         element->GetInputValue(&vx,gauss,VxEnum);
    625                         vmag=sqrt(vx*vx);
    626                         break;
    627                 case 2:
    628                         element->GetInputValue(&vx,gauss,VxEnum);
    629                         element->GetInputValue(&vy,gauss,VyEnum);
    630                         vmag=sqrt(vx*vx+vy*vy);
    631                         break;
    632                 case 3:
    633                         element->GetInputValue(&vx,gauss,VxEnum);
    634                         element->GetInputValue(&vy,gauss,VyEnum);
    635                         element->GetInputValue(&vz,gauss,VzEnum);
    636                         vmag=sqrt(vx*vx+vy*vy+vz*vz);
    637                         break;
    638                 default:
    639                         _error_("not supported");
    640         }
     496        /*Get velocity magnitude*/
     497        IssmDouble vmag = VelMag(gauss);
    641498
    642499        /*Check to prevent dividing by zero if vmag==0*/
    … …  
    741598        _assert_(!xIsNan<IssmDouble>(alpha2));
    742599        _assert_(!xIsInf<IssmDouble>(alpha2));
     600
     601        /*Assign output pointers:*/
     602        *palpha2=alpha2;
     603}/*}}}*/
     604void Friction::GetAlpha2Schoof(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
     605
     606        /*This routine calculates the basal friction coefficient
     607         *
     608         *               C |u_b|^(m-1)
     609         * alpha2= __________________________
     610         *          (1+(C/(Cmax N))^1/m |u_b| )^m
     611         *
     612         * */
     613
     614        /*diverse: */
     615        IssmDouble  C,Cmax,m;
     616
     617        /*Recover parameters: */
     618        element->GetInputValue(&Cmax,gauss,FrictionCmaxEnum);
     619        element->GetInputValue(&C,gauss,FrictionCEnum);
     620        element->GetInputValue(&m,FrictionMEnum);
     621
     622        /*Get effective pressure and velocity magnitude*/
     623        IssmDouble N  = EffectivePressure(gauss);
     624        IssmDouble ub = VelMag(gauss);
     625
     626        /*Compute alpha^2*/
     627        IssmDouble alpha2= (C*pow(ub,m-1.)) / pow(1.+  pow(C/(Cmax*N),1./m)*ub,m);
     628
     629        /*Checks*/
     630        _assert_(!xIsNan<IssmDouble>(alpha2));
     631        _assert_(alpha2>=0);
    743632
    744633        /*Assign output pointers:*/
    … …  
    810699
    811700}/*}}}*/
     701IssmDouble Friction::VelMag(Gauss* gauss){/*{{{*/
     702        /*Get effective pressure as a function of  flag */
     703
     704
     705        /*diverse*/
     706        IssmDouble vx,vy,vz,vmag;
     707
     708        switch(dim){
     709                case 1:
     710                        element->GetInputValue(&vx,gauss,VxEnum);
     711                        vmag=sqrt(vx*vx);
     712                        break;
     713                case 2:
     714                        element->GetInputValue(&vx,gauss,VxEnum);
     715                        element->GetInputValue(&vy,gauss,VyEnum);
     716                        vmag=sqrt(vx*vx+vy*vy);
     717                        break;
     718                case 3:
     719                        element->GetInputValue(&vx,gauss,VxEnum);
     720                        element->GetInputValue(&vy,gauss,VyEnum);
     721                        element->GetInputValue(&vz,gauss,VzEnum);
     722                        vmag=sqrt(vx*vx+vy*vy+vz*vz);
     723                        break;
     724                default:
     725                        _error_("not supported");
     726        }
     727
     728        return vmag;
     729
     730}/*}}}*/

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

    r23644 r23839  
    4141                void  GetAlpha2WeertmanTemp(IssmDouble* palpha2,Gauss* gauss);
    4242                void  GetAlpha2PISM(IssmDouble* palpha2,Gauss* gauss);
     43                void  GetAlpha2Schoof(IssmDouble* palpha2,Gauss* gauss);
    4344
    4445                IssmDouble EffectivePressure(Gauss* gauss);
     46                IssmDouble VelMag(Gauss* gauss);
    4547};
    4648

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

    r23814 r23839  
    502502        FrictionAsEnum,
    503503        FrictionCEnum,
     504        FrictionCmaxEnum,
    504505        FrictionCoefficientcoulombEnum,
    505506        FrictionCoefficientEnum,

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

    r23814 r23839  
    508508                case FrictionAsEnum : return "FrictionAs";
    509509                case FrictionCEnum : return "FrictionC";
     510                case FrictionCmaxEnum : return "FrictionCmax";
    510511                case FrictionCoefficientcoulombEnum : return "FrictionCoefficientcoulomb";
    511512                case FrictionCoefficientEnum : return "FrictionCoefficient";

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

    r23814 r23839  
    520520              else if (strcmp(name,"FrictionAs")==0) return FrictionAsEnum;
    521521              else if (strcmp(name,"FrictionC")==0) return FrictionCEnum;
     522              else if (strcmp(name,"FrictionCmax")==0) return FrictionCmaxEnum;
    522523              else if (strcmp(name,"FrictionCoefficientcoulomb")==0) return FrictionCoefficientcoulombEnum;
    523524              else if (strcmp(name,"FrictionCoefficient")==0) return FrictionCoefficientEnum;
    … …  
    628629              else if (strcmp(name,"SmbBNeg")==0) return SmbBNegEnum;
    629630              else if (strcmp(name,"SmbBPos")==0) return SmbBPosEnum;
    630               else if (strcmp(name,"SmbC")==0) return SmbCEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"SmbDailysnowfall")==0) return SmbDailysnowfallEnum;
     634              if (strcmp(name,"SmbC")==0) return SmbCEnum;
     635              else if (strcmp(name,"SmbDailysnowfall")==0) return SmbDailysnowfallEnum;
    635636              else if (strcmp(name,"SmbDailyrainfall")==0) return SmbDailyrainfallEnum;
    636637              else if (strcmp(name,"SmbDailydsradiation")==0) return SmbDailydsradiationEnum;
    … …  
    751752              else if (strcmp(name,"VzMesh")==0) return VzMeshEnum;
    752753              else if (strcmp(name,"VzSSA")==0) return VzSSAEnum;
    753               else if (strcmp(name,"Watercolumn")==0) return WatercolumnEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"WaterColumnOld")==0) return WaterColumnOldEnum;
     757              if (strcmp(name,"Watercolumn")==0) return WatercolumnEnum;
     758              else if (strcmp(name,"WaterColumnOld")==0) return WaterColumnOldEnum;
    758759              else if (strcmp(name,"WaterfractionDrainage")==0) return WaterfractionDrainageEnum;
    759760              else if (strcmp(name,"WaterfractionDrainageIntegrated")==0) return WaterfractionDrainageIntegratedEnum;
    … …  
    874875              else if (strcmp(name,"FSvelocity")==0) return FSvelocityEnum;
    875876              else if (strcmp(name,"FullMeltOnPartiallyFloating")==0) return FullMeltOnPartiallyFloatingEnum;
    876               else if (strcmp(name,"GaussPenta")==0) return GaussPentaEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"GaussSeg")==0) return GaussSegEnum;
     880              if (strcmp(name,"GaussPenta")==0) return GaussPentaEnum;
     881              else if (strcmp(name,"GaussSeg")==0) return GaussSegEnum;
    881882              else if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum;
    882883              else if (strcmp(name,"GaussTria")==0) return GaussTriaEnum;
    … …  
    997998              else if (strcmp(name,"NoFrictionOnPartiallyFloating")==0) return NoFrictionOnPartiallyFloatingEnum;
    998999              else if (strcmp(name,"NoMeltOnPartiallyFloating")==0) return NoMeltOnPartiallyFloatingEnum;
    999               else if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"None")==0) return NoneEnum;
     1003              if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum;
     1004              else if (strcmp(name,"None")==0) return NoneEnum;
    10041005              else if (strcmp(name,"Numberedcostfunction")==0) return NumberedcostfunctionEnum;
    10051006              else if (strcmp(name,"Numericalflux")==0) return NumericalfluxEnum;
    … …  
    11201121              else if (strcmp(name,"P2bubble")==0) return P2bubbleEnum;
    11211122              else if (strcmp(name,"P2")==0) return P2Enum;
    1122               else if (strcmp(name,"P2xP1")==0) return P2xP1Enum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"P2xP4")==0) return P2xP4Enum;
     1126              if (strcmp(name,"P2xP1")==0) return P2xP1Enum;
     1127              else if (strcmp(name,"P2xP4")==0) return P2xP4Enum;
    11271128              else if (strcmp(name,"Paterson")==0) return PatersonEnum;
    11281129              else if (strcmp(name,"Pengrid")==0) return PengridEnum;
    … …  
    12431244              else if (strcmp(name,"MeshZ")==0) return MeshZEnum;
    12441245              else if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum;
    1245               else if (strcmp(name,"OutputdefinitionList")==0) return OutputdefinitionListEnum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"SealevelObs")==0) return SealevelObsEnum;
     1249              if (strcmp(name,"OutputdefinitionList")==0) return OutputdefinitionListEnum;
     1250              else if (strcmp(name,"SealevelObs")==0) return SealevelObsEnum;
    12501251              else if (strcmp(name,"SealevelWeights")==0) return SealevelWeightsEnum;
    12511252              else if (strcmp(name,"StrainRate")==0) return StrainRateEnum;
Note: See TracChangeset for help on using the changeset viewer.