Changeset 20945

Timestamp:

07/19/16 13:48:09 (9 years ago)

Author:

felicity

Message:

CHG: ESTAR vorticity correction revision

Location:

issm/trunk-jpl/src/c

Files:

• classes/Materials/Matestar.h (modified) (2 diffs)
• shared/Elements/EstarComponents.cpp (modified) (3 diffs)
• shared/FSanalyticals/fsanalyticals.cpp (modified) (3 diffs)

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

@@ -3 +3 @@
  */
 
-#ifndef MATEARL_H_
-#define MATEARL_H_
+#ifndef MATESTAR_H_
+#define MATESTAR_H_
 
 /*Headers:*/
@@ -85 +85 @@
 };
 
-#endif  /* _MATEARL_H_ */
+#endif  /* _MATESTAR_H_ */

issm/trunk-jpl/src/c/shared/Elements/EstarComponents.cpp

@@ -6 +6 @@
 
         /*Intermediaries*/
-        IssmDouble omega[3],omega_norm;
+        IssmDouble omega[3];
         IssmDouble nrsp[3],nrsp_norm;
         IssmDouble eps[3][3],epso;
@@ -81 +81 @@
         /*Intermediaries*/
         IssmDouble omega_norm;
+        IssmDouble omega_rigid[3];
 
         /*Create vorticity vector*/
@@ -87 +88 @@
 
         /*Create vorticity vector, corrected for rigid body rotation
-         * \overline{\omega} =\omega - \Omega
-         *                   =-\nabla\times{\bf v} - 2*\nabla U\times{\bf v}/U;
+         * \overline{\omega} =\omega - \omega_rigid
+         *                   =\nabla\times{\bf v} - 2*U*\kappa*\hat{b};
+         *                   =\nabla\times{\bf v} - (2*{\bf v}\times(({\bf v}\cdot\nabla)*{\bf v}))/U^2
          * check the magnitude of the second term -- if it is small, then the two
          * vorticities (omega and first term in omega) are approx. equal
          *
          * */
-        omega[0] = -(dvz[1] - dvy[2]) + 2*(dvmag[1]*vz - dvmag[2]*vy)/vmag;
-        omega[1] = -(dvx[2] - dvz[0]) + 2*(dvmag[2]*vx - dvmag[0]*vz)/vmag;
-        omega[2] = -(dvy[0] - dvx[1]) + 2*(dvmag[0]*vy - dvmag[1]*vx)/vmag;
+        omega_rigid[0] = 2*(vy*(vx*dvz[0]+vy*dvz[1]+vz*dvz[2]) - vz*(vx*dvy[0]+vy*dvy[1]+vz*dvy[2]))/(vmag*vmag);
+        omega_rigid[1] = 2*(vz*(vx*dvx[0]+vy*dvx[1]+vz*dvx[2]) - vx*(vx*dvz[0]+vy*dvz[1]+vz*dvz[2]))/(vmag*vmag);
+        omega_rigid[2] = 2*(vx*(vx*dvy[0]+vy*dvy[1]+vz*dvy[2]) - vy*(vx*dvx[0]+vy*dvx[1]+vz*dvx[2]))/(vmag*vmag);
+
+        omega[0] = (dvz[1] - dvy[2]) - omega_rigid[0];
+        omega[1] = (dvx[2] - dvz[0]) - omega_rigid[1];
+        omega[2] = (dvy[0] - dvx[1]) - omega_rigid[2];
 
         /*Take out vorticity component aligned with the velocity*/

issm/trunk-jpl/src/c/shared/FSanalyticals/fsanalyticals.cpp

@@ -9 +9 @@
  * case 6: 3d test with sinusoidal functions, non homogeneous Dirichlet conditions
  *
- * case 201: 3d test quadratic functions FS, EARL flow law
- * case 202: 3d test quadratic functions HO, EARL flow law
+ * case 201: 3d test quadratic functions FS, ESTAR flow law
+ * case 202: 3d test quadratic functions HO, ESTAR flow law
  */
 
