WARPXM Application Development Information

This document provides an explanation of adding a physics application through example. Here we discuss the augmentation of the ideal MHD equations for the inclusion of divergence cleaning, as per Dedner.

WARPXM is designed to solve partial differential equations describing flow a plasma in an from an Eulerian point of view, written in conservative form. Thus, for a conserved variable \(\boldsymbol{q}\), equations are written in the form

\begin{align*} \frac{\partial \boldsymbol{q}}{\partial t} + \nabla\cdot \mathcal{F}(\boldsymbol{q}) =& \boldsymbol{S} \end{align*}

where \(\mathcal{F}\) represents a flux and \(\boldsymbol{S}\) represents a source. In terms of actual index notation this is written

\begin{align*} \frac{\partial q_{i}}{\partial t} + \frac{\partial \mathcal{F}_{ij}}{\partial x_{j}} =& S_{i} \end{align*}

These equations are solved on an unstructured grid, using the discontinuous Galerkin finite elmement method. When implementing a physics application, intimate details of the method are not required, but instead the flux and terms of the particular model must be written in. Specifically one needs to implement the flux term, dubbed "internal flux", tthe source term, and the "numerical flux". The numerical flux, as the name implies is the part that requires some understanding of the numerical method. This is a resolution of internal flux at element boundaries which must be consistent between elements.

Development Process

1. Create a files generating internal flux, numerical flux, and sources
2. Register the files
3. Add to CMake
4. Add to WARPY
5. Add warpy function to init
6. Add warpy function to CMake

Example 1 - Linear Advection

In this case our system for a single component q linearly advected with no source, written as:

\begin{align*} \frac{\partial q}{\partial t} + \nabla\cdot \left(\boldsymbol{c}q\right) =& 0 \end{align*}

or

\begin{align*} \frac{\partial q}{\partial t} + \frac{\partial \left(\mathcal{c}_{1} q\right)}{\partial x} + \frac{\partial \left(\mathcal{c}_{2} q\right)}{\partial y} + \frac{\partial \left(\mathcal{c}_{3} q\right)}{\partial z} =& 0 \end{align*}

or

\begin{align*} \frac{\partial q}{\partial t} + \begin{pmatrix} c_{1}q & c_{2}q & c_{3}q \end{pmatrix} \begin{pmatrix} \frac{\partial}{\partial x} \\ \frac{\partial}{\partial y} \\ \frac{\partial}{\partial z} \end{pmatrix} =& 0 \end{align*}

Implementation of the internal flux, \(\mathcal{F}\), looks is written as a function developing the 2-dimensinal flux array of [number of dimensions] x [number of conserved variables] for a given set of conserved variables \(q\) and auxialliary variables \(\text{aux}\). NOTE: This is actually the transpose of the form as written in the flux equation (i is conserved variables, j is dimensions).

1. Internal and numerical flux

   (src/apps/simple/advection.h)

   #ifndef WXM_APPS_ADVECTION_H
   #define WXM_APPS_ADVECTION_H
    
   // Wm includes
   #include "apps/wmapplication.h"
    
   namespace wxm
   {
   namespace apps
   {
    
   class advection_t : public WmApplication
   {
   public:
       advection_t();
    
       ~advection_t();
    
       void setup(const WxCryptSet& wxc) override;
    
       const std::vector<int>& getInputVariableIndexes(int flag) const override
       {
           return _variables;
       }
    
       const std::vector<int>& getOutputVariableIndexes(int flag) const override
       {
           return _variables;
       }
    
       real numerical_flux(const real* q_l, const real* q_r, const real* aux_l,
           const real* aux_r, const solverVariables_t* pFV,
           real* numericalFlux) const override;
       real internal_flux(const real* q, const real* aux, const solverVariables_t* pSV,
           std::vector<std::vector<real>>& internalFlux) const override;
    
   protected:
       std::vector<real> _velocity;
    
       std::vector<int> _variables;
    
   private:
       advection_t& operator=(const advection_t& var);
       advection_t(const advection_t& var);
   };
   }
   }
    
   #endif // WXM_APPS_ADVECTION_H

   (src/apps/simple/advection.cc)

   #include "advection.h"
    
   // STL includes
   #include <cmath>
   #include <iostream>
    
   // Wm includes
   #include "lib/wmudgfunctions.h"
   #include "solver/wmsolver.h"
    
   namespace wxm
   {
   namespace apps
   {
    
   advection_t::advection_t()
   {
       _allowedFlags.push_back(WMAPPLICATIONFLAG_NUMERICALFLUX);
       _allowedFlags.push_back(WMAPPLICATIONFLAG_INTERNALFLUX);
   }
    
   advection_t::~advection_t()
   {
   }
    
   void advection_t::setup(const WxCryptSet& wxc)
   {
       WmApplication::setup(wxc);
    
       _velocity = wxc.getVec<real>("Velocity");
       _velocity.resize(3, 0.);
    
       _variables = WmApplication_Variable(wxc, "Density", 1).getIndexes();
   }
    
   real advection_t::numerical_flux(const real* q_l, const real* q_r, const real* aux_l,
       const real* aux_r, const solverVariables_t* pSV, real* numericalFlux) const
   {
       real velocity_x = _velocity[0];
       real velocity_y = _velocity[1];
       real velocity_z = _velocity[2];
    
       const real nc =
           velocity_x * pSV->R[0][0] + velocity_y * pSV->R[0][1] + velocity_z * pSV->R[0][2];
       const real absc = fabs(nc);
    
       // upwinding
       numericalFlux[0] = 0.5 * (q_r[0] * (nc - absc) + q_l[0] * (nc + absc));
    
       // Rusanov
       // numericalFlux[0] =
       //    0.5 * nc * (q_l[0] + q_r[0]) - 0.5 * absc * (q_l[0] - q_r[0]);
    
       // // central difference
       // numericalFlux[0] = 0.5 * (q_r[0] + q_l[0] );
    
    
       return pSV->cfl_max * pSV->dxn / absc;
   }
    
   real advection_t::internal_flux(const real* q, const real* aux,
       const solverVariables_t* pSV, std::vector<std::vector<real>>& internalFlux) const
   {
       for (size_t i = 0; i < pSV->num_dims; ++i)
       {
           internalFlux[i][0] = _velocity[i] * q[0];
       }
       // Numerical flux should return the suggested dt
       return INFINITY;
   }
   }
   }
    
   // Add app to registry
   #include "lib/wxcreator.h"
   REGISTER_CREATOR("advection", wxm::apps::advection_t, WmApplication);

2. Registration - The last lines in the advection.cc file registers the scheme:

   (src/apps/simple/advection.cc)

   // Add app to registry
   #include "lib/wxcreator.h"
   REGISTER_CREATOR("advection", wxm::apps::advection_t, WmApplication);

3. Add to Cmake

   (src/apps/simple/CMakeLists.txt)

   add_sources(
          advection.cc
   )

4. Add to warpy

   (tools/warpy/apps/simple.py)

   from .application import application
    
   class advection(application):
       '''
       @brief Advection
       @param name Name of module
       @param variable Variable
       @param velocity Velocity (list of 1, 2, or 3 values)
       @param components Components of variable
       '''
       def __init__(self, name, velocity, variable, components=None):
           super(advection, self).__init__(
               name,
               req_attrs=[('Kind', 'advection'),
                          ('Velocity', velocity)],
               variables=[('Density', variable, components)])

5. Add warpy function to init

   (tools/warpy/apps/__init__.py)

   from .simple import advection

6. Add warpy function to CMake

   (tools/warpy/apps/CMakeLists.txt)

   WARPY_FILES(__init__.py
       application.py
       simple.py
       )

Example 2 - Ideal MHD

WARPXM has an Ideal MHD module upon which additional terms can be added, such as resistive diffusion, Hall terms, divergence cleaning, viscous terms, etc. The base Ideal MHD equations can be written as as conservation of:

Continuity

\begin{align*} \frac{\partial \rho}{\partial t} + \nabla\cdot\left(\rho \boldsymbol{u}\right) =& 0 \end{align*}

Momentum

\begin{align*} \frac{\partial\left(\rho \boldsymbol{u}\right)}{\partial t} + \nabla\cdot \left[\rho \boldsymbol{u}\boldsymbol{u} + p \mathcal{I} - \left(\boldsymbol{B}\boldsymbol{B}-\frac{1}{2}B^{2}\mathcal{I}\right) \right] =& 0 \end{align*}

Energy

\begin{align*} \frac{\partial e_{t}}{\partial t} + \nabla\cdot \left[\left(e_{t}+p + \frac{1}{2}B^{2}\right)\boldsymbol{u}-\left(\boldsymbol{u}\cdot\boldsymbol{B}\right)\boldsymbol{B} \right] =& 0 \end{align*}

Magnetic Field

\begin{align*} \frac{\partial \boldsymbol{B}}{\partial t} + \nabla\cdot\left(\boldsymbol{u}\boldsymbol{B} - \boldsymbol{B}\boldsymbol{u} \right) =& 0 \end{align*}

where \(\boldsymbol{u}=\left(u,v,w\right)\) and \(\boldsymbol{B}=\left(B_{x},B_{y},B_{z}\right)\), total energy \(e_{t}\) is related to pressure, \(p\), by

\begin{align*} e_{t} =& \frac{p}{\gamma-1} + \frac{1}{2}\rho \boldsymbol{u}\cdot\boldsymbol{u} + \frac{1}{2}\boldsymbol{B}\cdot\boldsymbol{B} = \frac{p}{\gamma-1} + \frac{1}{2}\rho \left(u^{2}+v^{2}+w^{2}\right) + \frac{1}{2}\left(B_{x}^{2}+B_{y}^{2}+B_{z}^2\right) \end{align*}

Altogether this is written

\begin{align*} \frac{\partial}{\partial t} \begin{pmatrix} \rho \\ \rho u \\ \rho v \\ \rho w \\ e_{t} \\ B_{x} \\ B_{y} \\ B_{z} \end{pmatrix} + \frac{\partial}{\partial x} \begin{pmatrix} \rho u \\ \rho u u - (B_{x} B_{x} - \frac{1}{2}B^{2}) + p \\ \rho u v - B_{x} B_{y} \\ \rho u w - B_{x} B_{z} \\ \left(e_{t} + p + \frac{1}{2}B^{2}\right)u - \left(\boldsymbol{B}\cdot\boldsymbol{u}\right)B_{x} \\ 0 \\ u B_{y} - B_{x} v \\ u B_{z} - B_{x} w \\ \end{pmatrix} + \frac{\partial}{\partial y} \begin{pmatrix} \rho v \\ \rho v u - B_{y} B_{x} \\ \rho v v - (B_{y} B_{y} - \frac{1}{2}B^{2}) + p \\ \rho v w - B_{y} B_{z} \\ \left(e_{t} + p + \frac{1}{2}B^{2}\right)v - \left(\boldsymbol{B}\cdot\boldsymbol{u}\right)B_{y} \\ v B_{x} - B_{y} u \\ 0 \\ v B_{z} - B_{y} w \\ \end{pmatrix} + \frac{\partial}{\partial z} \begin{pmatrix} \rho w \\ \rho w u - B_{z} B_{x} \\ \rho w v - B_{z} B_{y} \\ \rho w w - (B_{z} B_{z} - \frac{1}{2}B^{2}) + p \\ \left(e_{t} + p + \frac{1}{2}B^{2}\right)w - \left(\boldsymbol{B}\cdot\boldsymbol{u}\right)B_{z} \\ w B_{x} - B_{z} u \\ w B_{y} - B_{z} v \\ 0 \end{pmatrix} =& \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} \end{align*}

or

\begin{align*} \begin{pmatrix} \rho \\ \rho u \\ \rho v \\ \rho w \\ e_{t} \\ B_{x} \\ B_{y} \\ B_{z} \end{pmatrix} + \begin{pmatrix} \rho u & \rho v & \rho w \\ \rho u u - (B_{x} B_{x} - \frac{1}{2}B^{2}) + p & \rho v u - B_{y} B_{x} & \rho w u - B_{z} B_{x} \\ \rho u v - B_{x} B_{y} & \rho v v - (B_{y} B_{y} - \frac{1}{2}B^{2}) + p & \rho w v - B_{z} B_{y} \\ \rho u w - B_{x} B_{z} & \rho v w - B_{y} B_{z} & \rho w w - (B_{z} B_{z} - \frac{1}{2}B^{2}) + p \\ \left(e_{t} + p + \frac{1}{2}B^{2}\right)u - \left(\boldsymbol{B}\cdot\boldsymbol{u}\right)B_{x} & \left(e_{t} + p + \frac{1}{2}B^{2}\right)v - \left(\boldsymbol{B}\cdot\boldsymbol{u}\right)B_{y} & \left(e_{t} + p + \frac{1}{2}B^{2}\right)w - \left(\boldsymbol{B}\cdot\boldsymbol{u}\right)B_{z} \\ 0 & v B_{x} - B_{y} u & w B_{x} - B_{z} u \\ u B_{y} - B_{x} v & 0 & w B_{y} - B_{z} v \\ u B_{z} - B_{x} w & v B_{z} - B_{y} w & 0 \end{pmatrix} \begin{pmatrix} \frac{\partial}{\partial x} \\ \frac{\partial}{\partial y} \\ \frac{\partial}{\partial z} \end{pmatrix} =& \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} \end{align*}

NOTE: When we write the flux tensor in WARPXM we write its transpose (see code below).

