From b31e133fdd0abe79b184bad6ab005b0d5a6cbb91 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Tue, 4 Nov 2025 17:30:33 -0800 Subject: [PATCH] freeze_profiles: Add option to fix the Z-average profiles Removes the DC, i.e. Z-averaged, component of the time-derivatives. This effectively adds artificial, time-varying, sources and sinks to keep the profiles fixed. Not recommended for nonlinear simulations, but useful for linear analysis. --- include/evolve_density.hxx | 2 ++ include/evolve_energy.hxx | 2 ++ include/evolve_momentum.hxx | 2 ++ include/evolve_pressure.hxx | 2 ++ include/vorticity.hxx | 2 ++ .../implement_new_cooling_curves.ipynb | 13 ++----------- src/evolve_density.cxx | 6 ++++++ src/evolve_energy.cxx | 6 ++++++ src/evolve_momentum.cxx | 6 ++++++ src/evolve_pressure.cxx | 6 ++++++ src/vorticity.cxx | 6 ++++++ 11 files changed, 42 insertions(+), 11 deletions(-) diff --git a/include/evolve_density.hxx b/include/evolve_density.hxx index c2a7c15cc..a5d6f6fb4 100644 --- a/include/evolve_density.hxx +++ b/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) diff --git a/include/evolve_energy.hxx b/include/evolve_energy.hxx index 9eb0c903e..53e30e404 100644 --- a/include/evolve_energy.hxx +++ b/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 diff --git a/include/evolve_momentum.hxx b/include/evolve_momentum.hxx index c13f91571..9ff4181f4 100644 --- a/include/evolve_momentum.hxx +++ b/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 diff --git a/include/evolve_pressure.hxx b/include/evolve_pressure.hxx index 0c8f0869d..699927fbf 100644 --- a/include/evolve_pressure.hxx +++ b/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] diff --git a/include/vorticity.hxx b/include/vorticity.hxx index acf2c2815..6495de9bd 100644 --- a/include/vorticity.hxx +++ b/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 diff --git a/scripts/impurity_curves/implement_new_cooling_curves.ipynb b/scripts/impurity_curves/implement_new_cooling_curves.ipynb index 52e1bbafd..fac068317 100644 --- a/scripts/impurity_curves/implement_new_cooling_curves.ipynb +++ b/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", diff --git a/src/evolve_density.cxx b/src/evolve_density.cxx index 672d1ea6b..2fa5dada6 100644 --- a/src/evolve_density.cxx +++ b/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(false); + evolve_log = options["evolve_log"] .doc("Evolve the logarithm of density?") .withDefault(false); @@ -311,6 +313,10 @@ void EvolveDensity::finally(const Options& state) { flow_ylow += get(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) { diff --git a/src/evolve_energy.cxx b/src/evolve_energy.cxx index c0363b4a5..f2d580358 100644 --- a/src/evolve_energy.cxx +++ b/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(false); + diagnose = options["diagnose"] .doc("Save additional output diagnostics") .withDefault(false); @@ -465,6 +467,10 @@ void EvolveEnergy::finally(const Options& state) { ddt(E) *= get(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; } diff --git a/src/evolve_momentum.cxx b/src/evolve_momentum.cxx index 24df824dc..9e054821c 100644 --- a/src/evolve_momentum.cxx +++ b/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(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) { diff --git a/src/evolve_pressure.cxx b/src/evolve_pressure.cxx index 428679330..a261da8e7 100644 --- a/src/evolve_pressure.cxx +++ b/src/evolve_pressure.cxx @@ -88,6 +88,8 @@ EvolvePressure::EvolvePressure(std::string name, Options& alloptions, Solver* so .doc("4th-order dissipation of temperature") .withDefault(-1.0); + freeze_profiles = options["freeze_profiles"].doc("Subtract Z average from time derivatives?").withDefault(false); + diagnose = options["diagnose"] .doc("Save additional output diagnostics") .withDefault(false); @@ -558,6 +560,10 @@ void EvolvePressure::finally(const Options& state) { ddt(P) *= get(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; } diff --git a/src/vorticity.cxx b/src/vorticity.cxx index bc2c42d3a..b791ca77b 100644 --- a/src/vorticity.cxx +++ b/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(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) {