Simplify multiphase test setup + add Coefficient VTK post-processing

kernel/include/Coefficients/MfemCoefficient.hpp

@@ -27,55 +27,114 @@
 #include <span>
 #include <vector>
 
+#include "Coefficients/Coefficient.hpp"
 #include "Coefficients/Coefficients.hpp"
 #include "mfem.hpp"  // NOLINT [no include the directory when naming mfem include file]
 
+/**
+ * @brief Wraps a SLOTH ::Coefficient as an mfem::Coefficient, so it can be
+ *        projected onto a grid function or plugged into an mfem
+ *        integrator.
+ */
 class MfemCoefficient : public mfem::Coefficient {
  private:
   int id_ = 0;
-  Coefficients sloth_coefficient_;
-  const mfem::ParGridFunction& gf_;
+  ::Coefficient sloth_coefficient_;
+
+  const std::vector<mfem::ParGridFunction>& gf_;
+  const std::vector<mfem::ParGridFunction>& gfn_;
   const std::vector<mfem::ParGridFunction>& vaux_gf_;
 
   mutable std::vector<double> vaux_gf_at_ip_;
+  mutable std::vector<double> gf_at_ip_;
+  mutable std::vector<double> gfn_at_ip_;
 
  public:
-  MfemCoefficient(int id, const Coefficients& coef, const mfem::ParGridFunction& vun,
+  MfemCoefficient(::Coefficient coef, const std::vector<mfem::ParGridFunction>& vu,
+                  const std::vector<mfem::ParGridFunction>& vun,
                   const std::vector<mfem::ParGridFunction>& vauxn);
   double Eval(mfem::ElementTransformation& T, const mfem::IntegrationPoint& ip) override;
+
+  mfem::ParGridFunction ProjectToGridFunction(mfem::ParFiniteElementSpace& fespace);
+  mfem::ParGridFunction ProjectToGridFunction();
+
   virtual ~MfemCoefficient() = default;
 };
 
 /**
- * @brief Construct a new Mfem Coefficient:: Mfem Coefficient object
+ * @brief Construct a new MfemCoefficient object.
  *
- * @param coef
- * @param vun
- * @param vauxn
+ * @param coef SLOTH coefficient to wrap.
+ * @param vu Current-time vector of grid function for the coefficient's own variable.
+ * @param vun Previous-time vector of grid function for the coefficient's own variable.
+ * @param vauxn Grid functions of the auxiliary variables.
  */
-MfemCoefficient::MfemCoefficient(int id, const Coefficients& coef, const mfem::ParGridFunction& vun,
+MfemCoefficient::MfemCoefficient(::Coefficient coef, const std::vector<mfem::ParGridFunction>& vu,
+                                 const std::vector<mfem::ParGridFunction>& vun,
                                  const std::vector<mfem::ParGridFunction>& vauxn)
-    : id_(id), sloth_coefficient_(coef), gf_(vun), vaux_gf_(vauxn) {}
+    : sloth_coefficient_(coef), gf_(vu), gfn_(vun), vaux_gf_(vauxn) {
+  MFEM_VERIFY(this->gf_.size() == this->gfn_.size(),
+              "Error of consistency between Current-time and Previous-time vectors of grid "
+              "function. Please check their sizes.");
+}
 
 /**
- * @brief Eval MfemCoefficient at the integration point
+ * @brief Evaluate the wrapped coefficient at an integration point.
+ *
+ * Dispatches on the coefficient's scheme (scalar/semi-implicit/
+ * explicit/implicit) to call the matching ::Coefficient::compute overload.
  *
- * @param T
- * @param ip
- * @return double
+ * @param T Element transformation at the integration point.
+ * @param ip Integration point.
+ * @return Coefficient value.
  */
 double MfemCoefficient::Eval(mfem::ElementTransformation& T, const mfem::IntegrationPoint& ip) {
-  if (this->sloth_coefficient_[id_].is_scalar()) {
-    return this->sloth_coefficient_[id_].compute();
+  if (this->sloth_coefficient_.is_scalar()) {
+    return this->sloth_coefficient_.compute();
   }
 
   this->vaux_gf_at_ip_.resize(this->vaux_gf_.size());
   for (size_t k = 0; k < this->vaux_gf_.size(); ++k) {
     this->vaux_gf_at_ip_[k] = this->vaux_gf_[k].GetValue(T, ip);
   }
 
-  double var_at_ip = this->gf_.GetValue(T, ip);
-  std::array<double, 1> var = {var_at_ip};
-  return this->sloth_coefficient_[id_].compute(std::span<const double>(var),
-                                               std::span<const double>(this->vaux_gf_at_ip_));
+  this->gf_at_ip_.resize(this->gf_.size());
+  this->gfn_at_ip_.resize(this->gfn_.size());
+  for (size_t k = 0; k < this->gf_.size(); ++k) {
+    this->gf_at_ip_[k] = this->gf_[k].GetValue(T, ip);
+    this->gfn_at_ip_[k] = this->gfn_[k].GetValue(T, ip);
+  }
+
+  if (this->sloth_coefficient_.is_semi_implicit()) {
+    return this->sloth_coefficient_.compute(std::span<const double>(this->gf_at_ip_),
+                                            std::span<const double>(this->gfn_at_ip_),
+                                            std::span<const double>(this->vaux_gf_at_ip_));
+  } else if (this->sloth_coefficient_.is_explicit()) {
+    return this->sloth_coefficient_.compute(std::span<const double>(this->gfn_at_ip_),
+                                            std::span<const double>(this->vaux_gf_at_ip_));
+  } else {
+    return this->sloth_coefficient_.compute(std::span<const double>(this->gf_at_ip_),
+                                            std::span<const double>(this->vaux_gf_at_ip_));
+  }
+}
+/**
+ * @brief Project the coefficient onto a grid function on a given space.
+ *
+ * @param fespace Finite element space to project onto.
+ * @return Grid function holding the projected values.
+ */
+mfem::ParGridFunction MfemCoefficient::ProjectToGridFunction(mfem::ParFiniteElementSpace& fespace) {
+  mfem::ParGridFunction result(&fespace);
+  result.ProjectCoefficient(*this);
+  return result;
+}
+
+/**
+ * @brief Project the coefficient onto a grid function, using the finite
+ *        element space of the first index in `gf_` (assumption of homogeneous FESpace)
+ *
+ * @return Grid function holding the projected values.
+ */
+mfem::ParGridFunction MfemCoefficient::ProjectToGridFunction() {
+  return ProjectToGridFunction(*gf_[0].ParFESpace());
 }

kernel/include/Coefficients/SlothBaseCoefficient.hpp

@@ -44,6 +44,7 @@
  */
 class SlothBaseCoefficient {
  private:
+  std::string name_;
   GlossaryQuantity coefficient_quantity_;
   std::shared_ptr<FunctionCoefficient> coefficient_;
   std::vector<int> bdr_index_;
@@ -116,6 +117,9 @@ class SlothBaseCoefficient {
   bool is_scalar() const;
 
   void set_time(double time);
+
+  std::string get_name();
+  void set_name(const std::string& name);
   void set_bdr_index_coef(std::vector<int> ids);
   std::vector<int> get_bdr_index_coef() const;
 };

kernel/include/Couplings/Coupling.hpp

@@ -69,6 +69,7 @@ class Coupling {
   void initialize(const int& iter, const double& initial_time, const double time_step,
                   bool vtk_unified, std::vector<VAR*> all_vars);
   void collect_vtk_fields(std::map<std::string, mfem::ParGridFunction*>& all_fields);
+  void collect_vtk_coefficients(std::map<std::string, mfem::ParGridFunction*>& all_coefficients);
   auto get_shared_dc() { return std::get<0>(this->problems_).get_shared_dc(); }
   void check_duplicate_post_processing_directory(std::size_t nb_couplings,
                                                  const std::map<std::string, int>& directory_count);

kernel/include/Couplings/Coupling.tpp

@@ -240,6 +240,23 @@ void Coupling<Args...>::collect_vtk_fields(
              this->problems_);
 }
 
+/**
+ * @brief Collect the grid functions of every Problem in this coupling
+ *        into a shared coefficient map.
+ *
+ * @tparam Args Problem types in this coupling.
+ * @param all_coefficients Output map, accumulated across every Problem, for a
+ *                   single unified VTK save (see
+ *                   `TimeDiscretization::save_vtk_unified()`).
+ */
+template <class... Args>
+void Coupling<Args...>::collect_vtk_coefficients(
+    std::map<std::string, mfem::ParGridFunction*>& all_coefficients) {
+  std::apply([&all_coefficients](
+                 auto&... problem) { (problem.collect_vtk_coefficients(all_coefficients), ...); },
+             this->problems_);
+}
+
 /**
  * @brief Solve all the problems inside the coupling
  *

kernel/include/Factory/BoundaryConditionsFactory.hpp

@@ -0,0 +1,61 @@
+/**
+ * @file BoundaryConditionsFactory.hpp
+ * @author Clément Introïni (clement.introini@cea.fr)
+ * @brief BoundaryConditions factory to manage boundary conditions
+ * @version 2.1
+ * @date 2026-09-05
+ *
+ * @anchor bcs
+ *
+ * @copyright CEA (C) 2026
+ *
+ * This file is part of SLOTH.
+ *
+ * SLOTH 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.
+ *
+ * SLOTH 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+#pragma once
+
+#include <cstddef>
+#include <vector>
+
+#include "BCs/BoundaryConditions.hpp"
+/**
+ * @brief Build N BoundaryConditions<T, DIM>, one per element of `spatials`,
+ *        all sharing the same `boundaries...` list.
+ *
+ * @tparam T Finite element collection type.
+ * @tparam DIM Spatial dimension.
+ * @tparam BoundaryArgs Deduced argument types forwarded to every
+ *                      BoundaryConditions<T, DIM> constructor call.
+ * @param N Number of BoundaryConditions to build. Must match spatials.size().
+ * @param spatials One SpatialDiscretization<T, DIM>* per grain (SPAS).
+ * @param boundaries Forwarded as-is to each BoundaryConditions constructor.
+ * @return std::vector<BoundaryConditions<T, DIM>> of size N.
+ */
+template <class T, int DIM, class... BoundaryArgs>
+std::vector<BoundaryConditions<T, DIM>> setBoundaryConditions(
+    std::size_t N, const std::vector<SpatialDiscretization<T, DIM>*>& spatials,
+    const BoundaryArgs&... boundaries) {
+  MFEM_VERIFY(N >= 1, "setBoundaryConditions: N must be >= 1");
+  MFEM_VERIFY(spatials.size() == N,
+              "setBoundaryConditions: spatials must contain exactly N elements");
+
+  std::vector<BoundaryConditions<T, DIM>> bcs;
+  bcs.reserve(N);
+  for (std::size_t i = 0; i < N; ++i) {
+    bcs.emplace_back(spatials[i], boundaries...);
+  }
+  return bcs;
+}

kernel/include/Factory/CouplingFactory.hpp

@@ -0,0 +1,75 @@
+/**
+ * @file CouplingFactory.hpp
+ * @brief Helper building Coupling<PB, PB, ..., PB> (N times) from
+ *        a std::vector<PB>, with N given explicitly as a template parameter.
+ * @version 0.1
+ *
+ * @anchor couplings
+ *
+ * @copyright CEA (C) 2026
+ *
+ * This file is part of SLOTH.
+ *
+ * SLOTH 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.
+ *
+ * SLOTH 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+#pragma once
+#include <cassert>
+#include <cstddef>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include "Couplings/Coupling.hpp"
+
+/**
+ * @brief Maps an index to PB, whatever the index -> used to repeat PB N times
+ *        in a template parameter pack (Coupling<RepeatedProblemType<PB,0>,
+ * RepeatedProblemType<PB,1>, ...> = Coupling<PB, PB, ...>).
+ *
+ * @tparam PB Problem type to repeat.
+ */
+template <class PB, std::size_t>
+using RepeatedProblemType = PB;
+
+/**
+ * @brief Implementation detail of HomogeneousCoupling: expands the index sequence
+ *        to build the Coupling<PB, PB, ..., PB> return type and forward each
+ *        vector element by index.
+ */
+template <class PB, std::size_t N, std::size_t... I>
+Coupling<RepeatedProblemType<PB, I>...> HomogeneousCoupling(const std::string& name,
+                                                            std::vector<PB> vect,
+                                                            std::index_sequence<I...>) {
+  return Coupling<RepeatedProblemType<PB, I>...>(name, std::move(vect[I])...);
+}
+
+/**
+ * @brief Build a Coupling<PB, PB, ..., PB> (N times) from a std::vector<PB>
+ *        whose size is only known at runtime.
+ *
+ *
+ * @tparam N Number of problems in `vect`. Must match `vect.size()` exactly
+ *           (checked with assert in debug builds).
+ * @tparam PB Problem type shared by every element of `vect`.
+ * @param name Name of the coupling.
+ * @param vect Problems to couple, all of type PB, size exactly N.
+ * @return Coupling<PB, PB, ..., PB> (N times), ready to pass to
+ *         TimeDiscretization.
+ */
+template <std::size_t N, class PB>
+auto setCoupling(const std::string& name, std::vector<PB> vect) {
+  MFEM_VERIFY(vect.size() == N, "setCoupling<N>: vect must contain exactly N elements");
+  return HomogeneousCoupling<PB, N>(name, std::move(vect), std::make_index_sequence<N>{});
+}

kernel/include/Factory/ListFactory.hpp

@@ -0,0 +1,31 @@
+/**
+ * @file ListFactory.hpp
+ * @author ci230846  (clement.introini@cea.fr)
+ * @brief List of Factory used to simplify tests setup
+ * @version 0.1
+ * @date 2026-09-05
+ *
+ * @copyright CEA (C) 2026
+ *
+ * This file is part of SLOTH.
+ *
+ * SLOTH 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.
+ *
+ * SLOTH 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#pragma once
+
+#include "Factory/BoundaryConditionsFactory.hpp"
+#include "Factory/CouplingFactory.hpp"
+#include "Factory/SpatialDiscretizationFactory.hpp"