From d4d582a86d165a0d616d5483b94e35a123244e62 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Sat, 8 Oct 2022 12:34:43 +0100 Subject: [PATCH] v2.1.0 Use BoutEquation to allow individual terms in evolution equations to be saved. Better default values for input options. --- .gitignore | 5 +- .gitmodules | 3 + build_tools/save_git_version.py | 42 +++ shared/BoutEquation | 1 + storm3d/D-vpar.cxx | 62 ++-- storm3d/D-vpar.hxx | 9 +- storm3d/boundaries.cxx | 266 +++++++------- storm3d/data/BOUT.inp | 17 +- storm3d/initialise.cxx | 15 +- storm3d/makefile | 16 +- storm3d/neutral-model.cxx | 4 +- storm3d/neutral-model.hxx | 6 +- storm3d/storm.cxx | 339 +++++++++++------- storm3d/storm.hxx | 21 +- storm3d/utilities.cxx | 24 +- tests3d/run_3d_filament_test.py | 2 +- tests3d/test-3d/data/BOUT.inp | 54 +-- tests3d/test-equilib/data/BOUT.inp | 54 +-- .../data/expected_equilib/T_eq.dat | Bin 160 -> 160 bytes .../data/expected_equilib/U_eq.dat | Bin 160 -> 160 bytes .../data/expected_equilib/V_eq.dat | Bin 160 -> 160 bytes .../data/expected_equilib/equilibrium.nc | Bin 44432 -> 44432 bytes .../data/expected_equilib/n_eq.dat | Bin 160 -> 160 bytes .../data/expected_equilib/phi_eq.dat | 3 +- 24 files changed, 496 insertions(+), 447 deletions(-) create mode 100644 .gitmodules create mode 100755 build_tools/save_git_version.py create mode 160000 shared/BoutEquation diff --git a/.gitignore b/.gitignore index aace01c..58db75c 100644 --- a/.gitignore +++ b/.gitignore @@ -59,8 +59,9 @@ BOUT.settings *.bmp *.pdf -storm3d/storm -*/make.config +make.config +*storm data_* *.txt storm2d/storm2d +storm3d/storm_version.hxx diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..845e9e3 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "shared/BoutEquation"] + path = shared/BoutEquation + url = git@github.com:johnomotani/BoutEquation diff --git a/build_tools/save_git_version.py b/build_tools/save_git_version.py new file mode 100755 index 0000000..914e9b5 --- /dev/null +++ b/build_tools/save_git_version.py @@ -0,0 +1,42 @@ +#!/usr/bin/env python3 + +from pathlib import Path +from string import ascii_uppercase + +try: + from git import Repo +except ImportError: + # print a return code to stdout so we can capture and test it in the makefile + print(11) + raise + +scriptdir = Path(__file__).parent + +repo = Repo(scriptdir.parent) + +diff = repo.git.diff() + +# Find a string to use for the delimiter +delimiter = "_STORM_DIFF_" +for letter in ascii_uppercase: + if ")" + delimiter not in diff: + delimiter_success = True + break + delimiter = "_STORM_DIFF_" + letter + "_" +if ")" + delimiter in diff: + # print a return code to stdout so we can capture and test it in the makefile + print(12) + raise ValueError( + "save_git_version.py failed to find a delimiter that is not in the git diff" + ) + +with open("storm_version.hxx", "w") as f: + f.write("constexpr auto storm_git_hash = \"") + f.write(repo.git.describe(abbrev=40, dirty=True, always=True, tags=True)) + f.write("\";\n") + f.write("constexpr auto storm_git_diff = R\"" + delimiter + "(") + f.write(diff) + f.write(")" + delimiter + "\";\n") + +# print a return code to stdout so we can capture and test it in the makefile +print(0) diff --git a/shared/BoutEquation b/shared/BoutEquation new file mode 160000 index 0000000..e29baa1 --- /dev/null +++ b/shared/BoutEquation @@ -0,0 +1 @@ +Subproject commit e29baa101661921472dc343579768fb60ae02f64 diff --git a/storm3d/D-vpar.cxx b/storm3d/D-vpar.cxx index 68c287b..2eb8163 100644 --- a/storm3d/D-vpar.cxx +++ b/storm3d/D-vpar.cxx @@ -5,8 +5,8 @@ #include #include -NeutralDVpar::NeutralDVpar(Solver *solver, Options &options) - : NeutralModel(options), my_output_monitor(this) { +NeutralDVpar::NeutralDVpar(Solver *solver, Options &options, Datafile &dump) + : NeutralModel(options, dump), my_output_monitor(this) { mesh = bout::globals::mesh; coordinates_centre = mesh->getCoordinates(CELL_CENTRE); @@ -165,8 +165,7 @@ void NeutralDVpar::update(const Field3D& n, const Field3D& n_stag, const Field3D const Field3D& V, const Field3D& T, const Field3D& T_stag, BoutReal time) { TRACE("NeutralDVpar::update"); - ddt(lognn) = 0.; - ddt(nvn) = 0.; + rhs_counter++; Tn = Tn0/T0; Tn_stag = Tn0/T0; @@ -278,24 +277,31 @@ void NeutralDVpar::update(const Field3D& n, const Field3D& n_stag, const Field3D } Field3D nn_vperp_GradPerpvn = nn_stag*vnperp_GradPerpvn; - ddt(lognn) += - vn_centre*myGrad_par(lognn_aligned, CELL_CENTRE) - - myDiv_par(vn_aligned, CELL_CENTRE) - + Dn*LaplPerpnn/nn - + GradPerpDn_GradPerpnn/nn - + S/nn - + Snn/nn; - - ddt(nvn) += - vn*myGrad_par(nvn_aligned, CELL_YLOW) - - nvn*myDiv_par(vn_aligned, CELL_YLOW) - - Tn_stag*myGrad_par(nn_aligned, CELL_YLOW) - - vn*myInterp_to(div_vnperpnn_aligned, CELL_YLOW) - - nn_vperp_GradPerpvn - + myDiv_par(Dn_nn_Gradpar_vn_aligned, CELL_YLOW) - + Fi - + Fe/mu - + Snvn; - - if (munvn > 0.) ddt(nvn) += munvn*Delp2(nvn); + neutral_density_equation["parallel_advection"] = + - vn_centre*myGrad_par(lognn_aligned, CELL_CENTRE) + - myDiv_par(vn_aligned, CELL_CENTRE); + neutral_density_equation["perpendicular_diffusion"] = + Dn*LaplPerpnn/nn + GradPerpDn_GradPerpnn/nn; + neutral_density_equation["ionisation_recombination"] = S/nn; + neutral_density_equation["neutral_source"] = Snn/nn; + + neutral_momentum_equation["parallel_advection"] = + - vn*myGrad_par(nvn_aligned, CELL_YLOW) + - nvn*myDiv_par(vn_aligned, CELL_YLOW); + neutral_momentum_equation["density_gradient"] = + - Tn_stag*myGrad_par(nn_aligned, CELL_YLOW); + neutral_momentum_equation["perpendicular_diffusion"] = + - vn*myInterp_to(div_vnperpnn_aligned, CELL_YLOW) + - nn_vperp_GradPerpvn; + neutral_momentum_equation["parallel_diffusion"] = + myDiv_par(Dn_nn_Gradpar_vn_aligned, CELL_YLOW); + neutral_momentum_equation["ion_friction"] = Fi; + neutral_momentum_equation["electron_friction"] = Fe/mu; + neutral_momentum_equation["neutral_momentum_source"] = Snvn; + + if (munvn > 0.) { + neutral_momentum_equation["neutral_viscosity"] = munvn*Delp2(nvn); + } if (mesh->hasBndryLowerY()) { int j = mesh->ystart; @@ -494,14 +500,18 @@ void NeutralDVpar::recycleFluxes(BoutReal time) { if (bndry_neutrals == 2 || bndry_neutrals == 3) { // Sources to recycle density and momentum in the first cell + auto neutral_density_recycling = + neutral_density_equation["recycling_source"].localAccessor(); + auto neutral_momentum_recycling = + neutral_momentum_equation["recycling_source"].localAccessor(); j = mesh->ystart; for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { i = r.ind; BoutReal dd = 1./( sqrt(coordinates_stag->g_22(i,j)) * coordinates_stag->dy(i,j) ); for (int k = 0; k < mesh->LocalNz; k++) { BoutReal dnndt = recycled_lower(i,k)*dd; - ddt(lognn)(i,j,k) += dnndt/nn(i,j,k); - ddt(nvn)(i,j+1,k) += dnndt * vnth; + neutral_density_recycling(i,j,k) += dnndt/nn(i,j,k); + neutral_momentum_recycling(i,j+1,k) += dnndt * vnth; } } @@ -511,8 +521,8 @@ void NeutralDVpar::recycleFluxes(BoutReal time) { BoutReal dd = 1./( sqrt(coordinates_stag->g_22(i,j+1)) * coordinates_stag->dy(i,j+1) ); for (int k = 0; k < mesh->LocalNz; k++) { BoutReal dnndt = recycled_upper(i,k)*dd; - ddt(lognn)(i,j,k) += dnndt/nn(i,j,k); - ddt(nvn)(i,j,k) -= dnndt * vnth; + neutral_density_recycling(i,j,k) += dnndt/nn(i,j,k); + neutral_momentum_recycling(i,j,k) -= dnndt * vnth; } } } diff --git a/storm3d/D-vpar.hxx b/storm3d/D-vpar.hxx index d86e1df..429587e 100644 --- a/storm3d/D-vpar.hxx +++ b/storm3d/D-vpar.hxx @@ -7,12 +7,13 @@ #define __NEUTRAL_DVPAR_H__ #include "neutral-model.hxx" +#include "../shared/BoutEquation/equation.hxx" #include #include class NeutralDVpar : public NeutralModel { public: - NeutralDVpar(Solver *solver, Options &options); + NeutralDVpar(Solver *solver, Options &options, Datafile &dump); ~NeutralDVpar() {} /// Update plasma quantities @@ -27,6 +28,12 @@ public: void precon(BoutReal t, BoutReal gamma, BoutReal delta); private: + int rhs_counter = 0; + Equation neutral_density_equation{lognn, "lognn", Options::root()["save_equations"], + dump, rhs_counter}; + Equation neutral_momentum_equation{nvn, "nvn", Options::root()["save_equations"], dump, + rhs_counter}; + Field3D nn, nvn; Field3D Tn; Field3D lognn; diff --git a/storm3d/boundaries.cxx b/storm3d/boundaries.cxx index 5dbe30e..cc8e060 100644 --- a/storm3d/boundaries.cxx +++ b/storm3d/boundaries.cxx @@ -32,243 +32,241 @@ int i, j, k ; // Read and assign values here so that the default is taken as an explicitly set value, so // it won't be replaced by a different default later. void STORM::setBoundaryConditionsOptions() { - auto& opt = *globalOptions; - // density if (!average_radial_boundaries_core_SOL) { - opt["logn"]["bndry_xin"] = opt["logn"]["bndry_xin"].withDefault("neumann_o2"); - opt["logn"]["bndry_xout"] = opt["logn"]["bndry_xout"].withDefault("neumann_o2"); + globalOptions["logn"]["bndry_xin"] = globalOptions["logn"]["bndry_xin"].withDefault("neumann_o2"); + globalOptions["logn"]["bndry_xout"] = globalOptions["logn"]["bndry_xout"].withDefault("neumann_o2"); } else { - opt["logn"]["bndry_xin"] = opt["logn"]["bndry_xin"].withDefault("none"); - opt["logn"]["bndry_xout"] = opt["logn"]["bndry_xout"].withDefault("none"); + globalOptions["logn"]["bndry_xin"] = globalOptions["logn"]["bndry_xin"].withDefault("none"); + globalOptions["logn"]["bndry_xout"] = globalOptions["logn"]["bndry_xout"].withDefault("none"); } // y-boundary guard cells only used from n_aligned - opt["logn"]["bndry_ydown"] = opt["logn"]["bndry_ydown"].withDefault("none"); - opt["logn"]["bndry_yup"] = opt["logn"]["bndry_yup"].withDefault("none"); + globalOptions["logn"]["bndry_ydown"] = globalOptions["logn"]["bndry_ydown"].withDefault("none"); + globalOptions["logn"]["bndry_yup"] = globalOptions["logn"]["bndry_yup"].withDefault("none"); // x-boundary guard cells of n_aligned not used - opt["logn_aligned"]["bndry_xin"] = opt["logn_aligned"]["bndry_xin"].withDefault("none"); - opt["logn_aligned"]["bndry_xout"] = opt["logn_aligned"]["bndry_xout"].withDefault("none"); + globalOptions["logn_aligned"]["bndry_xin"] = globalOptions["logn_aligned"]["bndry_xin"].withDefault("none"); + globalOptions["logn_aligned"]["bndry_xout"] = globalOptions["logn_aligned"]["bndry_xout"].withDefault("none"); if (symmetry_plane) { - opt["logn_aligned"]["bndry_ydown"] = opt["logn_aligned"]["bndry_ydown"].withDefault("neumann_o2"); + globalOptions["logn_aligned"]["bndry_ydown"] = globalOptions["logn_aligned"]["bndry_ydown"].withDefault("neumann_o2"); } else { - opt["logn_aligned"]["bndry_ydown"] = opt["logn_aligned"]["bndry_ydown"].withDefault("free_o3"); + globalOptions["logn_aligned"]["bndry_ydown"] = globalOptions["logn_aligned"]["bndry_ydown"].withDefault("free_o3"); } - opt["logn_aligned"]["bndry_yup"] = opt["logn_aligned"]["bndry_yup"].withDefault("free_o3"); + globalOptions["logn_aligned"]["bndry_yup"] = globalOptions["logn_aligned"]["bndry_yup"].withDefault("free_o3"); // vorticity if (!average_radial_boundaries_core_SOL) { - opt["vort"]["bndry_xin"] = opt["vort"]["bndry_xin"].withDefault("neumann_o2"); - opt["vort"]["bndry_xout"] = opt["vort"]["bndry_xout"].withDefault("neumann_o2"); + globalOptions["vort"]["bndry_xin"] = globalOptions["vort"]["bndry_xin"].withDefault("neumann_o2"); + globalOptions["vort"]["bndry_xout"] = globalOptions["vort"]["bndry_xout"].withDefault("neumann_o2"); } else { - opt["vort"]["bndry_xin"] = opt["vort"]["bndry_xin"].withDefault("none"); - opt["vort"]["bndry_xout"] = opt["vort"]["bndry_xout"].withDefault("none"); + globalOptions["vort"]["bndry_xin"] = globalOptions["vort"]["bndry_xin"].withDefault("none"); + globalOptions["vort"]["bndry_xout"] = globalOptions["vort"]["bndry_xout"].withDefault("none"); } // y-boundary guard cells only used from vort_aligned - opt["vort"]["bndry_ydown"] = opt["vort"]["bndry_ydown"].withDefault("none"); - opt["vort"]["bndry_yup"] = opt["vort"]["bndry_yup"].withDefault("none"); + globalOptions["vort"]["bndry_ydown"] = globalOptions["vort"]["bndry_ydown"].withDefault("none"); + globalOptions["vort"]["bndry_yup"] = globalOptions["vort"]["bndry_yup"].withDefault("none"); // x-boundary guard cells of vort_aligned not used - opt["vort_aligned"]["bndry_xin"] = opt["vort_aligned"]["bndry_xin"].withDefault("none"); - opt["vort_aligned"]["bndry_xout"] = opt["vort_aligned"]["bndry_xout"].withDefault("none"); + globalOptions["vort_aligned"]["bndry_xin"] = globalOptions["vort_aligned"]["bndry_xin"].withDefault("none"); + globalOptions["vort_aligned"]["bndry_xout"] = globalOptions["vort_aligned"]["bndry_xout"].withDefault("none"); if (symmetry_plane) { - opt["vort_aligned"]["bndry_ydown"] = opt["vort_aligned"]["bndry_ydown"].withDefault("neumann_o2"); + globalOptions["vort_aligned"]["bndry_ydown"] = globalOptions["vort_aligned"]["bndry_ydown"].withDefault("neumann_o2"); } else { - opt["vort_aligned"]["bndry_ydown"] = opt["vort_aligned"]["bndry_ydown"].withDefault("free_o3"); + globalOptions["vort_aligned"]["bndry_ydown"] = globalOptions["vort_aligned"]["bndry_ydown"].withDefault("free_o3"); } - opt["vort_aligned"]["bndry_yup"] = opt["vort_aligned"]["bndry_yup"].withDefault("free_o3"); + globalOptions["vort_aligned"]["bndry_yup"] = globalOptions["vort_aligned"]["bndry_yup"].withDefault("free_o3"); // electrostatic potential // x-boundary guard cells set by Laplacian solver - opt["phi"]["bndry_xin"] = opt["phi"]["bndry_xin"].withDefault("none"); - opt["phi"]["bndry_xout"] = opt["phi"]["bndry_xout"].withDefault("none"); + globalOptions["phi"]["bndry_xin"] = globalOptions["phi"]["bndry_xin"].withDefault("none"); + globalOptions["phi"]["bndry_xout"] = globalOptions["phi"]["bndry_xout"].withDefault("none"); // y-boundary guard cells only used from phi_aligned - opt["phi"]["bndry_ydown"] = opt["phi"]["bndry_ydown"].withDefault("none"); - opt["phi"]["bndry_yup"] = opt["phi"]["bndry_yup"].withDefault("none"); + globalOptions["phi"]["bndry_ydown"] = globalOptions["phi"]["bndry_ydown"].withDefault("none"); + globalOptions["phi"]["bndry_yup"] = globalOptions["phi"]["bndry_yup"].withDefault("none"); // x-boundary guard cells of phi_aligned not used - opt["phi_aligned"]["bndry_xin"] = opt["phi_aligned"]["bndry_xin"].withDefault("none"); - opt["phi_aligned"]["bndry_xout"] = opt["phi_aligned"]["bndry_xout"].withDefault("none"); + globalOptions["phi_aligned"]["bndry_xin"] = globalOptions["phi_aligned"]["bndry_xin"].withDefault("none"); + globalOptions["phi_aligned"]["bndry_xout"] = globalOptions["phi_aligned"]["bndry_xout"].withDefault("none"); if (symmetry_plane) { - opt["phi_aligned"]["bndry_ydown"] = opt["phi_aligned"]["bndry_ydown"].withDefault("neumann_o2"); + globalOptions["phi_aligned"]["bndry_ydown"] = globalOptions["phi_aligned"]["bndry_ydown"].withDefault("neumann_o2"); } else { - opt["phi_aligned"]["bndry_ydown"] = opt["phi_aligned"]["bndry_ydown"].withDefault("free_o3"); + globalOptions["phi_aligned"]["bndry_ydown"] = globalOptions["phi_aligned"]["bndry_ydown"].withDefault("free_o3"); } - opt["phi_aligned"]["bndry_yup"] = opt["phi_aligned"]["bndry_yup"].withDefault("free_o3"); + globalOptions["phi_aligned"]["bndry_yup"] = globalOptions["phi_aligned"]["bndry_yup"].withDefault("free_o3"); if (run_1d) { // Not enough x-points to use free_o3 boundary conditions, so use neumann instead - opt["phi_stag"]["bndry_xin"] = opt["phi_stag"]["bndry_xin"].withDefault("neumann_o2"); - opt["phi_stag"]["bndry_xout"] = opt["phi_stag"]["bndry_xout"].withDefault("neumann_o2"); + globalOptions["phi_stag"]["bndry_xin"] = globalOptions["phi_stag"]["bndry_xin"].withDefault("neumann_o2"); + globalOptions["phi_stag"]["bndry_xout"] = globalOptions["phi_stag"]["bndry_xout"].withDefault("neumann_o2"); } else { - opt["phi_stag"]["bndry_xin"] = opt["phi_stag"]["bndry_xin"].withDefault("free_o3"); - opt["phi_stag"]["bndry_xout"] = opt["phi_stag"]["bndry_xout"].withDefault("free_o3"); + globalOptions["phi_stag"]["bndry_xin"] = globalOptions["phi_stag"]["bndry_xin"].withDefault("free_o3"); + globalOptions["phi_stag"]["bndry_xout"] = globalOptions["phi_stag"]["bndry_xout"].withDefault("free_o3"); } // y-boundary guard cells not used - opt["phi_stag"]["bndry_ydown"] = opt["phi_stag"]["bndry_ydown"].withDefault("none"); - opt["phi_stag"]["bndry_yup"] = opt["phi_stag"]["bndry_yup"].withDefault("none"); + globalOptions["phi_stag"]["bndry_ydown"] = globalOptions["phi_stag"]["bndry_ydown"].withDefault("none"); + globalOptions["phi_stag"]["bndry_yup"] = globalOptions["phi_stag"]["bndry_yup"].withDefault("none"); // electron temperature if (!isothermal) { if (!average_radial_boundaries_core_SOL) { - opt["logp"]["bndry_xin"] = opt["logp"]["bndry_xin"].withDefault("neumann_o2"); - opt["logp"]["bndry_xout"] = opt["logp"]["bndry_xout"].withDefault("neumann_o2"); + globalOptions["logp"]["bndry_xin"] = globalOptions["logp"]["bndry_xin"].withDefault("neumann_o2"); + globalOptions["logp"]["bndry_xout"] = globalOptions["logp"]["bndry_xout"].withDefault("neumann_o2"); } else { - opt["logp"]["bndry_xin"] = opt["logp"]["bndry_xin"].withDefault("none"); - opt["logp"]["bndry_xout"] = opt["logp"]["bndry_xout"].withDefault("none"); + globalOptions["logp"]["bndry_xin"] = globalOptions["logp"]["bndry_xin"].withDefault("none"); + globalOptions["logp"]["bndry_xout"] = globalOptions["logp"]["bndry_xout"].withDefault("none"); } // y-boundary guard cells only used from T_aligned - opt["logp"]["bndry_ydown"] = opt["logp"]["bndry_ydown"].withDefault("none"); - opt["logp"]["bndry_yup"] = opt["logp"]["bndry_yup"].withDefault("none"); + globalOptions["logp"]["bndry_ydown"] = globalOptions["logp"]["bndry_ydown"].withDefault("none"); + globalOptions["logp"]["bndry_yup"] = globalOptions["logp"]["bndry_yup"].withDefault("none"); // x-boundary guard cells of T_aligned not used - opt["logT_aligned"]["bndry_xin"] = opt["logT_aligned"]["bndry_xin"].withDefault("none"); - opt["logT_aligned"]["bndry_xout"] = opt["logT_aligned"]["bndry_xout"].withDefault("none"); + globalOptions["logT_aligned"]["bndry_xin"] = globalOptions["logT_aligned"]["bndry_xin"].withDefault("none"); + globalOptions["logT_aligned"]["bndry_xout"] = globalOptions["logT_aligned"]["bndry_xout"].withDefault("none"); if (symmetry_plane) { - opt["logT_aligned"]["bndry_ydown"] = opt["logT_aligned"]["bndry_ydown"].withDefault("neumann_o2"); + globalOptions["logT_aligned"]["bndry_ydown"] = globalOptions["logT_aligned"]["bndry_ydown"].withDefault("neumann_o2"); } else { - opt["logT_aligned"]["bndry_ydown"] = opt["logT_aligned"]["bndry_ydown"].withDefault("free_o3"); + globalOptions["logT_aligned"]["bndry_ydown"] = globalOptions["logT_aligned"]["bndry_ydown"].withDefault("free_o3"); } - opt["logT_aligned"]["bndry_yup"] = opt["logT_aligned"]["bndry_yup"].withDefault("free_o3"); + globalOptions["logT_aligned"]["bndry_yup"] = globalOptions["logT_aligned"]["bndry_yup"].withDefault("free_o3"); } // ion velocity if (!average_radial_boundaries_core_SOL) { - opt["U"]["bndry_xin"] = opt["U"]["bndry_xin"].withDefault("neumann_o2"); - opt["U"]["bndry_xout"] = opt["U"]["bndry_xout"].withDefault("neumann_o2"); + globalOptions["U"]["bndry_xin"] = globalOptions["U"]["bndry_xin"].withDefault("neumann_o2"); + globalOptions["U"]["bndry_xout"] = globalOptions["U"]["bndry_xout"].withDefault("neumann_o2"); } else { - opt["U"]["bndry_xin"] = opt["U"]["bndry_xin"].withDefault("none"); - opt["U"]["bndry_xout"] = opt["U"]["bndry_xout"].withDefault("none"); + globalOptions["U"]["bndry_xin"] = globalOptions["U"]["bndry_xin"].withDefault("none"); + globalOptions["U"]["bndry_xout"] = globalOptions["U"]["bndry_xout"].withDefault("none"); } // y-boundary guard cells only used from U_aligned - opt["U"]["bndry_ydown"] = opt["U"]["bndry_ydown"].withDefault("none"); - opt["U"]["bndry_yup"] = opt["U"]["bndry_yup"].withDefault("none"); + globalOptions["U"]["bndry_ydown"] = globalOptions["U"]["bndry_ydown"].withDefault("none"); + globalOptions["U"]["bndry_yup"] = globalOptions["U"]["bndry_yup"].withDefault("none"); // x-boundary guard cells of U_aligned not used - opt["U_aligned"]["bndry_xin"] = opt["U_aligned"]["bndry_xin"].withDefault("none"); - opt["U_aligned"]["bndry_xout"] = opt["U_aligned"]["bndry_xout"].withDefault("none"); + globalOptions["U_aligned"]["bndry_xin"] = globalOptions["U_aligned"]["bndry_xin"].withDefault("none"); + globalOptions["U_aligned"]["bndry_xout"] = globalOptions["U_aligned"]["bndry_xout"].withDefault("none"); if (symmetry_plane) { - opt["U_aligned"]["bndry_ydown"] = opt["U_aligned"]["bndry_ydown"].withDefault("dirichlet_o2"); + globalOptions["U_aligned"]["bndry_ydown"] = globalOptions["U_aligned"]["bndry_ydown"].withDefault("dirichlet_o2"); } else { // y-boundary guard cells set by sheath boundary conditions function - opt["U_aligned"]["bndry_ydown"] = opt["U_aligned"]["bndry_ydown"].withDefault("none"); + globalOptions["U_aligned"]["bndry_ydown"] = globalOptions["U_aligned"]["bndry_ydown"].withDefault("none"); } // y-boundary guard cells set by sheath boundary conditions function - opt["U_aligned"]["bndry_yup"] = opt["U_aligned"]["bndry_yup"].withDefault("none"); + globalOptions["U_aligned"]["bndry_yup"] = globalOptions["U_aligned"]["bndry_yup"].withDefault("none"); // electron velocity if (!average_radial_boundaries_core_SOL) { - opt["V"]["bndry_xin"] = opt["V"]["bndry_xin"].withDefault("neumann_o2"); - opt["V"]["bndry_xout"] = opt["V"]["bndry_xout"].withDefault("neumann_o2"); + globalOptions["V"]["bndry_xin"] = globalOptions["V"]["bndry_xin"].withDefault("neumann_o2"); + globalOptions["V"]["bndry_xout"] = globalOptions["V"]["bndry_xout"].withDefault("neumann_o2"); } else { - opt["V"]["bndry_xin"] = opt["V"]["bndry_xin"].withDefault("none"); - opt["V"]["bndry_xout"] = opt["V"]["bndry_xout"].withDefault("none"); + globalOptions["V"]["bndry_xin"] = globalOptions["V"]["bndry_xin"].withDefault("none"); + globalOptions["V"]["bndry_xout"] = globalOptions["V"]["bndry_xout"].withDefault("none"); } // y-boundary guard cells only used from V_aligned - opt["V"]["bndry_ydown"] = opt["V"]["bndry_ydown"].withDefault("none"); - opt["V"]["bndry_yup"] = opt["V"]["bndry_yup"].withDefault("none"); + globalOptions["V"]["bndry_ydown"] = globalOptions["V"]["bndry_ydown"].withDefault("none"); + globalOptions["V"]["bndry_yup"] = globalOptions["V"]["bndry_yup"].withDefault("none"); // x-boundary guard cells of V_aligned not used - opt["V_aligned"]["bndry_xin"] = opt["V_aligned"]["bndry_xin"].withDefault("none"); - opt["V_aligned"]["bndry_xout"] = opt["V_aligned"]["bndry_xout"].withDefault("none"); + globalOptions["V_aligned"]["bndry_xin"] = globalOptions["V_aligned"]["bndry_xin"].withDefault("none"); + globalOptions["V_aligned"]["bndry_xout"] = globalOptions["V_aligned"]["bndry_xout"].withDefault("none"); if (symmetry_plane) { - opt["V_aligned"]["bndry_ydown"] = opt["V_aligned"]["bndry_ydown"].withDefault("dirichlet_o2"); + globalOptions["V_aligned"]["bndry_ydown"] = globalOptions["V_aligned"]["bndry_ydown"].withDefault("dirichlet_o2"); } else { // y-boundary guard cells set by sheath boundary conditions function - opt["V_aligned"]["bndry_ydown"] = opt["V_aligned"]["bndry_ydown"].withDefault("none"); + globalOptions["V_aligned"]["bndry_ydown"] = globalOptions["V_aligned"]["bndry_ydown"].withDefault("none"); } // y-boundary guard cells set by sheath boundary conditions function - opt["V_aligned"]["bndry_yup"] = opt["V_aligned"]["bndry_yup"].withDefault("none"); + globalOptions["V_aligned"]["bndry_yup"] = globalOptions["V_aligned"]["bndry_yup"].withDefault("none"); // chiU and chiV do not need boundary conditions - we never take derivatives // of either of them, only of psi, U and V - opt["chiU"]["bndry_xin"] = opt["chiU"]["bndry_xin"].withDefault("none"); - opt["chiU"]["bndry_xout"] = opt["chiU"]["bndry_xout"].withDefault("none"); - opt["chiU"]["bndry_ydown"] = opt["chiU"]["bndry_ydown"].withDefault("none"); - opt["chiU"]["bndry_yup"] = opt["chiU"]["bndry_yup"].withDefault("none"); - opt["chiV"]["bndry_xin"] = opt["chiV"]["bndry_xin"].withDefault("none"); - opt["chiV"]["bndry_xout"] = opt["chiV"]["bndry_xout"].withDefault("none"); - opt["chiV"]["bndry_ydown"] = opt["chiV"]["bndry_ydown"].withDefault("none"); - opt["chiV"]["bndry_yup"] = opt["chiV"]["bndry_yup"].withDefault("none"); + globalOptions["chiU"]["bndry_xin"] = globalOptions["chiU"]["bndry_xin"].withDefault("none"); + globalOptions["chiU"]["bndry_xout"] = globalOptions["chiU"]["bndry_xout"].withDefault("none"); + globalOptions["chiU"]["bndry_ydown"] = globalOptions["chiU"]["bndry_ydown"].withDefault("none"); + globalOptions["chiU"]["bndry_yup"] = globalOptions["chiU"]["bndry_yup"].withDefault("none"); + globalOptions["chiV"]["bndry_xin"] = globalOptions["chiV"]["bndry_xin"].withDefault("none"); + globalOptions["chiV"]["bndry_xout"] = globalOptions["chiV"]["bndry_xout"].withDefault("none"); + globalOptions["chiV"]["bndry_ydown"] = globalOptions["chiV"]["bndry_ydown"].withDefault("none"); + globalOptions["chiV"]["bndry_yup"] = globalOptions["chiV"]["bndry_yup"].withDefault("none"); // electron parallel heat flux // x-boundaries are not needed if (!isothermal) { - opt["qpar_aligned"]["bndry_xin"] = opt["qpar_aligned"]["bndry_xin"].withDefault("none"); - opt["qpar_aligned"]["bndry_xout"] = opt["qpar_aligned"]["bndry_xout"].withDefault("none"); + globalOptions["qpar_aligned"]["bndry_xin"] = globalOptions["qpar_aligned"]["bndry_xin"].withDefault("none"); + globalOptions["qpar_aligned"]["bndry_xout"] = globalOptions["qpar_aligned"]["bndry_xout"].withDefault("none"); if (symmetry_plane) { - opt["qpar_aligned"]["bndry_ydown"] = opt["qpar_aligned"]["bndry_ydown"].withDefault("dirichlet_o2"); + globalOptions["qpar_aligned"]["bndry_ydown"] = globalOptions["qpar_aligned"]["bndry_ydown"].withDefault("dirichlet_o2"); } else { // y-boundary guard cells set by sheath boundary conditions function - opt["qpar_aligned"]["bndry_ydown"] = opt["qpar_aligned"]["bndry_ydown"].withDefault("none"); + globalOptions["qpar_aligned"]["bndry_ydown"] = globalOptions["qpar_aligned"]["bndry_ydown"].withDefault("none"); } // y-boundary guard cells set by sheath boundary conditions function - opt["qpar_aligned"]["bndry_yup"] = opt["qpar_aligned"]["bndry_yup"].withDefault("none"); + globalOptions["qpar_aligned"]["bndry_yup"] = globalOptions["qpar_aligned"]["bndry_yup"].withDefault("none"); } if (electromagnetic && !isothermal) { - opt["qpar_centre"]["bndry_xin"] = opt["qpar_centre"]["bndry_xin"].withDefault("neumann_o2"); - opt["qpar_centre"]["bndry_xout"] = opt["qpar_centre"]["bndry_xout"].withDefault("neumann_o2"); + globalOptions["qpar_centre"]["bndry_xin"] = globalOptions["qpar_centre"]["bndry_xin"].withDefault("neumann_o2"); + globalOptions["qpar_centre"]["bndry_xout"] = globalOptions["qpar_centre"]["bndry_xout"].withDefault("neumann_o2"); // y-boundaries not needed - opt["qpar_centre"]["bndry_ydown"] = opt["qpar_centre"]["bndry_ydown"].withDefault("none"); - opt["qpar_centre"]["bndry_yup"] = opt["qpar_centre"]["bndry_yup"].withDefault("none"); + globalOptions["qpar_centre"]["bndry_ydown"] = globalOptions["qpar_centre"]["bndry_ydown"].withDefault("none"); + globalOptions["qpar_centre"]["bndry_yup"] = globalOptions["qpar_centre"]["bndry_yup"].withDefault("none"); } // ExB speed squared if (boussinesq == 0) { if (run_1d) { // Not enough x-points to use free_o3 boundary conditions, so use neumann instead - opt["uE2"]["bndry_xin"] = opt["uE2"]["bndry_xin"].withDefault("neumann_o2"); - opt["uE2"]["bndry_xout"] = opt["uE2"]["bndry_xout"].withDefault("neumann_o2"); + globalOptions["uE2"]["bndry_xin"] = globalOptions["uE2"]["bndry_xin"].withDefault("neumann_o2"); + globalOptions["uE2"]["bndry_xout"] = globalOptions["uE2"]["bndry_xout"].withDefault("neumann_o2"); } else { - opt["uE2"]["bndry_xin"] = opt["uE2"]["bndry_xin"].withDefault("free_o3"); - opt["uE2"]["bndry_xout"] = opt["uE2"]["bndry_xout"].withDefault("free_o3"); + globalOptions["uE2"]["bndry_xin"] = globalOptions["uE2"]["bndry_xin"].withDefault("free_o3"); + globalOptions["uE2"]["bndry_xout"] = globalOptions["uE2"]["bndry_xout"].withDefault("free_o3"); } - opt["uE2"]["bndry_ydown"] = opt["uE2"]["bndry_ydown"].withDefault("none"); - opt["uE2"]["bndry_yup"] = opt["uE2"]["bndry_yup"].withDefault("none"); + globalOptions["uE2"]["bndry_ydown"] = globalOptions["uE2"]["bndry_ydown"].withDefault("none"); + globalOptions["uE2"]["bndry_yup"] = globalOptions["uE2"]["bndry_yup"].withDefault("none"); } // parallel component of magnetic vector potential if (electromagnetic) { // x-boundary guard cells set by Laplacian solver - opt["psi"]["bndry_xin"] = opt["psi"]["bndry_xin"].withDefault("none"); - opt["psi"]["bndry_xout"] = opt["psi"]["bndry_xout"].withDefault("none"); + globalOptions["psi"]["bndry_xin"] = globalOptions["psi"]["bndry_xin"].withDefault("none"); + globalOptions["psi"]["bndry_xout"] = globalOptions["psi"]["bndry_xout"].withDefault("none"); // y-boundary guard cells only used from psi_aligned - opt["psi"]["bndry_ydown"] = opt["psi"]["bndry_ydown"].withDefault("none"); - opt["psi"]["bndry_yup"] = opt["psi"]["bndry_yup"].withDefault("none"); + globalOptions["psi"]["bndry_ydown"] = globalOptions["psi"]["bndry_ydown"].withDefault("none"); + globalOptions["psi"]["bndry_yup"] = globalOptions["psi"]["bndry_yup"].withDefault("none"); // x-boundary guard cells of psi_aligned not used - opt["psi_aligned"]["bndry_xin"] = opt["psi_aligned"]["bndry_xin"].withDefault("none"); - opt["psi_aligned"]["bndry_xout"] = opt["psi_aligned"]["bndry_xout"].withDefault("none"); + globalOptions["psi_aligned"]["bndry_xin"] = globalOptions["psi_aligned"]["bndry_xin"].withDefault("none"); + globalOptions["psi_aligned"]["bndry_xout"] = globalOptions["psi_aligned"]["bndry_xout"].withDefault("none"); if (symmetry_plane) { - opt["psi_aligned"]["bndry_ydown"] = - opt["psi_aligned"]["bndry_ydown"].withDefault("dirichlet_o2"); + globalOptions["psi_aligned"]["bndry_ydown"] = + globalOptions["psi_aligned"]["bndry_ydown"].withDefault("dirichlet_o2"); } else { if (!use_psi_boundary_solver) { - opt["psi_aligned"]["bndry_ydown"] = - opt["psi_aligned"]["bndry_ydown"].withDefault("free_o3"); + globalOptions["psi_aligned"]["bndry_ydown"] = + globalOptions["psi_aligned"]["bndry_ydown"].withDefault("free_o3"); } else { // y-boundary guard cells set by sheath boundary conditions function - opt["psi_aligned"]["bndry_ydown"] = - opt["psi_aligned"]["bndry_ydown"].withDefault("none"); + globalOptions["psi_aligned"]["bndry_ydown"] = + globalOptions["psi_aligned"]["bndry_ydown"].withDefault("none"); } } if (!use_psi_boundary_solver) { - opt["psi_aligned"]["bndry_yup"] = - opt["psi_aligned"]["bndry_yup"].withDefault("free_o3"); + globalOptions["psi_aligned"]["bndry_yup"] = + globalOptions["psi_aligned"]["bndry_yup"].withDefault("free_o3"); } else { // y-boundary guard cells set by sheath boundary conditions function - opt["psi_aligned"]["bndry_yup"] = - opt["psi_aligned"]["bndry_yup"].withDefault("none"); + globalOptions["psi_aligned"]["bndry_yup"] = + globalOptions["psi_aligned"]["bndry_yup"].withDefault("none"); } - opt["psi_centre"]["bndry_xin"] = - opt["psi_centre"]["bndry_xin"].withDefault("dirichlet_o2"); - opt["psi_centre"]["bndry_xout"] = - opt["psi_centre"]["bndry_xout"].withDefault("dirichlet_o2"); + globalOptions["psi_centre"]["bndry_xin"] = + globalOptions["psi_centre"]["bndry_xin"].withDefault("dirichlet_o2"); + globalOptions["psi_centre"]["bndry_xout"] = + globalOptions["psi_centre"]["bndry_xout"].withDefault("dirichlet_o2"); // y-boundary guard cells not needed - opt["psi_centre"]["bndry_ydown"] = - opt["psi_centre"]["bndry_ydown"].withDefault("none"); - opt["psi_centre"]["bndry_yup"] = opt["psi_centre"]["bndry_yup"].withDefault("none"); + globalOptions["psi_centre"]["bndry_ydown"] = + globalOptions["psi_centre"]["bndry_ydown"].withDefault("none"); + globalOptions["psi_centre"]["bndry_yup"] = globalOptions["psi_centre"]["bndry_yup"].withDefault("none"); } } @@ -493,6 +491,9 @@ FieldPerp STORM::extrap_sheath_upper(const Field3D &var){ ASSERT1(var.getDirectionY() == YDirectionType::Aligned); ASSERT1(var.getLocation() == CELL_CENTRE); + // Should only be called when there is an upper sheath boundary + ASSERT1(mesh->hasBndryUpperY()); + //3rd order extrapolation of a field to find the sheath value FieldPerp result; result = 0.375*sliceXZ(var, mesh->yend + 1) + 0.75*sliceXZ(var, mesh->yend) - 0.125*sliceXZ(var, mesh->yend - 1) ; @@ -506,6 +507,9 @@ FieldPerp STORM::extrap_sheath_lower(const Field3D &var){ ASSERT1(var.getDirectionY() == YDirectionType::Aligned); ASSERT1(var.getLocation() == CELL_CENTRE); + // Should only be called when there is a lower sheath boundary + ASSERT1(mesh->hasBndryLowerY()); + //3rd order extrapolation of a field to find the sheath value FieldPerp result; result = 0.375*sliceXZ(var, mesh->ystart - 1) + 0.75*sliceXZ(var, mesh->ystart) - 0.125*sliceXZ(var, mesh->ystart + 1) ; @@ -523,17 +527,22 @@ FieldPerp STORM::extrap_sheath_lower(const Field3D &var){ void STORM::phi_bc_initialise(bool restarting){ bool append; OPTION(globalOptions, append, false); - if(restarting && append){ + OPTION(options, reset_evolving_bcs, false); + // replicate setting of dump_on_restart from BOUT++'s Solver::solve() + bool dump_on_restart; + OPTION(globalOptions, dump_on_restart, (not restarting) or (not append)); + if ((not dump_on_restart) or reset_evolving_bcs) { + // Note: always need to increment cstep_SOL, etc. on first step when + // reset_evolving_bcs=true to avoid divide-by-zero first_step = false; }else{ first_step = true; } // If monitor_timestep == false, the evolving boundary conditions are not working - Options* solver_options = globalOptions->getSection("solver"); - bool monitor_timestep; - solver_options->get("monitor_timestep", monitor_timestep, false); - if(!monitor_timestep) throw BoutException("Evolving boundary phi: monitor_timstep = false."); + Options& solver_options = *globalOptions.getSection("solver"); + bool monitor_timestep = solver_options["monitor_timestep"].withDefault(false); + if(!monitor_timestep) throw BoutException("Evolving boundary phi: monitor_timestep = false."); restart.add(phi_bc,"phi_bc",0); restart.add(time_last_SOL,"time_last_SOL",0); @@ -542,18 +551,21 @@ void STORM::phi_bc_initialise(bool restarting){ restart.add(cstep_SOL,"cstep_SOL",0); restart.add(cstep_PF,"cstep_PF",0); restart.add(cstep_core,"cstep_core",0); - if(!restarting){ + SAVE_REPEAT(phi_bc, time_last_SOL, time_last_PF, time_last_core, cstep_SOL, + cstep_PF, cstep_core); + if ((not restarting) or reset_evolving_bcs) { //Start the simulation from scratch - phi_bc = 0.; + phi_bc = 0.0; time_last_SOL = 0.; time_last_PF = 0.; time_last_core = 0.; cstep_SOL = 0.; cstep_PF = 0.; cstep_core = 0.; + auto phi_dc = DC(phi); for(int i = 0; iLocalNy; ++i){ - phi_bc(mesh->xstart,i) = phi(mesh->xstart,i,0); - phi_bc(mesh->xend,i) = phi(mesh->xend,i,0); + phi_bc(mesh->xstart,i) = phi_dc(mesh->xstart,i); + phi_bc(mesh->xend,i) = phi_dc(mesh->xend,i); } } } diff --git a/storm3d/data/BOUT.inp b/storm3d/data/BOUT.inp index 8c07f99..85a5876 100755 --- a/storm3d/data/BOUT.inp +++ b/storm3d/data/BOUT.inp @@ -5,20 +5,13 @@ MZ = 128 # Number of Z points zmin = 0 zmax = 15.9155 # 1.5915*2pi = 10, z is fracs of 2pi -[mesh:ddx] -first = C2 -second = C2 -upwind = C2 [mesh:ddy] -first = C2 -second = C2 upwind = U2 [mesh:ddz] first = C2 second = C2 -upwind = C2 [mesh] StaggerGrids = true # Enables staggered grids @@ -32,12 +25,6 @@ dy = Ly/ny ixseps1 = -1 # Set x location of separatrix 1 ixseps2 = 0 # Set x location of separatrix 2 -[laplace] -global_flags = 0 -inner_boundary_flags = 16 -outer_boundary_flags = 16 -include_yguards = false - [solver] mxstep = 100000000 # max steps before result is deemed not to converge @@ -50,7 +37,6 @@ q = 7 # Dimensionless R_c = 1.5 # m n_0 = 0.8e19 # m^-3 loglambda = -1 # Dimensionless -phi_wall = 0 # If these parameters are specified, they will be used instead of the values calculated from the primary parameters above. # mu_n0 = 0.01 @@ -66,7 +52,6 @@ normalise_lengths = false bracket = 2 # 0 = std, 1 = simple, 2 = arakawa add_blob = true symmetry_plane = true -run_1d = false [blob] delta_perp = 10.0 @@ -111,4 +96,4 @@ function = exp(-10)*exp( 10*(y/(2.0*pi))) #if symmetry_plane=true function = exp(-10)*exp( 10*(y/(2.0*pi))) #if symmetry_plane=true [phi] -function = (storm:phi_wall) + 3.18 +function = 3.18 diff --git a/storm3d/initialise.cxx b/storm3d/initialise.cxx index d094d4c..cfa3174 100644 --- a/storm3d/initialise.cxx +++ b/storm3d/initialise.cxx @@ -417,7 +417,7 @@ void STORM::initialise_blob(const char * imp_section){ BoutReal delta_front = -1.; // Parameter used to control gradient of the density blob front. 0 = step function BoutReal delta_front_T = -1.; // Parameter used to control gradient of the temperature blob front. 0 = step function bool A_relative_to_bg; // Switch for controlling whether the amplitude of the blob is relative to the midplane density or not. - Options *blob_options = Options::getRoot()->getSection(imp_section); + Options &blob_options = *globalOptions.getSection(imp_section); OPTION(blob_options, delta_perp, 10) ; OPTION(blob_options, elongation, 1) ; @@ -571,17 +571,16 @@ void STORM::add_noise(){ // Add a function given in the input file to a field void STORM::add_perturbation(Field3D &f, const std::string name) { - Options *options = Options::getRoot()->getSection("filaments"); - Field3D perturbation = FieldFactory::get()->create3D(name, options, mesh, f.getLocation(), 0); + Options &options = *globalOptions.getSection("storm"); + Field3D perturbation = FieldFactory::get()->create3D(name, &options, mesh, f.getLocation(), 0); f += perturbation; } // For consistency, U and V should have the same radial boundary conditions, as // V's boundary conditions are applied to UmV_centre void STORM::check_U_V_x_boundary_conditions() { - auto options = *globalOptions; - ASSERT0(options["U"]["bndry_xin"].as() - == options["V"]["bndry_xin"].as()); - ASSERT0(options["U"]["bndry_xout"].as() - == options["V"]["bndry_xout"].as()); + ASSERT0(globalOptions["U"]["bndry_xin"].as() + == globalOptions["V"]["bndry_xin"].as()); + ASSERT0(globalOptions["U"]["bndry_xout"].as() + == globalOptions["V"]["bndry_xout"].as()); } diff --git a/storm3d/makefile b/storm3d/makefile index 61c30ba..fa8acfd 100644 --- a/storm3d/makefile +++ b/storm3d/makefile @@ -24,10 +24,20 @@ else $(error "The 'make.config' file does not exist. You must create it for STORM to compile. See the example in 'make.config.example'") endif -SOURCEC = storm.cxx boundaries.cxx initialise.cxx monitors.cxx operators.cxx utilities.cxx ../shared/fast_output.cxx neutral-rates.cxx D-vpar.cxx neutral-model.cxx +SOURCEC = storm.cxx boundaries.cxx initialise.cxx monitors.cxx operators.cxx utilities.cxx ../shared/fast_output.cxx ../shared/BoutEquation/equation.cxx neutral-rates.cxx D-vpar.cxx neutral-model.cxx TARGET = storm -GIT_VERSION := $(shell git describe --abbrev=40 --dirty --always --tags) -CXXFLAGS += -DSTORM_VERSION=\"$(GIT_VERSION)\" +save_git_version_status := $(shell ../build_tools/save_git_version.py) +ifeq ($(save_git_version_status), 11) + $(error "gitpython not installed. You can install with 'conda install gitpython' or 'pip3 install --user gitpython'") +endif +ifeq ($(save_git_version_status), 12) + $(error "save_git_version.py failed to find a delimiter that is not in the git diff") +endif +ifneq ($(save_git_version_status), 0) + $(error "save_git_version.py failed with an unrecognised error $(save_git_version_status)") +endif + +$(shell [ -e storm.o ] && rm storm.o) include $(BOUT_TOP)/make.config diff --git a/storm3d/neutral-model.cxx b/storm3d/neutral-model.cxx index cbdfcba..b6ebf1a 100644 --- a/storm3d/neutral-model.cxx +++ b/storm3d/neutral-model.cxx @@ -5,7 +5,7 @@ #include "neutral-model.hxx" #include "D-vpar.hxx" -NeutralModel *NeutralModel::create(Solver *solver, Options &options) { +NeutralModel *NeutralModel::create(Solver *solver, Options &options, Datafile &dump) { // Decide which neutral model to use std::string type = options["type"].withDefault("none"); @@ -14,7 +14,7 @@ NeutralModel *NeutralModel::create(Solver *solver, Options &options) { return NULL; } else if (type == "d-vpar") { // Diffusive in X-Z, fluid in Y - return new NeutralDVpar(solver, options); + return new NeutralDVpar(solver, options, dump); } throw BoutException("Unrecognised neutral model '%s'", type.c_str()); } diff --git a/storm3d/neutral-model.hxx b/storm3d/neutral-model.hxx index f2dd427..9cac133 100644 --- a/storm3d/neutral-model.hxx +++ b/storm3d/neutral-model.hxx @@ -11,7 +11,7 @@ class NeutralModel { public: - NeutralModel(Options &options) { + NeutralModel(Options &options, Datafile &dump_in) : dump(dump_in) { // Better to initialize with NaN S = 0.; @@ -58,7 +58,7 @@ public: /*! * Creates an instance of NeutralModel, based on given options */ - static NeutralModel* create(Solver *solver, Options &options); + static NeutralModel* create(Solver *solver, Options &options, Datafile &dump); /*! * Set normalisations for temperature [eV], density [m^-3], @@ -118,6 +118,8 @@ protected: const Field3D& nn, const Field3D& nn_stag, const Field3D& vn, bool updaterates = true); + Datafile& dump; + private: NeutralModel(); }; diff --git a/storm3d/storm.cxx b/storm3d/storm.cxx index 382a3da..06b27e9 100644 --- a/storm3d/storm.cxx +++ b/storm3d/storm.cxx @@ -35,9 +35,11 @@ class STORM; BOUTMAIN(STORM); int STORM::init(bool restarting) { + + // Add version information to output files + dump.setAttribute("", "storm_git_hash", storm_git_hash); + dump.setAttribute("", "storm_git_diff", storm_git_diff); - globalOptions = Options::getRoot(); - options = Options::getRoot()->getSection("storm"); coordinates_centre = mesh->getCoordinates(CELL_CENTRE); OPTION(options, mu_n0, -1.0) ; @@ -86,8 +88,7 @@ int STORM::init(bool restarting) { // The user has not explicitly set equilibrium_file_path, use path relative to simulation output directory instead std::string equilibrium_directory; OPTION(options, equilibrium_directory, "equilibrium"); - std::string data_dir; - Options::getRoot()->get("datadir", data_dir, "data"); + std::string data_dir = globalOptions["datadir"].withDefault("data"); equilibrium_file_path = data_dir + "/" + equilibrium_directory; } OPTION(options, equilibrium_data_file, ""); @@ -100,7 +101,7 @@ int STORM::init(bool restarting) { OPTION(options, sources_realisticgeometry_background, false) ; OPTION(options, increased_dissipation_xbndries, false) ; OPTION(options, increased_resistivity_core, false) ; - OPTION(options, normalise_sources, true) ; + OPTION(options, normalise_sources, realistic_geometry == "none") ; // Set default values for boundary conditions for all variables. // Can be overridden in input file, but should never need to be. @@ -108,7 +109,7 @@ int STORM::init(bool restarting) { setBoundaryConditionsOptions(); // Check if we're using ShiftedMetric to run simulation in x-z orthogonal coordinates - auto& mesh_options = (*globalOptions)["mesh"]; + auto& mesh_options = globalOptions["mesh"]; isshifted = false; if (mesh_options["paralleltransform"] == "shifted") { // As we use staggered grids, we must transform to/from field-aligned coordinates @@ -268,49 +269,52 @@ int STORM::init(bool restarting) { } // Normalise grid mesh. + coordinates_stag = mesh->getCoordinates(CELL_YLOW); + if (normalise_lengths){ - coordinates_centre->g_11 /= SQ(rho_s) ; - coordinates_centre->g_22 /= SQ(rho_s) ; - coordinates_centre->g_33 /= SQ(rho_s) ; - coordinates_centre->g_12 /= SQ(rho_s) ; - coordinates_centre->g_13 /= SQ(rho_s) ; - coordinates_centre->g_23 /= SQ(rho_s) ; - coordinates_centre->g11 *= SQ(rho_s) ; - coordinates_centre->g22 *= SQ(rho_s) ; - coordinates_centre->g33 *= SQ(rho_s) ; - coordinates_centre->g12 *= SQ(rho_s) ; - coordinates_centre->g13 *= SQ(rho_s) ; - coordinates_centre->g23 *= SQ(rho_s) ; - //mesh->dx /= rho_s ; - //mesh->dz /= rho_s ; - //mesh->dy /= rho_s ; - //mesh->zlength/= rho_s ; - coordinates_centre->geometry(); + for (auto* coords : {coordinates_centre, coordinates_stag}) { + coords->g_11 /= SQ(rho_s); + coords->g_22 /= SQ(rho_s); + coords->g_33 /= SQ(rho_s); + coords->g_12 /= SQ(rho_s); + coords->g_13 /= SQ(rho_s); + coords->g_23 /= SQ(rho_s); + coords->g11 *= SQ(rho_s); + coords->g22 *= SQ(rho_s); + coords->g33 *= SQ(rho_s); + coords->g12 *= SQ(rho_s); + coords->g13 *= SQ(rho_s); + coords->g23 *= SQ(rho_s); + coords->geometry(); + } if (realistic_geometry != "none") { G3 = coordinates_centre->G3; } } if (normalise_all) { - if (normalise_lengths) + if (normalise_lengths) { throw BoutException("Error: not possible to use normalise_lengths = true and normalise_all = true."); - coordinates_centre->Bxy /= B_0; + } + for (auto* coords : {coordinates_centre, coordinates_stag}) { + coords->Bxy /= B_0; + coords->g_11 *= SQ(rho_s*B_0); + coords->g_22 /= SQ(rho_s); + coords->g_33 /= SQ(rho_s); + coords->g_12 *= B_0; + coords->g_13 *= B_0; + coords->g_23 /= SQ(rho_s); + coords->g11 /= SQ(rho_s*B_0); + coords->g22 *= SQ(rho_s); + coords->g33 *= SQ(rho_s); + coords->g12 /= B_0; + coords->g13 /= B_0; + coords->g23 *= SQ(rho_s); + coords->dx /= SQ(rho_s)*B_0; + coords->J *= B_0/rho_s; + coords->geometry(); + } B2 = SQ(coordinates_centre->Bxy); - coordinates_centre->g_11 *= SQ(rho_s*B_0) ; - coordinates_centre->g_22 /= SQ(rho_s) ; - coordinates_centre->g_33 /= SQ(rho_s) ; - coordinates_centre->g_12 *= B_0 ; - coordinates_centre->g_13 *= B_0 ; - coordinates_centre->g_23 /= SQ(rho_s) ; - coordinates_centre->g11 /= SQ(rho_s*B_0) ; - coordinates_centre->g22 *= SQ(rho_s) ; - coordinates_centre->g33 *= SQ(rho_s) ; - coordinates_centre->g12 /= B_0 ; - coordinates_centre->g13 /= B_0 ; - coordinates_centre->g23 *= SQ(rho_s) ; - coordinates_centre->dx /= SQ(rho_s)*B_0; - coordinates_centre->J *= B_0/rho_s; - coordinates_centre->geometry(); if (realistic_geometry != "none") { // Load curvature operator bxcv.x /= B_0; @@ -320,14 +324,26 @@ int STORM::init(bool restarting) { G3 = coordinates_centre->G3; } } + coordinates_stag->outputVars(dump); + // Save additional coordinate objects + // Note: Rxy, Zxy and psixy are *not* normalised + if (mesh->get(Rxy, "Rxy") == 0) { + mesh->communicate(Rxy); + SAVE_ONCE(Rxy); + } + if (mesh->get(Zxy, "Zxy") == 0) { + mesh->communicate(Zxy); + SAVE_ONCE(Zxy); + } + if (mesh->get(psixy, "psixy") == 0) { + mesh->communicate(psixy); + SAVE_ONCE(psixy); + } if (realistic_geometry != "none") { - // Save additional coordinate objects dump.addOnce(bxcv.x, "bxcvx"); dump.addOnce(bxcv.y, "bxcvy"); dump.addOnce(bxcv.z, "bxcvz"); - coordinates_stag = mesh->getCoordinates(CELL_YLOW); - coordinates_stag->outputVars(dump); } int bracket; @@ -380,10 +396,10 @@ int STORM::init(bool restarting) { // ******** Initialise constant fields ******** if (!sources_realistic_geometry) { - S = FieldFactory::get()->create3D("S:function", Options::getRoot(), mesh, CELL_CENTRE, 0); - S_stag = FieldFactory::get()->create3D("S:function", Options::getRoot(), mesh, CELL_YLOW, 0); + S = FieldFactory::get()->create3D("S:function", &globalOptions, mesh, CELL_CENTRE, 0); + S_stag = FieldFactory::get()->create3D("S:function", &globalOptions, mesh, CELL_YLOW, 0); if (!isothermal) { - S_E = FieldFactory::get()->create3D("S_E:function", Options::getRoot(), mesh, CELL_CENTRE, 0); + S_E = FieldFactory::get()->create3D("S_E:function", &globalOptions, mesh, CELL_CENTRE, 0); } if (normalise_sources) { @@ -424,8 +440,10 @@ int STORM::init(bool restarting) { logT_stag = 0.; T_stag = exp(logT_stag, "RGN_NOY"); sqrtT = 1.; - T_sheath_upper = exp(extrap_sheath_upper(logT_aligned)); - if (!symmetry_plane) { + if (mesh->hasBndryUpperY()) { + T_sheath_upper = exp(extrap_sheath_upper(logT_aligned)); + } + if (!symmetry_plane and mesh->hasBndryLowerY()) { T_sheath_lower = exp(extrap_sheath_lower(logT_aligned)); } SAVE_ONCE(T); @@ -456,13 +474,22 @@ int STORM::init(bool restarting) { } } - Options* phiOptions; - if (globalOptions->isSection("phiSolver")) { - phiOptions = globalOptions->getSection("phiSolver"); + Options& phiOptions = globalOptions.isSection("phiSolver") + ? *globalOptions.getSection("phiSolver") + : *globalOptions.getSection("laplace"); + if (split_n0) { + // Set DC part to zero - this part will be solved by phiSolverxy + phiOptions["global_flags"] = phiOptions["global_flags"].withDefault(1); + // keep boundary flags set to 0 (zero-value Dirichlet), which is the BOUT++ default } else { - phiOptions = globalOptions->getSection("laplace"); - } - phiSolver = Laplacian::create(phiOptions); + // Set the radial boundary conditions to non-zero Dirichlet, with the value passed + // in the boundary cells of the 'initial guess' + phiOptions["inner_boundary_flags"] + = phiOptions["inner_boundary_flags"].withDefault(16); + phiOptions["outer_boundary_flags"] + = phiOptions["outer_boundary_flags"].withDefault(16); + } + phiSolver = Laplacian::create(&phiOptions); phi = phi_wall ; // Sets phi field to 0 for inital 'guess' of Laplace solver phi.setBoundary("phi") ; phi_aligned.setBoundary("phi_aligned") ; @@ -488,11 +515,12 @@ int STORM::init(bool restarting) { } phi2D = 0.; restart.addOnce(phi2D,"phi2D"); + SAVE_REPEAT(phi2D); } if (electromagnetic) { - auto psi_options = globalOptions->getSection("psiSolver"); - psiSolver = Laplacian::create(psi_options, CELL_YLOW); + Options& psi_options = *globalOptions.getSection("psiSolver"); + psiSolver = Laplacian::create(&psi_options, CELL_YLOW); // Check that the Laplacian solver uses 3D coefficients ASSERT0(psiSolver->uses3DCoefs()); @@ -503,7 +531,7 @@ int STORM::init(bool restarting) { psi_aligned.setBoundary("psi_aligned"); psi_centre.setBoundary("psi_centre"); - use_psi_boundary_solver = (*psi_options)["use_psi_boundary_solver"].withDefault(true); + use_psi_boundary_solver = psi_options["use_psi_boundary_solver"].withDefault(true); if (use_psi_boundary_solver) { // Check there isn't a limiter-type boundary where we can't use a solver on the // boundary point @@ -515,20 +543,20 @@ int STORM::init(bool restarting) { or (ixseps1 <= 0 and ixseps2 <= 0) or (ixseps1 >= mesh->GlobalNx and ixseps2 >= mesh->GlobalNx)); - std::string psi_solver_type = (*globalOptions)["psiSolver"]["type"]; + std::string psi_solver_type = globalOptions["psiSolver"]["type"]; if (psi_solver_type == "naulin") { // Actually only need the FFT sub-solver as at the boundary we solve Delp2(psi) = // n*(U-V) which has constant coefficients - psiBoundarySolver = Laplacian::create(psi_options->getSection("delp2solver"), CELL_YLOW); + psiBoundarySolver = Laplacian::create(psi_options.getSection("delp2solver"), CELL_YLOW); // Set boundary flags from psiSolver options, like LaplaceNaulin does in its // constructor psiBoundarySolver->setInnerBoundaryFlags( - (*psi_options)["inner_boundary_flags"].as()); + psi_options["inner_boundary_flags"].as()); psiBoundarySolver->setOuterBoundaryFlags( - (*psi_options)["outer_boundary_flags"].as()); + psi_options["outer_boundary_flags"].as()); } else { - psiBoundarySolver = Laplacian::create(psi_options, CELL_YLOW); + psiBoundarySolver = Laplacian::create(&psi_options, CELL_YLOW); } if (mesh->hasBndryLowerY() and not symmetry_plane) { @@ -660,7 +688,7 @@ int STORM::init(bool restarting) { //************** Output of Summary ********************** output.write("\n*******************************************************************") ; - output.write("\nGit Version of STORM: %s", STORM_VERSION) ; + output.write("\nGit Version of STORM: %s", storm_git_hash) ; output.write("\n*******************************************************************") ; output.write("\nCalculated Parameters") ; output.write("\n*******************************************************************") ; @@ -721,12 +749,12 @@ int STORM::init(bool restarting) { int i = 0; BoutReal xpos, ypos, zpos; int ix, iy, iz; - Options* fast_output_options = globalOptions->getSection("fast_output"); + Options& fast_output_options = *globalOptions.getSection("fast_output"); while (true) { // Add more points if explicitly set in input file - fast_output_options->get("xpos"+std::to_string(i), xpos, -1.); - fast_output_options->get("ypos"+std::to_string(i), ypos, -1.); - fast_output_options->get("zpos"+std::to_string(i), zpos, -1.); + xpos = fast_output_options["xpos"+std::to_string(i)].withDefault(-1.); + ypos = fast_output_options["ypos"+std::to_string(i)].withDefault(-1.); + zpos = fast_output_options["zpos"+std::to_string(i)].withDefault(-1.); if (xpos<0. || ypos<0. || zpos<0.) { output.write("\tAdded %i fast_output points\n", i); break; @@ -758,7 +786,7 @@ int STORM::init(bool restarting) { // Neutral models TRACE("Initialising neutral models"); - neutrals = NeutralModel::create(solver, Options::root()["neutral"]); + neutrals = NeutralModel::create(solver, globalOptions["neutral"], dump); // Set normalisations if (neutrals != NULL) { @@ -779,6 +807,8 @@ int STORM::rhs(BoutReal time) { if(verbose){ output<<"\r\t\t\t\t\t\t\t t = "<communicate(comms) ; if (average_radial_boundaries_core_SOL) { average_Z_bndry(logn); @@ -855,13 +885,17 @@ int STORM::rhs(BoutReal time) { //********** Extrapolate BC for density and temperature to calculate BC on V, U and q_par ********** // Extrapolate in field-aligned coordinates since these are used for parallel boundary // conditions - n_sheath_upper = exp(extrap_sheath_upper(logn_aligned)); - if(!symmetry_plane){ + if (mesh->hasBndryUpperY()) { + n_sheath_upper = exp(extrap_sheath_upper(logn_aligned)); + } + if(!symmetry_plane and mesh->hasBndryLowerY()){ n_sheath_lower = exp(extrap_sheath_lower(logn_aligned)); } if (!isothermal) { - T_sheath_upper = exp(extrap_sheath_upper(logT_aligned)); - if(!symmetry_plane){ + if (mesh->hasBndryUpperY()) { + T_sheath_upper = exp(extrap_sheath_upper(logT_aligned)); + } + if(!symmetry_plane and mesh->hasBndryLowerY()){ T_sheath_lower = exp(extrap_sheath_lower(logT_aligned)); } } @@ -962,8 +996,12 @@ int STORM::rhs(BoutReal time) { } // Extrapolate phi at boundary for V boundary condition - phi_sheath_lower = extrap_sheath_lower(phi_aligned) ; - phi_sheath_upper = extrap_sheath_upper(phi_aligned) ; + if (not symmetry_plane and mesh->hasBndryLowerY()) { + phi_sheath_lower = extrap_sheath_lower(phi_aligned) ; + } + if (mesh->hasBndryUpperY()) { + phi_sheath_upper = extrap_sheath_upper(phi_aligned) ; + } if (electromagnetic) { //********* Laplace inversion: calculation of psi, U and V ********* @@ -1025,12 +1063,14 @@ int STORM::rhs(BoutReal time) { V_aligned.applyBoundary(time); //********** Sheath BC on the ion and electron velocity, U and V ******** - if (!symmetry_plane) { + if (!symmetry_plane and mesh->hasBndryLowerY()) { Vsheath_ydown_staggered(V_aligned, phi_sheath_lower, phi_wall, T_sheath_lower, Vsheath_BC_prefactor) ; Usheath_ydown_staggered(U_aligned, sqrt(T_sheath_lower), Usheath_BC_prefactor); } - Vsheath_yup_staggered(V_aligned, phi_sheath_upper, phi_wall, T_sheath_upper, Vsheath_BC_prefactor); - Usheath_yup_staggered(U_aligned, sqrt(T_sheath_upper), Usheath_BC_prefactor); + if (mesh->hasBndryUpperY()) { + Vsheath_yup_staggered(V_aligned, phi_sheath_upper, phi_wall, T_sheath_upper, Vsheath_BC_prefactor); + Usheath_yup_staggered(U_aligned, sqrt(T_sheath_upper), Usheath_BC_prefactor); + } if (mesh->yend - mesh->ystart + 1 < 3) { // Need to communicate mesh->ystart point that has just been set by boundary conditions mesh->communicate(U_aligned, V_aligned); @@ -1107,70 +1147,79 @@ int STORM::rhs(BoutReal time) { // ***Vorticity Equation*** if (hydrodynamic) { - ddt(vort) = 0.; + vorticity_equation["all"] = 0.; } else { - ddt(vort) = - bracket(phi, vort, bm, CELL_CENTRE) - - Vpar_Grad_par_EM(U_aligned, vort_aligned, U_centre, vort, - psi_centre, CELL_CENTRE); + vorticity_equation["ExB_advection"] = - bracket(phi, vort, bm, CELL_CENTRE); + vorticity_equation["parallel_advection"] = + - Vpar_Grad_par_EM(U_aligned, vort_aligned, U_centre, vort, psi_centre, + CELL_CENTRE); if (boussinesq == 1 || boussinesq == 3 || boussinesq == 4) { // Grad(n*Grad_perp(phi)) = n*Delp2(phi) - ddt(vort) += Div_par_EM(UmV_aligned, UmV_centre, psi_centre, CELL_CENTRE) - + UmV_centre*Grad_par_EM(logn_aligned, logn, psi_centre, CELL_CENTRE); + vorticity_equation["parallel_current"] = + Div_par_EM(UmV_aligned, UmV_centre, psi_centre, CELL_CENTRE) + + UmV_centre*Grad_par_EM(logn_aligned, logn, psi_centre, CELL_CENTRE); } else if (boussinesq == 0 || boussinesq == 2) { // Grad(n*Grad_perp(phi)) = Delp2(phi) or non-boussinesq - ddt(vort) += n*( Div_par_EM(UmV_aligned, UmV_centre, psi_centre, CELL_CENTRE) - + UmV_centre*Grad_par_EM(logn_aligned, logn, psi_centre, CELL_CENTRE) ); + vorticity_equation["parallel_current"] = + n*( Div_par_EM(UmV_aligned, UmV_centre, psi_centre, CELL_CENTRE) + + UmV_centre*Grad_par_EM(logn_aligned, logn, psi_centre, CELL_CENTRE) ); } else { throw BoutException("Boussinesq option unrecognized."); } if (boussinesq == 0) { - ddt(vort) -= bracket(uE2/2.,n,bm); //check normalisation of last term + vorticity_equation["bracket_uE2"] = -bracket(uE2/2.,n,bm); //check normalisation of last term } if (boussinesq == 3) { - ddt(vort) -= S*vort + coordinates_centre->g11*DDX(S)*DDX(phi)/B2; // assumes S axisymmetric + // assumes S axisymmetric + vorticity_equation["density_source"] = - S*vort + - coordinates_centre->g11*DDX(S)*DDX(phi)/B2; } if(uniform_diss_paras){ - ddt(vort) += mu_vort0*Delp2(vort) ; + vorticity_equation["perpendicular_diffusion"] = mu_vort0*Delp2(vort); } else{ - ddt(vort) += mu_vort*Delp2(vort) + Grad_perp_dot_Grad_perp(mu_vort, vort); + vorticity_equation["perpendicular_diffusion"] = + mu_vort*Delp2(vort) + Grad_perp_dot_Grad_perp(mu_vort, vort); } // Curvature terms for vorticity if (boussinesq == 1 || boussinesq == 3 || boussinesq == 4) { - ddt(vort) += Curv_p/n; + vorticity_equation["curvature"] = Curv_p/n; } else { - ddt(vort) += Curv_p; + vorticity_equation["curvature"] = Curv_p; } } Field3D div_par_V = Div_par_EM(V_aligned, V_centre, psi_centre, CELL_CENTRE); // ***Density Equation*** - ddt(logn) = - bracket(phi, logn, bm, CELL_CENTRE) - - div_par_V - - Vpar_Grad_par_EM(V_aligned, logn_aligned, V_centre, logn, psi_centre, CELL_CENTRE) - + S/n; + density_equation["ExB_advection"] = - bracket(phi, logn, bm, CELL_CENTRE); + density_equation["parallel_advection"] = + - div_par_V + - Vpar_Grad_par_EM(V_aligned, logn_aligned, V_centre, logn, psi_centre, CELL_CENTRE); + density_equation["density_source"] = S/n; if (uniform_diss_paras) { - ddt(logn) += mu_n0*Delp2(n)/n; + density_equation["perpendicular_diffusion"] = mu_n0*Delp2(n)/n; } else { // Include y derivatives?? - ddt(logn) += mu_n*Delp2(n)/n + Grad_perp_dot_Grad_perp(mu_n, logn) ; + density_equation["perpendicular_diffusion"] = mu_n*Delp2(n)/n + + Grad_perp_dot_Grad_perp(mu_n, logn); } // Curvature terms for density - ddt(logn) += -Curv_phi + Curv_p/n; + density_equation["ExB_curvature"] = -Curv_phi; + density_equation["diamagnetic_curvature"] = Curv_p/n; if (boussinesq == 4) { - ddt(vort) -= vort*ddt(logn); + vorticity_equation["vort*dndt"] = -vort*ddt(logn); } // ***Ion parallel velocity Equation*** @@ -1185,34 +1234,37 @@ int STORM::rhs(BoutReal time) { } else { Grad_par_phi_stag = Grad_par_EM(phi_aligned, phi_stag, psi, CELL_YLOW); } - ddt(chiU) = - bracket(phi_stag, U, bm, CELL_YLOW) - - Vpar_Grad_par_EM(U_aligned, U_aligned, U, U, psi, CELL_YLOW) - - Grad_par_phi_stag - - (nu_parallel/mu)*UmV - - (U*S_stag)/n_stag - + 0.71*Grad_par_T_stag; + ion_momentum_equation["ExB_advection"] = - bracket(phi_stag, U, bm, CELL_YLOW); + ion_momentum_equation["parallel_advection"] = + - Vpar_Grad_par_EM(U_aligned, U_aligned, U, U, psi, CELL_YLOW); + ion_momentum_equation["grad_phi"] = - Grad_par_phi_stag; + ion_momentum_equation["resistivity"] = - (nu_parallel/mu)*UmV; + ion_momentum_equation["thermal_force"] = 0.71*Grad_par_T_stag; + ion_momentum_equation["density_source"] = - (U*S_stag)/n_stag; if (diff_perp_U > 0.) { - ddt(chiU) += diff_perp_U*Delp2(U); + ion_momentum_equation["perpendicular_diffusion"] = diff_perp_U*Delp2(U); } set_lower_ddt_zero(chiU) ; // ***Electron parallel velocity Equation*** if (hydrodynamic) { - ddt(chiV) = ddt(chiU); + electron_momentum_equation["all"] = ddt(chiU); } else { - ddt(chiV) = - bracket(phi_stag, V, bm, CELL_YLOW) - - Vpar_Grad_par_EM(V_aligned, V_aligned, V, V, psi, CELL_YLOW) - + mu*Grad_par_phi_stag - + nu_parallel*UmV - - mu*T_stag*Grad_par_EM(logn_aligned, logn_stag, psi, CELL_YLOW) - - (1.71*mu)*Grad_par_T_stag - - (V*S_stag)/n_stag ; + electron_momentum_equation["ExB_advection"] = - bracket(phi_stag, V, bm, CELL_YLOW); + electron_momentum_equation["parallel_advection"] = + - Vpar_Grad_par_EM(V_aligned, V_aligned, V, V, psi, CELL_YLOW); + electron_momentum_equation["grad_phi"] = mu*Grad_par_phi_stag; + electron_momentum_equation["resistivity"] = nu_parallel*UmV; + electron_momentum_equation["density_gradient"] = + - mu*T_stag*Grad_par_EM(logn_aligned, logn_stag, psi, CELL_YLOW); + electron_momentum_equation["temperature_gradient"] = - (1.71*mu)*Grad_par_T_stag; + electron_momentum_equation["density_source"] = - (V*S_stag)/n_stag; if (diff_perp_V > 0.) { - ddt(chiV) += diff_perp_V*Delp2(V); + electron_momentum_equation["perpendicular_diffusion"] = diff_perp_V*Delp2(V); } set_lower_ddt_zero(chiV) ; @@ -1245,36 +1297,44 @@ int STORM::rhs(BoutReal time) { //***Electron temperature Equation*** Tcoef = (2.0/3.0)/p; - ddt(logp) = - bracket(phi,logp,bm,CELL_CENTRE) - - Vpar_Grad_par_EM(V_aligned, logn_aligned + logT_aligned, - V_centre, logp, psi_centre, CELL_CENTRE) - - Tcoef*Div_par_EM(qpar_aligned, qpar_centre, psi_centre,CELL_CENTRE) - - 2./3.*0.71*UmV_centre*Grad_par_EM(logT_aligned, logT, psi_centre, CELL_CENTRE) - - (5./3.)*div_par_V - + 2.0/(3.0*mu)*nu_parallel0*(pow(T, -2.5, "RGN_NOBNDRY"))*n*SQ(UmV_centre) - + Tcoef*S_E; + electron_pressure_equation["ExB_advection"] = - bracket(phi,logp,bm,CELL_CENTRE); + electron_pressure_equation["parallel_advection"] = + - Vpar_Grad_par_EM(V_aligned, logn_aligned + logT_aligned, V_centre, logp, + psi_centre, CELL_CENTRE); + electron_pressure_equation["parallel_conduction"] = + - Tcoef*Div_par_EM(qpar_aligned, qpar_centre, psi_centre,CELL_CENTRE); + electron_pressure_equation["thermal_force"] = + - 2./3.*0.71*UmV_centre*Grad_par_EM(logT_aligned, logT, psi_centre, CELL_CENTRE); + electron_pressure_equation["compression"] = - (5./3.)*div_par_V; + electron_pressure_equation["resistive_heating"] = + 2.0/(3.0*mu)*nu_parallel0*(pow(T, -2.5, "RGN_NOBNDRY"))*n*SQ(UmV_centre); + electron_pressure_equation["energy_source"] = Tcoef*S_E; if(uniform_diss_paras){ - ddt(logp) += Tcoef*kappa0_perp*Delp2(T); + electron_pressure_equation["perpendicular_diffusion"] = Tcoef*kappa0_perp*Delp2(T); } else{ - ddt(logp) += Tcoef*(kappa_perp*Delp2(T) + Grad_perp_dot_Grad_perp(kappa_perp, T)); + electron_pressure_equation["perpendicular_diffusion"] = + Tcoef*(kappa_perp*Delp2(T) + Grad_perp_dot_Grad_perp(kappa_perp, T)); } if (S_in_peq) { - ddt(logp) += S*SQ(V_centre)/(3.*mu*p); + electron_pressure_equation["S*V^2/3/mu/p"] = S*SQ(V_centre)/(3.*mu*p); } //Curvature terms for electron temperature - ddt(logp) += -5./3.*Curv_phi //Divergence of ExB - + 5./3.*(Curv_p/n + Curv_T) //Divergence of diamagnetic flow and diamagnetic heat flux - - SQ(V_centre)*Curv_p/(3.*mu*p); //gyro-viscous energy transfer term + // Divergence of ExB + electron_pressure_equation["ExB_curvature"] = -5./3.*Curv_phi; + // Divergence of diamagnetic flow and diamagnetic heat flux + electron_pressure_equation["diamagnetic_curvature"] = 5./3.*(Curv_p/n + Curv_T); + // gyro-viscous energy transfer term + electron_pressure_equation["gyroviscous_curvature"] = - SQ(V_centre)*Curv_p/(3.*mu*p); if (run_1d) { // assume we are only setting up equilibrium: // slow down the T evolution so we can take longer timesteps // assume that we solve for ddt(logp) = 0 - ddt(logp) = ddt(logp)/run_1d_T_slowdown + ddt(logn); + ddt(logp) = (ddt(logp) - ddt(logn))/run_1d_T_slowdown + ddt(logn); } } @@ -1291,25 +1351,30 @@ int STORM::rhs(BoutReal time) { neutrals->update(n, n_stag, U, V, T, T_stag, time); // Add sources/sinks to plasma equations - ddt(logn) -= neutrals->S/n; + density_equation["ionisation_recombination"] = -neutrals->S/n; Field3D Son_stag = neutrals->S_stag/n_stag; - ddt(chiU) += -neutrals->Fi/n_stag + Son_stag*U; + ion_momentum_equation["neutral_friction"] = -neutrals->Fi/n_stag; + ion_momentum_equation["neutral_source"] = Son_stag*U; - ddt(chiV) += -neutrals->Fe/n_stag + Son_stag*V; + electron_momentum_equation["neutral_friction"] = -neutrals->Fe/n_stag; + electron_momentum_equation["neutral_source"] = Son_stag*V; if (boussinesq == 2) { - ddt(vort) -= neutrals->Fperp*vort + Grad_perp_dot_Grad_perp(neutrals->Fperp, phi)/B2; + vorticity_equation["neutral_friction"] = + - neutrals->Fperp*vort - Grad_perp_dot_Grad_perp(neutrals->Fperp, phi)/B2; } else if (boussinesq == 1 || boussinesq == 3 || boussinesq == 4) { - ddt(vort) -= (neutrals->Fperp*vort + Grad_perp_dot_Grad_perp(neutrals->Fperp, phi)/B2)/n; + vorticity_equation["neutral_friction"] = + - (neutrals->Fperp*vort + Grad_perp_dot_Grad_perp(neutrals->Fperp, phi)/B2)/n; } else { - ddt(vort) -= neutrals->Fperp*vort/n + Grad_perp_dot_Grad_perp(neutrals->Fperp/n, phi)*n/B2; + vorticity_equation["neutral_friction"] = + - neutrals->Fperp*vort/n - Grad_perp_dot_Grad_perp(neutrals->Fperp/n, phi)*n/B2; } if (!isothermal) { - ddt(logp) -= Tcoef*neutrals->Rp; + electron_pressure_equation["neutral_radiation"] = - Tcoef*neutrals->Rp; } } diff --git a/storm3d/storm.hxx b/storm3d/storm.hxx index ed22650..eb6ffd7 100644 --- a/storm3d/storm.hxx +++ b/storm3d/storm.hxx @@ -29,7 +29,9 @@ #include #include #include "../shared/fast_output.hxx" +#include "../shared/BoutEquation/equation.hxx" #include "neutral-model.hxx" +#include "storm_version.hxx" class STORM : public PhysicsModel{ public: @@ -71,8 +73,8 @@ protected: int precon(BoutReal t, BoutReal gamma, BoutReal delta); private: - Options* globalOptions; - Options* options; + Options& globalOptions{Options::root()}; + Options& options{*globalOptions.getSection("storm")}; Coordinates* coordinates_centre; Coordinates* coordinates_stag; @@ -100,6 +102,19 @@ private: void read_equilibrium_file(BoutReal * data, const char * fname); + // rhs_counter is used by Equation objects to know if they have to reset for a new rhs evaluation + int rhs_counter = 0; + Equation vorticity_equation{vort, "vort", globalOptions["save_equations"], dump, + rhs_counter}; + Equation density_equation{logn, "logn", globalOptions["save_equations"], dump, + rhs_counter}; + Equation ion_momentum_equation{chiU, "chiU", globalOptions["save_equations"], dump, + rhs_counter}; + Equation electron_momentum_equation{chiV, "chiV", globalOptions["save_equations"], dump, + rhs_counter}; + Equation electron_pressure_equation{logp, "logp", globalOptions["save_equations"], dump, + rhs_counter}; + Field3D logn; Field3D n; // density Field3D U; // ion parallel velocity @@ -205,6 +220,7 @@ private: bool initial_perturbation ; // Switch to add perturbations specified in input file onto the equilibrium bool verbose ; // Add extra output bool evolving_bcs; // Switch to activate averaged Neumann boundary conditions for phi + bool reset_evolving_bcs; // Re-initialise, ignoring history of averaged boundary condition for phi when restarting bool use_psi_boundary_solver;// Use a Laplace solver to invert for psi on the boundary points, rather than extrapolating bool save_aligned_fields; // Switch to save field-aligned versions of fields to dump files bool average_radial_boundaries_core_SOL; // Ad hoc boundary conditions in radial direction @@ -221,6 +237,7 @@ private: Vector2D bxcv, Curlb_B; // Vectors for realistic curvature operator Field2D B2; // Bxy*Bxy Field2D G3; // Useful to temporary overwrite G3 + Field2D Rxy, Zxy, psixy; // Used to save extra grid info to dump files // BoutReal *phi_x_boundary ; FILE *file ; diff --git a/storm3d/utilities.cxx b/storm3d/utilities.cxx index 36d27ef..1776c82 100644 --- a/storm3d/utilities.cxx +++ b/storm3d/utilities.cxx @@ -261,14 +261,14 @@ Field3D STORM::this_Grad_perp2(const Field3D &f) { // If the x and z derivative schemes are not C2, then this operator is not // the inverse of the multigrid Laplacian solver std::string first; - OPTION(globalOptions->getSection("mesh:ddx"), first, "C2"); + OPTION(globalOptions.getSection("mesh:ddx"), first, "C2"); ASSERT1(first == "C2"); std::string second; - OPTION(globalOptions->getSection("mesh:ddx"), second, "C2"); + OPTION(globalOptions.getSection("mesh:ddx"), second, "C2"); ASSERT1(second == "C2"); - OPTION(globalOptions->getSection("mesh:ddz"), first, "FFT"); + OPTION(globalOptions.getSection("mesh:ddz"), first, "FFT"); ASSERT1(first == "C2"); - OPTION(globalOptions->getSection("mesh:ddz"), second, "FFT"); + OPTION(globalOptions.getSection("mesh:ddz"), second, "FFT"); ASSERT1(second == "C2"); #endif Coordinates* coords = f.getCoordinates(); @@ -284,14 +284,14 @@ Field3D STORM::this_Grad_perp_dot_Grad_perp(const Field3D &f, const Field3D &g) // If the x and z derivative schemes are not C2, then this operator is not // the inverse of the multigrid Laplacian solver std::string first; - OPTION(globalOptions->getSection("mesh:ddx"), first, "C2"); + OPTION(globalOptions.getSection("mesh:ddx"), first, "C2"); ASSERT1(first == "C2"); std::string second; - OPTION(globalOptions->getSection("mesh:ddx"), second, "C2"); + OPTION(globalOptions.getSection("mesh:ddx"), second, "C2"); ASSERT1(second == "C2"); - OPTION(globalOptions->getSection("mesh:ddz"), first, "FFT"); + OPTION(globalOptions.getSection("mesh:ddz"), first, "FFT"); ASSERT1(first == "C2"); - OPTION(globalOptions->getSection("mesh:ddz"), second, "FFT"); + OPTION(globalOptions.getSection("mesh:ddz"), second, "FFT"); ASSERT1(second == "C2"); #endif Coordinates* coords = f.getCoordinates(); @@ -387,7 +387,7 @@ void STORM::set_sources_realistic_geometry() { //PF region between jyseps2_2 and NyGlobal for(int iy = mesh->ystart; iy <= mesh->yend; ++iy){ for(int iz = 0; iz < mesh->LocalNz; ++iz){ - dy = 0.5/((double)(mesh->GlobalNy - jyseps2_2 - 1)); + dy = 0.5/((double)(mesh->GlobalNy -2*mesh->ystart - jyseps2_2 - 1)); y = ((double)(mesh->getGlobalYIndexNoBoundaries(iy) - jyseps2_2) - 0.5)*dy + 0.5; S(ix,iy,iz) += exp(10.*std::abs(y-0.5))/(exp(5.0)-1.0); } @@ -425,7 +425,7 @@ void STORM::set_sources_realistic_geometry() { else if (mesh->getGlobalYIndexNoBoundaries(mesh->ystart) > jyseps1_2-1) { for(int iy = mesh->ystart; iy <= mesh->yend; ++iy){ for(int iz = 0; iz < mesh->LocalNz; ++iz){ - dy = 0.5/((double)(mesh->GlobalNy - jyseps1_2 - 1)); + dy = 0.5/((double)(mesh->GlobalNy -2*mesh->ystart - jyseps1_2 - 1)); y = ((double)(mesh->getGlobalYIndexNoBoundaries(iy) - jyseps1_2) - 0.5)*dy + 0.5; S(ix,iy,iz) += exp(10.*std::abs(y-0.5))/(exp(5.0)-1.0); } @@ -434,7 +434,7 @@ void STORM::set_sources_realistic_geometry() { } else if (mesh->getGlobalXIndex(ix) >= ixseps2) { if (mesh->getGlobalYIndexNoBoundaries(mesh->yend) < ny_inner) { - for(int iy = mesh->xstart; iy <= mesh->yend; ++iy){ + for(int iy = mesh->ystart; iy <= mesh->yend; ++iy){ for(int iz = 0; iz < mesh->LocalNz; ++iz){ dy = 1./((double)(ny_inner)); y = ((double)(mesh->getGlobalYIndexNoBoundaries(iy)) + 0.5)*dy; @@ -445,7 +445,7 @@ void STORM::set_sources_realistic_geometry() { else { for(int iy = mesh->ystart; iy <= mesh->yend; ++iy){ for(int iz = 0; iz < mesh->LocalNz; ++iz){ - dy = 1./((double)(mesh->GlobalNy - ny_inner)); + dy = 1./((double)(mesh->GlobalNy -2*mesh->ystart - ny_inner)); y = ((double)(mesh->getGlobalYIndexNoBoundaries(iy)) - ny_inner + 0.5)*dy; S(ix,iy,iz) += exp(10.*std::abs(y-0.5))/(exp(5.0)-1.0); } diff --git a/tests3d/run_3d_filament_test.py b/tests3d/run_3d_filament_test.py index fab57a3..bce2bc4 100755 --- a/tests3d/run_3d_filament_test.py +++ b/tests3d/run_3d_filament_test.py @@ -40,7 +40,7 @@ # x-boundary cells are not set or used for these variables, so don't check them noXBoundaryVariables = ['qpar', 'chiU', 'chiV'] # y-boundary cells are not set or used for these variables, so don't check them -noYBoundaryVariables = ['uE2', 'chiU', 'chiV'] +noYBoundaryVariables = ['n', 'T', 'uE2', 'chiU', 'chiV'] def test(numProcs,retestOutput=False): diff --git a/tests3d/test-3d/data/BOUT.inp b/tests3d/test-3d/data/BOUT.inp index 8fec2f2..cd09c8a 100644 --- a/tests3d/test-3d/data/BOUT.inp +++ b/tests3d/test-3d/data/BOUT.inp @@ -3,38 +3,19 @@ nout = 1 # Number of outputted timesteps restart = true dump_on_restart = false - MZ = 128 # Number of Z points zmin = 0 zmax = 15.9155 # 1.5915*2pi = 10, z is fracs of 2pi #NXPE = 6 # Sets the number of processors in the x direction MXG = 2 # Number of X guard cells -MYG = 2 # Number of Y guard cells - -periodicX = 0 # Sets periodicity in the x direction 1=on, 0 =off -#ballooning = true # Not sure what this does - -[mesh:ddx] -first = C2 -second = C2 -upwind = C2 -flux = SPLIT [mesh:ddy] -first = C2 -second = C2 upwind = U2 -flux = SPLIT [mesh:ddz] first = C2 second = C2 -upwind = C2 -flux = SPLIT - -[output] -floats = false # Sets outputs to floats/doubles [mesh] StaggerGrids = true # Enables staggered grids @@ -47,23 +28,8 @@ dy = Ly/ny ixseps1 = -1 # Set x location of separatrix 1 ixseps2 = 0 # Set x location of separatrix 2 -#ny_inner = -#jyseps1_1 = -1 -#jyseps1_2 = -1 -#jyseps2_1 = -1 -#jyseps2_2 = 63 - -symmetricGlobalY = true - -[phiSolver] -type = cyclic -global_flags = 0 -inner_boundary_flags = 16 -outer_boundary_flags = 16 -include_yguards = false [solver] -mms = false type = rk4 timestep = 0.01 @@ -77,7 +43,6 @@ R_c = 1.5 # m n_0 = 0.8e19 # m^-3 Z = 1 # Dimensionless loglambda = -1 # Dimensionless -phi_wall = 0 # If these parameters are specified, they will be used instead of the values calculated from the primary parameters above. # mu_n0 = 0.01 @@ -86,7 +51,6 @@ phi_wall = 0 # nu_parallel0 = 1000 # g0 = 0.0025 -old_phi_wall_value = false isothermal = false boussinesq = 1 uniform_diss_paras = false @@ -104,7 +68,6 @@ curv_T_gyro = true init_bg = true add_blob = true symmetry_plane = true -run_1d = false [blob] delta_perp = 10.0 @@ -127,15 +90,7 @@ zoffset_T = 0.5 angle_blob_T = 0 delta_front_T = 0.1 -[turb] -max_xmode = 25 -max_ymode = 25 -max_zmode = 25 -amplitude_n = 0.001 -amplitude_T = 0.001 - [S] -scale = 1.0 # function = 10*exp(10.0*y/(2*pi))/(exp(10.0)-1.0) # Sheath loc # function = 10*exp(10.0*(y-pi)/(pi))/(0.5*(exp(10.0)-1.0))+10*exp(10.0*(pi-y)/(pi))/(0.5*(exp(10.0)-1.0)) # Sheath loc, no sym plane # function = 2*sqrt(100/pi)*exp(-100*(0.5*y/pi)^2) # Midplane loc @@ -145,32 +100,25 @@ function = 0.7*10*exp(10.0*y/(2*pi))/(exp(10.0)-1.0) # function = 0.7*10*(exp(10.0*(y-pi)/pi)+exp(-10.0*(y-pi)/pi))/(exp(10.0)-1.0)/0.5 # Sheath loc, no sym plane, norm-ed so n = 1 at midplane [S_E] -scale = 1.0; function = 16.0*exp(-5.0*y/(2*pi)) # function = 16.0*(exp(-5.0*(y-pi)/pi)*H(y-pi)+exp(5.0*(y-pi)/pi)*H(-y+pi))/0.5 # no sym plane [n] -scale = 1.0 function = 1.0 [T] -scale = 1.0 function = 1.0 [vort] -scale = 1.0 function = 0.0 [U] -scale = 1.0 function = exp(-10)*exp( 10*(y/(2.0*pi))) #if symmetry_plane=true # function = exp(-10)*(exp( 10*((y-pi)/pi)) - exp( 10*((pi-y)/pi))) #if symmetry_plane=false [V] -scale = 1.0 function = exp(-10)*exp( 10*(y/(2.0*pi))) #if symmetry_plane=true # function = exp(-10)*(exp( 10*((y-pi)/pi)) - exp( 10*((pi-y)/pi))) #if symmetry_plane=false [phi] -scale = 1.0 -function = (storm:phi_wall) + 3.18 +function = 3.18 diff --git a/tests3d/test-equilib/data/BOUT.inp b/tests3d/test-equilib/data/BOUT.inp index 79c950f..b19d591 100644 --- a/tests3d/test-equilib/data/BOUT.inp +++ b/tests3d/test-equilib/data/BOUT.inp @@ -1,37 +1,18 @@ timestep = 25 # Timestep length of outputted data nout = 20 # Number of outputted timesteps - MZ = 128 # Number of Z points zmin = 0 zmax = 15.9155 # 1.5915*2pi = 10, z is fracs of 2pi MXG = 2 # Number of X guard cells -MYG = 2 # Number of Y guard cells - -periodicX = 0 # Sets periodicity in the x direction 1=on, 0 =off -#ballooning = true # Not sure what this does - -[mesh:ddx] -first = C2 -second = C2 -upwind = C2 -flux = SPLIT [mesh:ddy] -first = C2 -second = C2 upwind = U2 -flux = SPLIT [mesh:ddz] first = C2 second = C2 -upwind = C2 -flux = SPLIT - -[output] -floats = false # Sets outputs to floats/doubles [mesh] StaggerGrids = true # Enables staggered grids @@ -44,23 +25,8 @@ dy = Ly/ny ixseps1 = -1 # Set x location of separatrix 1 ixseps2 = 0 # Set x location of separatrix 2 -#ny_inner = -#jyseps1_1 = -1 -#jyseps1_2 = -1 -#jyseps2_1 = -1 -#jyseps2_2 = 63 - -symmetricGlobalY = true - -[phiSolver] -type = cyclic -global_flags = 0 -inner_boundary_flags = 16 -outer_boundary_flags = 16 -include_yguards = false [solver] -mms = false type = rk4 timestep = 10. @@ -74,7 +40,6 @@ R_c = 1.5 # m n_0 = 0.8e19 # m^-3 Z = 1 # Dimensionless loglambda = -1 # Dimensionless -phi_wall = 0 # If these parameters are specified, they will be used instead of the values calculated from the primary parameters above. # mu_n0 = 0.01 @@ -83,7 +48,6 @@ phi_wall = 0 # nu_parallel0 = 1000 # g0 = 0.0025 -old_phi_wall_value = false isothermal = false boussinesq = 1 uniform_diss_paras = false @@ -101,7 +65,6 @@ curv_T_gyro = true init_bg = true add_blob = true symmetry_plane = true -run_1d = false run_1d_T_slowdown = 20. [blob] @@ -125,15 +88,7 @@ zoffset_T = 0.5 angle_blob_T = 0 delta_front_T = 0.1 -[turb] -max_xmode = 25 -max_ymode = 25 -max_zmode = 25 -amplitude_n = 0.001 -amplitude_T = 0.001 - [S] -scale = 1.0 # function = 10*exp(10.0*y/(2*pi))/(exp(10.0)-1.0) # Sheath loc # function = 10*exp(10.0*(y-pi)/(pi))/(0.5*(exp(10.0)-1.0))+10*exp(10.0*(pi-y)/(pi))/(0.5*(exp(10.0)-1.0)) # Sheath loc, no sym plane # function = 2*sqrt(100/pi)*exp(-100*(0.5*y/pi)^2) # Midplane loc @@ -143,32 +98,25 @@ function = 0.7*10*exp(10.0*y/(2*pi))/(exp(10.0)-1.0) # function = 0.7*10*(exp(10.0*(y-pi)/pi)+exp(-10.0*(y-pi)/pi))/(exp(10.0)-1.0)/0.5 # Sheath loc, no sym plane, norm-ed so n = 1 at midplane [S_E] -scale = 1.0; function = 16.0*exp(-5.0*y/(2*pi)) # function = 16.0*(exp(-5.0*(y-pi)/pi)*H(y-pi)+exp(5.0*(y-pi)/pi)*H(-y+pi))/0.5 # no sym plane [n] -scale = 1.0 function = 1.0 [T] -scale = 1.0 function = 1.0 [vort] -scale = 1.0 function = 0.0 [U] -scale = 1.0 function = exp(-10)*exp( 10*(y/(2.0*pi))) #if symmetry_plane=true # function = exp(-10)*(exp( 10*((y-pi)/pi)) - exp( 10*((pi-y)/pi))) #if symmetry_plane=false [V] -scale = 1.0 function = exp(-10)*exp( 10*(y/(2.0*pi))) #if symmetry_plane=true # function = exp(-10)*(exp( 10*((y-pi)/pi)) - exp( 10*((pi-y)/pi))) #if symmetry_plane=false [phi] -scale = 1.0 -function = (storm:phi_wall) + 3.18 +function = 3.18 diff --git a/tests3d/test-equilib/data/expected_equilib/T_eq.dat b/tests3d/test-equilib/data/expected_equilib/T_eq.dat index c3b3c857304a2e79c4bfb56f8e6ca88b96691f1e..d22860cf1f01d4b3c39a74feca4eb61991139f15 100644 GIT binary patch literal 160 zcmZQzKm}3dsvoYzzO!$R_G&v2_|ATp#+1!FY~I<2Y>fW*L*bqM5xJ|YWf*|+=ADau zZ@;y_%A2d6y7{er>f-~|%$;xT&o9`;vMA=QedM`=+ZXA*wcq_p>rVXtH}>-JnyL=x t-`K~zZQOfx-Wz-SaNb0b{5STeD_-o`?)b+3Qq#&q3Q}+EQC-8p000&?Ln#0N literal 160 zcmZQz0D}ey#n4b2*7@N|>^pmwDeKw}1irJk{rzL}4x4xOVZI0d{ZM#kKV!M|8X1On z_M47>SnPZIt$pLptLmwn-`YRfsa3<=`PTm2T2mdAr1W+y(tB&KS~KHL{Qo!h ymn>$eI-GxFpUL-W@6~y4>|J-%C5q&~v46-Gy?49g8~bM>vWFF<-q^$30iyw3RYgev diff --git a/tests3d/test-equilib/data/expected_equilib/U_eq.dat b/tests3d/test-equilib/data/expected_equilib/U_eq.dat index bfb4c00d24cbb7d8af49ad42d36e6907c55e1f1b..11e4a8cfa1a2648a7e25c2cce37a3f36e98fee4f 100644 GIT binary patch literal 160 zcmZ?Q`(Mj2i)X((nEtbm0ScgeyKjX;x)Nm~_TQgff5tys&Az?cDLU}1nSJox^YQMY zZuU}P&$Qych1mD_$#CUtOt!b!vbLQyt;pVe>*ahib|Ns?2T6^Bqo1)Yk&E9&`RH%f9(N$ Co<(y2 literal 160 zcmZQkp6|~vi)TLznEtbm0ScgeyM_f!ViIK{_OsNJukz1Uvo8={?jLy8%wCS6Ey7*Y z&3QydpX__Wo*`_ry07XD`&o*ahib|Ns?2T6^Bqo1)Yk&E9&`RH%f9(N$ Co<(y2 literal 160 zcmZQkp6|~vi)TLznEtbm0ScgeyM_f!ViIK{_OsNJukz1Uvo8={?jLy8%wCS6Ey7*Y z&3QydpX__Wo*`_ry07XD`&o3wR4R(7q*T(zbCRt+p_C$tHbvSfi7?Sdp+!?Bv{4e7kRqwk zMr9|dv{^2fB*mQ|9{OjGso+t$GOyf&vVXmdfm^QbDL~pCf=!2 zCmwvt!lxvBimiuV$cLCG^Mw=)VqVPGiumiz%-3SfK>sQSzjlUC3HRkI@B{eI!XNl& z@N-DOr#yV(#}%BMU0oeq-JBfPuW(xf55Nch8M(=()2BO3Hl1lb#}q#%w{(SvvlIRq z2m0CZ&&XJ(`Hx|S|JD5S==q;NebzkYH~-1=FV_wG=l*5=T}A+1a-*9E zKD7RV$4~d!c>b;b^Z6I?`c-si=-+;&cF@4ph2OZdG1Ue&1+#all# z#=>r@sK>T_qTR3U&tJ#%8P%B^F5XeI-RSd{(aU$w8jl3b42yUaLoLs3s(h=nINyJ? zKe+k(za7WF&tK-CiCu*pWgGKiAlReKaRCF}`Vmz?iEoxjG`Cu>^+#nG2% zyDG}Q{`!aWi}XFOl%?KBI?Dd-=ss?J{LR)cw?6(~&(FfYFSwDL0d5A`p8*lS?r0at z`IsqcYaZI0z_xzhpY`}(-A8rqX|q&Qxl&YPd-KP*l$=F1vVEvPvZuL>=@m-0*YI&) zM(_SM_WWUuuEuam_J+E?HFY0te zsbet3n|Rt=vM7;~{v|YUMxRZYMN@o!sxryTtSNfj&8-#pRr-UHwk&>7@9HJ$QRrch zqIL@A%px!L)*(A7}SP0;iOZou;)Sf%mqzzOf2O5Mga4p1%zV;*V&&I<1NX zzXW(h6nP@S!TpomCc>Yev)Fdc@;j*J$r|xd{+4VHdz<#nwv3~FY9TjhU(r_;+UI{v zkM^bX_oaO$S_5cbt}l4O_jJ56_BFNopQn9(@rsOZ^JK;s=tBDr^ccqY^pY9hxPgq% z@Au=4?GkQtKG?s1qxmy>SZpgey`t7An;C`s`6}5qYFH#zG^u0|Qo6Eb#oI*_kxsKg zKSxPhG)Y}W`LI_lvV5BSAXw}*a{RdJu!eyoF^@Ma{5m_fbv^dw`{ihWZ%l2^ffKM#{rOjK?7Mfkd?@xs@_N+MzOxG`+IL}7 z74{t(;ku0Woz$H``>uD3rhR+U?=U`P2gWyjMu+z`^fG8#c)%RBiucR(k85tL55fNJ z81u*e{|ltAbr>JFKDhNE%JWMHUmwErF*kBEz|BDWGk_X?O7y6>b{jQ>76ga|&qfX2 zL#LR(9gZ3nzdre_b|Gq@qNeE2utg23>&`do???O(6O{BOTtoaj`yVWORe|_XZ+_Dr zJRk8rQq%Poxgov@I&Gb0jrfXn@p0)}QT>bQo@>NUrY+1&)Q=Rd>dBphY)-k^64UBJ$C-Gfie6JDlT_rvz;>#yKSK{kOe3yu? z2jeU4!};j_`;VPJO6pPM!3P^XMK!i}{1I2Ms7AIAJ#5(N`u^b#iuZ(a)IA{eZG5b{ z^Kw=^#VcQu6Iz%<@rpdhI_T+Ay!;andxh#zyjyA#{rnRsUbcKoc#|!~OaJlXv^Vf` zNk?sZ?|Db@;-;o-8t|Ut9iKCD#GPDs=n{>_P*D*+-`s)-mel!w{2!FKU7d{f0P2LoE zdm|Do3K;OiJ6$B$*uCyr%UUD|*`2c1P8JEysQ4a!T#5u)U2oV0s3XChaPjIG2P8QA zHp@)<6KZVebWMNRPE;!)6}`ykEZf80W_*4JXy3THNXF+6oecYYN76oZgL1~VHk9`L zD7lmIJx7 z^QOJNS3fbMaDN+-ZPX*%Ht0t69;7)*>D};{D5NiN&3`kXA2JOravSoYJF<>`Y1IGO zdgSD5HguGPEAp87Wye6PACaG8^5SJ#K`7)w=!a=Jo6*t3Ha8OOQ&D_yMXG6M=tCEG z&sQJ380F4N__ffp8kK76Wvf~*V)r51C)46a`(XUe2ivqS``S9%w{~0v?K?2?Jnh>y z%;K; zbtE6ZDf(%6xfE)t%qh#+|AgJII@7|%xywoz}MSl&YolcVFTbpud?mLN;w!W`5ftr1iD!xhATwk-N`)(6B--ydS| zvjgIXn$GWD^&Ih?K5H7hGDZ9`T`!g#TZH(XtJRKJIH3CCvezr^V^RIdz&A4_koa~H-+bbOo`!wb#E8$H_`pm1Oo$Kq z8tro;z9iy%Kz#d%&yV;_INv{e{wS&9Q4)$nhKg!z@AwmyOZdEOA9`3ny1&#pcj#2I ztCA#Szl|j?j^DC)K=BHPjF$8{3Ej#t$ZPEt=vK?aD8r-Bt)58lx}in!F5jx1{z-=7 zU8L?_>l;Av&RV6viY}vgF)f!?mK>*eN5ldLm`$U2p)ME2o9+3fjhai8}4tE9`rG+G+a~EBYXaY z+oE;}pI_){YebD@2ZOYAgHYp3GcRAwMAX>4?z(G&30&9ou4mfTG9>61_-L(n4HAr8 za(-Lu&q!dQel^Th2MHY3@6br=iUgiYisi*8kzilU^W$^J!F5fOb0V|fBSH4e*OEh| zkl^0)r2#FeNN|4oq^jIGsPS$1KC!cVQElgfyn>cnZMlSxGd`UN#y8ZM@tw$HeB0r= zFYMdoPWyZ%9y7jW#~5G6Va8VoeG2>LEhWAz#y1$Q%fh}^=wY|jXa!W`<=DhkNnh+KmBS6{pzv(qJb78(9vTL)#a?SQ2egh zru!z0LuttcawVtzQ0|PpsD+nxQK{xJ9mkC~{?M<$cXrO6EbPPMzxwbt?c+)0)4n*> zp|sETeiQ90AEitCWChC^U!RGz@0?K~?Q6X{f$QwgLsaLNj=nyG=VNZ`g=s>KSrD7qn1Ar>EOo3_2nHBSrVE8%`tsMP0=u z^>E!)aNbGETOaWkx%p3Qfo?V2<4FIhzaW06kt2G<)uHIvn_qx9Z|ho zgY?5vc-%PbWQi2?tlq812M4}yPssAAw&4>ezEi~anB{|>Mf)xhUp?{dX8GW8&IkLK zn?G`3Qc^|RQYwzyiE3d|tK>J*<11HCU&H;*~Aj*?*zNx3R?FM5}2l#VhdU zjYyZLczGcMj>yGOyxg=G4`%m*ex)ammQc{IdPnYzy-D$s_YG9uiYQ)u%*E@g^r2tP znCYyfNbw?PUwW+_OYuTdeDxK^P`quQ9>|(rr+D_knpLwlK))L1e=694lJ-iQ@mb9g z`jyKRZ3msqn(*9%;6Uyj^r48+`YA`C56x9dyk`%6D6U}f9Jv0#o`@{zw3;T>*iLVp!JtIE2FOc>@m%={FI^w%ce7%Y9GVw_h-#p?gB0hKG z^Jjb5+nkTyzki$g%ajip;&*L-x2w&}DBRyhWE&~E-->>yQ;h~4XigP-V~@t%jtwr3 z$wp>Bgqv1P8;-0$ENuGGbO&-uEY_OvP7QfDCxyvhUXA<)%63)^{TYQ+-59ks_B=Wo zwd(26_;lz)o{ua%uR$NWn%mTz41FlF_nC?!=tFZQ2d{(Ul@|Pd{~+x<zV5^qMtrrz z2iI42a3BAz`_j`onb5E7=8arFx$W1#u)jON{B?96w?6Rt<<^I&&MzH(eF)FT+{n!U zHv{d@0BTrN^Lpy&O4MLkv`UzqBXbTAb0b<@`8*t|ydUpH26OZb55^XIMF{1aTybc}Bi zJ5m$+mCo=jn>V8RN%^-gG#5g@8Yivq+R&bm&8ajU^gK_j%8N zezp9!msY9Jue`hXuNn{i>g4V5$KZXgvirVzsxS%qmGtT#f~G^i+FEyWo;CEV&~@Vk zupKjAe5|Ys^ee-~;WqF-Ss~GVZ>v;7zk1>967&=Fs|}ZLy1j#b)qARMh%NN1Q$x1) zZGnC@ccVdzYw3p@%(++9iB`wu`+D{pHIYW5E*WLvb}- zcZNbADj%Cz{R;Y!Rzi#>ygu!h(apg#0Qyj5`i;kt(1+~uBT*OVLy0*(x2=FaWEC@X zpEmTN=Vg-i3DAdj=&jl~3i?p^-X^WF`p)SMIweQ0VQ z&*DvN4||*V7874M@jWFzP2z+1+oF9@#K$AP)x@Vke9+ZsAKZ7y`C$KY^T#^iGUX2y zFRv1pyV}f*!fzXqZKSxt&;G{4QD{(PRDk{S4`|Fislzv4>LD{FmD#1Y1CVvoq(Z^O zYUFgm?)t8e-JlOOyWVXa2YpEIa90m^=tK36>Dpn?hrE*p?Ja{ov~`x2lsfdGjL^D( zZP16p?UE(SpbyRLzuE?_Tcid5(yzev)KiA`!T6mIwrO8Zl0NNoh!SIberb&F^CQML zcrN2>sbhS62jYt(zK6tjg81}^&zJa0i0@K|_wnB@z2|5=0s7S}t)|KGZNL6y=2g`5 z*U^34`k>b@w?6RsPXyUncP-6QA4v(5FRwDx8nrzyH|zgD$kyHfqr#D*WW7x`jvHhd4eqTtK{zKCKH8T=rbPg-7WTudK7xt zqo`fN=V!av+t7!8tF(`mfIbx8x+V#(b84=zI@045^dW8ZObYJ9>z8qoHwUiA8hOMf zv;eMavRLTWI2G=Xxcn?O{K5`jfFn6(lp7r z40*Un?(ROZ0{W1Cjqe0W=tF{Y;xVhB5BUUDneE4Y=uVWO4(>y9islDG9|}vqb*< z&T7>K2^qNsLJ@pzbC%c z#1}z)L;tcq{#&P@8x=30U(I-!5#ZkT>tAMGMLmBV-N&sDdi`?igFgTL`RA98zCMKK zV{YVTfSZB#X8`(8(eUHq(1k3sZnP$zhdz}0E^Z-oEVa)@Zv>{$hgA7u74ZI26*FEe znsX2OmGg=ZHoc$^Nn}d*34lHnZE@2u8v0Ol^Pwg1zE)ostKYY|1O19x6%~5``c<)H zZtz{`SBhu)CoYG6l`55f+#C9ky1e?jg!Y6iAKceS`@l>4%!m&<+IPO*#OF_ZF~qmF z4c|TDTSk2FzVn<9_U}J_{%{w%WRkmDRAYPVgSb|TYGnIxze?YF9#w7uYOtCcLMHLEzRWnbhrP}DxcU1tk3V>Q z-ZFLUx$e-1=uvn)LcgN7p>HUM?Kj&FeW*)Pe|s6|L$mi7>K8&EGJAYI$`$(1qRNXt zQ=kv6s@!8O3w=mxP4Vm#(1*q>Idpgl^r1JJJ=b)_eJH~vb}#NjWrp>qaUYt_^Tzk* zg`FAh>WKT*kYRWKUB3cH!c)eV@br6M!c)d)62tf$ycpl;uf$hNe4U9;gZQF|&x!aH zh%bcrPW@$l(65{)4|;|B)w$&}eg2!*zx2F{dj2}Pk6Rz~`sLOKeg6CN&o3Q)eF)FT z+{n!UHv{d@0Q8~4rsbWWTUlo9-tV&p_o1{GZ}Ig_pC!}#48(mX$=Ctz53HCG*ZK1- z+^>c;hZW&IwB1Q#9PUGM?zVezAG*3L82)|1uMX8bx$U@L>8PpX_TLVT9QS3rE@S-va82k*~L`!tB}B=H?1zTw0- zl=FSJfB$jwhr3YXg~lGD8rwVmh^w2ZMz#<4t5g;34Y*qkcTD^6Z7dmJ{=}pK_p1~8 zSBK$#b;-pS{(ZvSSPxA(=v~=;_Fhk#kNZ`Re5>xbUxm)`r2Ccf_Lc*X^H!l z*vJ^z4vD(LI|uh2ZmX6EOoq{Z-2!Eu(YRkV-2eC^?pK0JRE7IhVE4jS+^_QQJb5Pc z;y#nRN%n%MN1=y3irOW7ezuFfjr&k*%m!!Nuhf)NKI4A1;oSo0QVQCZeq(0gew8{p zz7OtK%Wbo_<9_AueAWo}E1%OV;QljD>!nqL)o{Ox9=_BZ`c;tlm|qX#eiasa8or-5 zY{a}Nhj72rO;%9E{VJkB>?H12!9|}O*&g;b=i}z@&piGzM}=?pLRb-N)d5r4@6$DBd~2H+-^-80w}besh_5&C!TZb8zOlq-&hqIIpT%F+2mQ)%ar`LUug=VQ zRnYe9UuIrKJ%1hD$E}a=)-Sg{nDgJCc7EyT>qB@x=0^GV_p2>g)h@UXNv8jt zkNef7YkBi=A9^&T7CKfwr8(H4Kkh?U{P>G-A97Tfhj1Utcqsi_+n^E&@cjmSdzWp0>yBRa+Ut7VivhXGCywm~Tf&W?fFZ_G( zV@Si70({~7TU*;MUv9D7$=YhQgVRd51OCFlr_>wY#$|{--gI!83BOU2Zef(1QFq!7SZ_ z{;rhFU6{>?=NI(L39;Xd{9%~!zgmADz5es2n#^JT=D&IUTeZV}^Iz8AW)RS2);hW1 zm)0-%`e{D{ufOSEUw;wrUnOUT{^M7Q=lUYe_k*FR<`2m~z9lLb`?&vDB_&1iR%=I0 zTrgSGvvHqjpVy|>Pdt746~=m&Tgx};*LnA-ABoG@hr=A8!W2X)cP^5wLYSnTpzkGQ`1ke)_{`lI%&00 zbju&J%Ybg{Ww%iB9epF)t1Y19+r{ggt2j!@%fH{Ea{ekMFSAwo(u&_HUcG03kb1dIFVTuh^QW&OjBymzN~C1$_P*>{oR<(2M#K30w5T^(O2j5nZoN3%1U zY4}q-{TudO+=D5(B)e_+3R+ALAUGw?Bo13RX+pe8wQ5+7XS! zxh_aJEK#S+izFnRGCh8AMhFsGYcI~<_znp@1c}S)JdrRkwd0+UACNG%UP*qB2@)oj znny>ULBgQ&B2}YVsOqXxyZA6CRQ`0XV|hH6^QU9lXJKMV`}(Z`C-}76E}(s~OX2as zx3+f|+IJ!09PR6rx`Xyv?#QEk4~E^QeU6uV(!SgeMU2mF72~@T0R;WpM-R}-fXEC za089&B0qSe@jEp0Rg_l`-J!_p;}FB>fv?dj$Gyv5?757*!)Mf_Hs6B$7G4i{R5ucZ z+6y*+J$x6%ZSc9X??nhozNE6{yp0pew94GpYepg}q*Bkl+Xk;^TF}t1z*p|v`b2y1 z(en>J=KiXfj#5+R9sr+J->geQ@ZESKm>G$Eqy27sVV}p7H8I#X!F_xS+V^bOV%isU zza91kcG$Xv_Pv~6LHiCkJf?jYtn+AJ&x|9C&m8u-nEM1O>nqp$$e?f1cKzkM9&4-* z!SgE~>&N}?3#6}!jZd;aB>O{@=NB>GAME*9G9@b@S%Ick013jb1Y8;6g#`P`T*q0h zKmyNkv(D9*Ac1xBjJ&yvkYIdT-#PR9Awl;8r3vGnBK{|v&s25<;^$;sIG21G@sDY# zELpw*@mGJ)tL{7#@r^c}urpCd{7%D%6~ro|%GZZIZ$|i|%A}ayQSGOp%3TU;Q=PL= zrERy7i>J+QN=SS`jBiOC^@wFzt2aNB;J>v5t zzFEW<$oQrXmiXxB_rG@is1_yE`u0Z`i)uFY`XjETqMBSE$`fSXK5?vs;*~ygJ(bw; z$9$4{{)wRi#d|EHBIU13@g8=`8X7s9;@ve4Pj0u6;$>y1r<%u5yo^&;>7IUP|AK z1NsRlUIJy^f5?4`7c0<;TH_0y>O-wtNontAU-SFRqR+lKA z^o@{#ieD%>@0XztZBM7?t3>OYJq%1QKkQPr^z?bwOUY@64Ia9ugQzF#;Z9Mz*vIF( zxW}mKp~_$-b2n7=JT36So=2$abJ1^(vcDl=t1hQs$vYxpx6ns+-3pO#kVAc(o+1(& zDHjKhRY$@_W<`33bdYe(vnbW?=aJBF{_Z!KeULEbCGzNRi-hS#J8ey%Tissp`ObJX zBs`fMQgT-TRq_3+hIHLx^>B|FpN%Qwv!@u}_-?fCW!XN)H+~Z1Q&C}j z^DG#jPXg`hu5gs`U4`#spY{pHH|q}LqdG9YJCTgf+kx@jozL}fk0n0*{Ql>xUwX?& zi|uXqO-%dDO!oQKC&#EmgjD%Vv<7Jm*|X(cu__w=EO6(JPv&TXS0BYe!7|7+G=0v^ ziVetmL(ULYxgE%5(wde>Yl_fD#ZQjcd|o2IyQQDD##*A#Ln^P2+H6K~+n@L69e9nB z&z(zl(VK=cP0ID<60V^_Ew>ql-`;Th5bZOTYeoBD{^WyW+GkKa2m7|F6i=poa}T=H zzD4t;XKA%0swC~uOp^R_01LM>1Vtf~jXy1h9gT>y* z|0Gk87FUMqpY$83Jwlb+uf%Qr#N8*^AFzKV`$J^s7jfSo?D<48B`Y9Vfu>gg32djD zU9odVf_eG}+IC4q0>km+Z_R-oq_sCFa$-0V$lVZbjPHl|FIxF&T=78sl#yflI37iO z|H|UE9oiwj&9QwZPvQ_?H>Hb|8bW+If8W>lzoW{6?K1}{97mPWBNna=qEO|!-LFq) zZA6umjjujbTa7AvUK^d0mEV+*_`Hel67jVrz8K=uA-)~N7fO5@#1}<;ri|}FGvZ4i zzUjobllV3f-*c|do7}-pQTtd2I|;90h}zcxY2=2u`=}NVXPP?iixSmrWZ4?-Ca!2v z%?1d$f$}PmjrFC?6z@@=LVw}1AG1(eaeK#&6z_p?d~0Vrig#y$g6-o>=vG^L#+GVB zx2lhf41j)>@_K*K*D8v4Zt>RM{d6ea$yJ+6zm27MkzMz7f4Y+5g{ZCCI_@*Y^Y=Y{ zYIO&S=bgDE!V$i2F*owrWr}K9_oe1lmL(b@0BX3Lnqg zsYuw}v!z!|G!hP(zh5tAC=yQWV!p&H1_>=3Z}BRgBB7gtRHwPtNVqrH`AO7nB#awX z&0p9G2{Q~YsC0Hf!knT~nLbmH@YJ~73e&YvRkc^td$UjIt!#$vZ%5y7J=|mB>(2OG zn=!tw-x*&HJRbHPzsUGv;c>CgvmN8(L$AWVIRV61&G;64CBDhTx0K^6CB8rX`$OWx z&+mWk`lT!AgzJtyo|RU|O!oQKC&x$$b*(wn<|Wc}(Q#8ca3ATOxo}@~&;Vp4quFt~ zzA`ecY}3c^&~RiO=e4EScMNh_az1IsoN}~Lef9bApBEs%CqrwkPH#k^VOL80XUC$r zozs`e`rkpx358!vKNh0QX@$`pV!ctJM(jOhKc_$SA?#cE#*y~H{K*H$w9jX_GwoBK zcZ>FQ+NeeQX0!-reCysYJ{>K_7oW`dGTt%1XW5L;uruS^w~O(8DkZ*A#Aht_K7Mtj zL64jVsJ^IVVMdsLqd))R=T{8YPuzWy{ekzdWPgb2{37oAgFPQhrep;qE70@`Ai?4@ z+rF1&Bf+fDmu`k_kU;N%ZT6Q0Bv3C};raa-5=ed5C_ER5_{Br4m6KZ_{<*1h-)q_- z{&xPjA*#0#f6<D=b7 z@}{S#a{RgLcWs?eW!E0^Z*O`uCFJgx4@c?Q4KEazosGREvDyl_edFL^T^(wuZZj%UD#i0YYw|yowsn zcxo!eE3#6XVFaJ=V`stGG55XTx>)zUk1Xl}9V#~Mr^eMk{Ft7vw5*e_eQbJpU=DRvz6a~2vBBeP}9gVE;DIhl<}He((hPkmvUn?KPke?b2Iw65gkk z+iUL?_2tlqb~*I(tVLDt+wTY%pp4$ijXBq20=ZNpSM}qV_$CrxH1TaAz9{0GPJH&n z=S+NWi4Qs!?R!XkTEqvRC!~GwJ~P_)h4?y2eDw4CU$=hg3cD>0Mn_FZt79hnO4KLE zNXaStEHdbXG=rw^=MPpyy0@*eJDM|?TjAup^&h%JElV) zn)9hJ6g~$_3mQH@Nc)CMo6$a)Kl$L8_KnMRq!sB7uFuSLjzWA1PhmQi=p)^BiW+=nj3T!)y14 zhoKMY-ZPb5gZPh^D(HmCAbx^Z;k@&z(1*IamFmbL{=(FmJZX48)_^p#-&EeC%I_Z8 zMH=H#<(-8N!Af6IWmvV@hIkiL=`?t9(h6HtIqt#xwJ$YMrHcNT=%uJBA;&EdNX~b98gpcC-yvZHx6t$0au#@l_hNyiFkVbBZyAQh1@*5u7R-&4X zEL+3f#AP9>*#IFopkGxKULOP3rz)DPV;QvR$1Kn=GR^gae&rCEmzM|qYW;;%hH%}g z%xKr+{T!fQb>ADB4;?F|WvA}fA3(qI4AP#h1^sHT>gV10(62_vo5eMQex%&t z>R49ntfkPe%Gw83S3tj7YhSuR8~T+>%~Q>R(66FzkL}YL`qiA~nv-@zze+5QUscO` zq0g+ow=%at)RXmar>I@*<8xiyW9UO))~zg@0(~epka8FaeW>VpYZD*nLwz<4a6Juu zDCPOOL-(N%1@AQ3l?HuiLD;+xbD$5M_W%41I+xq*l`U<3p%1-0nq>u_!;9PUSz(bY z^r1knhe;LChx#i~gW&Ufr}no|Js<~tsHe4k_YTm9CY>MToyPTWk2$_k#J7p~;PZ(; z`O=BcocN$;(LVU|{3TyH*TX%Q`2Okh`^&9gy28O-_x9ATOsiui`+Vz@W2CfZQBoLx z4AKntc_cNf8`8~wepA6l=kc4H9qp$#Kz2CjlWbSb{q_A2N@f%~>=KZ8Cr-R_kz8vps1Syxf7pSb%Z`-9%UlKnxS|Nj2-i@5I(_Ixavk`<7w zK+`J#{h*-KFncugtE=wv0eaAfQeq!}vqXX(U%&b%^o2g8cDfmQ0exu3snyGsp${z` z^Lq_+qV2K=bNahOzl!Qt?qULcDC)IYD=p|()*WZ)sYAaqu>ZVdKJ=j?Ud>pTPMywI zw#|ormC`jz2%pze>D{-z?CYk4jrgn@@pU9V=xYC#&xH8OxIS-k2RlXWV;$@yyoMoa zUjw9(8{+PRF0`yw|0KBnm1t%o%hqr=aV-+nY=Dp((68PWX-hqZZe_f+=dB(;Wz>9zzk0Dq^-Xi=SFQ=i*JMGz>TK+}6t3Uy z9%Zs6f!Y8`P8_p9D!HiEIx zucnw9m_>0t+~Y=k$0R<<`u%^r{?Zi=<-h2??{ZolGuhXpJ~>88Ynwj|3Yiam=zEcK z=r*LAV|1kU{XFPH)QO%s^PvyT*L)s37y6JxQYOy}x!B3*IW4k-J~aGtmlxfj4+(u& zt%`s?nY#N+aUU|@RC5*jP(a2A-@eeVrcG5-g6qrCft?AjH?_;v;T z?AsN@_(H!hz7QYASF?`sWtS3P0P!s)zU9P+`_&(P=$H5LtF40yugOEdn)Yno0m%w9y#mk=@{LSlIzhj>Tzuu&0q9rB z^8S;CLLaI%Z*g=2^ea`FDdvlzUrnof9<~MgknP44>E+Oeqzl@`w1R$h%u0E>CiJU_ zj)R0_pgpN*PZB-$K7iQcHWZ9{P~Fre$1{z<<$a z(S#4KgY{2+o&QyzHxaW_)IQe1PQq&#qV_dF8o43vKIlRY7gvO3i)uEqYz=o4SB|J= z1BBeb{VF9IUr)Wrcx0Ys*pFGzbLK3IK-{lR$;oR#zjD9)#L^V^D~}`B;5zRaT@MQHHx7lR_aKAeA_*H-0uX?oi4#E8@{(4z3?pN}+t9wDe3cnGv!3FoL)_vs=?pM}H zk=LMKb^5$O_678-$l8+iUbtWFs<C-CUF_p? zUEE{bhidbrW*;O}sEACg$wS{n9tJYf6-}~m`ewDZ{8Q#w{B*Zp6_XzG+-T|wi zPg$JtxG2cO{mSY=?s43&!evWu;(qn5YRw?rul%+9UB&&Xe}R1>?pNBqYj@#(wdd8e zi@0C;D)wsn`iT;sWc~ie*B|_R{O;k65nRxTg&miA->taybt=7y?W0Q+^UD1)%*>0 z`M6(IjYzwO`_()BQQo*;ZML&c!~N>+Q19KW7x$UWOS;vfo~(yEMeSlApX=ft<9_wc z*;s)4mD=;SmvO&ZljfC#`_)Xl_WrnEC3TJ$;C{7qOg8-g`X(;)F^+-jMp>Nk+>ZL< zer1!E6^{GWvDWGDpkHlku2@iy`&Gbl*M7KP4I2E>8~3XLySv@O{c8Vq=@i_rwjW$^ zitFJXOMH^``#WEM=?cM8CgHte)9RSXz8>|-G4zej%FV;~LLX|IZ>=>3`p`_1Xhn7C zL&kHrNR5C#v{-ZZvN-5N%jFI|vcY|5%S6FU+=rZ$PYQ9r8tkhz4)>uems;rJKJ+No zVJPlbQ|7y_!~H6HNrDgVSA7!9rT*NnXqOD}$#8r!jPHmyvR3lM{#_^epw&%E32@GVBD`x1l*p}_|L!0x{7-J#N8*^AN2l}><{|<_xGP)#C?CT z=VQr~tbk+%nqC3i4}{yzH{gD?<@dlAxDO?K{yG-=)#qJ_q8P(3H=c Z0NE*OAM0Qz;WZ3V`x+pP+z>b4{{T)*Mt%SQ diff --git a/tests3d/test-equilib/data/expected_equilib/n_eq.dat b/tests3d/test-equilib/data/expected_equilib/n_eq.dat index d0c4d26b5a9319fe6273b3da82b6ca35e153b654..74f8a32d6fda6aaabc492f82d89967742bba3df1 100644 GIT binary patch literal 160 zcmZQzKn1C#i~s$;`QHBLRxyP&&wzB)=hEAM-rM_|H3n+%ez2EbZL(vJ+z0!nLlGIX zjXv1htx0W;@%Uh`nz^EGWBdnu1`4tl;Pd)U( reuihbf$q5v_9DES^`@==V84v<@>G@2@9h(RTcB@Z@)fwMyEIXD|=MeFfafBV^~Dw diff --git a/tests3d/test-equilib/data/expected_equilib/phi_eq.dat b/tests3d/test-equilib/data/expected_equilib/phi_eq.dat index 4e77160..49066ab 100644 --- a/tests3d/test-equilib/data/expected_equilib/phi_eq.dat +++ b/tests3d/test-equilib/data/expected_equilib/phi_eq.dat @@ -1,2 +1 @@ -ʪ>?dt/G @dt/G @FH" @Rf @̷y @ ->r @* @5 @={ @~ @M7 @; @Xq @ @}Zs^ @ľ7@9>y@cfTտ3в \ No newline at end of file +u(?a/G @a/G @DG" @mf @fy @ %r @YÎ @S @VxR{ @ @Q27 @2; @ @ @4,Zs^ @@?y@Tտϲ \ No newline at end of file