1. Writing Internal Flux and Numerical Flux

   (src/apps/imhd/wmapplication_imhd_flux.h)

   #ifndef WMAPPLICATION_IMHD_FLUX_H
   #define WMAPPLICATION_IMHD_FLUX_H
    
   // WarpM includes
   #include "apps/wmapplication.h"
    
   // IMHD includes
   #include "apps/imhd/imhd.h"
    
   class WmApplication_IMHD_Flux : public WmApplication
   {
   public:
       WmApplication_IMHD_Flux();
    
       ~WmApplication_IMHD_Flux();
    
       void setup(const WxCryptSet& wxc) override;
    
       const std::vector<int>& getInputVariableIndexes(int flag) const override
       {
           return _inputVariables;
       }
       const std::vector<int>& getOutputVariableIndexes(int flag) const override
       {
           return _outputVariables;
       }
    
       // redefinition of fluxes from parent
       real numerical_flux(const real* q_l, const real* q_r, const real* aux_l,
           const real* aux_r, const solverVariables_t* pFV,
           real* numericalFlux) const override;
       real internal_flux(const real* q, const real* aux,
           const solverVariables_t* pSV,
           std::vector<std::vector<real>>& internalFlux) const override;
    
   protected:
       real _gas_gamma;
       
       std::string _numerical_flux_type; // this is the type of numerical flux you want to use (rusanov, roe, hll, hlld, etc.) 
    
       std::vector<int> _inputVariables;
       std::vector<int> _outputVariables;
    
    
       real _min_density_floor;
       real _min_pressure_floor;
       
    
    
    
   private:
       WmApplication_IMHD_Flux& operator=(const WmApplication_IMHD_Flux& var);
       WmApplication_IMHD_Flux(const WmApplication_IMHD_Flux& var);
   };
    
   #endif // WMAPPLICATION_IMHD_FLUX_H

   (src/apps/imhd/wmapplication_imhd_flux.cc)

   #include "wmapplication_imhd_flux.h"
    
   // STD libraries
   #include <iostream>
   #include <cmath>
    
   // Wm libraries
   #include "lib/wmudgfunctions.h"
    
   WmApplication_IMHD_Flux::WmApplication_IMHD_Flux()
   {
       _allowedFlags.push_back(WMAPPLICATIONFLAG_NUMERICALFLUX);
       _allowedFlags.push_back(WMAPPLICATIONFLAG_INTERNALFLUX);
   }
    
   WmApplication_IMHD_Flux::~WmApplication_IMHD_Flux()
   {
    
   }
    
   void
   WmApplication_IMHD_Flux::setup(const WxCryptSet& wxc)
   {
       WmApplication::setup(wxc);
    
       WmApplication_Variable field(wxc, "Fluid", 8);
    
    
    
       _gas_gamma = wxc.get<real>(GAS_GAMMA);
    
    
       _numerical_flux_type = "Rusanov";
       if (wxc.has("NumericalFluxType"))
       {
           _numerical_flux_type = wxc.get<std::string>("NumericalFluxType");
       }
    
    
       _min_density_floor = 0.0;
       if (wxc.has("MinDensity"))
       {
           _min_density_floor = wxc.get<real>("MinDensity");
       }
       
       _min_pressure_floor = 0.0;
       if (wxc.has("MinPressure"))
       {
           _min_pressure_floor = wxc.get<real>("MinPressure");
       }   
       
       
    
       _inputVariables = field.getIndexes();
       _outputVariables = field.getIndexes();
    
   }
    
   real
   WmApplication_IMHD_Flux::numerical_flux(const real * q_l,
                                           const real * q_r,
                                           const real * aux_l,
                                           const real * aux_r,
                                           const solverVariables_t * pFV,
                                           real * numericalFlux) const
   {
       // Define the constants
       constants_imhd_t pC;
       populate_constants_imhd_t(&pC, _gas_gamma);
       // ensure positivity of density and pressure
       real dens_q_l[8], dens_q_r[8];
       real mod_q_l[8], mod_q_r[8];    
       ensure_floor_density(q_l,8,_min_density_floor,dens_q_l);
       ensure_floor_density(q_r,8,_min_density_floor,dens_q_r);    
       ensure_floor_pressure(dens_q_l,8,_gas_gamma,_min_pressure_floor,mod_q_l);
       ensure_floor_pressure(dens_q_r,8,_gas_gamma,_min_pressure_floor,mod_q_r);   
       
       real c = imhd_numerical_flux(&pC,mod_q_l, mod_q_r, aux_l, aux_r, pFV, _numerical_flux_type, numericalFlux);
    
       return pFV->cfl_max * pFV->dxn / fabs(c);
   }
    
   real
   WmApplication_IMHD_Flux::internal_flux(const real * q,
                                          const real * aux,
                                          // const elementGeometry_t * pEG,
                                          const solverVariables_t * pSV,
                                          // real *internalFlux[8]
                                          std::vector< std::vector<real> >& internalFlux
       ) const
   {
       constants_imhd_t pC;
       populate_constants_imhd_t(&pC, _gas_gamma);
    
       // ensure positivity of density and pressure
       real dens_q[8]; 
       real mod_q[8];
       ensure_floor_density(q,8,_min_density_floor,dens_q);
       ensure_floor_pressure(dens_q,8,_gas_gamma,_min_pressure_floor,mod_q);
    
       imhd_internal_flux(&pC, pSV, mod_q, internalFlux);
       
       // Speeds should be calculated in numerical flux
       return INFINITY;
    
   }
    
    
   #include "lib/wxcreator.h"
    
   // IMHD flux
   REGISTER_CREATOR("iMHDFlux", WmApplication_IMHD_Flux, WmApplication);

   where we are using helper functions for internal and numerical fluxes. For the ideal mhd internal flux we have from the file (src/apps/imhd/imhd.
   

   void imhd_internal_flux(const constants_imhd_t * pC, const solverVariables_t * pSV, const real* q, std::vector< std::vector<real> >& internalFlux)
   {
    
       
       real gas_gamma1 = pC->gas_gamma - 1;
       
       real rho = q[0];
       real u = q[1]/rho;
       real v = q[2]/rho;
       real w = q[3]/rho;
       real E = q[4];
       real bx = q[5];
       real by = q[6];
       real bz = q[7];
    
       real B2 = dot_product(q+5);
       real V2 = dot_product(q+1)/rho/rho;                 
       real BdotV = dot_product(q+1,q+5)/rho;  
       real p = get_p(pC->gas_gamma,rho,E,V2,B2);
    
       real Eterm1 = E + p + 0.5*B2;
       real Eterm2 = BdotV;
    
       internalFlux[0][0] = rho*u;
       internalFlux[0][1] = rho*u*u - (bx*bx-0.5*B2) + p;
       internalFlux[0][2] = rho*u*v - bx*by;
       internalFlux[0][3] = rho*u*w - bx*bz;
       internalFlux[0][4] = Eterm1*u - Eterm2*bx;
       internalFlux[0][5] = 0;     // u*bx - bx*u 
       internalFlux[0][6] = u*by - bx*v;
       internalFlux[0][7] = u*bz - bx*w;
    
       if (pSV->num_dims > 1)
       {
           internalFlux[1][0] = rho*v;
           internalFlux[1][1] = rho*v*u - by*bx;
           internalFlux[1][2] = rho*v*v - (by*by-0.5*B2) + p;
           internalFlux[1][3] = rho*v*w - by*bz;
           internalFlux[1][4] = Eterm1*v - Eterm2*by;
           internalFlux[1][5] = v*bx - by*u;
           internalFlux[1][6] = 0;     // v*by - by*v
           internalFlux[1][7] = v*bz - by*w;
       }   // NUM_DIMS > 1
    
    
       if (pSV->num_dims > 2)
       {
           internalFlux[2][0] = rho*w;
           internalFlux[2][1] = rho*w*u - bz*bx;
           internalFlux[2][2] = rho*w*v - bz*by;
           internalFlux[2][3] = rho*w*w - (bz*bz-0.5*B2) + p;
           internalFlux[2][4] = Eterm1*w - Eterm2*bz;
           internalFlux[2][5] = w*bx - bz*u;
           internalFlux[2][6] = w*by - bz*v;
           internalFlux[2][7] = 0;     // w*bz - bz*w
       }   // NUM_DIMS > 2
    
    
   }

   The numerical flux for ideal mhd is given also with a helper
   
   (src/apps/imhd/imhd.cc)
   

   real imhd_numerical_flux(
       const constants_imhd_t * pC,
       const real * q_l,
       const real * q_r,
       const real * aux_l,
       const real * aux_r,
       const solverVariables_t * pFV,
       const std::string& flux_type,
       real * numericalFlux)
   {
       real r_q_l[8];
       real r_q_r[8];
       real r_flux[8];
       real r_flux_l[8];
       real r_flux_r[8];
       real c;
    
       /* initialize scalars rho and e. No need for vectors p and B since they get rotated in next step */
       r_q_l[0] = q_l[0];          /* initialize rho on left */
       r_q_l[4] = q_l[4];          /* initialize e on left */
       r_q_r[0] = q_r[0];          /* ditto rho left */
       r_q_r[4] = q_r[4];          /* ditto e right */
    
       /* local face frame rotation */
       rotate_vector(q_l+1, pFV, r_q_l+1);
       rotate_vector(q_l+5, pFV, r_q_l+5);
       rotate_vector(q_r+1, pFV, r_q_r+1);
       rotate_vector(q_r+5, pFV, r_q_r+5);
    
       /* calculate flux left and right */
       imhd_numerical_flux_local_frame(pC->gas_gamma, r_q_l, r_flux_l);
       imhd_numerical_flux_local_frame(pC->gas_gamma, r_q_r, r_flux_r);
    
       // decide which numerical flux to use
       // switch doesn't work readily with strings in c/c++
       if (flux_type == "Rusanov")
       {
           // rusanov_flux
           c = rusanov_flux(pC, r_q_l, r_q_r, r_flux_l, r_flux_r, r_flux);
       }
    
       else if (flux_type == "HLL")
       {
           // hll flux
           c = hll_flux(pC, r_q_l, r_q_r, r_flux_l, r_flux_r, r_flux);
       }
    
       else if (flux_type == "HLLD")
       {
           // hlld flux
           c = hlld_flux(pC, r_q_l, r_q_r, r_flux_l, r_flux_r, r_flux);
       }
    
       else if (flux_type == "Roe")
       {
           // roe flux
           c = roe_flux(pC, r_q_l, r_q_r, r_flux_l, r_flux_r, r_flux);
       }
    
       else
       {
           WxExcept wxe("Got bad IMHD Numerical Flux Name\n");
           throw wxe;
       }
       
    
    
       
       /* rotate back momentum and mag field */
       antirotate_vector(r_flux+1, pFV, numericalFlux+1);
       antirotate_vector(r_flux+5, pFV, numericalFlux+5);
    
       /* apply calculation of scalars rho and e */
       numericalFlux[0] = r_flux[0];
       numericalFlux[4] = r_flux[4];
    
       return c;
    
   }

   where the specific flux functions are also given in imhd.cc, which can be looked up. Options are Rusanov, HLL, HLLD, and Roe (which requires the addition of a source term). Note that numerical fluxes are written in terms of a local face frame, which requires a dot with a normal or equivalently, a rotation operation, which is what is going on in this function. The procedure is:

   1. Rotate vector componants (may need to copy scalars if the vectors are rotated to different variables) from global frame to local face frame
   2. Apply numerical flux function based on left and right rotated states
   3. Antirotate vector componants back to global frame (copy any scalars if required)

2. Registration - this is at the bottom of wmapplication_imhd_flux.cc

   (src/apps/imhd/wmapplication_imhd_flux.cc)

   #include "lib/wxcreator.h"
    
   // IMHD flux
   REGISTER_CREATOR("iMHDFlux", WmApplication_IMHD_Flux, WmApplication);

3. Add to CMakeLists.txt

   (src/apps/imhd/CMakeLists.txt)

   add_sources(
       imhd.cc
       wmapplication_imhd_flux.cc
   )

4. Add to warpy

   (tools/warpy/apps/imhd/imhd.py)

   class imhd(application):
       '''
       @brief Ideal MHD
       '''
       def __init__(self, name, gamma,
                    variable, numerical_flux_type='Rusanov',floor_density=None, floor_pressure=None,components=None):
           super(imhd, self).__init__(
               name,
               req_attrs=[('Kind', 'iMHDFlux'),
                          ('gamma', gamma),
                          ('NumericalFluxType',numerical_flux_type)],
               opt_attrs=[('MinDensity', floor_density),
                          ('MinPressure',floor_pressure)],            
               variables=[('Fluid', variable, components)])

5. Add warpy function to init.py

   (tools/warpy/apps/__init__.py)

   from .imhd import imhd

6. Add warpy function to CMake

   (tools/warpy/apps/CMakeLists)

   WARPY_FILES(__init__.py
       application.py
       )
    
   add_subdirectory(imhd)

Example 3 Resistive MHD Formulation

Now we would like to add resistive diffusion terms. These appear in Ohm's law as an addition to the ideal term:

\begin{align*} \boldsymbol{E} =& -\boldsymbol{u}\times\boldsymbol{B} + \color{red}{\left(\frac{\delta_{p}}{L}\right)\left(\nu_{p}\tau\right)\eta \boldsymbol{j}} \end{align*}

where

\begin{align*} \boldsymbol{j} =& \left(\frac{\delta_{p}}{L}\right)\nabla\times\boldsymbol{B} \end{align*}

Also

\begin{align*} \eta =& K\frac{A_{e}^{\frac{1}{2}}}{T_{e}^{\frac{3}{2}}} \end{align*}

where for now \( KA_{e}^{\frac{1}{2}}\) is being lumped into the \(\left(\nu_{p}\tau\right)\) and

\begin{align*} T_{e} =& \frac{p}{2\rho} \end{align*}

where assumptions that \(T_{e}=T_{i}\) and singly-charged ions have been made.

Next, to actually incorporate these terms into the equation set we must notice where Ohm's Law fits into the MHD equations. This happens where \(\boldsymbol{E}\) appears, shows up in the energy and induction equation. The energy equation for the the electromagnetic terms

\begin{align*} \frac{\partial e_{em}}{\partial t} - \nabla \cdot\left(\boldsymbol{B}\times\boldsymbol{E}\right) =& 0 \end{align*}

and induction equation

\begin{align*} \frac{\partial \boldsymbol{B}}{\partial t} + \nabla \times \boldsymbol{E} =& 0 \end{align*}

where here we are adding the term \(\color{red}{\boldsymbol{E}=\left(\frac{\delta_{p}}{L}\right)\left(\nu_{p}\tau\right)\eta \boldsymbol{j}}\).

We can write the energy equation as

\begin{align*} \frac{\partial e_{em}}{\partial t} -\frac{\partial}{\partial x} \begin{pmatrix} B_{y}E_{z} - B_{z}E_{y} \end{pmatrix} -\frac{\partial}{\partial y} \begin{pmatrix} B_{z}E_{x} - B_{x}E_{z} \end{pmatrix} -\frac{\partial}{\partial z} \begin{pmatrix} B_{x}E_{y} - B_{y}E_{x} \end{pmatrix} =& 0 \end{align*}

We can write the induction equation in divergence form

\begin{align*} \frac{\partial}{\partial t} \begin{pmatrix} B_{x}\\ B_{y}\\ B_{z} \end{pmatrix} + \begin{pmatrix} \partial_{y}E_{z} - \partial_{z}E_{y}\\ -\partial_{x}E_{z} + \partial_{z}E_{x}\\ \partial_{x}E_{y} - \partial_{y}E_{x} \end{pmatrix} =& 0 \nonumber \\ \frac{\partial}{\partial t} \begin{pmatrix} B_{x}\\ B_{y}\\ B_{z} \end{pmatrix} + \partial_{x} \begin{pmatrix} 0 \\ -E_{z} \\ +E_{y} \end{pmatrix} + \partial_{y} \begin{pmatrix} +E_{z} \\ 0 \\ -E_{x} \end{pmatrix} + \partial_{z} \begin{pmatrix} -E_{y} \\ +E_{x} \\ 0 \end{pmatrix} =& 0 \nonumber \\ \frac{\partial}{\partial t} \begin{pmatrix} B_{x}\\ B_{y}\\ B_{z} \end{pmatrix} + \begin{pmatrix} 0 & +E_{z} & -E_{y} \\ -E_{z} & 0 & +E_{x} \\ +E_{y} & -E_{x} & 0 \end{pmatrix} \begin{pmatrix} \frac{\partial}{\partial x} \\ \frac{\partial}{\partial y} \\ \frac{\partial}{\partial z} \end{pmatrix} =& 0 \end{align*}

Altogether this is written for resitive MHD

\begin{align*} \frac{\partial}{\partial t} \begin{pmatrix} \rho \\ \rho u \\ \rho v \\ \rho w \\ e_{t} \\ B_{x} \\ B_{y} \\ B_{z} \end{pmatrix} + \frac{\partial}{\partial x} \begin{pmatrix} \rho u \\ \rho u u - (B_{x} B_{x} - \frac{1}{2}B^{2}) + p \\ \rho u v - B_{x} B_{y} \\ \rho u w - B_{x} B_{z} \\ \left(e_{t} + p + \frac{1}{2}B^{2}\right)u - \left(\boldsymbol{B}\cdot\boldsymbol{u}\right)B_{x} \color{red}{- \left(\frac{\delta_{p}}{L}\right)\left(\nu_{p}\tau\right)\eta \left(B_{y}j_{z} - B_{z}j_{y}\right)}\\ 0 \color{red}{+ 0} \\ u B_{y} - B_{x} v \color{red}{-\left(\frac{\delta_{p}}{L}\right)\left(\nu_{p}\tau\right)\eta j_{z}}\\ u B_{z} - B_{x} w \color{red}{+\left(\frac{\delta_{p}}{L}\right)\left(\nu_{p}\tau\right)\eta j_{y}}\\ \end{pmatrix} + \frac{\partial}{\partial y} \begin{pmatrix} \rho v \\ \rho v u - B_{y} B_{x} \\ \rho v v - (B_{y} B_{y} - \frac{1}{2}B^{2}) + p \\ \rho v w - B_{y} B_{z} \\ \left(e_{t} + p + \frac{1}{2}B^{2}\right)v - \left(\boldsymbol{B}\cdot\boldsymbol{u}\right)B_{y} \color{red}{- \left(\frac{\delta_{p}}{L}\right)\left(\nu_{p}\tau\right)\eta \left(B_{z}j_{x} - B_{x}j_{z}\right)}\\ v B_{x} - B_{y} u \color{red}{+\left(\frac{\delta_{p}}{L}\right)\left(\nu_{p}\tau\right)\eta j_{z}}\\ 0 \color{red}{+ 0}\\ v B_{z} - B_{y} w \color{red}{-\left(\frac{\delta_{p}}{L}\right)\left(\nu_{p}\tau\right)\eta j_{x}}\\ \end{pmatrix} + \frac{\partial}{\partial z} \begin{pmatrix} \rho w \\ \rho w u - B_{z} B_{x} \\ \rho w v - B_{z} B_{y} \\ \rho w w - (B_{z} B_{z} - \frac{1}{2}B^{2}) + p \\ \left(e_{t} + p + \frac{1}{2}B^{2}\right)w - \left(\boldsymbol{B}\cdot\boldsymbol{u}\right)B_{z} \color{red}{- \left(\frac{\delta_{p}}{L}\right)\left(\nu_{p}\tau\right)\eta \left(B_{x}j_{y} - B_{y}j_{x}\right)}\\ w B_{x} - B_{z} u \color{red}{-\left(\frac{\delta_{p}}{L}\right)\left(\nu_{p}\tau\right)\eta j_{y}}\\ w B_{y} - B_{z} v \color{red}{+\left(\frac{\delta_{p}}{L}\right)\left(\nu_{p}\tau\right)\eta j_{x}}\\ 0 \color{red}{+ 0} \end{pmatrix} =& \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} \end{align*}