@@ -83 +83 @@
                 case 212:
                         return (16.0L/9.0L)*x*(4*x + 1)/(pow(2*(4*pow(fabs((2*x + 1)*pow(y + 3*z, 3)/pow(fabs(y + 3*z), 3) - (2*x + 1)*(y + 3*z)/fabs(y + 3*z)), 2) + 9)/(8*pow(x, 2) + 2*pow(2*x + 1, 2) + 11) + 1, 1.0L/3.0L)*pow((4.0L/3.0L)*pow(x, 2) + (1.0L/3.0L)*pow(2*x + 1, 2) + 11.0L/6.0L, 4.0L/3.0L)) + (64.0L/3.0L)*x*(-(4*x + 1)*(4*pow(fabs((2*x + 1)*pow(y + 3*z, 3)/pow(fabs(y + 3*z), 3) - (2*x + 1)*(y + 3*z)/fabs(y + 3*z)), 2) + 9)/pow(8*pow(x, 2) + 2*pow(2*x + 1, 2) + 11, 2) + 2*(pow(y + 3*z, 3)/pow(fabs(y + 3*z), 3) - (y + 3*z)/fabs(y + 3*z))*((2*x + 1)*pow(y + 3*z, 3)/pow(fabs(y + 3*z), 3) - (2*x + 1)*(y + 3*z)/fabs(y + 3*z))/(8*pow(x, 2) + 2*pow(2*x + 1, 2) + 11))/(pow(2*(4*pow(fabs((2*x + 1)*pow(y + 3*z, 3)/pow(fabs(y + 3*z), 3) - (2*x + 1)*(y + 3*z)/fabs(y + 3*z)), 2) + 9)/(8*pow(x, 2) + 2*pow(2*x + 1, 2) + 11) + 1, 4.0L/3.0L)*pow((4.0L/3.0L)*pow(x, 2) + (1.0L/3.0L)*pow(2*x + 1, 2) + 11.0L/6.0L, 1.0L/3.0L)) + (4.0L/9.0L)*(4*x + 1)/(pow(2*(4*pow(fabs((2*x + 1)*pow(y + 3*z, 3)/pow(fabs(y + 3*z), 3) - (2*x + 1)*(y + 3*z)/fabs(y + 3*z)), 2) + 9)/(8*pow(x, 2) + 2*pow(2*x + 1, 2) + 11) + 1, 1.0L/3.0L)*pow((4.0L/3.0L)*pow(x, 2) + (1.0L/3.0L)*pow(2*x + 1, 2) + 11.0L/6.0L, 4.0L/3.0L)) - 4/(pow(2*(4*pow(fabs((2*x + 1)*pow(y + 3*z, 3)/pow(fabs(y + 3*z), 3) - (2*x + 1)*(y + 3*z)/fabs(y + 3*z)), 2) + 9)/(8*pow(x, 2) + 2*pow(2*x + 1, 2) + 11) + 1, 1.0L/3.0L)*pow((4.0L/3.0L)*pow(x, 2) + (1.0L/3.0L)*pow(2*x + 1, 2) + 11.0L/6.0L, 1.0L/3.0L)) + (16.0L/3.0L)*(-(4*x + 1)*(4*pow(fabs((2*x + 1)*pow(y + 3*z, 3)/pow(fabs(y + 3*z), 3) - (2*x + 1)*(y + 3*z)/fabs(y + 3*z)), 2) + 9)/pow(8*pow(x, 2) + 2*pow(2*x + 1, 2) + 11, 2) + 2*(pow(y + 3*z, 3)/pow(fabs(y + 3*z), 3) - (y + 3*z)/fabs(y + 3*z))*((2*x + 1)*pow(y + 3*z, 3)/pow(fabs(y + 3*z), 3) - (2*x + 1)*(y + 3*z)/fabs(y + 3*z))/(8*pow(x, 2) + 2*pow(2*x + 1, 2) + 11))/(pow(2*(4*pow(fabs((2*x + 1)*pow(y + 3*z, 3)/pow(fabs(y + 3*z), 3) - (2*x + 1)*(y + 3*z)/fabs(y + 3*z)), 2) + 9)/(8*pow(x, 2) + 2*pow(2*x + 1, 2) + 11) + 1, 4.0L/3.0L)*pow((4.0L/3.0L)*pow(x, 2) + (1.0L/3.0L)*pow(2*x + 1, 2) + 11.0L/6.0L, 1.0L/3.0L));
+
                 default:
                         _error_("FS analytical solution"<<testid<<" not implemented yet");
@@ -151 +152 @@
 
                 case 212:
