diff --git a/kernel/include/Time/Time.hpp b/kernel/include/Time/Time.hpp index 318b0221..e68e0ca6 100644 --- a/kernel/include/Time/Time.hpp +++ b/kernel/include/Time/Time.hpp @@ -66,9 +66,13 @@ class TimeDiscretization { void time_management(); void time_info(const int& iter); + std::function time_step_function_; + public: // explicit TimeDiscretization(const Parameters& params, Args&&... couplings); explicit TimeDiscretization(const Parameters& params, Args... couplings); + explicit TimeDiscretization(const std::function& given_time_step, + const Parameters& params, Args... couplings); void solve(); void get_tree(); ~TimeDiscretization(); diff --git a/kernel/include/Time/Time.tpp b/kernel/include/Time/Time.tpp index d52cf76c..f56bbb2f 100644 --- a/kernel/include/Time/Time.tpp +++ b/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 TimeDiscretization::TimeDiscretization(const Parameters& params, Args... couplings) : params_(params), couplings_(std::make_tuple(std::forward(couplings)...)) { this->get_parameters(); + + MFEM_VERIFY(this->time_step_ > 0., "Error: time_step Parameter must be defined."); + + this->time_step_function_ = + std::function([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 +TimeDiscretization::TimeDiscretization( + const std::function& given_time_step, const Parameters& params, + Args... couplings) + : params_(params), couplings_(std::make_tuple(std::forward(couplings)...)) { + this->get_parameters(); + this->time_step_function_ = given_time_step; } /** @@ -61,8 +83,8 @@ void TimeDiscretization::get_parameters() { this->initial_time_ = this->params_.template get_param_value_or_default("initial_time", 0.); this->final_time_ = this->params_.template get_param_value("final_time"); - this->time_step_ = this->params_.template get_param_value("time_step"); - this->current_time_step_ = this->initial_time_; + this->time_step_ = this->params_.template get_param_value_or_default("time_step", -1.0); + this->current_time_step_ = this->time_step_; } /** @@ -194,7 +216,8 @@ void TimeDiscretization::update() { */ template void TimeDiscretization::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_; diff --git a/tests/CahnHilliard/2D/CMakeLists.txt b/tests/CahnHilliard/2D/CMakeLists.txt index 3aa83e60..ffcdbe53 100644 --- a/tests/CahnHilliard/2D/CMakeLists.txt +++ b/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) diff --git a/tests/CahnHilliard/2D/test1_dt_variable/CMakeLists.txt b/tests/CahnHilliard/2D/test1_dt_variable/CMakeLists.txt new file mode 100644 index 00000000..b46b89c6 --- /dev/null +++ b/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") diff --git a/tests/CahnHilliard/2D/test1_dt_variable/main.cpp b/tests/CahnHilliard/2D/test1_dt_variable/main.cpp new file mode 100644 index 00000000..f7c0f535 --- /dev/null +++ b/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 +#include +#include +#include +#include +#include +#include + +#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::FECollection; + using VARS = Test::VARS; + using VAR = Test::VAR; + using PST = Test::PST; + using SPA = Test::SPA; + using BCS = Test::BCS; + ///////////////////////// + using OPE = TransientOperator; + using PB = Problem; + // ########################################### + // ########################################### + // 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& 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( + [](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(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 iterations_list = {1, 3, 5}; + const std::vector times_list = {0.35, 0.45}; + + std::string calculation_path = "CahnHilliard"; + std::map> 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 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 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; +} diff --git a/tests/CahnHilliard/2D/test1_dt_variable/ref/time_specialized.csv b/tests/CahnHilliard/2D/test1_dt_variable/ref/time_specialized.csv new file mode 100644 index 00000000..8ce10ae0 --- /dev/null +++ b/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