Correct tildeB sources and fluxes in ValenciaDivClean

docs/References.bib

@@ -2014,6 +2014,24 @@ @article{Palenzuela2010
   year      = {2010}
 }
 
+@article{Palenzuela2018,
+   title = "A Simflowny-based finite-difference code for high-performance
+           computing in numerical relativity",
+   volume = "35",
+   ISSN = "1361-6382",
+   url = "http://dx.doi.org/10.1088/1361-6382/aad7f6",
+   DOI = "10.1088/1361-6382/aad7f6",
+   number = "18",
+   journal = "Classical and Quantum Gravity",
+   publisher = "IOP Publishing",
+   author = "Palenzuela, Carlos and Miñano, Borja and Viganò, Daniele and
+            Arbona, Antoni and Bona-Casas, Carles and Rigo, Andreu and
+            Bezares, Miguel and Bona, Carles and Massó, Joan",
+   year = "2018",
+   month = "aug",
+   pages="185007"
+}
+
 @article{Pareschi2005,
   author        = {Pareschi, Lorenzo and Russo, Giovanni},
   title         = {Implicit-Explicit {Runge-Kutta} Schemes and Applications to

src/Evolution/Systems/GrMhd/ValenciaDivClean/Fluxes.cpp

@@ -65,9 +65,9 @@ void fluxes_impl(
       tilde_s_flux->get(i, j) = tilde_s.get(j) * transport_velocity->get(i) -
                                 lapse_b_over_w.get(j) * tilde_b.get(i);
       tilde_b_flux->get(i, j) =
-          tilde_b.get(j) * transport_velocity->get(i) +
-          get(lapse) * (get(tilde_phi) * inv_spatial_metric.get(i, j) -
-                        spatial_velocity.get(j) * tilde_b.get(i));
+          tilde_b.get(j) * transport_velocity->get(i) -
+          transport_velocity->get(j) * tilde_b.get(i) +
+          get(lapse) * (get(tilde_phi) * inv_spatial_metric.get(i, j));
     }
     tilde_s_flux->get(i, i) += get(pressure_star_lapse_sqrt_det_spatial_metric);
   }

src/Evolution/Systems/GrMhd/ValenciaDivClean/Fluxes.hpp

@@ -66,7 +66,7 @@ void fluxes_impl(
  * B^m v_m\\
  * F^i({\tilde \tau}) = &~  {\tilde \tau} v^i_{tr} + \sqrt{\gamma} \alpha \left(
  * p + p_m \right) v^i - \alpha {\tilde B}^i B^m v_m \\
- * F^i({\tilde B}^j) = &~  {\tilde B}^j v^i_{tr} - \alpha v^j {\tilde B}^i +
+ * F^i({\tilde B}^j) = &~  {\tilde B}^j v^i_{tr} - v^j_{tr} {\tilde B}^i +
  * \alpha \gamma^{ij} {\tilde \Phi} \\
  * F^i({\tilde \Phi}) = &~ \alpha {\tilde B^i} - \beta^i {\tilde \Phi}
  * \f}

src/Evolution/Systems/GrMhd/ValenciaDivClean/Sources.cpp

@@ -157,9 +157,6 @@ void sources_impl(
     source_tilde_b->get(i) *= get(tilde_phi);
     source_tilde_b->get(i) -=
         get(lapse) * get(tilde_phi) * trace_spatial_christoffel_second.get(i);
-    for (size_t m = 0; m < 3; ++m) {
-      source_tilde_b->get(i) -= tilde_b.get(m) * d_shift.get(m, i);
-    }
   }
 
   trace(source_tilde_phi, extrinsic_curvature, inv_spatial_metric);

src/Evolution/Systems/GrMhd/ValenciaDivClean/Sources.hpp

@@ -91,7 +91,7 @@ void sources_impl(
  * \alpha \\
  * S({\tilde \tau}) = & \alpha {\tilde S}^{mn} K_{mn} - {\tilde S}^m \partial_m
  * \alpha \\
- * S({\tilde B}^i) = & -{\tilde B}^m \partial_m \beta^i + {\tilde \Phi}
+ * S({\tilde B}^i) = & {\tilde \Phi}
  * \gamma^{im} \partial_m \alpha + \alpha {\tilde \Phi} \left( \frac{1}{2}
  * \gamma^{il} \gamma^{jk} - \gamma^{ij} \gamma^{lk} \right) \partial_l
  * \gamma_{jk} \\ S({\tilde \Phi}) = & {\tilde B}^k \partial_k \alpha - \alpha K

src/Evolution/Systems/GrMhd/ValenciaDivClean/System.hpp

@@ -32,6 +32,8 @@ namespace grmhd {
 /// - Black hole-neutron star mergers with a hot nuclear equation of state:
 /// outflow and neutrino-cooled disk for a low-mass, high-spin case
 /// \cite Deaton2013
+/// - A Simflowny-based finite-difference code for high-performance computing
+/// in Relativity \cite Palenzuela2018
 ///
 namespace ValenciaDivClean {
 
@@ -134,7 +136,7 @@ namespace ValenciaDivClean {
  *      \delta^i_j - \alpha b_j \tilde{B}^i / W \\
  * F^i({\tilde \tau}) &= {\tilde \tau} v^i_{tr} + \alpha \sqrt{\gamma} p^* v^i
  *      - \alpha^2 b^0 \tilde{B}^i / W \\
- * F^i({\tilde B}^j) &= {\tilde B}^j v^i_{tr} - \alpha v^j {\tilde B}^i
+ * F^i({\tilde B}^j) &= {\tilde B}^j v^i_{tr} - v^j_{tr} {\tilde B}^i
  *      + \alpha \gamma^{ij} {\tilde \Phi} \\
  * F^i({\tilde \Phi}) &= \alpha {\tilde B^i} - {\tilde \Phi} \beta^i
  * \f}
@@ -149,8 +151,7 @@ namespace ValenciaDivClean {
  *      + {\tilde \tau}) \partial_j \alpha \\
  * S({\tilde \tau}) &= \alpha {\tilde S}^{mn} K_{mn}
  *      - {\tilde S}^m \partial_m \alpha \\
- * S({\tilde B}^j) &= -{\tilde B}^m \partial_m \beta^j
- *      + \Phi \partial_k (\alpha \sqrt{\gamma}\gamma^{jk}) \\
+ * S({\tilde B}^j) &= \Phi \partial_k (\alpha \sqrt{\gamma}\gamma^{jk}) \\
  * S({\tilde \Phi}) &= {\tilde B}^k \partial_k \alpha - \alpha K
  *      {\tilde \Phi} - \alpha \kappa {\tilde \Phi}
  * \f}
@@ -166,6 +167,9 @@ namespace ValenciaDivClean {
  * divergence-free (no-monopole) condition \f$\Phi = \partial_i {\tilde B}^i =
  * 0\f$.
  *
+ * \note The fluxes and sources of ${\tilde B}$ follow eq. 32 of
+ * \cite Palenzuela2018 rather than eq. 87 of \cite Moesta2014
+ *
  * \note On the electron fraction side, the source term is currently set to
  * \f$S(\tilde{Y}_e) = 0\f$. Implementing the source term using neutrino scheme
  * is in progress (Last update : Oct 2022).

tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/Fluxes.py

@@ -159,8 +159,8 @@ def tilde_b_flux(
     magnetic_field,
 ):
     result = np.outer(lapse * spatial_velocity - shift, tilde_b)
+    result -= np.outer(tilde_b, lapse * spatial_velocity - shift)
     result += lapse * inv_spatial_metric * tilde_phi
-    result -= lapse * np.outer(tilde_b, spatial_velocity)
     return result
 
 

tests/Unit/Evolution/Systems/GrMhd/ValenciaDivClean/Sources.py

@@ -159,15 +159,11 @@ def source_tilde_b(
 ):
     term_one = np.einsum("ab, iab", inv_spatial_metric, d_spatial_metric)
     term_two = np.einsum("ab, aib", inv_spatial_metric, d_spatial_metric)
-    return (
-        tilde_phi * np.einsum("a, ia", d_lapse, inv_spatial_metric)
-        - np.einsum("a, ai", tilde_b, d_shift)
-        + lapse
-        * tilde_phi
-        * (
-            0.5 * np.einsum("ia, a", inv_spatial_metric, term_one)
-            - np.einsum("ia, a", inv_spatial_metric, term_two)
-        )
+    return tilde_phi * np.einsum(
+        "a, ia", d_lapse, inv_spatial_metric
+    ) + lapse * tilde_phi * (
+        0.5 * np.einsum("ia, a", inv_spatial_metric, term_one)
+        - np.einsum("ia, a", inv_spatial_metric, term_two)
     )