Add conservative flux divergence operator (FV::ConservativeFluxDiv)

CMakeLists.txt

@@ -115,6 +115,7 @@ set(BOUT_SOURCES
   ./include/bout/field_factory.hxx
   ./include/bout/fieldgroup.hxx
   ./include/bout/fieldperp.hxx
+  ./include/bout/facefield3d.hxx
   ./include/bout/fv_ops.hxx
   ./include/bout/generic_factory.hxx
   ./include/bout/globalfield.hxx
@@ -139,6 +140,7 @@ set(BOUT_SOURCES
   ./include/bout/macro_for_each.hxx
   ./include/bout/mask.hxx
   ./include/bout/mesh.hxx
+  ./include/bout/mesh_facefield_comm.hxx
   ./include/bout/monitor.hxx
   ./include/bout/mpi_wrapper.hxx
   ./include/bout/msg_stack.hxx
@@ -199,6 +201,7 @@ set(BOUT_SOURCES
   ./src/field/fieldgenerators.hxx
   ./src/field/fieldgroup.cxx
   ./src/field/fieldperp.cxx
+  ./src/field/facefield3d.cxx
   ./src/field/generated_fieldops.cxx
   ./src/field/globalfield.cxx
   ./src/field/initialprofiles.cxx
@@ -256,6 +259,7 @@ set(BOUT_SOURCES
   ./src/mesh/boundary_standard.cxx
   ./src/mesh/coordinates.cxx
   ./src/mesh/coordinates_accessor.cxx
+  ./src/mesh/conservative_flux_div.cxx
   ./src/mesh/data/gridfromfile.cxx
   ./src/mesh/data/gridfromoptions.cxx
   ./src/mesh/difops.cxx
@@ -272,6 +276,7 @@ set(BOUT_SOURCES
   ./src/mesh/interpolation/monotonic_hermite_spline_xz.cxx
   ./src/mesh/invert3x3.hxx
   ./src/mesh/mesh.cxx
+  ./src/mesh/mesh_facefield_comm.cxx
   ./src/mesh/parallel/fci.cxx
   ./src/mesh/parallel/fci.hxx
   ./src/mesh/parallel/identity.cxx

include/bout/conservative_flux_div.hxx

@@ -0,0 +1,82 @@
+/*!
+ * \file conservative_flux_div.hxx
+ *
+ * \brief Conservative flux divergence operator for finite-volume methods
+ *
+ * \author BOUT++ Team
+ *
+ **************************************************************************
+ * Copyright 2024 BOUT++ Contributors
+ *
+ * Contact: Ben Dudson, bd512@york.ac.uk
+ * 
+ * This file is part of BOUT++.
+ *
+ * BOUT++ is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * BOUT++ is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with BOUT++.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#pragma once
+#ifndef BOUT_CONSERVATIVE_FLUX_DIV_H
+#define BOUT_CONSERVATIVE_FLUX_DIV_H
+
+#include "bout/facefield3d.hxx"
+#include "bout/field3d.hxx"
+#include "bout/mesh.hxx"
+#include "bout/coordinates.hxx"
+#include "bout/region.hxx"
+
+namespace FV {
+
+/*!
+ * \brief Compute conservative divergence of face-centered fluxes
+ * 
+ * This function implements a finite-volume divergence operator that
+ * ensures exact conservation. Given fluxes on cell faces, it computes
+ * the divergence using:
+ * 
+ * div(F) = 1/V * sum(F·n·dA)
+ * 
+ * where the sum is over all cell faces, V is the cell volume,
+ * n is the outward normal, and dA is the face area.
+ * 
+ * \param F Face-centered flux field with components on cell faces
+ * \param bndry_flux If true, use fluxes at domain boundaries (default: false)
+ * \return Cell-centered divergence field
+ * 
+ * \note This operator ensures exact conservation: the sum of div(F)*V
+ *       over all cells equals the net flux through domain boundaries
+ */
+Field3D ConservativeFluxDiv(const FaceField3D& F, bool bndry_flux = false);
+
+/*!
+ * \brief Compute conservative divergence with specified boundary flux
+ * 
+ * This variant allows specifying the flux values at domain boundaries
+ * rather than using zero flux (default) or extrapolation.
+ * 
+ * \param F Face-centered flux field
+ * \param bndry_flux_x Flux at X boundaries (optional)
+ * \param bndry_flux_y Flux at Y boundaries (optional) 
+ * \param bndry_flux_z Flux at Z boundaries (optional)
+ * \return Cell-centered divergence field
+ */
+Field3D ConservativeFluxDiv(const FaceField3D& F,
+                           const Field3D* bndry_flux_x,
+                           const Field3D* bndry_flux_y,
+                           const Field3D* bndry_flux_z);
+
+} // namespace FV
+
+#endif // BOUT_CONSERVATIVE_FLUX_DIV_H