[docs]
class ShearAlfvenCurlCurl(StruphyModel):
"""Shear-Alfvén part of linear MHD with zero-flow equilibrium, with the magnetic field in :math:`H(\\textnormal{curl})`.
Same equations as :class:`~struphy.models.shear_alfven.ShearAlfven`, but with a weak induction equation, such that
the implicit time step leads to a curl-curl system for the magnetic field, which is preconditioned by the
Hiptmair-Xu preconditioner.
Parameters
----------
base_units: BaseUnits
Base units for normalization (default: BaseUnits())
mass_number: float
Mass number (in units of Proton mass) of the ion species (default: 1.0)
"""
@classmethod
def model_type(cls) -> LiteralOptions.ModelTypes:
return "Fluid"
## species
class EMFields(FieldSpecies):
def __init__(self):
self.b_field = FEECVariable(space="Hcurl")
self.init_variables()
class MHD(FluidSpecies):
def __init__(self, mass_number: float = 1.0):
self.velocity = FEECVariable(space="Hdiv")
self.init_variables(mass_number=mass_number)
class Propagators:
def __init__(self) -> None:
self.shear_alf = ShearAlfvenCurlCurlPropagator()
@property
def bulk_species(self):
return self.mhd
@property
def velocity_scale(self):
return "alfvén"
def post_allocate(self):
# project background magnetic field (1-form)
self._b_eq = Propagator.projected_equil.b1
# temporary vectors for scalar quantities
self._tmp_b1 = Propagator.derham.V1.zeros()
self._tmp_b2 = Propagator.derham.V1.zeros()
def __init__(self, base_units: BaseUnits = BaseUnits(), mass_number: float = 1.0):
# 0. store input parameters
self.params = copy.deepcopy(locals())
# 1. instantiate all species
self.em_fields = self.EMFields()
self.mhd = self.MHD(mass_number=mass_number)
# 2. derive units (must be done after instantiating species to access charge and mass numbers)
self.setup_equation_params(base_units=base_units)
# 3. instantiate all propagators
self.propagators = self.Propagators()
# 4. assign variables to propagators
self.propagators.shear_alf.variables.u = self.mhd.velocity
self.propagators.shear_alf.variables.b = self.em_fields.b_field
# 5. define scalars to be tracked during simulation
kinetic_energy = BilinearEnergyFEEC(self.mhd.velocity, bilinear_form_name="M2n")
magnetic_energy = BilinearEnergyFEEC(self.em_fields.b_field)
background_magnetic = FunctionScalarFEEC(self._compute_en_B_eq)
total_magnetic = FunctionScalarFEEC(self._compute_en_B_tot)
self.scalars = Scalars(
en_tot=kinetic_energy + magnetic_energy,
en_U=kinetic_energy,
en_B=magnetic_energy,
en_B_eq=background_magnetic,
en_B_tot=total_magnetic,
en_tot2=kinetic_energy + magnetic_energy + background_magnetic,
)
def _compute_en_B_eq(self):
Propagator.mass_ops.M1.dot(self._b_eq, apply_bc=False, out=self._tmp_b1)
return self._b_eq.inner(self._tmp_b1) / 2
def _compute_en_B_tot(self):
self._b_eq.copy(out=self._tmp_b1)
self._tmp_b1 += self.em_fields.b_field.spline.vector
Propagator.mass_ops.M1.dot(self._tmp_b1, apply_bc=False, out=self._tmp_b2)
return self._tmp_b1.inner(self._tmp_b2) / 2
[docs]
@classmethod
def doc_pde(cls):
r"""**PDEs solved by model:**
Momentum:
.. math::
\rho_0 \frac{\partial \tilde{\mathbf{U}}}{\partial t} = (\nabla \times \tilde{\mathbf{B}}) \times \mathbf{B}_0
Induction:
.. math::
\frac{\partial \tilde{\mathbf{B}}}{\partial t} - \nabla \times (\tilde{\mathbf{U}} \times \mathbf{B}_0) = 0
"""
[docs]
@classmethod
def doc_normalization(cls):
r"""All velocities are normalized by the Alfvén speed,
.. math::
\hat U = \hat v_A.
"""
[docs]
@classmethod
def doc_scalar_quantities(cls):
r"""**The following scalars are tracked during simulation:**
- Perturbed kinetic energy: ``en_U``
- Perturbed magnetic energy: ``en_B``
- Background magnetic energy: ``en_B_eq``
- Total magnetic energy including background: ``en_B_tot``
- Perturbation total energy: ``en_tot`` (not increased by the implicit time step: numerical damping that vanishes
with the resolution along the background field; exactly conserved in
:class:`~struphy.models.shear_alfven_curl_curl_with_proj.ShearAlfvenCurlCurlwithProj`)
- Total energy including background contribution: ``en_tot2``"""
[docs]
@classmethod
def doc_discretization(cls):
"""Time integration is performed by the following propagators (in sequence):
1. :class:`~struphy.propagators.shear_alfven_curl_curl_propagator.ShearAlfvenCurlCurlPropagator`
"""
doc = rf"""**1. ShearAlfvenCurlCurlPropagator:**
{ShearAlfvenCurlCurlPropagator.__doc__}
"""
return doc
[docs]
@classmethod
def doc_long_description(cls):
r"""ShearAlfvenCurlCurl solves the same equations as
:class:`~struphy.models.shear_alfven.ShearAlfven`, with the
perturbed magnetic field in :math:`H(\textnormal{curl})` instead of
:math:`H(\textnormal{div})`. The momentum equation is the
:math:`\rho_0`-weighted :math:`L^2` projection onto
:math:`H(\textnormal{div})`, and the induction equation is weak.
In the implicit time step, the velocity update is inserted into the
induction equation before the projection, which gives a field-aligned
curl-curl system for the magnetic field, preconditioned by the
Hiptmair-Xu auxiliary space preconditioner. This step is not exactly
energy conserving: it damps the waves numerically, with a damping that
vanishes quickly with the resolution along the background field (see
:class:`~struphy.models.shear_alfven_curl_curl_with_proj.ShearAlfvenCurlCurlwithProj`
for the exactly energy conserving version). With essential boundary
conditions in the de Rham sequence, the tangential components of the
magnetic field (a 1-form) vanish on that boundary (stronger than
:math:`\mathbf B \cdot \mathbf n = 0` of ShearAlfven)."""
[docs]
@classmethod
def doc_examples(cls):
r"""Create and initialize a shear-Alfvén curl-curl model:
.. code-block:: python
from struphy.models import ShearAlfvenCurlCurl
model = ShearAlfvenCurlCurl()
model.em_fields.b_field
model.mhd.velocity
"""
[docs]
@classmethod
def doc_use_cases(cls):
r"""This model is appropriate for:
- isolated shear-Alfvén wave propagation with large time steps
- testing curl-curl (Hiptmair-Xu) preconditioning of the shear-Alfvén step
- reduced linear-MHD studies without compressibility"""
[docs]
@classmethod
def doc_cannot_be_used_for(cls):
r"""This model is not suitable for:
- magnetosonic or compressible wave physics
- nonlinear MHD dynamics
- kinetic particle coupling
- resistive, viscous, or Hall-MHD effects
- polar splines with the Hiptmair-Xu preconditioner"""
