[WIP] Add freeze_profiles option to fix the Z-average profiles

include/evolve_density.hxx

@@ -76,6 +76,8 @@ private:
   bool low_p_diffuse_perp; ///< Add artificial cross-field diffusion at low pressure?
   BoutReal hyper_z;    ///< Hyper-diffusion in Z
 
+  bool freeze_profiles; ///< Subtract Z average from time derivatives?
+
   bool evolve_log; ///< Evolve logarithm of density?
   Field3D logN;    ///< Logarithm of density (if evolving)
 

include/evolve_energy.hxx

@@ -104,6 +104,8 @@ private:
 
   BoutReal hyper_z; ///< Hyper-diffusion
 
+  bool freeze_profiles; ///< Subtract Z average from time derivatives?
+
   bool diagnose;      ///< Output additional diagnostics?
   bool enable_precon; ///< Enable preconditioner?
   Field3D flow_xlow, flow_ylow; ///< Energy flow diagnostics

include/evolve_momentum.hxx

@@ -51,6 +51,8 @@ private:
 
   BoutReal hyper_z;  ///< Hyper-diffusion
 
+  bool freeze_profiles; ///< Subtract Z average from time derivatives?
+
   bool diagnose; ///< Output additional diagnostics?
   bool fix_momentum_boundary_flux; ///< Fix momentum flux to boundary condition?
   Field3D flow_xlow, flow_ylow; ///< Momentum flow diagnostics

include/evolve_pressure.hxx

@@ -109,6 +109,8 @@ private:
   BoutReal hyper_z; ///< Hyper-diffusion
   BoutReal hyper_z_T; ///< 4th-order dissipation in T
 
+  bool freeze_profiles; ///< Subtract Z average from time derivatives?
+
   bool diagnose; ///< Output additional diagnostics?
   bool enable_precon; ///< Enable preconditioner?
   BoutReal source_normalisation; ///< Normalisation factor [Pa/s]

include/vorticity.hxx

@@ -140,6 +140,8 @@ private:
   BoutReal hyper_z; ///< Hyper-viscosity in Z
   Field2D viscosity; ///< Kinematic viscosity
 
+  bool freeze_profiles; ///< Subtract Z average from time derivatives?
+
   // Diagnostic outputs
   Field3D DivJdia, DivJcol; // Divergence of diamagnetic and collisional current
 

scripts/impurity_curves/implement_new_cooling_curves.ipynb

@@ -2,18 +2,9 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 50,
+   "execution_count": 1,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "The autoreload extension is already loaded. To reload it, use:\n",
-      "  %reload_ext autoreload\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "%matplotlib inline\n",
     "import pandas as pd\n",

src/evolve_density.cxx

@@ -48,6 +48,8 @@ EvolveDensity::EvolveDensity(std::string name, Options& alloptions, Solver* solv
 
   hyper_z = options["hyper_z"].doc("Hyper-diffusion in Z").withDefault(-1.0);
 
+  freeze_profiles = options["freeze_profiles"].doc("Subtract Z average from time derivatives?").withDefault<bool>(false);
+
   evolve_log = options["evolve_log"]
                    .doc("Evolve the logarithm of density?")
                    .withDefault<bool>(false);
@@ -311,6 +313,10 @@ void EvolveDensity::finally(const Options& state) {
       flow_ylow += get<Field3D>(species["particle_flow_ylow"]);
     }
   }
+
+  if (freeze_profiles) {
+    ddt(N) -= DC(ddt(N)); // Remove the DC, i.e. Z-averaged, component.
+  }
 }
 
 void EvolveDensity::outputVars(Options& state) {

src/evolve_energy.cxx

@@ -81,6 +81,8 @@ EvolveEnergy::EvolveEnergy(std::string name, Options& alloptions, Solver* solver
 
   hyper_z = options["hyper_z"].doc("Hyper-diffusion in Z").withDefault(-1.0);
 
+  freeze_profiles = options["freeze_profiles"].doc("Subtract Z average from time derivatives?").withDefault<bool>(false);
+
   diagnose = options["diagnose"]
                  .doc("Save additional output diagnostics")
                  .withDefault<bool>(false);
@@ -465,6 +467,10 @@ void EvolveEnergy::finally(const Options& state) {
     ddt(E) *= get<Field3D>(state["scale_timederivs"]);
   }
 
+  if (freeze_profiles) {
+    ddt(E) -= DC(ddt(E)); // Remove the DC, i.e. Z-averaged, component.
+  }
+
   if (evolve_log) {
     ddt(logE) = ddt(E) / E;
   }

src/evolve_momentum.cxx

@@ -45,6 +45,8 @@ EvolveMomentum::EvolveMomentum(std::string name, Options &alloptions, Solver *so
 
   hyper_z = options["hyper_z"].doc("Hyper-diffusion in Z").withDefault(-1.0);
 
+  freeze_profiles = options["freeze_profiles"].doc("Subtract Z average from time derivatives?").withDefault<bool>(false);
+
   V.setBoundary(std::string("V") + name);
 
   diagnose = options["diagnose"]
@@ -245,6 +247,10 @@ void EvolveMomentum::finally(const Options &state) {
   // Note: Copy boundary condition so dump file has correct boundary.
   NV_solver.setBoundaryTo(NV);
   NV = NV_solver;
+
+  if (freeze_profiles) {
+    ddt(NV) -= DC(ddt(NV)); // Remove the DC, i.e. Z-averaged, component.
+  }
 }
 
 void EvolveMomentum::outputVars(Options &state) {

src/evolve_pressure.cxx

@@ -88,6 +88,8 @@ EvolvePressure::EvolvePressure(std::string name, Options& alloptions, Solver* so
     .doc("4th-order dissipation of temperature")
     .withDefault<BoutReal>(-1.0);
 
+  freeze_profiles = options["freeze_profiles"].doc("Subtract Z average from time derivatives?").withDefault<bool>(false);
+
   diagnose = options["diagnose"]
     .doc("Save additional output diagnostics")
     .withDefault<bool>(false);
@@ -558,6 +560,10 @@ void EvolvePressure::finally(const Options& state) {
     ddt(P) *= get<Field3D>(state["scale_timederivs"]);
   }
 
+  if (freeze_profiles) {
+    ddt(P) -= DC(ddt(P)); // Remove the DC, i.e. Z-averaged, component.
+  }
+
   if (evolve_log) {
     ddt(logP) = ddt(P) / P;
   }

src/vorticity.cxx

@@ -101,6 +101,8 @@ Vorticity::Vorticity(std::string name, Options& alloptions, Solver* solver) {
 
   hyper_z = options["hyper_z"].doc("Hyper-viscosity in Z. < 0 -> off").withDefault(-1.0);
 
+  freeze_profiles = options["freeze_profiles"].doc("Subtract Z average from time derivatives?").withDefault<bool>(false);
+
   // Numerical dissipation terms
   // These are required to suppress parallel zig-zags in
   // cell centred formulations. Essentially adds (hopefully small)
@@ -759,6 +761,10 @@ void Vorticity::finally(const Options& state) {
       }
     }
   }
+
+  if (freeze_profiles) {
+    ddt(Vort) -= DC(ddt(Vort)); // Remove the DC, i.e. Z-averaged, component.
+  }
 }
 
 void Vorticity::outputVars(Options& state) {