Move rho_dh into the density sub-structure. Update the pressure and...

Move rho_dh into the density sub-structure. Update the pressure and sound-speed after kick and drift

src/const.h

@@ -40,8 +40,8 @@
 #define const_isothermal_internal_energy 20.2615290634f
 
 /* Dimensionality of the problem */
-#define HYDRO_DIMENSION_3D
-//#define HYDRO_DIMENSION_2D
+//#define HYDRO_DIMENSION_3D
+#define HYDRO_DIMENSION_2D
 //#define HYDRO_DIMENSION_1D
 
 /* Hydrodynamical adiabatic index. */
@@ -63,8 +63,8 @@
 //#define WENDLAND_C6_KERNEL
 
 /* SPH variant to use */
-//#define MINIMAL_SPH
-#define GADGET2_SPH
+#define MINIMAL_SPH
+//#define GADGET2_SPH
 //#define HOPKINS_PE_SPH
 //#define DEFAULT_SPH
 //#define GIZMO_SPH

src/hydro/Gadget2/hydro.h

@@ -195,7 +195,7 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
   p->density.wcount = 0.f;
   p->density.wcount_dh = 0.f;
   p->rho = 0.f;
-  p->rho_dh = 0.f;
+  p->density.rho_dh = 0.f;
   p->density.div_v = 0.f;
   p->density.rot_v[0] = 0.f;
   p->density.rot_v[1] = 0.f;
@@ -222,20 +222,17 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
 
   /* Final operation on the density (add self-contribution). */
   p->rho += p->mass * kernel_root;
-  p->rho_dh -= hydro_dimension * p->mass * kernel_root;
+  p->density.rho_dh -= hydro_dimension * p->mass * kernel_root;
   p->density.wcount += kernel_root;
 
   /* Finish the calculation by inserting the missing h-factors */
   p->rho *= h_inv_dim;
-  p->rho_dh *= h_inv_dim_plus_one;
+  p->density.rho_dh *= h_inv_dim_plus_one;
   p->density.wcount *= kernel_norm;
   p->density.wcount_dh *= h_inv * kernel_gamma * kernel_norm;
 
   const float rho_inv = 1.f / p->rho;
 
-  /* Compute the derivative term */
-  p->rho_dh = 1.f / (1.f + hydro_dimension_inv * p->h * p->rho_dh * rho_inv);
-
   /* Finish calculation of the velocity curl components */
   p->density.rot_v[0] *= h_inv_dim_plus_one * rho_inv;
   p->density.rot_v[1] *= h_inv_dim_plus_one * rho_inv;
@@ -276,7 +273,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   const float rho_inv = 1.f / p->rho;
 
   /* Divide the pressure by the density and density gradient */
-  const float P_over_rho2 = pressure * rho_inv * rho_inv * p->rho_dh;
+  const float P_over_rho2 = pressure * rho_inv * rho_inv;
 
   /* Compute the sound speed */
   const float soundspeed = sqrtf(hydro_gamma * pressure * rho_inv);
@@ -285,7 +282,12 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   const float balsara =
       abs_div_v / (abs_div_v + curl_v + 0.0001f * soundspeed / fac_mu / p->h);
 
+  /* Compute the "grad h" term */
+  const float grad_h_term =
+      1.f / (1.f + hydro_dimension_inv * p->h * p->density.rho_dh * rho_inv);
+
   /* Update variables. */
+  p->force.f = grad_h_term;
   p->force.P_over_rho2 = P_over_rho2;
   p->force.soundspeed = soundspeed;
   p->force.balsara = balsara;
@@ -352,7 +354,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
   const float rho_inv = 1.f / p->rho;
 
   /* Divide the pressure by the density and density gradient */
-  const float P_over_rho2 = pressure * rho_inv * rho_inv * p->rho_dh;
+  const float P_over_rho2 = pressure * rho_inv * rho_inv;
 
   /* Compute the new sound speed */
   const float soundspeed = sqrtf(hydro_gamma * pressure * rho_inv);
@@ -400,6 +402,16 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
   /* Do not 'overcool' when timestep increases */
   if (p->entropy + p->entropy_dt * half_dt < 0.5f * p->entropy)
     p->entropy_dt = -0.5f * p->entropy / half_dt;
