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