Ignore:
Timestamp:
05/19/17 14:48:02 (8 years ago)
Author:
Mathieu Morlighem
Message:

merged trunk-jpl and trunk for revision 21727

Location:
issm/trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk

  • issm/trunk/src

  • issm/trunk/src/c/classes/Materials/Matestar.cpp

    r21341 r21729  
    134134/*}}}*/
    135135IssmDouble Matestar::GetA(){/*{{{*/
    136         _error_("not implemented yet");
     136        /*
     137         * A = 1/B^n
     138    */
     139
     140        IssmDouble B,n;
     141
     142        element->inputs->GetInputAverage(&B,MaterialsRheologyBEnum);
     143        n=this->GetN();
     144
     145        return pow(B,-n);
    137146}
    138147/*}}}*/
    139148IssmDouble Matestar::GetAbar(){/*{{{*/
    140         _error_("not implemented yet");
     149        /*
     150         * A = 1/B^n
     151         */
     152
     153        IssmDouble B,n;
     154
     155        element->inputs->GetInputAverage(&B,MaterialsRheologyBbarEnum);
     156        n=this->GetN();
     157
     158        return pow(B,-n);
    141159}
    142160/*}}}*/
    143161IssmDouble Matestar::GetB(){/*{{{*/
    144         _error_("not implemented yet");
     162
     163        /*Output*/
     164        IssmDouble B;
     165
     166        element->inputs->GetInputAverage(&B,MaterialsRheologyBEnum);
     167        return B;
    145168}
    146169/*}}}*/
    147170IssmDouble Matestar::GetBbar(){/*{{{*/
    148171
    149         _error_("not implemented yet");
     172        /*Output*/
     173        IssmDouble Bbar;
     174
     175        element->inputs->GetInputAverage(&Bbar,MaterialsRheologyBbarEnum);
     176        return Bbar;
    150177}
    151178/*}}}*/
    … …  
    159186}
    160187/*}}}*/
     188IssmDouble Matestar::GetEc(){/*{{{*/
     189
     190        /*Output*/
     191        IssmDouble Ec;
     192
     193        element->inputs->GetInputAverage(&Ec,MaterialsRheologyEcEnum);
     194        return Ec;
     195}
     196/*}}}*/
     197IssmDouble Matestar::GetEcbar(){/*{{{*/
     198
     199        /*Output*/
     200        IssmDouble Ecbar;
     201
     202        element->inputs->GetInputAverage(&Ecbar,MaterialsRheologyEcbarEnum);
     203        return Ecbar;
     204}
     205/*}}}*/
     206IssmDouble Matestar::GetEs(){/*{{{*/
     207
     208        /*Output*/
     209        IssmDouble Es;
     210
     211        element->inputs->GetInputAverage(&Es,MaterialsRheologyEsEnum);
     212        return Es;
     213}
     214/*}}}*/
     215IssmDouble Matestar::GetEsbar(){/*{{{*/
     216
     217        /*Output*/
     218        IssmDouble Esbar;
     219
     220        element->inputs->GetInputAverage(&Esbar,MaterialsRheologyEsbarEnum);
     221        return Esbar;
     222}
     223/*}}}*/
    161224IssmDouble Matestar::GetN(){/*{{{*/
    162         _error_("not implemented yet");
     225
     226        /*Output*/
     227        IssmDouble n=3.0;
     228        return n;
    163229}
    164230/*}}}*/
    165231void  Matestar::GetViscosity(IssmDouble* pviscosity,IssmDouble eps_eff){/*{{{*/
    166232        _error_("not implemented yet");
    167         /*From a string tensor and a material object, return viscosity, using Glen's flow law.
    168                                                                 (1-D) B
    169           viscosity= -------------------------
    170                                                   2 eps_eff ^[(n-1)/n]
    171 
    172           where viscosity is the viscotiy, B the flow law parameter , eps_eff is the effective strain rate
    173           and n the flow law exponent.
    174 
    175           If eps_eff = 0 , it means this is the first time SystemMatrices is being run, and we
    176           return 10^14, initial viscosity.
    177           */
     233}
     234/*}}}*/
     235void  Matestar::GetViscosityBar(IssmDouble* pviscosity,IssmDouble eps_eff){/*{{{*/
     236        _error_("not implemented yet");
     237}
     238/*}}}*/
     239void  Matestar::GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){/*{{{*/
     240        _error_("not implemented yet");
     241}
     242/*}}}*/
     243void  Matestar::GetViscosityDComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){/*{{{*/
     244        _error_("not implemented yet");
     245}
     246/*}}}*/
     247void  Matestar::GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* epsilon){/*{{{*/
     248        _error_("not implemented yet");
     249}
     250/*}}}*/
     251IssmDouble Matestar::GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,bool isdepthaveraged){/*{{{*/
    178252
    179253        /*output: */
    180254        IssmDouble viscosity;
    181255
    182         /*Intermediary: */
    183         IssmDouble B,D=0.,n;
    184 
    185         /*Get B and n*/
    186         B=GetB(); _assert_(B>0.);
    187         n=GetN(); _assert_(n>0.);
    188 
    189         if (n==1.){
    190                 /*Linear Viscous behavior (Newtonian fluid) viscosity=B/2: */
    191                 viscosity=(1.-D)*B/2.;
    192         }
    193         else{
    194 
    195                 /*if no strain rate, return maximum viscosity*/
    196                 if(eps_eff==0.){
    197                         viscosity = 1.e+14/2.;
    198                         //viscosity = B;
    199                         //viscosity=2.5*pow(10.,17);
    200                 }
    201 
    202                 else{
    203                         viscosity=(1.-D)*B/(2.*pow(eps_eff,(n-1.)/n));
    204                 }
    205         }
    206 
    207         /*Checks in debugging mode*/
    208         if(viscosity<=0) _error_("Negative viscosity");
    209 
    210         /*Return: */
    211         *pviscosity=viscosity;
    212 }
    213 /*}}}*/
    214 void  Matestar::GetViscosityBar(IssmDouble* pviscosity,IssmDouble eps_eff){/*{{{*/
    215         _error_("not implemented yet");
    216 }
    217 /*}}}*/
    218 void  Matestar::GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){/*{{{*/
    219         _error_("not implemented yet");
    220 }
    221 /*}}}*/
    222 void  Matestar::GetViscosityDComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){/*{{{*/
    223         _error_("not implemented yet");
    224 }
    225 /*}}}*/
    226 void  Matestar::GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* epsilon){/*{{{*/
    227         _error_("not implemented yet");
    228 }
    229 /*}}}*/
    230 IssmDouble Matestar::GetViscosityGeneral(IssmDouble ko,IssmDouble Ec, IssmDouble Es,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz){/*{{{*/
    231 
    232         /*Intermediaries*/
    233         IssmDouble viscosity;
    234         IssmDouble E,lambdas;
    235         IssmDouble epso,epsprime_norm;
     256        /*Intermediaries*/
     257        IssmDouble epseff,epsprime_norm;
     258        IssmDouble lambdas;
    236259        IssmDouble vmag,dvmag[3];
     260        IssmDouble B,Ec,Es,E,n;
    237261
    238262        /*Calculate velocity magnitude and its derivative*/
    … …  
    249273        }
    250274
    251         EstarStrainrateQuantities(&epso,&epsprime_norm,vx,vy,vz,vmag,dvx,dvy,dvz,&dvmag[0]);
    252         lambdas=EstarLambdaS(epso,epsprime_norm);
     275        EstarStrainrateQuantities(&epseff,&epsprime_norm,vx,vy,vz,vmag,dvx,dvy,dvz,&dvmag[0]);
     276        lambdas=EstarLambdaS(epseff,epsprime_norm);
     277
     278        /*Get B and enhancement*/
     279        n=GetN(); _assert_(n>0.);
     280        if (isdepthaveraged==0.){
     281                B=GetB(); _assert_(B>0.);
     282                Ec=GetEc(); _assert_(Ec>=0.);
     283                Es=GetEs(); _assert_(Es>=0.);
     284        }
     285        else{
     286                B=GetBbar(); _assert_(B>0.);
     287                Ec=GetEcbar(); _assert_(Ec>=0.);
     288                Es=GetEsbar(); _assert_(Es>=0.);
     289        }
    253290
    254291        /*Get total enhancement factor E(lambdas)*/
    255         E = Ec + (Es-Ec)*lambdas*lambdas;
     292        E = Ec + (Es-Ec)*lambdas*lambdas; _assert_(E>0.);
    256293
    257294        /*Compute viscosity*/
    258         _assert_(E>0.);
    259         _assert_(ko>0.);
    260         _assert_(epso>0.);
    261         viscosity = 1./(2*pow(ko*E*epso*epso,1./3.));
     295        /*if no strain rate, return maximum viscosity*/
     296        if(epseff==0.){
     297                viscosity = 1.e+14/2.;
     298                //viscosity = B;
     299                //viscosity=2.5*pow(10.,17);
     300        }
     301        else{
     302                viscosity = B/(2.*pow(E,1./n)*pow(epseff,2./n));
     303        }
    262304
    263305        /*Assign output pointer*/
    … …  
    265307}
    266308/*}}}*/
     309IssmDouble Matestar::GetViscosity_BGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,bool isdepthaveraged){/*{{{*/
     310
     311        /*Intermediaries*/
     312        IssmDouble dmudB;
     313        IssmDouble epseff,epsprime_norm;
     314        IssmDouble lambdas;
     315        IssmDouble vmag,dvmag[3];
     316        IssmDouble Ec,Es,E;
     317
     318        /*Calculate velocity magnitude and its derivative*/
     319        vmag = sqrt(vx*vx+vy*vy+vz*vz);
     320        if(vmag<1e-12){
     321                dvmag[0]=0;
     322                dvmag[1]=0;
     323                dvmag[2]=0;
     324        }
     325        else{
     326                dvmag[0]=1./(2*sqrt(vmag))*(2*vx*dvx[0]+2*vy*dvy[0]+2*vz*dvz[0]);
     327                dvmag[1]=1./(2*sqrt(vmag))*(2*vx*dvx[1]+2*vy*dvy[1]+2*vz*dvz[1]);
     328                dvmag[2]=1./(2*sqrt(vmag))*(2*vx*dvx[2]+2*vy*dvy[2]+2*vz*dvz[2]);
     329        }
     330
     331        EstarStrainrateQuantities(&epseff,&epsprime_norm,vx,vy,vz,vmag,dvx,dvy,dvz,&dvmag[0]);
     332        lambdas=EstarLambdaS(epseff,epsprime_norm);
     333
     334        /*Get enhancement*/
     335        if (isdepthaveraged==0.){
     336                Ec=GetEc(); _assert_(Ec>=0.);
     337                Es=GetEs(); _assert_(Es>=0.);
     338        }
     339        else{
     340                Ec=GetEcbar(); _assert_(Ec>=0.);
     341                Es=GetEsbar(); _assert_(Es>=0.);
     342        }
     343
     344        /*Get total enhancement factor E(lambdas)*/
     345        E = Ec + (Es-Ec)*lambdas*lambdas; _assert_(E>0.);
     346
     347        /*Compute dmudB*/
     348        if(epseff==0.) dmudB = 0.;
     349        else           dmudB = 1./(2.*pow(E,1./3.)*pow(epseff,2./3.));
     350
     351        /*Assign output*/
     352        return dmudB;
     353
     354}
     355/*}}}*/
    267356void  Matestar::GetViscosity_B(IssmDouble* pdmudB,IssmDouble eps_eff){/*{{{*/
    268         _error_("not implemented yet");
     357         _error_("not implemented yet");
    269358}
    270359/*}}}*/
    … …  
    319408}
    320409/*}}}*/
    321 void  Matestar::ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
     410void  Matestar::ViscosityBFS(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
    322411
    323412        /*Intermediaries*/
    324413        IssmDouble vx,vy,vz;
    325414        IssmDouble dvx[3],dvy[3],dvz[3];
    326         IssmDouble ko,Ec,Es;
     415        bool isdepthaveraged=0.;
    327416
    328417        /*Get velocity derivatives in all directions*/
    … …  
    344433        }
    345434
    346         /*Get enhancement factors and ko*/
    347         Input* ec_input = element->inputs->GetInput(MaterialsRheologyEcEnum); _assert_(ec_input);
    348         Input* es_input = element->inputs->GetInput(MaterialsRheologyEsEnum); _assert_(es_input);
    349         Input* ko_input = element->inputs->GetInput(MaterialsRheologyKoEnum); _assert_(ko_input);
    350         ec_input->GetInputValue(&Ec,gauss);
    351         es_input->GetInputValue(&Es,gauss);
    352         ko_input->GetInputValue(&ko,gauss);
    353 
    354         /*Compute viscosity*/
    355         *pviscosity=GetViscosityGeneral(ko,Ec,Es,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
    356 }
    357 /*}}}*/
    358 void  Matestar::ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
    359         this->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon);
    360 }/*}}}*/
    361 void  Matestar::ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
     435        /*Compute dmudB*/
     436        *pdmudB=GetViscosity_BGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged);
     437}
     438/*}}}*/
     439void  Matestar::ViscosityBHO(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
    362440
    363441        /*Intermediaries*/
    364442        IssmDouble vx,vy,vz;
    365443        IssmDouble dvx[3],dvy[3],dvz[3];
    366         IssmDouble ko,Ec,Es;
     444        bool isdepthaveraged=0.;
    367445
    368446        /*Get velocity derivatives in all directions*/
    … …  
    384462        dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];
    385463
    386         /*Get enhancement factors and ko*/
    387         Input* ec_input = element->inputs->GetInput(MaterialsRheologyEcEnum); _assert_(ec_input);
    388         Input* es_input = element->inputs->GetInput(MaterialsRheologyEsEnum); _assert_(es_input);
    389         Input* ko_input = element->inputs->GetInput(MaterialsRheologyKoEnum); _assert_(ko_input);
    390         ec_input->GetInputValue(&Ec,gauss);
    391         es_input->GetInputValue(&Es,gauss);
    392         ko_input->GetInputValue(&ko,gauss);
    393 
    394464        /*Compute viscosity*/
    395         *pviscosity=GetViscosityGeneral(ko,Ec,Es,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
    396 }/*}}}*/
    397 void  Matestar::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
    398         _error_("not implemented yet");
    399 }/*}}}*/
    400 void  Matestar::ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surface_input){/*{{{*/
    401         _error_("not implemented yet");
    402 }/*}}}*/
    403 void  Matestar::ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
     465        *pdmudB=GetViscosity_BGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged);
     466}/*}}}*/
     467void  Matestar::ViscosityBSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
    404468        /*Intermediaries*/
    405469        IssmDouble vx,vy,vz;
    406470        IssmDouble dvx[3],dvy[3],dvz[3];
    407         IssmDouble ko,Ec,Es;
     471        bool isdepthaveraged=1.;
    408472
    409473        /*Get velocity derivatives in all directions*/
    … …  
    428492        dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];
    429493
    430         /*Get enhancement factors and ko*/
    431         Input* ec_input = element->inputs->GetInput(MaterialsRheologyEcbarEnum); _assert_(ec_input);
    432         Input* es_input = element->inputs->GetInput(MaterialsRheologyEsbarEnum); _assert_(es_input);
    433         Input* ko_input = element->inputs->GetInput(MaterialsRheologyKobarEnum); _assert_(ko_input);
    434         ec_input->GetInputValue(&Ec,gauss);
    435         es_input->GetInputValue(&Es,gauss);
    436         ko_input->GetInputValue(&ko,gauss);
    437 
    438494        /*Compute viscosity*/
    439         *pviscosity=GetViscosityGeneral(ko,Ec,Es,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
     495        *pdmudB=GetViscosity_BGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged);
     496}/*}}}*/
     497void  Matestar::ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
     498
     499        /*Intermediaries*/
     500        IssmDouble vx,vy,vz;
     501        IssmDouble dvx[3],dvy[3],dvz[3];
     502        bool isdepthaveraged=0.;
     503
     504        /*Get velocity derivatives in all directions*/
     505        _assert_(dim>1);
     506        _assert_(vx_input);
     507        vx_input->GetInputValue(&vx,gauss);
     508        vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
     509        _assert_(vy_input);
     510        vy_input->GetInputValue(&vy,gauss);
     511        vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
     512        if(dim==3){
     513                _assert_(vz_input);
     514                vz_input->GetInputValue(&vz,gauss);
     515                vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
     516        }
     517        else{
     518                vz = 0.;
     519                dvz[0] = 0.; dvz[1] = 0.; dvz[2] = 0.;
     520        }
     521
     522        /*Compute viscosity*/
     523        *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged);
     524}
     525/*}}}*/
     526void  Matestar::ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
     527        this->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon);
     528}/*}}}*/
     529void  Matestar::ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
     530
     531        /*Intermediaries*/
     532        IssmDouble vx,vy,vz;
     533        IssmDouble dvx[3],dvy[3],dvz[3];
     534        bool isdepthaveraged=0.;
     535
     536        /*Get velocity derivatives in all directions*/
     537        _assert_(dim==2 || dim==3);
     538        _assert_(vx_input);
     539        vx_input->GetInputValue(&vx,gauss);
     540        vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
     541        if(dim==3){
     542                _assert_(vy_input);
     543                vy_input->GetInputValue(&vy,gauss);
     544                vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
     545        }
     546        else{
     547                dvx[2] = 0.;
     548                vy = 0.;
     549                dvy[0] = 0.; dvy[1] = 0.; dvy[2] = 0.;
     550        }
     551        vz = 0.;
     552        dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];
     553
     554        /*Compute viscosity*/
     555        *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged);
     556}/*}}}*/
     557void  Matestar::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
     558        _error_("not implemented yet");
     559}/*}}}*/
     560void  Matestar::ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surface_input){/*{{{*/
     561        _error_("not implemented yet");
     562}/*}}}*/
     563void  Matestar::ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
     564
     565        /*Intermediaries*/
     566        IssmDouble vx,vy,vz;
     567        IssmDouble dvx[3],dvy[3],dvz[3];
     568        bool isdepthaveraged=1.;
     569
     570        /*Get velocity derivatives in all directions*/
     571        _assert_(dim==1 || dim==2);
     572        _assert_(vx_input);
     573        vx_input->GetInputValue(&vx,gauss);
     574        vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
     575        if(dim==2){
     576                _assert_(vy_input);
     577                vy_input->GetInputValue(&vy,gauss);
     578                vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
     579        }
     580        else{
     581                dvx[1] = 0.;
     582                dvx[2] = 0.;
     583                vy = 0.;
     584                dvy[0] = 0.; dvy[1] = 0.; dvy[2] = 0.;
     585        }
     586        dvx[2] = 0.;
     587        dvy[2] = 0.;
     588        vz = 0.;
     589        dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];
     590
     591        /*Compute viscosity*/
     592        *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged);
    440593}/*}}}*/
    441594void  Matestar::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
Note: See TracChangeset for help on using the changeset viewer.