+
+  /* Compute the pressure */
+  const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+
+  /* Compute the sound speed from the pressure*/
+  const float rho_inv = 1.f / p->rho;
+  const float soundspeed = sqrtf(hydro_gamma * pressure * rho_inv);
+
+  p->force.soundspeed = soundspeed;
+  p->force.P_over_rho2 = pressure * rho_inv * rho_inv;
 }
 
 /**
@@ -414,6 +426,16 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
 
   /* We read u in the entropy field. We now get S from u */
   p->entropy = gas_entropy_from_internal_energy(p->rho, p->entropy);
+
+  /* Compute the pressure */
+  const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+
+  /* Compute the sound speed from the pressure*/
+  const float rho_inv = 1.f / p->rho;
+  const float soundspeed = sqrtf(hydro_gamma * pressure * rho_inv);
+
+  p->force.soundspeed = soundspeed;
+  p->force.P_over_rho2 = pressure * rho_inv * rho_inv;
 }
 
 #endif /* SWIFT_GADGET2_HYDRO_H */

src/hydro/Gadget2/hydro_debug.h

@@ -30,8 +30,8 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
       "v_sig=%e dh/dt=%.3e t_begin=%d, t_end=%d\n",
       p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
       xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
-      p->h, p->density.wcount, p->density.wcount_dh, p->mass, p->rho_dh, p->rho,
-      hydro_get_pressure(p, 0.), p->force.P_over_rho2, p->entropy,
+      p->h, p->density.wcount, p->density.wcount_dh, p->mass, p->density.rho_dh,
+      p->rho, hydro_get_pressure(p, 0.), p->force.P_over_rho2, p->entropy,
       p->entropy_dt, p->force.soundspeed, p->density.div_v, p->density.rot_v[0],
       p->density.rot_v[1], p->density.rot_v[2], p->force.balsara,
       p->force.v_sig, p->force.h_dt, p->ti_begin, p->ti_end);

src/hydro/Gadget2/hydro_iact.h

@@ -59,7 +59,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
 
   /* Compute contribution to the density */
   pi->rho += mj * wi;
-  pi->rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
+  pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
 
   /* Compute contribution to the number of neighbours */
   pi->density.wcount += wi;
@@ -72,7 +72,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
 
   /* Compute contribution to the density */
   pj->rho += mi * wj;
-  pj->rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx);
+  pj->density.rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx);
 
   /* Compute contribution to the number of neighbours */
   pj->density.wcount += wj;