or

\begin{align*} \begin{pmatrix} \rho \\ \rho u \\ \rho v \\ \rho w \\ e_{t} \\ B_{x} \\ B_{y} \\ B_{z} \end{pmatrix} + \begin{pmatrix} \rho u & \rho v & \rho w \\ \rho u u - (B_{x} B_{x} - \frac{1}{2}B^{2}) + p & \rho v u - B_{y} B_{x} & \rho w u - B_{z} B_{x} \\ \rho u v - B_{x} B_{y} & \rho v v - (B_{y} B_{y} - \frac{1}{2}B^{2}) + p & \rho w v - B_{z} B_{y} \\ \rho u w - B_{x} B_{z} & \rho v w - B_{y} B_{z} & \rho w w - (B_{z} B_{z} - \frac{1}{2}B^{2}) + p \\ \left(e_{t} + p + \frac{1}{2}B^{2}\right)u - \left(\boldsymbol{B}\cdot\boldsymbol{u}\right)B_{x} \color{red}{- \left(\frac{\delta_{p}}{L}\right)\left(\nu_{p}\tau\right)\eta \left(B_{y}j_{z} - B_{z}j_{y}\right)} & \left(e_{t} + p + \frac{1}{2}B^{2}\right)v - \left(\boldsymbol{B}\cdot\boldsymbol{u}\right)B_{y} \color{red}{- \left(\frac{\delta_{p}}{L}\right)\left(\nu_{p}\tau\right)\eta \left(B_{z}j_{x} - B_{x}j_{z}\right)} & \left(e_{t} + p + \frac{1}{2}B^{2}\right)w - \left(\boldsymbol{B}\cdot\boldsymbol{u}\right)B_{z} \color{red}{- \left(\frac{\delta_{p}}{L}\right)\left(\nu_{p}\tau\right)\eta \left(B_{x}j_{y} - B_{y}j_{x}\right)} \\ 0 \color{red}{+ 0} & v B_{x} - B_{y} u \color{red}{+\left(\frac{\delta_{p}}{L}\right)\left(\nu_{p}\tau\right)\eta j_{z}} & w B_{x} - B_{z} u \color{red}{-\left(\frac{\delta_{p}}{L}\right)\left(\nu_{p}\tau\right)\eta j_{y}} \\ u B_{y} - B_{x} v \color{red}{-\left(\frac{\delta_{p}}{L}\right)\left(\nu_{p}\tau\right)\eta j_{z}} & 0 \color{red}{+ 0} & w B_{y} - B_{z} v \color{red}{+\left(\frac{\delta_{p}}{L}\right)\left(\nu_{p}\tau\right)\eta j_{x}} \\ u B_{z} - B_{x} w \color{red}{+\left(\frac{\delta_{p}}{L}\right)\left(\nu_{p}\tau\right)\eta j_{y}} & v B_{z} - B_{y} w \color{red}{-\left(\frac{\delta_{p}}{L}\right)\left(\nu_{p}\tau\right)\eta j_{x}} & 0 \color{red}{+ 0} \end{pmatrix} \begin{pmatrix} \frac{\partial}{\partial x} \\ \frac{\partial}{\partial y} \\ \frac{\partial}{\partial z} \end{pmatrix} =& \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} \end{align*}

NOTE: When we write the flux tensor in WARPXM we write its transpose (see code below).

Now that we have the equations, we can add on these terms two ways:

1. Write a single resistive mhd module that writes the entire equation set
2. Make a module with only resistive diffusion and make a simulation solve both ideal mhd + resistive difffusion modules

Example 3.1 Resistive MHD in 1 module

In this example we write the resistive mhd module (wmapplication_rmhd.h/cc) using helpers rmhd.h/cc to write the entire resistive mhd equation set. It will consist of writing internal and numerical fluxes specifically for the resistive diffusion terms and adding them to the ideal mhd fluxes within the module.

