diff --git a/kernel/include/BCs/BoundaryConditions.hpp b/kernel/include/BCs/BoundaryConditions.hpp index 64e5ec82..8a8e2827 100644 --- a/kernel/include/BCs/BoundaryConditions.hpp +++ b/kernel/include/BCs/BoundaryConditions.hpp @@ -57,7 +57,7 @@ class BoundaryConditions { public: template BoundaryConditions(SpatialDiscretization* spatial, const Args&... boundaries); - void SetBoundaryConditions(mfem::Vector& u); + void SetDirichletBoundaryConditions(mfem::Vector& u, Coefficients& coefficients, std::vector auxvars_unk); mfem::Array GetEssentialDofs(); ~BoundaryConditions(); mfem::Array get_marker_array(const std::string& boundary_type); diff --git a/kernel/include/BCs/BoundaryConditions.tpp b/kernel/include/BCs/BoundaryConditions.tpp index 687a7d30..ca195dc6 100644 --- a/kernel/include/BCs/BoundaryConditions.tpp +++ b/kernel/include/BCs/BoundaryConditions.tpp @@ -163,23 +163,86 @@ mfem::Array BoundaryConditions::get_marker_array(const std::string& } /** - * @brief Set boundary conditions + * @brief Set dirichlet boundary conditions * * @param u unknown vector + * @param coefficients coefficients list for the variable + * @param auxvars_unk unknown vectors of the auxiliary variables + * */ template -void BoundaryConditions::SetBoundaryConditions(mfem::Vector& u) { +void BoundaryConditions::SetDirichletBoundaryConditions(mfem::Vector& u, Coefficients& coefficients, std::vector auxvars_unk) { + + const int nb_bdr = this->Dirichlet_bdr_.Size(); mfem::Array tmp_array_bdr(this->Dirichlet_bdr_.Size()); - for (auto i = 0; i < this->Dirichlet_bdr_.Size(); i++) { - tmp_array_bdr = 0; - mfem::Array dof; + const int coef_size = coefficients.size(); + Coefficient dirichlet_coef = Coefficient(Glossary::Default, 0.0); + + std::vector u_values; + std::vector vaux_values; + + // Inline function to compute Dirichlet coefficient + auto compute_dirichlet_coefficient = [&](Coefficient& coef, + const std::span& values, + const std::span& aux_values) -> double { + if (coef.is_scalar()) { + return coef.compute(); + } else { + return coef.compute(values, aux_values); + } + }; + + // Loop over boundaries + for (int i = 0; i < nb_bdr; i++) { + + // If the boundary is dirichlet if (this->Dirichlet_bdr_[i] > 0) { + + // Check if there is a coefficient given and store it + bool has_dirichlet_coef = false; + for (int l = 0; l < coef_size; l++) { + const auto& coef = coefficients[l]; + if (coef.get_type() == GlossaryType::Dirichlet) { + auto bdr_ids = coef.get_bdr_index_coef(); + if (std::find(bdr_ids.begin(), bdr_ids.end(), i) != bdr_ids.end()) { + dirichlet_coef = coef; + has_dirichlet_coef = true; + break; + } + } + } + + // Get the list of essential true dofs + tmp_array_bdr = 0; tmp_array_bdr[i] = 1; - this->fespace_->GetEssentialTrueDofs(tmp_array_bdr, dof); - u.SetSubVector(dof, this->Dirichlet_value_[i]); + mfem::Array dof_list; + this->fespace_->GetEssentialTrueDofs(tmp_array_bdr, dof_list); + + if (has_dirichlet_coef) { + mfem::Vector dirichlet_at_dofs(dof_list.Size()); + // Loop over essential dofs and compute the coefficient at the dofs + for (int j = 0; j < dof_list.Size(); j++) { + u_values.clear(); + vaux_values.clear(); + + int dof = dof_list[j]; + u_values.push_back(u(dof)); + for (auto aux : auxvars_unk) + vaux_values.emplace_back(aux(dof)); + + dirichlet_at_dofs[j] = compute_dirichlet_coefficient(dirichlet_coef, std::span(u_values), + std::span(vaux_values)); + } + + // SetSubVector with the calculated dirichlet values + u.SetSubVector(dof_list, dirichlet_at_dofs); + } else { + // If no dirichlet coefficient is given, use the constant dirichlet value + u.SetSubVector(dof_list, this->Dirichlet_value_[i]); + } } } -} // end of SetBoundaryConditions +} // end of SetDirichletBoundaryConditions /** * @brief Destroy the Boundary Conditions:: Boundary Conditions object diff --git a/kernel/include/Glossary/Glossary.hpp b/kernel/include/Glossary/Glossary.hpp index 6831fd6f..ca948859 100644 --- a/kernel/include/Glossary/Glossary.hpp +++ b/kernel/include/Glossary/Glossary.hpp @@ -120,7 +120,8 @@ enum class GlossaryType { Neumann, RobinA, RobinB, - ExplicitTime + ExplicitTime, + Dirichlet }; struct GlossaryQuantity { @@ -359,6 +360,13 @@ static const GlossaryQuantity Robin_a = static const GlossaryQuantity Robin_b = GlossaryQuantity(GlossaryType::RobinB, GlossaryUnit::None, "B-Robin boundary condition"); +/** + * @brief Quantity associated with the Dirichlet boundary condition + * + */ +static const GlossaryQuantity Dirichlet = + GlossaryQuantity(GlossaryType::Dirichlet, GlossaryUnit::None, "Dirichlet boundary condition"); + /** * @brief Quantity associated with the MPI rank * diff --git a/kernel/include/Operators/OperatorBase.tpp b/kernel/include/Operators/OperatorBase.tpp index ee7a391d..5665ba25 100644 --- a/kernel/include/Operators/OperatorBase.tpp +++ b/kernel/include/Operators/OperatorBase.tpp @@ -374,7 +374,15 @@ void OperatorBase::initialize([[maybe_unused]] const double& initial_tim this->bcs_.emplace_back(vv.get_boundary_conditions()); this->ess_tdof_list_.emplace_back(this->bcs_[iv]->GetEssentialDofs()); - this->bcs_[iv]->SetBoundaryConditions(u); + + // Get the unknown vectors of the auxiliary variables + std::vector auxvars_unk; + for (const auto& auxvar_vec : this->auxvariables_) { + for (const auto& auxvar : auxvar_vec->getVariables()) { + auxvars_unk.emplace_back(auxvar.get_unknown()); + } + } + this->bcs_[iv]->SetDirichletBoundaryConditions(u, this->coefficients_[iv], auxvars_unk); vv.update(u); u_vect.emplace_back(u); } diff --git a/kernel/include/Operators/SteadyOperator.tpp b/kernel/include/Operators/SteadyOperator.tpp index 990ed9b0..dbd979d9 100644 --- a/kernel/include/Operators/SteadyOperator.tpp +++ b/kernel/include/Operators/SteadyOperator.tpp @@ -168,9 +168,16 @@ void SteadyOperator::solve(std::vector>& v this->SetTransientParameters(dt, u_vect); /// Apply BCs: check if need to be uncomment + // Get the unknown vectors of the auxiliary variables + // std::vector auxvars_unk; + // for (const auto& auxvar_vec : this->auxvariables_) { + // for (const auto& auxvar : auxvar_vec->getVariables()) { + // auxvars_unk.emplace_back(auxvar.get_unknown()); + // } + // } // for (size_t i = 0; i < unk_size; i++) { // auto &unk_i = *(vect_unk[i]); - // this->bcs_[i]->SetBoundaryConditions(unk_i); + // this->bcs_[i]->SetDirichletBoundaryConditions(unk_i, this->coefficients_[i], auxvars_unk); // } // Source term diff --git a/kernel/include/Operators/TransientOperator.tpp b/kernel/include/Operators/TransientOperator.tpp index d1754389..3c47b4e8 100644 --- a/kernel/include/Operators/TransientOperator.tpp +++ b/kernel/include/Operators/TransientOperator.tpp @@ -534,11 +534,19 @@ void TransientOperator::ImplicitSolve(const double dt, const mfem::Vecto { MATools::MATrace::start(); Catch_Time_Section("ImplicitSolve::ApplyBCs"); + + // Get the unknown vectors of the auxiliary variables + std::vector auxvars_unk; + for (const auto& auxvar_vec : this->auxvariables_) { + for (const auto& auxvar : auxvar_vec->getVariables()) { + auxvars_unk.emplace_back(auxvar.get_unknown()); + } + } auto sc_1 = 0; auto sc_2 = sc / fes_size; for (int i = 0; i < fes_size; ++i) { mfem::Vector v_i(u.GetData() + sc_1, sc_2); - this->bcs_[i]->SetBoundaryConditions(v_i); + this->bcs_[i]->SetDirichletBoundaryConditions(v_i, this->coefficients_[i], auxvars_unk); sc_1 += sc_2; } reduced_oper->SetParameters(dt, &v); diff --git a/kernel/include/Problems/Problem.tpp b/kernel/include/Problems/Problem.tpp index 7a12667d..0a77bdab 100644 --- a/kernel/include/Problems/Problem.tpp +++ b/kernel/include/Problems/Problem.tpp @@ -394,7 +394,9 @@ void Problem::do_time_step( [[maybe_unused]] const std::vector>& unks_info) { const size_t unk_size = vect_unk.size(); - this->set_time_coefficients(next_time); + // Set time for coefficients, /!\ this is correct only for implicit time scheme + // TODO: use effective time step according to the time scheme + this->set_time_coefficients(next_time+current_time_step); this->oper_.setGeometry(this->geometry_); this->oper_.solve(vect_unk, next_time, current_time, current_time_step, iter); diff --git a/tests/HeatTransfer/2D/CMakeLists.txt b/tests/HeatTransfer/2D/CMakeLists.txt index 03456f5e..7db3404b 100644 --- a/tests/HeatTransfer/2D/CMakeLists.txt +++ b/tests/HeatTransfer/2D/CMakeLists.txt @@ -2,3 +2,4 @@ add_subdirectory(test1) add_subdirectory(test2) add_subdirectory(test3) add_subdirectory(test4) +add_subdirectory(test5) diff --git a/tests/HeatTransfer/2D/test5/CMakeLists.txt b/tests/HeatTransfer/2D/test5/CMakeLists.txt new file mode 100644 index 00000000..b4962b92 --- /dev/null +++ b/tests/HeatTransfer/2D/test5/CMakeLists.txt @@ -0,0 +1,6 @@ +execute_process(COMMAND python3 ${BUILD_SCRIPT_DIR}/GenerateCoefficient.py -r -f ${CMAKE_CURRENT_SOURCE_DIR}/coefficient.json + WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}) + +create_test("HeatTransfer2Dtest5" "HeatTransfer2Dtest5" FALSE "2D EXE HEAT" 1) +configure_file(${CMAKE_SOURCE_DIR}/tests/tools/convergence_study.py ${CMAKE_CURRENT_BINARY_DIR}/convergence_study.py COPYONLY) +create_col_comparison_convergence("CompareHeat2Dtest5Convergence" "convergence_output_ref.csv" "convergence_output.csv" -1 absolute 1e-16 FALSE "HeatTransfer2Dtest5" "2D Heat" 1. 100000) \ No newline at end of file diff --git a/tests/HeatTransfer/2D/test5/Coefficient.hpp b/tests/HeatTransfer/2D/test5/Coefficient.hpp new file mode 100644 index 00000000..8e4a9671 --- /dev/null +++ b/tests/HeatTransfer/2D/test5/Coefficient.hpp @@ -0,0 +1,190 @@ +/** + * + * 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 . + * + */ + +#include +#include +#include +#include +#include +#include + + +#include "Options/PhysicalPropertiesOptions.hpp" + +#include "Coefficients/FunctionCoefficient.hpp" + + +#pragma once + +/** + * + * @brief C++ function of the analytical expression + * + * F = -pi*t*cos(pi*x) + */ +class NeumannCoefficient : public FunctionCoefficient { + private: + double prefactor_; + protected: + std::function&,const std::span&, const std::span&, const unsigned int dimension)> F() final; + std::function(const std::span&,const std::span&, const std::span&, const unsigned int dimension)> GradientF() final; + std::function(const std::span&,const std::span&,const std::span&, const unsigned int dimension)> HessianF() final; + + public: + NeumannCoefficient() : prefactor_(1.0) {} + explicit NeumannCoefficient(const double prefactor): prefactor_(prefactor) {} + virtual ~NeumannCoefficient() = default; +}; + +/** + * + * @brief C++ function of the expression + * + * + * @return std::function&,const std::span&, const std::span&, const unsigned int dimension)> + */ + std::function&,const std::span&, const std::span&, const unsigned int dimension)> NeumannCoefficient::F() { + auto func = [&](const std::span& input_vector, [[maybe_unused]] const std::span&,const std::span& auxiliary_vector, [[maybe_unused]] const unsigned int dimension) { + double T = input_vector[0]; + double x = auxiliary_vector[0]; + double y = auxiliary_vector[1]; + double t = this->time_; + double F = -M_PI*t*std::cos(M_PI*x); + return this->prefactor_ * F; + }; + return func; +} + +/** + * + * @brief Gradient + * + * @return std::function(const std::span&,const std::span&, const unsigned int dimension)> + */ +std::function(const std::span&,const std::span&, const std::span&, const unsigned int dimension)> NeumannCoefficient::GradientF() { + auto func = [&](const std::span& input_vector, [[maybe_unused]] const std::span&,const std::span& auxiliary_vector, [[maybe_unused]] const unsigned int dimension) { + double T = input_vector[0]; + double x = auxiliary_vector[0]; + double y = auxiliary_vector[1]; + double t = this->time_; + std::vector gradient(1); + gradient[0] = this->prefactor_ * (0); + return gradient; + }; + return func; +} + +/** + * + * @brief Hessian + * @remark Hessian matrix stored in vector : H(i,j)->H(i*n+j) + * + * @return std::function(const std::span&,const std::span&, const std::span&, const unsigned int dimension)> + */ +std::function(const std::span&,const std::span&, const std::span&, const unsigned int dimension)> NeumannCoefficient::HessianF() { + auto func = [&](const std::span& input_vector, [[maybe_unused]] const std::span&,const std::span& auxiliary_vector, [[maybe_unused]] const unsigned int dimension) { + double T = input_vector[0]; + double x = auxiliary_vector[0]; + double y = auxiliary_vector[1]; + double t = this->time_; + std::vector hessian(1); + hessian[0] = this->prefactor_ * (0); + return hessian; + }; + return func; +} +/** + * + * @brief C++ function of the analytical expression + * + * F = t*sin(pi*y) + */ +class DirichletCoefficient : public FunctionCoefficient { + private: + double prefactor_; + protected: + std::function&,const std::span&, const std::span&, const unsigned int dimension)> F() final; + std::function(const std::span&,const std::span&, const std::span&, const unsigned int dimension)> GradientF() final; + std::function(const std::span&,const std::span&,const std::span&, const unsigned int dimension)> HessianF() final; + + public: + DirichletCoefficient() : prefactor_(1.0) {} + explicit DirichletCoefficient(const double prefactor): prefactor_(prefactor) {} + virtual ~DirichletCoefficient() = default; +}; + +/** + * + * @brief C++ function of the expression + * + * + * @return std::function&,const std::span&, const std::span&, const unsigned int dimension)> + */ + std::function&,const std::span&, const std::span&, const unsigned int dimension)> DirichletCoefficient::F() { + auto func = [&](const std::span& input_vector, [[maybe_unused]] const std::span&,const std::span& auxiliary_vector, [[maybe_unused]] const unsigned int dimension) { + double T = input_vector[0]; + double x = auxiliary_vector[0]; + double y = auxiliary_vector[1]; + double t = this->time_; + double F = t*std::sin(M_PI*y); + return this->prefactor_ * F; + }; + return func; +} + +/** + * + * @brief Gradient + * + * @return std::function(const std::span&,const std::span&, const unsigned int dimension)> + */ +std::function(const std::span&,const std::span&, const std::span&, const unsigned int dimension)> DirichletCoefficient::GradientF() { + auto func = [&](const std::span& input_vector, [[maybe_unused]] const std::span&,const std::span& auxiliary_vector, [[maybe_unused]] const unsigned int dimension) { + double T = input_vector[0]; + double x = auxiliary_vector[0]; + double y = auxiliary_vector[1]; + double t = this->time_; + std::vector gradient(1); + gradient[0] = this->prefactor_ * (0); + return gradient; + }; + return func; +} + +/** + * + * @brief Hessian + * @remark Hessian matrix stored in vector : H(i,j)->H(i*n+j) + * + * @return std::function(const std::span&,const std::span&, const std::span&, const unsigned int dimension)> + */ +std::function(const std::span&,const std::span&, const std::span&, const unsigned int dimension)> DirichletCoefficient::HessianF() { + auto func = [&](const std::span& input_vector, [[maybe_unused]] const std::span&,const std::span& auxiliary_vector, [[maybe_unused]] const unsigned int dimension) { + double T = input_vector[0]; + double x = auxiliary_vector[0]; + double y = auxiliary_vector[1]; + double t = this->time_; + std::vector hessian(1); + hessian[0] = this->prefactor_ * (0); + return hessian; + }; + return func; +} diff --git a/tests/HeatTransfer/2D/test5/coefficient.json b/tests/HeatTransfer/2D/test5/coefficient.json new file mode 100644 index 00000000..7327ccf0 --- /dev/null +++ b/tests/HeatTransfer/2D/test5/coefficient.json @@ -0,0 +1,18 @@ +[ + { + "expression":"-pi*t*cos(pi*x)", + "variables":"T", + "auxiliary_variables":"x,y", + "constants":"(t:T)", + "class_name":"NeumannCoefficient", + "outputfile":"Coefficient" + }, + { + "expression":"t*sin(pi*y)", + "variables":"T", + "auxiliary_variables":"x,y", + "constants":"(t:T)", + "class_name":"DirichletCoefficient", + "outputfile":"Coefficient" + } +] diff --git a/tests/HeatTransfer/2D/test5/main.cpp b/tests/HeatTransfer/2D/test5/main.cpp new file mode 100644 index 00000000..2648af43 --- /dev/null +++ b/tests/HeatTransfer/2D/test5/main.cpp @@ -0,0 +1,207 @@ +/** + * @file main.cpp + * @author mh286406 (marine.harel@cea.fr) + * @brief 2D heat transfer problem + * @version 0.1 + * @date 2026-07-17 + * + * @copyright Copyright (c) 2026 + * + */ +#include +#include +#include +#include +#include +#include +#include + +#include "./Coefficient.hpp" +#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 + 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 int refinement_level = 0; + + std::vector vect_elem{"InlineSquareWithQuadrangles"}; + std::vector vect_order{2, 1}; + std::vector vect_NN{160, 80, 40, 20}; + for (const auto& elem_type : vect_elem) { + for (const auto order_fe : vect_order) { + for (const auto NN : vect_NN) { + const int nx = NN; + const int ny = NN; + const double lx = 1.; + const double ly = 1.; + + const std::tuple& tuple_of_dimensions = std::make_tuple( + nx, ny, lx, ly); // Number of elements and maximum length in each direction + + SPA spatial(elem_type, order_fe, refinement_level, tuple_of_dimensions); + + // ############################## + // Boundary conditions // + // ############################## + + auto boundaries = {Boundary("lower", 0, "Neumann"), Boundary("right", 1, "Dirichlet"), + Boundary("upper", 2, "Neumann"), Boundary("left", 3, "Dirichlet")}; + auto bcs = BCS(&spatial, boundaries); + + auto Xboundaries = {Boundary("lower", 0, "Neumann"), Boundary("right", 1, "Neumann"), + Boundary("upper", 2, "Neumann"), Boundary("left", 3, "Neumann")}; + auto Xbcs = BCS(&spatial, Xboundaries); + + // ########################################### + // ########################################### + // Physical models // + // ########################################### + // ########################################### + // #################### + // parameters // + // #################### + // Heat + + Coefficient density(Glossary::Concentration, 1.); + Coefficient heat_capacity(Glossary::Cp, 1.); + Coefficient conductivity(Glossary::Conductivity, 1.); + Coefficient neumann(Glossary::Neumann, Scheme::Implicit, NeumannCoefficient()); + Coefficient dirichlet_left(Glossary::Dirichlet, Scheme::Implicit, DirichletCoefficient()); + Coefficient dirichlet_right(Glossary::Dirichlet, Scheme::Implicit, DirichletCoefficient(-1)); + neumann.set_bdr_index_coef(std::vector{0,2}); + dirichlet_left.set_bdr_index_coef(std::vector{3}); + dirichlet_right.set_bdr_index_coef(std::vector{1}); + + // #################### + // variables // + // #################### + + auto user_func = std::function( + [](const mfem::Vector& x, double time) { + const auto xx = x[0]; + const auto yy = x[1]; + const auto func = time * std::cos(M_PI*xx) * std::sin(M_PI*yy); + return func; + }); + auto T_analytical = AnalyticalFunctions(user_func); + + auto heat_vars = VARS(VAR(&spatial, bcs, "T", Glossary::Temperature, 2, 0, T_analytical)); + + // Coord + auto xcoord = std::function( + [](const mfem::Vector& vcoord, double time) { return vcoord[0]; }); + auto ycoord = std::function( + [](const mfem::Vector& vcoord, double time) { return vcoord[1]; }); + auto XC = VAR(&spatial, Xbcs, "XCOORD", Glossary::Coordinate, 2, + AnalyticalFunctions(xcoord)); + XC.set_additional_information("XCOORD"); + auto YC = VAR(&spatial, Xbcs, "YCOORD", Glossary::Coordinate, 2, + AnalyticalFunctions(ycoord)); + YC.set_additional_information("YCOORD"); + auto coord = VARS(XC, YC); + + // ########################################### + // ########################################### + // Post-processing // + // ########################################### + // ########################################### + const std::string& main_folder_path = + "Saves_order_" + std::to_string(order_fe) + "_Nx" + std::to_string(NN); + const auto& level_of_detail = 1; + const auto& frequency = 1; + // Heat + const std::string& calculation_path = "Problem1"; + auto p_pst = Parameters( + Parameter("main_folder_path", main_folder_path), + Parameter("calculation_path", calculation_path), Parameter("frequency", frequency), + Parameter("level_of_detail", level_of_detail), Parameter("enable_compute_energies", false)); + auto pst = PST(&spatial, p_pst); + + // #################### + // Problems // + // #################### + + auto user_func_source_term = std::function( + [](const mfem::Vector& x, [[maybe_unused]] double time) { + const auto xx = x[0]; + const auto yy = x[1]; + const auto func = (1+2*M_PI*M_PI*time) * std::cos(M_PI*xx) * std::sin(M_PI*yy); + return func; + }); + + std::vector > src_term; + src_term.emplace_back(AnalyticalFunctions(user_func_source_term)); + + // Heat + Coefficients coef_pb(density, heat_capacity, conductivity, neumann, dirichlet_left, dirichlet_right); + std::vector spatials{&spatial}; + OPE oper(spatials, {"Fourier"}, TimeScheme::EulerImplicit, "HeatTimeDerivative", src_term); + + 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-12))); + PB Heat_pb("Heat", oper, heat_vars, {coef_pb}, pst, coord); + + // Coupling 1 + auto cc = Coupling("Heat transfer", Heat_pb); + // ########################################### + // ########################################### + // Time-integration // + // ########################################### + // ########################################### + const auto& t_initial = 0.0; + const auto& t_final = 1.; + const auto& dt = 1.; + auto time_params = Parameters(Parameter("initial_time", t_initial), + Parameter("final_time", t_final), Parameter("time_step", dt)); + auto time = TimeDiscretization(time_params, cc); + + time.solve(); + //--------------------------------------- + // Profiling stop + //--------------------------------------- + Profiling::getInstance().print(); + } + } + } + //--------------------------------------- + // Finalize MPI + //--------------------------------------- + MPI_Finalize(); + //--------------------------------------- + return 0; +} diff --git a/tests/HeatTransfer/2D/test5/ref/convergence_output_ref.csv b/tests/HeatTransfer/2D/test5/ref/convergence_output_ref.csv new file mode 100644 index 00000000..a5a38a36 --- /dev/null +++ b/tests/HeatTransfer/2D/test5/ref/convergence_output_ref.csv @@ -0,0 +1,8 @@ +order 2,0.00625,4.17294e-08 +order 2,0.0125,2.9319e-07 +order 2,0.025,2.1636e-06 +order 2,0.05,1.65153e-05 +order 1,0.00625,2.73767e-05 +order 1,0.0125,0.000109681 +order 1,0.025,0.00044008 +order 1,0.05,0.00177063