diff --git a/.ci_fedora.sh b/.ci_fedora.sh index 77e7f45055..5f0a36efea 100755 --- a/.ci_fedora.sh +++ b/.ci_fedora.sh @@ -41,8 +41,9 @@ then # Ignore weak depencies echo "install_weak_deps=False" >> /etc/dnf/dnf.conf echo "minrate=10M" >> /etc/dnf/dnf.conf + export FORCE_COLUMNS=200 time dnf -y install dnf5 - time dnf5 -y install dnf5-plugins cmake python3-zoidberg python3-natsort + time dnf5 -y install dnf5-plugins cmake python3-zoidberg python3-natsort python3-boututils # Allow to override packages - see #2073 time dnf5 copr enable -y davidsch/fixes4bout || : time dnf5 -y upgrade diff --git a/.github/workflows/clang-format.yml b/.github/workflows/clang-format.yml index 87f1947802..49dfb31e25 100644 --- a/.github/workflows/clang-format.yml +++ b/.github/workflows/clang-format.yml @@ -22,7 +22,8 @@ jobs: - name: Run clang-format id: format - run: git clang-format origin/${{ github.base_ref }} + run: + while ! git clang-format origin/${{ github.base_ref }} ; do git add . ; done - name: Commit to the PR branch uses: stefanzweifel/git-auto-commit-action@v5 diff --git a/.github/workflows/clang-tidy-review.yml b/.github/workflows/clang-tidy-review.yml index 8ddcd4ae3f..1a9a0ea079 100644 --- a/.github/workflows/clang-tidy-review.yml +++ b/.github/workflows/clang-tidy-review.yml @@ -22,10 +22,9 @@ jobs: submodules: true - name: Run clang-tidy - uses: ZedThree/clang-tidy-review@v0.19.0 + uses: ZedThree/clang-tidy-review@v0.21.0 id: review with: - annotations: true build_dir: build apt_packages: "libfftw3-dev,libnetcdf-c++4-dev,libopenmpi-dev,petsc-dev,slepc-dev,liblapack-dev,libparpack2-dev,libsundials-dev,uuid-dev" config_file: ".clang-tidy" @@ -48,4 +47,4 @@ jobs: -DBOUT_UPDATE_GIT_SUBMODULE=OFF - name: Upload clang-tidy fixes - uses: ZedThree/clang-tidy-review/upload@v0.19.0 + uses: ZedThree/clang-tidy-review/upload@v0.21.0 diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 90c63e0233..4153ebaefc 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -208,8 +208,8 @@ jobs: - uses: actions/checkout@v4 with: submodules: true - - name: Build Fedora rawhide - run: ./.ci_fedora.sh setup openmpi rawhide + - name: Build Fedora + run: ./.ci_fedora.sh setup openmpi latest shell: bash env: TRAVIS_BUILD_DIR: ${{ github.workspace }} diff --git a/cmake/FindnetCDF.cmake b/cmake/FindnetCDF.cmake index 393c57549b..361095954e 100644 --- a/cmake/FindnetCDF.cmake +++ b/cmake/FindnetCDF.cmake @@ -32,6 +32,7 @@ if (NOT netCDF_ROOT AND EXISTS "${BOUT_USE_NETCDF}") set(netCDF_ROOT "${BOUT_USE_NETCDF}") endif() +enable_language(C) find_package(netCDF QUIET CONFIG) if (netCDF_FOUND) diff --git a/cmake/SetupBOUTThirdParty.cmake b/cmake/SetupBOUTThirdParty.cmake index 10942f8aa9..53ceb4351c 100644 --- a/cmake/SetupBOUTThirdParty.cmake +++ b/cmake/SetupBOUTThirdParty.cmake @@ -272,6 +272,7 @@ option(BOUT_DOWNLOAD_SUNDIALS "Download and build SUNDIALS" OFF) cmake_dependent_option(BOUT_USE_SUNDIALS "Enable support for SUNDIALS time solvers" OFF "NOT BOUT_DOWNLOAD_SUNDIALS" ON) if (BOUT_USE_SUNDIALS) + enable_language(C) if (BOUT_DOWNLOAD_SUNDIALS) message(STATUS "Downloading and configuring SUNDIALS") include(FetchContent) @@ -293,7 +294,6 @@ if (BOUT_USE_SUNDIALS) FetchContent_MakeAvailable(sundials) message(STATUS "SUNDIALS done configuring") else() - enable_language(C) find_package(SUNDIALS REQUIRED) if (SUNDIALS_VERSION VERSION_LESS 4.0.0) message(FATAL_ERROR "SUNDIALS_VERSION 4.0.0 or newer is required. Found version ${SUNDIALS_VERSION}.") diff --git a/examples/elm-pb/elm_pb.cxx b/examples/elm-pb/elm_pb.cxx index b8a8bd41b9..d981796200 100644 --- a/examples/elm-pb/elm_pb.cxx +++ b/examples/elm-pb/elm_pb.cxx @@ -5,6 +5,7 @@ * Can also include the Vpar compressional term *******************************************************************************/ +#include #include #include #include @@ -14,9 +15,11 @@ #include #include #include -#include #include +#include +#include + #include #if BOUT_HAS_HYPRE @@ -36,1860 +39,2131 @@ BOUT_OVERRIDE_DEFAULT_OPTION("phi:bndry_xout", "none"); /// 3-field ELM simulation class ELMpb : public PhysicsModel { private: - // 2D inital profiles - Field2D J0, P0; // Current and pressure - Vector2D b0xcv; // Curvature term - Field2D beta; // Used for Vpar terms - Coordinates::FieldMetric gradparB; - Field2D phi0; // When diamagnetic terms used - Field2D Psixy, x; - Coordinates::FieldMetric U0; // 0th vorticity of equilibrium flow, - // radial flux coordinate, normalized radial flux coordinate - - bool constn0; - // the total height, average width and center of profile of N0 - BoutReal n0_height, n0_ave, n0_width, n0_center, n0_bottom_x, Nbar, Tibar, Tebar; - - BoutReal Tconst; // the ampitude of constant temperature - - Field2D N0, Ti0, Te0, Ne0; // number density and temperature - Field2D Pi0, Pe0; - Field2D q95; - Field3D ubyn; - bool n0_fake_prof, T0_fake_prof; - - // B field vectors - Vector2D B0vec; // B0 field vector - - // V0 field vectors - Vector2D V0net; // net flow - - // 3D evolving variables - Field3D U, Psi, P, Vpar, Psi_loc; - - // Derived 3D variables - Field3D Jpar, phi; // Parallel current, electric potential - - Field3D Jpar2; // Delp2 of Parallel current - - Field3D tmpP2; // Grad2_par2new of pressure - Field3D tmpU2; // Grad2_par2new of Parallel vorticity - Field3D tmpA2; // Grad2_par2new of Parallel vector potential - - // Constraint - Field3D C_phi; - - // Parameters - BoutReal density; // Number density [m^-3] - BoutReal Bbar, Lbar, Tbar, Va; // Normalisation constants - BoutReal dnorm; // For diamagnetic terms: 1 / (2. * wci * Tbar) - BoutReal dia_fact; // Multiply diamagnetic term by this - BoutReal delta_i; // Normalized ion skin depth - BoutReal omega_i; // ion gyrofrequency - - BoutReal diffusion_p4; // xqx: parallel hyper-viscous diffusion for pressure - BoutReal diffusion_u4; // xqx: parallel hyper-viscous diffusion for vorticity - BoutReal diffusion_a4; // xqx: parallel hyper-viscous diffusion for vector potential - - BoutReal diffusion_par; // Parallel pressure diffusion - BoutReal heating_P; // heating power in pressure - BoutReal hp_width; // heating profile radial width in pressure - BoutReal hp_length; // heating radial domain in pressure - BoutReal sink_P; // sink in pressure - BoutReal sp_width; // sink profile radial width in pressure - BoutReal sp_length; // sink radial domain in pressure - - BoutReal sink_Ul; // left edge sink in vorticity - BoutReal su_widthl; // left edge sink profile radial width in vorticity - BoutReal su_lengthl; // left edge sink radial domain in vorticity - - BoutReal sink_Ur; // right edge sink in vorticity - BoutReal su_widthr; // right edge sink profile radial width in vorticity - BoutReal su_lengthr; // right edge sink radial domain in vorticity - - BoutReal viscos_par; // Parallel viscosity - BoutReal viscos_perp; // Perpendicular viscosity - BoutReal hyperviscos; // Hyper-viscosity (radial) - Field3D hyper_mu_x; // Hyper-viscosity coefficient - - Field3D Dperp2Phi0, Dperp2Phi, GradPhi02, - GradPhi2; // Temporary variables for gyroviscous - Field3D GradparPhi02, GradparPhi2, GradcPhi, GradcparPhi; - Field3D Dperp2Pi0, Dperp2Pi, bracketPhi0P, bracketPhiP0, bracketPhiP; - BoutReal Upara2; - - // options - bool include_curvature, include_jpar0, compress; - bool evolve_pressure, gyroviscous; - - BoutReal vacuum_pressure; - BoutReal vacuum_trans; // Transition width - Field3D vac_mask; - - bool nonlinear; - bool evolve_jpar; - BoutReal g; // Only if compressible - bool phi_curv; - - // Poisson brackets: b0 x Grad(f) dot Grad(g) / B = [f, g] - // Method to use: BRACKET_ARAKAWA, BRACKET_STD or BRACKET_SIMPLE - /* - * Bracket method - * - * BRACKET_STD - Same as b0xGrad_dot_Grad, methods in BOUT.inp - * BRACKET_SIMPLE - Subset of terms, used in BOUT-06 - * BRACKET_ARAKAWA - Arakawa central differencing (2nd order) - * BRACKET_CTU - 1st order upwind method - * - */ - - // Bracket method for advection terms - BRACKET_METHOD bm_exb; - BRACKET_METHOD bm_mag; - int bm_exb_flag; - int bm_mag_flag; - /* BRACKET_METHOD bm_ExB = BRACKET_STD; - BRACKET_METHOD bm_mflutter = BRACKET_STD; */ - - bool diamag; - bool diamag_grad_t; // Grad_par(Te) term in Psi equation - bool diamag_phi0; // Include the diamagnetic equilibrium phi0 - - bool eHall; - BoutReal AA; // ion mass in units of the proton mass; AA=Mi/Mp - - // net flow, Er=-R*Bp*Dphi0,Dphi0=-D_min-0.5*D_0*(1.0-tanh(D_s*(x-x0))) - Field2D V0; // net flow amplitude - Field2D Dphi0; // differential potential to flux - BoutReal D_0; // potential amplitude - BoutReal D_s; // shear parameter - BoutReal x0; // velocity peak location - BoutReal sign; // direction of flow - BoutReal Psiaxis, Psibndry; - bool withflow; - bool K_H_term; // Kelvin-Holmhotz term - Field2D perp; // for test - BoutReal D_min; // constant in flow - - // for C_mod - bool experiment_Er; // read in total Er from experiment - - bool nogradparj; - bool filter_z; - int filter_z_mode; - int low_pass_z; - bool zonal_flow; - bool zonal_field; - bool zonal_bkgd; - - bool split_n0; // Solve the n=0 component of potential + // 2D inital profiles + Field2D J0, P0; // Current and pressure + Vector2D b0xcv; // Curvature term + Field2D beta; // Used for Vpar terms + Coordinates::FieldMetric gradparB; + Field2D phi0; // When diamagnetic terms used + Field2D Psixy, x; + Coordinates::FieldMetric U0; // 0th vorticity of equilibrium flow, + // radial flux coordinate, normalized radial flux coordinate + + bool laplace_perp; // Use Laplace_perp or Delp2? + bool constn0; + // the total height, average width and center of profile of N0 + BoutReal n0_height, n0_ave, n0_width, n0_center, n0_bottom_x, Nbar, Tibar, Tebar; + + BoutReal Tconst; // the ampitude of constant temperature + + Field2D N0, Ti0, Te0, Ne0; // number density and temperature + Field2D Pi0, Pe0; + Field2D q95; + Field3D ubyn; + bool n0_fake_prof, T0_fake_prof; + + // B field vectors + Vector2D B0vec; // B0 field vector + + // V0 field vectors + Vector2D V0net; // net flow + + // 3D evolving variables + Field3D U, Psi, P, Vpar, Psi_loc; + + // Derived 3D variables + Field3D Jpar, phi; // Parallel current, electric potential + + Field3D Jpar2; // Delp2 of Parallel current + + Field3D tmpP2; // Grad2_par2new of pressure + Field3D tmpU2; // Grad2_par2new of Parallel vorticity + Field3D tmpA2; // Grad2_par2new of Parallel vector potential + + // Constraint + Field3D C_phi; + + // Parameters + BoutReal density; // Number density [m^-3] + BoutReal Bbar, Lbar, Tbar, Va; // Normalisation constants + BoutReal dnorm; // For diamagnetic terms: 1 / (2. * wci * Tbar) + BoutReal dia_fact; // Multiply diamagnetic term by this + BoutReal delta_i; // Normalized ion skin depth + BoutReal omega_i; // ion gyrofrequency + + BoutReal diffusion_p4; // xqx: parallel hyper-viscous diffusion for pressure + BoutReal diffusion_u4; // xqx: parallel hyper-viscous diffusion for vorticity + BoutReal diffusion_a4; // xqx: parallel hyper-viscous diffusion for vector potential + + BoutReal diffusion_par; // Parallel pressure diffusion + BoutReal heating_P; // heating power in pressure + BoutReal hp_width; // heating profile radial width in pressure + BoutReal hp_length; // heating radial domain in pressure + BoutReal sink_P; // sink in pressure + BoutReal sp_width; // sink profile radial width in pressure + BoutReal sp_length; // sink radial domain in pressure + + BoutReal sink_Ul; // left edge sink in vorticity + BoutReal su_widthl; // left edge sink profile radial width in vorticity + BoutReal su_lengthl; // left edge sink radial domain in vorticity + + BoutReal sink_Ur; // right edge sink in vorticity + BoutReal su_widthr; // right edge sink profile radial width in vorticity + BoutReal su_lengthr; // right edge sink radial domain in vorticity + + BoutReal viscos_par; // Parallel viscosity + BoutReal viscos_perp; // Perpendicular viscosity + BoutReal hyperviscos; // Hyper-viscosity (radial) + Field3D hyper_mu_x; // Hyper-viscosity coefficient + + Field3D Dperp2Phi0, Dperp2Phi, GradPhi02, + GradPhi2; // Temporary variables for gyroviscous + Field3D GradparPhi02, GradparPhi2, GradcPhi, GradcparPhi; + Field3D Dperp2Pi0, Dperp2Pi, bracketPhi0P, bracketPhiP0, bracketPhiP; + BoutReal Upara2; + + // options + bool include_curvature, include_jpar0, compress; + bool evolve_pressure, gyroviscous; + + BoutReal vacuum_pressure; + BoutReal vacuum_trans; // Transition width + Field3D vac_mask; + + bool nonlinear; + bool evolve_jpar; + BoutReal g; // Only if compressible + bool phi_curv; + + // Poisson brackets: b0 x Grad(f) dot Grad(g) / B = [f, g] + // Method to use: BRACKET_ARAKAWA, BRACKET_STD or BRACKET_SIMPLE + /* + * Bracket method + * + * BRACKET_STD - Same as b0xGrad_dot_Grad, methods in BOUT.inp + * BRACKET_SIMPLE - Subset of terms, used in BOUT-06 + * BRACKET_ARAKAWA - Arakawa central differencing (2nd order) + * BRACKET_CTU - 1st order upwind method + * + */ + + // Bracket method for advection terms + BRACKET_METHOD bm_exb; + BRACKET_METHOD bm_mag; + int bm_exb_flag; + int bm_mag_flag; + /* BRACKET_METHOD bm_ExB = BRACKET_STD; + BRACKET_METHOD bm_mflutter = BRACKET_STD; */ + + bool diamag; + bool diamag_grad_t; // Grad_par(Te) term in Psi equation + bool diamag_phi0; // Include the diamagnetic equilibrium phi0 + + bool eHall; + BoutReal AA; // ion mass in units of the proton mass; AA=Mi/Mp + + // net flow, Er=-R*Bp*Dphi0,Dphi0=-D_min-0.5*D_0*(1.0-tanh(D_s*(x-x0))) + Field2D V0; // net flow amplitude + Field2D Dphi0; // differential potential to flux + BoutReal D_0; // potential amplitude + BoutReal D_s; // shear parameter + BoutReal x0; // velocity peak location + BoutReal sign; // direction of flow + BoutReal Psiaxis, Psibndry; + bool withflow; + bool K_H_term; // Kelvin-Holmhotz term + Field2D perp; // for test + BoutReal D_min; // constant in flow + + // for C_mod + bool experiment_Er; // read in total Er from experiment + + bool nogradparj; + bool filter_z; + int filter_z_mode; + int low_pass_z; + bool zonal_flow; + bool zonal_field; + bool zonal_bkgd; + + bool split_n0; // Solve the n=0 component of potential #if BOUT_HAS_HYPRE - std::unique_ptr laplacexy{nullptr}; // Laplacian solver in X-Y (n=0) + std::unique_ptr laplacexy{nullptr}; // Laplacian solver in X-Y (n=0) #else - std::unique_ptr laplacexy{nullptr}; // Laplacian solver in X-Y (n=0) + std::unique_ptr laplacexy{nullptr}; // Laplacian solver in X-Y (n=0) #endif - Field2D phi2D; // Axisymmetric phi - - bool relax_j_vac; - BoutReal relax_j_tconst; // Time-constant for j relax - - bool smooth_j_x; // Smooth Jpar in the x direction - - int jpar_bndry_width; // Zero jpar in a boundary region - - bool sheath_boundaries; // Apply sheath boundaries in Y - - bool parallel_lr_diff; // Use left and right shifted stencils for parallel differences - - bool phi_constraint; // Solver for phi using a solver constraint - - bool include_rmp; // Include RMP coil perturbation - bool simple_rmp; // Just use a simple form for the perturbation - - BoutReal rmp_factor; // Multiply amplitude by this factor - BoutReal rmp_ramp; // Ramp-up time for RMP [s]. negative -> instant - BoutReal rmp_freq; // Amplitude oscillation frequency [Hz] (negative -> no oscillation) - BoutReal rmp_rotate; // Rotation rate [Hz] - bool rmp_vac_mask; - Field3D rmp_Psi0; // Parallel vector potential from Resonant Magnetic Perturbation (RMP) - // coils - Field3D rmp_Psi; // Value used in calculations - Field3D rmp_dApdt; // Time variation - - BoutReal vac_lund, core_lund; // Lundquist number S = (Tau_R / Tau_A). -ve -> infty - BoutReal vac_resist, core_resist; // The resistivities (just 1 / S) - Field3D eta; // Resistivity profile (1 / S) - bool spitzer_resist; // Use Spitzer formula for resistivity - BoutReal Zeff; // Z effective for resistivity formula - - BoutReal hyperresist; // Hyper-resistivity coefficient (in core only) - BoutReal ehyperviscos; // electron Hyper-viscosity coefficient - - int damp_width; // Width of inner damped region - BoutReal damp_t_const; // Timescale of damping - - bout::TokamakOptions tokamak_options = bout::TokamakOptions(*mesh); - - const BoutReal MU0 = 4.0e-7 * PI; - const BoutReal Mi = 2.0 * 1.6726e-27; // Ion mass - const BoutReal Me = 9.1094e-31; // Electron mass - const BoutReal mi_me = Mi / Me; - - // Communication objects - FieldGroup comms; - - /// Solver for inverting Laplacian - std::unique_ptr phiSolver{nullptr}; - std::unique_ptr aparSolver{nullptr}; - - const Field2D N0tanh(BoutReal n0_height, BoutReal n0_ave, BoutReal n0_width, - BoutReal n0_center, BoutReal n0_bottom_x) { - Field2D result; - result.allocate(); - - BoutReal Grid_NX, Grid_NXlimit; // the grid number on x, and the - BoutReal Jysep; - mesh->get(Grid_NX, "nx"); - mesh->get(Jysep, "jyseps1_1"); - Grid_NXlimit = n0_bottom_x * Grid_NX; - output.write("Jysep1_1 = {:d} Grid number = {:e}\n", int(Jysep), Grid_NX); - - if (Jysep > 0.) { // for single null geometry - BoutReal Jxsep, Jysep2; - mesh->get(Jxsep, "ixseps1"); - mesh->get(Jysep2, "jyseps2_2"); - - for (auto i: result) { - BoutReal mgx = mesh->GlobalX(i.x()); - BoutReal xgrid_num = (Jxsep + 1.) / Grid_NX; - - int globaly = mesh->getGlobalYIndex(i.y()); - - if (mgx > xgrid_num || (globaly <= int(Jysep) - 2) - || (globaly > int(Jysep2) + 2)) { - mgx = xgrid_num; - } - BoutReal rlx = mgx - n0_center; - BoutReal temp = exp(rlx / n0_width); - BoutReal dampr = ((temp - 1.0 / temp) / (temp + 1.0 / temp)); - result[i] = 0.5 * (1.0 - dampr) * n0_height + n0_ave; - } - } else { // circular geometry - for (auto i: result) { - BoutReal mgx = mesh->GlobalX(i.x()); - BoutReal xgrid_num = Grid_NXlimit / Grid_NX; - if (mgx > xgrid_num) { - mgx = xgrid_num; - } - BoutReal rlx = mgx - n0_center; - BoutReal temp = exp(rlx / n0_width); - BoutReal dampr = ((temp - 1.0 / temp) / (temp + 1.0 / temp)); - result[i] = 0.5 * (1.0 - dampr) * n0_height + n0_ave; - } - } - - mesh->communicate(result); - - return result; + Field2D phi2D; // Axisymmetric phi + + bool relax_j_vac; + BoutReal relax_j_tconst; // Time-constant for j relax + + bool smooth_j_x; // Smooth Jpar in the x direction + + int jpar_bndry_width; // Zero jpar in a boundary region + + bool sheath_boundaries; // Apply sheath boundaries in Y + + bool parallel_lr_diff; // Use left and right shifted stencils for parallel differences + + bool phi_constraint; // Solver for phi using a solver constraint + bool phi_boundary_relax; // Relax x boundaries of phi towards Neumann? + bool phi_core_averagey; // Average phi core boundary in Y? + BoutReal phi_boundary_timescale; // Relaxation timescale + BoutReal phi_boundary_last_update; // Time when last updated + + bool include_rmp; // Include RMP coil perturbation + bool simple_rmp; // Just use a simple form for the perturbation + + BoutReal rmp_factor; // Multiply amplitude by this factor + BoutReal rmp_ramp; // Ramp-up time for RMP [s]. negative -> instant + BoutReal rmp_freq; // Amplitude oscillation frequency [Hz] (negative -> no oscillation) + BoutReal rmp_rotate; // Rotation rate [Hz] + bool rmp_vac_mask; + Field3D rmp_Psi0; // Parallel vector potential from Resonant Magnetic Perturbation (RMP) + // coils + Field3D rmp_Psi; // Value used in calculations + Field3D rmp_dApdt; // Time variation + + BoutReal vac_lund, core_lund; // Lundquist number S = (Tau_R / Tau_A). -ve -> infty + BoutReal vac_resist, core_resist; // The resistivities (just 1 / S) + Field3D eta; // Resistivity profile (1 / S) + bool spitzer_resist; // Use Spitzer formula for resistivity + BoutReal Zeff; // Z effective for resistivity formula + + BoutReal hyperresist; // Hyper-resistivity coefficient (in core only) + BoutReal ehyperviscos; // electron Hyper-viscosity coefficient + + int damp_width; // Width of inner damped region + BoutReal damp_t_const; // Timescale of damping + + // Metric coefficients + Field2D Rxy, Bpxy, Btxy, B0, hthe; + Field2D I; // Shear factor + + const BoutReal MU0 = 4.0e-7 * PI; + const BoutReal Mi = 2.0 * 1.6726e-27; // Ion mass + const BoutReal Me = 9.1094e-31; // Electron mass + const BoutReal mi_me = Mi / Me; + + // Communication objects + FieldGroup comms; + + /// Solver for inverting Laplacian + std::unique_ptr phiSolver{nullptr}; + std::unique_ptr aparSolver{nullptr}; + + const Field2D N0tanh(BoutReal n0_height, BoutReal n0_ave, BoutReal n0_width, + BoutReal n0_center, BoutReal n0_bottom_x) { + Field2D result; + result.allocate(); + + BoutReal Grid_NX, Grid_NXlimit; // the grid number on x, and the + BoutReal Jysep; + mesh->get(Grid_NX, "nx"); + mesh->get(Jysep, "jyseps1_1"); + Grid_NXlimit = n0_bottom_x * Grid_NX; + output.write("Jysep1_1 = {:d} Grid number = {:e}\n", int(Jysep), Grid_NX); + + if (Jysep > 0.) { // for single null geometry + BoutReal Jxsep, Jysep2; + mesh->get(Jxsep, "ixseps1"); + mesh->get(Jysep2, "jyseps2_2"); + + for (auto i : result) { + BoutReal mgx = mesh->GlobalX(i.x()); + BoutReal xgrid_num = (Jxsep + 1.) / Grid_NX; + + int globaly = mesh->getGlobalYIndex(i.y()); + + if (mgx > xgrid_num || (globaly <= int(Jysep) - 2) + || (globaly > int(Jysep2) + 2)) { + mgx = xgrid_num; + } + BoutReal rlx = mgx - n0_center; + BoutReal temp = exp(rlx / n0_width); + BoutReal dampr = ((temp - 1.0 / temp) / (temp + 1.0 / temp)); + result[i] = 0.5 * (1.0 - dampr) * n0_height + n0_ave; + } + } else { // circular geometry + for (auto i : result) { + BoutReal mgx = mesh->GlobalX(i.x()); + BoutReal xgrid_num = Grid_NXlimit / Grid_NX; + if (mgx > xgrid_num) { + mgx = xgrid_num; + } + BoutReal rlx = mgx - n0_center; + BoutReal temp = exp(rlx / n0_width); + BoutReal dampr = ((temp - 1.0 / temp) / (temp + 1.0 / temp)); + result[i] = 0.5 * (1.0 - dampr) * n0_height + n0_ave; + } } -protected: - int init(bool restarting) override { + mesh->communicate(result); - bool noshear; + return result; + } - Coordinates *metric = mesh->getCoordinates(); +protected: + int init(bool restarting) override { + bool noshear; - output.write("Solving high-beta flute reduced equations\n"); - output.write("\tFile : {:s}\n", __FILE__); - output.write("\tCompiled: {:s} at {:s}\n", __DATE__, __TIME__); + Coordinates* metric = mesh->getCoordinates(); - ////////////////////////////////////////////////////////////// - // Load data from the grid + output.write("Solving high-beta flute reduced equations\n"); + output.write("\tFile : {:s}\n", __FILE__); + output.write("\tCompiled: {:s} at {:s}\n", __DATE__, __TIME__); - // Load 2D profiles - mesh->get(J0, "Jpar0"); // A / m^2 - mesh->get(P0, "pressure"); // Pascals + ////////////////////////////////////////////////////////////// + // Load data from the grid - // Load curvature term - b0xcv.covariant = false; // Read contravariant components - mesh->get(b0xcv, "bxcv"); // mixed units x: T y: m^-2 z: m^-2 + // Load 2D profiles + mesh->get(J0, "Jpar0"); // A / m^2 + mesh->get(P0, "pressure"); // Pascals - mesh->get(Psixy, "psixy"); // get Psi - mesh->get(Psiaxis, "psi_axis"); // axis flux - mesh->get(Psibndry, "psi_bndry"); // edge flux + // Load curvature term + b0xcv.covariant = false; // Read contravariant components + mesh->get(b0xcv, "bxcv"); // mixed units x: T y: m^-2 z: m^-2 - // Set locations of staggered variables - // Note, use of staggered grids in elm-pb is untested and may not be completely - // implemented. Parallel boundary conditions are especially likely to be wrong. - if (mesh->StaggerGrids) { - loc = CELL_YLOW; - } else { - loc = CELL_CENTRE; - } - Jpar.setLocation(loc); - Vpar.setLocation(loc); - Psi.setLocation(loc); - eta.setLocation(loc); - Psi_loc.setBoundary("Psi_loc"); - - ////////////////////////////////////////////////////////////// - auto &globalOptions = Options::root(); - auto &options = globalOptions["highbeta"]; - - constn0 = options["constn0"].withDefault(true); - // use the hyperbolic profile of n0. If both n0_fake_prof and - // T0_fake_prof are false, use the profiles from grid file - n0_fake_prof = options["n0_fake_prof"].withDefault(false); - // the total height of profile of N0, in percentage of Ni_x - n0_height = options["n0_height"].withDefault(0.4); - // the center or average of N0, in percentage of Ni_x - n0_ave = options["n0_ave"].withDefault(0.01); - // the width of the gradient of N0,in percentage of x - n0_width = options["n0_width"].withDefault(0.1); - // the grid number of the center of N0, in percentage of x - n0_center = options["n0_center"].withDefault(0.633); - // the start of flat region of N0 on SOL side, in percentage of x - n0_bottom_x = options["n0_bottom_x"].withDefault(0.81); - T0_fake_prof = options["T0_fake_prof"].withDefault(false); - // the amplitude of constant temperature, in percentage - Tconst = options["Tconst"].withDefault(-1.0); - - density = options["density"].doc("Number density [m^-3]").withDefault(1.0e19); - - evolve_jpar = options["evolve_jpar"] - .doc("If true, evolve J rather than Psi") - .withDefault(false); - phi_constraint = options["phi_constraint"] - .doc("Use solver constraint for phi?") - .withDefault(false); - - // Effects to include/exclude - include_curvature = options["include_curvature"].withDefault(true); - include_jpar0 = options["include_jpar0"].withDefault(true); - evolve_pressure = options["evolve_pressure"].withDefault(true); - nogradparj = options["nogradparj"].withDefault(false); - - compress = options["compress"] - .doc("Include compressibility effects (evolve Vpar)?") - .withDefault(false); - gyroviscous = options["gyroviscous"].withDefault(false); - nonlinear = options["nonlinear"].doc("Include nonlinear terms?").withDefault(false); - - // option for ExB Poisson Bracket - bm_exb_flag = options["bm_exb_flag"] - .doc("ExB Poisson bracket method. 0=standard;1=simple;2=arakawa") - .withDefault(0); - switch (bm_exb_flag) { - case 0: { - bm_exb = BRACKET_STD; - output << "\tBrackets for ExB: default differencing\n"; - break; - } - case 1: { - bm_exb = BRACKET_SIMPLE; - output << "\tBrackets for ExB: simplified operator\n"; - break; - } - case 2: { - bm_exb = BRACKET_ARAKAWA; - output << "\tBrackets for ExB: Arakawa scheme\n"; - break; - } - case 3: { - bm_exb = BRACKET_CTU; - output << "\tBrackets for ExB: Corner Transport Upwind method\n"; - break; - } - default: - throw BoutException("Invalid choice of bracket method. Must be 0 - 3\n"); - } + // Load metrics + if (mesh->get(Rxy, "Rxy")) { // m + throw BoutException("Error: Cannot read Rxy from grid\n"); + } + if (mesh->get(Bpxy, "Bpxy")) { // T + throw BoutException("Error: Cannot read Bpxy from grid\n"); + } + mesh->get(Btxy, "Btxy"); // T + mesh->get(B0, "Bxy"); // T + mesh->get(hthe, "hthe"); // m + mesh->get(I, "sinty"); // m^-2 T^-1 + mesh->get(Psixy, "psixy"); // get Psi + mesh->get(Psiaxis, "psi_axis"); // axis flux + mesh->get(Psibndry, "psi_bndry"); // edge flux + + // Set locations of staggered variables + // Note, use of staggered grids in elm-pb is untested and may not be completely + // implemented. Parallel boundary conditions are especially likely to be wrong. + if (mesh->StaggerGrids) { + loc = CELL_YLOW; + } else { + loc = CELL_CENTRE; + } + Jpar.setLocation(loc); + Vpar.setLocation(loc); + Psi.setLocation(loc); + eta.setLocation(loc); + Psi_loc.setBoundary("Psi_loc"); + + ////////////////////////////////////////////////////////////// + auto& globalOptions = Options::root(); + auto& options = globalOptions["highbeta"]; + laplace_perp = options["laplace_perp"].withDefault(false); + // Use Laplace_perp rather than Delp2 + constn0 = options["constn0"].withDefault(true); + // use the hyperbolic profile of n0. If both n0_fake_prof and + // T0_fake_prof are false, use the profiles from grid file + n0_fake_prof = options["n0_fake_prof"].withDefault(false); + // the total height of profile of N0, in percentage of Ni_x + n0_height = options["n0_height"].withDefault(0.4); + // the center or average of N0, in percentage of Ni_x + n0_ave = options["n0_ave"].withDefault(0.01); + // the width of the gradient of N0,in percentage of x + n0_width = options["n0_width"].withDefault(0.1); + // the grid number of the center of N0, in percentage of x + n0_center = options["n0_center"].withDefault(0.633); + // the start of flat region of N0 on SOL side, in percentage of x + n0_bottom_x = options["n0_bottom_x"].withDefault(0.81); + T0_fake_prof = options["T0_fake_prof"].withDefault(false); + // the amplitude of constant temperature, in percentage + Tconst = options["Tconst"].withDefault(-1.0); + + density = options["density"].doc("Number density [m^-3]").withDefault(1.0e19); + + evolve_jpar = options["evolve_jpar"] + .doc("If true, evolve J rather than Psi") + .withDefault(false); + phi_constraint = options["phi_constraint"] + .doc("Use solver constraint for phi?") + .withDefault(false); + phi_boundary_relax = options["phi_boundary_relax"] + .doc("Relax x boundaries of phi towards Neumann?") + .withDefault(false); + + // Effects to include/exclude + include_curvature = options["include_curvature"].withDefault(true); + include_jpar0 = options["include_jpar0"].withDefault(true); + evolve_pressure = options["evolve_pressure"].withDefault(true); + nogradparj = options["nogradparj"].withDefault(false); + + compress = options["compress"] + .doc("Include compressibility effects (evolve Vpar)?") + .withDefault(false); + gyroviscous = options["gyroviscous"].withDefault(false); + nonlinear = options["nonlinear"].doc("Include nonlinear terms?").withDefault(false); + + // option for ExB Poisson Bracket + bm_exb_flag = options["bm_exb_flag"] + .doc("ExB Poisson bracket method. 0=standard;1=simple;2=arakawa") + .withDefault(0); + switch (bm_exb_flag) { + case 0: { + bm_exb = BRACKET_STD; + output << "\tBrackets for ExB: default differencing\n"; + break; + } + case 1: { + bm_exb = BRACKET_SIMPLE; + output << "\tBrackets for ExB: simplified operator\n"; + break; + } + case 2: { + bm_exb = BRACKET_ARAKAWA; + output << "\tBrackets for ExB: Arakawa scheme\n"; + break; + } + case 3: { + bm_exb = BRACKET_CTU; + output << "\tBrackets for ExB: Corner Transport Upwind method\n"; + break; + } + default: + throw BoutException("Invalid choice of bracket method. Must be 0 - 3\n"); + } - bm_mag_flag = - options["bm_mag_flag"].doc("magnetic flutter Poisson Bracket").withDefault(0); - switch (bm_mag_flag) { - case 0: { - bm_mag = BRACKET_STD; - output << "\tBrackets: default differencing\n"; - break; - } - case 1: { - bm_mag = BRACKET_SIMPLE; - output << "\tBrackets: simplified operator\n"; - break; - } - case 2: { - bm_mag = BRACKET_ARAKAWA; - output << "\tBrackets: Arakawa scheme\n"; - break; - } - case 3: { - bm_mag = BRACKET_CTU; - output << "\tBrackets: Corner Transport Upwind method\n"; - break; - } - default: - throw BoutException("Invalid choice of bracket method. Must be 0 - 3\n"); - } + bm_mag_flag = + options["bm_mag_flag"].doc("magnetic flutter Poisson Bracket").withDefault(0); + switch (bm_mag_flag) { + case 0: { + bm_mag = BRACKET_STD; + output << "\tBrackets: default differencing\n"; + break; + } + case 1: { + bm_mag = BRACKET_SIMPLE; + output << "\tBrackets: simplified operator\n"; + break; + } + case 2: { + bm_mag = BRACKET_ARAKAWA; + output << "\tBrackets: Arakawa scheme\n"; + break; + } + case 3: { + bm_mag = BRACKET_CTU; + output << "\tBrackets: Corner Transport Upwind method\n"; + break; + } + default: + throw BoutException("Invalid choice of bracket method. Must be 0 - 3\n"); + } - eHall = options["eHall"] + eHall = options["eHall"] .doc("electron Hall or electron parallel pressue gradient effects?") .withDefault(false); - AA = options["AA"].doc("ion mass in units of proton mass").withDefault(1.0); - - diamag = options["diamag"].doc("Diamagnetic effects?").withDefault(false); - diamag_grad_t = options["diamag_grad_t"] - .doc("Grad_par(Te) term in Psi equation") - .withDefault(diamag); - diamag_phi0 = - options["diamag_phi0"].doc("Include equilibrium phi0").withDefault(diamag); - dia_fact = options["dia_fact"] - .doc("Scale diamagnetic effects by this factor") - .withDefault(1.0); - - // withflow or not - withflow = options["withflow"].withDefault(false); - // keep K-H term - K_H_term = options["K_H_term"].withDefault(true); - // velocity magnitude - D_0 = options["D_0"].withDefault(0.0); - // flowshear - D_s = options["D_s"].withDefault(0.0); - // flow location - x0 = options["x0"].withDefault(0.0); - // flow direction, -1 means negative electric field - sign = options["sign"].withDefault(1.0); - // a constant - D_min = options["D_min"].withDefault(3000.0); - - experiment_Er = options["experiment_Er"].withDefault(false); - - noshear = options["noshear"].withDefault(false); - - relax_j_vac = options["relax_j_vac"] - .doc("Relax vacuum current to zero") - .withDefault(false); - relax_j_tconst = options["relax_j_tconst"] - .doc("Time constant for relaxation of vacuum current. Alfven " - "(normalised) units") - .withDefault(0.1); - - // Toroidal filtering - filter_z = options["filter_z"] - .doc("Filter a single toroidal mode number? The mode to keep is " - "filter_z_mode.") - .withDefault(false); - filter_z_mode = options["filter_z_mode"] - .doc("Single toroidal mode number to keep") - .withDefault(1); - low_pass_z = options["low_pass_z"].doc("Low-pass filter").withDefault(-1); - zonal_flow = options["zonal_flow"] - .doc("Keep zonal (n=0) component of potential?") - .withDefault(false); - zonal_field = options["zonal_field"] - .doc("Keep zonal (n=0) component of magnetic potential?") - .withDefault(false); - zonal_bkgd = options["zonal_bkgd"] - .doc("Evolve zonal (n=0) pressure profile?") - .withDefault(false); - - // n = 0 electrostatic potential solve - split_n0 = options["split_n0"] - .doc("Solve zonal (n=0) component of potential using LaplaceXY?") - .withDefault(false); - if (split_n0) { - // Create an XY solver for n=0 component + AA = options["AA"].doc("ion mass in units of proton mass").withDefault(1.0); + + diamag = options["diamag"].doc("Diamagnetic effects?").withDefault(false); + diamag_grad_t = options["diamag_grad_t"] + .doc("Grad_par(Te) term in Psi equation") + .withDefault(diamag); + diamag_phi0 = + options["diamag_phi0"].doc("Include equilibrium phi0").withDefault(diamag); + dia_fact = options["dia_fact"] + .doc("Scale diamagnetic effects by this factor") + .withDefault(1.0); + + // withflow or not + withflow = options["withflow"].withDefault(false); + // keep K-H term + K_H_term = options["K_H_term"].withDefault(true); + // velocity magnitude + D_0 = options["D_0"].withDefault(0.0); + // flowshear + D_s = options["D_s"].withDefault(0.0); + // flow location + x0 = options["x0"].withDefault(0.0); + // flow direction, -1 means negative electric field + sign = options["sign"].withDefault(1.0); + // a constant + D_min = options["D_min"].withDefault(3000.0); + + experiment_Er = options["experiment_Er"].withDefault(false); + + noshear = options["noshear"].withDefault(false); + + relax_j_vac = options["relax_j_vac"] + .doc("Relax vacuum current to zero") + .withDefault(false); + relax_j_tconst = options["relax_j_tconst"] + .doc("Time constant for relaxation of vacuum current. Alfven " + "(normalised) units") + .withDefault(0.1); + + // Toroidal filtering + filter_z = options["filter_z"] + .doc("Filter a single toroidal mode number? The mode to keep is " + "filter_z_mode.") + .withDefault(false); + filter_z_mode = options["filter_z_mode"] + .doc("Single toroidal mode number to keep") + .withDefault(1); + low_pass_z = options["low_pass_z"].doc("Low-pass filter").withDefault(-1); + zonal_flow = options["zonal_flow"] + .doc("Keep zonal (n=0) component of potential?") + .withDefault(false); + zonal_field = options["zonal_field"] + .doc("Keep zonal (n=0) component of magnetic potential?") + .withDefault(false); + zonal_bkgd = options["zonal_bkgd"] + .doc("Evolve zonal (n=0) pressure profile?") + .withDefault(false); + + // n = 0 electrostatic potential solve + split_n0 = options["split_n0"] + .doc("Solve zonal (n=0) component of potential using LaplaceXY?") + .withDefault(false); + if (split_n0) { + // Create an XY solver for n=0 component #if BOUT_HAS_HYPRE - laplacexy = bout::utils::make_unique(mesh); + laplacexy = bout::utils::make_unique(mesh); #else - laplacexy = bout::utils::make_unique(mesh); + laplacexy = bout::utils::make_unique(mesh); #endif - // Set coefficients for Boussinesq solve - laplacexy->setCoefs(1.0, 0.0); - phi2D = 0.0; // Starting guess - phi2D.setBoundary("phi"); - } - - // Radial smoothing - smooth_j_x = options["smooth_j_x"].doc("Smooth Jpar in x").withDefault(false); - - // Jpar boundary region - jpar_bndry_width = options["jpar_bndry_width"] - .doc("Number of cells near the boundary where jpar = 0") - .withDefault(-1); - - sheath_boundaries = options["sheath_boundaries"] - .doc("Apply sheath boundaries in Y?") - .withDefault(false); + // Set coefficients for Boussinesq solve + laplacexy->setCoefs(1.0, 0.0); + phi2D = 0.0; // Starting guess + phi2D.setBoundary("phi"); + } - // Parallel differencing - parallel_lr_diff = - options["parallel_lr_diff"] - .doc("Use left and right shifted stencils for parallel differences?") - .withDefault(false); - - // RMP-related options - include_rmp = options["include_rmp"] - .doc("Read RMP field rmp_A from grid?") - .withDefault(false); - - simple_rmp = - options["simple_rmp"].doc("Include a simple RMP model?").withDefault(false); - rmp_factor = options["rmp_factor"].withDefault(1.0); - rmp_ramp = options["rmp_ramp"].withDefault(-1.0); - rmp_freq = options["rmp_freq"].withDefault(-1.0); - rmp_rotate = options["rmp_rotate"].withDefault(0.0); - - // Vacuum region control - vacuum_pressure = - options["vacuum_pressure"] - .doc("Fraction of peak pressure, below which is considered vacuum.") - .withDefault(0.02); - vacuum_trans = - options["vacuum_trans"] - .doc("Vacuum boundary transition width, as fraction of peak pressure.") - .withDefault(0.005); - - // Resistivity and hyper-resistivity options - vac_lund = - options["vac_lund"].doc("Lundquist number in vacuum region").withDefault(0.0); - core_lund = - options["core_lund"].doc("Lundquist number in core region").withDefault(0.0); - hyperresist = options["hyperresist"].withDefault(-1.0); - ehyperviscos = options["ehyperviscos"].withDefault(-1.0); - spitzer_resist = options["spitzer_resist"] - .doc("Use Spitzer resistivity?") - .withDefault(false); - Zeff = options["Zeff"].withDefault(2.0); // Z effective - - // Inner boundary damping - damp_width = options["damp_width"] - .doc("Width of the radial damping regions, in grid cells") - .withDefault(0); - damp_t_const = - options["damp_t_const"] - .doc("Time constant for damping in radial regions. Normalised time units.") - .withDefault(0.1); - - // Viscosity and hyper-viscosity - viscos_par = options["viscos_par"].doc("Parallel viscosity").withDefault(-1.0); - viscos_perp = options["viscos_perp"].doc("Perpendicular viscosity").withDefault(-1.0); - hyperviscos = options["hyperviscos"].doc("Radial hyperviscosity").withDefault(-1.0); - - diffusion_par = - options["diffusion_par"].doc("Parallel pressure diffusion").withDefault(-1.0); - diffusion_p4 = options["diffusion_p4"] - .doc("parallel hyper-viscous diffusion for pressure") - .withDefault(-1.0); - diffusion_u4 = options["diffusion_u4"] - .doc("parallel hyper-viscous diffusion for vorticity") - .withDefault(-1.0); - diffusion_a4 = options["diffusion_a4"] - .doc("parallel hyper-viscous diffusion for vector potential") - .withDefault(-1.0); - - // heating factor in pressure - // heating power in pressure - heating_P = options["heating_P"].withDefault(-1.0); - // the percentage of radial grid points for heating profile radial - // width in pressure - hp_width = options["hp_width"].withDefault(0.1); - // the percentage of radial grid points for heating profile radial - // domain in pressure - hp_length = options["hp_length"].withDefault(0.04); - - // sink factor in pressure - // sink in pressure - sink_P = options["sink_P"].withDefault(-1.0); - // the percentage of radial grid points for sink profile radial - // width in pressure - sp_width = options["sp_width"].withDefault(0.05); - // the percentage of radial grid points for sink profile radial - // domain in pressure - sp_length = options["sp_length"].withDefault(0.04); - - // left edge sink factor in vorticity - // left edge sink in vorticity - sink_Ul = options["sink_Ul"].withDefault(-1.0); - // the percentage of left edge radial grid points for sink profile - // radial width in vorticity - su_widthl = options["su_widthl"].withDefault(0.06); - // the percentage of left edge radial grid points for sink profile - // radial domain in vorticity - su_lengthl = options["su_lengthl"].withDefault(0.15); - - // right edge sink factor in vorticity - // right edge sink in vorticity - sink_Ur = options["sink_Ur"].withDefault(-1.0); - // the percentage of right edge radial grid points for sink profile - // radial width in vorticity - su_widthr = options["su_widthr"].withDefault(0.06); - // the percentage of right edge radial grid points for sink profile - // radial domain in vorticity - su_lengthr = options["su_lengthr"].withDefault(0.15); - - // Compressional terms - phi_curv = - options["phi_curv"].doc("ExB compression in P equation?").withDefault(true); - g = options["gamma"].doc("Ratio of specific heats").withDefault(5.0 / 3.0); - - x = (Psixy - Psiaxis) / (Psibndry - Psiaxis); - - if (experiment_Er) { // get er from experiment - mesh->get(Dphi0, "Epsi"); - diamag_phi0 = false; - K_H_term = false; - } else { - Dphi0 = -D_min - 0.5 * D_0 * (1.0 - tanh(D_s * (x - x0))); - } + // Radial smoothing + smooth_j_x = options["smooth_j_x"].doc("Smooth Jpar in x").withDefault(false); + + // Jpar boundary region + jpar_bndry_width = options["jpar_bndry_width"] + .doc("Number of cells near the boundary where jpar = 0") + .withDefault(-1); + + sheath_boundaries = options["sheath_boundaries"] + .doc("Apply sheath boundaries in Y?") + .withDefault(false); + + // Parallel differencing + parallel_lr_diff = + options["parallel_lr_diff"] + .doc("Use left and right shifted stencils for parallel differences?") + .withDefault(false); + + // RMP-related options + include_rmp = options["include_rmp"] + .doc("Read RMP field rmp_A from grid?") + .withDefault(false); + + simple_rmp = + options["simple_rmp"].doc("Include a simple RMP model?").withDefault(false); + rmp_factor = options["rmp_factor"].withDefault(1.0); + rmp_ramp = options["rmp_ramp"].withDefault(-1.0); + rmp_freq = options["rmp_freq"].withDefault(-1.0); + rmp_rotate = options["rmp_rotate"].withDefault(0.0); + + // Vacuum region control + vacuum_pressure = + options["vacuum_pressure"] + .doc("Fraction of peak pressure, below which is considered vacuum.") + .withDefault(0.02); + vacuum_trans = + options["vacuum_trans"] + .doc("Vacuum boundary transition width, as fraction of peak pressure.") + .withDefault(0.005); + + // Resistivity and hyper-resistivity options + vac_lund = + options["vac_lund"].doc("Lundquist number in vacuum region").withDefault(0.0); + core_lund = + options["core_lund"].doc("Lundquist number in core region").withDefault(0.0); + hyperresist = options["hyperresist"].withDefault(-1.0); + ehyperviscos = options["ehyperviscos"].withDefault(-1.0); + spitzer_resist = options["spitzer_resist"] + .doc("Use Spitzer resistivity?") + .withDefault(false); + Zeff = options["Zeff"].withDefault(2.0); // Z effective + + // Inner boundary damping + damp_width = options["damp_width"] + .doc("Width of the radial damping regions, in grid cells") + .withDefault(0); + damp_t_const = + options["damp_t_const"] + .doc("Time constant for damping in radial regions. Normalised time units.") + .withDefault(0.1); + + // Viscosity and hyper-viscosity + viscos_par = options["viscos_par"].doc("Parallel viscosity").withDefault(-1.0); + viscos_perp = options["viscos_perp"].doc("Perpendicular viscosity").withDefault(-1.0); + hyperviscos = options["hyperviscos"].doc("Radial hyperviscosity").withDefault(-1.0); + + diffusion_par = + options["diffusion_par"].doc("Parallel pressure diffusion").withDefault(-1.0); + diffusion_p4 = options["diffusion_p4"] + .doc("parallel hyper-viscous diffusion for pressure") + .withDefault(-1.0); + diffusion_u4 = options["diffusion_u4"] + .doc("parallel hyper-viscous diffusion for vorticity") + .withDefault(-1.0); + diffusion_a4 = options["diffusion_a4"] + .doc("parallel hyper-viscous diffusion for vector potential") + .withDefault(-1.0); + + // heating factor in pressure + // heating power in pressure + heating_P = options["heating_P"].withDefault(-1.0); + // the percentage of radial grid points for heating profile radial + // width in pressure + hp_width = options["hp_width"].withDefault(0.1); + // the percentage of radial grid points for heating profile radial + // domain in pressure + hp_length = options["hp_length"].withDefault(0.04); + + // sink factor in pressure + // sink in pressure + sink_P = options["sink_P"].withDefault(-1.0); + // the percentage of radial grid points for sink profile radial + // width in pressure + sp_width = options["sp_width"].withDefault(0.05); + // the percentage of radial grid points for sink profile radial + // domain in pressure + sp_length = options["sp_length"].withDefault(0.04); + + // left edge sink factor in vorticity + // left edge sink in vorticity + sink_Ul = options["sink_Ul"].withDefault(-1.0); + // the percentage of left edge radial grid points for sink profile + // radial width in vorticity + su_widthl = options["su_widthl"].withDefault(0.06); + // the percentage of left edge radial grid points for sink profile + // radial domain in vorticity + su_lengthl = options["su_lengthl"].withDefault(0.15); + + // right edge sink factor in vorticity + // right edge sink in vorticity + sink_Ur = options["sink_Ur"].withDefault(-1.0); + // the percentage of right edge radial grid points for sink profile + // radial width in vorticity + su_widthr = options["su_widthr"].withDefault(0.06); + // the percentage of right edge radial grid points for sink profile + // radial domain in vorticity + su_lengthr = options["su_lengthr"].withDefault(0.15); + + // Compressional terms + phi_curv = + options["phi_curv"].doc("ExB compression in P equation?").withDefault(true); + g = options["gamma"].doc("Ratio of specific heats").withDefault(5.0 / 3.0); + + x = (Psixy - Psiaxis) / (Psibndry - Psiaxis); + + if (experiment_Er) { // get er from experiment + mesh->get(Dphi0, "Epsi"); + diamag_phi0 = false; + K_H_term = false; + } else { + Dphi0 = -D_min - 0.5 * D_0 * (1.0 - tanh(D_s * (x - x0))); + } - if (sign < 0) { // change flow direction - Dphi0 *= -1; - } + if (sign < 0) { // change flow direction + Dphi0 *= -1; + } - auto Bpxy = tokamak_options.Bpxy; - auto hthe = tokamak_options.hthe; - auto Rxy = tokamak_options.Rxy; - auto Btxy = tokamak_options.Btxy; - auto B0 = tokamak_options.Bxy; + V0 = -Rxy * Bpxy * Dphi0 / B0; - V0 = -Rxy * Bpxy * Dphi0 / B0; + if (simple_rmp) { + include_rmp = true; + } - if (simple_rmp) { - include_rmp = true; - } + // toroidal and poloidal mode numbers + const int rmp_n = + options["rmp_n"].doc("Simple RMP toroidal mode number").withDefault(3); + const int rmp_m = + options["rmp_m"].doc("Simple RMP poloidal mode number").withDefault(9); + const int rmp_polwid = options["rmp_polwid"] + .doc("Poloidal width (-ve -> full, fraction of 2pi)") + .withDefault(-1.0); + const int rmp_polpeak = options["rmp_polpeak"] + .doc("Peak poloidal location (fraction of 2pi)") + .withDefault(0.5); + rmp_vac_mask = + options["rmp_vac_mask"].doc("Should a vacuum mask be applied?").withDefault(true); + // Divide n by the size of the domain + const int zperiod = globalOptions["zperiod"].withDefault(1); + + if (include_rmp) { + // Including external field coils. + if (simple_rmp) { + // Use a fairly simple form for the perturbation + Field2D pol_angle; + if (mesh->get(pol_angle, "pol_angle")) { + output_warn.write(" ***WARNING: need poloidal angle for simple RMP\n"); + include_rmp = false; + } else { + if ((rmp_n % zperiod) != 0) { + output_warn.write( + " ***WARNING: rmp_n ({:d}) not a multiple of zperiod ({:d})\n", rmp_n, + zperiod); + } + + output.write("\tMagnetic perturbation: n = {:d}, m = {:d}, magnitude {:e} Tm\n", + rmp_n, rmp_m, rmp_factor); + + rmp_Psi0 = 0.0; + if (mesh->lastX()) { + // Set the outer boundary + for (int jx = mesh->LocalNx - 4; jx < mesh->LocalNx; jx++) { + for (int jy = 0; jy < mesh->LocalNy; jy++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { - // toroidal and poloidal mode numbers - const int rmp_n = - options["rmp_n"].doc("Simple RMP toroidal mode number").withDefault(3); - const int rmp_m = - options["rmp_m"].doc("Simple RMP poloidal mode number").withDefault(9); - const int rmp_polwid = options["rmp_polwid"] - .doc("Poloidal width (-ve -> full, fraction of 2pi)") - .withDefault(-1.0); - const int rmp_polpeak = options["rmp_polpeak"] - .doc("Peak poloidal location (fraction of 2pi)") - .withDefault(0.5); - rmp_vac_mask = - options["rmp_vac_mask"].doc("Should a vacuum mask be applied?").withDefault(true); - // Divide n by the size of the domain - const int zperiod = globalOptions["zperiod"].withDefault(1); - - if (include_rmp) { - // Including external field coils. - if (simple_rmp) { - // Use a fairly simple form for the perturbation - Field2D pol_angle; - if (mesh->get(pol_angle, "pol_angle")) { - output_warn.write(" ***WARNING: need poloidal angle for simple RMP\n"); - include_rmp = false; - } else { - if ((rmp_n % zperiod) != 0) { - output_warn.write( - " ***WARNING: rmp_n ({:d}) not a multiple of zperiod ({:d})\n", rmp_n, - zperiod); - } - - output.write("\tMagnetic perturbation: n = {:d}, m = {:d}, magnitude {:e} Tm\n", - rmp_n, rmp_m, rmp_factor); - - rmp_Psi0 = 0.0; - if (mesh->lastX()) { - // Set the outer boundary - for (int jx = mesh->LocalNx - 4; jx < mesh->LocalNx; jx++) { - for (int jy = 0; jy < mesh->LocalNy; jy++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { - - BoutReal angle = - rmp_m * pol_angle(jx, jy) - + rmp_n * ((BoutReal) jz) * mesh->getCoordinates()->dz(jx, jy, jz); - rmp_Psi0(jx, jy, jz) = - (((BoutReal) (jx - 4)) / ((BoutReal) (mesh->LocalNx - 5))) - * rmp_factor * cos(angle); - if (rmp_polwid > 0.0) { - // Multiply by a Gaussian in poloidal angle - BoutReal gx = - ((pol_angle(jx, jy) / (2. * PI)) - rmp_polpeak) / rmp_polwid; - rmp_Psi0(jx, jy, jz) *= exp(-gx * gx); - } - } - } - } - } - - // Now have a simple model for Psi due to coils at the outer boundary - // Need to calculate Psi inside the domain, enforcing j = 0 - - Jpar = 0.0; - auto psiLap = std::unique_ptr{Laplacian::create(nullptr, loc)}; - psiLap->setInnerBoundaryFlags(INVERT_AC_GRAD); // Zero gradient inner BC - psiLap->setOuterBoundaryFlags(INVERT_SET); // Set to rmp_Psi0 on outer boundary - rmp_Psi0 = psiLap->solve(Jpar, rmp_Psi0); - mesh->communicate(rmp_Psi0); + BoutReal angle = + rmp_m * pol_angle(jx, jy) + + rmp_n * ((BoutReal)jz) * mesh->getCoordinates()->dz(jx, jy, jz); + rmp_Psi0(jx, jy, jz) = + (((BoutReal)(jx - 4)) / ((BoutReal)(mesh->LocalNx - 5))) + * rmp_factor * cos(angle); + if (rmp_polwid > 0.0) { + // Multiply by a Gaussian in poloidal angle + BoutReal gx = + ((pol_angle(jx, jy) / (2. * PI)) - rmp_polpeak) / rmp_polwid; + rmp_Psi0(jx, jy, jz) *= exp(-gx * gx); + } } - } else { - // Load perturbation from grid file. - include_rmp = !mesh->get(rmp_Psi0, "rmp_A"); // Only include if found - if (!include_rmp) { - output_warn.write("WARNING: Couldn't read 'rmp_A' from grid file\n"); - } - // Multiply by factor - rmp_Psi0 *= rmp_factor; + } } - } + } + + // Now have a simple model for Psi due to coils at the outer boundary + // Need to calculate Psi inside the domain, enforcing j = 0 + + Jpar = 0.0; + auto psiLap = std::unique_ptr{Laplacian::create(nullptr, loc)}; + psiLap->setInnerBoundaryFlags(INVERT_AC_GRAD); // Zero gradient inner BC + psiLap->setOuterBoundaryFlags(INVERT_SET); // Set to rmp_Psi0 on outer boundary + rmp_Psi0 = psiLap->solve(Jpar, rmp_Psi0); + mesh->communicate(rmp_Psi0); + } + } else { + // Load perturbation from grid file. + include_rmp = !mesh->get(rmp_Psi0, "rmp_A"); // Only include if found + if (!include_rmp) { + output_warn.write("WARNING: Couldn't read 'rmp_A' from grid file\n"); + } + // Multiply by factor + rmp_Psi0 *= rmp_factor; + } + } - if (!include_curvature) { - b0xcv = 0.0; - } + if (!include_curvature) { + b0xcv = 0.0; + } - if (!include_jpar0) { - J0 = 0.0; - } + if (!include_jpar0) { + J0 = 0.0; + } - ////////////////////////////////////////////////////////////// - // SHIFTED RADIAL COORDINATES + if (noshear) { + if (include_curvature) { + b0xcv.z += I * b0xcv.x; + } + I = 0.0; + } - BoutReal shearFactor = 1.0; - if (!noshear && mesh->IncIntShear) { - // BOUT-06 style, using d/dx = d/dpsi + I * d/dz - metric->setIntShiftTorsion(tokamak_options.I); + ////////////////////////////////////////////////////////////// + // SHIFTED RADIAL COORDINATES - } else { - // Dimits style, using local coordinate system - if (include_curvature) { - b0xcv.z += tokamak_options.I * b0xcv.x; - } - shearFactor = 0.0; // I disappears from metric - } + if (mesh->IncIntShear) { + // BOUT-06 style, using d/dx = d/dpsi + I * d/dz + metric->IntShiftTorsion = I; - set_tokamak_coordinates_on_mesh(tokamak_options, *mesh, Lbar, Bbar, shearFactor); + } else { + // Dimits style, using local coordinate system + if (include_curvature) { + b0xcv.z += I * b0xcv.x; + } + I = 0.0; // I disappears from metric + } - ////////////////////////////////////////////////////////////// - // NORMALISE QUANTITIES + ////////////////////////////////////////////////////////////// + // NORMALISE QUANTITIES - if (mesh->get(Bbar, "bmag")) { // Typical magnetic field - Bbar = 1.0; - } - if (mesh->get(Lbar, "rmag")) { // Typical length scale - Lbar = 1.0; - } + if (mesh->get(Bbar, "bmag")) { // Typical magnetic field + Bbar = 1.0; + } + if (mesh->get(Lbar, "rmag")) { // Typical length scale + Lbar = 1.0; + } - Va = sqrt(Bbar * Bbar / (MU0 * density * Mi)); + Va = sqrt(Bbar * Bbar / (MU0 * density * Mi)); - Tbar = Lbar / Va; + Tbar = Lbar / Va; - dnorm = dia_fact * Mi / (2. * 1.602e-19 * Bbar * Tbar); + dnorm = dia_fact * Mi / (2. * 1.602e-19 * Bbar * Tbar); - delta_i = AA * 60.67 * 5.31e5 / sqrt(density / 1e6) / (Lbar * 100.0); + delta_i = AA * 60.67 * 5.31e5 / sqrt(density / 1e6) / (Lbar * 100.0); - output.write("Normalisations: Bbar = {:e} T Lbar = {:e} m\n", Bbar, Lbar); - output.write(" Va = {:e} m/s Tbar = {:e} s\n", Va, Tbar); - output.write(" dnorm = {:e}\n", dnorm); - output.write(" Resistivity\n"); + output.write("Normalisations: Bbar = {:e} T Lbar = {:e} m\n", Bbar, Lbar); + output.write(" Va = {:e} m/s Tbar = {:e} s\n", Va, Tbar); + output.write(" dnorm = {:e}\n", dnorm); + output.write(" Resistivity\n"); - if (gyroviscous) { - omega_i = 9.58e7 * Zeff * Bbar; - Upara2 = 0.5 / (Tbar * omega_i); - output.write("Upara2 = {:e} Omega_i = {:e}\n", Upara2, omega_i); - } + if (gyroviscous) { + omega_i = 9.58e7 * Zeff * Bbar; + Upara2 = 0.5 / (Tbar * omega_i); + output.write("Upara2 = {:e} Omega_i = {:e}\n", Upara2, omega_i); + } - if (eHall) { - output.write(" delta_i = {:e} AA = {:e} \n", delta_i, AA); - } + if (eHall) { + output.write(" delta_i = {:e} AA = {:e} \n", delta_i, AA); + } - if (vac_lund > 0.0) { - output.write(" Vacuum Tau_R = {:e} s eta = {:e} Ohm m\n", vac_lund * Tbar, - MU0 * Lbar * Lbar / (vac_lund * Tbar)); - vac_resist = 1. / vac_lund; - } else { - output.write(" Vacuum - Zero resistivity -\n"); - vac_resist = 0.0; - } - if (core_lund > 0.0) { - output.write(" Core Tau_R = {:e} s eta = {:e} Ohm m\n", - core_lund * Tbar, MU0 * Lbar * Lbar / (core_lund * Tbar)); - core_resist = 1. / core_lund; - } else { - output.write(" Core - Zero resistivity -\n"); - core_resist = 0.0; - } + if (vac_lund > 0.0) { + output.write(" Vacuum Tau_R = {:e} s eta = {:e} Ohm m\n", vac_lund * Tbar, + MU0 * Lbar * Lbar / (vac_lund * Tbar)); + vac_resist = 1. / vac_lund; + } else { + output.write(" Vacuum - Zero resistivity -\n"); + vac_resist = 0.0; + } + if (core_lund > 0.0) { + output.write(" Core Tau_R = {:e} s eta = {:e} Ohm m\n", + core_lund * Tbar, MU0 * Lbar * Lbar / (core_lund * Tbar)); + core_resist = 1. / core_lund; + } else { + output.write(" Core - Zero resistivity -\n"); + core_resist = 0.0; + } - if (hyperresist > 0.0) { - output.write(" Hyper-resistivity coefficient: {:e}\n", hyperresist); - } + if (hyperresist > 0.0) { + output.write(" Hyper-resistivity coefficient: {:e}\n", hyperresist); + } - if (ehyperviscos > 0.0) { - output.write(" electron Hyper-viscosity coefficient: {:e}\n", ehyperviscos); - } + if (ehyperviscos > 0.0) { + output.write(" electron Hyper-viscosity coefficient: {:e}\n", ehyperviscos); + } - if (hyperviscos > 0.0) { - output.write(" Hyper-viscosity coefficient: {:e}\n", hyperviscos); - SAVE_ONCE(hyper_mu_x); - } + if (hyperviscos > 0.0) { + output.write(" Hyper-viscosity coefficient: {:e}\n", hyperviscos); + SAVE_ONCE(hyper_mu_x); + } - if (diffusion_par > 0.0) { - output.write(" diffusion_par: {:e}\n", diffusion_par); - SAVE_ONCE(diffusion_par); - } + if (diffusion_par > 0.0) { + output.write(" diffusion_par: {:e}\n", diffusion_par); + SAVE_ONCE(diffusion_par); + } - // xqx: parallel hyper-viscous diffusion for pressure - if (diffusion_p4 > 0.0) { - output.write(" diffusion_p4: {:e}\n", diffusion_p4); - SAVE_ONCE(diffusion_p4); - } + // xqx: parallel hyper-viscous diffusion for pressure + if (diffusion_p4 > 0.0) { + output.write(" diffusion_p4: {:e}\n", diffusion_p4); + SAVE_ONCE(diffusion_p4); + } - // xqx: parallel hyper-viscous diffusion for vorticity - if (diffusion_u4 > 0.0) { - output.write(" diffusion_u4: {:e}\n", diffusion_u4); - SAVE_ONCE(diffusion_u4) - } + // xqx: parallel hyper-viscous diffusion for vorticity + if (diffusion_u4 > 0.0) { + output.write(" diffusion_u4: {:e}\n", diffusion_u4); + SAVE_ONCE(diffusion_u4) + } - // xqx: parallel hyper-viscous diffusion for vector potential - if (diffusion_a4 > 0.0) { - output.write(" diffusion_a4: {:e}\n", diffusion_a4); - SAVE_ONCE(diffusion_a4); - } + // xqx: parallel hyper-viscous diffusion for vector potential + if (diffusion_a4 > 0.0) { + output.write(" diffusion_a4: {:e}\n", diffusion_a4); + SAVE_ONCE(diffusion_a4); + } - if (heating_P > 0.0) { - output.write(" heating_P(watts): {:e}\n", heating_P); + if (heating_P > 0.0) { + output.write(" heating_P(watts): {:e}\n", heating_P); - output.write(" hp_width(%%): {:e}\n", hp_width); + output.write(" hp_width(%%): {:e}\n", hp_width); - output.write(" hp_length(%%): {:e}\n", hp_length); + output.write(" hp_length(%%): {:e}\n", hp_length); - SAVE_ONCE(heating_P, hp_width, hp_length); - } + SAVE_ONCE(heating_P, hp_width, hp_length); + } - if (sink_P > 0.0) { - output.write(" sink_P(rate): {:e}\n", sink_P); - output.write(" sp_width(%%): {:e}\n", sp_width); - output.write(" sp_length(%%): {:e}\n", sp_length); + if (sink_P > 0.0) { + output.write(" sink_P(rate): {:e}\n", sink_P); + output.write(" sp_width(%%): {:e}\n", sp_width); + output.write(" sp_length(%%): {:e}\n", sp_length); - SAVE_ONCE(sink_P, sp_width, sp_length); - } + SAVE_ONCE(sink_P, sp_width, sp_length); + } - if (K_H_term) { - output.write(" keep K-H term\n"); - } else { - output.write(" drop K-H term\n"); - } + if (K_H_term) { + output.write(" keep K-H term\n"); + } else { + output.write(" drop K-H term\n"); + } - Field2D Te; - Te = P0 / (2.0 * density * 1.602e-19); // Temperature in eV + Field2D Te; + Te = P0 / (2.0 * density * 1.602e-19); // Temperature in eV + + J0 = -MU0 * Lbar * J0 / B0; + P0 = 2.0 * MU0 * P0 / (Bbar * Bbar); + V0 = V0 / Va; + Dphi0 *= Tbar; + + b0xcv.x /= Bbar; + b0xcv.y *= Lbar * Lbar; + b0xcv.z *= Lbar * Lbar; + + Rxy /= Lbar; + Bpxy /= Bbar; + Btxy /= Bbar; + B0 /= Bbar; + hthe /= Lbar; + metric->dx /= Lbar * Lbar * Bbar; + I *= Lbar * Lbar * Bbar; + + if (constn0) { + T0_fake_prof = false; + n0_fake_prof = false; + } else { + Nbar = 1.0; + Tibar = 1000.0; + Tebar = 1000.0; + + if ((!T0_fake_prof) && n0_fake_prof) { + N0 = N0tanh(n0_height * Nbar, n0_ave * Nbar, n0_width, n0_center, n0_bottom_x); + + Ti0 = P0 / N0 / 2.0; + Te0 = Ti0; + } else if (T0_fake_prof) { + Ti0 = Tconst; + Te0 = Ti0; + N0 = P0 / (Ti0 + Te0); + } else { + if (mesh->get(N0, "Niexp")) { // N_i0 + throw BoutException("Error: Cannot read Ni0 from grid\n"); + } + + if (mesh->get(Ti0, "Tiexp")) { // T_i0 + throw BoutException("Error: Cannot read Ti0 from grid\n"); + } + + if (mesh->get(Te0, "Teexp")) { // T_e0 + throw BoutException("Error: Cannot read Te0 from grid\n"); + } + N0 /= Nbar; + Ti0 /= Tibar; + Te0 /= Tebar; + } + } - J0 = -MU0 * Lbar * J0 / B0; - P0 = 2.0 * MU0 * P0 / (Bbar * Bbar); - V0 = V0 / Va; - Dphi0 *= Tbar; + if (gyroviscous) { + Dperp2Phi0.setLocation(CELL_CENTRE); + Dperp2Phi0.setBoundary("phi"); + Dperp2Phi.setLocation(CELL_CENTRE); + Dperp2Phi.setBoundary("phi"); + GradPhi02.setLocation(CELL_CENTRE); + GradPhi02.setBoundary("phi"); + GradcPhi.setLocation(CELL_CENTRE); + GradcPhi.setBoundary("phi"); + Dperp2Pi0.setLocation(CELL_CENTRE); + Dperp2Pi0.setBoundary("P"); + Dperp2Pi.setLocation(CELL_CENTRE); + Dperp2Pi.setBoundary("P"); + bracketPhi0P.setLocation(CELL_CENTRE); + bracketPhi0P.setBoundary("P"); + bracketPhiP0.setLocation(CELL_CENTRE); + bracketPhiP0.setBoundary("P"); + if (nonlinear) { + GradPhi2.setLocation(CELL_CENTRE); + GradPhi2.setBoundary("phi"); + bracketPhiP.setLocation(CELL_CENTRE); + bracketPhiP.setBoundary("P"); + } + } - b0xcv.x /= Bbar; - b0xcv.y *= Lbar * Lbar; - b0xcv.z *= Lbar * Lbar; + BoutReal pnorm = max(P0, true); // Maximum over all processors + + vacuum_pressure *= pnorm; // Get pressure from fraction + vacuum_trans *= pnorm; + + // Transitions from 0 in core to 1 in vacuum + vac_mask = (1.0 - tanh((P0 - vacuum_pressure) / vacuum_trans)) / 2.0; + + if (spitzer_resist) { + // Use Spitzer resistivity + output.write("\tTemperature: {:e} -> {:e} [eV]\n", min(Te), max(Te)); + eta = 0.51 * 1.03e-4 * Zeff * 20. + * pow(Te, -1.5); // eta in Ohm-m. NOTE: ln(Lambda) = 20 + output.write("\tSpitzer resistivity: {:e} -> {:e} [Ohm m]\n", min(eta), max(eta)); + eta /= MU0 * Va * Lbar; + output.write("\t -> Lundquist {:e} -> {:e}\n", 1.0 / max(eta), 1.0 / min(eta)); + } else { + // transition from 0 for large P0 to resistivity for small P0 + eta = core_resist + (vac_resist - core_resist) * vac_mask; + } - if (constn0) { - T0_fake_prof = false; - n0_fake_prof = false; - } else { - Nbar = 1.0; - Tibar = 1000.0; - Tebar = 1000.0; - - if ((!T0_fake_prof) && n0_fake_prof) { - N0 = N0tanh(n0_height * Nbar, n0_ave * Nbar, n0_width, n0_center, n0_bottom_x); - - Ti0 = P0 / N0 / 2.0; - Te0 = Ti0; - } else if (T0_fake_prof) { - Ti0 = Tconst; - Te0 = Ti0; - N0 = P0 / (Ti0 + Te0); - } else { - if (mesh->get(N0, "Niexp")) { // N_i0 - throw BoutException("Error: Cannot read Ni0 from grid\n"); - } + eta = interp_to(eta, loc); + + SAVE_ONCE(eta); + + if (include_rmp) { + // Normalise RMP quantities + + rmp_Psi0 /= Bbar * Lbar; + + rmp_ramp /= Tbar; + rmp_freq *= Tbar; + rmp_rotate *= Tbar; + + rmp_Psi = rmp_Psi0; + rmp_dApdt = 0.0; + + bool apar_changing = false; + + output.write("Including magnetic perturbation\n"); + if (rmp_ramp > 0.0) { + output.write("\tRamping up over period t = {:e} ({:e} ms)\n", rmp_ramp, + rmp_ramp * Tbar * 1000.); + apar_changing = true; + } + if (rmp_freq > 0.0) { + output.write("\tOscillating with frequency f = {:e} ({:e} kHz)\n", rmp_freq, + rmp_freq / Tbar / 1000.); + apar_changing = true; + } + if (rmp_rotate != 0.0) { + output.write("\tRotating with a frequency f = {:e} ({:e} kHz)\n", rmp_rotate, + rmp_rotate / Tbar / 1000.); + apar_changing = true; + } + + if (apar_changing) { + SAVE_REPEAT(rmp_Psi, rmp_dApdt); + } else { + SAVE_ONCE(rmp_Psi); + } + } else { + rmp_Psi = 0.0; + } - if (mesh->get(Ti0, "Tiexp")) { // T_i0 - throw BoutException("Error: Cannot read Ti0 from grid\n"); - } + /**************** CALCULATE METRICS ******************/ - if (mesh->get(Te0, "Teexp")) { // T_e0 - throw BoutException("Error: Cannot read Te0 from grid\n"); - } - N0 /= Nbar; - Ti0 /= Tibar; - Te0 /= Tebar; - } - } + metric->g11 = SQ(Rxy * Bpxy); + metric->g22 = 1.0 / SQ(hthe); + metric->g33 = SQ(I) * metric->g11 + SQ(B0) / metric->g11; + metric->g12 = 0.0; + metric->g13 = -I * metric->g11; + metric->g23 = -Btxy / (hthe * Bpxy * Rxy); - if (gyroviscous) { - Dperp2Phi0.setLocation(CELL_CENTRE); - Dperp2Phi0.setBoundary("phi"); - Dperp2Phi.setLocation(CELL_CENTRE); - Dperp2Phi.setBoundary("phi"); - GradPhi02.setLocation(CELL_CENTRE); - GradPhi02.setBoundary("phi"); - GradcPhi.setLocation(CELL_CENTRE); - GradcPhi.setBoundary("phi"); - Dperp2Pi0.setLocation(CELL_CENTRE); - Dperp2Pi0.setBoundary("P"); - Dperp2Pi.setLocation(CELL_CENTRE); - Dperp2Pi.setBoundary("P"); - bracketPhi0P.setLocation(CELL_CENTRE); - bracketPhi0P.setBoundary("P"); - bracketPhiP0.setLocation(CELL_CENTRE); - bracketPhiP0.setBoundary("P"); - if (nonlinear) { - GradPhi2.setLocation(CELL_CENTRE); - GradPhi2.setBoundary("phi"); - bracketPhiP.setLocation(CELL_CENTRE); - bracketPhiP.setBoundary("P"); - } - } + metric->J = hthe / Bpxy; + metric->Bxy = B0; - BoutReal pnorm = max(P0, true); // Maximum over all processors + metric->g_11 = 1.0 / metric->g11 + SQ(I * Rxy); + metric->g_22 = SQ(B0 * hthe / Bpxy); + metric->g_33 = Rxy * Rxy; + metric->g_12 = Btxy * hthe * I * Rxy / Bpxy; + metric->g_13 = I * Rxy * Rxy; + metric->g_23 = Btxy * hthe * Rxy / Bpxy; - vacuum_pressure *= pnorm; // Get pressure from fraction - vacuum_trans *= pnorm; + metric->geometry(); // Calculate quantities from metric tensor - // Transitions from 0 in core to 1 in vacuum - vac_mask = (1.0 - tanh((P0 - vacuum_pressure) / vacuum_trans)) / 2.0; + // Set B field vector - if (spitzer_resist) { - // Use Spitzer resistivity - output.write("\tTemperature: {:e} -> {:e} [eV]\n", min(Te), max(Te)); - eta = 0.51 * 1.03e-4 * Zeff * 20. - * pow(Te, -1.5); // eta in Ohm-m. NOTE: ln(Lambda) = 20 - output.write("\tSpitzer resistivity: {:e} -> {:e} [Ohm m]\n", min(eta), max(eta)); - eta /= MU0 * Va * Lbar; - output.write("\t -> Lundquist {:e} -> {:e}\n", 1.0 / max(eta), 1.0 / min(eta)); - } else { - // transition from 0 for large P0 to resistivity for small P0 - eta = core_resist + (vac_resist - core_resist) * vac_mask; - } + B0vec.covariant = false; + B0vec.x = 0.; + B0vec.y = Bpxy / hthe; + B0vec.z = 0.; - eta = interp_to(eta, loc); + V0net.covariant = false; // presentation for net flow + V0net.x = 0.; + V0net.y = Rxy * Btxy * Bpxy / (hthe * B0 * B0) * Dphi0; + V0net.z = -Dphi0; - SAVE_ONCE(eta); + U0 = B0vec * Curl(V0net) / B0; // get 0th vorticity for Kelvin-Holmholtz term - if (include_rmp) { - // Normalise RMP quantities + /**************** SET EVOLVING VARIABLES *************/ - rmp_Psi0 /= Bbar * Lbar; + // Tell BOUT which variables to evolve + SOLVE_FOR(U, P); - rmp_ramp /= Tbar; - rmp_freq *= Tbar; - rmp_rotate *= Tbar; + if (evolve_jpar) { + output.write("Solving for jpar: Inverting to get Psi\n"); + SOLVE_FOR(Jpar); + SAVE_REPEAT(Psi); + } else { + output.write("Solving for Psi, Differentiating to get jpar\n"); + SOLVE_FOR(Psi); + SAVE_REPEAT(Jpar); + } - rmp_Psi = rmp_Psi0; - rmp_dApdt = 0.0; + if (compress) { + output.write("Including compression (Vpar) effects\n"); - bool apar_changing = false; + SOLVE_FOR(Vpar); + comms.add(Vpar); - output.write("Including magnetic perturbation\n"); - if (rmp_ramp > 0.0) { - output.write("\tRamping up over period t = {:e} ({:e} ms)\n", rmp_ramp, - rmp_ramp * Tbar * 1000.); - apar_changing = true; - } - if (rmp_freq > 0.0) { - output.write("\tOscillating with frequency f = {:e} ({:e} kHz)\n", rmp_freq, - rmp_freq / Tbar / 1000.); - apar_changing = true; - } - if (rmp_rotate != 0.0) { - output.write("\tRotating with a frequency f = {:e} ({:e} kHz)\n", rmp_rotate, - rmp_rotate / Tbar / 1000.); - apar_changing = true; - } + beta = B0 * B0 / (0.5 + (B0 * B0 / (g * P0))); + gradparB = Grad_par(B0) / B0; - if (apar_changing) { - SAVE_REPEAT(rmp_Psi, rmp_dApdt); - } else { - SAVE_ONCE(rmp_Psi); - } - } else { - rmp_Psi = 0.0; - } + output.write("Beta in range {:e} -> {:e}\n", min(beta), max(beta)); + } else { + Vpar = 0.0; + } - // Set B field vector + if (phi_constraint) { + // Implicit Phi solve using IDA - B0vec.covariant = false; - B0vec.x = 0.; + if (!solver->constraints()) { + throw BoutException("Cannot constrain. Run again with phi_constraint=false.\n"); + } - B0vec.y = Bpxy / hthe; - B0vec.z = 0.; + solver->constraint(phi, C_phi, "phi"); - V0net.covariant = false; // presentation for net flow - V0net.x = 0.; - V0net.y = Rxy * Btxy * Bpxy / (hthe * B0 * B0) * Dphi0; - V0net.z = -Dphi0; + // Set preconditioner + setPrecon(&ELMpb::precon_phi); - U0 = B0vec * Curl(V0net) / B0; // get 0th vorticity for Kelvin-Holmholtz term + } else { + // Phi solved in RHS (explicitly) + SAVE_REPEAT(phi); - /**************** SET EVOLVING VARIABLES *************/ + // Set preconditioner + setPrecon(&ELMpb::precon); - // Tell BOUT which variables to evolve - SOLVE_FOR(U, P); + // Set Jacobian + setJacobian((jacobianfunc)&ELMpb::jacobian); + } - if (evolve_jpar) { - output.write("Solving for jpar: Inverting to get Psi\n"); - SOLVE_FOR(Jpar); - SAVE_REPEAT(Psi); - } else { - output.write("Solving for Psi, Differentiating to get jpar\n"); - SOLVE_FOR(Psi); - SAVE_REPEAT(Jpar); - } + // Diamagnetic phi0 + if (diamag_phi0) { + if (constn0) { + phi0 = -0.5 * dnorm * P0 / B0; + } else { + // Stationary equilibrium plasma. ExB velocity balances diamagnetic drift + phi0 = -0.5 * dnorm * P0 / B0 / N0; + } + SAVE_ONCE(phi0); + } - if (compress) { - output.write("Including compression (Vpar) effects\n"); + // Add some equilibrium quantities and normalisations + // everything needed to recover physical units + SAVE_ONCE(J0, P0); + SAVE_ONCE(density, Lbar, Bbar, Tbar); + SAVE_ONCE(Va, B0); + SAVE_ONCE(Dphi0, U0); + SAVE_ONCE(V0); + if (!constn0) { + SAVE_ONCE(Ti0, Te0, N0); + } - SOLVE_FOR(Vpar); - comms.add(Vpar); + // Create a solver for the Laplacian + phiSolver = Laplacian::create(&globalOptions["phiSolver"]); + // Save performance metrics to output, using the + // given name as the prefix. + phiSolver->savePerformance(*solver, "phiSolver"); + + aparSolver = Laplacian::create(&globalOptions["aparSolver"], loc); + + if (phi_boundary_relax) { + // Set the last update time to -1, so it will reset + // the first time RHS function is called + phi_boundary_last_update = -1.; + phi_core_averagey = options["phi_core_averagey"] + .doc("Average phi core boundary in Y?") + .withDefault(false) + and mesh->periodicY(mesh->xstart); + + phi_boundary_timescale = options["phi_boundary_timescale"] + .doc("Timescale for phi boundary relaxation [seconds]") + .withDefault(1e-7) + / Tbar; // Normalise time units to Tbar + + phiSolver->setInnerBoundaryFlags(INVERT_SET); + phiSolver->setOuterBoundaryFlags(INVERT_SET); + } - beta = B0 * B0 / (0.5 + (B0 * B0 / (g * P0))); - gradparB = Grad_par(B0) / B0; + /////////////// CHECK VACUUM /////////////////////// + // In vacuum region, initial vorticity should equal zero - output.write("Beta in range {:e} -> {:e}\n", min(beta), max(beta)); - } else { - Vpar = 0.0; - } + if (!restarting) { + // Only if not restarting: Check initial perturbation - if (phi_constraint) { - // Implicit Phi solve using IDA + // Set U to zero where P0 < vacuum_pressure + U = where(P0 - vacuum_pressure, U, 0.0); - if (!solver->constraints()) { - throw BoutException("Cannot constrain. Run again with phi_constraint=false.\n"); - } + if (constn0) { + ubyn = U; + // Phi should be consistent with U + phi = phiSolver->solve(ubyn); + } else { + ubyn = U / N0; + phiSolver->setCoefC(N0); + phi = phiSolver->solve(ubyn); + } - solver->constraint(phi, C_phi, "phi"); + // if(diamag) { + // phi -= 0.5*dnorm * P / B0; + //} + } - // Set preconditioner - setPrecon(&ELMpb::precon_phi); + /************** SETUP COMMUNICATIONS **************/ - } else { - // Phi solved in RHS (explicitly) - SAVE_REPEAT(phi); + comms.add(U, P); - // Set preconditioner - setPrecon(&ELMpb::precon); + phi.setBoundary("phi"); // Set boundary conditions + tmpU2.setBoundary("U"); + tmpP2.setBoundary("P"); + tmpA2.setBoundary("J"); - // Set Jacobian - setJacobian((jacobianfunc) & ELMpb::jacobian); - } + if (evolve_jpar) { + comms.add(Jpar); + } else { + comms.add(Psi); + // otherwise Need to communicate Jpar separately + Jpar.setBoundary("J"); + } + Jpar2.setBoundary("J"); - // Diamagnetic phi0 - if (diamag_phi0) { - if (constn0) { - phi0 = -0.5 * dnorm * P0 / B0; - } else { - // Stationary equilibrium plasma. ExB velocity balances diamagnetic drift - phi0 = -0.5 * dnorm * P0 / B0 / N0; - } - SAVE_ONCE(phi0); - } + return 0; + } - // Add some equilibrium quantities and normalisations - // everything needed to recover physical units - SAVE_ONCE(J0, P0); - SAVE_ONCE(density, Lbar, Bbar, Tbar); - SAVE_ONCE(Va, B0); - SAVE_ONCE(Dphi0, U0); - SAVE_ONCE(V0); - if (!constn0) { - SAVE_ONCE(Ti0, Te0, N0); - } + // Parallel gradient along perturbed field-line + Field3D Grad_parP(const Field3D& f, CELL_LOC loc = CELL_DEFAULT) const { - // Create a solver for the Laplacian - phiSolver = Laplacian::create(&globalOptions["phiSolver"]); + if (loc == CELL_DEFAULT) { + loc = f.getLocation(); + } - aparSolver = Laplacian::create(&globalOptions["aparSolver"], loc); + Field3D result = Grad_par(f, loc); - /////////////// CHECK VACUUM /////////////////////// - // In vacuum region, initial vorticity should equal zero + if (nonlinear) { + result -= bracket(interp_to(Psi, loc), f, bm_mag) * B0; - if (!restarting) { - // Only if not restarting: Check initial perturbation + if (include_rmp) { + result -= bracket(interp_to(rmp_Psi, loc), f, bm_mag) * B0; + } + } - // Set U to zero where P0 < vacuum_pressure - U = where(P0 - vacuum_pressure, U, 0.0); + return result; + } - if (constn0) { - ubyn = U; - // Phi should be consistent with U - phi = phiSolver->solve(ubyn); - } else { - ubyn = U / N0; - phiSolver->setCoefC(N0); - phi = phiSolver->solve(ubyn); - } + bool first_run = true; // For printing out some diagnostics first time around - // if(diamag) { - // phi -= 0.5*dnorm * P / B0; - //} - } + int rhs(BoutReal t) override { + // Perform communications + mesh->communicate(comms); - /************** SETUP COMMUNICATIONS **************/ + Coordinates* metric = mesh->getCoordinates(); - comms.add(U, P); + //////////////////////////////////////////// + // Transitions from 0 in core to 1 in vacuum + if (nonlinear) { + vac_mask = (1.0 - tanh(((P0 + P) - vacuum_pressure) / vacuum_trans)) / 2.0; - phi.setBoundary("phi"); // Set boundary conditions - tmpU2.setBoundary("U"); - tmpP2.setBoundary("P"); - tmpA2.setBoundary("J"); + // Update resistivity + if (spitzer_resist) { + // Use Spitzer formula + Field3D Te; + Te = (P0 + P) * Bbar * Bbar / (4. * MU0) / (density * 1.602e-19); // eV - if (evolve_jpar) { - comms.add(Jpar); - } else { - comms.add(Psi); - // otherwise Need to communicate Jpar separately - Jpar.setBoundary("J"); - } - Jpar2.setBoundary("J"); + // eta in Ohm-m. ln(Lambda) = 20 + eta = interp_to(0.51 * 1.03e-4 * Zeff * 20. * pow(Te, -1.5), loc); - return 0; + // Normalised eta + eta /= MU0 * Va * Lbar; + } else { + // Use specified core and vacuum Lundquist numbers + eta = core_resist + (vac_resist - core_resist) * vac_mask; + } + eta = interp_to(eta, loc); } - // Parallel gradient along perturbed field-line - Field3D Grad_parP(const Field3D &f, CELL_LOC loc = CELL_DEFAULT) const { + //////////////////////////////////////////// + // Resonant Magnetic Perturbation code - if (loc == CELL_DEFAULT) { - loc = f.getLocation(); - } + if (include_rmp) { - Field3D result = Grad_par(f, loc); + if ((rmp_ramp > 0.0) || (rmp_freq > 0.0) || (rmp_rotate != 0.0)) { + // Need to update the RMP terms - auto B0 = tokamak_options.Bxy; + if ((rmp_ramp > 0.0) && (t < rmp_ramp)) { + // Still in ramp phase - if (nonlinear) { - result -= bracket(interp_to(Psi, loc), f, bm_mag) * B0; + rmp_Psi = (t / rmp_ramp) * rmp_Psi0; // Linear ramp - if (include_rmp) { - result -= bracket(interp_to(rmp_Psi, loc), f, bm_mag) * B0; - } + rmp_dApdt = rmp_Psi0 / rmp_ramp; + } else { + rmp_Psi = rmp_Psi0; + rmp_dApdt = 0.0; } - return result; - } - - bool first_run = true; // For printing out some diagnostics first time around - - int rhs(BoutReal t) override { - // Perform communications - mesh->communicate(comms); - - Coordinates *metric = mesh->getCoordinates(); + if (rmp_freq > 0.0) { + // Oscillating the amplitude - auto B0 = tokamak_options.Bxy; + rmp_dApdt = rmp_dApdt * sin(2. * PI * rmp_freq * t) + + rmp_Psi * (2. * PI * rmp_freq) * cos(2. * PI * rmp_freq * t); - //////////////////////////////////////////// - // Transitions from 0 in core to 1 in vacuum - if (nonlinear) { - vac_mask = (1.0 - tanh(((P0 + P) - vacuum_pressure) / vacuum_trans)) / 2.0; - - // Update resistivity - if (spitzer_resist) { - // Use Spitzer formula - Field3D Te; - Te = (P0 + P) * Bbar * Bbar / (4. * MU0) / (density * 1.602e-19); // eV - - // eta in Ohm-m. ln(Lambda) = 20 - eta = interp_to(0.51 * 1.03e-4 * Zeff * 20. * pow(Te, -1.5), loc); - - // Normalised eta - eta /= MU0 * Va * Lbar; - } else { - // Use specified core and vacuum Lundquist numbers - eta = core_resist + (vac_resist - core_resist) * vac_mask; - } - eta = interp_to(eta, loc); + rmp_Psi *= sin(2. * PI * rmp_freq * t); } - //////////////////////////////////////////// - // Resonant Magnetic Perturbation code - - if (include_rmp) { - - if ((rmp_ramp > 0.0) || (rmp_freq > 0.0) || (rmp_rotate != 0.0)) { - // Need to update the RMP terms + if (rmp_rotate != 0.0) { + // Rotate toroidally at given frequency - if ((rmp_ramp > 0.0) && (t < rmp_ramp)) { - // Still in ramp phase + shiftZ(rmp_Psi, 2 * PI * rmp_rotate * t); + shiftZ(rmp_dApdt, 2 * PI * rmp_rotate * t); - rmp_Psi = (t / rmp_ramp) * rmp_Psi0; // Linear ramp + // Add toroidal rotation term. CHECK SIGN - rmp_dApdt = rmp_Psi0 / rmp_ramp; - } else { - rmp_Psi = rmp_Psi0; - rmp_dApdt = 0.0; - } - - if (rmp_freq > 0.0) { - // Oscillating the amplitude - - rmp_dApdt = rmp_dApdt * sin(2. * PI * rmp_freq * t) - + rmp_Psi * (2. * PI * rmp_freq) * cos(2. * PI * rmp_freq * t); - - rmp_Psi *= sin(2. * PI * rmp_freq * t); - } - - if (rmp_rotate != 0.0) { - // Rotate toroidally at given frequency - - shiftZ(rmp_Psi, 2 * PI * rmp_rotate * t); - shiftZ(rmp_dApdt, 2 * PI * rmp_rotate * t); - - // Add toroidal rotation term. CHECK SIGN - - rmp_dApdt += DDZ(rmp_Psi) * 2 * PI * rmp_rotate; - } - - // Set to zero in the core - if (rmp_vac_mask) { - rmp_Psi *= vac_mask; - } - } else { - // Set to zero in the core region - if (rmp_vac_mask) { - // Only in vacuum -> skin current -> diffuses inwards - rmp_Psi = rmp_Psi0 * vac_mask; - } - } - - mesh->communicate(rmp_Psi); + rmp_dApdt += DDZ(rmp_Psi) * 2 * PI * rmp_rotate; } - //////////////////////////////////////////// - // Inversion - - if (evolve_jpar) { - // Invert laplacian for Psi - Psi = aparSolver->solve(Jpar); - mesh->communicate(Psi); + // Set to zero in the core + if (rmp_vac_mask) { + rmp_Psi *= vac_mask; } + } else { + // Set to zero in the core region + if (rmp_vac_mask) { + // Only in vacuum -> skin current -> diffuses inwards + rmp_Psi = rmp_Psi0 * vac_mask; + } + } - if (phi_constraint) { - // Phi being solved as a constraint - - Field3D Ctmp = phi; - Ctmp.setBoundary("phi"); // Look up boundary conditions for phi - Ctmp.applyBoundary(); - Ctmp -= phi; // Now contains error in the boundary - - C_phi = Delp2(phi) - U; // Error in the bulk - C_phi.setBoundaryTo(Ctmp); - - } else { - - if (constn0) { - if (split_n0) { - //////////////////////////////////////////// - // Boussinesq, split - // Split into axisymmetric and non-axisymmetric components - Field2D Vort2D = DC(U); // n=0 component - - // Applies boundary condition for "phi". - phi2D.applyBoundary(t); - - // Solve axisymmetric (n=0) part - phi2D = laplacexy->solve(Vort2D, phi2D); + mesh->communicate(rmp_Psi); + } - // Solve non-axisymmetric part - phi = phiSolver->solve(U - Vort2D); + //////////////////////////////////////////// + // Inversion - phi += phi2D; // Add axisymmetric part - } else { - phi = phiSolver->solve(U); - } + if (evolve_jpar) { + // Invert laplacian for Psi + Psi = aparSolver->solve(Jpar); + mesh->communicate(Psi); + } - if (diamag) { - phi -= 0.5 * dnorm * P / B0; - } - } else { - ubyn = U / N0; - if (diamag) { - ubyn -= 0.5 * dnorm / (N0 * B0) * Delp2(P); - mesh->communicate(ubyn); + if (phi_constraint) { + // Phi being solved as a constraint + + Field3D Ctmp = phi; + Ctmp.setBoundary("phi"); // Look up boundary conditions for phi + Ctmp.applyBoundary(); + Ctmp -= phi; // Now contains error in the boundary + + if (laplace_perp) { + C_phi = Laplace_perp(phi) - U; // Error in the bulk + } else { + C_phi = Delp2(phi) - U; // Error in the bulk + } + C_phi.setBoundaryTo(Ctmp); + + } else { + if (phi_boundary_relax) { + // Update the boundary regions by relaxing towards zero gradient + // on a given timescale. + + if (phi_boundary_last_update < 0.0) { + // First time this has been called. + phi_boundary_last_update = t; + } else if (t > phi_boundary_last_update) { + // Only update if simulation time has advanced + // Uses an exponential decay of the weighting of the value in the boundary + // so that the solution is well behaved for arbitrary steps + BoutReal const weight = exp(-(t - phi_boundary_last_update) / phi_boundary_timescale); + phi_boundary_last_update = t; + + if (mesh->firstX()) { + BoutReal phivalue = 0.0; + if (phi_core_averagey) { + // Calculate a single phi boundary value for all Y slices + BoutReal philocal = 0.0; + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + philocal += phi(mesh->xstart, j, k); } - // Invert laplacian for phi - phiSolver->setCoefC(N0); - phi = phiSolver->solve(ubyn); - } - // Apply a boundary condition on phi for target plates - phi.applyBoundary(); - mesh->communicate(phi); - } - - if (!evolve_jpar) { - // Get J from Psi - Jpar = Delp2(Psi); - if (include_rmp) { - Jpar += Delp2(rmp_Psi); + } + MPI_Comm comm_inner = mesh->getYcomm(0); + int np = 0; + MPI_Comm_size(comm_inner, &np); + MPI_Allreduce(&philocal, &phivalue, 1, MPI_DOUBLE, MPI_SUM, comm_inner); + phivalue /= (np * mesh->LocalNz * mesh->LocalNy); } - - Jpar.applyBoundary(); - mesh->communicate(Jpar); - - if (jpar_bndry_width > 0) { - // Zero j in boundary regions. Prevents vorticity drive - // at the boundary - - for (int i = 0; i < jpar_bndry_width; i++) { - for (int j = 0; j < mesh->LocalNy; j++) { - for (int k = 0; k < mesh->LocalNz; k++) { - if (mesh->firstX()) { - Jpar(i, j, k) = 0.0; - } - if (mesh->lastX()) { - Jpar(mesh->LocalNx - 1 - i, j, k) = 0.0; - } - } - } + for (int j = mesh->ystart; j <= mesh->yend; j++) { + if (!phi_core_averagey) { + phivalue = 0.0; // Calculate phi boundary for each Y index separately + for (int k = 0; k < mesh->LocalNz; k++) { + phivalue += phi(mesh->xstart, j, k); } - } + phivalue /= mesh->LocalNz; // Average in Z of point next to boundary + } - // Smooth j in x - if (smooth_j_x) { - Jpar = smooth_x(Jpar); - Jpar.applyBoundary(); + // Old value of phi at boundary. Note: this is constant in Z + BoutReal const oldvalue = + 0.5 * (phi(mesh->xstart - 1, j, 0) + phi(mesh->xstart, j, 0)); - // Recommunicate now smoothed - mesh->communicate(Jpar); - } + // New value of phi at boundary, relaxing towards phivalue + BoutReal const newvalue = weight * oldvalue + (1. - weight) * phivalue; - // Get Delp2(J) from J - Jpar2 = Delp2(Jpar); - - Jpar2.applyBoundary(); - mesh->communicate(Jpar2); - - if (jpar_bndry_width > 0) { - // Zero jpar2 in boundary regions. Prevents vorticity drive - // at the boundary - - for (int i = 0; i < jpar_bndry_width; i++) { - for (int j = 0; j < mesh->LocalNy; j++) { - for (int k = 0; k < mesh->LocalNz; k++) { - if (mesh->firstX()) { - Jpar2(i, j, k) = 0.0; - } - if (mesh->lastX()) { - Jpar2(mesh->LocalNx - 1 - i, j, k) = 0.0; - } - } - } - } + // Set phi at the boundary to this value + for (int k = 0; k < mesh->LocalNz; k++) { + phi(mesh->xstart - 1, j, k) = 2. * newvalue - phi(mesh->xstart, j, k); + phi(mesh->xstart - 2, j, k) = phi(mesh->xstart - 1, j, k); + } } - } - - //////////////////////////////////////////////////// - // Sheath boundary conditions - // Normalised and linearised, since here we have only pressure - // rather than density and temperature. Applying a boundary - // to Jpar so that Jpar = sqrt(mi/me)/(2*pi) * phi - // + } + + if (mesh->lastX()) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + BoutReal phivalue = 0.0; + for (int k = 0; k < mesh->LocalNz; k++) { + phivalue += phi(mesh->xend, j, k); + } + phivalue /= mesh->LocalNz; // Average in Z of point next to boundary + + // Old value of phi at boundary. Note: this is constant in Z + BoutReal oldvalue = + 0.5 * (phi(mesh->xend + 1, j, 0) + phi(mesh->xend, j, 0)); + + // New value of phi at boundary, relaxing towards phivalue + BoutReal const newvalue = weight * oldvalue + (1. - weight) * phivalue; + + // Set phi at the boundary to this value + for (int k = 0; k < mesh->LocalNz; k++) { + phi(mesh->xend + 1, j, k) = 2. * newvalue - phi(mesh->xend, j, k); + phi(mesh->xend + 2, j, k) = phi(mesh->xend + 1, j, k); + } + } + } + } + } + + Field3D phi_shift = phi; + if (constn0 and diamag) { + // Solving for phi + ion pressure term + phi_shift += 0.5 * dnorm * P / B0; + } else { + // Ensure that memory is not shared between phi and phi_shift + phi_shift.allocate(); + } + + // Update boundary conditions. + // The INVERT_SET flag takes the value in the guard (boundary) cell + // and sets the boundary between cells to this value. + // This shift by 1/2 grid cell is important. + + if (mesh->firstX()) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + // Average phi + Pi at the boundary, and set the boundary cell + // to this value. The phi solver will then put the value back + // onto the cell mid-point + phi_shift(mesh->xstart - 1, j, k) = + 0.5 * (phi_shift(mesh->xstart - 1, j, k) + phi_shift(mesh->xstart, j, k)); + } + } + } + + if (mesh->lastX()) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + phi_shift(mesh->xend + 1, j, k) = + 0.5 * (phi_shift(mesh->xend + 1, j, k) + phi_shift(mesh->xend, j, k)); + } + } + } + + if (constn0) { + if (split_n0) { + //////////////////////////////////////////// + // Boussinesq, split + // Split into axisymmetric and non-axisymmetric components + Field2D Vort2D = DC(U); // n=0 component + Field2D phi_shift_2d = phi2D; - if (sheath_boundaries) { + if (phi_boundary_relax) { + phi_shift_2d = DC(phi_shift); + } + phi_shift -= phi_shift_2d; - // Need to shift into field-aligned coordinates before applying - // parallel boundary conditions + // Applies boundary condition for "phi". + phi2D.applyBoundary(t); - auto phi_fa = toFieldAligned(phi); - auto P_fa = toFieldAligned(P); - auto Jpar_fa = toFieldAligned(Jpar); + // Solve axisymmetric (n=0) part + phi2D = laplacexy->solve(Vort2D, phi_shift_2d); - // At y = ystart (lower boundary) + // Solve non-axisymmetric part + phi = phiSolver->solve(U - Vort2D, phi_shift); - for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { + phi += phi2D; // Add axisymmetric part + } else { + if (phi_boundary_relax) { + phi = phiSolver->solve(U, phi_shift); + } else { + phi = phiSolver->solve(U); + } + } + if (diamag) { + phi -= 0.5 * dnorm * P / B0; + } + } else { + ubyn = U / N0; + if (diamag) { + if (laplace_perp) { + ubyn -= 0.5 * dnorm / (N0 * B0) * Laplace_perp(P); + } else { + ubyn -= 0.5 * dnorm / (N0 * B0) * Delp2(P); + } + mesh->communicate(ubyn); + } + // Invert laplacian for phi + phiSolver->setCoefC(N0); + phi = phiSolver->solve(ubyn); + } + mesh->communicate(phi); + } - // Zero-gradient potential - BoutReal const phisheath = phi_fa(r.ind, mesh->ystart, jz); + if (mesh->firstX()) { + for (int i = mesh->xstart - 2; i >= 0; --i) { + for (int j = mesh->ystart; j <= mesh->yend; ++j) { + for (int k = 0; k < mesh->LocalNz; ++k) { + phi(i, j, k) = phi(i + 1, j, k); + } + } + } + } - BoutReal jsheath = -(sqrt(mi_me) / (2. * sqrt(PI))) * phisheath; + if (mesh->lastX()) { + for (int i = mesh->xend + 2; i < mesh->LocalNx; ++i) { + for (int j = mesh->ystart; j <= mesh->yend; ++j) { + for (int k = 0; k < mesh->LocalNz; ++k) { + phi(i, j, k) = phi(i - 1, j, k); + } + } + } + } - // Apply boundary condition half-way between cells - for (int jy = mesh->ystart - 1; jy >= 0; jy--) { - // Neumann conditions - P_fa(r.ind, jy, jz) = P_fa(r.ind, mesh->ystart, jz); - phi_fa(r.ind, jy, jz) = phisheath; - // Dirichlet condition on Jpar - Jpar_fa(r.ind, jy, jz) = 2. * jsheath - Jpar_fa(r.ind, mesh->ystart, jz); - } - } + if (!evolve_jpar) { + // Get J from Psi + if (laplace_perp) { + Jpar = Laplace_perp(Psi); + } else { + Jpar = Delp2(Psi); + } + + if (include_rmp) { + if (laplace_perp) { + Jpar += Laplace_perp(rmp_Psi); + } else { + Jpar += Delp2(rmp_Psi); + } + } + + Jpar.applyBoundary(); + mesh->communicate(Jpar); + + if (jpar_bndry_width > 0) { + // Zero j in boundary regions. Prevents vorticity drive + // at the boundary + + for (int i = 0; i < jpar_bndry_width; i++) { + for (int j = 0; j < mesh->LocalNy; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + if (mesh->firstX()) { + Jpar(i, j, k) = 0.0; + } + if (mesh->lastX()) { + Jpar(mesh->LocalNx - 1 - i, j, k) = 0.0; + } } - - // At y = yend (upper boundary) - - for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { - - // Zero-gradient potential - BoutReal const phisheath = phi_fa(r.ind, mesh->yend, jz); - - BoutReal jsheath = (sqrt(mi_me) / (2. * sqrt(PI))) * phisheath; - - // Apply boundary condition half-way between cells - for (int jy = mesh->yend + 1; jy < mesh->LocalNy; jy++) { - // Neumann conditions - P_fa(r.ind, jy, jz) = P_fa(r.ind, mesh->yend, jz); - phi_fa(r.ind, jy, jz) = phisheath; - // Dirichlet condition on Jpar - // WARNING: this is not correct if staggered grids are used - ASSERT3(not mesh->StaggerGrids); - Jpar_fa(r.ind, jy, jz) = 2. * jsheath - Jpar_fa(r.ind, mesh->yend, jz); - } - } + } + } + } + + // Smooth j in x + if (smooth_j_x) { + Jpar = smooth_x(Jpar); + Jpar.applyBoundary(); + + // Recommunicate now smoothed + mesh->communicate(Jpar); + } + + // Get Delp2(J) from J + if (laplace_perp) { + Jpar2 = Laplace_perp(Jpar); + } else { + Jpar2 = Delp2(Jpar); + } + Jpar2.applyBoundary(); + mesh->communicate(Jpar2); + + if (jpar_bndry_width > 0) { + // Zero jpar2 in boundary regions. Prevents vorticity drive + // at the boundary + + for (int i = 0; i < jpar_bndry_width; i++) { + for (int j = 0; j < mesh->LocalNy; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + if (mesh->firstX()) { + Jpar2(i, j, k) = 0.0; + } + if (mesh->lastX()) { + Jpar2(mesh->LocalNx - 1 - i, j, k) = 0.0; + } } - - // Shift back from field aligned coordinates - phi = fromFieldAligned(phi_fa); - P = fromFieldAligned(P_fa); - Jpar = fromFieldAligned(Jpar_fa); + } } + } + } - //////////////////////////////////////////////////// - // Parallel electric field + //////////////////////////////////////////////////// + // Sheath boundary conditions + // Normalised and linearised, since here we have only pressure + // rather than density and temperature. Applying a boundary + // to Jpar so that Jpar = sqrt(mi/me)/(2*pi) * phi + // - if (evolve_jpar) { - // Jpar - Field3D B0U = B0 * U; - mesh->communicate(B0U); - ddt(Jpar) = -Grad_parP(B0U, loc) / B0 + eta * Delp2(Jpar); + if (sheath_boundaries) { - if (relax_j_vac) { - // Make ddt(Jpar) relax to zero. + // Need to shift into field-aligned coordinates before applying + // parallel boundary conditions - ddt(Jpar) -= vac_mask * Jpar / relax_j_tconst; - } - } else { - // Vector potential - ddt(Psi) = -Grad_parP(phi * B0, loc) / B0 + eta * Jpar; + auto phi_fa = toFieldAligned(phi); + auto P_fa = toFieldAligned(P); + auto Jpar_fa = toFieldAligned(Jpar); - if (eHall) { // electron parallel pressure - ddt(Psi) += 0.25 * delta_i - * (Grad_parP(B0 * P, loc) / B0 - + bracket(interp_to(P0, loc), Psi, bm_mag) * B0); - } + // At y = ystart (lower boundary) - if (diamag_phi0) { // Equilibrium flow - ddt(Psi) -= bracket(interp_to(phi0, loc), Psi, bm_exb) * B0; - } + for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { - if (withflow) { // net flow - ddt(Psi) -= V_dot_Grad(V0net, Psi); - } + // Zero-gradient potential + BoutReal const phisheath = phi_fa(r.ind, mesh->ystart, jz); - if (diamag_grad_t) { // grad_par(T_e) correction - ddt(Psi) += 1.71 * dnorm * 0.5 * Grad_parP(P, loc) / B0; - } - - if (hyperresist > 0.0) { // Hyper-resistivity - ddt(Psi) -= eta * hyperresist * Delp2(Jpar); - } + BoutReal jsheath = -(sqrt(mi_me) / (2. * sqrt(PI))) * phisheath; - if (ehyperviscos > 0.0) { // electron Hyper-viscosity coefficient - ddt(Psi) -= eta * ehyperviscos * Delp2(Jpar2); - } + // Apply boundary condition half-way between cells + for (int jy = mesh->ystart - 1; jy >= 0; jy--) { + // Neumann conditions + P_fa(r.ind, jy, jz) = P_fa(r.ind, mesh->ystart, jz); + phi_fa(r.ind, jy, jz) = phisheath; + // Dirichlet condition on Jpar + Jpar_fa(r.ind, jy, jz) = 2. * jsheath - Jpar_fa(r.ind, mesh->ystart, jz); + } + } + } - // Parallel hyper-viscous diffusion for vector potential - if (diffusion_a4 > 0.0) { - tmpA2 = D2DY2(Psi); - mesh->communicate(tmpA2); - tmpA2.applyBoundary(); - ddt(Psi) -= diffusion_a4 * D2DY2(tmpA2); - } + // At y = yend (upper boundary) - // Vacuum solution - if (relax_j_vac) { - // Calculate the J and Psi profile we're aiming for - Field3D Jtarget = Jpar * (1.0 - vac_mask); // Zero in vacuum + for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { - // Invert laplacian for Psi - Field3D Psitarget = aparSolver->solve(Jtarget); + // Zero-gradient potential + BoutReal const phisheath = phi_fa(r.ind, mesh->yend, jz); - // Add a relaxation term in the vacuum - ddt(Psi) = - ddt(Psi) * (1. - vac_mask) - (Psi - Psitarget) * vac_mask / relax_j_tconst; - } - } + BoutReal jsheath = (sqrt(mi_me) / (2. * sqrt(PI))) * phisheath; - //////////////////////////////////////////////////// - // Vorticity equation - Psi_loc = interp_to(Psi, CELL_CENTRE, "RGN_ALL"); - Psi_loc.applyBoundary(); - // Grad j term - ddt(U) = SQ(B0) * b0xGrad_dot_Grad(Psi_loc, J0, CELL_CENTRE); - if (include_rmp) { - ddt(U) += SQ(B0) * b0xGrad_dot_Grad(rmp_Psi, J0, CELL_CENTRE); + // Apply boundary condition half-way between cells + for (int jy = mesh->yend + 1; jy < mesh->LocalNy; jy++) { + // Neumann conditions + P_fa(r.ind, jy, jz) = P_fa(r.ind, mesh->yend, jz); + phi_fa(r.ind, jy, jz) = phisheath; + // Dirichlet condition on Jpar + // WARNING: this is not correct if staggered grids are used + ASSERT3(not mesh->StaggerGrids); + Jpar_fa(r.ind, jy, jz) = 2. * jsheath - Jpar_fa(r.ind, mesh->yend, jz); + } } + } - ddt(U) += b0xcv * Grad(P); // curvature term + // Shift back from field aligned coordinates + phi = fromFieldAligned(phi_fa); + P = fromFieldAligned(P_fa); + Jpar = fromFieldAligned(Jpar_fa); + } - if (!nogradparj) { // Parallel current term - ddt(U) -= SQ(B0) * Grad_parP(Jpar, CELL_CENTRE); // b dot grad j + //////////////////////////////////////////////////// + // Parallel electric field + + if (evolve_jpar) { + // Jpar + Field3D B0U = B0 * U; + mesh->communicate(B0U); + if (laplace_perp) { + ddt(Jpar) = -Grad_parP(B0U, loc) / B0 + eta * Laplace_perp(Jpar); + } else { + ddt(Jpar) = -Grad_parP(B0U, loc) / B0 + eta * Delp2(Jpar); + } + + if (relax_j_vac) { + // Make ddt(Jpar) relax to zero. + + ddt(Jpar) -= vac_mask * Jpar / relax_j_tconst; + } + } else { + // Vector potential + ddt(Psi) = -Grad_parP(phi * B0, loc) / B0 + eta * Jpar; + + if (eHall) { // electron parallel pressure + ddt(Psi) += 0.25 * delta_i + * (Grad_parP(B0 * P, loc) / B0 + + bracket(interp_to(P0, loc), Psi, bm_mag) * B0); + } + + if (diamag_phi0) { // Equilibrium flow + ddt(Psi) -= bracket(interp_to(phi0, loc), Psi, bm_exb) * B0; + } + + if (withflow) { // net flow + ddt(Psi) -= V_dot_Grad(V0net, Psi); + } + + if (diamag_grad_t) { // grad_par(T_e) correction + ddt(Psi) += 1.71 * dnorm * 0.5 * Grad_parP(P, loc) / B0; + } + + if (hyperresist > 0.0) { // Hyper-resistivity + if (laplace_perp) { + ddt(Psi) -= eta * hyperresist * Laplace_perp(Jpar); + } else { + ddt(Psi) -= eta * hyperresist * Delp2(Jpar); } + } - if (withflow && K_H_term) { // K_H_term - ddt(U) -= b0xGrad_dot_Grad(phi, U0); - } + if (ehyperviscos > 0.0) { // electron Hyper-viscosity coefficient + if (laplace_perp) { + ddt(Psi) -= eta * ehyperviscos * Laplace_perp(Jpar2); + } else { + ddt(Psi) -= eta * ehyperviscos * Delp2(Jpar2); + } + } + + // Parallel hyper-viscous diffusion for vector potential + if (diffusion_a4 > 0.0) { + tmpA2 = D2DY2(Psi); + mesh->communicate(tmpA2); + tmpA2.applyBoundary(); + ddt(Psi) -= diffusion_a4 * D2DY2(tmpA2); + } + + // Vacuum solution + if (relax_j_vac) { + // Calculate the J and Psi profile we're aiming for + Field3D Jtarget = Jpar * (1.0 - vac_mask); // Zero in vacuum + + // Invert laplacian for Psi + Field3D Psitarget = aparSolver->solve(Jtarget); + + // Add a relaxation term in the vacuum + ddt(Psi) = + ddt(Psi) * (1. - vac_mask) - (Psi - Psitarget) * vac_mask / relax_j_tconst; + } + } - if (diamag_phi0) { // Equilibrium flow - ddt(U) -= b0xGrad_dot_Grad(phi0, U); - } + //////////////////////////////////////////////////// + // Vorticity equation + Psi_loc = interp_to(Psi, CELL_CENTRE, "RGN_ALL"); + Psi_loc.applyBoundary(); + // Grad j term + ddt(U) = SQ(B0) * b0xGrad_dot_Grad(Psi_loc, J0, CELL_CENTRE); + if (include_rmp) { + ddt(U) += SQ(B0) * b0xGrad_dot_Grad(rmp_Psi, J0, CELL_CENTRE); + } - if (withflow) { // net flow - ddt(U) -= V_dot_Grad(V0net, U); - } + ddt(U) += b0xcv * Grad(P); // curvature term - if (nonlinear) { // Advection - ddt(U) -= bracket(phi, U, bm_exb) * B0; - } + if (!nogradparj) { // Parallel current term + ddt(U) -= SQ(B0) * Grad_parP(Jpar, CELL_CENTRE); // b dot grad j + } - // Viscosity terms + if (withflow && K_H_term) { // K_H_term + ddt(U) -= b0xGrad_dot_Grad(phi, U0); + } - if (viscos_par > 0.0) { // Parallel viscosity - ddt(U) += viscos_par * Grad2_par2(U); - } + if (diamag_phi0) { // Equilibrium flow + ddt(U) -= b0xGrad_dot_Grad(phi0, U); + } - if (diffusion_u4 > 0.0) { - tmpU2 = D2DY2(U); - mesh->communicate(tmpU2); - tmpU2.applyBoundary(); - ddt(U) -= diffusion_u4 * D2DY2(tmpU2); - } + if (withflow) { // net flow + ddt(U) -= V_dot_Grad(V0net, U); + } - if (viscos_perp > 0.0) { - ddt(U) += viscos_perp * Delp2(U); // Perpendicular viscosity - } + if (nonlinear) { // Advection + ddt(U) -= bracket(phi, U, bm_exb) * B0; + } - // Hyper-viscosity - if (hyperviscos > 0.0) { - // Calculate coefficient. + // Viscosity terms - hyper_mu_x = hyperviscos * metric->g_11() * SQ(metric->dx()) - * abs(metric->g11() * D2DX2(U)) / (abs(U) + 1e-3); - hyper_mu_x.applyBoundary("dirichlet"); // Set to zero on all boundaries + if (viscos_par > 0.0) { // Parallel viscosity + ddt(U) += viscos_par * Grad2_par2(U); + } - ddt(U) += hyper_mu_x * metric->g11() * D2DX2(U); + if (diffusion_u4 > 0.0) { + tmpU2 = D2DY2(U); + mesh->communicate(tmpU2); + tmpU2.applyBoundary(); + ddt(U) -= diffusion_u4 * D2DY2(tmpU2); + } - if (first_run) { // Print out maximum values of viscosity used on this processor - output.write(" Hyper-viscosity values:\n"); - output.write(" Max mu_x = {:e}, Max_DC mu_x = {:e}\n", max(hyper_mu_x), - max(DC(hyper_mu_x))); - } - } + if (viscos_perp > 0.0) { + if (laplace_perp) { + ddt(U) += viscos_perp * Laplace_perp(U); // Perpendicular viscosity + } else { + ddt(U) += viscos_perp * Delp2(U); // Perpendicular viscosity + } + } - if (gyroviscous) { - - Field3D Pi; - Field2D Pi0; - Pi = 0.5 * P; - Pi0 = 0.5 * P0; - - Dperp2Phi0 = Field3D(Delp2(B0 * phi0)); - Dperp2Phi0.applyBoundary(); - mesh->communicate(Dperp2Phi0); - - Dperp2Phi = Delp2(B0 * phi); - Dperp2Phi.applyBoundary(); - mesh->communicate(Dperp2Phi); - - Dperp2Pi0 = Field3D(Delp2(Pi0)); - Dperp2Pi0.applyBoundary(); - mesh->communicate(Dperp2Pi0); - - Dperp2Pi = Delp2(Pi); - Dperp2Pi.applyBoundary(); - mesh->communicate(Dperp2Pi); - - bracketPhi0P = bracket(B0 * phi0, Pi, bm_exb); - bracketPhi0P.applyBoundary(); - mesh->communicate(bracketPhi0P); - - bracketPhiP0 = bracket(B0 * phi, Pi0, bm_exb); - bracketPhiP0.applyBoundary(); - mesh->communicate(bracketPhiP0); - - ddt(U) -= 0.5 * Upara2 * bracket(Pi, Dperp2Phi0, bm_exb) / B0; - ddt(U) -= 0.5 * Upara2 * bracket(Pi0, Dperp2Phi, bm_exb) / B0; - Field3D B0phi = B0 * phi; - mesh->communicate(B0phi); - Field3D B0phi0 = B0 * phi0; - mesh->communicate(B0phi0); - ddt(U) += 0.5 * Upara2 * bracket(B0phi, Dperp2Pi0, bm_exb) / B0; - ddt(U) += 0.5 * Upara2 * bracket(B0phi0, Dperp2Pi, bm_exb) / B0; - ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhi0P) / B0; - ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhiP0) / B0; - - if (nonlinear) { - Field3D B0phi = B0 * phi; - mesh->communicate(B0phi); - bracketPhiP = bracket(B0phi, Pi, bm_exb); - bracketPhiP.applyBoundary(); - mesh->communicate(bracketPhiP); - - ddt(U) -= 0.5 * Upara2 * bracket(Pi, Dperp2Phi, bm_exb) / B0; - ddt(U) += 0.5 * Upara2 * bracket(B0phi, Dperp2Pi, bm_exb) / B0; - ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhiP) / B0; - } - } + // Hyper-viscosity + if (hyperviscos > 0.0) { + // Calculate coefficient. - // left edge sink terms - if (sink_Ul > 0.0) { - ddt(U) -= sink_Ul * sink_tanhxl(P0, U, su_widthl, su_lengthl); // core sink - } + hyper_mu_x = hyperviscos * metric->g_11 * SQ(metric->dx) + * abs(metric->g11 * D2DX2(U)) / (abs(U) + 1e-3); + hyper_mu_x.applyBoundary("dirichlet"); // Set to zero on all boundaries - // right edge sink terms - if (sink_Ur > 0.0) { - ddt(U) -= sink_Ur * sink_tanhxr(P0, U, su_widthr, su_lengthr); // sol sink - } + ddt(U) += hyper_mu_x * metric->g11 * D2DX2(U); - //////////////////////////////////////////////////// - // Pressure equation + if (first_run) { // Print out maximum values of viscosity used on this processor + output.write(" Hyper-viscosity values:\n"); + output.write(" Max mu_x = {:e}, Max_DC mu_x = {:e}\n", max(hyper_mu_x), + max(DC(hyper_mu_x))); + } + } - ddt(P) = 0.0; - if (evolve_pressure) { - ddt(P) -= b0xGrad_dot_Grad(phi, P0); + if (gyroviscous) { + + Field3D Pi; + Field2D Pi0; + Pi = 0.5 * P; + Pi0 = 0.5 * P0; + + if (laplace_perp) { + Dperp2Phi0 = Field3D(Laplace_perp(B0 * phi0)); + Dperp2Phi0.applyBoundary(); + mesh->communicate(Dperp2Phi0); + + Dperp2Phi = Laplace_perp(B0 * phi); + Dperp2Phi.applyBoundary(); + mesh->communicate(Dperp2Phi); + + Dperp2Pi0 = Field3D(Laplace_perp(Pi0)); + Dperp2Pi0.applyBoundary(); + mesh->communicate(Dperp2Pi0); + + Dperp2Pi = Laplace_perp(Pi); + Dperp2Pi.applyBoundary(); + mesh->communicate(Dperp2Pi); + } else { + Dperp2Phi0 = Field3D(Delp2(B0 * phi0)); + Dperp2Phi0.applyBoundary(); + mesh->communicate(Dperp2Phi0); + + Dperp2Phi = Delp2(B0 * phi); + Dperp2Phi.applyBoundary(); + mesh->communicate(Dperp2Phi); + + Dperp2Pi0 = Field3D(Delp2(Pi0)); + Dperp2Pi0.applyBoundary(); + mesh->communicate(Dperp2Pi0); + + Dperp2Pi = Delp2(Pi); + Dperp2Pi.applyBoundary(); + mesh->communicate(Dperp2Pi); + } + + bracketPhi0P = bracket(B0 * phi0, Pi, bm_exb); + bracketPhi0P.applyBoundary(); + mesh->communicate(bracketPhi0P); + + bracketPhiP0 = bracket(B0 * phi, Pi0, bm_exb); + bracketPhiP0.applyBoundary(); + mesh->communicate(bracketPhiP0); + + ddt(U) -= 0.5 * Upara2 * bracket(Pi, Dperp2Phi0, bm_exb) / B0; + ddt(U) -= 0.5 * Upara2 * bracket(Pi0, Dperp2Phi, bm_exb) / B0; + Field3D B0phi = B0 * phi; + mesh->communicate(B0phi); + Field3D B0phi0 = B0 * phi0; + mesh->communicate(B0phi0); + ddt(U) += 0.5 * Upara2 * bracket(B0phi, Dperp2Pi0, bm_exb) / B0; + ddt(U) += 0.5 * Upara2 * bracket(B0phi0, Dperp2Pi, bm_exb) / B0; + if (laplace_perp) { + ddt(U) -= 0.5 * Upara2 * Laplace_perp(bracketPhi0P) / B0; + ddt(U) -= 0.5 * Upara2 * Laplace_perp(bracketPhiP0) / B0; + } else { + ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhi0P) / B0; + ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhiP0) / B0; + } + + if (nonlinear) { + Field3D B0phi = B0 * phi; + mesh->communicate(B0phi); + bracketPhiP = bracket(B0phi, Pi, bm_exb); + bracketPhiP.applyBoundary(); + mesh->communicate(bracketPhiP); + + ddt(U) -= 0.5 * Upara2 * bracket(Pi, Dperp2Phi, bm_exb) / B0; + ddt(U) += 0.5 * Upara2 * bracket(B0phi, Dperp2Pi, bm_exb) / B0; + if (laplace_perp) { + ddt(U) -= 0.5 * Upara2 * Laplace_perp(bracketPhiP) / B0; + } else { + ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhiP) / B0; + } + } + } - if (diamag_phi0) { // Equilibrium flow - ddt(P) -= b0xGrad_dot_Grad(phi0, P); - } + // left edge sink terms + if (sink_Ul > 0.0) { + ddt(U) -= sink_Ul * sink_tanhxl(P0, U, su_widthl, su_lengthl); // core sink + } - if (withflow) { // net flow - ddt(P) -= V_dot_Grad(V0net, P); - } + // right edge sink terms + if (sink_Ur > 0.0) { + ddt(U) -= sink_Ur * sink_tanhxr(P0, U, su_widthr, su_lengthr); // sol sink + } - if (nonlinear) { // Advection - ddt(P) -= bracket(phi, P, bm_exb) * B0; - } - } + //////////////////////////////////////////////////// + // Pressure equation - // Parallel diffusion terms + ddt(P) = 0.0; + if (evolve_pressure) { + ddt(P) -= b0xGrad_dot_Grad(phi, P0); - if (diffusion_par > 0.0) { // Parallel diffusion - ddt(P) += diffusion_par * Grad2_par2(P); - } + if (diamag_phi0) { // Equilibrium flow + ddt(P) -= b0xGrad_dot_Grad(phi0, P); + } - if (diffusion_p4 > 0.0) { // parallel hyper-viscous diffusion for pressure - tmpP2 = D2DY2(P); - mesh->communicate(tmpP2); - tmpP2.applyBoundary(); - ddt(P) = diffusion_p4 * D2DY2(tmpP2); - } + if (withflow) { // net flow + ddt(P) -= V_dot_Grad(V0net, P); + } - if (heating_P > 0.0) { // heating source terms - BoutReal pnorm = P0(0, 0); - ddt(P) += heating_P * source_expx2(P0, 2. * hp_width, 0.5 * hp_length) - * (Tbar / pnorm); // heat source - ddt(P) += (100. * source_tanhx(P0, hp_width, hp_length) + 0.01) * metric->g11() - * D2DX2(P) * (Tbar / Lbar / Lbar); // radial diffusion - } + if (nonlinear) { // Advection + ddt(P) -= bracket(phi, P, bm_exb) * B0; + } + } - if (sink_P > 0.0) { // sink terms - ddt(P) -= sink_P * sink_tanhxr(P0, P, sp_width, sp_length) * Tbar; // sink - } + // Parallel diffusion terms - //////////////////////////////////////////////////// - // Compressional effects + if (diffusion_par > 0.0) { // Parallel diffusion + ddt(P) += diffusion_par * Grad2_par2(P); + } - if (compress) { + if (diffusion_p4 > 0.0) { // parallel hyper-viscous diffusion for pressure + tmpP2 = D2DY2(P); + mesh->communicate(tmpP2); + tmpP2.applyBoundary(); + ddt(P) = diffusion_p4 * D2DY2(tmpP2); + } - ddt(P) -= beta * Div_par(Vpar, CELL_CENTRE); + if (heating_P > 0.0) { // heating source terms + BoutReal pnorm = P0(0, 0); + ddt(P) += heating_P * source_expx2(P0, 2. * hp_width, 0.5 * hp_length) + * (Tbar / pnorm); // heat source + ddt(P) += (100. * source_tanhx(P0, hp_width, hp_length) + 0.01) * metric->g11 + * D2DX2(P) * (Tbar / Lbar / Lbar); // radial diffusion + } - if (phi_curv) { - ddt(P) -= 2. * beta * b0xcv * Grad(phi); - } + if (sink_P > 0.0) { // sink terms + ddt(P) -= sink_P * sink_tanhxr(P0, P, sp_width, sp_length) * Tbar; // sink + } - // Vpar equation + //////////////////////////////////////////////////// + // Compressional effects - ddt(Vpar) = -0.5 * (Grad_par(P, loc) + Grad_par(P0, loc)); + if (compress) { - if (nonlinear) { - ddt(Vpar) -= bracket(interp_to(phi, loc), Vpar, bm_exb) * B0; // Advection - } - } + ddt(P) -= beta * Div_par(Vpar, CELL_CENTRE); - if (filter_z) { - // Filter out all except filter_z_mode + if (phi_curv) { + ddt(P) -= 2. * beta * b0xcv * Grad(phi); + } - if (evolve_jpar) { - ddt(Jpar) = filter(ddt(Jpar), filter_z_mode); - } else { - ddt(Psi) = filter(ddt(Psi), filter_z_mode); - } + // Vpar equation - ddt(U) = filter(ddt(U), filter_z_mode); - ddt(P) = filter(ddt(P), filter_z_mode); - } + ddt(Vpar) = -0.5 * (Grad_par(P, loc) + Grad_par(P0, loc)); - if (low_pass_z > 0) { - // Low-pass filter, keeping n up to low_pass_z - if (evolve_jpar) { - ddt(Jpar) = lowPass(ddt(Jpar), low_pass_z, zonal_field); - } else { - ddt(Psi) = lowPass(ddt(Psi), low_pass_z, zonal_field); - } - ddt(U) = lowPass(ddt(U), low_pass_z, zonal_flow); - ddt(P) = lowPass(ddt(P), low_pass_z, zonal_bkgd); - } + if (nonlinear) { + ddt(Vpar) -= bracket(interp_to(phi, loc), Vpar, bm_exb) * B0; // Advection + } + } - if (damp_width > 0) { - for (int i = 0; i < damp_width; i++) { - for (int j = 0; j < mesh->LocalNy; j++) { - for (int k = 0; k < mesh->LocalNz; k++) { - if (mesh->firstX()) { - ddt(U)(i, j, k) -= U(i, j, k) / damp_t_const; - } - if (mesh->lastX()) { - ddt(U)(mesh->LocalNx - 1 - i, j, k) -= - U(mesh->LocalNx - 1 - i, j, k) / damp_t_const; - } - } - } - } - } + if (filter_z) { + // Filter out all except filter_z_mode - first_run = false; + if (evolve_jpar) { + ddt(Jpar) = filter(ddt(Jpar), filter_z_mode); + } else { + ddt(Psi) = filter(ddt(Psi), filter_z_mode); + } - return 0; + ddt(U) = filter(ddt(U), filter_z_mode); + ddt(P) = filter(ddt(P), filter_z_mode); } - /******************************************************************************* - * Preconditioner - * - * o System state in variables (as in rhs function) - * o Values to be inverted in time derivatives - * - * o Return values should be in time derivatives - * - * enable by setting solver / use_precon = true in BOUT.inp - *******************************************************************************/ - - int precon(BoutReal UNUSED(t), BoutReal gamma, BoutReal UNUSED(delta)) { - // First matrix, applying L - mesh->communicate(ddt(Psi)); - Field3D Jrhs = Delp2(ddt(Psi)); - Jrhs.applyBoundary("neumann"); + if (low_pass_z > 0) { + // Low-pass filter, keeping n up to low_pass_z + if (evolve_jpar) { + ddt(Jpar) = lowPass(ddt(Jpar), low_pass_z, zonal_field); + } else { + ddt(Psi) = lowPass(ddt(Psi), low_pass_z, zonal_field); + } + ddt(U) = lowPass(ddt(U), low_pass_z, zonal_flow); + ddt(P) = lowPass(ddt(P), low_pass_z, zonal_bkgd); + } - if (jpar_bndry_width > 0) { - // Boundary in jpar + if (damp_width > 0) { + for (int i = 0; i < damp_width; i++) { + for (int j = 0; j < mesh->LocalNy; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { if (mesh->firstX()) { - for (int i = jpar_bndry_width; i >= 0; i--) { - for (int j = 0; j < mesh->LocalNy; j++) { - for (int k = 0; k < mesh->LocalNz; k++) { - Jrhs(i, j, k) = 0.5 * Jrhs(i + 1, j, k); - } - } - } + ddt(U)(i, j, k) -= U(i, j, k) / damp_t_const; } if (mesh->lastX()) { - for (int i = mesh->LocalNx - jpar_bndry_width - 1; i < mesh->LocalNx; i++) { - for (int j = 0; j < mesh->LocalNy; j++) { - for (int k = 0; k < mesh->LocalNz; k++) { - Jrhs(i, j, k) = 0.5 * Jrhs(i - 1, j, k); - } - } - } + ddt(U)(mesh->LocalNx - 1 - i, j, k) -= + U(mesh->LocalNx - 1 - i, j, k) / damp_t_const; } + } } + } + } - mesh->communicate(Jrhs, ddt(P)); - - Field3D U1 = ddt(U); - - auto B0 = tokamak_options.Bxy; - - U1 += (gamma * B0 * B0) * Grad_par(Jrhs, CELL_CENTRE) + (gamma * b0xcv) * Grad(P); - - // Second matrix, solving Alfven wave dynamics - static std::unique_ptr invU{nullptr}; - if (!invU) { - invU = InvertPar::create(); + first_run = false; + + return 0; + } + + /******************************************************************************* + * Preconditioner + * + * o System state in variables (as in rhs function) + * o Values to be inverted in time derivatives + * + * o Return values should be in time derivatives + * + * enable by setting solver / use_precon = true in BOUT.inp + *******************************************************************************/ + + int precon(BoutReal UNUSED(t), BoutReal gamma, BoutReal UNUSED(delta)) { + // First matrix, applying L + mesh->communicate(ddt(Psi)); + ddt(Psi).applyBoundary("neumann"); + Field3D Jrhs; + if (laplace_perp) { + Jrhs = Laplace_perp(ddt(Psi)); + } else { + Jrhs = Delp2(ddt(Psi)); + } + Jrhs.applyBoundary("neumann"); + + if (jpar_bndry_width > 0) { + // Boundary in jpar + if (mesh->firstX()) { + for (int i = jpar_bndry_width; i >= 0; i--) { + for (int j = 0; j < mesh->LocalNy; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + Jrhs(i, j, k) = 0.5 * Jrhs(i + 1, j, k); + } + } + } + } + if (mesh->lastX()) { + for (int i = mesh->LocalNx - jpar_bndry_width - 1; i < mesh->LocalNx; i++) { + for (int j = 0; j < mesh->LocalNy; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + Jrhs(i, j, k) = 0.5 * Jrhs(i - 1, j, k); + } + } } - - invU->setCoefA(1.); - invU->setCoefB(-SQ(gamma) * B0 * B0); - ddt(U) = invU->solve(U1); - ddt(U).applyBoundary(); - - // Third matrix, applying U - Field3D phi3 = phiSolver->solve(ddt(U)); - mesh->communicate(phi3); - phi3.applyBoundary("neumann"); - Field3D B0phi3 = B0 * phi3; - mesh->communicate(B0phi3); - ddt(Psi) = ddt(Psi) - gamma * Grad_par(B0phi3, loc) / B0; - ddt(Psi).applyBoundary(); - - return 0; + } } - /******************************************************************************* - * Jacobian-vector multiply - * - * Input - * System state is in (P, Psi, U) - * Vector v is in (F_P, F_Psi, F_U) - * Output - * Jacobian-vector multiplied Jv should be in (P, Psi, U) - * - * NOTE: EXPERIMENTAL - * enable by setting solver / use_jacobian = true in BOUT.inp - *******************************************************************************/ - - int jacobian(BoutReal UNUSED(t)) { - // Communicate - mesh->communicate(ddt(P), ddt(Psi), ddt(U)); - - phi = phiSolver->solve(ddt(U)); + mesh->communicate(Jrhs, ddt(P)); + ddt(P).applyBoundary("neumann"); - Jpar = Delp2(ddt(Psi)); + Field3D U1 = ddt(U); + U1 += (gamma * B0 * B0) * Grad_par(Jrhs, CELL_CENTRE) + (gamma * b0xcv) * Grad(ddt(P)); - mesh->communicate(phi, Jpar); - - Field3D JP = -b0xGrad_dot_Grad(phi, P0); - JP.setBoundary("P"); - JP.applyBoundary(); - - auto B0 = tokamak_options.Bxy; - - Field3D B0phi = B0 * phi; - mesh->communicate(B0phi); - Field3D JPsi = -Grad_par(B0phi, loc) / B0; - JPsi.setBoundary("Psi"); - JPsi.applyBoundary(); - - Field3D JU = b0xcv * Grad(ddt(P)) - SQ(B0) * Grad_par(Jpar, CELL_CENTRE) - + SQ(B0) * b0xGrad_dot_Grad(ddt(Psi), J0, CELL_CENTRE); - JU.setBoundary("U"); - JU.applyBoundary(); - - // Put result into time-derivatives - - ddt(P) = JP; - ddt(Psi) = JPsi; - ddt(U) = JU; - - return 0; + // Second matrix, solving Alfven wave dynamics + static std::unique_ptr invU{nullptr}; + if (!invU) { + invU = InvertPar::create(); } - /******************************************************************************* - * Preconditioner for when phi solved as a constraint - * Currently only possible with the IDA solver - * - * o System state in variables (as in rhs function) - * o Values to be inverted in F_vars - * - * o Return values should be in vars (overwriting system state) - *******************************************************************************/ - - int precon_phi(BoutReal UNUSED(t), BoutReal UNUSED(cj), BoutReal UNUSED(delta)) { - ddt(phi) = phiSolver->solve(C_phi - ddt(U)); - return 0; + invU->setCoefA(1.); + invU->setCoefB(-SQ(gamma) * B0 * B0); + ddt(U) = invU->solve(U1); + ddt(U).applyBoundary(); + + // Third matrix, applying U + Field3D phi3 = phiSolver->solve(ddt(U)); + mesh->communicate(phi3); + phi3.applyBoundary("neumann"); + Field3D B0phi3 = B0 * phi3; + mesh->communicate(B0phi3); + ddt(Psi) = ddt(Psi) - gamma * Grad_par(B0phi3, loc) / B0; + ddt(Psi).applyBoundary(); + + return 0; + } + + /******************************************************************************* + * Jacobian-vector multiply + * + * Input + * System state is in (P, Psi, U) + * Vector v is in (F_P, F_Psi, F_U) + * Output + * Jacobian-vector multiplied Jv should be in (P, Psi, U) + * + * NOTE: EXPERIMENTAL + * enable by setting solver / use_jacobian = true in BOUT.inp + *******************************************************************************/ + + int jacobian(BoutReal UNUSED(t)) { + // Communicate + mesh->communicate(ddt(P), ddt(Psi), ddt(U)); + + phi = phiSolver->solve(ddt(U)); + + if (laplace_perp) { + Jpar = Laplace_perp(ddt(Psi)); + } else { + Jpar = Delp2(ddt(Psi)); } + mesh->communicate(phi, Jpar); + + Field3D JP = -b0xGrad_dot_Grad(phi, P0); + JP.setBoundary("P"); + JP.applyBoundary(); + Field3D B0phi = B0 * phi; + mesh->communicate(B0phi); + Field3D JPsi = -Grad_par(B0phi, loc) / B0; + JPsi.setBoundary("Psi"); + JPsi.applyBoundary(); + + Field3D JU = b0xcv * Grad(ddt(P)) - SQ(B0) * Grad_par(Jpar, CELL_CENTRE) + + SQ(B0) * b0xGrad_dot_Grad(ddt(Psi), J0, CELL_CENTRE); + JU.setBoundary("U"); + JU.applyBoundary(); + + // Put result into time-derivatives + + ddt(P) = JP; + ddt(Psi) = JPsi; + ddt(U) = JU; + + return 0; + } + + /******************************************************************************* + * Preconditioner for when phi solved as a constraint + * Currently only possible with the IDA solver + * + * o System state in variables (as in rhs function) + * o Values to be inverted in F_vars + * + * o Return values should be in vars (overwriting system state) + *******************************************************************************/ + + int precon_phi(BoutReal UNUSED(t), BoutReal UNUSED(cj), BoutReal UNUSED(delta)) { + ddt(phi) = phiSolver->solve(C_phi - ddt(U)); + return 0; + } }; BOUTMAIN(ELMpb); diff --git a/examples/elm-pb/relaxing/BOUT.inp b/examples/elm-pb/relaxing/BOUT.inp new file mode 100644 index 0000000000..d9489611e6 --- /dev/null +++ b/examples/elm-pb/relaxing/BOUT.inp @@ -0,0 +1,299 @@ +# settings file for BOUT++ +# High-Beta reduced MHD case + +################################################## +# Global settings used by the core code + +nout = 1000 # number of time-steps +timestep = 1 # time between outputs + +zperiod = 5 # Fraction of a torus to simulate +MZ = 64 # Number of points in Z + +grid = "cbm18_dens8.grid_nx68ny64.nc" #"cbm18_dens3_0.5BS_516nx64ny.grid.nc" + +[mesh] + +staggergrids = false # Use staggered grids + +[mesh:paralleltransform] +#type = shifted # Use shifted metric method +type = shiftedinterp + +################################################## +# derivative methods + +[mesh:ddx] +first = C4 # order of first x derivatives +second = C4 # order of second x derivatives +upwind = W3 # order of upwinding method W3 = Weno3 + +[mesh:ddy] +first = C4 +second = C4 +upwind = W3 +flux = C2 + +[mesh:ddz] +first = C4 # Z derivatives can be done using FFT +second = C4 +upwind = W3 + +################################################## +# FFTs + +[fft] + +fft_measurement_flag = measure # If using FFTW, perform tests to determine fastest method + +[output] + +#type = adios +#shiftoutput = true # Put the output into field-aligned coordinates + +[restart_files] + +#type = adios + +[laplace] + +#type = hypre3d +#rtol = 1.e-9 +#atol = 1.e-14 +#pctype = sor # Preconditioner + +#atol = 1e-12 +#rtol = 1e-08 + + +################################################## +# Solver settings + +[solver] + +# mudq, mldq, mukeep, mlkeep preconditioner options +atol = 1.0e-8 # absolute tolerance +rtol = 1.0e-5 # relative tolerance + +use_precon = false # Use preconditioner: User-supplied or BBD + +mxstep = 50000 # Number of internal steps between outputs + +################################################## +# settings for high-beta reduced MHD + +[highbeta] + +density = 1.0e19 # number density of deuterium [m^-3] + # used to produce output normalisations +constn0 = true +n0_fake_prof = false + +sheath_boundaries = false + +evolve_jpar = false # If true, evolve J raher than Psi + +evolve_pressure = true # If false, switch off all pressure evolution + +phi_constraint = false # Solve phi as a constraint (DAE system, needs IDA) + +## Effects to include/exclude + +include_jpar0 = true # determines whether to include jpar0 terms +include_curvature = true # include curvature drive term? + +compress = false # set compressible (evolve Vpar) +nonlinear = true # include non-linear terms? + +diamag = true # Include diamagnetic effects? +diamag_grad_t = false # Include Grad_par(Te) term in Psi equation +diamag_phi0 = true # Balance ExB against Vd for stationary equilibrium + +split_n0 = false + +laplace_perp = true + +################################################## +# BRACKET_METHOD flags: +# 0:BRACKET_STD; derivative methods will be determined +# by the choices C or W in this input file +# 1:BRACKET_SIMPLE; 2:BRACKET_ARAKAWA; 3:BRACKET_CTU. + +bm_exb_flag = 0 +bm_mag_flag = 0 +################################################################## +withflow = false # With flow or not +D_0 = 130000 # differential potential +D_s = 20 # shear parameter +K_H_term = false # Contain K-H term +sign = -1 # flow direction +x0 = 0.855 # peak location +D_min = 3000 # constant +################################################################## + +eHall = false # Include electron pressure effects in Ohm's law? +AA = 2.0 # ion mass in units of proton mass +Zeff = 1.0 +#Zi = 1.0 + +noshear = false # zero all shear + +relax_j_vac = false # Relax to zero-current in the vacuum +relax_j_tconst = 1e-2 # Time constant for vacuum relaxation + +## Toroidal filtering +filter_z = false # remove all except one mode +filter_z_mode = 1 # Specify which harmonic to keep (1 = fundamental) +low_pass_z = 16 # Keep up to and including this harmonic (-1 = keep all) +zonal_flow = true # keep zonal component of vorticity? +zonal_field = false # keep zonal component of Psi? +zonal_bkgd = true # keep zonal component of P? + +## Jpar smoothing +smooth_j_x = true # Filter Jpar in the X direction + +## Magnetic perturbations +include_rmp = false # Read RMP data from grid file + +simple_rmp = false # Enable/disable a simple model of RMP +rmp_n = 3 # Toroidal mode number +rmp_m = 6 # Poloidal mode number +rmp_factor = 1.e-4 # Amplitude of Apar [Tm] +rmp_ramp = 1.e-4 # Timescale [s] of ramp +rmp_polwid = -1.0 # Width of Gaussian factor (< 0 = No Gaussian) +rmp_polpeak = 0.5 # Y location of maximum (fraction) + +## Vacuum region control + +vacuum_pressure = 0.02 # the pressure below which it is considered vacuum + # fraction of peak pressure +vacuum_trans = 0.01 # transition width (fraction of P) + +## Resistivity and Hyper-resistivity + +vac_lund = 1.0e8 # Lundquist number in vacuum (negative -> infinity) +core_lund = 1.0e8 # Lundquist number in core (negative -> infinity) +hyperresist = 1.e-4 # Hyper-resistivity coefficient (like 1 / Lundquist number) + +## Inner boundary damping + +damp_width = -1 # Width of damping region (grid cells) +damp_t_const = 1e-2 # Damping time constant + +## Parallel pressure diffusion + +diffusion_par = 1.0e-2 # Parallel pressure diffusion (< 0 = none) +diffusion_p4 = -1e-05 # parallel hyper-viscous diffusion for pressure (< 0 = none) +diffusion_u4 = -1e-05 # parallel hyper-viscous diffusion for vorticity (< 0 = none) +diffusion_a4 = -1e-05 # parallel hyper-viscous diffusion for vector potential (< 0 = none) + +## heat source in pressure in watts + +heating_P = -1 # heat power in watts (< 0 = none) +hp_width = 0.1 # heat width, in percentage of nx (< 0 = none) +hp_length = 0.3 # heat length in percentage of nx (< 0 = none) + +## sink rate in pressure + +sink_P = -1 # sink rate in pressure (< 0 = none) +sp_width = 0.04 # sink width, in percentage of nx (< 0 = none) +sp_length = 0.15 # sink length in percentage of nx (< 0 = none) + + +## left edge sink rate in vorticity +sink_Ul = 10.0 # left edge sink rate in vorticity (< 0 = none) +su_widthl = 0.06 # left edge sink width, in percentage of nx (< 0 = none) +su_lengthl = 0.1 # left edge sink length in percentage of nx (< 0 = none) + +## right edge sink rate in vorticity +sink_Ur = 10.0 # right edge sink rate in vorticity (< 0 = none) +su_widthr = 0.06 # right edge sink width, in percentage of nx (< 0 = none) +su_lengthr = 0.1 # right edge sink length in percentage of nx (< 0 = none) + +## Viscosity and Hyper-viscosity + +viscos_par = 0.1 # Parallel viscosity (< 0 = none) +viscos_perp = 1.0e-7 # Perpendicular viscosity (< 0 = none) +hyperviscos = -1.0 # Radial hyper viscosity + +## Compressional terms (only when compress = true) +phi_curv = true # Include curvature*Grad(phi) in P equation +# gamma = 1.6666 + +phi_boundary_relax = true # Relax phi at radial boundaries towards zero gradient +phi_core_averagey = false +phi_boundary_timescale = 1.0e-6 # In seconds + +[phiSolver] +# INVERT_DC_GRAD = 1 // Zero-gradient for DC (constant in Z) component. Default is zero value +# INVERT_AC_GRAD = 2 // Zero-gradient for AC (non-constant in Z) component. Default is zero value +# INVERT_AC_LAP = 4 // Use zero-laplacian (decaying solution) to AC component +# INVERT_DC_LAP = 64 // Use zero-laplacian solution for DC component +#type = hypre3d +#rtol = 1.e-9 +#atol = 1.e-14 +#inner_boundary_flags = 0 # 0 for Dirichlet; 2 for Neumann +#outer_boundary_flags = 0 # 0 for Dirichlet; 2 for Neumann +#hypre_print_level = 1 # print information on the matrices, solver settings, and iterations. + +[aparSolver] +#type = hypre3d +#rtol = 1.e-9 +#atol = 1.e-14 +#inner_boundary_flags = 4 # INVERT_AC_LAP +#outer_boundary_flags = 1 + 4 # INVERT_DC_GRAD + INVERT_AC_LAP + +################################################## +# settings for individual variables +# The section "All" defines default settings for all variables +# These can be overridden for individual variables in +# a section of that name. + +[all] +scale = 0.0 # default size of initial perturbations + +# boundary conditions +# ------------------- +# dirichlet - Zero value +# neumann - Zero gradient +# zerolaplace - Laplacian = 0, decaying solution +# constlaplace - Laplacian = const, decaying solution +# +# relax( ) - Make boundary condition relaxing + +bndry_all = dirichlet_o2 # Default to zero-value + +[U] # vorticity +scale = 1e-05 +function = ballooning(gauss(x-0.5,0.1)*gauss(y-pi,0.6*pi)*sin(z),3) + +[P] # pressure +bndry_core = dirichlet +#scale = 1.0e-5 + +[Psi] # Vector potential + +# zero laplacian +bndry_xin = zerolaplace +bndry_xout = zerolaplace + +[Psi_loc] # for staggering + +bndry_xin = zerolaplace +bndry_xout = zerolaplace + +bndry_yup = free_o3 +bndry_ydown = free_o3 + +[J] # parallel current + +# Zero gradient in the core +bndry_core = neumann + +[phi] + +#bndry_core = neumann +bndry_xin = none +bndry_xout = none +#bndry_target = neumann + diff --git a/examples/laplacexy/laplace_perp/square/BOUT.inp b/examples/laplacexy/laplace_perp/square/BOUT.inp index 95ecb7119f..38fb175d90 100644 --- a/examples/laplacexy/laplace_perp/square/BOUT.inp +++ b/examples/laplacexy/laplace_perp/square/BOUT.inp @@ -1,5 +1,5 @@ -input = sin(pi*x - y) +input_field = sin(pi*x - y) non_uniform = true # Include corrections to second derivatives diff --git a/examples/laplacexy/laplace_perp/test.cxx b/examples/laplacexy/laplace_perp/test.cxx index 87d4b5ecb0..3c6c1a979d 100644 --- a/examples/laplacexy/laplace_perp/test.cxx +++ b/examples/laplacexy/laplace_perp/test.cxx @@ -1,4 +1,6 @@ #include +#include +#include #include #include @@ -11,24 +13,37 @@ int main(int argc, char** argv) { BoutInitialise(argc, argv); /////////////////////////////////////// - bool calc_metric; - calc_metric = Options::root()["calc_metric"].withDefault(false); + const bool calc_metric = Options::root()["calc_metric"].withDefault(false); if (calc_metric) { auto tokamak_options = bout::TokamakOptions(*mesh); set_tokamak_coordinates_on_mesh(tokamak_options, *mesh, 1.0, 1.0); } + // Read metric tensor + Field2D Rxy; + Field2D Btxy; + Field2D Bpxy; + Field2D B0; + Field2D hthe; + Field2D I; + mesh->get(Rxy, "Rxy"); // m + mesh->get(Btxy, "Btxy"); // T + mesh->get(Bpxy, "Bpxy"); // T + mesh->get(B0, "Bxy"); // T + mesh->get(hthe, "hthe"); // m + mesh->get(I, "sinty"); // m^-2 T^-1 Coordinates* coord = mesh->getCoordinates(); /////////////////////////////////////// // Read an analytic input - Field2D input = FieldFactory::get()->create2D("input", Options::getRoot(), mesh); + const Field2D input = + FieldFactory::get()->create2D("input_field", Options::getRoot(), mesh); // Create a LaplaceXY solver - LaplaceXY* laplacexy = new LaplaceXY(mesh); + LaplaceXY laplacexy{mesh}; // Solve, using 0.0 as starting guess - Field2D solved = laplacexy->solve(input, 0.0); + Field2D solved = laplacexy.solve(input, 0.0); // Need to communicate guard cells mesh->communicate(solved); diff --git a/examples/laplacexy/laplace_perp/torus/BOUT.inp b/examples/laplacexy/laplace_perp/torus/BOUT.inp index 365294174f..4109b44bab 100644 --- a/examples/laplacexy/laplace_perp/torus/BOUT.inp +++ b/examples/laplacexy/laplace_perp/torus/BOUT.inp @@ -1,5 +1,5 @@ -input = sin(pi*x - y) +input_field = sin(pi*x - y) calc_metric = true # Read Rxy, Bpxy etc and calculate metric non_uniform = true # Include corrections to second derivatives diff --git a/include/bout/adios_object.hxx b/include/bout/adios_object.hxx index b14316f1ba..12da8ce4c7 100755 --- a/include/bout/adios_object.hxx +++ b/include/bout/adios_object.hxx @@ -57,11 +57,17 @@ public: } template - adios2::Variable GetArrayVariable(const std::string& varname, adios2::Dims& shape) { + adios2::Variable GetArrayVariable(const std::string& varname, adios2::Dims& shape, + const std::vector& dimNames, + int rank) { adios2::Variable v = io.InquireVariable(varname); if (!v) { adios2::Dims start(shape.size()); v = io.DefineVariable(varname, shape, start, shape); + if (!rank && dimNames.size()) { + io.DefineAttribute("__xarray_dimensions__", dimNames.data(), + dimNames.size(), varname, "/", true); + } } else { v.SetShape(shape); } diff --git a/include/bout/bout.hxx b/include/bout/bout.hxx index 8588fe56d9..a08696ec66 100644 --- a/include/bout/bout.hxx +++ b/include/bout/bout.hxx @@ -34,6 +34,8 @@ #ifndef BOUT_H #define BOUT_H +#include // std::filesystem (C++17) + // IWYU pragma: begin_keep, begin_export #include "bout/build_defines.hxx" @@ -106,10 +108,10 @@ void setupGetText(); struct CommandLineArgs { int verbosity{4}; bool color_output{false}; - std::string data_dir{"data"}; ///< Directory for data input/output - std::string opt_file{"BOUT.inp"}; ///< Filename for the options file - std::string set_file{"BOUT.settings"}; ///< Filename for the options file - std::string log_file{"BOUT.log"}; ///< File name for the log file + std::filesystem::path data_dir{"data"}; ///< Directory for data input/output + std::filesystem::path opt_file{"BOUT.inp"}; ///< Filename for the options file + std::filesystem::path set_file{"BOUT.settings"}; ///< Filename for the options file + std::filesystem::path log_file{"BOUT.log"}; ///< File name for the log file /// The original set of command line arguments std::vector original_argv; /// The "canonicalised" command line arguments, with single-letter diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index 4ae3ac0740..cf29f928ee 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -9,9 +9,9 @@ * * ************************************************************************** - * Copyright 2014 B.D.Dudson + * Copyright 2014-2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -269,9 +269,10 @@ public: transform = std::move(pt); } + bool hasParallelTransform() const { return transform != nullptr; } /// Return the parallel transform - ParallelTransform& getParallelTransform() { - ASSERT1(transform != nullptr); + ParallelTransform& getParallelTransform() const { + ASSERT1(hasParallelTransform()); return *transform; } diff --git a/include/bout/field.hxx b/include/bout/field.hxx index 188d529ef0..61edff6723 100644 --- a/include/bout/field.hxx +++ b/include/bout/field.hxx @@ -79,6 +79,8 @@ public: std::string name; + bool isFci() const; + #if CHECK > 0 // Routines to test guard/boundary cells set diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx index a4e01d4572..8160876d42 100644 --- a/include/bout/mesh.hxx +++ b/include/bout/mesh.hxx @@ -17,9 +17,9 @@ * * Incorporates code from topology.cpp and Communicator * ************************************************************************** - * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2010-2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -650,9 +650,22 @@ public: return inserted.first->second; } + std::shared_ptr + getCoordinatesConst(const CELL_LOC location = CELL_CENTRE) const { + ASSERT1(location != CELL_DEFAULT); + ASSERT1(location != CELL_VSHIFT); + + auto found = coords_map.find(location); + if (found != coords_map.end()) { + // True branch most common, returns immediately + return found->second; + } + throw BoutException("Coordinates not yet set. Use non-const version!"); + } + /// Returns the non-CELL_CENTRE location /// allowed as a staggered location - CELL_LOC getAllowedStaggerLoc(DIRECTION direction) const { + static CELL_LOC getAllowedStaggerLoc(DIRECTION direction) { AUTO_TRACE(); switch (direction) { case (DIRECTION::X): @@ -841,6 +854,16 @@ public: ASSERT1(RegionID.has_value()); return region3D[RegionID.value()]; } + bool isFci() const { + const auto coords = this->getCoordinatesConst(); + if (coords == nullptr) { + return false; + } + if (not coords->hasParallelTransform()) { + return false; + } + return not coords->getParallelTransform().canToFromFieldAligned(); + } private: /// Allocates default Coordinates objects diff --git a/include/bout/petsc_interface.hxx b/include/bout/petsc_interface.hxx index 1378419872..90dd188d8a 100644 --- a/include/bout/petsc_interface.hxx +++ b/include/bout/petsc_interface.hxx @@ -566,6 +566,26 @@ PetscVector operator*(const PetscMatrix& mat, const PetscVector& vec) { return PetscVector(vec, result); } +// Compatibility wrappers +// For < 3.24 +#if PETSC_VERSION_GE(3, 24, 0) \ + || (PETSC_VERSION_GE(3, 23, 0) && PETSC_VERSION_RELEASE == 0) +namespace bout { +template +constexpr auto cast_MatFDColoringFn(T func) { + return func; +} +} // namespace bout +#else +using MatFDColoringFn = PetscErrorCode (*)(); +namespace bout { +template +constexpr auto cast_MatFDColoringFn(T func) { + return reinterpret_cast(func); // NOLINT(*-reinterpret-cast) +} +} // namespace bout +#endif + #endif // BOUT_HAS_PETSC #endif // BOUT_PETSC_INTERFACE_H diff --git a/include/bout/tokamak_coordinates.hxx b/include/bout/tokamak_coordinates.hxx index f13e718787..0f3f2ea65f 100644 --- a/include/bout/tokamak_coordinates.hxx +++ b/include/bout/tokamak_coordinates.hxx @@ -11,18 +11,31 @@ using FieldMetric = MetricTensor::FieldMetric; namespace bout { + struct Coordinates3D { + + Field3D x; + Field3D y; + Field3D z; + + Coordinates3D(Field3D& x, Field3D& y, Field3D& z) : x(x), y(y), z(z) {} + }; + struct TokamakOptions { Field2D Rxy; + Field2D Zxy; Field2D Bpxy; Field2D Btxy; Field2D Bxy; Field2D hthe; FieldMetric I; FieldMetric dx; + std::vector toroidal_angles; - TokamakOptions(Mesh &mesh); + TokamakOptions(Mesh& mesh); void normalise(BoutReal Lbar, BoutReal Bbar, BoutReal ShearFactor); + + Coordinates3D cylindrical_coordinates_to_cartesian(); }; BoutReal get_sign_of_bp(const Field2D &Bpxy); diff --git a/manual/RELEASE_HOWTO.md b/manual/RELEASE_HOWTO.md index 4103e6f197..d8ced2afb5 100644 --- a/manual/RELEASE_HOWTO.md +++ b/manual/RELEASE_HOWTO.md @@ -51,6 +51,7 @@ Before merging PR: - [ ] [`manual/doxygen/Doxyfile_readthedocs`][Doxyfile_readthedocs]: `PROJECT_NUMBER` - [ ] [`manual/doxygen/Doxyfile`][Doxyfile]: `PROJECT_NUMBER` - [ ] [`CMakeLists.txt`][CMakeLists]: `_bout_previous_version`, `_bout_next_version` +- [ ] Update what version of PETSc and SUNDIALS we support (upper bound) After PR is merged: diff --git a/manual/sphinx/user_docs/advanced_install.rst b/manual/sphinx/user_docs/advanced_install.rst index 048a26a6e3..f8e594be07 100644 --- a/manual/sphinx/user_docs/advanced_install.rst +++ b/manual/sphinx/user_docs/advanced_install.rst @@ -355,7 +355,7 @@ BOUT++ can use PETSc https://www.mcs.anl.gov/petsc/ for time-integration and for solving elliptic problems, such as inverting Poisson and Helmholtz equations. -Currently, BOUT++ supports PETSc versions 3.7 - 3.19. More recent versions may +Currently, BOUT++ supports PETSc versions 3.7 - 3.23. More recent versions may well work, but the PETSc API does sometimes change in backward-incompatible ways, so this is not guaranteed. To install PETSc version 3.19, use the following steps:: diff --git a/manual/sphinx/user_docs/time_integration.rst b/manual/sphinx/user_docs/time_integration.rst index 7d0b56abcc..c1526c14e2 100644 --- a/manual/sphinx/user_docs/time_integration.rst +++ b/manual/sphinx/user_docs/time_integration.rst @@ -165,6 +165,9 @@ nonlinear solvers: `CVodeSetEpsLin `_. +The linear solver type can be set using the ``linear_solver`` option. +Valid choices include ``gmres`` (the default), ``fgmres``, ``tfqmr``, ``bcgs``. + IMEX-BDF2 --------- @@ -393,7 +396,9 @@ iterations within a given range. +---------------------------+---------------+----------------------------------------------------+ | predictor | true | Use linear predictor? | +---------------------------+---------------+----------------------------------------------------+ -| matrix_free | false | Use matrix free Jacobian-vector product? | +| matrix_free | false | Matrix-free preconditioning? | ++---------------------------+---------------+----------------------------------------------------+ +| matrix_free_operator | false | Use matrix free Jacobian-vector product? | +---------------------------+---------------+----------------------------------------------------+ | use_coloring | true | If ``matrix_free=false``, use coloring to speed up | | | | calculation of the Jacobian elements. | @@ -402,11 +407,18 @@ iterations within a given range. +---------------------------+---------------+----------------------------------------------------+ | kspsetinitialguessnonzero | false | If true, Use previous solution as KSP initial | +---------------------------+---------------+----------------------------------------------------+ -| use_precon | false | Use user-supplied preconditioner? | +| use_precon | false | If ``matrix_free=true``, use user-supplied | +| | | preconditioner? | | | | If false, the default PETSc preconditioner is used | +---------------------------+---------------+----------------------------------------------------+ | diagnose | false | Print diagnostic information every iteration | +---------------------------+---------------+----------------------------------------------------+ +| stencil:cross | 0 | If ``matrix_free=false`` and ``use_coloring=true`` | +| stencil:square | 0 | Set the size and shape of the Jacobian coloring | +| stencil:taxi | 2 | stencil. | ++---------------------------+---------------+----------------------------------------------------+ +| force_symmetric_coloring | false | Ensure that the Jacobian coloring is symmetric | ++---------------------------+---------------+----------------------------------------------------+ The predictor is linear extrapolation from the last two timesteps. It seems to be effective, but can be disabled by setting ``predictor = false``. @@ -444,6 +456,51 @@ Preconditioner types: Enable with command-line args ``-pc_type hypre -pc_hypre_type euclid -pc_hypre_euclid_levels k`` where ``k`` is the level (1-8 typically). +Jacobian coloring stencil +~~~~~~~~~~~~~~~~~~~~~~~~~ + +The stencil used to create the Jacobian colouring can be varied, +depending on which numerical operators are in use. + + +``solver:stencil:cross = N`` +e.g. for N == 2 + +.. code-block:: bash + + * + * + * * x * * + * + * + + +``solver:stencil:square = N`` +e.g. for N == 2 + +.. code-block:: bash + + * * * * * + * * * * * + * * x * * + * * * * * + * * * * * + +``solver:stencil:taxi = N`` +e.g. for N == 2 + +.. code-block:: bash + + * + * * * + * * x * * + * * * + * + +Setting ``solver:force_symmetric_coloring = true``, will make sure +that the jacobian colouring matrix is symmetric. This will often +include a few extra non-zeros that the stencil will miss otherwise + ODE integration --------------- diff --git a/src/bout++.cxx b/src/bout++.cxx index 7f23cf5f91..6de1c5713b 100644 --- a/src/bout++.cxx +++ b/src/bout++.cxx @@ -4,7 +4,7 @@ * Adapted from the BOUT code by B.Dudson, University of York, Oct 2007 * ************************************************************************** - * Copyright 2010-2023 BOUT++ contributors + * Copyright 2010-2025 BOUT++ contributors * * Contact Ben Dudson, dudson2@llnl.gov * @@ -84,7 +84,7 @@ const char DEFAULT_DIR[] = "data"; // Define S_ISDIR if not defined by system headers (that is, MSVC) // Taken from https://github.com/curl/curl/blob/e59540139a398dc70fde6aec487b19c5085105af/lib/curl_setup.h#L748-L751 #if !defined(S_ISDIR) && defined(S_IFMT) && defined(S_IFDIR) -#define S_ISDIR(m) (((m)&S_IFMT) == S_IFDIR) +#define S_ISDIR(m) (((m) & S_IFMT) == S_IFDIR) #endif #ifdef _MSC_VER @@ -180,7 +180,7 @@ int BoutInitialise(int& argc, char**& argv) { // `optionfile` here, but we'd need to call parseCommandLine // _first_ in order to do that and set the source, etc., but we // need to call that _second_ in order to override the input file - reader->read(Options::getRoot(), "{}/{}", args.data_dir, args.opt_file); + reader->read(Options::getRoot(), "{}", (args.data_dir / args.opt_file).string()); // Get options override from command-line reader->parseCommandLine(Options::getRoot(), args.argv); diff --git a/src/field/field.cxx b/src/field/field.cxx index e48a8f3ef7..797df6c405 100644 --- a/src/field/field.cxx +++ b/src/field/field.cxx @@ -39,3 +39,14 @@ int Field::getNx() const { return getMesh()->LocalNx; } int Field::getNy() const { return getMesh()->LocalNy; } int Field::getNz() const { return getMesh()->LocalNz; } + +bool Field::isFci() const { + const auto* coords = this->getCoordinates(); + if (coords == nullptr) { + return false; + } + if (not coords->hasParallelTransform()) { + return false; + } + return not coords->getParallelTransform().canToFromFieldAligned(); +} diff --git a/src/mesh/fv_ops.cxx b/src/mesh/fv_ops.cxx index 96a95f4cd9..6ee782370d 100644 --- a/src/mesh/fv_ops.cxx +++ b/src/mesh/fv_ops.cxx @@ -444,8 +444,8 @@ void communicateFluxes(Field3D& f) { comm_handle xin, xout; // Cache results to silence spurious compiler warning about xin, // xout possibly being uninitialised when used - bool const not_first = !mesh->firstX(); - bool const not_last = !mesh->lastX(); + const bool not_first = mesh->periodicX || !mesh->firstX(); + const bool not_last = mesh->periodicX || !mesh->lastX(); if (not_first) { xin = mesh->irecvXIn(f(0, 0), size, 0); } diff --git a/src/mesh/impls/bout/boutmesh.cxx b/src/mesh/impls/bout/boutmesh.cxx index faa02a4c33..82ab004f4d 100644 --- a/src/mesh/impls/bout/boutmesh.cxx +++ b/src/mesh/impls/bout/boutmesh.cxx @@ -2,19 +2,10 @@ * Implementation of the Mesh class, handling input files compatible with * BOUT / BOUT-06. * - * Changelog - * --------- - * - * 2015-01 Ben Dudson - * * - * - * 2010-05 Ben Dudson - * * Initial version, adapted from grid.cpp and topology.cpp - * ************************************************************************** - * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2010-2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -1530,42 +1521,66 @@ bool BoutMesh::firstX() const { return PE_XIND == 0; } bool BoutMesh::lastX() const { return PE_XIND == NXPE - 1; } int BoutMesh::sendXOut(BoutReal* buffer, int size, int tag) { + Timer timer("comms"); + + int proc {-1}; if (PE_XIND == NXPE - 1) { - return 1; + if (periodicX) { + // Wrap around to first processor in X + proc = PROC_NUM(0, PE_YIND); + } else { + return 1; + } + } else { + proc = PROC_NUM(PE_XIND + 1, PE_YIND); } - Timer timer("comms"); - - mpi->MPI_Send(buffer, size, PVEC_REAL_MPI_TYPE, PROC_NUM(PE_XIND + 1, PE_YIND), tag, + mpi->MPI_Send(buffer, size, PVEC_REAL_MPI_TYPE, proc, tag, BoutComm::get()); return 0; } int BoutMesh::sendXIn(BoutReal* buffer, int size, int tag) { + Timer timer("comms"); + + int proc {-1}; if (PE_XIND == 0) { - return 1; + if (periodicX) { + // Wrap around to last processor in X + proc = PROC_NUM(NXPE - 1, PE_YIND); + } else { + return 1; + } + } else { + proc = PROC_NUM(PE_XIND - 1, PE_YIND); } - Timer timer("comms"); - - mpi->MPI_Send(buffer, size, PVEC_REAL_MPI_TYPE, PROC_NUM(PE_XIND - 1, PE_YIND), tag, + mpi->MPI_Send(buffer, size, PVEC_REAL_MPI_TYPE, proc, tag, BoutComm::get()); return 0; } comm_handle BoutMesh::irecvXOut(BoutReal* buffer, int size, int tag) { + Timer timer("comms"); + + int proc {-1}; if (PE_XIND == NXPE - 1) { - return nullptr; + if (periodicX) { + // Wrap around to first processor in X + proc = PROC_NUM(0, PE_YIND); + } else { + return nullptr; + } + } else { + proc = PROC_NUM(PE_XIND + 1, PE_YIND); } - Timer timer("comms"); - // Get a communications handle. Not fussy about size of arrays CommHandle* ch = get_handle(0, 0); - mpi->MPI_Irecv(buffer, size, PVEC_REAL_MPI_TYPE, PROC_NUM(PE_XIND + 1, PE_YIND), tag, + mpi->MPI_Irecv(buffer, size, PVEC_REAL_MPI_TYPE, proc, tag, BoutComm::get(), ch->request); ch->in_progress = true; @@ -1574,16 +1589,24 @@ comm_handle BoutMesh::irecvXOut(BoutReal* buffer, int size, int tag) { } comm_handle BoutMesh::irecvXIn(BoutReal* buffer, int size, int tag) { + Timer timer("comms"); + + int proc {-1}; if (PE_XIND == 0) { - return nullptr; + if (periodicX) { + // Wrap around to last processor in X + proc = PROC_NUM(NXPE - 1, PE_YIND); + } else { + return nullptr; + } + } else { + proc = PROC_NUM(PE_XIND - 1, PE_YIND); } - Timer timer("comms"); - // Get a communications handle. Not fussy about size of arrays CommHandle* ch = get_handle(0, 0); - mpi->MPI_Irecv(buffer, size, PVEC_REAL_MPI_TYPE, PROC_NUM(PE_XIND - 1, PE_YIND), tag, + mpi->MPI_Irecv(buffer, size, PVEC_REAL_MPI_TYPE, proc, tag, BoutComm::get(), ch->request); ch->in_progress = true; diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx index acd4954f97..659df8600d 100644 --- a/src/mesh/parallel/fci.cxx +++ b/src/mesh/parallel/fci.cxx @@ -346,3 +346,9 @@ void FCITransform::integrateParallelSlices(Field3D& f) { f.ynext(map.offset) = map.integrate(f); } } + +void FCITransform::outputVars(Options& output_options) { + // Real-space coordinates of grid points + output_options["R"].force(R, "FCI"); + output_options["Z"].force(Z, "FCI"); +} diff --git a/src/mesh/parallel/fci.hxx b/src/mesh/parallel/fci.hxx index 3ec3321a6a..1a02f558e1 100644 --- a/src/mesh/parallel/fci.hxx +++ b/src/mesh/parallel/fci.hxx @@ -73,11 +73,15 @@ public: FCITransform() = delete; FCITransform(Mesh& mesh, const Coordinates::FieldMetric& dy, bool zperiodic = true, Options* opt = nullptr) - : ParallelTransform(mesh, opt) { + : ParallelTransform(mesh, opt), R{&mesh}, Z{&mesh} { // check the coordinate system used for the grid data source FCITransform::checkInputGrid(); + // Real-space coordinates of grid cells + mesh.get(R, "R", 0.0, false); + mesh.get(Z, "Z", 0.0, false); + auto forward_boundary_xin = std::make_shared("FCI_forward", BNDRY_PAR_FWD_XIN, +1, &mesh); auto backward_boundary_xin = std::make_shared( @@ -142,6 +146,10 @@ public: bool canToFromFieldAligned() const override { return false; } + /// Save mesh variables to output + /// If R and Z(x,y,z) coordinates are in the input then these are saved to output. + void outputVars(Options& output_options) override; + bool requiresTwistShift(bool UNUSED(twist_shift_enabled), [[maybe_unused]] YDirectionType ytype) override { // No Field3Ds require twist-shift, because they cannot be field-aligned @@ -156,6 +164,9 @@ protected: private: /// FCI maps for each of the parallel slices std::vector field_line_maps; + + /// Real-space coordinates of grid points + Field3D R, Z; }; #endif // BOUT_FCITRANSFORM_H diff --git a/src/mesh/tokamak_coordinates.cxx b/src/mesh/tokamak_coordinates.cxx index e5a7636b7d..d90ec43078 100644 --- a/src/mesh/tokamak_coordinates.cxx +++ b/src/mesh/tokamak_coordinates.cxx @@ -4,6 +4,7 @@ #include "bout/bout_types.hxx" #include "bout/field2d.hxx" #include "bout/utils.hxx" +#include "bout/constants.hxx" namespace bout { @@ -17,7 +18,7 @@ namespace bout { TokamakOptions::TokamakOptions(Mesh &mesh) { mesh.get(Rxy, "Rxy"); - // mesh->get(Zxy, "Zxy"); + mesh.get(Zxy, "Zxy"); mesh.get(Bpxy, "Bpxy"); mesh.get(Btxy, "Btxy"); mesh.get(Bxy, "Bxy"); @@ -26,6 +27,34 @@ namespace bout { if (mesh.get(dx, "dpsi")) { dx = mesh.getCoordinates()->dx(); } +// mesh.get(toroidal_angle, "z"); + const auto d_phi = TWOPI / mesh.LocalNz; + auto current_phi = 0.0; + for (int k = 0; k < mesh.LocalNz; k++) { + toroidal_angles.push_back(current_phi); + current_phi += d_phi; + } + } + + Coordinates3D TokamakOptions::cylindrical_coordinates_to_cartesian() { + + auto* mesh = Rxy.getMesh(); + Field3D x = Field3D(0.0, mesh); + Field3D y = Field3D(0.0, mesh); + Field3D z = Field3D(0.0, mesh); + + for (int i = 0; i < Rxy.getNx(); i++) { + for (int j = 0; j < Rxy.getNy(); j++) { + int k = 0; + for (auto angle : toroidal_angles) { + x(i, j, k) = Rxy(i, j) * std::cos(angle); + y(i, j, k) = Rxy(i, j) * std::sin(angle); + z(i, j, k) = Zxy(i, j); + k++; + } + } + } + return Coordinates3D(x, y, z); } void TokamakOptions::normalise(BoutReal Lbar, BoutReal Bbar, BoutReal ShearFactor) { @@ -47,7 +76,7 @@ namespace bout { const BoutReal sign_of_bp = get_sign_of_bp(tokamak_options.Bpxy); - auto *coord = mesh.getCoordinates(); + auto* coord = mesh.getCoordinates(); const auto g11 = SQ(tokamak_options.Rxy * tokamak_options.Bpxy); const auto g22 = 1.0 / SQ(tokamak_options.hthe); diff --git a/src/solver/impls/arkode/arkode.cxx b/src/solver/impls/arkode/arkode.cxx index 23883cc043..6f20ce11a6 100644 --- a/src/solver/impls/arkode/arkode.cxx +++ b/src/solver/impls/arkode/arkode.cxx @@ -229,9 +229,11 @@ int ArkodeSolver::init() { throw BoutException("ARKodeSetUserData failed\n"); } - if (ARKodeSetLinear(arkode_mem, static_cast(set_linear)) - != ARK_SUCCESS) { - throw BoutException("ARKodeSetLinear failed\n"); + if (set_linear) { + constexpr bool is_time_dep = false; + if (ARKodeSetLinear(arkode_mem, is_time_dep) != ARK_SUCCESS) { + throw BoutException("ARKodeSetLinear failed\n"); + } } if (fixed_step) { diff --git a/src/solver/impls/cvode/cvode.cxx b/src/solver/impls/cvode/cvode.cxx index f0d42d39bc..8b81a66f6d 100644 --- a/src/solver/impls/cvode/cvode.cxx +++ b/src/solver/impls/cvode/cvode.cxx @@ -51,6 +51,10 @@ #include #include +#include +#include +#include + #include #include #include @@ -59,6 +63,8 @@ BOUT_ENUM_CLASS(positivity_constraint, none, positive, non_negative, negative, non_positive); +BOUT_ENUM_CLASS(linear_solver, gmres, fgmres, tfqmr, bcgs); + // NOLINTBEGIN(readability-identifier-length) namespace { int cvode_linear_rhs(BoutReal t, N_Vector u, N_Vector du, void* user_data); @@ -344,7 +350,24 @@ int CvodeSolver::init() { const auto prectype = use_precon ? (rightprec ? SUN_PREC_RIGHT : SUN_PREC_LEFT) : SUN_PREC_NONE; - sun_solver = callWithSUNContext(SUNLinSol_SPGMR, suncontext, uvec, prectype, maxl); + + switch ((*options)["linear_solver"] + .doc("Set linear solver type. Default is gmres.") + .withDefault(linear_solver::gmres)) { + case linear_solver::gmres: + sun_solver = callWithSUNContext(SUNLinSol_SPGMR, suncontext, uvec, prectype, maxl); + break; + case linear_solver::fgmres: + sun_solver = callWithSUNContext(SUNLinSol_SPFGMR, suncontext, uvec, prectype, maxl); + break; + case linear_solver::tfqmr: + sun_solver = + callWithSUNContext(SUNLinSol_SPTFQMR, suncontext, uvec, prectype, maxl); + break; + case linear_solver::bcgs: + sun_solver = callWithSUNContext(SUNLinSol_SPBCGS, suncontext, uvec, prectype, maxl); + break; + }; if (sun_solver == nullptr) { throw BoutException("Creating SUNDIALS linear solver failed\n"); } diff --git a/src/solver/impls/imex-bdf2/imex-bdf2.cxx b/src/solver/impls/imex-bdf2/imex-bdf2.cxx index 853d8eaf12..07188110ab 100644 --- a/src/solver/impls/imex-bdf2/imex-bdf2.cxx +++ b/src/solver/impls/imex-bdf2/imex-bdf2.cxx @@ -9,6 +9,7 @@ #include #include #include +#include #include #include @@ -113,7 +114,7 @@ static PetscErrorCode FormFunctionForDifferencing(void* ctx, Vec x, Vec f) { * * This can be a linearised and simplified form of FormFunction */ -static PetscErrorCode FormFunctionForColoring(SNES UNUSED(snes), Vec x, Vec f, +static PetscErrorCode FormFunctionForColoring(void* UNUSED(snes), Vec x, Vec f, void* ctx) { return static_cast(ctx)->snes_function(x, f, true); } @@ -651,9 +652,8 @@ void IMEXBDF2::constructSNES(SNES* snesIn) { // Create data structure for SNESComputeJacobianDefaultColor MatFDColoringCreate(Jmf, iscoloring, &fdcoloring); // Set the function to difference - MatFDColoringSetFunction( - fdcoloring, reinterpret_cast(FormFunctionForColoring), - this); + MatFDColoringSetFunction(fdcoloring, + bout::cast_MatFDColoringFn(FormFunctionForColoring), this); MatFDColoringSetFromOptions(fdcoloring); MatFDColoringSetUp(Jmf, iscoloring, fdcoloring); ISColoringDestroy(&iscoloring); @@ -677,12 +677,7 @@ void IMEXBDF2::constructSNES(SNES* snesIn) { 0, // Number of nonzeros per row in off-diagonal portion of local submatrix nullptr, &Jmf); -#if PETSC_VERSION_GE(3, 4, 0) SNESSetJacobian(*snesIn, Jmf, Jmf, SNESComputeJacobianDefault, this); -#else - // Before 3.4 - SNESSetJacobian(*snesIn, Jmf, Jmf, SNESDefaultComputeJacobian, this); -#endif MatSetOption(Jmf, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); } diff --git a/src/solver/impls/petsc/petsc.cxx b/src/solver/impls/petsc/petsc.cxx index 090087d0de..dd63811f0f 100644 --- a/src/solver/impls/petsc/petsc.cxx +++ b/src/solver/impls/petsc/petsc.cxx @@ -514,9 +514,7 @@ int PetscSolver::init() { CHKERRQ(ierr); ierr = ISColoringDestroy(&iscoloring); CHKERRQ(ierr); - ierr = MatFDColoringSetFunction(matfdcoloring, - reinterpret_cast(solver_f), this); - CHKERRQ(ierr); + throw BoutException("Coloring is not working"); ierr = SNESSetJacobian(snes, J, J, SNESComputeJacobianDefaultColor, matfdcoloring); CHKERRQ(ierr); diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx index e382bbd3f8..2bb163f324 100644 --- a/src/solver/impls/snes/snes.cxx +++ b/src/solver/impls/snes/snes.cxx @@ -1,4 +1,4 @@ -#include "bout/build_defines.hxx" +#include "bout/build_config.hxx" #if BOUT_HAS_PETSC @@ -6,18 +6,57 @@ #include #include +#include #include +#include #include -#include #include #include #include -#include "petscmat.h" #include "petscsnes.h" +class ColoringStencil { +private: + bool static isInSquare(int const i, int const j, int const n_square) { + return std::abs(i) <= n_square && std::abs(j) <= n_square; + } + bool static isInCross(int const i, int const j, int const n_cross) { + if (i == 0) { + return std::abs(j) <= n_cross; + } + if (j == 0) { + return std::abs(i) <= n_cross; + } + return false; + } + bool static isInTaxi(int const i, int const j, int const n_taxi) { + return std::abs(i) + std::abs(j) <= n_taxi; + } + +public: + auto static getOffsets(int n_square, int n_taxi, int n_cross) { + ASSERT2(n_square >= 0 && n_cross >= 0 && n_taxi >= 0 + && n_square + n_cross + n_taxi > 0); + auto inside = [&](int i, int j) { + return isInSquare(i, j, n_square) || isInTaxi(i, j, n_taxi) + || isInCross(i, j, n_cross); + }; + std::vector> xy_offsets; + auto loop_bound = std::max({n_square, n_taxi, n_cross}); + for (int i = -loop_bound; i <= loop_bound; ++i) { + for (int j = -loop_bound; j <= loop_bound; ++j) { + if (inside(i, j)) { + xy_offsets.emplace_back(i, j); + } + } + } + return xy_offsets; + } +}; + /* * PETSc callback function, which evaluates the nonlinear * function to be solved by SNES. @@ -43,7 +82,7 @@ static PetscErrorCode FormFunctionForDifferencing(void* ctx, Vec x, Vec f) { * * This can be a linearised and simplified form of FormFunction */ -static PetscErrorCode FormFunctionForColoring(SNES UNUSED(snes), Vec x, Vec f, +static PetscErrorCode FormFunctionForColoring(void* UNUSED(snes), Vec x, Vec f, void* ctx) { return static_cast(ctx)->snes_function(x, f, true); } @@ -79,12 +118,32 @@ SNESSolver::SNESSolver(Options* opts) maxits((*options)["max_nonlinear_iterations"] .doc("Maximum number of nonlinear iterations per SNES solve") .withDefault(50)), + maxf((*options)["maxf"] + .doc("Maximum number of function evaluations per SNES solve") + .withDefault(10000)), lower_its((*options)["lower_its"] .doc("Iterations below which the next timestep is increased") .withDefault(static_cast(maxits * 0.5))), upper_its((*options)["upper_its"] .doc("Iterations above which the next timestep is reduced") .withDefault(static_cast(maxits * 0.8))), + timestep_factor_on_failure((*options)["timestep_factor_on_failure"] + .doc("Multiply timestep on convergence failure") + .withDefault(0.5)), + timestep_factor_on_upper_its( + (*options)["timestep_factor_on_upper_its"] + .doc("Multiply timestep if iterations exceed upper_its") + .withDefault(0.9)), + timestep_factor_on_lower_its( + (*options)["timestep_factor_on_lower_its"] + .doc("Multiply timestep if iterations are below lower_its") + .withDefault(1.4)), + pid_controller( + (*options)["pid_controller"].doc("Use PID controller?").withDefault(false)), + target_its((*options)["target_its"].doc("Target snes iterations").withDefault(7)), + kP((*options)["kP"].doc("Proportional PID parameter").withDefault(0.7)), + kI((*options)["kI"].doc("Integral PID parameter").withDefault(0.3)), + kD((*options)["kD"].doc("Derivative PID parameter").withDefault(0.2)), diagnose( (*options)["diagnose"].doc("Print additional diagnostics").withDefault(false)), diagnose_failures((*options)["diagnose_failures"] @@ -122,7 +181,7 @@ SNESSolver::SNESSolver(Options* opts) .withDefault(false)), matrix_free_operator((*options)["matrix_free_operator"] .doc("Use matrix free Jacobian-vector operator?") - .withDefault(true)), + .withDefault(false)), lag_jacobian((*options)["lag_jacobian"] .doc("Re-use the Jacobian this number of SNES iterations") .withDefault(50)), @@ -259,13 +318,21 @@ int SNESSolver::init() { SNESSetJacobian(snes, Jmf, Jmf, MatMFFDComputeJacobian, this); } else { - // Calculate the Jacobian using finite differences. - // The finite difference Jacobian (Jfd) may be used for both operator - // and preconditioner or, if matrix_free_operator, in only the preconditioner. + // Calculate the Jacobian using finite differences. The finite + // difference Jacobian (Jfd) may be used for both operator and + // preconditioner or, if matrix_free_operator, in only the + // preconditioner. + + // Create a vector to store interpolated output solution + // Used so that the timestep does not have to be adjusted, + // because that would require updating the preconditioner. + ierr = VecDuplicate(snes_x, &output_x); + CHKERRQ(ierr); + if (use_coloring) { - // Use matrix coloring - // This greatly reduces the number of times the rhs() function needs - // to be evaluated when calculating the Jacobian. + // Use matrix coloring. + // This greatly reduces the number of times the rhs() function + // needs to be evaluated when calculating the Jacobian. // Use global mesh for now Mesh* mesh = bout::globals::mesh; @@ -279,225 +346,165 @@ int SNESSolver::init() { output_progress.write("Setting Jacobian matrix sizes\n"); - int n2d = f2d.size(); - int n3d = f3d.size(); + const int n2d = f2d.size(); + const int n3d = f3d.size(); // Set size of Matrix on each processor to nlocal x nlocal MatCreate(BoutComm::get(), &Jfd); + MatSetOption(Jfd, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE); MatSetSizes(Jfd, nlocal, nlocal, PETSC_DETERMINE, PETSC_DETERMINE); MatSetFromOptions(Jfd); - - std::vector d_nnz(nlocal); - std::vector o_nnz(nlocal); - - // Set values for most points - const int ncells_x = (mesh->LocalNx > 1) ? 2 : 0; - const int ncells_y = (mesh->LocalNy > 1) ? 2 : 0; - const int ncells_z = (mesh->LocalNz > 1) ? 2 : 0; - - const auto star_pattern = (1 + ncells_x + ncells_y) * (n3d + n2d) + ncells_z * n3d; - - // Offsets. Start with the central cell - std::vector> xyoffsets{{0, 0}}; - if (ncells_x != 0) { - // Stencil includes points in X - xyoffsets.push_back({-1, 0}); - xyoffsets.push_back({1, 0}); - } - if (ncells_y != 0) { - // Stencil includes points in Y - xyoffsets.push_back({0, -1}); - xyoffsets.push_back({0, 1}); + // Determine which row/columns of the matrix are locally owned + int Istart, Iend; + MatGetOwnershipRange(Jfd, &Istart, &Iend); + // Convert local into global indices + // Note: Not in the boundary cells, to keep -1 values + for (const auto& i : mesh->getRegion3D("RGN_NOBNDRY")) { + index[i] += Istart; } + // Now communicate to fill guard cells + mesh->communicate(index); - output_info.write("Star pattern: {} non-zero entries\n", star_pattern); - for (int i = 0; i < nlocal; i++) { - // Non-zero elements on this processor - d_nnz[i] = star_pattern; - // Non-zero elements on neighboring processor - o_nnz[i] = 0; + // Non-zero elements on this processor + std::vector d_nnz; + std::vector o_nnz; + auto n_square = (*options)["stencil:square"] + .doc("Extent of stencil (square)") + .withDefault(0); + auto n_taxi = (*options)["stencil:taxi"] + .doc("Extent of stencil (taxi-cab norm)") + .withDefault(0); + auto n_cross = (*options)["stencil:cross"] + .doc("Extent of stencil (cross)") + .withDefault(0); + // Set n_taxi 2 if nothing else is set + // Probably a better way to do this + if (n_square == 0 && n_taxi == 0 && n_cross == 0) { + output_info.write("Setting solver:stencil:taxi = 2\n"); + n_taxi = 2; } - // X boundaries - if (ncells_x != 0) { - if (mesh->firstX()) { - // Lower X boundary - for (int y = mesh->ystart; y <= mesh->yend; y++) { - for (int z = 0; z < mesh->LocalNz; z++) { - const int localIndex = ROUND(index(mesh->xstart, y, z)); - ASSERT2((localIndex >= 0) && (localIndex < nlocal)); - const int num_fields = (z == 0) ? n2d + n3d : n3d; - for (int i = 0; i < num_fields; i++) { - d_nnz[localIndex + i] -= (n3d + n2d); - } - } - } - } else { - // On another processor - for (int y = mesh->ystart; y <= mesh->yend; y++) { - for (int z = 0; z < mesh->LocalNz; z++) { - const int localIndex = ROUND(index(mesh->xstart, y, z)); - ASSERT2((localIndex >= 0) && (localIndex < nlocal)); - const int num_fields = (z == 0) ? n2d + n3d : n3d; - for (int i = 0; i < num_fields; i++) { - d_nnz[localIndex + i] -= (n3d + n2d); - o_nnz[localIndex + i] += (n3d + n2d); - } - } - } - } - if (mesh->lastX()) { - // Upper X boundary - for (int y = mesh->ystart; y <= mesh->yend; y++) { - for (int z = 0; z < mesh->LocalNz; z++) { - const int localIndex = ROUND(index(mesh->xend, y, z)); - ASSERT2((localIndex >= 0) && (localIndex < nlocal)); - const int num_fields = (z == 0) ? n2d + n3d : n3d; - for (int i = 0; i < num_fields; i++) { - d_nnz[localIndex + i] -= (n3d + n2d); - } - } - } - } else { - // On another processor + auto const xy_offsets = ColoringStencil::getOffsets(n_square, n_taxi, n_cross); + { + // This is ugly but can't think of a better and robust way to + // count the non-zeros for some arbitrary stencil + // effectively the same loop as the one that sets the non-zeros below + std::vector> d_nnz_map2d(nlocal); + std::vector> o_nnz_map2d(nlocal); + std::vector> d_nnz_map3d(nlocal); + std::vector> o_nnz_map3d(nlocal); + // Loop over every element in 2D to count the *unique* non-zeros + for (int x = mesh->xstart; x <= mesh->xend; x++) { for (int y = mesh->ystart; y <= mesh->yend; y++) { - for (int z = 0; z < mesh->LocalNz; z++) { - const int localIndex = ROUND(index(mesh->xend, y, z)); - ASSERT2((localIndex >= 0) && (localIndex < nlocal)); - const int num_fields = (z == 0) ? n2d + n3d : n3d; - for (int i = 0; i < num_fields; i++) { - d_nnz[localIndex + i] -= (n3d + n2d); - o_nnz[localIndex + i] += (n3d + n2d); - } - } - } - } - } - // Y boundaries - if (ncells_y != 0) { - for (int x = mesh->xstart; x <= mesh->xend; x++) { - // Default to no boundary - // NOTE: This assumes that communications in Y are to other - // processors. If Y is communicated with this processor (e.g. NYPE=1) - // then this will result in PETSc warnings about out of range allocations - - // z = 0 case - int localIndex = ROUND(index(x, mesh->ystart, 0)); - ASSERT2(localIndex >= 0); - - // All 2D and 3D fields - for (int i = 0; i < n2d + n3d; i++) { - o_nnz[localIndex + i] += (n3d + n2d); - d_nnz[localIndex + i] -= (n3d + n2d); - } + const int ind0 = ROUND(index(x, y, 0)) - Istart; - for (int z = 1; z < mesh->LocalNz; z++) { - localIndex = ROUND(index(x, mesh->ystart, z)); + // 2D fields + for (int i = 0; i < n2d; i++) { + const PetscInt row = ind0 + i; + // Loop through each point in the stencil + for (const auto& [x_off, y_off] : xy_offsets) { + const int xi = x + x_off; + const int yi = y + y_off; + if ((xi < 0) || (yi < 0) || (xi >= mesh->LocalNx) + || (yi >= mesh->LocalNy)) { + continue; + } - // Only 3D fields - for (int i = 0; i < n3d; i++) { - o_nnz[localIndex + i] += (n3d + n2d); - d_nnz[localIndex + i] -= (n3d + n2d); + const int ind2 = ROUND(index(xi, yi, 0)); + if (ind2 < 0) { + continue; // A boundary point + } + + // Depends on all variables on this cell + for (int j = 0; j < n2d; j++) { + const PetscInt col = ind2 + j; + if (col >= Istart && col < Iend) { + d_nnz_map2d[row].insert(col); + } else { + o_nnz_map2d[row].insert(col); + } + } + } } - } + // 3D fields + for (int z = 0; z < mesh->LocalNz; z++) { + const int ind = ROUND(index(x, y, z)) - Istart; - // z = 0 case - localIndex = ROUND(index(x, mesh->yend, 0)); - // All 2D and 3D fields - for (int i = 0; i < n2d + n3d; i++) { - o_nnz[localIndex + i] += (n3d + n2d); - d_nnz[localIndex + i] -= (n3d + n2d); - } + for (int i = 0; i < n3d; i++) { + PetscInt row = ind + i; + if (z == 0) { + row += n2d; + } - for (int z = 1; z < mesh->LocalNz; z++) { - localIndex = ROUND(index(x, mesh->yend, z)); + // Depends on 2D fields + for (int j = 0; j < n2d; j++) { + const PetscInt col = ind0 + j; + if (col >= Istart && col < Iend) { + d_nnz_map2d[row].insert(col); + } else { + o_nnz_map2d[row].insert(col); + } + } - // Only 3D fields - for (int i = 0; i < n3d; i++) { - o_nnz[localIndex + i] += (n3d + n2d); - d_nnz[localIndex + i] -= (n3d + n2d); - } - } - } + // Star pattern + for (const auto& [x_off, y_off] : xy_offsets) { + const int xi = x + x_off; + const int yi = y + y_off; - for (RangeIterator it = mesh->iterateBndryLowerY(); !it.isDone(); it++) { - // A boundary, so no communication + if ((xi < 0) || (yi < 0) || (xi >= mesh->LocalNx) + || (yi >= mesh->LocalNy)) { + continue; + } - // z = 0 case - int localIndex = ROUND(index(it.ind, mesh->ystart, 0)); - if (localIndex < 0) { - // This can occur because it.ind includes values in x boundary e.g. x=0 - continue; - } - // All 2D and 3D fields - for (int i = 0; i < n2d + n3d; i++) { - o_nnz[localIndex + i] -= (n3d + n2d); - } + int ind2 = ROUND(index(xi, yi, 0)); + if (ind2 < 0) { + continue; // Boundary point + } - for (int z = 1; z < mesh->LocalNz; z++) { - int localIndex = ROUND(index(it.ind, mesh->ystart, z)); + if (z == 0) { + ind2 += n2d; + } - // Only 3D fields - for (int i = 0; i < n3d; i++) { - o_nnz[localIndex + i] -= (n3d + n2d); + // 3D fields on this cell + for (int j = 0; j < n3d; j++) { + const PetscInt col = ind2 + j; + if (col >= Istart && col < Iend) { + d_nnz_map3d[row].insert(col); + } else { + o_nnz_map3d[row].insert(col); + } + } + } + } } } } - for (RangeIterator it = mesh->iterateBndryUpperY(); !it.isDone(); it++) { - // A boundary, so no communication - - // z = 0 case - const int localIndex = ROUND(index(it.ind, mesh->yend, 0)); - if (localIndex < 0) { - continue; // Out of domain - } - - // All 2D and 3D fields - for (int i = 0; i < n2d + n3d; i++) { - o_nnz[localIndex + i] -= (n3d + n2d); - } - - for (int z = 1; z < mesh->LocalNz; z++) { - const int localIndex = ROUND(index(it.ind, mesh->yend, z)); + d_nnz.reserve(nlocal); + d_nnz.reserve(nlocal); - // Only 3D fields - for (int i = 0; i < n3d; i++) { - o_nnz[localIndex + i] -= (n3d + n2d); - } - } + for (int i = 0; i < nlocal; ++i) { + // Assume all elements in the z direction are potentially coupled + d_nnz.emplace_back(d_nnz_map3d[i].size() * mesh->LocalNz + + d_nnz_map2d[i].size()); + o_nnz.emplace_back(o_nnz_map3d[i].size() * mesh->LocalNz + + o_nnz_map2d[i].size()); } } output_progress.write("Pre-allocating Jacobian\n"); - // Pre-allocate MatMPIAIJSetPreallocation(Jfd, 0, d_nnz.data(), 0, o_nnz.data()); MatSeqAIJSetPreallocation(Jfd, 0, d_nnz.data()); MatSetUp(Jfd); MatSetOption(Jfd, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE); - // Determine which row/columns of the matrix are locally owned - int Istart, Iend; - MatGetOwnershipRange(Jfd, &Istart, &Iend); - - // Convert local into global indices - // Note: Not in the boundary cells, to keep -1 values - for (const auto& i : mesh->getRegion3D("RGN_NOBNDRY")) { - index[i] += Istart; - } - - // Now communicate to fill guard cells - mesh->communicate(index); - ////////////////////////////////////////////////// // Mark non-zero entries output_progress.write("Marking non-zero Jacobian entries\n"); - - const PetscScalar val = 1.0; - + PetscScalar val = 1.0; for (int x = mesh->xstart; x <= mesh->xend; x++) { for (int y = mesh->ystart; y <= mesh->yend; y++) { @@ -507,35 +514,31 @@ int SNESSolver::init() { for (int i = 0; i < n2d; i++) { const PetscInt row = ind0 + i; - // Loop through each point in the 5-point stencil - for (const auto& xyoffset : xyoffsets) { - const int xi = x + xyoffset.first; - const int yi = y + xyoffset.second; - + // Loop through each point in the stencil + for (const auto& [x_off, y_off] : xy_offsets) { + const int xi = x + x_off; + const int yi = y + y_off; if ((xi < 0) || (yi < 0) || (xi >= mesh->LocalNx) || (yi >= mesh->LocalNy)) { continue; } - const int ind2 = ROUND(index(xi, yi, 0)); - + int ind2 = ROUND(index(xi, yi, 0)); if (ind2 < 0) { continue; // A boundary point } // Depends on all variables on this cell for (int j = 0; j < n2d; j++) { - const PetscInt col = ind2 + j; + PetscInt col = ind2 + j; ierr = MatSetValues(Jfd, 1, &row, 1, &col, &val, INSERT_VALUES); CHKERRQ(ierr); } } } - // 3D fields for (int z = 0; z < mesh->LocalNz; z++) { - - const int ind = ROUND(index(x, y, z)); + int ind = ROUND(index(x, y, z)); for (int i = 0; i < n3d; i++) { PetscInt row = ind + i; @@ -545,73 +548,48 @@ int SNESSolver::init() { // Depends on 2D fields for (int j = 0; j < n2d; j++) { - const PetscInt col = ind0 + j; + PetscInt col = ind0 + j; ierr = MatSetValues(Jfd, 1, &row, 1, &col, &val, INSERT_VALUES); CHKERRQ(ierr); } // Star pattern - for (const auto& xyoffset : xyoffsets) { - const int xi = x + xyoffset.first; - const int yi = y + xyoffset.second; + for (const auto& [x_off, y_off] : xy_offsets) { + int xi = x + x_off; + int yi = y + y_off; if ((xi < 0) || (yi < 0) || (xi >= mesh->LocalNx) || (yi >= mesh->LocalNy)) { continue; } - - int ind2 = ROUND(index(xi, yi, z)); - if (ind2 < 0) { - continue; // Boundary point - } - - if (z == 0) { - ind2 += n2d; - } - - // 3D fields on this cell - for (int j = 0; j < n3d; j++) { - const PetscInt col = ind2 + j; - ierr = MatSetValues(Jfd, 1, &row, 1, &col, &val, INSERT_VALUES); - if (ierr != 0) { - output.write("ERROR: {} : ({}, {}) -> ({}, {}) : {} -> {}\n", row, x, - y, xi, yi, ind2, ind2 + n3d - 1); + for (int zi = 0; zi < mesh->LocalNz; ++zi) { + int ind2 = ROUND(index(xi, yi, zi)); + if (ind2 < 0) { + continue; // Boundary point } - CHKERRQ(ierr); - } - } - int nz = mesh->LocalNz; - if (nz > 1) { - // Multiple points in z + if (z == 0) { + ind2 += n2d; + } - int zp = (z + 1) % nz; + // 3D fields on this cell + for (int j = 0; j < n3d; j++) { + PetscInt col = ind2 + j; + ierr = MatSetValues(Jfd, 1, &row, 1, &col, &val, INSERT_VALUES); - int ind2 = ROUND(index(x, y, zp)); - if (zp == 0) { - ind2 += n2d; - } - for (int j = 0; j < n3d; j++) { - const PetscInt col = ind2 + j; - ierr = MatSetValues(Jfd, 1, &row, 1, &col, &val, INSERT_VALUES); - CHKERRQ(ierr); - } - - int zm = (z - 1 + nz) % nz; - ind2 = ROUND(index(x, y, zm)); - if (zm == 0) { - ind2 += n2d; - } - for (int j = 0; j < n3d; j++) { - const PetscInt col = ind2 + j; - ierr = MatSetValues(Jfd, 1, &row, 1, &col, &val, INSERT_VALUES); - CHKERRQ(ierr); + if (ierr != 0) { + output.write("ERROR: {} {} : ({}, {}) -> ({}, {}) : {} -> {}\n", + row, x, y, xi, yi, ind2, ind2 + n3d - 1); + } + CHKERRQ(ierr); + } } } } } } } + // Finished marking non-zero entries output_progress.write("Assembling Jacobian matrix\n"); @@ -620,6 +598,30 @@ int SNESSolver::init() { MatAssemblyBegin(Jfd, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(Jfd, MAT_FINAL_ASSEMBLY); + { + // Test if the matrix is symmetric + // Values are 0 or 1 so tolerance (1e-5) shouldn't matter + PetscBool symmetric; + ierr = MatIsSymmetric(Jfd, 1e-5, &symmetric); + CHKERRQ(ierr); + if (!symmetric) { + output_warn.write("Jacobian pattern is not symmetric\n"); + } + } + + // The above can miss entries around the X-point branch cut: + // The diagonal terms are complicated because moving in X then Y + // is different from moving in Y then X at the X-point. + // Making sure the colouring matrix is symmetric does not + // necessarily give the correct stencil but may help. + if ((*options)["force_symmetric_coloring"] + .doc("Modifies coloring matrix to force it to be symmetric") + .withDefault(false)) { + Mat Jfd_T; + MatCreateTranspose(Jfd, &Jfd_T); + MatAXPY(Jfd, 1, Jfd_T, DIFFERENT_NONZERO_PATTERN); + } + output_progress.write("Creating Jacobian coloring\n"); updateColoring(); @@ -655,14 +657,21 @@ int SNESSolver::init() { // Note: If the 'Amat' Jacobian is matrix free, SNESComputeJacobian // always updates its reference 'u' vector every nonlinear iteration SNESSetLagJacobian(snes, lag_jacobian); - // Set Jacobian and preconditioner to persist across time steps - SNESSetLagJacobianPersists(snes, PETSC_TRUE); - SNESSetLagPreconditionerPersists(snes, PETSC_TRUE); + if (pid_controller) { + nl_its_prev = target_its; + nl_its_prev2 = target_its; + SNESSetLagJacobianPersists(snes, PETSC_FALSE); + SNESSetLagPreconditionerPersists(snes, PETSC_FALSE); + } else { + // Set Jacobian and preconditioner to persist across time steps + SNESSetLagJacobianPersists(snes, PETSC_TRUE); + SNESSetLagPreconditionerPersists(snes, PETSC_TRUE); + } SNESSetLagPreconditioner(snes, 1); // Rebuild when Jacobian is rebuilt } // Set tolerances - SNESSetTolerances(snes, atol, rtol, stol, maxits, PETSC_DEFAULT); + SNESSetTolerances(snes, atol, rtol, stol, maxits, maxf); // Force SNES to take at least one nonlinear iteration. // This may prevent the solver from getting stuck in false steady state conditions @@ -758,14 +767,18 @@ int SNESSolver::run() { CHKERRQ(ierr); } + BoutReal target = simtime; for (int s = 0; s < getNumberOutputSteps(); s++) { - BoutReal target = simtime + getOutputTimestep(); + target += getOutputTimestep(); bool looping = true; int snes_failures = 0; // Count SNES convergence failures + int steps_since_snes_failure = 0; int saved_jacobian_lag = 0; int loop_count = 0; do { + if (simtime >= target) + break; // Could happen if step over multiple outputs if (scale_vars) { // Individual variable scaling // Note: If variables are rescaled then the Jacobian columns @@ -843,9 +856,15 @@ int SNESSolver::run() { dt = timestep; looping = true; if (simtime + dt >= target) { - // Ensure that the timestep goes to the next output time and then stops + // Note: When the timestep is changed the preconditioner needs to be updated + // => Step over the output time and interpolate if not matrix free + + if (matrix_free) { + // Ensure that the timestep goes to the next output time and then stops. + // This avoids the need to interpolate + dt = target - simtime; + } looping = false; - dt = target - simtime; } if (predictor and (time1 > 0.0)) { @@ -856,6 +875,10 @@ int SNESSolver::run() { VecAXPBY(snes_x, -beta, (1. + beta), x1); } + if (pid_controller) { + SNESSetLagJacobian(snes, lag_jacobian); + } + // Run the solver PetscErrorCode ierr = SNESSolve(snes, nullptr, snes_x); @@ -892,9 +915,10 @@ int SNESSolver::run() { } ++snes_failures; + steps_since_snes_failure = 0; // Try a smaller timestep - timestep /= 2.0; + timestep *= timestep_factor_on_failure; // Restore state VecCopy(x0, snes_x); @@ -912,7 +936,9 @@ int SNESSolver::run() { updateColoring(); jacobian_pruned = false; // Reset flag. Will be set after pruning. } + if (saved_jacobian_lag == 0) { + // This triggers a Jacobian recalculation SNESGetLagJacobian(snes, &saved_jacobian_lag); SNESSetLagJacobian(snes, 1); } @@ -984,6 +1010,7 @@ int SNESSolver::run() { } simtime += dt; + ++steps_since_snes_failure; if (diagnose) { // Gather and print diagnostic information @@ -993,7 +1020,6 @@ int SNESSolver::run() { simtime, timestep, nl_its, lin_its, static_cast(reason)); if (snes_failures > 0) { output.write(", SNES failures: {}", snes_failures); - snes_failures = 0; } output.write("\n"); } @@ -1044,24 +1070,91 @@ int SNESSolver::run() { #endif // PETSC_VERSION_GE(3,20,0) if (looping) { - if (nl_its <= lower_its) { - // Increase timestep slightly - timestep *= 1.1; - if (timestep > max_timestep) { - timestep = max_timestep; + if (pid_controller) { + // Changing the timestep. + // Note: The preconditioner depends on the timestep, + // so we recalculate the jacobian and the preconditioner + // every time the timestep changes + + timestep = pid(timestep, nl_its); + + // NOTE(malamast): Do we really need this? + // Recompute Jacobian (for now) + if (saved_jacobian_lag == 0) { + SNESGetLagJacobian(snes, &saved_jacobian_lag); + SNESSetLagJacobian(snes, 1); + } + + } else { + + // Consider changing the timestep. + // Note: The preconditioner depends on the timestep, + // so if it is not recalculated the it will be less + // effective. + if ((nl_its <= lower_its) && (timestep < max_timestep) + && (steps_since_snes_failure > 2)) { + // Increase timestep slightly + timestep *= timestep_factor_on_lower_its; + + timestep = std::min(timestep, max_timestep); + + // Note: Setting the SNESJacobianFn to NULL retains + // previously set evaluation function. + // + // The SNES Jacobian is a combination of the RHS Jacobian + // and a factor involving the timestep. + // Depends on equation_form + // -> Probably call SNESSetJacobian(snes, Jfd, Jfd, NULL, fdcoloring); + + if (static_cast(lin_its) / nl_its > 4) { + // Recompute Jacobian (for now) + if (saved_jacobian_lag == 0) { + SNESGetLagJacobian(snes, &saved_jacobian_lag); + SNESSetLagJacobian(snes, 1); + } + } + + } else if (nl_its >= upper_its) { + // Reduce timestep slightly + timestep *= timestep_factor_on_upper_its; + + // Recompute Jacobian + if (saved_jacobian_lag == 0) { + SNESGetLagJacobian(snes, &saved_jacobian_lag); + SNESSetLagJacobian(snes, 1); + } } - } else if (nl_its >= upper_its) { - // Reduce timestep slightly - timestep *= 0.9; } } + snes_failures = 0; } while (looping); + if (!matrix_free) { + ASSERT2(simtime >= target); + ASSERT2(simtime - dt < target); + // Stepped over output timestep => Interpolate + // snes_x is the solution at t = simtime + // x0 is the solution at t = simtime - dt + // Calculate output_x at t = target + VecCopy(snes_x, output_x); + + // Note: If simtime = target then alpha = 0 + // and output_x = snes_x + BoutReal alpha = (simtime - target) / dt; + + // output_x <- alpha * x0 + (1 - alpha) * output_x + VecAXPBY(output_x, alpha, 1. - alpha, x0); + + } else { + // Timestep was adjusted to hit target output time + output_x = snes_x; + } + // Put the result into variables if (scale_vars) { - // scaled_x <- snes_x * var_scaling_factors - int ierr = VecPointwiseMult(scaled_x, snes_x, var_scaling_factors); + // scaled_x <- output_x * var_scaling_factors + int ierr = VecPointwiseMult(scaled_x, output_x, var_scaling_factors); CHKERRQ(ierr); const BoutReal* xdata = nullptr; @@ -1072,15 +1165,15 @@ int SNESSolver::run() { CHKERRQ(ierr); } else { const BoutReal* xdata = nullptr; - int ierr = VecGetArrayRead(snes_x, &xdata); + int ierr = VecGetArrayRead(output_x, &xdata); CHKERRQ(ierr); load_vars(const_cast(xdata)); - ierr = VecRestoreArrayRead(snes_x, &xdata); + ierr = VecRestoreArrayRead(output_x, &xdata); CHKERRQ(ierr); } - run_rhs(simtime); // Run RHS to calculate auxilliary variables + run_rhs(target); // Run RHS to calculate auxilliary variables - if (call_monitors(simtime, s, getNumberOutputSteps()) != 0) { + if (call_monitors(target, s, getNumberOutputSteps()) != 0) { break; // User signalled to quit } } @@ -1321,7 +1414,10 @@ void SNESSolver::updateColoring() { // Re-calculate the coloring MatColoring coloring = NULL; MatColoringCreate(Jfd, &coloring); - MatColoringSetType(coloring, MATCOLORINGSL); + // MatColoringSetType(coloring, MATCOLORINGSL); // Serial algorithm. Better for smale-to-medium size problems. + MatColoringSetType( + coloring, MATCOLORINGGREEDY); // Parallel algorith. Better for large parallel runs + // MatColoringSetType(coloring, MATCOLORINGJP); // This didn't work MatColoringSetFromOptions(coloring); // Calculate new index sets @@ -1332,8 +1428,8 @@ void SNESSolver::updateColoring() { // Replace the old coloring with the new one MatFDColoringDestroy(&fdcoloring); MatFDColoringCreate(Jfd, iscoloring, &fdcoloring); - MatFDColoringSetFunction( - fdcoloring, reinterpret_cast(FormFunctionForColoring), this); + MatFDColoringSetFunction(fdcoloring, + bout::cast_MatFDColoringFn(FormFunctionForColoring), this); MatFDColoringSetFromOptions(fdcoloring); MatFDColoringSetUp(Jfd, iscoloring, fdcoloring); ISColoringDestroy(&iscoloring); @@ -1347,4 +1443,25 @@ void SNESSolver::updateColoring() { } } +BoutReal SNESSolver::pid(BoutReal timestep, int nl_its) { + + /* ---------- multiplicative PID factors ---------- */ + const BoutReal facP = std::pow(double(target_its) / double(nl_its), kP); + const BoutReal facI = std::pow(double(nl_its_prev) / double(nl_its), kI); + const BoutReal facD = std::pow(double(nl_its_prev) * double(nl_its_prev) + / double(nl_its) / double(nl_its_prev2), + kD); + + // clamp groth factor to avoid huge changes + const BoutReal fac = std::clamp(facP * facI * facD, 0.2, 5.0); + + /* ---------- update timestep and history ---------- */ + const BoutReal dt_new = std::min(timestep * fac, max_timestep); + + nl_its_prev2 = nl_its_prev; + nl_its_prev = nl_its; + + return dt_new; +} + #endif // BOUT_HAS_PETSC diff --git a/src/solver/impls/snes/snes.hxx b/src/solver/impls/snes/snes.hxx index 17050ad775..bd942f09ff 100644 --- a/src/solver/impls/snes/snes.hxx +++ b/src/solver/impls/snes/snes.hxx @@ -4,7 +4,7 @@ * using PETSc for the SNES interface * ************************************************************************** - * Copyright 2015-2024 BOUT++ contributors + * Copyright 2015-2025 BOUT++ contributors * * Contact: Ben Dudson, dudson2@llnl.gov * @@ -98,10 +98,28 @@ private: BoutReal atol; ///< Absolute tolerance BoutReal rtol; ///< Relative tolerance BoutReal stol; ///< Convergence tolerance + int maxf; ///< Maximum number of function evaluations allowed in the solver (default: 10000) int maxits; ///< Maximum nonlinear iterations int lower_its, upper_its; ///< Limits on iterations for timestep adjustment + BoutReal timestep_factor_on_failure; + BoutReal timestep_factor_on_upper_its; + BoutReal timestep_factor_on_lower_its; + + ///< PID controller parameters + bool pid_controller; ///< Use PID controller? + int target_its; ///< Target number of nonlinear iterations for the PID controller. + ///< Use with caution! Not tested values. + BoutReal kP; ///< (0.6 - 0.8) Proportional parameter (main response to current step) + BoutReal kI; ///< (0.2 - 0.4) Integral parameter (smooths history of changes) + BoutReal kD; ///< (0.1 - 0.3) Derivative (dampens oscillation - optional) + + int nl_its_prev; + int nl_its_prev2; + + BoutReal pid(BoutReal timestep, int nl_its); ///< Updates the timestep + bool diagnose; ///< Output additional diagnostics bool diagnose_failures; ///< Print diagnostics on SNES failures @@ -116,14 +134,15 @@ private: Vec snes_x; ///< Result of SNES Vec x0; ///< Solution at start of current timestep Vec delta_x; ///< Change in solution + Vec output_x; ///< Solution to output. Used if interpolating. bool predictor; ///< Use linear predictor? Vec x1; ///< Previous solution BoutReal time1{-1.0}; ///< Time of previous solution - SNES snes; ///< SNES context - Mat Jmf; ///< Matrix Free Jacobian - Mat Jfd; ///< Finite Difference Jacobian + SNES snes; ///< SNES context + Mat Jmf; ///< Matrix Free Jacobian + Mat Jfd; ///< Finite Difference Jacobian MatFDColoring fdcoloring{nullptr}; ///< Matrix coloring context ///< Jacobian evaluation @@ -135,10 +154,10 @@ private: std::string pc_hypre_type; ///< Hypre preconditioner type std::string line_search_type; ///< Line search type - bool matrix_free; ///< Use matrix free Jacobian - bool matrix_free_operator; ///< Use matrix free Jacobian in the operator? - int lag_jacobian; ///< Re-use Jacobian - bool use_coloring; ///< Use matrix coloring + bool matrix_free; ///< Use matrix free Jacobian + bool matrix_free_operator; ///< Use matrix free Jacobian in the operator? + int lag_jacobian; ///< Re-use Jacobian + bool use_coloring; ///< Use matrix coloring bool jacobian_recalculated; ///< Flag set when Jacobian is recalculated bool prune_jacobian; ///< Remove small elements in the Jacobian? diff --git a/src/sys/options/options_adios.cxx b/src/sys/options/options_adios.cxx index b3acbaada6..09797647ed 100644 --- a/src/sys/options/options_adios.cxx +++ b/src/sys/options/options_adios.cxx @@ -307,6 +307,12 @@ void OptionsADIOS::verifyTimesteps() const { return; } +const std::vector DIMS_NONE; +const std::vector DIMS_X = {"x"}; +const std::vector DIMS_XY = {"x", "y"}; +const std::vector DIMS_XZ = {"x", "z"}; +const std::vector DIMS_XYZ = {"x", "y", "z"}; + /// Visit a variant type, and put the data into a NcVar struct ADIOSPutVarVisitor { ADIOSPutVarVisitor(const std::string& name, ADIOSStream& stream) @@ -388,7 +394,8 @@ void ADIOSPutVarVisitor::operator()(const Field2D& value) { adios2::Dims memCount = {static_cast(value.getNx()), static_cast(value.getNy())}; - adios2::Variable var = stream.GetArrayVariable(varname, shape); + adios2::Variable var = + stream.GetArrayVariable(varname, shape, DIMS_XY, BoutComm::rank()); var.SetSelection({start, count}); var.SetMemorySelection({memStart, memCount}); stream.engine.Put(var, &value(0, 0)); @@ -425,7 +432,8 @@ void ADIOSPutVarVisitor::operator()(const Field3D& value) { static_cast(value.getNy()), static_cast(value.getNz())}; - adios2::Variable var = stream.GetArrayVariable(varname, shape); + adios2::Variable var = + stream.GetArrayVariable(varname, shape, DIMS_XYZ, BoutComm::rank()); var.SetSelection({start, count}); var.SetMemorySelection({memStart, memCount}); stream.engine.Put(var, &value(0, 0, 0)); @@ -457,7 +465,8 @@ void ADIOSPutVarVisitor::operator()(const FieldPerp& value) { adios2::Dims memCount = {static_cast(value.getNx()), static_cast(value.getNz())}; - adios2::Variable var = stream.GetArrayVariable(varname, shape); + adios2::Variable var = + stream.GetArrayVariable(varname, shape, DIMS_XZ, BoutComm::rank()); var.SetSelection({start, count}); var.SetMemorySelection({memStart, memCount}); stream.engine.Put(var, &value(0, 0)); @@ -469,7 +478,8 @@ void ADIOSPutVarVisitor::operator()>(const Array& valu adios2::Dims shape = {(size_t)BoutComm::size(), (size_t)value.size()}; adios2::Dims start = {(size_t)BoutComm::rank(), 0}; adios2::Dims count = {1, shape[1]}; - adios2::Variable var = stream.GetArrayVariable(varname, shape); + adios2::Variable var = + stream.GetArrayVariable(varname, shape, DIMS_NONE, BoutComm::rank()); var.SetSelection({start, count}); stream.engine.Put(var, value.begin()); } @@ -482,7 +492,8 @@ void ADIOSPutVarVisitor::operator()>(const Matrix& va (size_t)std::get<1>(s)}; adios2::Dims start = {(size_t)BoutComm::rank(), 0, 0}; adios2::Dims count = {1, shape[1], shape[2]}; - adios2::Variable var = stream.GetArrayVariable(varname, shape); + adios2::Variable var = + stream.GetArrayVariable(varname, shape, DIMS_NONE, BoutComm::rank()); var.SetSelection({start, count}); stream.engine.Put(var, value.begin()); } @@ -495,7 +506,8 @@ void ADIOSPutVarVisitor::operator()>(const Tensor& va (size_t)std::get<1>(s), (size_t)std::get<2>(s)}; adios2::Dims start = {(size_t)BoutComm::rank(), 0, 0, 0}; adios2::Dims count = {1, shape[1], shape[2], shape[3]}; - adios2::Variable var = stream.GetArrayVariable(varname, shape); + adios2::Variable var = + stream.GetArrayVariable(varname, shape, DIMS_NONE, BoutComm::rank()); var.SetSelection({start, count}); stream.engine.Put(var, value.begin()); } diff --git a/tests/MMS/laplace/data/BOUT.inp b/tests/MMS/laplace/data/BOUT.inp index 1965b1e8f7..5b3cc37b09 100644 --- a/tests/MMS/laplace/data/BOUT.inp +++ b/tests/MMS/laplace/data/BOUT.inp @@ -5,7 +5,7 @@ MXG = 1 solution = sin(pi*x) -input = -pi^2*sin(pi*x) +input_field = -pi^2*sin(pi*x) [mesh] diff --git a/tests/MMS/laplace/laplace.cxx b/tests/MMS/laplace/laplace.cxx index 475725f269..28fbfd903f 100644 --- a/tests/MMS/laplace/laplace.cxx +++ b/tests/MMS/laplace/laplace.cxx @@ -1,16 +1,21 @@ #include - #include +#include #include #include +#include + +#include using bout::globals::mesh; +using namespace std::string_literals; int main(int argc, char** argv) { int init_err = BoutInitialise(argc, argv); if (init_err < 0) { return 0; - } else if (init_err > 0) { + } + if (init_err > 0) { return init_err; } @@ -31,16 +36,14 @@ int main(int argc, char** argv) { FieldFactory fact(mesh); - std::shared_ptr gen = fact.parse("input"); - output << "GEN = " << gen->str() << endl; - - Field3D input = fact.create3D("input"); - - Field3D result = lap->solve(input); - - Field3D solution = fact.create3D("solution"); + const auto input_name = "input_field"s; + const auto gen = fact.parse(input_name); + output.write("GEN = {}\n", gen->str()); - Field3D error = result - solution; + const Field3D input = fact.create3D(input_name); + const Field3D result = lap->solve(input); + const Field3D solution = fact.create3D("solution"); + const Field3D error = result - solution; Options dump; dump["input"] = input; diff --git a/tests/MMS/spatial/d2dx2/data/BOUT.inp b/tests/MMS/spatial/d2dx2/data/BOUT.inp index 5c1deaf0af..cdc0d96139 100644 --- a/tests/MMS/spatial/d2dx2/data/BOUT.inp +++ b/tests/MMS/spatial/d2dx2/data/BOUT.inp @@ -5,7 +5,7 @@ MXG = 1 MYG = 1 # No guard cells in Y -input = sin(0.5*pi*x) +input_field = sin(0.5*pi*x) solution = -sin(0.5*pi*x) * 0.25*pi*pi [mesh] diff --git a/tests/MMS/spatial/d2dx2/test_d2dx2.cxx b/tests/MMS/spatial/d2dx2/test_d2dx2.cxx index 6b35b719e1..7f7127cea9 100644 --- a/tests/MMS/spatial/d2dx2/test_d2dx2.cxx +++ b/tests/MMS/spatial/d2dx2/test_d2dx2.cxx @@ -4,7 +4,9 @@ #include #include +#include #include +#include using bout::globals::mesh; @@ -12,8 +14,9 @@ int main(int argc, char** argv) { BoutInitialise(argc, argv); - Field3D input = FieldFactory::get()->create3D("input", Options::getRoot(), mesh); - Field3D solution = FieldFactory::get()->create3D("solution", Options::getRoot(), mesh); + Field3D input = FieldFactory::get()->create3D("input_field", Options::getRoot(), mesh); + const Field3D solution = + FieldFactory::get()->create3D("solution", Options::getRoot(), mesh); // At this point the boundary cells are set to the analytic solution input.setBoundary("bndry"); @@ -21,7 +24,7 @@ int main(int argc, char** argv) { // Boundaries of input now set using extrapolation around mid-point boundary - Field3D result = D2DX2(input); + const Field3D result = D2DX2(input); // At this point result is not set in the boundary cells Options dump; diff --git a/tests/MMS/spatial/d2dz2/data/BOUT.inp b/tests/MMS/spatial/d2dz2/data/BOUT.inp index 6618d6f574..640411c112 100644 --- a/tests/MMS/spatial/d2dz2/data/BOUT.inp +++ b/tests/MMS/spatial/d2dz2/data/BOUT.inp @@ -5,7 +5,7 @@ zperiod = 1 MXG = 0 # No guard cells in X MYG = 1 # No guard cells in Y -input = sin(z) +input_field = sin(z) solution = -sin(z) [mesh] diff --git a/tests/MMS/spatial/d2dz2/test_d2dz2.cxx b/tests/MMS/spatial/d2dz2/test_d2dz2.cxx index 0b7aa2cde4..68a12b062f 100644 --- a/tests/MMS/spatial/d2dz2/test_d2dz2.cxx +++ b/tests/MMS/spatial/d2dz2/test_d2dz2.cxx @@ -4,7 +4,9 @@ #include #include +#include #include +#include using bout::globals::mesh; @@ -12,8 +14,10 @@ int main(int argc, char** argv) { BoutInitialise(argc, argv); - Field3D input = FieldFactory::get()->create3D("input", Options::getRoot(), mesh); - Field3D solution = FieldFactory::get()->create3D("solution", Options::getRoot(), mesh); + const Field3D input = + FieldFactory::get()->create3D("input_field", Options::getRoot(), mesh); + const Field3D solution = + FieldFactory::get()->create3D("solution", Options::getRoot(), mesh); Field3D result = D2DZ2(input); diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index 47253c508f..e76a15dd57 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -67,6 +67,7 @@ set(serial_tests_source ./mesh/test_boundary_factory.cxx ./mesh/test_boutmesh.cxx ./mesh/test_coordinates.cxx + ./mesh/test_change_coordinate_system.cxx ./mesh/test_coordinates_accessor.cxx ./mesh/test_interpolation.cxx ./mesh/test_invert3x3.cxx diff --git a/tests/unit/fake_mesh_fixture.hxx b/tests/unit/fake_mesh_fixture.hxx index 4c33093aee..3c2215acba 100644 --- a/tests/unit/fake_mesh_fixture.hxx +++ b/tests/unit/fake_mesh_fixture.hxx @@ -25,121 +25,134 @@ /// alias to make a new test: /// /// using MyTest = FakeMeshFixture; -class FakeMeshFixture : public ::testing::Test { +/// +/// Type alias `FakeMeshFixture = FakeMeshFixture_tmpl<3, 5, 7>` +/// is used as a shim to allow FakeMeshFixture to be used with default values for nx, ny, nz. +/// Use this template class directly to use different sized grid: +/// +/// using MyTest = FakeMeshFixture_tmpl<7, 9, 11>; +template +class FakeMeshFixture_tmpl : public ::testing::Test { public: - FakeMeshFixture() { - WithQuietOutput quiet_info{output_info}; - WithQuietOutput quiet_warn{output_warn}; - - delete bout::globals::mesh; - bout::globals::mpi = new MpiWrapper(); - bout::globals::mesh = new FakeMesh(nx, ny, nz); - bout::globals::mesh->createDefaultRegions(); - static_cast(bout::globals::mesh)->setCoordinates(nullptr); - test_coords = std::make_shared( - bout::globals::mesh, Field2D{1.0}, Field2D{1.0}, Field2D{1.0}, Field2D{1.0}, - Field2D{1.0}, Field2D{1.0}, Field2D{1.0}, Field2D{1.0}, Field2D{0.0}, - Field2D{0.0}, Field2D{0.0}, Field2D{1.0}, Field2D{1.0}, Field2D{1.0}, - Field2D{0.0}, Field2D{0.0}, Field2D{0.0}, Field2D{0.0}, Field2D{0.0}); - static_cast(bout::globals::mesh)->setCoordinates(test_coords); - - // Set nonuniform corrections - test_coords->setNon_uniform(true); - test_coords->setD1_dx(0.2); - test_coords->setD1_dy(0.2); - test_coords->setD1_dz(0.0); + FakeMeshFixture_tmpl() + : mesh_m(NX, NY, NZ, mpi), mesh_staggered_m(NX, NY, NZ, mpi), + mesh_staggered(&mesh_staggered_m) { + + bout::globals::mpi = &mpi; + bout::globals::mesh = &mesh_m; + bout::globals::mesh->createDefaultRegions(); + static_cast(bout::globals::mesh)->setCoordinates(nullptr); + test_coords = std::make_shared( + bout::globals::mesh, Field2D{1.0}, Field2D{1.0}, Field2D{1.0}, Field2D{1.0}, + Field2D{1.0}, Field2D{1.0}, Field2D{1.0}, Field2D{1.0}, Field2D{0.0}, + Field2D{0.0}, Field2D{0.0}, Field2D{1.0}, Field2D{1.0}, Field2D{1.0}, + Field2D{0.0}, Field2D{0.0}, Field2D{0.0}, Field2D{0.0}, Field2D{0.0}); + + // Set nonuniform corrections + test_coords->setNon_uniform(true); + test_coords->setD1_dx(0.2); + test_coords->setD1_dy(0.2); + test_coords->setD1_dz(0.0); + #if BOUT_USE_METRIC_3D - FieldMetric mutable_Bxy = test_coords->Bxy(); - mutable_Bxy.splitParallelSlices(); - test_coords->setBxy(mutable_Bxy); + FieldMetric mutable_Bxy = test_coords->Bxy(); + mutable_Bxy.splitParallelSlices(); + test_coords->setBxy(mutable_Bxy); - mutable_Bxy = test_coords->Bxy(); - mutable_Bxy.yup() = test_coords->Bxy(); - mutable_Bxy.ydown() = test_coords->Bxy(); - test_coords->setBxy(mutable_Bxy); + mutable_Bxy = test_coords->Bxy(); + mutable_Bxy.yup() = test_coords->Bxy(); + mutable_Bxy.ydown() = test_coords->Bxy(); + test_coords->setBxy(mutable_Bxy); #endif - static_cast(bout::globals::mesh)->setCoordinates(test_coords); - static_cast(bout::globals::mesh) - ->setGridDataSource(new FakeGridDataSource()); - // May need a ParallelTransform to create fields, because create3D calls - // fromFieldAligned - test_coords->setParallelTransform( - bout::utils::make_unique(*bout::globals::mesh)); - dynamic_cast(bout::globals::mesh)->createBoundaryRegions(); - - delete mesh_staggered; - mesh_staggered = new FakeMesh(nx, ny, nz); - mesh_staggered->StaggerGrids = true; - dynamic_cast(mesh_staggered)->setCoordinates(nullptr); - dynamic_cast(mesh_staggered)->setCoordinates(nullptr, CELL_XLOW); - dynamic_cast(mesh_staggered)->setCoordinates(nullptr, CELL_YLOW); - dynamic_cast(mesh_staggered)->setCoordinates(nullptr, CELL_ZLOW); - mesh_staggered->createDefaultRegions(); - - test_coords_staggered = std::make_shared( - mesh_staggered, Field2D{1.0, mesh_staggered}, Field2D{1.0, mesh_staggered}, - Field2D{1.0, mesh_staggered}, Field2D{1.0, mesh_staggered}, - Field2D{1.0, mesh_staggered}, Field2D{1.0, mesh_staggered}, - Field2D{1.0, mesh_staggered}, Field2D{1.0, mesh_staggered}, - Field2D{0.0, mesh_staggered}, Field2D{0.0, mesh_staggered}, - Field2D{0.0, mesh_staggered}, Field2D{1.0, mesh_staggered}, - Field2D{1.0, mesh_staggered}, Field2D{1.0, mesh_staggered}, - Field2D{0.0, mesh_staggered}, Field2D{0.0, mesh_staggered}, - Field2D{0.0, mesh_staggered}, Field2D{0.0, mesh_staggered}, - Field2D{0.0, mesh_staggered}); - static_cast(mesh_staggered)->setCoordinates(test_coords_staggered); - - // Set nonuniform corrections - test_coords_staggered->setNon_uniform(true); - test_coords_staggered->setD1_dx(0.2); - test_coords_staggered->setD1_dy(0.2); - test_coords_staggered->setD1_dz(0.0); + mesh_m.setCoordinates(test_coords); + mesh_m.setGridDataSource(new FakeGridDataSource()); + // May need a ParallelTransform to create fields, because create3D calls + // fromFieldAligned + test_coords->setParallelTransform( + bout::utils::make_unique(*bout::globals::mesh)); + mesh_m.createBoundaryRegions(); + + mesh_staggered_m.StaggerGrids = true; + mesh_staggered_m.setCoordinates(nullptr); + mesh_staggered_m.setCoordinates(nullptr, CELL_XLOW); + mesh_staggered_m.setCoordinates(nullptr, CELL_YLOW); + mesh_staggered_m.setCoordinates(nullptr, CELL_ZLOW); + mesh_staggered_m.createDefaultRegions(); + + test_coords_staggered = std::make_shared( + &mesh_staggered_m, Field2D{1.0, &mesh_staggered_m}, + Field2D{1.0, &mesh_staggered_m}, Field2D{1.0, &mesh_staggered_m}, + Field2D{1.0, &mesh_staggered_m}, Field2D{1.0, &mesh_staggered_m}, + Field2D{1.0, &mesh_staggered_m}, Field2D{1.0, &mesh_staggered_m}, + Field2D{1.0, &mesh_staggered_m}, Field2D{0.0, &mesh_staggered_m}, + Field2D{0.0, &mesh_staggered_m}, Field2D{0.0, &mesh_staggered_m}, + Field2D{1.0, &mesh_staggered_m}, Field2D{1.0, &mesh_staggered_m}, + Field2D{1.0, &mesh_staggered_m}, Field2D{0.0, &mesh_staggered_m}, + Field2D{0.0, &mesh_staggered_m}, Field2D{0.0, &mesh_staggered_m}, + Field2D{0.0, &mesh_staggered_m}, Field2D{0.0, &mesh_staggered_m}); + + // Set nonuniform corrections + test_coords_staggered->setNon_uniform(true); + test_coords_staggered->setD1_dx(0.2); + test_coords_staggered->setD1_dy(0.2); + test_coords_staggered->setD1_dz(0.0); + #if BOUT_USE_METRIC_3D - mutable_Bxy = test_coords_staggered->Bxy(); - mutable_Bxy.splitParallelSlices(); - test_coords_staggered->setBxy(mutable_Bxy); + mutable_Bxy = test_coords_staggered->Bxy(); + mutable_Bxy.splitParallelSlices(); + test_coords_staggered->setBxy(mutable_Bxy); - mutable_Bxy = test_coords_staggered->Bxy(); - mutable_Bxy.yup() = test_coords_staggered->Bxy(); - mutable_Bxy.ydown() = test_coords_staggered->Bxy(); - test_coords_staggered->setBxy(mutable_Bxy); + mutable_Bxy = test_coords_staggered->Bxy(); + mutable_Bxy.yup() = test_coords_staggered->Bxy(); + mutable_Bxy.ydown() = test_coords_staggered->Bxy(); + test_coords_staggered->setBxy(mutable_Bxy); #endif - test_coords_staggered->setParallelTransform( - bout::utils::make_unique(*mesh_staggered)); - - // Set all coordinates to the same Coordinates object for now - dynamic_cast(mesh_staggered)->setCoordinates(test_coords_staggered); - dynamic_cast(mesh_staggered) - ->setCoordinates(test_coords_staggered, CELL_XLOW); - dynamic_cast(mesh_staggered) - ->setCoordinates(test_coords_staggered, CELL_YLOW); - dynamic_cast(mesh_staggered) - ->setCoordinates(test_coords_staggered, CELL_ZLOW); - } - - ~FakeMeshFixture() override { - delete bout::globals::mesh; - bout::globals::mesh = nullptr; - delete mesh_staggered; - mesh_staggered = nullptr; - delete bout::globals::mpi; - bout::globals::mpi = nullptr; - - Options::cleanup(); - } - - static constexpr int nx = 3; - static constexpr int ny = 5; - static constexpr int nz = 7; - - Mesh* mesh_staggered = nullptr; - - std::shared_ptr test_coords{nullptr}; - std::shared_ptr test_coords_staggered{nullptr}; -}; \ No newline at end of file + test_coords_staggered->setParallelTransform( + bout::utils::make_unique(mesh_staggered_m)); + + // Set all coordinates to the same Coordinates object for now + mesh_staggered_m.setCoordinates(test_coords_staggered); + mesh_staggered_m.setCoordinates(test_coords_staggered, CELL_XLOW); + mesh_staggered_m.setCoordinates(test_coords_staggered, CELL_YLOW); + mesh_staggered_m.setCoordinates(test_coords_staggered, CELL_ZLOW); + } + + FakeMeshFixture_tmpl(const FakeMeshFixture_tmpl&) = delete; + FakeMeshFixture_tmpl& operator=(const FakeMeshFixture_tmpl&) = delete; + FakeMeshFixture_tmpl(FakeMeshFixture_tmpl&&) = delete; + FakeMeshFixture_tmpl& operator=(FakeMeshFixture_tmpl&&) = delete; + + ~FakeMeshFixture_tmpl() override { + bout::globals::mesh = nullptr; + bout::globals::mpi = nullptr; + + Options::cleanup(); + } + + static constexpr int nx = NX; + static constexpr int ny = NY; + static constexpr int nz = NZ; + +private: + std::shared_ptr test_coords{nullptr}; + std::shared_ptr test_coords_staggered{nullptr}; + + WithQuietOutput quiet_info{output_info}; + WithQuietOutput quiet_warn{output_warn}; + MpiWrapper mpi; + + FakeMesh mesh_m; + FakeMesh mesh_staggered_m; + +public: + // Public pointer to our staggered mesh + Mesh* mesh_staggered; // NOLINT +}; + +using FakeMeshFixture = FakeMeshFixture_tmpl<3, 5, 7>; diff --git a/tests/unit/mesh/test_change_coordinate_system.cxx b/tests/unit/mesh/test_change_coordinate_system.cxx new file mode 100644 index 0000000000..d513007d65 --- /dev/null +++ b/tests/unit/mesh/test_change_coordinate_system.cxx @@ -0,0 +1,77 @@ +#include "gtest/gtest.h" +#include "bout/coordinates.hxx" +#include "bout/mesh.hxx" +#include "fake_mesh_fixture.hxx" +#include "bout/constants.hxx" +#include + +static constexpr int NX = 3; +static constexpr int NY = 5; +static constexpr int NZ = 8; + +using bout::globals::mesh; + +class CoordinateTransformTest : public FakeMeshFixture_tmpl { +public: + using FieldMetric = Coordinates::FieldMetric; + + CoordinateTransformTest() : FakeMeshFixture_tmpl() {} +}; + +TEST_F(CoordinateTransformTest, CylindricalToCartesian) { + + // arrange + + // Set up test values + // Calculate cylindrical coordinates (Rxy, Zxy) + // from (2D) orthogonal poloidal coordinates (r, theta) + + const double R0 = 2.0; // major radius + const std::array r_values = {0.1, 0.2, 0.3}; // minor radius + const std::array theta_values = { // poloidal angle + 0.0, + PI / 2, + PI, + 3 * PI / 2, + 2 * PI}; + + auto tokamak_options = bout::TokamakOptions(*mesh); + + int i = 0; + for (auto r: r_values) { + int j = 0; + for (auto theta: theta_values) { + tokamak_options.Rxy(i, j) = R0 + r * std::cos(theta); + tokamak_options.Zxy(i, j) = r * std::sin(theta); + j++; + } + i++; + } + + // act + bout::Coordinates3D cartesian_coords = tokamak_options.cylindrical_coordinates_to_cartesian(); + + // assert + const auto max_r = *std::max_element(begin(r_values), end(r_values)); + const auto expected_max_x = R0 + max_r; + const auto expected_max_y = (R0 + max_r); + const auto expected_max_z = max_r; + + const auto expected_min_x = -1 * (R0 + max_r); + const auto expected_min_y = -1 * (R0 + max_r); + const auto expected_min_z = -1 * expected_max_z; + + const auto actual_max_x = max(cartesian_coords.x, false, "RGN_ALL"); + const auto actual_max_y = max(cartesian_coords.y, false, "RGN_ALL"); + const auto actual_max_z = max(cartesian_coords.z, false, "RGN_ALL"); + const auto actual_min_x = min(cartesian_coords.x, false, "RGN_ALL"); + const auto actual_min_y = min(cartesian_coords.y, false, "RGN_ALL"); + const auto actual_min_z = min(cartesian_coords.z, false, "RGN_ALL"); + + EXPECT_EQ(expected_max_x, actual_max_x); + EXPECT_EQ(expected_max_y, actual_max_y); + EXPECT_EQ(expected_max_z, actual_max_z); + EXPECT_EQ(expected_min_x, actual_min_x); + EXPECT_EQ(expected_min_y, actual_min_y); + EXPECT_EQ(expected_min_z, actual_min_z); +}