1. Internal Flux, Numerical Flux

   (src/apps/rmhd/wmapplication_rmhd.h)

   #ifndef WMAPPLICATION_RMHD_H
   #define WMAPPLICATION_RMHD_H
    
   // WarpM includes
   #include "apps/wmapplication.h"
    
   class WmApplication_RMHD : public WmApplication
   {
   public:
       WmApplication_RMHD();
    
       ~WmApplication_RMHD();
    
       void setup(const WxCryptSet& wxc) override;
    
       const std::vector<int>& getInputVariableIndexes(
           int flag) const override
       {
           return _inputVariables;
       }
       const std::vector<int>& getAuxiliaryVariableIndexes(
           int flag) const override
       {
           return _auxVariables;
       }
       const std::vector<int>& getOutputVariableIndexes(
           int flag) const override
       {
           return _outputVariables;
       }
       const std::vector<int>& getCrossVariableIndexes(
           int flag) const override
       {
           return _inputVariables;
       }
    
       // redefinition of fluxes from parent
       real numerical_flux(const real* q_l, const real* q_r, const real* aux_l,
           const real* aux_r, const solverVariables_t* pFV,
           real* numericalFlux) const override;
       real internal_flux(const real* q, const real* aux,
           const solverVariables_t* pSV,
           std::vector<std::vector<real>>& internalFlux) const override;
    
   protected:
       real _gas_gamma;
       real _skin_depth_norm;
       real _nu_p_norm;
    
       std::string _imhd_numerical_flux_type; // this is the type of numerical flux you want to use (rusanov, roe, hll, hlld, etc.)        
    
       std::vector<int> _inputVariables;
       std::vector<int> _auxVariables;
       std::vector<int> _outputVariables;
    
    
       real _min_density_floor;
       real _min_pressure_floor;
    
    
   private:
       WmApplication_RMHD& operator=(const WmApplication_RMHD& var);
       WmApplication_RMHD(const WmApplication_RMHD& var);
   };
    
   #endif // WMAPPLICATION_RMHD

   (src/apps/rmhd/wmapplication_rmhd.cc)
   

   #include "wmapplication_rmhd.h"
   #include "rmhd.h"
    
   // STD libraries
   #include <iostream>
    
   // Wm libraries
   #include "lib/wmudgfunctions.h"
    
   WmApplication_RMHD::WmApplication_RMHD()
   {
       _allowedFlags.push_back(WMAPPLICATIONFLAG_NUMERICALFLUX);
       _allowedFlags.push_back(WMAPPLICATIONFLAG_INTERNALFLUX);
   }
    
   WmApplication_RMHD::~WmApplication_RMHD()
   {
    
   }
    
   void
   WmApplication_RMHD::setup(const WxCryptSet& wxc)
   {
       WmApplication::setup(wxc);
    
       WmApplication_Variable field(wxc, "Fluid", 8);
       WmApplication_Variable gradient_field(wxc, "FluidGradient", 24);
    
    
       _gas_gamma = wxc.get<real>(GAS_GAMMA);
       _skin_depth_norm = wxc.get<real>(SKIN_DEPTH_NORM);
       _nu_p_norm = wxc.get<real>(NU_P_NORM);
    
    
       _min_density_floor = 0.0;
       if (wxc.has("MinDensity"))
       {
           _min_density_floor = wxc.get<real>("MinDensity");
       }
       
       _min_pressure_floor = 0.0;
       if (wxc.has("MinPressure"))
       {
           _min_pressure_floor = wxc.get<real>("MinPressure");
       }
       
    
       _imhd_numerical_flux_type = "Rusanov";
       if (wxc.has("ImhdNumericalFluxType"))
       {
           _imhd_numerical_flux_type = wxc.get<std::string>("ImhdNumericalFluxType");
       }
       
    
       _inputVariables = field.getIndexes();
       _outputVariables = field.getIndexes();
       _auxVariables = gradient_field.getIndexes();
   }
    
   real
   WmApplication_RMHD::internal_flux(const real * q,
                              const real * aux,
                              const solverVariables_t * pSV,
                              // real *internalFlux[8]
                              std::vector< std::vector<real> >& internalFlux
       ) const
   {
    
    
       // q = [rho, px, py, pz, e, Bx, By, Bz]
       // aux = [rho_x, rho_y, rho_z,
       // px_x, px_y, px_z, py_x, py_y, py_z, pz_x, pz_y, pz_z,
       // e_x, e_y, e_z
       // Bx_x, Bx_y, Bx_z, By_x, By_y, By_z, Bz_x, Bz_y, Bz_z]
    
       // Initialize constants
       constants_rmhd_t pC;
       populate_rmhd_t(&pC, _gas_gamma, _skin_depth_norm, _nu_p_norm);
    
       // ensure positivity of density and pressure
       real dens_q[8]; 
       real mod_q[8];
       ensure_floor_density(q,8,_min_density_floor,dens_q);
       ensure_floor_pressure(dens_q,8,_gas_gamma,_min_pressure_floor,mod_q);
       
    
       // Advective (IDEAL MHD) portion
       // output vector allocation is [num_dims] x [_num_output_Q_max]
       std::vector<real> N_output_dimension_vector(8, 0.);
       std::vector< std::vector<real> > internal_advective_flux(pSV->num_dims, N_output_dimension_vector);
       // calculate internal flux from ideal mhd, put into internal_advective_flux
       imhd_internal_flux(&(pC.imhd), pSV, q, internal_advective_flux);
    
    
    
       
       // Diffusive (Resistive Terms) portion 
       // output vector allocation is [num_dims] x [_num_output_Q_max]
       N_output_dimension_vector.assign(8, 0.);
       std::vector< std::vector<real> > internal_diffusive_flux(pSV->num_dims, N_output_dimension_vector);
    
       // calculate resistivity
       real eta = get_resistivity(&pC, q);
    
    
       // calculate J
       // const real * dp[3] = {aux+3, aux+6, aux+9}; // dpx, dpy, dpz not needed
       // const real * de[3] = {aux+12}; // de not needed
       const real * dB[3] = {aux+15, aux+18, aux+21};
       real J[3];
       real abs_j_coeff = current_density(&pC,dB,J);
       /* calculate internal diffusion flux */
       real abs_diff_coeff = rmhd_internal_diffusion_flux(&pC,pSV,mod_q,J,eta,internal_diffusive_flux);
    
    
       
       /* Add the 2 fluxes together */
       for(int i=0;i<pSV->num_dims;i++)
       {
           for(int j=0;j<8;j++)
           {
               internalFlux[i][j] = internal_advective_flux[i][j] + internal_diffusive_flux[i][j];
               /* printf("advective_flux[%d][%d]=%f\n",i,j,internal_advective_flux[i][j]); */
           }
    
       }
       
       return set_time_step(pSV, abs_j_coeff, abs_diff_coeff, eta);
       
       /* return pSV->cfl_max * pSV->cfl_max * pSV->dxn * pSV->dxn / (abs_j_coeff * abs_diff_coeff * fabs(eta)); */
    
   }
    
   real
   WmApplication_RMHD::numerical_flux(const real * q_l,
                                           const real * q_r,
                                           const real * aux_l,
                                           const real * aux_r,
                                           const solverVariables_t * pFV,
                                           real * numericalFlux) const
   {
       // q = [rho, px, py, pz, e, Bx, By, Bz]
       // aux = [rho_x, rho_y, rho_z,
       // px_x, px_y, px_z, py_x, py_y, py_z, pz_x, pz_y, pz_z,
       // e_x, e_y, e_z
       // Bx_x, Bx_y, Bx_z, By_x, By_y, By_z, Bz_x, Bz_y, Bz_z]
    
    
    
       /* advective part */
       // Define the constants
       constants_imhd_t pCi;
       populate_constants_imhd_t(&pCi, _gas_gamma);    
       // pCi.gas_gamma = _gas_gamma;
       // ensure positivity of density and pressure
       real dens_q_l[8], dens_q_r[8];
       real mod_q_l[8], mod_q_r[8];    
       ensure_floor_density(q_l,8,_min_density_floor,dens_q_l);
       ensure_floor_density(q_r,8,_min_density_floor,dens_q_r);    
       ensure_floor_pressure(dens_q_l,8,_gas_gamma,_min_pressure_floor,mod_q_l);
       ensure_floor_pressure(dens_q_r,8,_gas_gamma,_min_pressure_floor,mod_q_r);
       /* allocate numerical advective flux */ 
       real numerical_advective_flux[8];
       for(int i=0;i<8;i++) { numerical_advective_flux[i] = 0; }
       real c = imhd_numerical_flux(&pCi, q_l, q_r, aux_l, aux_r, pFV, _imhd_numerical_flux_type, numerical_advective_flux);
    
    
       /* diffusive part */
       const real * dB_l[3] = {aux_l+15, aux_l+18, aux_l+21};
       const real * dB_r[3] = {aux_r+15, aux_r+18, aux_r+21};
    
       real J_l[3];
       real J_r[3];
    
       // Initialize constants
       constants_rmhd_t pC;
       populate_rmhd_t(&pC, _gas_gamma, _skin_depth_norm, _nu_p_norm);
    
    
       /* allocate numerical diffusive flux */
       real numerical_diffusive_flux[8];
       for(int i=0;i<8;i++) { numerical_diffusive_flux[i] = 0; }
    
    
       // calcualte resistivity
       real eta_l = get_resistivity(&pC, mod_q_l);
       real eta_r = get_resistivity(&pC, mod_q_r);
    
    
       /* calculate J */
       real abs_j_coeff_l = current_density(&pC,dB_l,J_l); // careful:  we might have modded density and pressure but not their gradients...
       real abs_j_coeff_r = current_density(&pC,dB_r,J_r);
       
       
    
       /* calculate numerical diffusion flux */
       real abs_diff_coeff = rmhd_numerical_diffusion_flux(&pC,pFV,mod_q_l,mod_q_r,J_l,J_r,eta_l,eta_r,numerical_diffusive_flux);
    
       /* add the 2 fluxes together */
       for(int i=0;i<8;i++)
       {
           numericalFlux[i] =  numerical_advective_flux[i] +  numerical_diffusive_flux[i];
       }
    
    
    
       return set_time_step(pFV,fmax(abs_j_coeff_l,abs_j_coeff_r), abs_diff_coeff, fmax(eta_l,eta_r));
       
       /* return pFV->cfl_max * pFV->cfl_max * pFV->dxn * pFV->dxn / (max(abs_j_coeff_l,abs_j_coeff_r) * abs_diff_coeff * fabs(eta)); */
    
    
    
   }
    
    
   #include "lib/wxcreator.h"
    
   // RMHD
   REGISTER_CREATOR("rMHD", WmApplication_RMHD, WmApplication);

   the resistive diffusion flux functions are written in rmhd.cc.  For internal flux of diffusive terms we have
   (src/apps/rmhd/rmhd.cc)
   

   real rmhd_internal_diffusion_flux(
       const constants_rmhd_t * pC,
       const solverVariables_t * pSV,
       const real * q,
       const real * J,
       const real eta,
       // real * internal_diffusion_flux[8]
       std::vector< std::vector<real> >& internal_diffusion_flux
       )
   {
       real Bx = q[5], By = q[6], Bz = q[7];
       real Jx = J[0], Jy = J[1], Jz = J[2];
       // real norm_factor_energy = pC->nu_p_norm * pC->imhd.omega_c_norm * pC->imhd.skin_depth_norm * pC->imhd.skin_depth_norm;
       real norm_factor_energy = pC->nu_p_norm * pC->skin_depth_norm;
       // real norm_factor_induction = pC->nu_p_norm / pC->imhd.omega_c_norm; /* setting omega_c tau to 0 is gonna give problems */
       real norm_factor_induction = pC->nu_p_norm * pC->skin_depth_norm; // actually same factor
    
    
       
       /* clear the tensor */
       for (int i=0;i<pSV->num_dims;i++)
       {
           for (int j=0;j<8;j++)
           {
               internal_diffusion_flux[i][j] = 0.;
               
           }
           
       }
       /* add to specific components */
       /* solving dq_j/dt + dF_ij/dx_i = 0  */
    
       /* dx terms */
       internal_diffusion_flux[0][4] += -norm_factor_energy * eta * (By*Jz - Bz*Jy); /* energy */
       internal_diffusion_flux[0][5] += 0; /* Bx */
       internal_diffusion_flux[0][6] += -norm_factor_induction * eta * Jz; /* By */
       internal_diffusion_flux[0][7] += +norm_factor_induction * eta * Jy; /* Bz */
    
       /* printf("got to line 107 of internal diffusion flux\n"); */
       if (pSV->num_dims > 1)
       {
           /* dy terms */
           internal_diffusion_flux[1][4] += -norm_factor_energy * eta * (Bz*Jx - Bx*Jz); /* energy */;
           internal_diffusion_flux[1][5] += +norm_factor_induction * eta * Jz; /* Bx */;
           internal_diffusion_flux[1][6] += 0; /* By */
           internal_diffusion_flux[1][7] += -norm_factor_induction * eta * Jx; /* Bz */;
    
       }
    
       if (pSV->num_dims > 2)
       {
           /* dz terms */
           internal_diffusion_flux[2][4] += -norm_factor_energy * eta * (Bx*Jy - By*Jx); /* energy */;
           internal_diffusion_flux[2][5] += -norm_factor_induction * eta * Jy; /* Bx */;
           internal_diffusion_flux[2][6] += +norm_factor_induction * eta * Jx; /* By */
           internal_diffusion_flux[2][7] += 0; /* Bz */;
    
       }
    
    
       /* printf("got to line 123 of internal diffusion flux\n"); */
       return fmax(fabs(norm_factor_energy), fabs(norm_factor_induction));
   }

   These terms also have a numerical flux defined, according to the Local Discontinuous Galerkin formulation
   
     (src/apps/rmhd/rmhd.cc)
   

   real rmhd_numerical_diffusion_flux(
           const constants_rmhd_t * pC,
           const solverVariables_t * pSV,
           const real * q_l,
           const real * q_r,
           const real * J_l,
           const real * J_r,
           const real eta_l,
           const real eta_r,
           real * numerical_diffusion_flux
       )
   {
       // iman's version LDG
       // The numerical flux will be given by
       // F = (1/2-beta_times_common_n_dot_n) * (F_l dot n) + (1/2+beta_times_common_n_dot_n) * (F_r dot n) - eta * (q_l - q_r)
       // where l represents "inside" or "lambda" element, and r reppresents "outside" or gamma element
       // let alpha_l = (1/2-beta_times_common_n_dot_n)
       // let alpha_r = (1/2+beta_times_common_n_dot_n)
    
    
       const real * norm = pSV->R[0]; /* normal vector */
       
    
       const real alpha_l = 0.5 - pSV->beta_times_common_n_dot_n;  
       const real alpha_r = 0.5 + pSV->beta_times_common_n_dot_n;
    
       const real eta = pSV->eta;
    
       /* calculate internal diffusion flux */
       /* memory allocation */
    
       // output vector allocation is [num_dims] x [_num_output_Q_max]
       std::vector<real> N_output_dimension_vector(8, 0.);
       std::vector< std::vector<real> > internal_diffusion_flux(pSV->num_dims, N_output_dimension_vector);
           
    
    
       
       /* Calculate left state */
       real abs_diff_coeff_l = rmhd_internal_diffusion_flux(pC,pSV,q_l,J_l,eta_l,internal_diffusion_flux);
       // Initialize the numerical flux vector with the left state
    
       for(int i = 0; i < 8; i++){
           numerical_diffusion_flux[i] = 0.;
           for(int j = 0; j < pSV->num_dims; j++)
           {
               numerical_diffusion_flux[i] += alpha_l * internal_diffusion_flux[j][i] * norm[j];
           }
       }
    
    
       /* Right state */
       N_output_dimension_vector.assign(8, 0.);
       internal_diffusion_flux.assign(pSV->num_dims, N_output_dimension_vector);
       
       real abs_diff_coeff_r = rmhd_internal_diffusion_flux(pC,pSV,q_r,J_r,eta_r,internal_diffusion_flux);
       // Initialize the numerical flux vector with the left state
       for(int i = 0; i < 8; i++)
       {
           for(int j = 0; j < pSV->num_dims; j++){
               numerical_diffusion_flux[i] += alpha_r * internal_diffusion_flux[j][i] * norm[j];
           }
       }
    
    
       // Next we do the final penalty (don't apply to density)
       for(int i = 0; i < 8; i++)
       {
           // if( i != 0)
           {
               numerical_diffusion_flux[i] -= eta * (q_l[i] - q_r[i]);
           }
    
       }
    
    
       // This probably needs to be a bit more restrictive
       return fmax(abs_diff_coeff_l,abs_diff_coeff_r);
   }

2. Registration - this is at the bottom of wmapplication_rmhd.cc (src/apps/imhd/wmapplication_rmhd.cc)

   #include "lib/wxcreator.h"
    
   // RMHD
   REGISTER_CREATOR("rMHD", WmApplication_RMHD, WmApplication);

3. Add to CMakeLists.txt

   (src/apps/rmhd/CMakeLists.txt)

   add_sources(
       rmhd.cc
       wmapplication_rmhd.cc
   )

4. Add to warpy

   (tools/warpy/apps/rmhd/rmhd.py)

   from ..application import application
    
   class rmhd(application):
       '''
       @brief Resistive MHD
       '''
       def __init__(self, name, skin_depth_norm,nu_p_norm, gamma,
                    variable, gradient_variable, imhd_numerical_flux_type=None, floor_density=None, floor_pressure=None, components=None, gradient_variable_components=None):
           super(rmhd, self).__init__(
               name,
               req_attrs=[('Kind', 'rMHD'),
                          ('skin_depth_norm', skin_depth_norm),
                          ('nu_p_norm', nu_p_norm),
                          ('gamma', gamma),
                          ('ImhdNumericalFluxType',imhd_numerical_flux_type)],
               opt_attrs=[('MinDensity', floor_density),
                          ('MinPressure',floor_pressure)],            
               variables=[('Fluid', variable, components),
                          ('FluidGradient', gradient_variable, gradient_variable_components)])

5. Add warpy function to init

   (tools/warpy/apps/__init__.py)

   from . import rmhd

   (tools/warpy/apps/rmhd/__init__.py)

   from .rmhd import rmhd

6. Add warpy function to CMake

   (tools/warpy/apps/CMakeLists.txt)

   WARPY_FILES(__init__.py
       application.py
       )
   add_subdirectory(rmhd)

   (tools/warpy/apps/rmhd/CMakeLists.txt)

   WARPY_FILES(
       __init__.py
       rmhd.py
   )

Example 3.1.1 Initial Condition

In order to actually run a problem you'll want to set some intial condition and boundary conditions (unless periodic, which you set in the mesh). There is a library of initial conditions, many of which are for single variables, with various options, such as "set_to", "Heaviside", and "polynomial", while others are for equation sets, such as mhd. Here we will develop wmicfunction_mhd_custom. The steps are very similar to writing an application, but instead of writing functions for fluxes and sources, we write an "applyFunction", which for a given coordinate position \((x,y,z)\).

1. Writing a function class for initial condition. Here I make a custom file, defining

   • Constructor
   • Destructor
   • setup
   • applyFunction
   • an required variables

   (src/initialConditions/wmicfunction_mhd_custom.h)

   #ifndef WMICFUNCTION_MHD_CUSTOM_H
   #define WMICFUNCTION_MHD_CUSTOM_H
    
   // std includes
   #include <vector>
   #include <string>
    
   // Wm includes
   #include "initialConditions/wmicfunction.h"
    
   class WmICFunction_Mhd_Custom: public WmICFunction
   {
   public:
    
       WmICFunction_Mhd_Custom();
    
       ~WmICFunction_Mhd_Custom();
    
       void setup(const WxCryptSet& wxc);
    
       void applyFunction(const real * position, real * value) const;
    
   private:
    
       real _gas_gamma;
       real _skin_depth_param;
    
   };
    
    
   #endif // WMICFUNCTION_MHD_CUSTOM_H

   Set the number of arguments, initialize variables from cryptset in setup, and define function in applyFunction. (src/initialConditions/wmicfunction_mhd_custom.cc)

   #include "wmicfunction_mhd_custom.h"
    
   // STL includes
   #include <cmath>
    
    
   // Wm includes
   // #include "lib/wmpearsonIVfunctions.h"
    
   WmICFunction_Mhd_Custom::WmICFunction_Mhd_Custom()
   {
       _numArgs = 8;
   }
    
   WmICFunction_Mhd_Custom::~WmICFunction_Mhd_Custom()
   {
    
   }
    
   void
   WmICFunction_Mhd_Custom::setup(const WxCryptSet& wxc)
   {
       WmICFunction::setup(wxc);
       
       _gas_gamma = wxc.get<real>(GAS_GAMMA);
       _skin_depth_param = wxc.get<real>(SKIN_DEPTH_NORM);     
    
   }
    
   void
   WmICFunction_Mhd_Custom::applyFunction(const real * position, real * qc) const
   {
    
       const real x = position[0];
       const real y = position[1];
       const real z = position[2];
    
       // Gaussian
       real min = 0.0;
       real max = 1.0;
       std::vector<real> average = std::vector<real>(3,0.);
       int numDims = average.size();
       real standardDeviation = 0.125;
           
       real C0 = _min;
       real C1 = (max - min);
       real C3 = 0.;
       for(int i = 0; i < numDims; i++){
           C3 += (position[i] - average[i]) * (position[i] - average[i]);
       }
       real C2 = - C3 / (2. * standardDeviation * standardDeviation);
    
       real gaussian_value = C0 + C1 * exp(C2);
    
    
       // constant
       real constant_value1 = 1.0;
       real constant_value2 = 1.0;
    
       // Sine wave
       // real sinewave = std::cos(4*wxm::pi*x);
       real sinewave = std::sin(2*wxm::pi*x);  
    
       real rho = constant_value1;
       real vx = 0.0;
       real vy = 0.0;
       real vz = 0.0;
       real P = constant_value2;
       real Bx = 0.0;
       real By = constant_value1;//gaussian_value; //sinewave;// 
       real Bz = 0.0;
    
    
       // Convert to the output conserved variables
       real gamma1 = _gas_gamma - 1.0;
       real V2 = vx*vx + vy*vy + vz*vz;
       real B2 = Bx*Bx + By*By + Bz*Bz;
       real E = P/gamma1 + 0.5*rho*V2 + 0.5*B2;
       
       qc[0] = rho;
       qc[1] = rho*vx;
       qc[2] = rho*vy;
       qc[3] = rho*vz;
       qc[4] = E;
       qc[5] = Bx;
       qc[6] = By;
       qc[7] = Bz;
    
   }
    
   // Register
   #include "lib/wxcreator.h"
   REGISTER_CREATOR("custom_mhd", WmICFunction_Mhd_Custom, WmICFunction);

2. Registration - The last lines in the wmicfunction_mhd_custom.cc file registers the scheme:

   (src/initialConditions/wmicfunction_mhd_custom.cc)

   // Register
   #include "lib/wxcreator.h"
   REGISTER_CREATOR("custom_mhd", WmICFunction_Mhd_Custom, WmICFunction);

3. Add to CMakeLists.txt

   (src/apps/imhd/CMakeLists.txt)

   add_sources(
       wmicfunction.cc
       wmicfunction_mhd_custom.cc
   )

4. Add to warpy

   (tools/warpy/functions/custom_mhd.py)

   from .. import generator as gen
    
   class custom_mhd(gen.warpy_obj):
       '''
       '''
    
       def __init__(self, name, skin_depth_norm, gamma):
           '''
           '''
           super(custom_mhd, self).__init__(name,
                                            req_attrs=[('Type','function'),
                                                       ('Kind', 'custom_mhd'),
                                                       ('skin_depth_norm', skin_depth_norm),
                                                       ('gamma', gamma)])

5. Add warpy function to init

   (tools/warpy/functions/__init__.py)

   from .applicator import applicator
   from .custom_mhd import custom_mhd

6. Add warpy function to CMake

   (tools/warpy/functions/CMakeLists.txt)

   WARPY_FILES(__init__.py
       applicator.py
       custom_mhd.py
       )

Example 3.1.2 Boundary Condition

Boundary conditions are just another application, specific to a particular to a particular set of variables. Often they make sense physically for a specific equation set, e.g. no slip wall for Euler or Ideal MHD. One can also apply Dirichlet or Neumann conditions to particular variables. Here, I show a no-slip wall condition used by resistive mhd.

1. Write boundary condition

   (src/apps/rmhd/wmapplication_rmhd_bc_noslipwall.h)

   #ifndef WMAPPLICATION_RMHD_BC_NOSLIPWALL_H
   #define WMAPPLICATION_RMHD_BC_NOSLIPWALL_H
    
   // WarpM includes
   #include "apps/wmapplication.h"
    
   class WmApplication_Rmhd_BC_NoslipWall : public WmApplication
   {
   public:
       WmApplication_Rmhd_BC_NoslipWall();
    
       ~WmApplication_Rmhd_BC_NoslipWall();
    
       void setup(const WxCryptSet& wxc) override;
    
       const std::vector<int>& getInputVariableIndexes(int flag) const override
       {
           return _inputVariables;
       }
    
       const std::vector<int>& getOutputVariableIndexes(int flag) const override
       {
           return _outputVariables;
       }
    
       // redefinition of bc
       void bc_q(const real* q_in, const solverVariables_t* pFV,
           real* q_out) const override;
    
   protected:
       real _gas_gamma;
       real _skin_depth_norm;
       real _nu_p_norm;
    
       std::vector<real> _wallVelocity;
       std::vector<real> _wallMagneticField;
    
       std::vector<int> _inputVariables;
       std::vector<int> _outputVariables;
    
   private:
       WmApplication_Rmhd_BC_NoslipWall& operator=(
           const WmApplication_Rmhd_BC_NoslipWall& var);
       WmApplication_Rmhd_BC_NoslipWall(
           const WmApplication_Rmhd_BC_NoslipWall& var);
   };
    
   #endif // WMAPPLICATION_RMHD_BC_NOSLIPWALL_H

   (src/apps/rmhd/wmapplication_rmhd_bc_noslipwall.cc)
   

   #include "wmapplication_rmhd_bc_noslipwall.h"
   #include "rmhd.h"
    
   // STD libraries
   #include <iostream>
    
   // Wm includes
   #include "lib/wmudgfunctions.h"
    
   WmApplication_Rmhd_BC_NoslipWall::WmApplication_Rmhd_BC_NoslipWall()
   {
       _allowedFlags.push_back(WMAPPLICATIONFLAG_BC_Q);
   }
    
   WmApplication_Rmhd_BC_NoslipWall::~WmApplication_Rmhd_BC_NoslipWall()
   {
    
   }
    
   void
   WmApplication_Rmhd_BC_NoslipWall::setup(const WxCryptSet& wxc)
   {
    
       WmApplication::setup(wxc);
    
       WmApplication_Variable fluid(wxc, "Fluid", 8);
    
       _gas_gamma = wxc.get<real>(GAS_GAMMA);
       _skin_depth_norm = wxc.get<real>(SKIN_DEPTH_NORM);
       _nu_p_norm = wxc.get<real>(NU_P_NORM);
    
       
       _wallVelocity.resize(3,0.);
       if(wxc.has("WallVelocity")){
           _wallVelocity = wxc.getVec<real>("WallVelocity");
       }
       _wallMagneticField.resize(3,0.);
       if(wxc.has("WallMagneticField")){
           _wallMagneticField = wxc.getVec<real>("WallMagneticField");
       }
    
       _inputVariables = fluid.getIndexes();
       _outputVariables = _inputVariables;
   }
    
   void
   WmApplication_Rmhd_BC_NoslipWall::bc_q(const real * q_in,
           /* const real * restrict aux, */
           const solverVariables_t * pFV,
           real * q_out) const
   {
       /* calcualte numerical diffusive flux from rmhd */
       // Define the constants
       constants_rmhd_t pC;
       populate_rmhd_t(&pC, _gas_gamma, _skin_depth_norm, _nu_p_norm);
    
    
       real wall_vx = _wallVelocity[0], wall_vy = _wallVelocity[1], wall_vz = _wallVelocity[2];
       real wall_bx = _wallMagneticField[0], wall_by = _wallMagneticField[1], wall_bz = _wallMagneticField[2];
    
       // Apply boundary conditions
       noslipWallCondition_MHD(&pC, pFV, q_in, wall_vx, wall_vy, wall_vz, wall_bx, wall_by, wall_bz, q_out);
    
   }
    
    
   // Register the class with wxcreator
   #include "lib/wxcreator.h"
   REGISTER_CREATOR("rMHDNoSlipWall", WmApplication_Rmhd_BC_NoslipWall, WmApplication);

   where noslipWallCondition_MHD is defined with helper in rmhd.cc
   
   (src/apps/rmhd/rmhd.cc)
   

   void
   noslipWallCondition_MHD(
           const constants_rmhd_t * pC,
           const solverVariables_t * pSV,
           const real * QC,
           const real vx_w,
           const real vy_w,
           const real vz_w,
           const real bx_w,
           const real by_w,
           const real bz_w,
           real * QC_w)
   {
    
       // Allocate space for the nonconserved variables
       real QNC[8];
    
       // Calculate nonconserved variables
       // CtoNC_MHD(&(pC->imhd),QC,QNC);
       CtoNC_MHD(&(pC->imhd),QC,QNC);  
    
       // We are assuming that the pressure and the density of the fluid in the wall is the same
    
       // FORM 0: This is the simple form - it works, but I can do better
       // The velocity in the wall is assumed to be the velocity of the wall
      QNC[1] = vx_w;
      QNC[2] = vy_w;
      QNC[3] = vz_w;
      QNC[5] = bx_w;
      QNC[6] = by_w;
      QNC[7] = bz_w;
    
       // // FORM 1: Slightly more complex, but it is derived from the kinetic expression when reflecting the distribution across the wall
       // QNC[1] = 2. * vx_w - QNC[1];
       // QNC[2] = 2. * vy_w - QNC[2];
       // QNC[3] = 2. * vz_w - QNC[3];
       // QNC[5] = 2. * bx_w - QNC[5];
       // QNC[6] = 2. * by_w - QNC[6];
       // QNC[7] = 2. * bz_w - QNC[7];
    
    
    
       // Convert the nonconserved variables into conserved variables
       // NCtoC_MHD(&(pC->imhd),QNC,QC_w);
      NCtoC_MHD(&(pC->imhd),QNC,QC_w); 
       /* printf("rho= %f, px = %f, py = %f, pz = %f, e = %f, bx = %f, by = %f, bz = %f\n",QC_w[0],QC_w[1],QC_w[2],QC_w[3],QC_w[4],QC_w[5],QC_w[6],QC_w[7]); */
    
      
   }

2. Registration - this is at the bottom of wmapplication_rmhd_bc_noslipwall.cc

   (src/apps/rmhd/wmapplication_rmhd_noslipwall.cc)

   // Register the class with wxcreator
   #include "lib/wxcreator.h"
   REGISTER_CREATOR("rMHDNoSlipWall", WmApplication_Rmhd_BC_NoslipWall, WmApplication);

3. Add to CMakeLists.txt

   (src/apps/rmhd/CMakeLists.txt)

   add_sources(
       rmhd.cc
       wmapplication_rmhd.cc
       wmapplication_resistive_diffusion.cc
       wmapplication_rmhd_bc_noslipwall.cc
   )

4. Add to warpy

   (tools/warpy/apps/rmhd/rmhd_noslip_wall.py)

   from ..application import application
    
   class rmhd_noslip_wall(application):
       def __init__(self, name, on_boundaries, skin_depth_norm, nu_p_norm,gamma,
                    variable, gradient_variable, wall_velocity=[0.0,0.0,0.0], wall_magnetic_field=[0.0,0.0,0.0], components=None, gradient_variable_components=None):
           '''
           @brief
           @param name
           @param boundary_subdomains list of boundary subdomains this BC should be applied to
           @param skin_depth_norm
           @param gamma
           @param variable
           @param components
           '''
           super(rmhd_noslip_wall, self).__init__(
               name,
               req_attrs=[('Kind', 'rMHDNoSlipWall'),
                          ('OnBoundaries', on_boundaries),
                          ('skin_depth_norm', skin_depth_norm),
                          ('nu_p_norm', nu_p_norm),
                          ('gamma', gamma),
                          ("WallVelocity",wall_velocity),
                          ("WallMagneticField",wall_magnetic_field)],
               variables=[('Fluid', variable, components),
                          ('FluidGradient', gradient_variable, gradient_variable_components)])

5. Add warpy function to init

   (tools/warpy/apps/__init__.py)

   from . import rmhd

   (tools/warpy/apps/rmhd/__init__.py)

   from .rmhd_noslip_wall import rmhd_noslip_wall

