Template Struct ResistiveSource

Defined in File mhd_resistivity.hpp

Struct Documentation

template<typename Real, typename Backend> struct ResistiveSource

Explicit resistivity (magnetic diffusion) as a source term.

Adds -curl(eta curl B) to the induction equation and the corresponding energy flux div((1/4pi) eta (curl B) x B) to the energy equation (the Joule heating goes to the internal energy). In index notation, dB_i/dt = d_j G_ij, dE/dt = d_j ((1/4pi) B_i G_ij), G_ij = eta (d_j B_i - d_i B_j), discretized in the flux form with second order accuracy (see resistive_flux).

Not included by <miso/mhd.hpp>; include <miso/mhd_resistivity.hpp> explicitly.

This struct holds views: the diffusivity array eta (at cell centers, including ghost cells) is owned by the caller and must outlive this object. eta may be changed during the run (e.g., time-dependent eta); dt_limit() recomputes the limit every step.

struct Model : public mhd::ModelBase<Model, Real, Backend> {
  ...
  Array3D<Real, Backend> eta;
  mhd::ResistiveSource<Real, Backend> src;
  Model(Config &config)
      : ModelBase(config), ..., eta(make_eta(...)),
        src(config, mhd.grid, eta) {}
};

To combine it with other source terms, hold it in a user-defined source term and forward bx, by, bz, ei and dt_limit(). Without dt_limit(), neither the stability limit of the diffusion nor the check of eta is applied, and no error is reported.

struct MySource {
  mhd::ResistiveSource<Real, Backend> resistivity;
  __host__ __device__ Real vz(FieldsView<const Real> qq, int i, int j,
int k) const { return -g * qq.ro(i, j, k); }
  __host__ __device__ Real bx(FieldsView<const Real> qq, int i, int j,
int k) const {
return resistivity.bx(qq, i, j, k);
  }
  ...  // by, bz, ei in the same way
  Real dt_limit() const { return resistivity.dt_limit(); }
};

Note

Ghost cells of the MHD fields including edges and corners are used.

Public Functions

inline ResistiveSource(const Grid<Real, Backend> &grid_, const Array3D<Real, Backend> &eta_, Real cfl_number_)

Throws:: std::runtime_error – if eta is negative or not finite somewhere on this rank

inline ResistiveSource(const Config &config, const Grid<Real, Backend> &grid_, const Array3D<Real, Backend> &eta_): Construct with mhd.resistivity.cfl_number in the config.

inline Real bx(const FieldsView<const Real> &qq, int i, int j, int k) const: x induction: d_y G_xy + d_z G_xz

inline Real by(const FieldsView<const Real> &qq, int i, int j, int k) const: y induction: d_x G_yx + d_z G_yz

inline Real bz(const FieldsView<const Real> &qq, int i, int j, int k) const: z induction: d_x G_zx + d_y G_zy

inline Real ei(const FieldsView<const Real> &qq, int i, int j, int k) const: Energy: divergence of the energy flux (1/4pi) B_m G_md.

inline Real dt_limit() const

Upper limit of the time step on this rank: cfl_number / max(eta*sum(1/ds^2)).

Recomputed at every call, so that a change of eta is followed. ModelBase::update takes the minimum over all ranks.

Returns:: -1 if eta is negative or not finite somewhere on this rank (ModelBase::update then throws on every rank)

inline Real max_rate() const: max(eta*sum(1/ds^2)) over the interior cells of this rank, or -1 if eta is negative or not finite in any cell (including ghost cells)

Note

Public because nvcc does not allow extended device lambdas in private member functions.

Public Members

GridView<const Real> grid: Grid (of the backend)

Array3DView<const Real> eta: Magnetic diffusivity at cell centers (including ghost cells)

Real cfl_number

Safety factor of the diffusive time step.

Stability of RK4 requires lambda*dt < 2.78 with lambda <= 4*eta*sum(1/ds^2), i.e., cfl_number < 0.69.

impl_resistivity::ReduceWorkspace<Real, Backend> workspace: Workspace of the reduction in dt_limit() (used on the host only)