diff --git a/.gitignore b/.gitignore index 2e9c6ec..2b7943e 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,8 @@ doc/html doc/latex +experiments/ + .aider* .vscode diff --git a/backends/dealii/include/laplace.h b/backends/dealii/include/laplace.h new file mode 100644 index 0000000..12c2296 --- /dev/null +++ b/backends/dealii/include/laplace.h @@ -0,0 +1,681 @@ +/* ------------------------------------------------------------------------ + * + * SPDX-License-Identifier: LGPL-2.1-or-later + * Copyright (C) 2010 - 2024 by the deal.II authors + * + * This file is part of the deal.II library. + * + * Part of the source code is dual licensed under Apache-2.0 WITH + * LLVM-exception OR LGPL-2.1-or-later. Detailed license information + * governing the source code and code contributions can be found in + * LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. + * + * ------------------------------------------------------------------------ + * + * Authors: Wolfgang Bangerth, Texas A&M University, 2009, 2010 + * Timo Heister, University of Goettingen, 2009, 2010 + */ + +// @sect3{Include files} +// +// Most of the include files we need for this program have already been +// discussed in previous programs. In particular, all of the following should +// already be familiar friends: +#include +#include +#include + +#include + +// This program can use either PETSc or Trilinos for its parallel +// algebra needs. By default, if deal.II has been configured with +// PETSc, it will use PETSc. Otherwise, the following few lines will +// check that deal.II has been configured with Trilinos and take that. +// +// But there may be cases where you want to use Trilinos, even though +// deal.II has *also* been configured with PETSc, for example to +// compare the performance of these two libraries. To do this, +// add the following \#define to the source code: +// @code +// #define FORCE_USE_OF_TRILINOS +// @endcode +// +// Using this logic, the following lines will then import either the +// PETSc or Trilinos wrappers into the namespace `LA` (for linear +// algebra). In the former case, we are also defining the macro +// `USE_PETSC_LA` so that we can detect if we are using PETSc (see +// solve() for an example where this is necessary). +namespace LA +{ +#if defined(DEAL_II_WITH_PETSC) && !defined(DEAL_II_PETSC_WITH_COMPLEX) && \ + !(defined(DEAL_II_WITH_TRILINOS) && defined(FORCE_USE_OF_TRILINOS)) + // # pragma message("Compiling with PETSc") + using namespace dealii::LinearAlgebraPETSc; +# define USE_PETSC_LA +#elif defined(DEAL_II_WITH_TRILINOS) + // # pragma message("Compiling with Trilinos") + using namespace dealii::LinearAlgebraTrilinos; +#else +# error DEAL_II_WITH_PETSC or DEAL_II_WITH_TRILINOS required +#endif +} // namespace LA + +#include +#include + +#include +#include + +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +// The following, however, will be new or be used in new roles. Let's walk +// through them. The first of these will provide the tools of the +// Utilities::System namespace that we will use to query things like the +// number of processors associated with the current MPI universe, or the +// number within this universe the processor this job runs on is: +#include +// The next one provides a class, ConditionOStream that allows us to write +// code that would output things to a stream (such as std::cout) +// on every processor but throws the text away on all but one of them. We +// could achieve the same by simply putting an if statement in +// front of each place where we may generate output, but this doesn't make the +// code any prettier. In addition, the condition whether this processor should +// or should not produce output to the screen is the same every time -- and +// consequently it should be simple enough to put it into the statements that +// generate output itself. +#include +// After these preliminaries, here is where it becomes more interesting. As +// mentioned in the @ref distributed topic, one of the fundamental truths of +// solving problems on large numbers of processors is that there is no way for +// any processor to store everything (e.g. information about all cells in the +// mesh, all degrees of freedom, or the values of all elements of the solution +// vector). Rather, every processor will own a few of each of these +// and, if necessary, may know about a few more, for example the ones +// that are located on cells adjacent to the ones this processor owns +// itself. We typically call the latter ghost cells, ghost nodes +// or ghost elements of a vector. The point of this discussion here is +// that we need to have a way to indicate which elements a particular +// processor owns or need to know of. This is the realm of the IndexSet class: +// if there are a total of $N$ cells, degrees of freedom, or vector elements, +// associated with (non-negative) integral indices $[0,N)$, then both the set +// of elements the current processor owns as well as the (possibly larger) set +// of indices it needs to know about are subsets of the set $[0,N)$. IndexSet +// is a class that stores subsets of this set in an efficient format: +#include +// The next header file is necessary for a single function, +// SparsityTools::distribute_sparsity_pattern. The role of this function will +// be explained below. +#include +// The final two, new header files provide the class +// parallel::distributed::Triangulation that provides meshes distributed +// across a potentially very large number of processors, while the second +// provides the namespace parallel::distributed::GridRefinement that offers +// functions that can adaptively refine such distributed meshes: +#include +#include + +#include +#include +#include +#include + +using namespace dealii; +using MPIHandle = Utilities::MPI::MPI_InitFinalize; + +// @sect3{The LaplaceProblem class template} + +// Next let's declare the main class of this program. Its structure is +// almost exactly that of the step-6 tutorial program. The only significant +// differences are: +// - The mpi_communicator variable that +// describes the set of processors we want this code to run on. In practice, +// this will be MPI_COMM_WORLD, i.e. all processors the batch scheduling +// system has assigned to this particular job. +// - The presence of the pcout variable of type ConditionOStream. +// - The obvious use of parallel::distributed::Triangulation instead of +// Triangulation. +// - The presence of two IndexSet objects that denote which sets of degrees of +// freedom (and associated elements of solution and right hand side vectors) +// we own on the current processor and which we need (as ghost elements) for +// the algorithms in this program to work. +// - The fact that all matrices and vectors are now distributed. We use +// either the PETSc or Trilinos wrapper classes so that we can use one of +// the sophisticated preconditioners offered by Hypre (with PETSc) or ML +// (with Trilinos). Note that as part of this class, we store a solution +// vector that does not only contain the degrees of freedom the current +// processor owns, but also (as ghost elements) all those vector elements +// that correspond to "locally relevant" degrees of freedom (i.e. all +// those that live on locally owned cells or the layer of ghost cells that +// surround it). +template +class LaplaceProblem +{ +public: + LaplaceProblem(bool mpi_initialize = true); + + void + run(const std::string &dir); + +private: + void + setup_system(); + void + assemble_system(); + void + solve(); + void + refine_grid(); + void + output_results(const unsigned int cycle, std::filesystem::path dir); + + std::unique_ptr mpi_handle{}; + + MPI_Comm mpi_communicator; + + parallel::distributed::Triangulation triangulation; + + const FE_Q fe; + DoFHandler dof_handler; + + IndexSet locally_owned_dofs; + IndexSet locally_relevant_dofs; + + AffineConstraints constraints; + + LA::MPI::SparseMatrix system_matrix; + LA::MPI::Vector locally_relevant_solution; + LA::MPI::Vector system_rhs; + + ConditionalOStream pcout; + TimerOutput computing_timer; +}; + +// @sect3{The LaplaceProblem class implementation} + +// @sect4{Constructor} + +inline std::unique_ptr +init_mpi(bool mpi_handle_init) +{ + if (mpi_handle_init) + { + int argc = 0; + char **argv = nullptr; + return std::move(std::make_unique(argc, argv, 1)); + } + else + { + return std::move(std::unique_ptr(nullptr)); + } +} + +// Constructors and destructors are rather trivial. In addition to what we +// do in step-6, we set the set of processors we want to work on to all +// machines available (MPI_COMM_WORLD); ask the triangulation to ensure that +// the mesh remains smooth and free to refined islands, for example; and +// initialize the pcout variable to only allow processor zero +// to output anything. The final piece is to initialize a timer that we +// use to determine how much compute time the different parts of the program +// take: +template +LaplaceProblem::LaplaceProblem(bool mpi_handle_init) + : mpi_handle(init_mpi(mpi_handle_init)) + , mpi_communicator(MPI_COMM_WORLD) + , triangulation(mpi_communicator, + typename Triangulation::MeshSmoothing( + Triangulation::smoothing_on_refinement | + Triangulation::smoothing_on_coarsening)) + , fe(2) + , dof_handler(triangulation) + , pcout(std::cout, (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)) + , computing_timer(mpi_communicator, + pcout, + TimerOutput::never, + TimerOutput::wall_times) +{} + +// @sect4{LaplaceProblem::setup_system} + +// The following function is, arguably, the most interesting one in the +// entire program since it goes to the heart of what distinguishes %parallel +// step-40 from sequential step-6. +// +// At the top we do what we always do: tell the DoFHandler object to +// distribute degrees of freedom. Since the triangulation we use here is +// distributed, the DoFHandler object is smart enough to recognize that on +// each processor it can only distribute degrees of freedom on cells it +// owns; this is followed by an exchange step in which processors tell each +// other about degrees of freedom on ghost cell. The result is a DoFHandler +// that knows about the degrees of freedom on locally owned cells and ghost +// cells (i.e. cells adjacent to locally owned cells) but nothing about +// cells that are further away, consistent with the basic philosophy of +// distributed computing that no processor can know everything. +template +void +LaplaceProblem::setup_system() +{ + TimerOutput::Scope t(computing_timer, "setup"); + + dof_handler.distribute_dofs(fe); + + pcout << " Number of active cells: " + << triangulation.n_global_active_cells() << std::endl + << " Number of degrees of freedom: " << dof_handler.n_dofs() + << std::endl; + + // The next two lines extract some information we will need later on, + // namely two index sets that provide information about which degrees of + // freedom are owned by the current processor (this information will be + // used to initialize solution and right hand side vectors, and the system + // matrix, indicating which elements to store on the current processor and + // which to expect to be stored somewhere else); and an index set that + // indicates which degrees of freedom are locally relevant (i.e. live on + // cells that the current processor owns or on the layer of ghost cells + // around the locally owned cells; we need all of these degrees of + // freedom, for example, to estimate the error on the local cells). + locally_owned_dofs = dof_handler.locally_owned_dofs(); + locally_relevant_dofs = DoFTools::extract_locally_relevant_dofs(dof_handler); + + // Next, let us initialize the solution and right hand side vectors. As + // mentioned above, the solution vector we seek does not only store + // elements we own, but also ghost entries; on the other hand, the right + // hand side vector only needs to have the entries the current processor + // owns since all we will ever do is write into it, never read from it on + // locally owned cells (of course the linear solvers will read from it, + // but they do not care about the geometric location of degrees of + // freedom). + locally_relevant_solution.reinit(locally_owned_dofs, + locally_relevant_dofs, + mpi_communicator); + system_rhs.reinit(locally_owned_dofs, mpi_communicator); + + // The next step is to compute hanging node and boundary value + // constraints, which we combine into a single object storing all + // constraints. + // + // As with all other things in %parallel, the mantra must be that no + // processor can store all information about the entire universe. As a + // consequence, we need to tell the AffineConstraints object for which + // degrees of freedom it can store constraints and for which it may not + // expect any information to store. In our case, as explained in the + // @ref distributed topic, the degrees of freedom we need to care about on + // each processor are the locally relevant ones, so we pass this to the + // AffineConstraints::reinit() function as a second argument. A further + // optimization, AffineConstraint can avoid certain operations if you also + // provide it with the set of locally owned degrees of freedom -- the + // first argument to AffineConstraints::reinit(). + // + // (What would happen if we didn't pass this information to + // AffineConstraints, for example if we called the argument-less version of + // AffineConstraints::reinit() typically used in non-parallel codes? In that + // case, the AffineConstraints class will allocate an array + // with length equal to the largest DoF index it has seen so far. For + // processors with large numbers of MPI processes, this may be very large -- + // maybe on the order of billions. The program would then allocate more + // memory than for likely all other operations combined for this single + // array. Fortunately, recent versions of deal.II would trigger an assertion + // that tells you that this is considered a bug.) + constraints.clear(); + constraints.reinit(locally_owned_dofs, locally_relevant_dofs); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + VectorTools::interpolate_boundary_values(dof_handler, + 0, + Functions::ZeroFunction(), + constraints); + constraints.close(); + + // The last part of this function deals with initializing the matrix with + // accompanying sparsity pattern. As in previous tutorial programs, we use + // the DynamicSparsityPattern as an intermediate with which we + // then initialize the system matrix. To do so, we have to tell the sparsity + // pattern its size, but as above, there is no way the resulting object will + // be able to store even a single pointer for each global degree of + // freedom; the best we can hope for is that it stores information about + // each locally relevant degree of freedom, i.e., all those that we may + // ever touch in the process of assembling the matrix (the + // @ref distributed_paper "distributed computing paper" has a long + // discussion why one really needs the locally relevant, and not the small + // set of locally active degrees of freedom in this context). + // + // So we tell the sparsity pattern its size and what DoFs to store + // anything for and then ask DoFTools::make_sparsity_pattern to fill it + // (this function ignores all cells that are not locally owned, mimicking + // what we will do below in the assembly process). After this, we call a + // function that exchanges entries in these sparsity pattern between + // processors so that in the end each processor really knows about all the + // entries that will exist in that part of the finite element matrix that + // it will own. The final step is to initialize the matrix with the + // sparsity pattern. + DynamicSparsityPattern dsp(locally_relevant_dofs); + + DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false); + SparsityTools::distribute_sparsity_pattern(dsp, + dof_handler.locally_owned_dofs(), + mpi_communicator, + locally_relevant_dofs); + + system_matrix.reinit(locally_owned_dofs, + locally_owned_dofs, + dsp, + mpi_communicator); +} + +// @sect4{LaplaceProblem::assemble_system} + +// The function that then assembles the linear system is comparatively +// boring, being almost exactly what we've seen before. The points to watch +// out for are: +// - Assembly must only loop over locally owned cells. There +// are multiple ways to test that; for example, we could compare a cell's +// subdomain_id against information from the triangulation as in +// cell->subdomain_id() == +// triangulation.locally_owned_subdomain(), or skip all cells for +// which the condition cell->is_ghost() || +// cell->is_artificial() is true. The simplest way, however, is to +// simply ask the cell whether it is owned by the local processor. +// - Copying local contributions into the global matrix must include +// distributing constraints and boundary values not just from the local +// matrix and vector into the global ones, but in the process +// also -- possibly -- from one MPI process to other processes if the +// entries we want to write to are not stored on the current process. +// Interestingly, this requires essentially no additional work: The +// AffineConstraints class we already used in step-6 is perfectly +// capable to also do this in parallel, and the only difference in this +// regard is that at the very end of the function, we have to call a +// `compress()` function on the global matrix and right hand side vector +// objects (see the description of what this does just before these calls). +// - The way we compute the right hand side (given the +// formula stated in the introduction) may not be the most elegant but will +// do for a program whose focus lies somewhere entirely different. +template +void +LaplaceProblem::assemble_system() +{ + TimerOutput::Scope t(computing_timer, "assembly"); + + const QGauss quadrature_formula(fe.degree + 1); + + FEValues fe_values(fe, + quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + const unsigned int dofs_per_cell = fe.n_dofs_per_cell(); + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); + Vector cell_rhs(dofs_per_cell); + + std::vector local_dof_indices(dofs_per_cell); + + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + { + fe_values.reinit(cell); + + cell_matrix = 0.; + cell_rhs = 0.; + + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + const double rhs_value = + (fe_values.quadrature_point(q_point)[1] > + 0.5 + + 0.25 * std::sin(4.0 * numbers::PI * + fe_values.quadrature_point(q_point)[0]) ? + 1. : + -1.); + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < dofs_per_cell; ++j) + cell_matrix(i, j) += fe_values.shape_grad(i, q_point) * + fe_values.shape_grad(j, q_point) * + fe_values.JxW(q_point); + + cell_rhs(i) += rhs_value * // + fe_values.shape_value(i, q_point) * // + fe_values.JxW(q_point); + } + } + + cell->get_dof_indices(local_dof_indices); + constraints.distribute_local_to_global( + cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs); + } + + // In the operations above, specifically the call to + // `distribute_local_to_global()` in the last line, every MPI + // process was only working on its local data. If the operation + // required adding something to a matrix or vector entry that is + // not actually stored on the current process, then the matrix or + // vector object keeps track of this for a later data exchange, + // but for efficiency reasons, this part of the operation is only + // queued up, rather than executed right away. But now that we got + // here, it is time to send these queued-up additions to those + // processes that actually own these matrix or vector entries. In + // other words, we want to "finalize" the global data + // structures. This is done by invoking the function `compress()` + // on both the matrix and vector objects. See + // @ref GlossCompress "Compressing distributed objects" + // for more information on what `compress()` actually does. + system_matrix.compress(VectorOperation::add); + system_rhs.compress(VectorOperation::add); +} + +// @sect4{LaplaceProblem::solve} + +// Even though solving linear systems on potentially tens of thousands of +// processors is by far not a trivial job, the function that does this is -- +// at least at the outside -- relatively simple. Most of the parts you've +// seen before. There are really only two things worth mentioning: +// - Solvers and preconditioners are built on the deal.II wrappers of PETSc +// and Trilinos functionality. It is relatively well known that the +// primary bottleneck of massively %parallel linear solvers is not +// actually the communication between processors, but the fact that it is +// difficult to produce preconditioners that scale well to large numbers +// of processors. Over the second half of the first decade of the 21st +// century, it has become clear that algebraic multigrid (AMG) methods +// turn out to be extremely efficient in this context, and we will use one +// of them -- either the BoomerAMG implementation of the Hypre package +// that can be interfaced to through PETSc, or a preconditioner provided +// by ML, which is part of Trilinos -- for the current program. The rest +// of the solver itself is boilerplate and has been shown before. Since +// the linear system is symmetric and positive definite, we can use the CG +// method as the outer solver. +// - Ultimately, we want a vector that stores not only the elements +// of the solution for degrees of freedom the current processor owns, but +// also all other locally relevant degrees of freedom. On the other hand, +// the solver itself needs a vector that is uniquely split between +// processors, without any overlap. We therefore create a vector at the +// beginning of this function that has these properties, use it to solve the +// linear system, and only assign it to the vector we want at the very +// end. This last step ensures that all ghost elements are also copied as +// necessary. +template +void +LaplaceProblem::solve() +{ + TimerOutput::Scope t(computing_timer, "solve"); + + LA::MPI::Vector completely_distributed_solution(locally_owned_dofs, + mpi_communicator); + + SolverControl solver_control(dof_handler.n_dofs(), + 1e-6 * system_rhs.l2_norm()); + LA::SolverCG solver(solver_control); + + LA::MPI::PreconditionAMG::AdditionalData data; +#ifdef USE_PETSC_LA + data.symmetric_operator = true; +#else + /* Trilinos defaults are good */ +#endif + LA::MPI::PreconditionAMG preconditioner; + preconditioner.initialize(system_matrix, data); + + solver.solve(system_matrix, + completely_distributed_solution, + system_rhs, + preconditioner); + + pcout << " Solved in " << solver_control.last_step() << " iterations." + << std::endl; + + constraints.distribute(completely_distributed_solution); + + locally_relevant_solution = completely_distributed_solution; +} + +// @sect4{LaplaceProblem::refine_grid} + +// The function that estimates the error and refines the grid is again +// almost exactly like the one in step-6. The only difference is that the +// function that flags cells to be refined is now in namespace +// parallel::distributed::GridRefinement -- a namespace that has functions +// that can communicate between all involved processors and determine global +// thresholds to use in deciding which cells to refine and which to coarsen. +// +// Note that we didn't have to do anything special about the +// KellyErrorEstimator class: we just give it a vector with as many elements +// as the local triangulation has cells (locally owned cells, ghost cells, +// and artificial ones), but it only fills those entries that correspond to +// cells that are locally owned. +template +void +LaplaceProblem::refine_grid() +{ + TimerOutput::Scope t(computing_timer, "refine"); + + Vector estimated_error_per_cell(triangulation.n_active_cells()); + KellyErrorEstimator::estimate( + dof_handler, + QGauss(fe.degree + 1), + std::map *>(), + locally_relevant_solution, + estimated_error_per_cell); + parallel::distributed::GridRefinement::refine_and_coarsen_fixed_number( + triangulation, estimated_error_per_cell, 0.3, 0.03); + triangulation.execute_coarsening_and_refinement(); +} + +// @sect4{LaplaceProblem::output_results} + +// Compared to the corresponding function in step-6, the one here is +// a tad more complicated. There are two reasons: the first one is +// that we do not just want to output the solution but also for each +// cell which processor owns it (i.e. which "subdomain" it is +// in). Secondly, as discussed at length in step-17 and step-18, +// generating graphical data can be a bottleneck in +// parallelizing. In those two programs, we simply generate one +// output file per process. That worked because the +// parallel::shared::Triangulation cannot be used with large numbers +// of MPI processes anyway. But this doesn't scale: Creating a +// single file per processor will overwhelm the filesystem with a +// large number of processors. +// +// We here follow a more sophisticated approach that uses +// high-performance, parallel IO routines using MPI I/O to write to +// a small, fixed number of visualization files (here 8). We also +// generate a .pvtu record referencing these .vtu files, which can +// be opened directly in visualizatin tools like ParaView and VisIt. +// +// To start, the top of the function looks like it usually does. In addition +// to attaching the solution vector (the one that has entries for all locally +// relevant, not only the locally owned, elements), we attach a data vector +// that stores, for each cell, the subdomain the cell belongs to. This is +// slightly tricky, because of course not every processor knows about every +// cell. The vector we attach therefore has an entry for every cell that the +// current processor has in its mesh (locally owned ones, ghost cells, and +// artificial cells), but the DataOut class will ignore all entries that +// correspond to cells that are not owned by the current processor. As a +// consequence, it doesn't actually matter what values we write into these +// vector entries: we simply fill the entire vector with the number of the +// current MPI process (i.e. the subdomain_id of the current process); this +// correctly sets the values we care for, i.e. the entries that correspond +// to locally owned cells, while providing the wrong value for all other +// elements -- but these are then ignored anyway. +template +void +LaplaceProblem::output_results(const unsigned int cycle, + std::filesystem::path dir) +{ + TimerOutput::Scope t(computing_timer, "output"); + + DataOut data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(locally_relevant_solution, "u"); + + Vector subdomain(triangulation.n_active_cells()); + for (unsigned int i = 0; i < subdomain.size(); ++i) + subdomain(i) = triangulation.locally_owned_subdomain(); + data_out.add_data_vector(subdomain, "subdomain"); + + data_out.build_patches(); + + // The final step is to write this data to disk. We write up to 8 VTU files + // in parallel with the help of MPI-IO. Additionally a PVTU record is + // generated, which groups the written VTU files. + data_out.write_vtu_with_pvtu_record( + dir.string() + "/", "solution", cycle, mpi_communicator, 2, 8); +} + +// @sect4{LaplaceProblem::run} + +// The function that controls the overall behavior of the program is again +// like the one in step-6. The minor difference are the use of +// pcout instead of std::cout for output to the +// console (see also step-17). +// +// A functional difference to step-6 is the use of a square domain and that +// we start with a slightly finer mesh (5 global refinement cycles) -- there +// just isn't much of a point showing a massively %parallel program starting +// on 4 cells (although admittedly the point is only slightly stronger +// starting on 1024). +template +void +LaplaceProblem::run(const std::string &dir) +{ + pcout << "Running with " +#ifdef USE_PETSC_LA + << "PETSc" +#else + << "Trilinos" +#endif + << " on " << Utilities::MPI::n_mpi_processes(mpi_communicator) + << " MPI rank(s)..." << std::endl; + + const unsigned int n_cycles = 8; + for (unsigned int cycle = 0; cycle < n_cycles; ++cycle) + { + pcout << "Cycle " << cycle << ':' << std::endl; + + if (cycle == 0) + { + GridGenerator::hyper_cube(triangulation); + triangulation.refine_global(5); + } + else + refine_grid(); + + setup_system(); + assemble_system(); + solve(); + output_results(cycle, dir); + + computing_timer.print_summary(); + computing_timer.reset(); + + pcout << std::endl; + } +} diff --git a/backends/dealii/include/register_types.h b/backends/dealii/include/register_types.h index 2ccad6f..d31371a 100644 --- a/backends/dealii/include/register_types.h +++ b/backends/dealii/include/register_types.h @@ -18,6 +18,7 @@ #include "coral.h" #include "coral_network.h" +#include "laplace.h" #include "poisson.h" /** \cond INTERNAL */ @@ -164,6 +165,13 @@ namespace coral &PoissonSolver::solve, {"PoissonSolver::solve<" + Utilities::dim_string(dim, spacedim) + ">", "poisson_solver"}); + + NodeObject::register_type, bool>("mpi_initialize"); + NodeObject::register_method, void, const std::string &>( + &LaplaceProblem::run, + {"LaplaceProblem::run<" + Utilities::dim_string(dim, spacedim) + ">", + "laplace_problem", + "output_dir"}); } diff --git a/backends/dealii/tests/CMakeLists.txt b/backends/dealii/tests/CMakeLists.txt index 2888fce..d0f000e 100644 --- a/backends/dealii/tests/CMakeLists.txt +++ b/backends/dealii/tests/CMakeLists.txt @@ -1,4 +1,5 @@ find_package(GTest) +find_package(MPI) if(NOT GTest_FOUND) message(WARNING "GTest not found; skipping deal.II backend tests") return() @@ -11,7 +12,10 @@ include(GoogleTest) # Keep plugin smoke test in its own executable so it doesn't get "help" from # backend tests that register types directly. set(_plugin_test_file ${CMAKE_CURRENT_LIST_DIR}/plugin_smoke_test.cc) -list(REMOVE_ITEM _test_files ${_plugin_test_file}) +set(_mpi_test_file ${CMAKE_CURRENT_LIST_DIR}/dealii_mpi.cc) +list(REMOVE_ITEM _test_files ${_plugin_test_file} ${_mpi_test_file}) + +# BACKEND TESTS add_executable(dealii_backend_tests ${_test_files}) @@ -47,6 +51,70 @@ gtest_discover_tests(dealii_backend_tests WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} PROPERTIES LABELS "dealii;backend") +# MPI TESTS + +add_executable(dealii_mpi_tests ${_mpi_test_file}) + +target_include_directories(dealii_mpi_tests + PRIVATE + ${CMAKE_CURRENT_LIST_DIR}/../include) + +target_link_libraries(dealii_mpi_tests + coral_core + ${GTEST_LIBRARIES} + ${GTEST_MAIN_LIBRARIES}) + +if(GTEST_INCLUDE_DIRS) + target_include_directories(dealii_mpi_tests SYSTEM PRIVATE ${GTEST_INCLUDE_DIRS}) +endif() + +deal_ii_setup_target(dealii_mpi_tests) + +if(MSVC) + target_compile_options(dealii_mpi_tests PRIVATE /W0) +else() + target_compile_options(dealii_mpi_tests PRIVATE -w) +endif() + +# Prefer running the gtest binary directly (faster iteration than ctest). +add_custom_target(run_mpi_backend_tests + COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 4 "$" + DEPENDS dealii_mpi_tests + USES_TERMINAL) + +# Register MPI tests with CTest. +# gtest_discover_tests(EXECUTOR) requires CMake >= 3.29. For 3.28 we set the +# CROSSCOMPILING_EMULATOR target property, which gtest_discover_tests respects. +# A wrapper script is used so discovery (--gtest_list_tests) runs with -n 1 +# (single clean output) while actual test execution uses -n 4. +if(MPI_FOUND) + set(_mpi_wrapper ${CMAKE_CURRENT_BINARY_DIR}/mpi_gtest_wrapper.sh) + file(WRITE ${_mpi_wrapper} +"#!/bin/bash +if [[ \"$*\" == *\"--gtest_list_tests\"* ]]; then + exec \"$@\" +else + exec ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 4 \"$@\" +fi +") + file(CHMOD ${_mpi_wrapper} + PERMISSIONS OWNER_READ OWNER_WRITE OWNER_EXECUTE + GROUP_READ GROUP_EXECUTE WORLD_READ WORLD_EXECUTE) + + set_target_properties(dealii_mpi_tests PROPERTIES + CROSSCOMPILING_EMULATOR ${_mpi_wrapper}) + + gtest_discover_tests(dealii_mpi_tests + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + PROPERTIES + LABELS "dealii;mpi" + SKIP_REGULAR_EXPRESSION "\\[ SKIPPED") +else() + message(WARNING "MPI not found; skipping MPI tests registration with ctest") +endif() + +# PLUGIN TESTS + add_executable(dealii_plugin_tests ${_plugin_test_file}) target_link_libraries(dealii_plugin_tests coral_core dl ${GTEST_LIBRARIES} ${GTEST_MAIN_LIBRARIES}) if(GTEST_INCLUDE_DIRS) diff --git a/backends/dealii/tests/dealii_mpi.cc b/backends/dealii/tests/dealii_mpi.cc new file mode 100644 index 0000000..51d8d23 --- /dev/null +++ b/backends/dealii/tests/dealii_mpi.cc @@ -0,0 +1,162 @@ +#include + +#include +#include + +#include "coral.h" +#include "coral_network.h" +#include "register_types.h" +#include "test_utils.h" + +using namespace dealii; +using namespace coral; +using coral_test::ScopedTestOutputDir; + +/** + * MPI-aware wrapper around ScopedTestOutputDir. + * + * Only the master rank (rank 0) constructs and destructs the underlying + * ScopedTestOutputDir, so all filesystem operations (create/delete) happen + * exactly once. MPI barriers ensure: + * - the directory exists before non-master ranks proceed after construction, + * - all ranks have finished writing before master rank deletes on destruction. + */ +class MPIScopedTestOutputDir +{ +public: + MPIScopedTestOutputDir( + const std::string &test_name, + const std::filesystem::path &base_dir = "./test_output") + : path_(base_dir / test_name) + { + int rank; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + if (rank == 0) + scoped_dir_ = std::make_unique(test_name, base_dir); + MPI_Barrier(MPI_COMM_WORLD); // ensure dir exists before all ranks proceed + } + + ~MPIScopedTestOutputDir() + { + MPI_Barrier(MPI_COMM_WORLD); // wait for all ranks to finish writing + // scoped_dir_ destructs only on rank 0; other ranks hold nullptr + } + + const std::filesystem::path & + path() const + { + return path_; + } + + operator const std::filesystem::path &() const + { + return path_; + } + +private: + std::filesystem::path path_; + std::unique_ptr scoped_dir_; +}; + +class DealiiMPITest : public ::testing::Test +{ +private: + using MPI_Session = Utilities::MPI::MPI_InitFinalize; + + inline static std::unique_ptr m_mpi_session = nullptr; + +protected: + void static SetUpTestSuite() + { + int argc = 0; + char **argv = nullptr; + m_mpi_session = std::make_unique(argc, argv, 1); + + if (Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD) == 1) + GTEST_SKIP() << "Please run MPI test with non trivial world size."; + } + + void static TearDownTestSuite() + { + m_mpi_session.reset(); + } +}; + +TEST_F(DealiiMPITest, Setup) +{ + int size; + ASSERT_EQ(MPI_Comm_size(MPI_COMM_WORLD, &size), MPI_SUCCESS); + ASSERT_GE(size, 0); + + int rank; + ASSERT_EQ(MPI_Comm_rank(MPI_COMM_WORLD, &rank), MPI_SUCCESS); + ASSERT_GE(rank, 0); +} + +TEST_F(DealiiMPITest, LaplaceProblem) +{ + MARK_LONG_TEST(); + + MPIScopedTestOutputDir output_dir("DealiiMPITest_PoissonSolver"); + + register_non_dimensional_types(); + register_dimensional_types<2, 2>(); + + auto laplace = make_node>(); + auto mpi_init = make_node(false); + laplace->set_arguments({mpi_init}); + ASSERT_TRUE((*laplace)()); + ASSERT_TRUE(laplace->ready()); + + auto run_method = make_node("LaplaceProblem::run<2>"); + auto output_dir_string = make_node(std::string(output_dir.path())); + run_method->set_arguments({laplace, output_dir_string}); + (*run_method)(); +} + +TEST_F(DealiiMPITest, LaplaceProblemNetwork) +{ + MARK_LONG_TEST(); + + MPIScopedTestOutputDir output_dir("dealiiMPITest_LaplaceTransformNetwork"); + + register_non_dimensional_types(); + register_dimensional_types<2, 2>(); + + Network network; + network.set_touch_file_base_path(output_dir); + network.clear_network(); + + std::ifstream file(SOURCE_DIR "/test_files/laplace.json"); + ASSERT_TRUE(file.is_open()) << "Failed to open laplace.json file."; + + nlohmann::json json_data; + file >> json_data; + file.close(); + + ASSERT_FALSE(json_data.empty()) << "JSON data is empty."; + + // Update the initialize constant + // Node 5 contains the bool value + if (json_data["workflow"]["nodes"].contains("5") && + json_data["workflow"]["nodes"]["5"].contains("value")) + { + json_data["workflow"]["nodes"]["5"]["value"] = "false"; + } + + // Update the output file path to use the test output directory + // Node 2 contains the output filename + if (json_data["workflow"]["nodes"].contains("2") && + json_data["workflow"]["nodes"]["2"].contains("value")) + { + json_data["workflow"]["nodes"]["2"]["value"] = output_dir.path().string(); + } + + network.from_json(json_data); + + slog_debug("Network has %u nodes and %u connections", + network.n_nodes(), + network.n_connections()); + + network.run(); +} diff --git a/backends/dealii/tests/test_utils.h b/backends/dealii/tests/test_utils.h index d7c7f38..aae1618 100644 --- a/backends/dealii/tests/test_utils.h +++ b/backends/dealii/tests/test_utils.h @@ -158,4 +158,17 @@ namespace coral_test } // namespace coral_test + +#define MARK_LONG_TEST() \ + do \ + { \ + if (std::getenv("CORAL_KEEP_LONG_TESTS") == nullptr) \ + { \ + std::cout << "Skip On" << std::endl; \ + GTEST_SKIP() \ + << "Long test and CORAL_KEEP_LONG_TEST disable. Skipping."; \ + } \ + } \ + while (0) + #endif // GTESTS_TEST_UTILS_H diff --git a/test_files/laplace.json b/test_files/laplace.json new file mode 100644 index 0000000..a800db7 --- /dev/null +++ b/test_files/laplace.json @@ -0,0 +1,36 @@ +{ + "workflow": { + "nodes": { + "0": { + "type": "LaplaceProblem::run<2>", + "position": { "x": 293.7562932408905, "y": 327.37380674085716 }, + "qualified_id": "0" + }, + "1": { + "type": "LaplaceProblem<2>", + "position": { "x": -183.0478681433714, "y": 226.5708997085626 }, + "qualified_id": "1" + }, + "2": { + "type": "std::string", + "position": { "x": -164.89097078018403, "y": 417.02495316403406 }, + "value": "/app/build/output/", + "qualified_id": "2" + }, + "5": { + "type": "bool", + "position": { "x": -515.9703517123914, "y": 221.60561502104412 }, + "value": "true", + "qualified_id": "5" + } + }, + "edges": { + "0": { "source": 1, "target": 0, "source_output": 0, "target_input": 0 }, + "1": { "source": 2, "target": 0, "source_output": 0, "target_input": 1 }, + "2": { "source": 5, "target": 1, "source_output": 0, "target_input": 0 } + } + }, + "version": 1, + "author": "dealiix-platform", + "date_time_utc": "2026-03-11T11:35:34.762Z" +}