@@ -209,13 +209,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_density(
   /* Update particles. */
   for (k = 0; k < VEC_SIZE; k++) {
     pi[k]->rho += rhoi.f[k];
-    pi[k]->rho_dh -= rhoi_dh.f[k];
+    pi[k]->density.rho_dh -= rhoi_dh.f[k];
     pi[k]->density.wcount += wcounti.f[k];
     pi[k]->density.wcount_dh -= wcounti_dh.f[k];
     pi[k]->density.div_v -= div_vi.f[k];
     for (j = 0; j < 3; j++) pi[k]->density.rot_v[j] += curl_vi[j].f[k];
     pj[k]->rho += rhoj.f[k];
-    pj[k]->rho_dh -= rhoj_dh.f[k];
+    pj[k]->density.rho_dh -= rhoj_dh.f[k];
     pj[k]->density.wcount += wcountj.f[k];
     pj[k]->density.wcount_dh -= wcountj_dh.f[k];
     pj[k]->density.div_v -= div_vj.f[k];
@@ -254,7 +254,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
 
   /* Compute contribution to the density */
   pi->rho += mj * wi;
-  pi->rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
+  pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
 
   /* Compute contribution to the number of neighbours */
   pi->density.wcount += wi;
@@ -366,7 +366,7 @@ runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj,
   /* Update particles. */
   for (k = 0; k < VEC_SIZE; k++) {
     pi[k]->rho += rhoi.f[k];
-    pi[k]->rho_dh -= rhoi_dh.f[k];
+    pi[k]->density.rho_dh -= rhoi_dh.f[k];
     pi[k]->density.wcount += wcounti.f[k];
     pi[k]->density.wcount_dh -= wcounti_dh.f[k];
     pi[k]->density.div_v -= div_vi.f[k];
@@ -415,7 +415,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   kernel_deval(xj, &wj, &wj_dx);
   const float wj_dr = hjd_inv * wj_dx;
 
-  /* Compute gradient terms */
+  /* Compute h-gradient terms */
+  const float f_i = pi->force.f;
+  const float f_j = pj->force.f;
+
+  /* Compute pressure terms */
   const float P_over_rho2_i = pi->force.P_over_rho2;
   const float P_over_rho2_j = pj->force.P_over_rho2;
 
@@ -447,7 +451,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   /* Now, convolve with the kernel */
   const float visc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
   const float sph_term =
-      (P_over_rho2_i * wi_dr + P_over_rho2_j * wj_dr) * r_inv;
+      (f_i * P_over_rho2_i * wi_dr + f_j * P_over_rho2_j * wj_dr) * r_inv;
 
   /* Eventually got the acceleration */
   const float acc = visc_term + sph_term;
@@ -690,7 +694,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   kernel_deval(xj, &wj, &wj_dx);
   const float wj_dr = hjd_inv * wj_dx;
 
-  /* Compute gradient terms */
+  /* Compute h-gradient terms */
+  const float f_i = pi->force.f;
+  const float f_j = pj->force.f;
+
+  /* Compute pressure terms */
   const float P_over_rho2_i = pi->force.P_over_rho2;
   const float P_over_rho2_j = pj->force.P_over_rho2;
 
@@ -722,7 +730,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   /* Now, convolve with the kernel */
   const float visc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
   const float sph_term =
-      (P_over_rho2_i * wi_dr + P_over_rho2_j * wj_dr) * r_inv;
+      (f_i * P_over_rho2_i * wi_dr + f_j * P_over_rho2_j * wj_dr) * r_inv;
 
   /* Eventually got the acceleration */
   const float acc = visc_term + sph_term;

src/hydro/Gadget2/hydro_part.h

@@ -80,9 +80,6 @@ struct part {
   /* Entropy time derivative */
   float entropy_dt;
 
-  /* Derivative of the density with respect to smoothing length. */
-  float rho_dh;
-
   union {
 
     struct {
@@ -93,6 +90,9 @@ struct part {
       /* Number of neighbours spatial derivative. */
       float wcount_dh;
 
+      /* Derivative of the density with respect to h. */
+      float rho_dh;
+
       /* Particle velocity curl. */
       float rot_v[3];
 
@@ -106,6 +106,9 @@ struct part {
       /* Balsara switch */
       float balsara;
 
+      /*! "Grad h" term */
+      float f;
+
       /* Pressure over density squared (including drho/dh term) */
       float P_over_rho2;
 

src/hydro/Minimal/hydro.h

@@ -214,7 +214,7 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
   p->density.wcount = 0.f;
   p->density.wcount_dh = 0.f;
   p->rho = 0.f;
-  p->rho_dh = 0.f;
+  p->density.rho_dh = 0.f;
 }
 
 /**
@@ -240,19 +240,14 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
 
   /* Final operation on the density (add self-contribution). */
   p->rho += p->mass * kernel_root;
-  p->rho_dh -= hydro_dimension * p->mass * kernel_root;
+  p->density.rho_dh -= hydro_dimension * p->mass * kernel_root;
   p->density.wcount += kernel_root;
 
   /* Finish the calculation by inserting the missing h-factors */
   p->rho *= h_inv_dim;
-  p->rho_dh *= h_inv_dim_plus_one;
+  p->density.rho_dh *= h_inv_dim_plus_one;
   p->density.wcount *= kernel_norm;
   p->density.wcount_dh *= h_inv * kernel_gamma * kernel_norm;
-
-  const float irho = 1.f / p->rho;
-
-  /* Compute the derivative term */
-  p->rho_dh = 1.f / (1.f + hydro_dimension_inv * p->h * p->rho_dh * irho);
 }
 
 /**
@@ -276,7 +271,14 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
 
   const float half_dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
   const float pressure = hydro_get_pressure(p, half_dt);
+  const float rho_inv = 1.f / p->rho;
+
+  /* Compute the "grad h" term */
+  const float grad_h_term =
+      1.f / (1.f + hydro_dimension_inv * p->h * p->density.rho_dh * rho_inv);
 
+  /* Update variables. */
+  p->force.f = grad_h_term;
   p->force.pressure = pressure;
 }
 
@@ -379,6 +381,10 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
 
   /* Do not 'overcool' when timestep increases */
   if (p->u + p->u_dt * half_dt < 0.5f * p->u) p->u_dt = -0.5f * p->u / half_dt;
+
+  /* Compute the pressure */
+  const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
+  p->force.pressure = pressure;
 }
 
 /**
@@ -392,6 +398,11 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
  * @param p The particle to act upon
  */
 __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
-    struct part *restrict p) {}
+    struct part *restrict p) {
+
+  /* Compute the pressure */
+  const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
+  p->force.pressure = pressure;
+}
 
 #endif /* SWIFT_MINIMAL_HYDRO_H */

src/hydro/Minimal/hydro_debug.h

@@ -44,7 +44,7 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
       p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
       xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
       p->u, p->u_dt, p->force.v_sig, p->force.pressure, p->h, p->force.h_dt,
-      (int)p->density.wcount, p->mass, p->rho_dh, p->rho, p->ti_begin,
+      (int)p->density.wcount, p->mass, p->density.rho_dh, p->rho, p->ti_begin,
       p->ti_end);
 }
 

src/hydro/Minimal/hydro_iact.h

@@ -55,7 +55,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
   kernel_deval(xi, &wi, &wi_dx);
 
   pi->rho += mj * wi;
-  pi->rho_dh -= mj * (hydro_dimension * wi + xi * wi_dx);
+  pi->density.rho_dh -= mj * (hydro_dimension * wi + xi * wi_dx);
   pi->density.wcount += wi;
   pi->density.wcount_dh -= xi * wi_dx;
 
@@ -65,7 +65,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
   kernel_deval(xj, &wj, &wj_dx);
 
   pj->rho += mi * wj;
-  pj->rho_dh -= mi * (hydro_dimension * wj + xj * wj_dx);
+  pj->density.rho_dh -= mi * (hydro_dimension * wj + xj * wj_dx);
   pj->density.wcount += wj;
   pj->density.wcount_dh -= xj * wj_dx;
 }
@@ -100,7 +100,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
   kernel_deval(xi, &wi, &wi_dx);
 
   pi->rho += mj * wi;
-  pi->rho_dh -= mj * (hydro_dimension * wi + xi * wi_dx);
+  pi->density.rho_dh -= mj * (hydro_dimension * wi + xi * wi_dx);
   pi->density.wcount += wi;
   pi->density.wcount_dh -= xi * wi_dx;
 }
@@ -152,8 +152,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   const float wj_dr = hjd_inv * wj_dx;
 
   /* Compute gradient terms */
-  const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->rho_dh;
-  const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
+  const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->force.f;
+  const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->force.f;
 
   /* Compute dv dot r. */
   const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
@@ -263,8 +263,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   const float wj_dr = hjd_inv * wj_dx;
 
   /* Compute gradient terms */
-  const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->rho_dh;
-  const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
+  const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->force.f;
+  const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->force.f;
 
   /* Compute dv dot r. */
   const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +

src/hydro/Minimal/hydro_part.h

@@ -93,9 +93,6 @@ struct part {
   /*! Particle density. */
   float rho;
 
-  /*! Derivative of density with respect to h */
-  float rho_dh;
-
   /* Store density/force specific stuff. */
   union {
 
@@ -114,6 +111,9 @@ struct part {
       /*! Derivative of the neighbour number with respect to h. */
       float wcount_dh;
 
+      /*! Derivative of density with respect to h */
+      float rho_dh;
+
     } density;
 
     /**
@@ -125,6 +125,9 @@ struct part {
      */
     struct {
 
+      /*! "Grad h" term */
+      float f;
+
       /*! Particle pressure. */
       float pressure;
 

tests/test27cells.c

@@ -199,10 +199,8 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
             hydro_get_density(&main_cell->parts[pid]),
 #if defined(GIZMO_SPH)
             0.f,
-#elif defined(HOPKINS_PE_SPH)
-            main_cell->parts[pid].density.rho_dh,
 #else
-            main_cell->parts[pid].rho_dh,
+            main_cell->parts[pid].density.rho_dh,
 #endif
             main_cell->parts[pid].density.wcount,
             main_cell->parts[pid].density.wcount_dh,
@@ -238,10 +236,8 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
               cj->parts[pjd].v[2], hydro_get_density(&cj->parts[pjd]),
 #if defined(GIZMO_SPH)
               0.f,
-#elif defined(HOPKINS_PE_SPH)
-              main_cell->parts[pjd].density.rho_dh,
 #else
-              main_cell->parts[pjd].rho_dh,
+              main_cell->parts[pjd].density.rho_dh,
 #endif
               cj->parts[pjd].density.wcount, cj->parts[pjd].density.wcount_dh,
 #if defined(GADGET2_SPH) || defined(DEFAULT_SPH) || defined(HOPKINS_PE_SPH)

tests/testInteractions.c

@@ -111,13 +111,8 @@ void dump_indv_particle_fields(char *fileName, struct part *p) {
           "%13e %13e %13e %10f\n",
           p->id, p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2],
           p->a_hydro[0], p->a_hydro[1], p->a_hydro[2], p->rho,
-#ifdef HOPKINS_PE_SPH
-          p->density.rho_dh,
-#else
-          p->rho_dh,
-#endif
-          p->density.wcount, p->density.wcount_dh, p->force.h_dt,
-          p->force.v_sig,
+          p->density.rho_dh, p->density.wcount, p->density.wcount_dh,
+          p->force.h_dt, p->force.v_sig,
 #if defined(GADGET2_SPH)
           p->density.div_v, p->density.rot_v[0], p->density.rot_v[1],
           p->density.rot_v[2], p->entropy_dt

tests/testPair.c

@@ -138,10 +138,8 @@ void dump_particle_fields(char *fileName, struct cell *ci, struct cell *cj) {
             ci->parts[pid].v[2], hydro_get_density(&ci->parts[pid]),
 #if defined(GIZMO_SPH)
             0.f,
-#elif defined(HOPKINS_PE_SPH)
-            ci->parts[pid].density.rho_dh,
 #else
-            ci->parts[pid].rho_dh,
+            ci->parts[pid].density.rho_dh,
 #endif
             ci->parts[pid].density.wcount, ci->parts[pid].density.wcount_dh,
 #if defined(GADGET2_SPH) || defined(DEFAULT_SPH) || defined(HOPKINS_PE_SPH)
@@ -164,10 +162,8 @@ void dump_particle_fields(char *fileName, struct cell *ci, struct cell *cj) {
             cj->parts[pjd].v[2], hydro_get_density(&cj->parts[pjd]),
 #if defined(GIZMO_SPH)
             0.f,
-#elif defined(HOPKINS_PE_SPH)
-            cj->parts[pjd].density.rho_dh,
 #else
-            cj->parts[pjd].rho_dh,
+            cj->parts[pjd].density.rho_dh,
 #endif
             cj->parts[pjd].density.wcount, cj->parts[pjd].density.wcount_dh,
 #if defined(GADGET2_SPH) || defined(DEFAULT_SPH) || defined(HOPKINS_PE_SPH)

tests/tolerance_27_normal.dat

 #   ID      pos_x      pos_y      pos_z        v_x        v_y        v_z           rho        rho_dh        wcount     wcount_dh         div_v       curl_vx       curl_vy       curl_vz
     0	    1e-6       1e-6	  1e-6 	       1e-6 	  1e-6	     1e-6	   1e-6	      2e-6	    2e-5       2e-3		 2e-6	     2e-6	   2e-6		 2e-6
-    0	    1e-6       1e-6	  1e-6 	       1e-6 	  1e-6	     1e-6	   1e-6	      1e-6	    1e-5       1e-4		 2e-5	     2e-5	   2e-5	 	 2e-5
+    0	    1e-6       1e-6	  1e-6 	       1e-6 	  1e-6	     1e-6	   1e-6	      2e-6	    1e-5       1e-4		 2e-5	     2e-5	   2e-5	 	 2e-5

tests/tolerance_27_perturbed.dat

 #   ID      pos_x      pos_y      pos_z        v_x        v_y        v_z           rho        rho_dh        wcount     wcount_dh         div_v       curl_vx       curl_vy       curl_vz
     0	    1e-6       1e-6	  1e-6 	       1e-6 	  1e-6	     1e-6	   1.2e-6     2e-6	    2.1e-5     2e-3		 2.1e-6	     2e-6	   2e-6		 2e-6
-    0	    1e-6       1e-6	  1e-6 	       1e-6 	  1e-6	     1e-6	   1e-6	      1e-6	    1e-5       1e-4		 2e-5	     4e-4	   4e-4	 	 4e-4
+    0	    1e-6       1e-6	  1e-6 	       1e-6 	  1e-6	     1e-6	   1e-6	      1e-3	    1e-5       1e-4		 2e-5	     4e-4	   4e-4	 	 4e-4