source: issm/oecreview/Archive/21337-21723/ISSM-21699-21700.diff@ 21726

Last change on this file since 21726 was 21726, checked in by Mathieu Morlighem, 8 years ago

CHG added Archive/21337-21723
File size: 15.4 KB

  • ../trunk-jpl/src/c/classes/Materials/Matestar.cpp

     
    159159}
    160160/*}}}*/
    161161IssmDouble Matestar::GetB(){/*{{{*/
     162
    162163        /*Output*/
    163164        IssmDouble B;
    164165
    … …  
    167168}
    168169/*}}}*/
    169170IssmDouble Matestar::GetBbar(){/*{{{*/
     171
    170172        /*Output*/
    171173        IssmDouble Bbar;
    172174
    … …  
    192194        return Ec;
    193195}
    194196/*}}}*/
     197IssmDouble Matestar::GetEcbar(){/*{{{*/
     198
     199        /*Output*/
     200        IssmDouble Ecbar;
     201
     202        element->inputs->GetInputAverage(&Ecbar,MaterialsRheologyEcbarEnum);
     203        return Ecbar;
     204}
     205/*}}}*/
    195206IssmDouble Matestar::GetEs(){/*{{{*/
    196207
    197208        /*Output*/
    … …  
    201212        return Es;
    202213}
    203214/*}}}*/
     215IssmDouble Matestar::GetEsbar(){/*{{{*/
     216
     217        /*Output*/
     218        IssmDouble Esbar;
     219
     220        element->inputs->GetInputAverage(&Esbar,MaterialsRheologyEsbarEnum);
     221        return Esbar;
     222}
     223/*}}}*/
    204224IssmDouble Matestar::GetN(){/*{{{*/
    205225
    206226        /*Output*/
    … …  
    210230/*}}}*/
    211231void  Matestar::GetViscosity(IssmDouble* pviscosity,IssmDouble eps_eff){/*{{{*/
    212232        _error_("not implemented yet");
    213         /*From a string tensor and a material object, return viscosity, using Glen's flow law.
    214                                                                 (1-D) B
    215           viscosity= -------------------------
    216                                                   2 eps_eff ^[(n-1)/n]
    217 
    218           where viscosity is the viscotiy, B the flow law parameter , eps_eff is the effective strain rate
    219           and n the flow law exponent.
    220 
    221           If eps_eff = 0 , it means this is the first time SystemMatrices is being run, and we
    222           return 10^14, initial viscosity.
    223           */
    224 
    225         /*output: */
    226         IssmDouble viscosity;
    227 
    228         /*Intermediary: */
    229         IssmDouble B,D=0.,n;
    230 
    231         /*Get B and n*/
    232         B=GetB(); _assert_(B>0.);
    233         n=GetN(); _assert_(n>0.);
    234 
    235         if (n==1.){
    236                 /*Linear Viscous behavior (Newtonian fluid) viscosity=B/2: */
    237                 viscosity=(1.-D)*B/2.;
    238         }
    239         else{
    240 
    241                 /*if no strain rate, return maximum viscosity*/
    242                 if(eps_eff==0.){
    243                         viscosity = 1.e+14/2.;
    244                         //viscosity = B;
    245                         //viscosity=2.5*pow(10.,17);
    246                 }
    247 
    248                 else{
    249                         viscosity=(1.-D)*B/(2.*pow(eps_eff,(n-1.)/n));
    250                 }
    251         }
    252 
    253         /*Checks in debugging mode*/
    254         if(viscosity<=0) _error_("Negative viscosity");
    255 
    256         /*Return: */
    257         *pviscosity=viscosity;
    258233}
    259234/*}}}*/
    260235void  Matestar::GetViscosityBar(IssmDouble* pviscosity,IssmDouble eps_eff){/*{{{*/
    … …  
    273248        _error_("not implemented yet");
    274249}
    275250/*}}}*/
    276 IssmDouble Matestar::GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz){/*{{{*/
     251IssmDouble Matestar::GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,bool isdepthaveraged){/*{{{*/
    277252
     253        /*output: */
     254        IssmDouble viscosity;
     255
    278256        /*Intermediaries*/
    279         IssmDouble viscosity;
    280257        IssmDouble epseff,epsprime_norm;
    281258        IssmDouble lambdas;
    282259        IssmDouble vmag,dvmag[3];
    283         IssmDouble B,Ec,Es,E;
     260        IssmDouble B,Ec,Es,E,n;
    284261
    285262        /*Calculate velocity magnitude and its derivative*/
    286263        vmag = sqrt(vx*vx+vy*vy+vz*vz);
    … …  
    299276        lambdas=EstarLambdaS(epseff,epsprime_norm);
    300277
    301278        /*Get B and enhancement*/
    302         B=GetB(); _assert_(B>0.);
    303         Ec=GetEc(); _assert_(Ec>=0.);
    304         Es=GetEs(); _assert_(Es>=0.);
     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        }
    305290
    306291        /*Get total enhancement factor E(lambdas)*/
    307292        E = Ec + (Es-Ec)*lambdas*lambdas; _assert_(E>0.);
    … …  
    314299                //viscosity=2.5*pow(10.,17);
    315300        }
    316301        else{
    317                 viscosity = B/(2.*pow(E*epseff*epseff,1./3.));
     302                viscosity = B/(2.*pow(E,1./n)*pow(epseff,2./n));
    318303        }
    319304
    320305        /*Assign output pointer*/
    321306        return viscosity;
    322307}
    323308/*}}}*/
    324 void  Matestar::GetViscosity_B(IssmDouble* pdmudB,IssmDouble eps_eff){/*{{{*/
    325         /*output: */
     309IssmDouble Matestar::GetViscosity_BGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,bool isdepthaveraged){/*{{{*/
     310
     311        /*Intermediaries*/
    326312        IssmDouble dmudB;
     313        IssmDouble epseff,epsprime_norm;
     314        IssmDouble lambdas;
     315        IssmDouble vmag,dvmag[3];
     316        IssmDouble Ec,Es,E;
    327317
    328         /*Intermediary: */
    329         IssmDouble E=1.,n;
     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        }
    330330
    331         n=GetN(); _assert_(n>0.);
    332         if(n==1.){
    333                 /*Linear Viscous behavior (Newtonian fluid) dmudB=B/2E: */
    334                 dmudB=1./(2.*E);
     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.);
    335338        }
    336339        else{
    337                 if(eps_eff==0.) dmudB = 0.;
    338                 else            dmudB = 1./(2.*pow(E*eps_eff*eps_eff,1./3.));
     340                Ec=GetEcbar(); _assert_(Ec>=0.);
     341                Es=GetEsbar(); _assert_(Es>=0.);
    339342        }
    340343
    341         /*Return: */
    342         *pdmudB=dmudB;
     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
    343354}
    344355/*}}}*/
     356void  Matestar::GetViscosity_B(IssmDouble* pdmudB,IssmDouble eps_eff){/*{{{*/
     357         _error_("not implemented yet");
     358}
     359/*}}}*/
    345360void  Matestar::GetViscosity_D(IssmDouble* pdmudD,IssmDouble eps_eff){/*{{{*/
    346361         _error_("not implemented yet");
    347362}
    … …  
    392407
    393408}
    394409/*}}}*/
     410void  Matestar::ViscosityBFS(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
     411
     412        /*Intermediaries*/
     413        IssmDouble vx,vy,vz;
     414        IssmDouble dvx[3],dvy[3],dvz[3];
     415        bool isdepthaveraged=0.;
     416
     417        /*Get velocity derivatives in all directions*/
     418        _assert_(dim>1);
     419        _assert_(vx_input);
     420        vx_input->GetInputValue(&vx,gauss);
     421        vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
     422        _assert_(vy_input);
     423        vy_input->GetInputValue(&vy,gauss);
     424        vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
     425        if(dim==3){
     426                _assert_(vz_input);
     427                vz_input->GetInputValue(&vz,gauss);
     428                vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
     429        }
     430        else{
     431                vz = 0.;
     432                dvz[0] = 0.; dvz[1] = 0.; dvz[2] = 0.;
     433        }
     434
     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){/*{{{*/
     440
     441        /*Intermediaries*/
     442        IssmDouble vx,vy,vz;
     443        IssmDouble dvx[3],dvy[3],dvz[3];
     444        bool isdepthaveraged=0.;
     445
     446        /*Get velocity derivatives in all directions*/
     447        _assert_(dim==2 || dim==3);
     448        _assert_(vx_input);
     449        vx_input->GetInputValue(&vx,gauss);
     450        vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
     451        if(dim==3){
     452                _assert_(vy_input);
     453                vy_input->GetInputValue(&vy,gauss);
     454                vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
     455        }
     456        else{
     457                dvx[2] = 0.;
     458                vy = 0.;
     459                dvy[0] = 0.; dvy[1] = 0.; dvy[2] = 0.;
     460        }
     461        vz = 0.;
     462        dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];
     463
     464        /*Compute viscosity*/
     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){/*{{{*/
     468        /*Intermediaries*/
     469        IssmDouble vx,vy,vz;
     470        IssmDouble dvx[3],dvy[3],dvz[3];
     471        bool isdepthaveraged=1.;
     472
     473        /*Get velocity derivatives in all directions*/
     474        _assert_(dim==1 || dim==2);
     475        _assert_(vx_input);
     476        vx_input->GetInputValue(&vx,gauss);
     477        vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
     478        if(dim==2){
     479                _assert_(vy_input);
     480                vy_input->GetInputValue(&vy,gauss);
     481                vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
     482        }
     483        else{
     484                dvx[1] = 0.;
     485                dvx[2] = 0.;
     486                vy = 0.;
     487                dvy[0] = 0.; dvy[1] = 0.; dvy[2] = 0.;
     488        }
     489        dvx[2] = 0.;
     490        dvy[2] = 0.;
     491        vz = 0.;
     492        dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];
     493
     494        /*Compute viscosity*/
     495        *pdmudB=GetViscosity_BGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged);
     496}/*}}}*/
    395497void  Matestar::ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
    396498
    397499        /*Intermediaries*/
    398500        IssmDouble vx,vy,vz;
    399501        IssmDouble dvx[3],dvy[3],dvz[3];
    400         IssmDouble B,Ec,Es;
     502        bool isdepthaveraged=0.;
    401503
    402504        /*Get velocity derivatives in all directions*/
    403505        _assert_(dim>1);
    … …  
    418520        }
    419521
    420522        /*Compute viscosity*/
    421         *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
     523        *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged);
    422524}
    423525/*}}}*/
    424526void  Matestar::ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
    … …  
    429531        /*Intermediaries*/
    430532        IssmDouble vx,vy,vz;
    431533        IssmDouble dvx[3],dvy[3],dvz[3];
     534        bool isdepthaveraged=0.;
    432535
    433536        /*Get velocity derivatives in all directions*/
    434537        _assert_(dim==2 || dim==3);
    … …  
    449552        dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];
    450553
    451554        /*Compute viscosity*/
    452         *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
     555        *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged);
    453556}/*}}}*/
    454557void  Matestar::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
    455558        _error_("not implemented yet");
    … …  
    458561        _error_("not implemented yet");
    459562}/*}}}*/
    460563void  Matestar::ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
     564
    461565        /*Intermediaries*/
    462566        IssmDouble vx,vy,vz;
    463567        IssmDouble dvx[3],dvy[3],dvz[3];
    464         IssmDouble B,Ec,Es;
     568        bool isdepthaveraged=1.;
    465569
    466570        /*Get velocity derivatives in all directions*/
    467571        _assert_(dim==1 || dim==2);
    … …  
    485589        dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];
    486590
    487591        /*Compute viscosity*/
    488         *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
     592        *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged);
    489593}/*}}}*/
    490594void  Matestar::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
    491595        _error_("not implemented yet");

  • ../trunk-jpl/src/c/classes/Materials/Material.h

     
    5252                virtual void       ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf)=0;
    5353                virtual void       ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0;
    5454                virtual void       ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon)=0;
     55                virtual void       ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input)=0;
     56                virtual void       ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0;
     57                virtual void       ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0;
    5558
    5659};
    5760#endif

  • ../trunk-jpl/src/c/classes/Materials/Matestar.h

     
    6969                IssmDouble GetD();
    7070                IssmDouble GetDbar();
    7171                IssmDouble GetEc();
     72                IssmDouble GetEcbar();
    7273                IssmDouble GetEs();
     74                IssmDouble GetEsbar();
    7375                IssmDouble GetN();
    7476                bool       IsDamage();
    7577                bool       IsEnhanced(){_error_("not supported");};
    … …  
    8385                void       ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf);
    8486                void       ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
    8587                void       ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
     88                void       ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
     89                void       ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
     90                void       ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
    8691                /*}}}*/
    87                 IssmDouble GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz);
     92                IssmDouble GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,bool isdepthaveraged);
     93                IssmDouble GetViscosity_BGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,bool isdepthaveraged);
    8894};
    8995
    9096#endif  /* _MATESTAR_H_ */

  • ../trunk-jpl/src/c/classes/Materials/Matice.h

     
    8787                void       ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf);
    8888                void       ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
    8989                void       ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
     90                void       ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not supported");};
     91                void       ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");};
     92                void       ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");};
    9093                /*}}}*/
    9194};
    9295

  • ../trunk-jpl/src/c/classes/Materials/Matpar.h

     
    124124                void       ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf){_error_("not supported");};
    125125                void       ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");};
    126126                void       ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not supported");};
     127                void       ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not supported");};
     128                void       ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");};
     129                void       ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");};
    127130                /*}}}*/
    128131                /*Numerics: {{{*/
    129132                void       EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure);
Note: See TracBrowser for help on using the repository browser.