Index: ../trunk-jpl/src/c/shared/Elements/elements.h =================================================================== --- ../trunk-jpl/src/c/shared/Elements/elements.h (revision 21825) +++ ../trunk-jpl/src/c/shared/Elements/elements.h (revision 21826) @@ -16,7 +16,7 @@ // IssmDouble LliboutryDuval(IssmDouble temperature, IssmDouble waterfraction, IssmDouble depth,IssmDouble n); IssmDouble EstarLambdaS(IssmDouble epseff, IssmDouble epsprime_norm); void EstarOmega(IssmDouble* omega,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz, IssmDouble* dvmag); -void EstarStrainrateQuantities(IssmDouble *pepseff, IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble* dvmag); +void EstarStrainrateQuantities(IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble* dvmag); IssmDouble PddSurfaceMassBalance(IssmDouble* monthlytemperatures, IssmDouble* monthlyprec, IssmDouble* pdds, IssmDouble* pds, IssmDouble* melt, IssmDouble* accu, IssmDouble signorm, IssmDouble yts, IssmDouble h, IssmDouble s, IssmDouble desfac,IssmDouble s0t, Index: ../trunk-jpl/src/c/shared/Elements/EstarComponents.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Elements/EstarComponents.cpp (revision 21825) +++ ../trunk-jpl/src/c/shared/Elements/EstarComponents.cpp (revision 21826) @@ -2,12 +2,12 @@ #include "../Numerics/types.h" #include "../Exceptions/exceptions.h" #include "./elements.h" -void EstarStrainrateQuantities(IssmDouble *pepseff, IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble* dvmag){/*{{{*/ +void EstarStrainrateQuantities(IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble* dvmag){/*{{{*/ /*Intermediaries*/ IssmDouble omega[3]; IssmDouble nrsp[3],nrsp_norm; - IssmDouble eps[3][3],epseff; + IssmDouble eps[3][3]; IssmDouble epsprime[3],epsprime_norm; /*Get omega, correction for rigid body rotation*/ @@ -36,11 +36,6 @@ eps[1][0] = .5*(dvx[1]+dvy[0]); eps[1][1] = dvy[1]; eps[1][2] = .5*(dvy[2]+dvz[1]); eps[2][0] = .5*(dvx[2]+dvz[0]); eps[2][1] = .5*(dvy[2]+dvz[1]); eps[2][2] = dvz[2]; - /*effective strain rate*/ - epseff = 0.; - /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */ - epseff = sqrt(eps[0][0]*eps[0][0] + eps[1][1]*eps[1][1] + eps[0][1]*eps[0][1] + eps[0][2]*eps[0][2] + eps[1][2]*eps[1][2] + eps[0][0]*eps[1][1]); - /*Compute the shear strain rate on the non rotating shear plane*/ epsprime[0]=0.; epsprime[1]=0.; @@ -72,7 +67,6 @@ epsprime_norm = sqrt(epsprime[0]*epsprime[0] + epsprime[1]*epsprime[1] + epsprime[2]*epsprime[2]); /*Assign output pointers*/ - *pepseff=epseff; *pepsprime_norm=epsprime_norm; }/*}}}*/ void EstarOmega(IssmDouble* omega,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz, IssmDouble* dvmag){/*{{{*/ Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 21825) +++ ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 21826) @@ -47,7 +47,7 @@ /*Intermediaries*/ IssmDouble vx,vy,vz,vmag; IssmDouble dvx[3],dvy[3],dvz[3],dvmag[3]; - IssmDouble epseff,epsprime; + IssmDouble eps[3][3],epseff,epsprime; int dim; IssmDouble *xyz_list = NULL; @@ -100,7 +100,18 @@ 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]); } - EstarStrainrateQuantities(&epseff,&epsprime,vx,vy,vz,vmag,&dvx[0],&dvy[0],&dvz[0],&dvmag[0]); + + /*Build strain rate tensor*/ + eps[0][0] = dvx[0]; eps[0][1] = .5*(dvx[1]+dvy[0]); eps[0][2] = .5*(dvx[2]+dvz[0]); + eps[1][0] = .5*(dvx[1]+dvy[0]); eps[1][1] = dvy[1]; eps[1][2] = .5*(dvy[2]+dvz[1]); + eps[2][0] = .5*(dvx[2]+dvz[0]); eps[2][1] = .5*(dvy[2]+dvz[1]); eps[2][2] = dvz[2]; + + /*effective strain rate*/ + epseff = 0.; + /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */ + epseff = sqrt(eps[0][0]*eps[0][0] + eps[1][1]*eps[1][1] + eps[0][1]*eps[0][1] + eps[0][2]*eps[0][2] + eps[1][2]*eps[1][2] + eps[0][0]*eps[1][1]); + + EstarStrainrateQuantities(&epsprime,vx,vy,vz,vmag,&dvx[0],&dvy[0],&dvz[0],&dvmag[0]); lambdas[iv]=EstarLambdaS(epseff,epsprime); } @@ -676,7 +687,7 @@ material->GetViscosity_B(&dmudB,eps_eff); break; case MatestarEnum: - material->ViscosityBFS(&dmudB,dim,xyz_list,gauss,vx_input,vy_input,vz_input); + material->ViscosityBFS(&dmudB,dim,xyz_list,gauss,vx_input,vy_input,vz_input,eps_eff); break; default: _error_("not supported"); } @@ -713,7 +724,7 @@ material->GetViscosity_B(&dmudB,eps_eff); break; case MatestarEnum: - material->ViscosityBHO(&dmudB,dim,xyz_list,gauss,vx_input,vy_input); + material->ViscosityBHO(&dmudB,dim,xyz_list,gauss,vx_input,vy_input,eps_eff); break; default: _error_("not supported"); } @@ -750,7 +761,7 @@ material->GetViscosity_B(&dmudB,eps_eff); break; case MatestarEnum: - material->ViscosityBSSA(&dmudB,dim,xyz_list,gauss,vx_input,vy_input); + material->ViscosityBSSA(&dmudB,dim,xyz_list,gauss,vx_input,vy_input,eps_eff); break; default: _error_("not supported"); } Index: ../trunk-jpl/src/c/classes/Materials/Matice.h =================================================================== --- ../trunk-jpl/src/c/classes/Materials/Matice.h (revision 21825) +++ ../trunk-jpl/src/c/classes/Materials/Matice.h (revision 21826) @@ -87,9 +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");}; + void ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input,IssmDouble eps_eff){_error_("not supported");}; + void ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble eps_eff){_error_("not supported");}; + void ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble eps_eff){_error_("not supported");}; /*}}}*/ }; Index: ../trunk-jpl/src/c/classes/Materials/Matpar.h =================================================================== --- ../trunk-jpl/src/c/classes/Materials/Matpar.h (revision 21825) +++ ../trunk-jpl/src/c/classes/Materials/Matpar.h (revision 21826) @@ -124,9 +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");}; + void ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input,IssmDouble eps_eff){_error_("not supported");}; + void ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble eps_eff){_error_("not supported");}; + void ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble eps_eff){_error_("not supported");}; /*}}}*/ /*Numerics: {{{*/ void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure); Index: ../trunk-jpl/src/c/classes/Materials/Matestar.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Materials/Matestar.cpp (revision 21825) +++ ../trunk-jpl/src/c/classes/Materials/Matestar.cpp (revision 21826) @@ -248,13 +248,13 @@ _error_("not implemented yet"); } /*}}}*/ -IssmDouble Matestar::GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,bool isdepthaveraged){/*{{{*/ +IssmDouble Matestar::GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble eps_eff,bool isdepthaveraged){/*{{{*/ /*output: */ IssmDouble viscosity; /*Intermediaries*/ - IssmDouble epseff,epsprime_norm; + IssmDouble epsprime_norm; IssmDouble lambdas; IssmDouble vmag,dvmag[3]; IssmDouble B,Ec,Es,E,n; @@ -272,8 +272,8 @@ dvmag[2]=1./(2*sqrt(vmag))*(2*vx*dvx[2]+2*vy*dvy[2]+2*vz*dvz[2]); } - EstarStrainrateQuantities(&epseff,&epsprime_norm,vx,vy,vz,vmag,dvx,dvy,dvz,&dvmag[0]); - lambdas=EstarLambdaS(epseff,epsprime_norm); + EstarStrainrateQuantities(&epsprime_norm,vx,vy,vz,vmag,dvx,dvy,dvz,&dvmag[0]); + lambdas=EstarLambdaS(eps_eff,epsprime_norm); /*Get B and enhancement*/ n=GetN(); _assert_(n>0.); @@ -293,24 +293,27 @@ /*Compute viscosity*/ /*if no strain rate, return maximum viscosity*/ - if(epseff==0.){ + if(eps_eff==0.){ viscosity = 1.e+14/2.; //viscosity = B; //viscosity=2.5*pow(10.,17); } else{ - viscosity = B/(2.*pow(E,1./n)*pow(epseff,2./n)); + viscosity = B/(2.*pow(E,1./n)*pow(eps_eff,2./n)); } + /*Checks in debugging mode*/ + if(viscosity<=0) _error_("Negative viscosity"); + /*Assign output pointer*/ return viscosity; } /*}}}*/ -IssmDouble Matestar::GetViscosity_BGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,bool isdepthaveraged){/*{{{*/ +IssmDouble Matestar::GetViscosity_BGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble eps_eff,bool isdepthaveraged){/*{{{*/ /*Intermediaries*/ IssmDouble dmudB; - IssmDouble epseff,epsprime_norm; + IssmDouble epsprime_norm; IssmDouble lambdas; IssmDouble vmag,dvmag[3]; IssmDouble Ec,Es,E; @@ -328,8 +331,8 @@ dvmag[2]=1./(2*sqrt(vmag))*(2*vx*dvx[2]+2*vy*dvy[2]+2*vz*dvz[2]); } - EstarStrainrateQuantities(&epseff,&epsprime_norm,vx,vy,vz,vmag,dvx,dvy,dvz,&dvmag[0]); - lambdas=EstarLambdaS(epseff,epsprime_norm); + EstarStrainrateQuantities(&epsprime_norm,vx,vy,vz,vmag,dvx,dvy,dvz,&dvmag[0]); + lambdas=EstarLambdaS(eps_eff,epsprime_norm); /*Get enhancement*/ if (isdepthaveraged==0.){ @@ -345,8 +348,8 @@ 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.)); + if(eps_eff==0.) dmudB = 0.; + else dmudB = 1./(2.*pow(E,1./3.)*pow(eps_eff,2./3.)); /*Assign output*/ return dmudB; @@ -407,7 +410,7 @@ } /*}}}*/ -void Matestar::ViscosityBFS(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/ +void Matestar::ViscosityBFS(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input,IssmDouble eps_eff){/*{{{*/ /*Intermediaries*/ IssmDouble vx,vy,vz; @@ -433,10 +436,10 @@ } /*Compute dmudB*/ - *pdmudB=GetViscosity_BGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged); + *pdmudB=GetViscosity_BGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],eps_eff,isdepthaveraged); } /*}}}*/ -void Matestar::ViscosityBHO(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ +void Matestar::ViscosityBHO(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble eps_eff){/*{{{*/ /*Intermediaries*/ IssmDouble vx,vy,vz; @@ -462,9 +465,9 @@ 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); + *pdmudB=GetViscosity_BGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],eps_eff,isdepthaveraged); }/*}}}*/ -void Matestar::ViscosityBSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ +void Matestar::ViscosityBSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble eps_eff){/*{{{*/ /*Intermediaries*/ IssmDouble vx,vy,vz; IssmDouble dvx[3],dvy[3],dvz[3]; @@ -492,7 +495,7 @@ 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); + *pdmudB=GetViscosity_BGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],eps_eff,isdepthaveraged); }/*}}}*/ void Matestar::ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/ @@ -499,6 +502,9 @@ /*Intermediaries*/ IssmDouble vx,vy,vz; IssmDouble dvx[3],dvy[3],dvz[3]; + IssmDouble epsilon3d[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + IssmDouble epsilon2d[3]; /* epsilon=[exx,eyy,exy];*/ + IssmDouble eps_eff,eps0=1.e-27; bool isdepthaveraged=0.; /*Get velocity derivatives in all directions*/ @@ -519,8 +525,19 @@ dvz[0] = 0.; dvz[1] = 0.; dvz[2] = 0.; } + if(dim==3){ + /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */ + element->StrainRateFS(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input,vz_input); + eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[3]*epsilon3d[3] + epsilon3d[4]*epsilon3d[4] + epsilon3d[5]*epsilon3d[5] + epsilon3d[0]*epsilon3d[1]+eps0*eps0); + } + else{ + /* eps_eff^2 = 1/2 ( exx^2 + eyy^2 + 2*exy^2 )*/ + element->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input); + eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + 2.*epsilon2d[2]*epsilon2d[2]); + } + /*Compute viscosity*/ - *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged); + *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],eps_eff,isdepthaveraged); } /*}}}*/ void Matestar::ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/ @@ -531,8 +548,22 @@ /*Intermediaries*/ IssmDouble vx,vy,vz; IssmDouble dvx[3],dvy[3],dvz[3]; + IssmDouble epsilon3d[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ + IssmDouble epsilon2d[5]; /* epsilon=[exx,exy];*/ + IssmDouble eps_eff; bool isdepthaveraged=0.; + if(dim==3){ + /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */ + element->StrainRateHO(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input); + eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[2]*epsilon3d[2] + epsilon3d[3]*epsilon3d[3] + epsilon3d[4]*epsilon3d[4] + epsilon3d[0]*epsilon3d[1]); + } + else{ + /* eps_eff^2 = 1/2 (2*exx^2 + 2*exy^2 ) (since eps_zz = - eps_xx)*/ + element->StrainRateHO2dvertical(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input); + eps_eff = 1./sqrt(2.)*sqrt(2*epsilon2d[0]*epsilon2d[0] + 2*epsilon2d[1]*epsilon2d[1]); + } + /*Get velocity derivatives in all directions*/ _assert_(dim==2 || dim==3); _assert_(vx_input); @@ -552,7 +583,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],isdepthaveraged); + *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],eps_eff,isdepthaveraged); }/*}}}*/ void Matestar::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/ _error_("not implemented yet"); @@ -565,6 +596,9 @@ /*Intermediaries*/ IssmDouble vx,vy,vz; IssmDouble dvx[3],dvy[3],dvz[3]; + IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy]; */ + IssmDouble epsilon1d; /* epsilon=[exx]; */ + IssmDouble eps_eff; bool isdepthaveraged=1.; /*Get velocity derivatives in all directions*/ @@ -588,8 +622,19 @@ vz = 0.; dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1]; + if(dim==2){ + /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy*/ + element->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input); + eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + epsilon2d[2]*epsilon2d[2] + epsilon2d[0]*epsilon2d[1]); + } + else{ + /* eps_eff^2 = exx^2*/ + element->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input); + eps_eff = fabs(epsilon1d); + } + /*Compute viscosity*/ - *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged); + *pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],eps_eff,isdepthaveraged); }/*}}}*/ void Matestar::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/ _error_("not implemented yet"); Index: ../trunk-jpl/src/c/classes/Materials/Matice.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 21825) +++ ../trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 21826) @@ -723,7 +723,7 @@ IssmDouble viscosity; IssmDouble epsilon3d[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/ IssmDouble epsilon2d[2];/* epsilon=[exx,exy]; */ - IssmDouble eps_eff,E=1.0; + IssmDouble eps_eff; if(dim==3){ /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */ Index: ../trunk-jpl/src/c/classes/Materials/Material.h =================================================================== --- ../trunk-jpl/src/c/classes/Materials/Material.h (revision 21825) +++ ../trunk-jpl/src/c/classes/Materials/Material.h (revision 21826) @@ -52,9 +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; + virtual void ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input,IssmDouble epseff)=0; + virtual void ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble epseff)=0; + virtual void ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble epseff)=0; }; #endif Index: ../trunk-jpl/src/c/classes/Materials/Matestar.h =================================================================== --- ../trunk-jpl/src/c/classes/Materials/Matestar.h (revision 21825) +++ ../trunk-jpl/src/c/classes/Materials/Matestar.h (revision 21826) @@ -85,12 +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); + void ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input,IssmDouble eps_eff); + void ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble eps_eff); + void ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble eps_eff); /*}}}*/ - 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); + IssmDouble GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble eps_eff,bool isdepthaveraged); + IssmDouble GetViscosity_BGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble eps_eff,bool isdepthaveraged); }; #endif /* _MATESTAR_H_ */