Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <cstddef>
       7             : 
       8             : #include "DataStructures/Tensor/TypeAliases.hpp"
       9             : #include "DataStructures/VariablesTag.hpp"
      10             : #include "Evolution/Systems/ForceFree/BoundaryConditions/BoundaryCondition.hpp"
      11             : #include "Evolution/Systems/ForceFree/Characteristics.hpp"
      12             : #include "Evolution/Systems/ForceFree/Tags.hpp"
      13             : #include "Evolution/Systems/ForceFree/TimeDerivativeTerms.hpp"
      14             : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
      15             : 
      16             : ///\cond
      17             : class DataVector;
      18             : ///\endcond
      19             : 
      20             : /*!
      21             :  * \ingroup EvolutionSystemsGroup
      22             :  * \brief Items related to evolving the GRFFE system with divergence cleaning
      23             :  *
      24             :  */
      25             : namespace ForceFree {
      26             : 
      27             : /*!
      28             :  * \brief General relativistic force-free electrodynamics (GRFFE) system with
      29             :  * divergence cleaning.
      30             :  *
      31             :  * For electromagnetism in a curved spacetime, Maxwell equations are given as
      32             :  * \f{align*}{
      33             :  *  \nabla_a F^{ab} &= -J^b,\\
      34             :  *  \nabla_a ^* F^{ab} &= 0.
      35             :  * \f}
      36             :  *
      37             :  * where \f$F^{ab}\f$ is the electromagnetic field tensor, \f$^*F^{ab}\f$ is its
      38             :  * dual \f$^*F^{ab} = \epsilon^{abcd}F_{cd} / 2\f$, and \f$J^a\f$ is the
      39             :  * 4-current.
      40             :  *
      41             :  * \note
      42             :  * - We are using the electromagnetic variables with the scaling convention
      43             :  * that the factor \f$4\pi\f$ does not appear in Maxwell equations and the
      44             :  * stress-energy tensor of the EM fields (geometrized Heaviside-Lorentz units).
      45             :  * - We adopt following definition of the Levi-Civita tensor by
      46             :  * \cite Misner1973
      47             :  * \f{align*}{
      48             :  *  \epsilon_{abcd} &= \sqrt{-g} \, [abcd] ,\\
      49             :  *  \epsilon^{abcd} &= -\frac{1}{\sqrt{-g}} \, [abcd] ,
      50             :  * \f}
      51             :  * where \f$g\f$ is the determinant of spacetime metric, and \f$[abcd]=\pm 1\f$
      52             :  * is the antisymmetric symbol with \f$[0123]=+1\f$.
      53             :  *
      54             :  * In SpECTRE, we evolve 'extended' (or augmented) version of Maxwell equations
      55             :  * with two divergence cleaning scalar fields \f$\psi\f$ and \f$\phi\f$ :
      56             :  *
      57             :  * \f{align*}{
      58             :  *  \nabla_a(F^{ab}+g^{ab}\psi) & = -J^b + \kappa_\psi n^b \psi \\
      59             :  *  \nabla_a(^* F^{ab}+ g^{ab}\phi) & = \kappa_\phi n^b \phi
      60             :  * \f}
      61             :  *
      62             :  * which reduce to the original Maxwell equations when \f$\psi=\phi=0\f$. For
      63             :  * damping constants $\kappa_{\psi, \phi} > 0$, Gauss constraint violations are
      64             :  * damped with timescales $\kappa_{\psi,\phi}^{-1}$ and propagated away.
      65             :  *
      66             :  * We decompose the EM field tensor as follows
      67             :  * \f{align*}{
      68             :  *  F^{ab} = n^a E^b - n^b E^a - \epsilon^{abcd}B_c n_d,
      69             :  * \f}
      70             :  *
      71             :  * where $n^a$ is the normal to spatial hypersurface, $E^a$ and $B^a$ are
      72             :  * electric and magnetic fields.
      73             :  *
      74             :  * Evolved variables are
      75             :  *
      76             :  * \f{align*}{
      77             :  * \mathbf{U} = \sqrt{\gamma}\left[\,\begin{matrix}
      78             :  *      E^i \\
      79             :  *      B^i \\
      80             :  *      \psi \\
      81             :  *      \phi \\
      82             :  *      q \\
      83             :  * \end{matrix}\,\right] \equiv \left[\,\,\begin{matrix}
      84             :  *      \tilde{E}^i \\
      85             :  *      \tilde{B}^i \\
      86             :  *      \tilde{\psi} \\
      87             :  *      \tilde{\phi} \\
      88             :  *      \tilde{q} \\
      89             :  * \end{matrix}\,\,\right] \f}
      90             :  *
      91             :  * where \f$E^i\f$ is electric field, \f$B^i\f$ is magnetic field, $\psi$ is
      92             :  * electric divergence cleaning field, $\phi$ is magnetic divergence cleaning
      93             :  * field, \f$q\equiv-n_aJ^a\f$ is electric charge density, and \f$\gamma\f$ is
      94             :  * the determinant of spatial metric.
      95             :  *
      96             :  * Corresponding fluxes \f$\mathbf{F}^j\f$ are
      97             :  *
      98             :  * \f{align*}{
      99             :  *  F^j(\tilde{E}^i) & = -\beta^j\tilde{E}^i + \alpha
     100             :  *      (\gamma^{ij}\tilde{\psi} - \epsilon^{ijk}_{(3)}\tilde{B}_k) \\
     101             :  *  F^j(\tilde{B}^i) & = -\beta^j\tilde{B}^i + \alpha (\gamma^{ij}\tilde{\phi} +
     102             :  *      \epsilon^{ijk}_{(3)}\tilde{E}_k) \\
     103             :  *  F^j(\tilde{\psi}) & = -\beta^j \tilde{\psi} + \alpha \tilde{E}^j \\
     104             :  *  F^j(\tilde{\phi}) & = -\beta^j \tilde{\phi} + \alpha \tilde{B}^j \\
     105             :  *  F^j(\tilde{q}) & = \tilde{J}^j - \beta^j \tilde{q}
     106             :  * \f}
     107             :  *
     108             :  * and source terms are
     109             :  *
     110             :  * \f{align*}{
     111             :  *  S(\tilde{E}^i) &= -\tilde{J}^i - \tilde{E}^j \partial_j \beta^i
     112             :  *    + \tilde{\psi} ( \gamma^{ij} \partial_j \alpha - \alpha \gamma^{jk}
     113             :  *      \Gamma^i_{jk} ) \\
     114             :  *  S(\tilde{B}^i) &= -\tilde{B}^j \partial_j \beta^i + \tilde{\phi} (
     115             :  *      \gamma^{ij} \partial_j \alpha - \alpha \gamma^{jk} \Gamma^i_{jk} ) \\
     116             :  *  S(\tilde{\psi}) &= \tilde{E}^k \partial_k \alpha + \alpha \tilde{q} -
     117             :  *      \alpha \tilde{\phi} ( K + \kappa_\phi ) \\
     118             :  *  S(\tilde{\phi}) &= \tilde{B}^k \partial_k \alpha - \alpha \tilde{\phi}
     119             :  *      (K + \kappa_\phi ) \\
     120             :  *  S(\tilde{q}) &= 0
     121             :  * \f}
     122             :  *
     123             :  * where $\tilde{J}^i \equiv \alpha \sqrt{\gamma}J^i$.
     124             :  *
     125             :  * See the documentation of Fluxes and Sources for further details.
     126             :  *
     127             :  * In addition to Maxwell equations, general relativistic force-free
     128             :  * electrodynamics (GRFFE) assumes the following which are called the force-free
     129             :  * (FF) conditions.
     130             :  *
     131             :  * \f{align*}{
     132             :  *  F^{ab}J_b & = 0, \\
     133             :  *  ^*F^{ab}F_{ab} & = 0, \\
     134             :  *  F^{ab}F_{ab} & > 0.
     135             :  * \f}
     136             :  *
     137             :  * In terms of electric and magnetic fields, the FF conditions above read
     138             :  *
     139             :  * \f{align*}{
     140             :  *  E_iJ^i & = 0 , \\
     141             :  *  qE^i + \epsilon_{(3)}^{ijk} J_jB_k & = 0 , \\
     142             :  *  B_iE^i & = 0 , \\
     143             :  *  B^2 - E^2 & > 0.
     144             :  * \f}
     145             :  *
     146             :  * where \f$B^2=B^aB_a\f$ and \f$E^2 = E^aE_a\f$. Also,
     147             :  * \f$\epsilon_{(3)}^{ijk}\f$ is the spatial Levi-Civita tensor defined as
     148             :  *
     149             :  * \f{align*}
     150             :  *  \epsilon_{(3)}^{ijk} \equiv n_\mu \epsilon^{\mu ijk}
     151             :  *   = -\frac{1}{\sqrt{-g}} n_\mu [\mu ijk] = \frac{1}{\sqrt{\gamma}} [ijk]
     152             :  * \f}
     153             :  *
     154             :  * where \f$n^\mu\f$ is the normal to spatial hypersurface and \f$[ijk]\f$ is
     155             :  * the antisymmetric symbol with \f$[123] = +1\f$.
     156             :  *
     157             :  * There are a number of different ways in literature to numerically treat the
     158             :  * FF conditions. For the constraint $B_iE^i = 0$, cleaning of the parallel
     159             :  * electric field after every time step (e.g. \cite Palenzuela2010) or adopting
     160             :  * analytically determined parallel current density \cite Komissarov2011
     161             :  * were explored. On the magnetic dominance condition $B^2 - E^2 > 0$, there
     162             :  * have been approaches with modification of the drift current
     163             :  * \cite Komissarov2006 or manual rescaling of the electric field
     164             :  * \cite Palenzuela2010.
     165             :  *
     166             :  * We take the strategy that introduces special driver terms in the electric
     167             :  * current density \f$J^i\f$ following \cite Alic2012 :
     168             :  *
     169             :  * \f{align}{
     170             :  *  J^i = J^i_\mathrm{drift} + J^i_\mathrm{parallel}
     171             :  * \f}
     172             :  *
     173             :  * with
     174             :  *
     175             :  * \f{align}{
     176             :  *  J^i_\mathrm{drift} & = q \frac{\epsilon^{ijk}_{(3)}E_jB_k}{B_lB^l}, \\
     177             :  *  J^i_\mathrm{parallel} & = \eta \left[ \frac{E_jB^j}{B_lB^l}B^i
     178             :  *          + \frac{\mathcal{R}(E_lE^l-B_lB^l)}{B_lB^l}E^i \right] .
     179             :  * \f}
     180             :  *
     181             :  * where \f$\eta\f$ is the parallel conductivity and \f$\eta\rightarrow\infty\f$
     182             :  * corresponds to the ideal force-free limit. \f$\mathcal{R}(x)\f$ is the ramp
     183             :  * (or rectifier) function defined as
     184             :  *
     185             :  * \f{align*}
     186             :  *  \mathcal{R}(x) = \left\{\begin{array}{lc}
     187             :  *          x, & \text{if } x \geq 0 \\
     188             :  *          0, & \text{if } x < 0 \\
     189             :  * \end{array}\right\} = \max (x, 0) .
     190             :  * \f}
     191             :  *
     192             :  * Internally we handle each pieces \f$\tilde{J}^i_\mathrm{drift} \equiv
     193             :  * \alpha\sqrt{\gamma}J^i_\mathrm{drift}\f$ and \f$\tilde{J}^i_\mathrm{parallel}
     194             :  * \equiv \alpha\sqrt{\gamma}J^i_\mathrm{parallel}\f$ as two separate Tags
     195             :  * since the latter term is stiff and needs to be evolved in conjunction with
     196             :  * implicit time steppers.
     197             :  *
     198             :  */
     199           1 : struct System {
     200           0 :   static constexpr bool is_in_flux_conservative_form = true;
     201           0 :   static constexpr bool has_primitive_and_conservative_vars = false;
     202           0 :   static constexpr size_t volume_dim = 3;
     203             : 
     204           0 :   using boundary_conditions_base = BoundaryConditions::BoundaryCondition;
     205             : 
     206           0 :   using variables_tag =
     207             :       ::Tags::Variables<tmpl::list<Tags::TildeE, Tags::TildeB, Tags::TildePsi,
     208             :                                    Tags::TildePhi, Tags::TildeQ>>;
     209             : 
     210           0 :   using flux_variables = tmpl::list<Tags::TildeE, Tags::TildeB, Tags::TildePsi,
     211             :                                     Tags::TildePhi, Tags::TildeQ>;
     212             : 
     213           0 :   using non_conservative_variables = tmpl::list<>;
     214           0 :   using gradient_variables = tmpl::list<>;
     215             : 
     216           0 :   using spacetime_variables_tag =
     217             :       ::Tags::Variables<gr::tags_for_hydro<volume_dim, DataVector>>;
     218             : 
     219           0 :   using flux_spacetime_variables_tag = ::Tags::Variables<
     220             :       tmpl::list<gr::Tags::Lapse<DataVector>, gr::Tags::Shift<DataVector, 3>,
     221             :                  gr::Tags::SqrtDetSpatialMetric<DataVector>,
     222             :                  gr::Tags::SpatialMetric<DataVector, 3>,
     223             :                  gr::Tags::InverseSpatialMetric<DataVector, 3>>>;
     224             : 
     225           0 :   using compute_volume_time_derivative_terms = TimeDerivativeTerms;
     226             : 
     227           0 :   using compute_largest_characteristic_speed =
     228             :       Tags::LargestCharacteristicSpeedCompute;
     229             : 
     230           0 :   using inverse_spatial_metric_tag =
     231             :       gr::Tags::InverseSpatialMetric<DataVector, volume_dim>;
     232             : };
     233             : }  // namespace ForceFree