6. Add warpy function to CMake

   (tools/warpy/apps/CMakeLists.txt)

   WARPY_FILES(__init__.py
       application.py
       )
   add_subdirectory(rmhd)

   (tools/warpy/apps/rmhd/CMakeLists.txt)

   WARPY_FILES(__init__.py
       rmhd.py
       rmhd_noslip_wall.py
       )

Example 3.1.3 Putting it all together

Armed with an equation set, initial conditions, and boundary conditions, we can now write a python input file. Here we setup a simple 1D problem with constant density, pressure and By and set By to 0 at the boundaries. This is effectively running the diffusion equatoin on magnetic field where the boundary conditions are pulling down the magnetic field in the middle of the domain. The input file is written:

(examples/dg/mhd/resistive_mhd_pure_diffusion.py)

# helps with making this script work with Python 2 or 3
from __future__ import division, print_function
# all WARPXM Python interface is found in the warpy package
import warpy
 
# helpers for temporal and spatial order used later
torder = 'RK4'
sorder = 'secondOrder'
# normalization params
skin_depth_norm = 1.0
nu_p_norm = 1.0# 0.1
gamma = 2.0
 
penalty_beta =  -0.5#0.0#-0.5#  -1.0#0.0#0.0#
penalty_eta = 0.0
 
imhd_numerical_flux_type = "Rusanov"
 
cfl = 1.0
 
 
# Generate a 1D mesh
# For periodic FVM simulations, the mesh must be specified with no external ghost layers (GenLayers=0)
mesh = warpy.mesh.block(Bounds=[-.5, .5],
                        NumCells=[64],
                        NodeSets=['leftWall', 'rightWall'],
                        NumLayers=1,
                        # PeriodicBoundaries=['leftWall','rightWall'],                        
                        basis_array_set=sorder) 
 
# Create a variable for the fluid
component_names = ['rho', 'px', 'py', 'pz', 'e', 'Bx', 'By', 'Bz']
fluid = warpy.variable(name='fluid',components=component_names,basis_array_set=sorder)
gradient_component_names = [prefix+"_"+suffix for prefix in component_names for suffix in ['x','y','z']]
fluid_gradient = warpy.variable(name='fluid_gradient', components=gradient_component_names,basis_array_set=sorder) 
 
# Initial conditions consist of a function and applicator
pure_diffusion_ic = warpy.functions.custom_mhd(name='pure_diffusion_ic',skin_depth_norm=skin_depth_norm,gamma=gamma)
# bandpass_ic = warpy.functions.bandpass(name='fluid_ic', width=[0.25])
# sinewave_ic = warpy.functions.fourier(name='sinewave',sine_amplitudes=1.0,sine_frequencies=2*3.1415)
# setto_ic = warpy.functions.set_to(name='setto',value=0.5)
 
# applicators:
ic_applicators = [warpy.applicator(spatial_order=sorder,
                                   fun=pure_diffusion_ic,
                                   var=fluid,components=None,
                                   spatial_scheme='Nodal')]
 
# ic_applicators = [warpy.applicator(spatial_order=sorder,fun=bandpass_ic,var=fluid,components=['By'],spatial_scheme='Nodal'),
#                   warpy.applicator(spatial_order=sorder,fun=setto_ic,var=fluid,components=['rho'],spatial_scheme='Nodal')]
 
 
# applications consist of equation sets and boundary conditions
apps = []
dirichlet_apps = []
apps.append(warpy.apps.rmhd.rmhd(name='rmhd',
                                 skin_depth_norm=skin_depth_norm,
                                 nu_p_norm=nu_p_norm,
                                 gamma=gamma,
                                 variable=fluid,
                                 gradient_variable=fluid_gradient,
                                 components=None,
                                 gradient_variable_components=None,
                                 imhd_numerical_flux_type=imhd_numerical_flux_type))
 
 
# bc
dirichlet_apps.append(warpy.apps.rmhd.rmhd_noslip_wall(name='noslip_wall_bc',
                                                       on_boundaries=['leftWall', 'rightWall'],
                                                       skin_depth_norm=skin_depth_norm,
                                                       nu_p_norm=nu_p_norm,
                                                       gamma=gamma,
                                                       wall_velocity=[0.0,0.0,0.0],
                                                       wall_magnetic_field=[0.0,0.0,0.0],
                                                       variable=fluid,
                                                       gradient_variable=fluid_gradient,
                                                       components=None,
                                                       gradient_variable_components=None))
 
# write out the fluid variable
writer = warpy.host_actions.writer(name='writer',
                                   ReadVars=[fluid,
                                             fluid_gradient] )
# Spatial Solver
spatial_solver = warpy.spatial_solvers.dg(name="dg",
                                          spatial_order=sorder,
                                          applications=apps,
                                          penalty_beta=penalty_beta,
                                          penalty_eta=penalty_eta,
                                          cfl=cfl)
 
 
 
# Variable adjusters
variable_adjusters = []
variable_adjusters.append(warpy.variable_adjusters.boundary_conditions(name='dirichlet_conditions',
                                                                       applications=dirichlet_apps,
                                                                       priority=0,
                                                                       spatial_order=sorder,
                                                                       copy_variables=[fluid],
                                                                       on_boundaries=['leftWall','rightWall']))
variable_adjusters.append(warpy.variable_adjusters.dg_gradient(name='gradient',
                                                               priority=2,
                                                               spatial_order=sorder,
                                                               penalty_beta=penalty_beta,
                                                               penalty_eta=penalty_eta,
                                                               input_variables=[fluid],
                                                               output_variables=[fluid_gradient]))
variable_adjusters.append(warpy.variable_adjusters.boundary_conditions(name='neumann',
                                                                       priority=3,
                                                                       copy_variables=[fluid_gradient],
                                                                       spatial_order=sorder,
                                                                       copy_all=True,
                                                                       on_boundaries=['leftWall','rightWall']))
 
adjust_bc_apps = []
adjust_bc_apps.append(warpy.apps.simple.bc_dirichlet(name='bc_dirichlet_grad_rhoL_grad_eL',
                                                     values=[0.0]*6,
                                                     on_boundaries=['leftWall'],
                                                     variable=fluid_gradient,
                                                     components=['rho_x','rho_y','rho_z','e_x','e_y','e_z']))
 
adjust_bc_apps.append(warpy.apps.simple.bc_dirichlet(name='bc_dirichlet_grad_rhoR_grad_eR',
                                                     values=[0.0]*6,
                                                     on_boundaries=['rightWall'],
                                                     variable=fluid_gradient,
                                                     components=['rho_x','rho_y','rho_z','e_x','e_y','e_z']))
 
variable_adjusters.append(warpy.variable_adjusters.boundary_conditions(name='bc_grad_q',
                                                                       priority=4,
                                                                       spatial_order=sorder,
                                                                       applications=adjust_bc_apps))
 
 
# Time integrator
ti = warpy.host_actions.erk(name="ti",
                            scheme=torder,
                            spatial_solvers=[spatial_solver],
                            variable_adjusters=variable_adjusters)
 
dt_controller = warpy.dt_calc.fixed_dt(init_dt=1.0e-6)
 
sim = warpy.dg_sim(name='test_rmhd_secondSorder',
# sim = warpy.dg_sim(name='test_imhd_plus_resistive_diffusion_secondSorder',                   
                   meshes=[mesh],
                   initial_conditions=ic_applicators,
                   temporal_solvers=[ti],
                   writers=[writer],
                   time=[0.,0.001],
                   dt_controller=dt_controller,
                   write_steps=100,
                   verbosity='info')
 
sim.run('mhd/resistive/jan26_2018', gen_xdmf=True, detect_nonscalar=True)

Example 3.2 Resistive MHD in 2 modules

Instead of running 1 resitive mhd application which adds the ideal and resistive diffusion terms for you, you could also run these as separate applications stacked together. You need to write a specific resistive diffusion application with only those terms:

1. Resistive Diffusion Application

   (src/apps/rmhd/wmapplication_resistive_diffusion.h)

   #ifndef WMAPPLICATION_RESISTIVE_DIFFUSION_H
   #define WMAPPLICATION_RESISTIVE_DIFFUSION_H
    
   // WarpM includes
   #include "apps/wmapplication.h"
    
   class WmApplication_Resistive_Diffusion : public WmApplication
   {
   public:
       WmApplication_Resistive_Diffusion();
    
       ~WmApplication_Resistive_Diffusion();
    
       void setup(const WxCryptSet& wxc) override;
    
       const std::vector<int>& getInputVariableIndexes(
           int flag) const override
       {
           return _inputVariables;
       }
       const std::vector<int>& getAuxiliaryVariableIndexes(
           int flag) const override
       {
           return _auxVariables;
       }
       const std::vector<int>& getOutputVariableIndexes(
           int flag) const override
       {
           return _outputVariables;
       }
       const std::vector<int>& getCrossVariableIndexes(
           int flag) const override
       {
           return _inputVariables;
       }
    
       // redefinition of fluxes from parent
       real numerical_flux(const real* q_l, const real* q_r, const real* aux_l,
           const real* aux_r, const solverVariables_t* pFV,
           real* numericalFlux) const override;
       // real internal_flux(const real * q, const real * aux, const
       // elementGeometry_t * pEG, real *internalFlux[8]) const;
       real internal_flux(const real* q, const real* aux,
           const solverVariables_t* pSV,
           std::vector<std::vector<real>>& internalFlux) const override;
    
   protected:
       real _gas_gamma;
       real _skin_depth_norm;
       real _nu_p_norm;
    
       // std::string _imhd_numerical_flux_type; // this is the type of numerical flux you want to use (rusanov, roe, hll, hlld, etc.)     
    
       std::vector<int> _inputVariables;
       std::vector<int> _auxVariables;
       std::vector<int> _outputVariables;
    
    
       real _min_density_floor;
       real _min_pressure_floor;
    
    
   private:
       WmApplication_Resistive_Diffusion& operator=(const WmApplication_Resistive_Diffusion& var);
       WmApplication_Resistive_Diffusion(const WmApplication_Resistive_Diffusion& var);
   };
    
   #endif // WMAPPLICATION_RESISTIVE_DIFFUSION

   (src/apps/rmhd/wmapplication_resistive_diffusion.cc)
   

   #include "wmapplication_resistive_diffusion.h"
   #include "rmhd.h"
    
   // STD libraries
   #include <iostream>
    
   // Wm libraries
   #include "lib/wmudgfunctions.h"
    
   WmApplication_Resistive_Diffusion::WmApplication_Resistive_Diffusion()
   {
       _allowedFlags.push_back(WMAPPLICATIONFLAG_NUMERICALFLUX);
       _allowedFlags.push_back(WMAPPLICATIONFLAG_INTERNALFLUX);
   }
    
   WmApplication_Resistive_Diffusion::~WmApplication_Resistive_Diffusion()
   {
    
   }
    
   void
   WmApplication_Resistive_Diffusion::setup(const WxCryptSet& wxc)
   {
       WmApplication::setup(wxc);
    
       WmApplication_Variable field(wxc, "Fluid", 8);
       WmApplication_Variable gradient_field(wxc, "FluidGradient", 24);
    
    
       _gas_gamma = wxc.get<real>(GAS_GAMMA);
       _skin_depth_norm = wxc.get<real>(SKIN_DEPTH_NORM);
       _nu_p_norm = wxc.get<real>(NU_P_NORM);
    
    
       _min_density_floor = 0.0;
       if (wxc.has("MinDensity"))
       {
           _min_density_floor = wxc.get<real>("MinDensity");
       }
       
       _min_pressure_floor = 0.0;
       if (wxc.has("MinPressure"))
       {
           _min_pressure_floor = wxc.get<real>("MinPressure");
       }
       
    
       // _imhd_numerical_flux_type = "Rusanov";
       // if (wxc.has("ImhdNumericalFluxType"))
       // {
       //  _imhd_numerical_flux_type = wxc.get<std::string>("ImhdNumericalFluxType");
       // }
       
    
       _inputVariables = field.getIndexes();
       _outputVariables = field.getIndexes();
       _auxVariables = gradient_field.getIndexes();
   }
    
   real
   WmApplication_Resistive_Diffusion::internal_flux(const real * q,
                              const real * aux,
                              const solverVariables_t * pSV,
                              // real *internalFlux[8]
                              std::vector< std::vector<real> >& internalFlux
       ) const
   {
    
    
       // q = [rho, px, py, pz, e, Bx, By, Bz]
       // aux = [rho_x, rho_y, rho_z,
       // px_x, px_y, px_z, py_x, py_y, py_z, pz_x, pz_y, pz_z,
       // e_x, e_y, e_z
       // Bx_x, Bx_y, Bx_z, By_x, By_y, By_z, Bz_x, Bz_y, Bz_z]
    
       // Initialize constants
       constants_rmhd_t pC;
       populate_rmhd_t(&pC, _gas_gamma, _skin_depth_norm, _nu_p_norm);
    
    
       // Diffusive portion only
    
       
       // calculate resistivity
       real eta = get_resistivity(&pC, q);
    
    
       // calculate J
       // const real * dp[3] = {aux+3, aux+6, aux+9};
       const real * dB[3] = {aux+15, aux+18, aux+21};
       real J[3];
       real abs_j_coeff = current_density(&pC,dB,J);
       // calculate internal diffusion flux
       real abs_diff_coeff = rmhd_internal_diffusion_flux(&pC,pSV,q,J,eta,internalFlux);
    
    
       
       return set_time_step(pSV, abs_j_coeff, abs_diff_coeff, eta);
       
    
   }
    
   real
   WmApplication_Resistive_Diffusion::numerical_flux(const real * q_l,
                                           const real * q_r,
                                           const real * aux_l,
                                           const real * aux_r,
                                           const solverVariables_t * pFV,
                                           real * numericalFlux) const
   {
       // q = [rho, px, py, pz, e, Bx, By, Bz]
       // aux = [rho_x, rho_y, rho_z,
       // px_x, px_y, px_z, py_x, py_y, py_z, pz_x, pz_y, pz_z,
       // e_x, e_y, e_z
       // Bx_x, Bx_y, Bx_z, By_x, By_y, By_z, Bz_x, Bz_y, Bz_z]
    
    
       // diffusive part only
       const real * dB_l[3] = {aux_l+15, aux_l+18, aux_l+21};
       const real * dB_r[3] = {aux_r+15, aux_r+18, aux_r+21};
    
       real J_l[3], J_r[3];
    
       // Initialize constants
       constants_rmhd_t pC;
       populate_rmhd_t(&pC, _gas_gamma, _skin_depth_norm, _nu_p_norm);
    
       // calculate resistivity
       real eta_l = get_resistivity(&pC, q_l);
       real eta_r = get_resistivity(&pC, q_r);
    
    
       // calculate J
       real abs_j_coeff_l = current_density(&pC,dB_l,J_l); // careful:  we might have modded density and pressure (in imhd part) but not their gradients...
       real abs_j_coeff_r = current_density(&pC,dB_r,J_r);
       
       
    
       // calculate numerical resistive diffusion flux
       real abs_diff_coeff = rmhd_numerical_diffusion_flux(&pC,pFV,q_l,q_r,J_l,J_r,eta_l,eta_r,numericalFlux);
       
    
       return set_time_step(pFV,fmax(abs_j_coeff_l,abs_j_coeff_r), abs_diff_coeff, fmax(eta_l,eta_r));
       
    
   }
    
    
   #include "lib/wxcreator.h"
    
   // Resistive_Diffusion
   REGISTER_CREATOR("resistive_diffusion", WmApplication_Resistive_Diffusion, WmApplication);