-                                return -8.0L/3.0L*x*(9*(y + 3*z)*(pow(y + 3*z, 2)/pow(fabs(y + 3*z), 3) - 1/fabs(y + 3*z))/fabs(y + 3*z) + 4*((2*x + 1)*pow(y + 3*z, 3)/pow(fabs(y + 3*z), 3) - (2*x + 1)*(y + 3*z)/fabs(y + 3*z))*(3*(2*x + 1)*pow(y + 3*z, 4)/pow(fabs(y + 3*z), 5) - 4*(2*x + 1)*pow(y + 3*z, 2)/pow(fabs(y + 3*z), 3) + (2*x + 1)/fabs(y + 3*z)))/(pow(2*(4*pow(fabs((2*x + 1)*pow(y + 3*z, 3)/pow(fabs(y + 3*z), 3) - (2*x + 1)*(y + 3*z)/fabs(y + 3*z)), 2) + 9)/(8*pow(x, 2) + 2*pow(2*x + 1, 2) + 11) + 1, 4.0L/3.0L)*pow((4.0L/3.0L)*pow(x, 2) + (1.0L/3.0L)*pow(2*x + 1, 2) + 11.0L/6.0L, 1.0L/3.0L)*(8*pow(x, 2) + 2*pow(2*x + 1, 2) + 11)) - 26.0L/3.0L*(9*(y + 3*z)*(pow(y + 3*z, 2)/pow(fabs(y + 3*z), 3) - 1/fabs(y + 3*z))/fabs(y + 3*z) + 4*((2*x + 1)*pow(y + 3*z, 3)/pow(fabs(y + 3*z), 3) - (2*x + 1)*(y + 3*z)/fabs(y + 3*z))*(3*(2*x + 1)*pow(y + 3*z, 4)/pow(fabs(y + 3*z), 5) - 4*(2*x + 1)*pow(y + 3*z, 2)/pow(fabs(y + 3*z), 3) + (2*x + 1)/fabs(y + 3*z)))/(pow(2*(4*pow(fabs((2*x + 1)*pow(y + 3*z, 3)/pow(fabs(y + 3*z), 3) - (2*x + 1)*(y + 3*z)/fabs(y + 3*z)), 2) + 9)/(8*pow(x, 2) + 2*pow(2*x + 1, 2) + 11) + 1, 4.0L/3.0L)*pow((4.0L/3.0L)*pow(x, 2) + (1.0L/3.0L)*pow(2*x + 1, 2) + 11.0L/6.0L, 1.0L/3.0L)*(8*pow(x, 2) + 2*pow(2*x + 1, 2) + 11));
+                        return -8.0L/3.0L*x*(9*(y + 3*z)*(pow(y + 3*z, 2)/pow(fabs(y + 3*z), 3) - 1/fabs(y + 3*z))/fabs(y + 3*z) + 4*((2*x + 1)*pow(y + 3*z, 3)/pow(fabs(y + 3*z), 3) - (2*x + 1)*(y + 3*z)/fabs(y + 3*z))*(3*(2*x + 1)*pow(y + 3*z, 4)/pow(fabs(y + 3*z), 5) - 4*(2*x + 1)*pow(y + 3*z, 2)/pow(fabs(y + 3*z), 3) + (2*x + 1)/fabs(y + 3*z)))/(pow(2*(4*pow(fabs((2*x + 1)*pow(y + 3*z, 3)/pow(fabs(y + 3*z), 3) - (2*x + 1)*(y + 3*z)/fabs(y + 3*z)), 2) + 9)/(8*pow(x, 2) + 2*pow(2*x + 1, 2) + 11) + 1, 4.0L/3.0L)*pow((4.0L/3.0L)*pow(x, 2) + (1.0L/3.0L)*pow(2*x + 1, 2) + 11.0L/6.0L, 1.0L/3.0L)*(8*pow(x, 2) + 2*pow(2*x + 1, 2) + 11)) - 26.0L/3.0L*(9*(y + 3*z)*(pow(y + 3*z, 2)/pow(fabs(y + 3*z), 3) - 1/fabs(y + 3*z))/fabs(y + 3*z) + 4*((2*x + 1)*pow(y + 3*z, 3)/pow(fabs(y + 3*z), 3) - (2*x + 1)*(y + 3*z)/fabs(y + 3*z))*(3*(2*x + 1)*pow(y + 3*z, 4)/pow(fabs(y + 3*z), 5) - 4*(2*x + 1)*pow(y + 3*z, 2)/pow(fabs(y + 3*z), 3) + (2*x + 1)/fabs(y + 3*z)))/(pow(2*(4*pow(fabs((2*x + 1)*pow(y + 3*z, 3)/pow(fabs(y + 3*z), 3) - (2*x + 1)*(y + 3*z)/fabs(y + 3*z)), 2) + 9)/(8*pow(x, 2) + 2*pow(2*x + 1, 2) + 11) + 1, 4.0L/3.0L)*pow((4.0L/3.0L)*pow(x, 2) + (1.0L/3.0L)*pow(2*x + 1, 2) + 11.0L/6.0L, 1.0L/3.0L)*(8*pow(x, 2) + 2*pow(2*x + 1, 2) + 11));
 
                 default: