Merge branch 'next' into next-update

Author: bshanahan

diffusion2d.cxx

@@ -23,12 +23,6 @@ Diffusion2D::Diffusion2D(Solver *solver, Mesh*, Options &options) : NeutralModel
   inv = 0;
 }
 
-Diffusion2D::~Diffusion2D() {
-  if(inv) {
-    delete inv;
-  }
-}
-
 void Diffusion2D::update(const Field3D &Ne, const Field3D &Te, const Field3D &UNUSED(Ti), const Field3D &UNUSED(Vi)) {
 
   mesh->communicate(Nn, Pn);

diffusion2d.hxx

@@ -8,7 +8,7 @@
 class Diffusion2D : public NeutralModel {
 public:
   Diffusion2D(Solver *solver, Mesh *mesh, Options &options);
-  ~Diffusion2D();
+  ~Diffusion2D() {}
 
   void update(const Field3D &Ne, const Field3D &Te, const Field3D &Ti, const Field3D &Vi);
 
@@ -20,7 +20,7 @@ private:
 
   BoutReal Lmax; // Maximum mean free path [m]
 
-  Laplacian *inv; // Laplacian inversion used for preconditioning
+  std::unique_ptr<Laplacian> inv; // Laplacian inversion used for preconditioning
 };
 
 #endif // __NEUTRAL_DIFFUSION2D_H__

hermes-2.cxx

@@ -359,6 +359,7 @@ int Hermes::init(bool restarting) {
 
   OPTION(optsc, ne_hyper_z, -1.0);
   OPTION(optsc, pe_hyper_z, -1.0);
+  OPTION(optsc, pi_hyper_z, -1.0);
 
   OPTION(optsc, low_n_diffuse, false);
   OPTION(optsc, low_n_diffuse_perp, false);
@@ -701,7 +702,6 @@ int Hermes::init(bool restarting) {
           .doc("Include a fixed fraction carbon impurity. < 0 means none.")
           .withDefault(-1.);
   if (carbon_fraction > 0.0) {
-    SAVE_REPEAT(Rzrad);
     SAVE_ONCE(carbon_fraction);
     carbon_rad = new HutchinsonCarbonRadiation();
   }
@@ -1870,6 +1870,7 @@ int Hermes::rhs(BoutReal t) {
           Vi(r.ind, jy, jz) = -Vi(r.ind, mesh->ystart, jz);
           NVi(r.ind, jy, jz) = -NVi(r.ind, mesh->ystart, jz);
           Ve(r.ind, jy, jz) = -Ve(r.ind, mesh->ystart, jz);
+          VePsi(r.ind, jy, jz) = -VePsi(r.ind, mesh->ystart, jz);
           Jpar(r.ind, jy, jz) = -Jpar(r.ind, mesh->ystart, jz);
         }
       }
@@ -3428,6 +3429,10 @@ int Hermes::rhs(BoutReal t) {
         (2. / 3) * FV::Div_a_Laplace_perp(anomalous_chi * DC(Ne), DC(Ti));
   }
 
+  if (pi_hyper_z > 0.0) {
+    ddt(Pi) -= pi_hyper_z * SQ(SQ(coord->dz)) * D4DZ4(Pi);
+  }
+  
   ///////////////////////////////////
   // Heat transmission through sheath
 
@@ -4007,7 +4012,7 @@ int Hermes::rhs(BoutReal t) {
  * @param[in] delta   Not used here
  */
 int Hermes::precon(BoutReal t, BoutReal gamma, BoutReal delta) {
-  static InvertPar *inv = NULL;
+  static std::unique_ptr<InvertPar> inv;
   if (!inv) {
     // Initialise parallel inversion class
     auto inv = InvertPar::create();

hermes-2.hxx

@@ -183,7 +183,7 @@ private:
   BoutReal z_hyper_viscos, x_hyper_viscos, y_hyper_viscos; // 4th-order derivatives
   bool low_n_diffuse; // Diffusion in parallel direction at low density
   bool low_n_diffuse_perp; // Diffusion in perpendicular direction at low density
-  BoutReal ne_hyper_z, pe_hyper_z; // Hyper-diffusion
+  BoutReal ne_hyper_z, pe_hyper_z, pi_hyper_z; // Hyper-diffusion
   BoutReal scale_num_cs; // Scale numerical sound speed
   BoutReal floor_num_cs; // Apply a floor to the numerical sound speed
   bool vepsi_dissipation; // Dissipation term in VePsi equation
@@ -230,8 +230,8 @@ private:
   // Electromagnetic solver for finite electron mass case
   bool split_n0_psi;   // Split the n=0 component of Apar (psi)?
   //Laplacian *aparSolver;
-  std::unique_ptr<LaplaceXZ> aparSolver{nullptr};
 
+  std::unique_ptr<LaplaceXZ> aparSolver{nullptr};
   LaplaceXY *aparXY;    // Solves n=0 component
   Field2D psi2D;        // Axisymmetric Psi
 
@@ -246,6 +246,7 @@ private:
   BoutReal phi_boundary_last_update; ///< The last time the boundary was updated
 
   bool newXZsolver; 
+
   std::unique_ptr<Laplacian> phiSolver{nullptr}; // Old Laplacian in X-Z
   std::unique_ptr<LaplaceXZ> newSolver{nullptr}; // New Laplacian in X-Z