2. Registration - this is at the bottom of wmapplication_resistive_diffusion.cc (src/apps/imhd/wmapplication_resistive_diffusion.cc)

   // Resistive_Diffusion
   REGISTER_CREATOR("resistive_diffusion", WmApplication_Resistive_Diffusion, WmApplication);

3. Add to CMakeLists.txt

   (src/apps/rmhd/CMakeLists.txt)

   add_sources(
       rmhd.cc
       wmapplication_rmhd.cc
       wmapplication_resistive_diffusion.cc
       wmapplication_rmhd_bc_noslipwall.cc 
   )

4. Add to warpy

   (tools/warpy/apps/rmhd/rmhd.py)

   from ..application import application
    
   class resistive_diffusion(application):
       '''
       @brief Resistive Diffusion Terms for MHD Only
       '''
       def __init__(self, name, skin_depth_norm,nu_p_norm, gamma,
                    variable, gradient_variable, floor_density=None, floor_pressure=None, components=None, gradient_variable_components=None):
           super(resistive_diffusion, self).__init__(
               name,
               req_attrs=[('Kind', 'resistive_diffusion'),
                          ('skin_depth_norm', skin_depth_norm),
                          ('nu_p_norm', nu_p_norm),
                          ('gamma', gamma)],
               opt_attrs=[('MinDensity', floor_density),
                          ('MinPressure',floor_pressure)],            
               variables=[('Fluid', variable, components),
                          ('FluidGradient', gradient_variable, gradient_variable_components)])

5. Add warpy function to init

   (tools/warpy/apps/__init__.py)

   from . import rmhd

   (tools/warpy/apps/rmhd/__init__.py)

   from .rmhd import rmhd
   from .rmhd import resistive_diffusion

6. Add warpy function to CMake

   (tools/warpy/apps/CMakeLists.txt)

   WARPY_FILES(__init__.py
       application.py
       )
   add_subdirectory(rmhd)

   (tools/warpy/apps/rmhd/CMakeLists.txt)

   WARPY_FILES(
       __init__.py
       rmhd.py
       rmhd_noslip_wall.py
   )

Example 3.2.1 Putting it together

Now we can use the same in put file as in the previous example, except that we swap the rmhd app:

(examples/dg/mhd/resistive_mhd_pure_diffusion.py)

apps.append(warpy.apps.rmhd.rmhd(name='rmhd',
                                 skin_depth_norm=skin_depth_norm,
                                 nu_p_norm=nu_p_norm,
                                 gamma=gamma,
                                 variable=fluid,
                                 gradient_variable=fluid_gradient,
                                 components=None,
                                 gradient_variable_components=None,
                                 imhd_numerical_flux_type=imhd_numerical_flux_type))

for an imhd app followed by a resistive_diffusion app: (examples/dg/mhd/resistive_mhd_pure_diffusion.py)

apps.append(warpy.apps.imhd.imhd(name='imhd',
                                 gamma=gamma,
                                 variable=fluid,
                                 numerical_flux_type=imhd_numerical_flux_type,
                                 components=None))
 
apps.append(warpy.apps.rmhd.resistive_diffusion(name='resistive_diffusion',
                                 skin_depth_norm=skin_depth_norm,
                                 nu_p_norm=nu_p_norm,gamma=gamma,
                                 variable=fluid,
                                 gradient_variable=fluid_gradient,
                                 components=None,
                                 gradient_variable_components=None))

Example 4 - Ideal MHD with Divergence Cleaning

The ideal mhd equations ammended with divergence cleaning terms (in green), as coded in WARPXM are as follows:

Continuity

\begin{align*} \frac{\partial \rho}{\partial t} + \nabla\cdot\left(\rho \boldsymbol{u}\right) =& 0 \end{align*}

Momentum

\begin{align*} \frac{\partial\left(\rho \boldsymbol{u}\right)}{\partial t} + \nabla\cdot \left[\rho \boldsymbol{u}\boldsymbol{u} + p \mathcal{I} - \left(\boldsymbol{B}\boldsymbol{B}-\frac{1}{2}B^{2}\mathcal{I}\right) \right] =& 0 \end{align*}

Energy

\begin{align*} \frac{\partial e_{t}}{\partial t} + \nabla\cdot \left[\left(e_{t}+p + \frac{1}{2}B^{2}\right)\boldsymbol{u}-\left(\boldsymbol{u}\cdot\boldsymbol{B}\right)\boldsymbol{B} \right] =& 0 \end{align*}

Induction

\begin{align*} \frac{\partial \boldsymbol{B}}{\partial t} + \nabla\cdot\left(\boldsymbol{u}\boldsymbol{B} - \boldsymbol{B}\boldsymbol{u} {\color{green}{+\psi\mathcal{I}}} \right) =& 0 \end{align*}

\(\boldsymbol{B}\) advection

\begin{align*} \color{green} { \frac{\partial \psi}{\partial t} + \nabla\cdot\left(c_{h}^{2}\boldsymbol{B}\right) } =& \color{green} { -\frac{c_{h}^{2}}{c_{p}^{2}}\psi } \end{align*}

Altogether this is written

\begin{align*} \frac{\partial}{\partial t} \begin{pmatrix} \rho \\ \rho u \\ \rho v \\ \rho w \\ e_{t} \\ B_{x} \\ B_{y} \\ B_{z} \\ \color{green}{\psi} \end{pmatrix} + \frac{\partial}{\partial x} \begin{pmatrix} \rho u \\ \rho u u - (B_{x} B_{x} - \frac{1}{2}B^{2}) + p \\ \rho u v - B_{x} B_{y} \\ \rho u w - B_{x} B_{z} \\ \left(e_{t} + p + \frac{1}{2}B^{2}\right)u - \left(\boldsymbol{B}\cdot\boldsymbol{u}\right)B_{x} \\ 0 \color{green}{+\psi}\\ u B_{y} - B_{x} v \\ u B_{z} - B_{x} w \\ \color{green}{c_{h}^2 B_{x}} \end{pmatrix} + \frac{\partial}{\partial y} \begin{pmatrix} \rho v \\ \rho v u - B_{y} B_{x} \\ \rho v v - (B_{y} B_{y} - \frac{1}{2}B^{2}) + p \\ \rho v w - B_{y} B_{z} \\ \left(e_{t} + p + \frac{1}{2}B^{2}\right)v - \left(\boldsymbol{B}\cdot\boldsymbol{u}\right)B_{y} \\ v B_{x} - B_{y} u \\ 0 \color{green}{+\psi} \\ v B_{z} - B_{y} w \\ \color{green}{c_{h}^2 B_{y}} \end{pmatrix} + \frac{\partial}{\partial z} \begin{pmatrix} \rho w \\ \rho w u - B_{z} B_{x} \\ \rho w v - B_{z} B_{y} \\ \rho w w - (B_{z} B_{z} - \frac{1}{2}B^{2}) + p \\ \left(e_{t} + p + \frac{1}{2}B^{2}\right)w - \left(\boldsymbol{B}\cdot\boldsymbol{u}\right)B_{z} \\ w B_{x} - B_{z} u \\ w B_{y} - B_{z} v \\ 0 \color{green}{+\psi} \\ \color{green}{c_{h}^2 B_{z}} \end{pmatrix} =& \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \color{green}{-\frac{c_{h}^{2}}{c_{p}^{2}}\psi} \end{pmatrix} \end{align*}

or

\begin{align*} \begin{pmatrix} \rho \\ \rho u \\ \rho v \\ \rho w \\ e_{t} \\ B_{x} \\ B_{y} \\ B_{z} \\ \color{green}{\psi} \end{pmatrix} + \begin{pmatrix} \rho u & \rho v & \rho w \\ \rho u u - (B_{x} B_{x} - \frac{1}{2}B^{2}) + p & \rho v u - B_{y} B_{x} & \rho w u - B_{z} B_{x} \\ \rho u v - B_{x} B_{y} & \rho v v - (B_{y} B_{y} - \frac{1}{2}B^{2}) + p & \rho w v - B_{z} B_{y} \\ \rho u w - B_{x} B_{z} & \rho v w - B_{y} B_{z} & \rho w w - (B_{z} B_{z} - \frac{1}{2}B^{2}) + p \\ \left(e_{t} + p + \frac{1}{2}B^{2}\right)u - \left(\boldsymbol{B}\cdot\boldsymbol{u}\right)B_{x} & \left(e_{t} + p + \frac{1}{2}B^{2}\right)v - \left(\boldsymbol{B}\cdot\boldsymbol{u}\right)B_{y} & \left(e_{t} + p + \frac{1}{2}B^{2}\right)w - \left(\boldsymbol{B}\cdot\boldsymbol{u}\right)B_{z} \\ 0 \color{green}{+\psi} & v B_{x} - B_{y} u & w B_{x} - B_{z} u \\ u B_{y} - B_{x} v & 0 \color{green}{+\psi} & w B_{y} - B_{z} v \\ u B_{z} - B_{x} w & v B_{z} - B_{y} w & 0 \color{green}{+\psi} \\ \color{green}{c_{h}^2 B_{x}} & \color{green}{c_{h}^2 B_{y}} & \color{green}{c_{h}^2 B_{z}} \end{pmatrix} \begin{pmatrix} \frac{\partial}{\partial x} \\ \frac{\partial}{\partial y} \\ \frac{\partial}{\partial z} \end{pmatrix} =& \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \color{green}{-\frac{c_{h}^{2}}{c_{p}^{2}}\psi} \end{pmatrix} \end{align*}

NOTE: When we write the flux tensor in WARPXM we write its transpose (see code below).

1. Internal Flux, Numerical Flux, and Source

   (src/apps/imhd/wmapplication_imhd_mixed_glm.h)

   #ifndef WMAPPLICATION_IMHD_MIXED_GLM_H
   #define WMAPPLICATION_IMHD_MIXED_GLM_H
    
   // WarpM includes
   #include "apps/wmapplication.h"
    
   // IMHD includes
   #include "apps/imhd/imhd.h"
    
   class WmApplication_IMHD_MIXED_GLM : public WmApplication
   {
   public:
       WmApplication_IMHD_MIXED_GLM();
    
       ~WmApplication_IMHD_MIXED_GLM();
    
       void setup(const WxCryptSet& wxc) override;
    
       const std::vector<int>& getInputVariableIndexes(int flag) const override
       {
           return _inputVariables;
       }
       const std::vector<int>& getOutputVariableIndexes(int flag) const override
       {
           return _outputVariables;
       }
    
       // const std::vector<int>& getAuxiliaryVariableIndexes(int flag) const override
       // {
       //     return _aux_variables;
       // }
       
    
       // redefinition of fluxes from parent
       real numerical_flux(const real* q_l, const real* q_r, const real* aux_l,
           const real* aux_r, const solverVariables_t* pFV,
           real* numericalFlux) const override;
       real internal_flux(const real* q, const real* aux,
           const solverVariables_t* pSV,
           std::vector<std::vector<real>>& internalFlux) const override;
       real source(const real* q, const real* aux, const elementGeometry_t* pEG,
           real* source) const override;
       
    
   protected:
       real _gas_gamma;
       
       std::string _numerical_flux_type; // this is the type of numerical flux you want to use (rusanov, roe, hll, hlld, etc.) 
    
       std::vector<int> _inputVariables;
       std::vector<int> _outputVariables;
       // std::vector<int> _aux_variables;
    
       real _c_h, _c_p, _c_r;
    
       real _min_density_floor;
       real _min_pressure_floor;
    
    
   private:
       WmApplication_IMHD_MIXED_GLM& operator=(const WmApplication_IMHD_MIXED_GLM& var);
       WmApplication_IMHD_MIXED_GLM(const WmApplication_IMHD_MIXED_GLM& var);
   };
    
   #endif // WMAPPLICATION_IMHD_MIXED_GLM_H

   (src/apps/imhd/wmapplication_imhd_mixed_glm.cc)

   #include "wmapplication_imhd_mixed_glm.h"
    
   #include "imhd_glm.h"
    
   // STD libraries
   #include <iostream>
   #include <cmath>
    
   // Wm libraries
   #include "lib/wmudgfunctions.h"
    
   WmApplication_IMHD_MIXED_GLM::WmApplication_IMHD_MIXED_GLM()
   {
       _allowedFlags.push_back(WMAPPLICATIONFLAG_NUMERICALFLUX);
       _allowedFlags.push_back(WMAPPLICATIONFLAG_INTERNALFLUX);
       _allowedFlags.push_back(WMAPPLICATIONFLAG_SOURCE);  
   }
    
   WmApplication_IMHD_MIXED_GLM::~WmApplication_IMHD_MIXED_GLM()
   {
    
   }
    
   void
   WmApplication_IMHD_MIXED_GLM::setup(const WxCryptSet& wxc)
   {
       WmApplication::setup(wxc);
    
       WmApplication_Variable field(wxc, "Fluid", 9);
    
    
       _gas_gamma = wxc.get<real>(GAS_GAMMA);
    
    
       _numerical_flux_type = "Rusanov";
       if (wxc.has("NumericalFluxType"))
       {
           _numerical_flux_type = wxc.get<std::string>("NumericalFluxType");
       }
       
       _min_density_floor = 0.0;
       if (wxc.has("MinDensity"))
       {
           _min_density_floor = wxc.get<real>("MinDensity");
       }
       
       _min_pressure_floor = 0.0;
       if (wxc.has("MinPressure"))
       {
           _min_pressure_floor = wxc.get<real>("MinPressure");
       }   
       
       _c_h = 1.0;
       if (wxc.has("C_H"))
       {
           _c_h = wxc.get<real>("C_H");
       }   
       _c_r = 0.18;
       _c_p = std::sqrt(_c_r*_c_h);
       
       _inputVariables = field.getIndexes();
       _outputVariables = field.getIndexes();
    
   }
    
   real
   WmApplication_IMHD_MIXED_GLM::numerical_flux(const real * q_l,
                                           const real * q_r,
                                           const real * aux_l,
                                           const real * aux_r,
                                           const solverVariables_t * pFV,
                                           real * numericalFlux) const
   {
       // Define the constants
       constants_imhd_glm_t pC;
       populate_constants_imhd_glm_t(&pC, _gas_gamma, _c_h, _c_p, _c_r);
       real initial_flux[8];
    
    
       real dens_q_l[9], dens_q_r[9];
       real positivized_q_l[9], positivized_q_r[9];    
       ensure_floor_density(q_l,9,_min_density_floor,dens_q_l);
       ensure_floor_density(q_r,9,_min_density_floor,dens_q_r);    
       ensure_floor_pressure(dens_q_l,9,_gas_gamma,_min_pressure_floor,positivized_q_l);
       ensure_floor_pressure(dens_q_r,9,_gas_gamma,_min_pressure_floor,positivized_q_r);   
       
    
    
       // lets assign Bxm as the value for the numerical flux to use
       // First need to rotate
       real r_q_l[9];              // rotated q_l
       real r_q_r[9];              // rotated q_r
       /* initialize scalars rho and e. No need for vectors p and B since they get rotated in next step */
       r_q_l[0] = positivized_q_l[0];          /* initialize rho on left */
       r_q_l[4] = positivized_q_l[4];          /* initialize e on left */
       r_q_l[8] = positivized_q_l[8];          // initialize psi on left
       r_q_r[0] = positivized_q_r[0];          /* ditto rho left */
       r_q_r[4] = positivized_q_r[4];          /* ditto e right */
       r_q_r[8] = positivized_q_r[8];          // initialize psi on right
       /* local face frame rotation */
       rotate_vector(q_l+1, pFV, r_q_l+1);
       rotate_vector(q_l+5, pFV, r_q_l+5);
       rotate_vector(q_r+1, pFV, r_q_r+1);
       rotate_vector(q_r+5, pFV, r_q_r+5);
       real Bxm = calc_Bxm(r_q_l,r_q_r,_c_h);
       r_q_l[5] = r_q_r[5] = Bxm;  // apply new Bx here
       real mod_q_l[9];
       real mod_q_r[9];
       mod_q_l[0] = r_q_l[0];
       mod_q_l[4] = r_q_l[4];
       mod_q_l[8] = r_q_l[8];      // not really going to be used
       mod_q_r[0] = r_q_r[0];
       mod_q_r[4] = r_q_r[4];
       mod_q_r[8] = r_q_r[8];      // not really going to be used
       antirotate_vector(r_q_l+1,pFV,mod_q_l+1);
       antirotate_vector(r_q_l+5,pFV,mod_q_l+5);   
       antirotate_vector(r_q_r+1,pFV,mod_q_r+1);
       antirotate_vector(r_q_r+5,pFV,mod_q_r+5);   
       
       real c = imhd_numerical_flux(&(pC.imhd),mod_q_l, mod_q_r, aux_l, aux_r, pFV, _numerical_flux_type, initial_flux);   
    
    
    
       real dedner_flux[9];
       imhd_glm_numerical_flux(&pC,positivized_q_l, positivized_q_r, aux_l, aux_r, pFV, dedner_flux); // should be based on unmodified q
    
       
       // add up initial imhd flux and Dedner correction flux in global frame
       for (int i =0;i<8;i++)
       {
           numericalFlux[i] = initial_flux[i] + dedner_flux[i];
       }
       numericalFlux[8] = dedner_flux[8];  // incorporate the last component as this is not in initial_flux
    
    
    
       return pFV->cfl_max * pFV->dxn / std::max(fabs(c),fabs(_c_h));  
   }
    
   real
   WmApplication_IMHD_MIXED_GLM::internal_flux(const real * q,
                                          const real * aux,
                                          // const elementGeometry_t * pEG,
                                          const solverVariables_t * pSV,
                                          // real *internalFlux[8]
                                          std::vector< std::vector<real> >& internalFlux
       ) const
   {
    
       // q = [rho, px, py, pz, e, Bx, By, Bz, psi]
       
       constants_imhd_glm_t pC;
       populate_constants_imhd_glm_t(&pC, _gas_gamma, _c_h, _c_p, _c_r);
    
       // ensure positivity of density and pressure
       real dens_q[9]; 
       real positivized_q[9];
       ensure_floor_density(q,9,_min_density_floor,dens_q);
       ensure_floor_pressure(dens_q,9,_gas_gamma,_min_pressure_floor,positivized_q);
    
       
       /* Ideal portion */
       // output vector allocation is [num_dims] x [_num_output_Q_max]
       std::vector<real> N_output_dimension_vector(9, 0.);
       std::vector< std::vector<real> > internal_ideal_flux(pSV->num_dims, N_output_dimension_vector);
       // apply flux (only first 8 components will be filled)
       imhd_internal_flux(&(pC.imhd), pSV, positivized_q, internal_ideal_flux);
    
       
    
    
       /* GLM portion */
       // output vector allocation is [num_dims] x [_num_output_Q_max]
       N_output_dimension_vector.assign(9, 0.);
       std::vector< std::vector<real> > internal_glm_flux(pSV->num_dims, N_output_dimension_vector);
       // apply glm flux
       imhd_glm_internal_flux(&pC, pSV, positivized_q, internal_glm_flux); // the postiived part doesnt affect the glm flux so doesnt matter if use q or positivized_q
    
    
       /* Add the 2 fluxes together */
       for(int i=0;i<pSV->num_dims;i++)
       {
           for(int j=0;j<9;j++)
           {
               internalFlux[i][j] = internal_ideal_flux[i][j] + internal_glm_flux[i][j];
           }
    
       }
       
       
       
       // Speeds should be calculated in numerical flux
       return INFINITY;
    
   }
    
   real
   WmApplication_IMHD_MIXED_GLM::source(
       const real* q, const real* aux, const elementGeometry_t* pEG, real* source) const
   {
       // q = [rho, px, py, pz, e, Bx, By, Bz, psi]
    
       real psi = q[8];
    
       source[0] = 0.0;
    
       source[1] = 0.0;
       source[2] = 0.0;
       source[3] = 0.0;
    
       source[4] = 0.0;
    
       source[5] = 0.0;
       source[6] = 0.0;
       source[7] = 0.0;
    
       source[8] = -_c_h*_c_h/(_c_p*_c_p)*psi;
    
    
       // Use the check frequency application
       return INFINITY;
   }
    
    
   #include "lib/wxcreator.h"
    
   // IMHD flux
   REGISTER_CREATOR("iMHD_MIXED_GLM", WmApplication_IMHD_MIXED_GLM, WmApplication);

   where we are using helper functions for internal fluxes for both the ideal and divergence cleaning corrections. For the ideal mhd internal flux we reuse the the function shown above for ideal mhd in (src/apps/imhd/imhd.cc). For the divergence cleaning we write another function specifically for those terms

   (src/apps/imhd/imhd_glm.cc)
   

   void imhd_glm_internal_flux(
       const constants_imhd_glm_t * pC,
       const solverVariables_t * pSV,
       const real* q,
       std::vector< std::vector<real> >& internalFlux
       )
   {
       // add div (psi * I) term and flux for psi equation which is c_h^2 * div (B)
       real psi = q[8];
       real bx = q[5];
       real by = q[6];
       real bz = q[7];
       real ch_factor = pC->c_h*pC->c_h;
    
       /* clear the tensor */
       for (int i=0;i<pSV->num_dims;i++)
       {
           for (int j=0;j<9;j++)
           {
               internalFlux[i][j] = 0.;
               
           }
           
       }
    
       // add components
       internalFlux[0][5] += psi;  // induction term
       internalFlux[0][8] = ch_factor * bx; // psi term
       if (pSV->num_dims > 1)
       {
           internalFlux[1][6] += psi; // induction term
           internalFlux[1][8] = ch_factor * by; // psi term        
       }
       if (pSV->num_dims > 2)
       {
           internalFlux[2][7] += psi;
           internalFlux[2][8] = ch_factor * bz; // psi term                
       }
    
   }

   The numerical flux for ideal mhd again done with a reuse from (src/apps/imhd/imhd.cc) where the specific flux functions are also given in (src/apps/imhd/imhd.cc), which can be looked up. The divergence cleaning addition is added in a file

     (src/apps/imhd/imhd_glm.cc)
   

   void imhd_glm_numerical_flux(
       const constants_imhd_glm_t * pC,
       const real * q_l,
       const real * q_r,
       const real * aux_l,
       const real * aux_r,
       const solverVariables_t * pFV,
       real * numericalFlux)
   {
       
       // add correction of Dedner:
       // First need to rotate
       real r_q_l[9];              // rotated q_l
       real r_q_r[9];              // rotated q_r
       real r_dedner_flux[9];              // rotated flux
       // real dedner_flux[9];             // antirotated flux
    
       /* initialize scalars rho and e. No need for vectors p and B since they get rotated in next step */
       r_q_l[0] = q_l[0];          /* initialize rho on left */
       r_q_l[4] = q_l[4];          /* initialize e on left */
       r_q_l[8] = q_l[8];          // initialize psi on left
       r_q_r[0] = q_r[0];          /* ditto rho left */
       r_q_r[4] = q_r[4];          /* ditto e right */
       r_q_r[8] = q_r[8];          // initialize psi on right
    
       /* local face frame rotation */
       rotate_vector(q_l+1, pFV, r_q_l+1);
       rotate_vector(q_l+5, pFV, r_q_l+5);
       rotate_vector(q_r+1, pFV, r_q_r+1);
       rotate_vector(q_r+5, pFV, r_q_r+5);
    
       // calculate intermediate states as per Dedner
       // real Bxm = r_q_l[5] + 0.5*(r_q_r[5]-r_q_l[5]) - 0.5/pC->c_h * (r_q_r[8]-r_q_l[8]);
       real Bxm = calc_Bxm(r_q_l,r_q_r,pC->c_h);
       // real Bxm = 0.5*(r_q_l[5]+r_q_r[5]);
       real psi_m = r_q_l[8] + 0.5*(r_q_r[8] - r_q_l[8]) - 0.5*pC->c_h*(r_q_r[5]-r_q_l[5]);
       
       // calculate in-frame flux terms, first initialize to 0
       for (int i=0;i<9;i++)
       {
           r_dedner_flux[i] = 0.0;
           numericalFlux[i] = 0.0; // this is the dedner flux after rotation
       }
       // apply Dedner numerical fluxes on local Bx and psi
       r_dedner_flux[5] = psi_m;
       r_dedner_flux[8] = pC->c_h*pC->c_h*Bxm;
    
       /* rotate back momentum and mag field */
       antirotate_vector(r_dedner_flux+1, pFV, numericalFlux+1);
       antirotate_vector(r_dedner_flux+5, pFV, numericalFlux+5);
    
       
    
       
   }

2. Registration - this is at the bottom of wmapplication_imhd_mixed_glm.cc

   (src/apps/imhd/wmapplication_imhd_mixed_glm.cc)

   #include "lib/wxcreator.h"
    
   // IMHD flux
   REGISTER_CREATOR("iMHD_MIXED_GLM", WmApplication_IMHD_MIXED_GLM, WmApplication);

3. Add to CMakeLists.txt

   (src/apps/imhd/CMakeLists.txt)

   add_sources(
       imhd.cc
       imhd_glm.cc
       wmapplication_imhd_mixed_glm.cc
   )

4. Add to warpy

   (tools/warpy/apps/imhd/imhd.py)

   class imhd_mixed_glm(application):
       '''
       @brief Ideal MHD with Mixed GLM formulation as per Dedner
       '''
       def __init__(self, name, gamma,
                    variable, c_h=1.0,numerical_flux_type='Rusanov',components=None, floor_density=None, floor_pressure=None):
           super(imhd_mixed_glm, self).__init__(
               name,
               req_attrs=[('Kind', 'iMHD_MIXED_GLM'),
                          ('gamma', gamma),
                          ('C_H', c_h),                       
                          ('NumericalFluxType',numerical_flux_type)],
               opt_attrs=[('MinDensity', floor_density),
                          ('MinPressure',floor_pressure)],            
               variables=[('Fluid', variable, components)])        

5. Add warpy function to init

   (tools/warpy/apps/__init__.py)

   from . import imhd

   (tools/warpy/apps/imhd/__init__.py)

   from .imhd import imhd_mixed_glm

6. Add warpy function to CMake

   (tools/warpy/apps/CMakeLists.txt)

   WARPY_FILES(__init__.py
       application.py
       )
   add_subdirectory(imhd)

   (tools/warpy/apps/imhd/CMakeLists.txt)

   WARPY_FILES(
       __init__.py
       imhd.py
   )