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 std::string name() { return "ForceFree"; }
201 :
202 0 : static constexpr bool is_in_flux_conservative_form = true;
203 0 : static constexpr bool has_primitive_and_conservative_vars = false;
204 0 : static constexpr size_t volume_dim = 3;
205 :
206 0 : using boundary_conditions_base = BoundaryConditions::BoundaryCondition;
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
|
