boutproject · ZedThree · Aug 13, 2026 · May 11, 2026 · Jul 8, 2026 · Jul 8, 2026
diff --git a/CHANGELOG.md b/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)

diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -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
@@ -135,18 +136,21 @@ set(BOUT_SOURCES
     ./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
@@ -163,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
@@ -215,8 +220,8 @@ set(BOUT_SOURCES
     ./include/bout/vector2d.hxx
     ./include/bout/vector3d.hxx
     ./include/bout/where.hxx
-    ./src/bout++.cxx
     ./src/bout++-time.hxx
+    ./src/bout++.cxx
     ./src/field/field.cxx
     ./src/field/field2d.cxx
     ./src/field/field3d.cxx
@@ -239,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
@@ -261,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
@@ -285,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
@@ -368,18 +376,17 @@ 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/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
@@ -389,7 +396,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
 )

diff --git a/examples/6field-simple/elm_6f.cxx b/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

diff --git a/examples/IMEX/drift-wave-constraint/test-drift.cxx b/examples/IMEX/drift-wave-constraint/test-drift.cxx
@@ -61,10 +61,10 @@ class DriftWave : public PhysicsModel {
       // This code results in ddt(Ne) depending on y+2, y-2
       // which are not (currently) included in the coloring
       // The result is that IMEX-BDF2 with coloring doesn't converge
-      
+
     Ve = ( Grad_par(phi) - Grad_par(Ne) ) / nu;
     mesh->communicate(Ve);
-    
+
     ddt(Ne) = -Div_par(Ve);
     ddt(Vort) = -Div_par(Ve);
     */
@@ -82,7 +82,7 @@ class DriftWave : public PhysicsModel {
     // ddt(phi) = Delp2(phi) - Vort;
 
     // This version uses central differencing for Delp2
-    ddt(phi) = (coord->g11 * D2DX2(phi) + coord->g33 * D2DZ2(phi)) - Vort;
+    ddt(phi) = (coord->g11() * D2DX2(phi) + coord->g33() * D2DZ2(phi)) - Vort;
 
     return 0;
   }

diff --git a/examples/conducting-wall-mode/cwm.cxx b/examples/conducting-wall-mode/cwm.cxx
@@ -5,14 +5,26 @@
  * Model version in the code created by M. Umansky and J. Myra.
  *******************************************************************************/
 
+#include <bout/bout_types.hxx>
+#include <bout/coordinates.hxx>
 #include <bout/derivs.hxx>
+#include <bout/difops.hxx>
+#include <bout/field2d.hxx>
+#include <bout/field3d.hxx>
+#include <bout/globals.hxx>
 #include <bout/initialprofiles.hxx>
 #include <bout/interpolation.hxx>
 #include <bout/invert_laplace.hxx>
+#include <bout/output.hxx>
 #include <bout/physicsmodel.hxx>
+#include <bout/sys/range.hxx>
 #include <bout/tokamak_coordinates.hxx>
+#include <bout/unused.hxx>
 #include <bout/utils.hxx>
 
+#include <cmath>
+#include <memory>
+
 class CWM : public PhysicsModel {
 private:
   // 2D initial profiles
@@ -319,7 +331,7 @@ class CWM : public PhysicsModel {
       result = VDDX(DDZ(p), f);
     } else {
       // Use full expression with all terms
-      result = b0xGrad_dot_Grad(p, f) / coord->Bxy;
+      result = b0xGrad_dot_Grad(p, f) / coord->Bxy();
     }
     return result;
   }
@@ -328,10 +340,10 @@ class CWM : public PhysicsModel {
     Field3D result;
     if (bout_exb) {
       // Use a subset of terms for comparison to BOUT-06
-      result = VDDZ(-DDX(p), f);
+      result = VDDZ(Field2D{-DDX(p)}, f);
     } else {
       // Use full expression with all terms
-      result = b0xGrad_dot_Grad(p, f) / coord->Bxy;
+      result = b0xGrad_dot_Grad(p, f) / coord->Bxy();
     }
     return result;
   }
@@ -343,7 +355,7 @@ class CWM : public PhysicsModel {
       result = VDDX(DDZ(p), f) + VDDZ(-DDX(p), f);
     } else {
       // Use full expression with all terms
-      result = b0xGrad_dot_Grad(p, f) / coord->Bxy;
+      result = b0xGrad_dot_Grad(p, f) / coord->Bxy();
     }
     return result;
   }