Index: ../trunk-jpl/src/c/classes/Materials/Matestar.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Materials/Matestar.cpp (revision 21699) +++ ../trunk-jpl/src/c/classes/Materials/Matestar.cpp (revision 21700) @@ -159,6 +159,7 @@ } /*}}}*/ IssmDouble Matestar::GetB(){/*{{{*/ + /*Output*/ IssmDouble B; @@ -167,6 +168,7 @@ } /*}}}*/ IssmDouble Matestar::GetBbar(){/*{{{*/ + /*Output*/ IssmDouble Bbar; @@ -192,6 +194,15 @@ return Ec; } /*}}}*/ +IssmDouble Matestar::GetEcbar(){/*{{{*/ + + /*Output*/ + IssmDouble Ecbar; + + element->inputs->GetInputAverage(&Ecbar,MaterialsRheologyEcbarEnum); + return Ecbar; +} +/*}}}*/ IssmDouble Matestar::GetEs(){/*{{{*/ /*Output*/ @@ -201,6 +212,15 @@ return Es; } /*}}}*/ +IssmDouble Matestar::GetEsbar(){/*{{{*/ + + /*Output*/ + IssmDouble Esbar; + + element->inputs->GetInputAverage(&Esbar,MaterialsRheologyEsbarEnum); + return Esbar; +} +/*}}}*/ IssmDouble Matestar::GetN(){/*{{{*/ /*Output*/ @@ -210,51 +230,6 @@ /*}}}*/ void Matestar::GetViscosity(IssmDouble* pviscosity,IssmDouble eps_eff){/*{{{*/ _error_("not implemented yet"); - /*From a string tensor and a material object, return viscosity, using Glen's flow law. - (1-D) B - viscosity= ------------------------- - 2 eps_eff ^[(n-1)/n] - - where viscosity is the viscotiy, B the flow law parameter , eps_eff is the effective strain rate - and n the flow law exponent. - - If eps_eff = 0 , it means this is the first time SystemMatrices is being run, and we - return 10^14, initial viscosity. - */ - - /*output: */ - IssmDouble viscosity; - - /*Intermediary: */ - IssmDouble B,D=0.,n; - - /*Get B and n*/ - B=GetB(); _assert_(B>0.); - n=GetN(); _assert_(n>0.); - - if (n==1.){ - /*Linear Viscous behavior (Newtonian fluid) viscosity=B/2: */ - viscosity=(1.-D)*B/2.; - } - else{ - - /*if no strain rate, return maximum viscosity*/ - if(eps_eff==0.){ - viscosity = 1.e+14/2.; - //viscosity = B; - //viscosity=2.5*pow(10.,17); - } - - else{ - viscosity=(1.-D)*B/(2.*pow(eps_eff,(n-1.)/n)); - } - } - - /*Checks in debugging mode*/ - if(viscosity<=0) _error_("Negative viscosity"); - - /*Return: */ - *pviscosity=viscosity; } /*}}}*/ void Matestar::GetViscosityBar(IssmDouble* pviscosity,IssmDouble eps_eff){/*{{{*/ @@ -273,14 +248,16 @@ _error_("not implemented yet"); } /*}}}*/ -IssmDouble Matestar::GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz){/*{{{*/ +IssmDouble Matestar::GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,bool isdepthaveraged){/*{{{*/ + /*output: */ + IssmDouble viscosity; + /*Intermediaries*/ - IssmDouble viscosity; IssmDouble epseff,epsprime_norm; IssmDouble lambdas; IssmDouble vmag,dvmag[3]; - IssmDouble B,Ec,Es,E; + IssmDouble B,Ec,Es,E,n; /*Calculate velocity magnitude and its derivative*/ vmag = sqrt(vx*vx+vy*vy+vz*vz); @@ -299,9 +276,17 @@ lambdas=EstarLambdaS(epseff,epsprime_norm); /*Get B and enhancement*/ - B=GetB(); _assert_(B>0.); - Ec=GetEc(); _assert_(Ec>=0.); - Es=GetEs(); _assert_(Es>=0.); + n=GetN(); _assert_(n>0.); + if (isdepthaveraged==0.){ + B=GetB(); _assert_(B>0.); + Ec=GetEc(); _assert_(Ec>=0.); + Es=GetEs(); _assert_(Es>=0.); + } + else{ + B=GetBbar(); _assert_(B>0.); + Ec=GetEcbar(); _assert_(Ec>=0.); + Es=GetEsbar(); _assert_(Es>=0.); + } /*Get total enhancement factor E(lambdas)*/ E = Ec + (Es-Ec)*lambdas*lambdas; _assert_(E>0.); @@ -314,34 +299,64 @@ //viscosity=2.5*pow(10.,17); } else{ - viscosity = B/(2.*pow(E*epseff*epseff,1./3.)); + viscosity = B/(2.*pow(E,1./n)*pow(epseff,2./n)); } /*Assign output pointer*/ return viscosity; } /*}}}*/ -void Matestar::GetViscosity_B(IssmDouble* pdmudB,IssmDouble eps_eff){/*{{{*/ - /*output: */ +IssmDouble Matestar::GetViscosity_BGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,bool isdepthaveraged){/*{{{*/ + + /*Intermediaries*/ IssmDouble dmudB; + IssmDouble epseff,epsprime_norm; + IssmDouble lambdas; + IssmDouble vmag,dvmag[3]; + IssmDouble Ec,Es,E; - /*Intermediary: */ - IssmDouble E=1.,n; + /*Calculate velocity magnitude and its derivative*/ + vmag = sqrt(vx*vx+vy*vy+vz*vz); + if(vmag<1e-12){ + dvmag[0]=0; + dvmag[1]=0; + dvmag[2]=0; + } + else{ + dvmag[0]=1./(2*sqrt(vmag))*(2*vx*dvx[0]+2*vy*dvy[0]+2*vz*dvz[0]); + dvmag[1]=1./(2*sqrt(vmag))*(2*vx*dvx[1]+2*vy*dvy[1]+2*vz*dvz[1]); + dvmag[2]=1./(2*sqrt(vmag))*(2*vx*dvx[2]+2*vy*dvy[2]+2*vz*dvz[2]); + } - n=GetN(); _assert_(n>0.); - if(n==1.){ - /*Linear Viscous behavior (Newtonian fluid) dmudB=B/2E: */ - dmudB=1./(2.*E); + EstarStrainrateQuantities(&epseff,&epsprime_norm,vx,vy,vz,vmag,dvx,dvy,dvz,&dvmag[0]); + lambdas=EstarLambdaS(epseff,epsprime_norm); + + /*Get enhancement*/ + if (isdepthaveraged==0.){ + Ec=GetEc(); _assert_(Ec>=0.); + Es=GetEs(); _assert_(Es>=0.); } else{ - if(eps_eff==0.) dmudB = 0.; - else dmudB = 1./(2.*pow(E*eps_eff*eps_eff,1./3.)); + Ec=GetEcbar(); _assert_(Ec>=0.); + Es=GetEsbar(); _assert_(Es>=0.); } - /*Return: */ - *pdmudB=dmudB; + /*Get total enhancement factor E(lambdas)*/ + E = Ec + (Es-Ec)*lambdas*lambdas; _assert_(E>0.); + + /*Compute dmudB*/ + if(epseff==0.) dmudB = 0.; + else dmudB = 1./(2.*pow(E,1./3.)*pow(epseff,2./3.)); + + /*Assign output*/ + return dmudB; + } /*}}}*/ +void Matestar::GetViscosity_B(IssmDouble* pdmudB,IssmDouble eps_eff){/*{{{*/ + _error_("not implemented yet"); +} +/*}}}*/ void Matestar::GetViscosity_D(IssmDouble* pdmudD,IssmDouble eps_eff){/*{{{*/ _error_("not implemented yet"); } @@ -392,12 +407,99 @@ } /*}}}*/ +void Matestar::ViscosityBFS(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/ + + /*Intermediaries*/ + IssmDouble vx,vy,vz; + IssmDouble dvx[3],dvy[3],dvz[3]; + bool isdepthaveraged=0.; + + /*Get velocity derivatives in all directions*/ + _assert_(dim>1); + _assert_(vx_input); + vx_input->GetInputValue(&vx,gauss); + vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); + _assert_(vy_input); + vy_input->GetInputValue(&vy,gauss); + vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); + if(dim==3){ + _assert_(vz_input); + vz_input->GetInputValue(&vz,gauss); + vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss); + } + else{ + vz = 0.; + dvz[0] = 0.; dvz[1] = 0.; dvz[2] = 0.; + } + + /*Compute dmudB*/ + *pdmudB=GetViscosity_BGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged); +} +/*}}}*/ +void Matestar::ViscosityBHO(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ + + /*Intermediaries*/ + IssmDouble vx,vy,vz; + IssmDouble dvx[3],dvy[3],dvz[3]; + bool isdepthaveraged=0.; + + /*Get velocity derivatives in all directions*/ + _assert_(dim==2 || dim==3); + _assert_(vx_input); + vx_input->GetInputValue(&vx,gauss); + vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); + if(dim==3){ + _assert_(vy_input); + vy_input->GetInputValue(&vy,gauss); + vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); + } + else{ + dvx[2] = 0.; + vy = 0.; + dvy[0] = 0.; dvy[1] = 0.; dvy[2] = 0.; + } + vz = 0.; + dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1]; + + /*Compute viscosity*/ + *pdmudB=GetViscosity_BGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged); +}/*}}}*/ +void Matestar::ViscosityBSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ + /*Intermediaries*/ + IssmDouble vx,vy,vz; + IssmDouble dvx[3],dvy[3],dvz[3]; + bool isdepthaveraged=1.; + + /*Get velocity derivatives in all directions*/ + _assert_(dim==1 || dim==2); + _assert_(vx_input); + vx_input->GetInputValue(&vx,gauss); + vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); + if(dim==2){ + _assert_(vy_input); + vy_input->GetInputValue(&vy,gauss); + vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); + } + else{ + dvx[1] = 0.; + dvx[2] = 0.; + vy = 0.; + dvy[0] = 0.; dvy[1] = 0.; dvy[2] = 0.; + } + dvx[2] = 0.; + dvy[2] = 0.; + vz = 0.; + dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1]; + + /*Compute viscosity*/ + *pdmudB=GetViscosity_BGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged); +}/*}}}*/ void Matestar::ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/ /*Intermediaries*/ IssmDouble vx,vy,vz; IssmDouble dvx[3],dvy[3],dvz[3]; - IssmDouble B,Ec,Es; + bool isdepthaveraged=0.; /*Get velocity derivatives in all directions*/ _assert_(dim>1); @@ -418,7 +520,7 @@ } /*Compute viscosity*/ - *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]); + *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged); } /*}}}*/ void Matestar::ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/ @@ -429,6 +531,7 @@ /*Intermediaries*/ IssmDouble vx,vy,vz; IssmDouble dvx[3],dvy[3],dvz[3]; + bool isdepthaveraged=0.; /*Get velocity derivatives in all directions*/ _assert_(dim==2 || dim==3); @@ -449,7 +552,7 @@ dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1]; /*Compute viscosity*/ - *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]); + *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged); }/*}}}*/ void Matestar::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/ _error_("not implemented yet"); @@ -458,10 +561,11 @@ _error_("not implemented yet"); }/*}}}*/ void Matestar::ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ + /*Intermediaries*/ IssmDouble vx,vy,vz; IssmDouble dvx[3],dvy[3],dvz[3]; - IssmDouble B,Ec,Es; + bool isdepthaveraged=1.; /*Get velocity derivatives in all directions*/ _assert_(dim==1 || dim==2); @@ -485,7 +589,7 @@ dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1]; /*Compute viscosity*/ - *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]); + *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged); }/*}}}*/ void Matestar::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/ _error_("not implemented yet"); Index: ../trunk-jpl/src/c/classes/Materials/Material.h =================================================================== --- ../trunk-jpl/src/c/classes/Materials/Material.h (revision 21699) +++ ../trunk-jpl/src/c/classes/Materials/Material.h (revision 21700) @@ -52,6 +52,9 @@ virtual void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf)=0; virtual void ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0; virtual void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon)=0; + virtual void ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input)=0; + virtual void ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0; + virtual void ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0; }; #endif Index: ../trunk-jpl/src/c/classes/Materials/Matestar.h =================================================================== --- ../trunk-jpl/src/c/classes/Materials/Matestar.h (revision 21699) +++ ../trunk-jpl/src/c/classes/Materials/Matestar.h (revision 21700) @@ -69,7 +69,9 @@ IssmDouble GetD(); IssmDouble GetDbar(); IssmDouble GetEc(); + IssmDouble GetEcbar(); IssmDouble GetEs(); + IssmDouble GetEsbar(); IssmDouble GetN(); bool IsDamage(); bool IsEnhanced(){_error_("not supported");}; @@ -83,8 +85,12 @@ void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf); void ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon); + void ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); + void ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); + void ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); /*}}}*/ - IssmDouble GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz); + IssmDouble GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,bool isdepthaveraged); + IssmDouble GetViscosity_BGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,bool isdepthaveraged); }; #endif /* _MATESTAR_H_ */ Index: ../trunk-jpl/src/c/classes/Materials/Matice.h =================================================================== --- ../trunk-jpl/src/c/classes/Materials/Matice.h (revision 21699) +++ ../trunk-jpl/src/c/classes/Materials/Matice.h (revision 21700) @@ -87,6 +87,9 @@ void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf); void ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon); + void ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not supported");}; + void ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");}; + void ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");}; /*}}}*/ }; Index: ../trunk-jpl/src/c/classes/Materials/Matpar.h =================================================================== --- ../trunk-jpl/src/c/classes/Materials/Matpar.h (revision 21699) +++ ../trunk-jpl/src/c/classes/Materials/Matpar.h (revision 21700) @@ -124,6 +124,9 @@ void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf){_error_("not supported");}; void ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");}; void ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not supported");}; + void ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not supported");}; + void ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");}; + void ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");}; /*}}}*/ /*Numerics: {{{*/ void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure);