Defining the time-step as a function of the current simulation time

kernel/include/Time/Time.hpp

@@ -66,9 +66,13 @@ class TimeDiscretization {
   void time_management();
   void time_info(const int& iter);
 
+  std::function<double(double)> time_step_function_;
+
  public:
   // explicit TimeDiscretization(const Parameters& params, Args&&... couplings);
   explicit TimeDiscretization(const Parameters& params, Args... couplings);
+  explicit TimeDiscretization(const std::function<double(double)>& given_time_step,
+                              const Parameters& params, Args... couplings);
   void solve();
   void get_tree();
   ~TimeDiscretization();

kernel/include/Time/Time.tpp

@@ -38,7 +38,7 @@
 #include "mfem.hpp"  // NOLINT [no include the directory when naming mfem include file]
 
 /**
- * @brief Construct a new Time Discretization< Args...>:: Time Discretization object
+ * @brief Construct a new TimeDiscretization< Args...>::TimeDiscretization object
  *
  * @tparam Args
  * @param params
@@ -48,6 +48,28 @@ template <class... Args>
 TimeDiscretization<Args...>::TimeDiscretization(const Parameters& params, Args... couplings)
     : params_(params), couplings_(std::make_tuple(std::forward<Args>(couplings)...)) {
   this->get_parameters();
+
+  MFEM_VERIFY(this->time_step_ > 0., "Error: time_step Parameter must be defined.");
+
+  this->time_step_function_ =
+      std::function<double(double)>([this](double) { return this->time_step_; });
+}
+
+/**
+ * @brief Construct a new TimeDiscretization< Args...>::TimeDiscretization object
+ *
+ * @tparam Args
+ * @param given_time_step
+ * @param params
+ * @param couplings
+ */
+template <class... Args>
+TimeDiscretization<Args...>::TimeDiscretization(
+    const std::function<double(double)>& given_time_step, const Parameters& params,
+    Args... couplings)
+    : params_(params), couplings_(std::make_tuple(std::forward<Args>(couplings)...)) {
+  this->get_parameters();
+  this->time_step_function_ = given_time_step;
 }
 
 /**
@@ -61,8 +83,8 @@ void TimeDiscretization<Args...>::get_parameters() {
   this->initial_time_ =
       this->params_.template get_param_value_or_default<double>("initial_time", 0.);
   this->final_time_ = this->params_.template get_param_value<double>("final_time");
-  this->time_step_ = this->params_.template get_param_value<double>("time_step");
-  this->current_time_step_ = this->initial_time_;
+  this->time_step_ = this->params_.template get_param_value_or_default<double>("time_step", -1.0);
+  this->current_time_step_ = this->time_step_;
 }
 
 /**
@@ -194,7 +216,8 @@ void TimeDiscretization<Args...>::update() {
  */
 template <class... Args>
 void TimeDiscretization<Args...>::time_management() {
-  auto dt_real = this->time_step_;
+  auto dt_real = this->time_step_function_(this->current_time_);
+  this->time_step_ = dt_real;
   if ((this->current_time_ + this->time_step_) + 0.5 * this->time_step_ > this->final_time_) {
     this->last_step_ = true;
     dt_real = this->final_time_ - this->current_time_;

tests/CahnHilliard/2D/CMakeLists.txt

@@ -1,4 +1,5 @@
 add_subdirectory(test1)
+add_subdirectory(test1_dt_variable)
 add_subdirectory(test2)
 add_subdirectory(test3)
 add_subdirectory(test4)

tests/CahnHilliard/2D/test1_dt_variable/CMakeLists.txt

@@ -0,0 +1,3 @@
+
+create_test("CahnHilliard2Dtest1_dt" "CahnHilliard2Dtest1_dt" FALSE "2D EXE CH" 1)
+create_col_comparison("CompareCahnHilliard2Dtest1_dt" "time_specialized.csv" "Saves/CahnHilliard/time_specialized.csv" -1 relativeEpsilon 1e-12 FALSE "CahnHilliard2Dtest1_dt" "2D CH")

tests/CahnHilliard/2D/test1_dt_variable/main.cpp

@@ -0,0 +1,218 @@
+/**
+ * @file main.cpp
+ * @author ci230846  (clement.introini@cea.fr)
+ * @brief 2D coalescence bubbles solved by Cahn-Hilliard equations
+ * @version 0.1
+ * @date 2025-07-04
+ *
+ * Copyright CEA (c) 2025
+ *
+ */
+#include <iostream>
+#include <map>
+#include <memory>
+#include <sstream>
+#include <string>
+#include <tuple>
+#include <vector>
+
+#include "Sloth/sloth.hpp"
+#include "Sloth/tests.hpp"
+
+///---------------
+/// Main program
+///---------------
+int main(int argc, char* argv[]) {
+  //---------------------------------------
+  // Initialize MPI and HYPRE
+  //---------------------------------------
+  setVerbosity(Verbosity::Verbose);
+
+  mfem::Mpi::Init(argc, argv);
+  mfem::Hypre::Init();
+  //
+  //---------------------------------------
+  // Profiling start
+  Profiling::getInstance().enable();
+  //---------------------------------------
+  /////////////////////////
+  const int DIM = 2;
+  using FECollection = Test<DIM>::FECollection;
+  using VARS = Test<DIM>::VARS;
+  using VAR = Test<DIM>::VAR;
+  using PST = Test<DIM>::PST;
+  using SPA = Test<DIM>::SPA;
+  using BCS = Test<DIM>::BCS;
+  /////////////////////////
+  using OPE = TransientOperator<FECollection, DIM>;
+  using PB = Problem<OPE, VARS, PST>;
+  // ###########################################
+  // ###########################################
+  //         Spatial Discretization           //
+  // ###########################################
+  // ###########################################
+  // ##############################
+  //           Meshing           //
+  // ##############################
+  const std::string mesh_type =
+      "InlineSquareWithQuadrangles";  // type of mesh // "InlineSquareWithTriangles"
+  const int order_fe = 1;             // finite element order
+
+  const int refinement_level = 0;  // number of levels of uniform refinement
+  const int nx = 128;
+  const int ny = 128;
+  const double lx = 2. * M_PI;
+  const double ly = 2. * M_PI;
+  const std::tuple<int, int, double, double>& tuple_of_dimensions =
+      std::make_tuple(nx, ny, lx, ly);  // Number of elements and maximum length in each direction
+
+  SPA spatial(mesh_type, order_fe, refinement_level, tuple_of_dimensions);
+  // ##############################
+  //     Boundary conditions     //
+  // ##############################
+  auto boundaries = {Boundary("lower", 0, "Neumann"), Boundary("right", 1, "Neumann"),
+                     Boundary("upper", 2, "Neumann"), Boundary("left", 3, "Neumann")};
+  auto bcs_phi = BCS(&spatial, boundaries);
+  auto boundaries_mu = {Boundary("lower", 0, "Neumann"), Boundary("right", 1, "Neumann"),
+                        Boundary("upper", 2, "Neumann"), Boundary("left", 3, "Neumann")};
+  auto bcs_mu = BCS(&spatial, boundaries_mu);
+
+  // ###########################################
+  // ###########################################
+  //            Physical models               //
+  // ###########################################
+  // ###########################################
+  // ####################
+  //     parameters    //
+  // ####################
+  //  Interface thickness
+  const double epsilon(0.02);
+  // Interfacial energy
+  const double sigma(1.);
+  // Two-phase mobility
+  const double mob(1.);
+  const double lambda = (epsilon * epsilon);
+  const double omega = 1.;
+  auto params = Parameters(Parameter("epsilon", epsilon), Parameter("sigma", sigma),
+                           Parameter("lambda", lambda), Parameter("omega", omega));
+  // ####################
+  //     coefficients  //
+  // ####################
+
+  Coefficient grad_energy(Glossary::GradEnergy, Scheme::Implicit, GradientEnergy(lambda));
+  Coefficient double_well(Glossary::FreeEnergy, Scheme::Implicit, Fw(omega));
+  Coefficient capillary(Glossary::Capillary, lambda);
+  Coefficient mobility(Glossary::Mobility, mob);
+  // ####################
+  //     variables     //
+  // ####################
+
+  auto user_func_solution = std::function<double(const mfem::Vector&, double)>(
+      [](const mfem::Vector& x, [[maybe_unused]] double time) {
+        const double xx = x[0];
+        const double yy = x[1];
+        const double r1 = (xx - M_PI + 1) * (xx - M_PI + 1) + (yy - M_PI) * (yy - M_PI);
+        const double r2 = (xx - M_PI - 1) * (xx - M_PI - 1) + (yy - M_PI) * (yy - M_PI);
+        double sol = 0.;
+        if (r1 < 1 || r2 < 1) {
+          sol = 1.;
+        } else {
+          sol = -1.;
+        }
+        return sol;
+      });
+
+  auto phi_initial_condition = AnalyticalFunctions<DIM>(user_func_solution);
+  double mu_initial_condition = 0.0;
+  const std::string& var_name_1 = "phi";
+  const std::string& var_name_2 = "mu";
+  auto v1 = VAR(&spatial, bcs_phi, var_name_1, Glossary::PhaseField, 2, phi_initial_condition);
+  auto v2 = VAR(&spatial, bcs_mu, var_name_2, Glossary::ChemicalPotential, 2, mu_initial_condition);
+  auto vars = VARS(v1, v2);
+
+  // ###########################################
+  // ###########################################
+  //      Post-processing                     //
+  // ###########################################
+  // ###########################################
+
+  const std::string& main_folder_path = "Saves";
+  const int level_of_detail = 1;
+
+  const std::vector<int> iterations_list = {1, 3, 5};
+  const std::vector<double> times_list = {0.35, 0.45};
+
+  std::string calculation_path = "CahnHilliard";
+  std::map<std::string, std::tuple<double, double>> map_threshold_integral = {
+      {var_name_1, {-1.1, 1.1}}};
+  bool enable_save_specialized_at_iter = true;
+  auto p_pst =
+      Parameters(Parameter("main_folder_path", main_folder_path),
+                 Parameter("calculation_path", calculation_path),
+                 Parameter("iterations_list", iterations_list), Parameter("times_list", times_list),
+
+                 Parameter("level_of_detail", level_of_detail),
+                 Parameter("integral_to_compute", map_threshold_integral),
+                 Parameter("enable_save_specialized_at_iter", enable_save_specialized_at_iter));
+  // ####################
+  //     operators     //
+  // ####################
+
+  // Problem 1:
+  Coefficients coef_pb1(double_well, capillary, mobility, grad_energy);
+  std::vector<SPA*> spatials{&spatial, &spatial};
+  OPE oper(spatials, {"CahnHilliard"}, params, TimeScheme::EulerImplicit, "SplitTimeDerivative");
+  oper.overload_nl_solver(
+      NLSolverType::NEWTON,
+      Parameters(Parameter("description", "Newton solver "), Parameter("print_level", -1),
+                 Parameter("rel_tol", 1.e-10), Parameter("abs_tol", 1.e-14)));
+  const auto& solver = HypreSolverType::HYPRE_GMRES;
+  const auto& precond = HyprePreconditionerType::HYPRE_ILU;
+  oper.overload_solver(solver);
+  oper.overload_preconditioner(precond);
+
+  auto pst = PST(&spatial, p_pst);
+  PB problem1(oper, vars, {coef_pb1, coef_pb1}, pst);
+
+  // Coupling 1
+  auto cc = Coupling("CahnHilliard Coupling", problem1);
+
+  // ###########################################
+  // ###########################################
+  //            Time-integration              //
+  // ###########################################
+  // ###########################################
+  const double t_initial = 0.0;
+  const double t_final = 0.5;
+
+  auto user_time_step = std::function<double(double)>([](double time) {
+    double dt;
+    if (time < 0.1) {
+      dt = 0.01;
+    } else if (time < 0.2) {
+      dt = 0.02;
+    } else if (time < 0.4) {
+      dt = 0.04;
+    } else {
+      dt = 0.05;
+    }
+    return dt;
+  });
+
+  auto time_params =
+      Parameters(Parameter("initial_time", t_initial), Parameter("final_time", t_final));
+  auto time = TimeDiscretization(user_time_step, time_params, cc);
+
+  time.solve();
+  //---------------------------------------
+  // Profiling stop
+  //---------------------------------------
+  Profiling::getInstance().print();
+
+  //---------------------------------------
+  // Finalize MPI
+  //---------------------------------------
+  MPI_Finalize();
+  //---------------------------------------
+  return 0;
+}

tests/CahnHilliard/2D/test1_dt_variable/ref/time_specialized.csv

@@ -0,0 +1,24 @@
+Iter[-],Dt[s],Time[s],phi_integral[-],phi_average[-],Density[J.m-3],Sigma[J.m-3]
+1,0.01,0.01,-26.9101,-0.681641,0.189534,0.219443
+2,0.01,0.02,-26.9101,-0.681641,0.187541,0.213779
+3,0.01,0.03,-26.9101,-0.681641,0.185774,0.213026
+4,0.01,0.04,-26.9101,-0.681641,0.184989,0.212594
+5,0.01,0.05,-26.9101,-0.681641,0.184745,0.212059
+6,0.01,0.06,-26.9101,-0.681641,0.184656,0.211524
+7,0.01,0.07,-26.9101,-0.681641,0.184605,0.210962
+8,0.01,0.08,-26.9101,-0.681641,0.184498,0.210312
+9,0.01,0.09,-26.9101,-0.681641,0.184063,0.209712
+10,0.01,0.1,-26.9101,-0.681641,0.183143,0.209662
+11,0.01,0.11,-26.9101,-0.681641,0.182538,0.209495
+12,0.02,0.13,-26.9101,-0.681641,0.181943,0.209004
+13,0.02,0.15,-26.9101,-0.681641,0.181812,0.208455
+14,0.02,0.17,-26.9101,-0.681641,0.181983,0.207802
+15,0.02,0.19,-26.9101,-0.681641,0.182237,0.206944
+16,0.02,0.21,-26.9101,-0.681641,0.181214,0.206568
+17,0.04,0.25,-26.9101,-0.681641,0.179187,0.206477
+18,0.04,0.29,-26.9101,-0.681641,0.177705,0.206197
+19,0.04,0.33,-26.9101,-0.681641,0.176544,0.205984
+20,0.04,0.37,-26.9101,-0.681641,0.175738,0.205705
+21,0.04,0.41,-26.9101,-0.681641,0.175407,0.204804
+22,0.05,0.46,-26.9101,-0.681641,0.174452,0.204718
+23,0.04,0.5,-26.9101,-0.681641,0.174151,0.204395