Changeset 17615

Timestamp:

04/01/14 09:03:21 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: added 1/sqrt(2) factor in effective straine rate for flow band models

File:

• issm/trunk-jpl/src/c/classes/Elements/Element.cpp (modified) (5 diffs)

issm/trunk-jpl/src/c/classes/Elements/Element.cpp

@@ -994 +994 @@
 /*}}}*/
 void Element::ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
+        /*The effective strain rate is defined in Paterson 3d Ed p 91 eq 9,
+         * and Cuffey p 303 eq 8.18:
+         *
+         *  2 eps_eff^2 = eps_xx^2 + eps_yy^2 + eps_zz^2 + 2(eps_xy^2 + eps_xz^2 + eps_yz^2)
+         *
+         *  or
+         *
+         *  eps_eff = 1/sqrt(2) sqrt( \sum_ij eps_ij^2 )
+         *
+         *          = 1/sqrt(2) ||eps||_F
+         *
+         *  where ||.||_F is the Frobenius norm */
 
         /*Intermediaries*/
@@ -1008 +1020 @@
         }
         else{
-                /* eps_eff^2 = exx^2 + eyy^2 + 2*exy^2 */
+                /* eps_eff^2 = 1/2 ( exx^2 + eyy^2 + 2*exy^2 )*/
                 this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
-                eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + 2.*epsilon2d[2]*epsilon2d[2]);
+                eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + 2.*epsilon2d[2]*epsilon2d[2]);
         }
 
@@ -1084 +1096 @@
         IssmDouble viscosity;
         IssmDouble epsilon3d[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/
-        IssmDouble epsilon2d[2]; /* epsilon=[exx,exy];           */
+        IssmDouble epsilon2d[2];/* epsilon=[exx,exy];            */
         IssmDouble eps_eff;
 
@@ -1093 +1105 @@
         }
         else{
-                /* eps_eff^2 = exx^2 + exy^2 */
+                /* eps_eff^2 = 1/2 (exx^2 + 2*exy^2 )*/
                 this->StrainRateHO2dvertical(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
-                eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1]);
+                eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + 2.*epsilon2d[1]*epsilon2d[1]);
         }
 
@@ -1113 +1125 @@
 
          if(dim==2){
-                 /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy+*/
+                 /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy*/
                  this->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*/
+                 /* eps_eff^2 = 1/2 exx^2*/
                  this->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input);
-                 eps_eff = sqrt(epsilon1d*epsilon1d);
+                 eps_eff = sqrt(epsilon1d*epsilon1d/2.);
          }