merge next-elimBE-new-mri-filter with next-elimBE-new-mri

.clang-format

@@ -112,6 +112,6 @@ StatementMacros:
   - BOUT_OMP
   - BOUT_OMP_PERF
   - BOUT_OMP_SAFE
-Standard:        c++14
+Standard:        c++20
 TabWidth:        8
 UseTab:          Never

.github/workflows/auto-formatting.yml

@@ -16,7 +16,7 @@ jobs:
           fetch-depth: 0
 
       - name: "Set up Python"
-        uses: actions/setup-python@v6
+        uses: actions/setup-python@v7
 
       - name: Install uv
         uses: astral-sh/setup-uv@v7

.github/workflows/python-package.yml

@@ -22,7 +22,7 @@ jobs:
         fetch-depth: 0
         submodules: true
     - name: Set up Python
-      uses: actions/setup-python@v6
+      uses: actions/setup-python@v7
     - name: Install dependencies
       run: python -m pip install --upgrade pip &&
            pip install --upgrade build &&
@@ -62,7 +62,7 @@ jobs:
         fetch-depth: 0
         submodules: true
     - name: Set up Python
-      uses: actions/setup-python@v6
+      uses: actions/setup-python@v7
     - name: Install dependencies
       run: python -m pip install --upgrade pip &&
            pip install --upgrade build &&
@@ -111,7 +111,7 @@ jobs:
         fetch-depth: 0
         submodules: true
     - name: Set up Python
-      uses: actions/setup-python@v6
+      uses: actions/setup-python@v7
     - name: Install dependencies
       run: python -m pip install --upgrade pip &&
            pip install --upgrade build &&

.github/workflows/tests.yml

@@ -158,11 +158,11 @@ jobs:
           submodules: true
 
       - name: Install uv
-        uses: astral-sh/setup-uv@v5
+        uses: astral-sh/setup-uv@v7
         with:
           enable-cache: true
 
-      - uses: actions/setup-python@v6
+      - uses: actions/setup-python@v7
         with:
           python-version: '3.x'
 

.pre-commit-config.yaml

@@ -24,14 +24,14 @@ repos:
 
 # C++ formatting
 - repo: https://github.com/pre-commit/mirrors-clang-format
-  rev: v22.1.1
+  rev: v22.1.8
   hooks:
   - id: clang-format
     types_or: [c++, c, cuda]
 
 # Python linting and formatting
 - repo: https://github.com/astral-sh/ruff-pre-commit
-  rev: v0.15.6
+  rev: v0.16.2
   hooks:
     # Run the linter.
     - id: ruff-check

CHANGELOG.md

@@ -1,5 +1,21 @@
 # Changelog
 
+## [v6.0.0](https://github.com/boutproject/BOUT-dev/tree/v6.0.0
+
+[Full Changelog](https://github.com/boutproject/BOUT-dev/compare/v5.2.0...)
+
+### Breaking changes
+
+- `Coordinates` has been refactored to with a more encapsulated interface. This
+  allows us to ensure that derived quantities are correctly updated when the
+  metric components are changed, as well as defer creation of the Christoffel
+  symbols, and other derivatives of the metric components. This change will
+  require most physics models to be updated. You can do this with the
+  `bout-upgrader` from `boutdata`:
+  ```console
+  $ bout-upgrader v6 v6_upgrader file/to/update.cxx
+  ````
+
 ## [v5.2.0](https://github.com/boutproject/BOUT-dev/tree/v5.2.0
 
 [Full Changelog](https://github.com/boutproject/BOUT-dev/compare/v5.1.1...v5.2.0)

CMakeLists.txt

@@ -31,7 +31,7 @@ endif()
 # Set the version here, strip any extra tags to use in `project`
 # We try to use git to get a full description, inspired by setuptools_scm
 set(_bout_previous_version "5.2.0")
-set(_bout_next_version "5.2.1")
+set(_bout_next_version "6.0.0")
 execute_process(
   COMMAND "git" describe --tags --match=v${_bout_previous_version}
   COMMAND sed -e s/${_bout_previous_version}-/${_bout_next_version}.dev/ -e
@@ -116,10 +116,11 @@ set(BOUT_SOURCES
     ./include/bout/bout.hxx
     ./include/bout/bout_enum_class.hxx
     ./include/bout/bout_types.hxx
-    ./include/bout/build_config.hxx
     ./include/bout/boutcomm.hxx
     ./include/bout/boutexception.hxx
+    ./include/bout/build_config.hxx
     ./include/bout/caliper_wrapper.hxx
+    ./include/bout/christoffel_symbols.hxx
     ./include/bout/constants.hxx
     ./include/bout/coordinates.hxx
     ./include/bout/coordinates_accessor.hxx
@@ -128,26 +129,28 @@ set(BOUT_SOURCES
     ./include/bout/deriv_store.hxx
     ./include/bout/derivs.hxx
     ./include/bout/difops.hxx
-    ./include/bout/expr.hxx
     ./include/bout/fft.hxx
     ./include/bout/field.hxx
     ./include/bout/field2d.hxx
     ./include/bout/field3d.hxx
     ./include/bout/field_accessor.hxx
     ./include/bout/field_data.hxx
     ./include/bout/field_factory.hxx
-    ./include/bout/fieldops.hxx
     ./include/bout/fieldgroup.hxx
+    ./include/bout/fieldops.hxx
     ./include/bout/fieldperp.hxx
     ./include/bout/fv_ops.hxx
     ./include/bout/fv_ops_impl.hxx
+    ./include/bout/g_values.hxx
     ./include/bout/generic_factory.hxx
+    ./include/bout/git_metadata.hxx
     ./include/bout/globalfield.hxx
     ./include/bout/globalindexer.hxx
     ./include/bout/globals.hxx
     ./include/bout/griddata.hxx
     ./include/bout/gyro_average.hxx
     ./include/bout/hypre_interface.hxx
+    ./include/bout/hyprelib.hxx
     ./include/bout/index_derivs.hxx
     ./include/bout/index_derivs_interface.hxx
     ./include/bout/initialprofiles.hxx
@@ -164,6 +167,7 @@ set(BOUT_SOURCES
     ./include/bout/macro_for_each.hxx
     ./include/bout/mask.hxx
     ./include/bout/mesh.hxx
+    ./include/bout/metric_tensor.hxx
     ./include/bout/monitor.hxx
     ./include/bout/mpi_wrapper.hxx
     ./include/bout/msg_stack.hxx
@@ -176,7 +180,6 @@ set(BOUT_SOURCES
     ./include/bout/output.hxx
     ./include/bout/output_bout_types.hxx
     ./include/bout/parallel_boundary_op.hxx
-    ./include/bout/parallel_boundary_region.hxx
     ./include/bout/paralleltransform.hxx
     ./include/bout/petsc_interface.hxx
     ./include/bout/petsc_operators.hxx
@@ -216,8 +219,9 @@ set(BOUT_SOURCES
     ./include/bout/vector2d.hxx
     ./include/bout/vector3d.hxx
     ./include/bout/where.hxx
-    ./src/bout++.cxx
+    ./include/bout/yboundary_regions.hxx
     ./src/bout++-time.hxx
+    ./src/bout++.cxx
     ./src/field/field.cxx
     ./src/field/field2d.cxx
     ./src/field/field3d.cxx
@@ -240,6 +244,8 @@ set(BOUT_SOURCES
     ./src/invert/laplace/common_transform.hxx
     ./src/invert/laplace/impls/cyclic/cyclic_laplace.cxx
     ./src/invert/laplace/impls/cyclic/cyclic_laplace.hxx
+    ./src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx
+    ./src/invert/laplace/impls/hypre3d/hypre3d_laplace.hxx
     ./src/invert/laplace/impls/iterative_parallel_tri/iterative_parallel_tri.cxx
     ./src/invert/laplace/impls/iterative_parallel_tri/iterative_parallel_tri.hxx
     ./src/invert/laplace/impls/multigrid/multigrid_alg.cxx
@@ -262,8 +268,6 @@ set(BOUT_SOURCES
     ./src/invert/laplace/impls/serial_tri/serial_tri.hxx
     ./src/invert/laplace/impls/spt/spt.cxx
     ./src/invert/laplace/impls/spt/spt.hxx
-    ./src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx
-    ./src/invert/laplace/impls/hypre3d/hypre3d_laplace.hxx
     ./src/invert/laplace/invert_laplace.cxx
     ./src/invert/laplacexy/impls/hypre/laplacexy-hypre.cxx
     ./src/invert/laplacexy/impls/hypre/laplacexy-hypre.hxx
@@ -286,23 +290,26 @@ set(BOUT_SOURCES
     ./src/mesh/boundary_factory.cxx
     ./src/mesh/boundary_region.cxx
     ./src/mesh/boundary_standard.cxx
+    ./src/mesh/christoffel_symbols.cxx
     ./src/mesh/coordinates.cxx
     ./src/mesh/coordinates_accessor.cxx
     ./src/mesh/data/gridfromfile.cxx
     ./src/mesh/data/gridfromoptions.cxx
     ./src/mesh/difops.cxx
     ./src/mesh/fv_ops.cxx
+    ./src/mesh/g_values.cxx
     ./src/mesh/impls/bout/boutmesh.cxx
     ./src/mesh/impls/bout/boutmesh.hxx
     ./src/mesh/index_derivs.cxx
-    ./src/mesh/interpolation_xz.cxx
     ./src/mesh/interpolation/bilinear_xz.cxx
     ./src/mesh/interpolation/hermite_spline_xz.cxx
     ./src/mesh/interpolation/hermite_spline_z.cxx
     ./src/mesh/interpolation/interpolation_z.cxx
     ./src/mesh/interpolation/lagrange_4pt_xz.cxx
+    ./src/mesh/interpolation_xz.cxx
     ./src/mesh/invert3x3.hxx
     ./src/mesh/mesh.cxx
+    ./src/mesh/metric_tensor.cxx
     ./src/mesh/parallel/fci.cxx
     ./src/mesh/parallel/fci.hxx
     ./src/mesh/parallel/fci_comm.cxx
@@ -316,6 +323,7 @@ set(BOUT_SOURCES
     ./src/mesh/petsc_operators.cxx
     ./src/mesh/surfaceiter.cxx
     ./src/mesh/tokamak_coordinates.cxx
+    ./src/mesh/yboundary_regions.cxx
     ./src/physics/gyro_average.cxx
     ./src/physics/physicsmodel.cxx
     ./src/physics/smoothing.cxx
@@ -371,19 +379,18 @@ set(BOUT_SOURCES
     ./src/sys/derivs.cxx
     ./src/sys/expressionparser.cxx
     ./src/sys/generator_context.cxx
-    ./include/bout/hyprelib.hxx
     ./src/sys/hyprelib.cxx
     ./src/sys/hypre_interface.cxx
     ./src/sys/msg_stack.cxx
     ./src/sys/options.cxx
     ./src/sys/options/optionparser.hxx
+    ./src/sys/options/options_adios.cxx
+    ./src/sys/options/options_adios.hxx
     ./src/sys/options/options_ini.cxx
     ./src/sys/options/options_ini.hxx
     ./src/sys/options/options_io.cxx
     ./src/sys/options/options_netcdf.cxx
     ./src/sys/options/options_netcdf.hxx
-    ./src/sys/options/options_adios.cxx
-    ./src/sys/options/options_adios.hxx
     ./src/sys/optionsreader.cxx
     ./src/sys/output.cxx
     ./src/sys/output_bout_types.cxx
@@ -393,7 +400,6 @@ set(BOUT_SOURCES
     ./src/sys/timer.cxx
     ./src/sys/type_name.cxx
     ./src/sys/utils.cxx
-    ./include/bout/git_metadata.hxx
     ${CMAKE_CURRENT_BINARY_DIR}/include/bout/revision.hxx
     ${CMAKE_CURRENT_BINARY_DIR}/include/bout/version.hxx
 )

examples/6field-simple/elm_6f.cxx

@@ -6,19 +6,32 @@
  * T. Xia
  *******************************************************************************/
 
+#include "bout/assert.hxx"
 #include "bout/bout.hxx"
+#include "bout/bout_types.hxx"
+#include "bout/build_defines.hxx"
 #include "bout/constants.hxx"
+#include "bout/coordinates.hxx"
 #include "bout/derivs.hxx"
+#include "bout/difops.hxx"
+#include "bout/field2d.hxx"
+#include "bout/field3d.hxx"
+#include "bout/fieldgroup.hxx"
 #include "bout/initialprofiles.hxx"
 #include "bout/interpolation_xz.hxx"
 #include "bout/invert_laplace.hxx"
-#include "bout/invert_parderiv.hxx"
 #include "bout/msg_stack.hxx"
+#include "bout/output.hxx"
 #include "bout/physicsmodel.hxx"
+#include "bout/smoothing.hxx"
 #include "bout/sourcex.hxx"
 #include "bout/tokamak_coordinates.hxx"
+#include "bout/unused.hxx"
+#include "bout/vecops.hxx"
+#include "bout/where.hxx"
 
 #include <cmath>
+#include <memory>
 
 constexpr BoutReal eV_K = 11605.0; // 1eV = 11605K
 
@@ -356,7 +369,7 @@ class Elm_6f : public PhysicsModel {
       result.allocate();
       for (auto i : result) {
         result[i] =
-            (fp[i.yp()] - fm[i.ym()]) / (2. * coord->dy[i] * sqrt(coord->g_22[i]));
+            (fp[i.yp()] - fm[i.ym()]) / (2. * coord->dy()[i] * sqrt(coord->g_22()[i]));
       }
     } else {
       result = Grad_par(f, loc);
@@ -696,7 +709,7 @@ class Elm_6f : public PhysicsModel {
 
     if (mesh->IncIntShear) {
       // BOUT-06 style, using d/dx = d/dpsi + I * d/dz
-      mesh->getCoordinates()->IntShiftTorsion = I;
+      mesh->getCoordinates()->setIntShiftTorsion(I);
     } else {
       // Dimits style, using local coordinate system
       if (include_curvature) {
@@ -929,7 +942,7 @@ class Elm_6f : public PhysicsModel {
     output.write("\tlog Lambda: {:e}\n", LnLambda);
 
     nu_e = 2.91e-6 * LnLambda * ((N0)*Nbar * density / 1.e6)
-           * pow(Te0 * Tebar, -1.5); // nu_e in 1/S.
+           * pow(Field2D{Te0 * Tebar}, -1.5); // nu_e in 1/S.
     output.write("\telectron collision rate: {:e} -> {:e} [1/s]\n", min(nu_e), max(nu_e));
     // nu_e.applyBoundary();
     // mesh->communicate(nu_e);
@@ -941,7 +954,8 @@ class Elm_6f : public PhysicsModel {
       // xqx addition, begin
       // Use Spitzer thermal conductivities
       nu_i = 4.80e-8 * (Zi * Zi * Zi * Zi / sqrt(AA)) * LnLambda
-             * ((N0)*Nbar * density / 1.e6) * pow(Ti0 * Tibar, -1.5); // nu_i in 1/S.
+             * ((N0)*Nbar * density / 1.e6)
+             * pow(Field2D{Ti0 * Tibar}, -1.5); // nu_i in 1/S.
       // output.write("\tCoulomb Logarithm: {:e} \n", max(LnLambda));
       output.write("\tion collision rate: {:e} -> {:e} [1/s]\n", min(nu_i), max(nu_i));
 
@@ -1006,8 +1020,9 @@ class Elm_6f : public PhysicsModel {
       // Use Spitzer resistivity
       output.write("\n\tSpizter parameters");
       // output.write("\tTemperature: {:e} -> {:e} [eV]\n", min(Te), max(Te));
-      eta_spitzer = 0.51 * 1.03e-4 * Zi * LnLambda
-                    * pow(Te0 * Tebar, -1.5); // eta in Ohm-m. NOTE: ln(Lambda) = 20
+      eta_spitzer =
+          0.51 * 1.03e-4 * Zi * LnLambda
+          * pow(Field2D{Te0 * Tebar}, -1.5); // eta in Ohm-m. NOTE: ln(Lambda) = 20
       output.write("\tSpitzer resistivity: {:e} -> {:e} [Ohm m]\n", min(eta_spitzer),
                    max(eta_spitzer));
       eta_spitzer /= SI::mu0 * Va * Lbar;
@@ -1114,7 +1129,7 @@ class Elm_6f : public PhysicsModel {
       // Only if not restarting: Check initial perturbation
 
       // Set U to zero where P0 < vacuum_pressure
-      U = where(P0 - vacuum_pressure, U, 0.0);
+      U = where(Field2D{P0 - vacuum_pressure}, U, 0.0);
 
       //    Field2D lap_temp = 0.0;
       Field2D logn0 = laplace_alpha * N0;
@@ -1241,24 +1256,25 @@ class Elm_6f : public PhysicsModel {
       // Update resistivity
       if (spitzer_resist) {
         // Use Spitzer formula
-        eta_spitzer = 0.51 * 1.03e-4 * Zi * LnLambda
-                      * pow(Te_tmp * Tebar, -1.5); // eta in Ohm-m. ln(Lambda) = 20
+        eta_spitzer =
+            0.51 * 1.03e-4 * Zi * LnLambda
+            * pow(Field3D{Te_tmp * Tebar}, -1.5); // eta in Ohm-m. ln(Lambda) = 20
         eta_spitzer /= SI::mu0 * Va * Lbar;
       } else {
         eta = core_resist + (vac_resist - core_resist) * vac_mask;
       }
 
       nu_e = 2.91e-6 * LnLambda * (N_tmp * Nbar * density / 1.e6)
-             * pow(Te_tmp * Tebar, -1.5); // nu_e in 1/S.
+             * pow(Field3D{Te_tmp * Tebar}, -1.5); // nu_e in 1/S.
 
       if (diffusion_par > 0.0) {
         // Use Spitzer thermal conductivities
 
         nu_i = 4.80e-8 * (Zi * Zi * Zi * Zi / sqrt(AA)) * LnLambda
                * (N_tmp * Nbar * density / 1.e6)
-               * pow(Ti_tmp * Tibar, -1.5);         // nu_i in 1/S.
-        vth_i = 9.79e3 * sqrt(Ti_tmp * Tibar / AA); // vth_i in m/S.
-        vth_e = 4.19e5 * sqrt(Te_tmp * Tebar);      // vth_e in m/S.
+               * pow(Field3D{Ti_tmp * Tibar}, -1.5); // nu_i in 1/S.
+        vth_i = 9.79e3 * sqrt(Ti_tmp * Tibar / AA);  // vth_i in m/S.
+        vth_e = 4.19e5 * sqrt(Te_tmp * Tebar);       // vth_e in m/S.
       }
 
       if (diffusion_par > 0.0) {
@@ -1399,11 +1415,11 @@ class Elm_6f : public PhysicsModel {
       if (hyperviscos > 0.0) {
         // Calculate coefficient.
 
-        hyper_mu_x = hyperviscos * coord->g_11 * SQ(coord->dx)
-                     * abs(coord->g11 * D2DX2(U)) / (abs(U) + 1e-3);
+        hyper_mu_x = hyperviscos * coord->g_11() * SQ(coord->dx())
+                     * abs(coord->g11() * D2DX2(U)) / (abs(U) + 1e-3);
         hyper_mu_x.applyBoundary("dirichlet"); // Set to zero on all boundaries
 
-        ddt(U) += hyper_mu_x * coord->g11 * D2DX2(U);
+        ddt(U) += hyper_mu_x * coord->g11() * D2DX2(U);
 
         if (first_run) {
           // Print out maximum values of viscosity used on this processor