diff --git a/.github/workflows/pages.yml b/.github/workflows/pages.yml index bede61d6d..7d11084b1 100644 --- a/.github/workflows/pages.yml +++ b/.github/workflows/pages.yml @@ -161,6 +161,9 @@ jobs: mkdir docs_out chmod a+rwx docs_out docker run -v ${PWD}:/src ghcr.io/cexa-project/ddc/doxygen:${GITHUB_SHA:0:7} bash /src/run.sh + if [ -f docs_out/doxygen.log ]; then + cat docs_out/doxygen.log + fi - name: Publish site if: ${{ github.event_name == 'push' && github.ref_name == 'main' && needs.id_repo.outputs.in_base_repo == 'true' }} run: | diff --git a/docs/CMakeLists.txt b/docs/CMakeLists.txt index 76974a080..411e35490 100644 --- a/docs/CMakeLists.txt +++ b/docs/CMakeLists.txt @@ -59,7 +59,8 @@ set(DOXYGEN_IMAGE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/images") doxygen_add_docs(doc "${CMAKE_CURRENT_SOURCE_DIR}/About.md" "${CMAKE_CURRENT_SOURCE_DIR}/first_steps.md" + "${CMAKE_CURRENT_SOURCE_DIR}/going_further.md" "${CMAKE_CURRENT_SOURCE_DIR}/uniform_heat_equation.md" - "${CMAKE_CURRENT_SOURCE_DIR}/heat_equation.md" + "${CMAKE_CURRENT_SOURCE_DIR}/non_uniform_heat_equation.md" "${DDC_SOURCE_DIR}/include/ddc/" ALL) diff --git a/docs/first_steps.md b/docs/first_steps.md index 1b30af633..3d10b84b6 100644 --- a/docs/first_steps.md +++ b/docs/first_steps.md @@ -5,7 +5,7 @@ Copyright (C) The DDC development team, see COPYRIGHT.md file SPDX-License-Identifier: MIT --> -In \subpage heat_equation "examples/uniform_heat_equation.cpp" is a DDC example implementing a forward +In \subpage uniform_heat_equation "examples/uniform_heat_equation.cpp" is a DDC example implementing a forward finite-difference solver for the heat equation over a rectangle 2D domain with periodic boundary conditions. @@ -175,6 +175,33 @@ And we display the initial data. \snippet uniform_heat_equation.cpp initial-display + +\snippet uniform_heat_equation.cpp time iteration + +To display the data, a chunk is created on the host. + + +\snippet uniform_heat_equation.cpp host-chunk + +We deepcopy the data from the `ghosted_last_temp` chunk to `ghosted_temp` on the host. + +\snippet uniform_heat_equation.cpp boundary conditions + +\snippet uniform_heat_equation.cpp initial-deepcopy + +And we display the initial data. + +\snippet uniform_heat_equation.cpp initial-display + +For the numerical scheme, two chunkspans are created: ++ `next_temp` a span excluding ghosts of the temperature at the time-step we will build. ++ `last_temp` a read-only view of the temperature at the previous time-step.Note that *span_cview* returns a read-only ChunkSpan. + +\snippet uniform_heat_equation.cpp manipulated views + +We then solve the equation. + +\snippet uniform_heat_equation.cpp numerical scheme # Time loop \snippet uniform_heat_equation.cpp time iteration @@ -182,11 +209,16 @@ And we display the initial data. ## Periodic conditions +\snippet uniform_heat_equation.cpp output + \snippet uniform_heat_equation.cpp boundary conditions ## Numerical scheme + +\snippet uniform_heat_equation.cpp swap + For the numerical scheme, two chunkspans are created: + `next_temp` a span excluding ghosts of the temperature at the time-step we will build. + `last_temp` a read-only view of the temperature at the previous time-step.Note that *span_cview* returns a read-only ChunkSpan. @@ -195,4 +227,7 @@ For the numerical scheme, two chunkspans are created: We then solve the equation. + +\snippet uniform_heat_equation.cpp final output + \snippet uniform_heat_equation.cpp numerical scheme diff --git a/docs/going_further.md b/docs/going_further.md new file mode 100644 index 000000000..3feec26cf --- /dev/null +++ b/docs/going_further.md @@ -0,0 +1,123 @@ +# Commented example: the non uniform heat equation {#going_further} + + +In \subpage non_uniform_heat_equation "examples/non_uniform_heat_equation.cpp" is a DDC example implementing a forward +finite-difference solver for the heat equation over a rectangle 2D domain with periodic boundary +conditions and non-uniform space discretization. + +As usual, the file starts with a few includes that will be used in the code. + +\snippet non_uniform_heat_equation.cpp includes + +# Differences with the uniform problem resolution + +For non-uniform discretization, differences arise compared to uniform discretization, particularly in terms of dimension labeling, domain construction, and even in the resolution of the problem. + +## Dimensions naming + +Just like the uniform case, we start by defining types that we later use to name our +dimensions. + +\snippet non_uniform_heat_equation.cpp X-dimension + +However, we then define the discretization as non uniform. + +\snippet non_uniform_heat_equation.cpp X-discretization + +We do the same thing for the second dimension `Y`. + +\snippet non_uniform_heat_equation.cpp Y-space + +And for this case, we kept an uniform discretization for the time dimension. + +\snippet non_uniform_heat_equation.cpp time-space + +## Domains + +### Dimension X + +Just like the uniform case, we define the main characteristics of the domain. + +\snippet non_uniform_heat_equation.cpp main-start-x-parameters + +The first step to create a `DiscreteDomain` with non-uniform spatial discretization is to create a C++ iterator containing the actual coordinates of the points along each dimension. For tutorial purposes, we have constructed a function `generate_random_vector` that generates a vector with random data ranging from the lower bound to the higher one to populate our vectors. + +For the `X` dimension, we want to build 4 domains: +* `x_domain`: the main domain from the start (included) to end (included) but excluding "ghost" + points. +* `ghosted_x_domain`: the domain including all "ghost" points. +* `x_pre_ghost`: the "ghost" points that come before the main domain. +* `x_post_ghost`: the "ghost" points that come after the main domain. + +To do so, we have to create the iterator for the main domain. + +\snippet non_uniform_heat_equation.cpp iterator_main-domain + +For the initialization of ghost points in the non-uniform case, it was necessary to explicitly describe the position of the pre-ghost and post-ghost points. For the pre-ghost, it was necessary to start from the first coordinate along x and subtract the difference between the last and the penultimate coordinate of the x dimension to maintain periodicity. For the post-ghost point, take the last point and add the difference between the first and second points along the x dimension. + +\snippet non_uniform_heat_equation.cpp ghost_points_x + +Then we can create our 4 domains using the `init_discretization` +function with the `init_ghosted` function that takes the vectors as parameters. + +\snippet non_uniform_heat_equation.cpp build-domains + +**To summarize,** in this section we saw how to build a domain with non-uniform space discretization: ++ Create an iterator for your domain. ++ Call the `init_discretization` function with `init_ghosted` if you need ghost points or with `init` for a regular domain. + +### Dimension Y + +For the `Y` dimension, we do the same. We start by defining our 3 vectors. + +\snippet non_uniform_heat_equation.cpp Y-vectors + +And we use them to build the 4 domains in the same way as for the X dimension. + +\snippet non_uniform_heat_equation.cpp build-Y-domain + +### Time dimension + +Then we handle the domains for the simulated time dimension. We first give the simulated time at which to stard and end the simulation. + +\snippet non_uniform_heat_equation.cpp main-start-t-parameters + +The CFL conditions are more challenging to achieve for the case of non-uniform discretization. +Since the spatial steps are not uniform, we first need to find the maximum of the inverse of the square of the spatial step for the `X` and `Y` dimensions. And then we obtain the minimum value for `dt`. + +\snippet non_uniform_heat_equation.cpp CFL-condition + +We can calculate the number of time steps and build the `DiscreteDomain` for the time dimension. + +\snippet non_uniform_heat_equation.cpp time-domain + +## Time loop + +Allocation and initialization are the same as for the uniform case. Let's focus on resolving the numerical scheme. +The main difference in solving the numerical equation is that we need to account for the fact that the values of dx and dy on the left and right sides are different. We use the functions `distance_at_left` and `distance_at_right` to solve the equation. + +\snippet non_uniform_heat_equation.cpp numerical scheme + + + + + + + + + + + + + + + + + + + + diff --git a/docs/heat_equation.md b/docs/heat_equation.md deleted file mode 100644 index 555ac22d5..000000000 --- a/docs/heat_equation.md +++ /dev/null @@ -1,8 +0,0 @@ -# examples/heat_equation.cpp {#heat_equation} - - -\include{lineno} heat_equation.cpp diff --git a/docs/non_uniform_heat_equation.md b/docs/non_uniform_heat_equation.md new file mode 100644 index 000000000..9e2bf9a35 --- /dev/null +++ b/docs/non_uniform_heat_equation.md @@ -0,0 +1,8 @@ +# examples/non_uniform_heat_equation.cpp {#non_uniform_heat_equation} + + +\include{lineno} non_uniform_heat_equation.cpp \ No newline at end of file diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index c0392b47d..c7b19a0d6 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -6,11 +6,18 @@ add_executable(heat_equation heat_equation.cpp) target_link_libraries(heat_equation PUBLIC DDC::DDC) add_executable(uniform_heat_equation uniform_heat_equation.cpp) +add_executable(non_uniform_heat_equation non_uniform_heat_equation.cpp) + +target_link_libraries(uniform_heat_equation PUBLIC DDC::DDC) +target_link_libraries(non_uniform_heat_equation PUBLIC DDC::DDC) target_link_libraries(uniform_heat_equation PUBLIC DDC::DDC) + if("${DDC_BUILD_PDI_WRAPPER}") target_link_libraries(heat_equation PUBLIC DDC::PDI_Wrapper) target_link_libraries(uniform_heat_equation PUBLIC DDC::PDI_Wrapper) + target_link_libraries(non_uniform_heat_equation PUBLIC DDC::PDI_Wrapper) + endif() if("${DDC_BUILD_KERNELS_FFT}") add_executable(heat_equation_spectral heat_equation_spectral.cpp) diff --git a/examples/non_uniform_heat_equation.cpp b/examples/non_uniform_heat_equation.cpp new file mode 100644 index 000000000..1fcdf5218 --- /dev/null +++ b/examples/non_uniform_heat_equation.cpp @@ -0,0 +1,401 @@ +// Copyright (C) The DDC development team, see COPYRIGHT.md file +// +// SPDX-License-Identifier: MIT + +//! [includes] +#include +#include +#include +#include +#include +#include +#include + +#include +//! [includes] + +//! [vector_generator] +std::vector generate_random_vector( + int n, + int lower_bound, + int higher_bound) +{ + std::random_device rd; + std::mt19937 gen(rd()); + std::uniform_real_distribution + dis(lower_bound, higher_bound); + + std::vector vec(n); + vec[0] = lower_bound; + vec[n - 1] = higher_bound; + + for (int i = 1; i < vec.size() - 1; ++i) { + vec[i] = dis(gen); + } + + std::sort(vec.begin(), vec.end()); + return vec; +} +//! [vector_generator] + +//! [X-dimension] +struct X; +//! [X-dimension] + +//! [X-discretization] +struct DDimX : ddc::NonUniformPointSampling +{ +}; +//! [X-discretization] + +//! [Y-space] +struct Y; +struct DDimY : ddc::NonUniformPointSampling +{ +}; +//! [Y-space] + +//! [time-space] +struct T; +struct DDimT : ddc::UniformPointSampling +{ +}; +//! [time-space] + +//! [display] +/** A function to pretty print the temperature + * @tparam ChunkType The type of chunk span. This way the template parameters are avoided, + * should be deduced by the compiler. + * @param time The time at which the output is made. + * @param temp The temperature at this time-step. + */ +template +void display(double time, ChunkType temp) +{ + double const mean_temp = ddc::transform_reduce( + temp.domain(), + 0., + ddc::reducer::sum(), + temp) + / temp.domain().size(); + std::cout << std::fixed << std::setprecision(3); + std::cout << "At t = " << time << ",\n"; + std::cout << " * mean temperature = " << mean_temp << "\n"; + ddc::ChunkSpan temp_slice + = temp[ddc::get_domain(temp).front() + + ddc::get_domain(temp).size() / 2]; + std::cout << " * temperature[y:" + << ddc::get_domain(temp).size() / 2 << "] = {"; + ddc::for_each( + ddc::get_domain(temp), + [=](ddc::DiscreteElement const ix) { + std::cout << std::setw(6) << temp_slice(ix); + }); + std::cout << " }" << std::endl; + +#ifdef DDC_BUILD_PDI_WRAPPER + ddc::PdiEvent("display") + .with("temp", temp) + .and_with("mean_temp", mean_temp) + .and_with("temp_slice", temp_slice); +#endif +} +//! [display] + + +//! [main-start] +//! [main-start-x-parameters] +int main(int argc, char** argv) +{ +#ifdef DDC_BUILD_PDI_WRAPPER + auto pdi_conf = PC_parse_string(""); + PDI_init(pdi_conf); +#endif + Kokkos::ScopeGuard const kokkos_scope(argc, argv); + ddc::ScopeGuard const ddc_scope(argc, argv); + + + //! [parameters] + double const x_start = -1.; + double const x_end = 1.; + std::size_t const nb_x_points = 10; + double const kx = .01; + //! [main-start-x-parameters] + //! [main-start-y-parameters] + double const y_start = -1.; + double const y_end = 1.; + std::size_t const nb_y_points = 100; + double const ky = .002; + //! [main-start-y-parameters] + //! [main-start-t-parameters] + double const start_time = 0.; + double const end_time = 10.; + //! [main-start-t-parameters] + std::ptrdiff_t const t_output_period = 10; + //! [parameters] + + //! [iterator_main-domain] + std::vector x_domain_vect + = generate_random_vector(nb_x_points, x_start, x_end); + //! [iterator_main-domain] + + std::size_t size_x = x_domain_vect.size(); + + //! [ghost_points_x] + std::vector x_pre_ghost_vect { + x_domain_vect.front() + - (x_domain_vect.back() - x_domain_vect[size_x - 2])}; + + std::vector x_post_ghost_vect { + x_domain_vect.back() + + (x_domain_vect[1] - x_domain_vect.front())}; + //! [ghost_points_x] + + //! [build-domains] + auto const [x_domain, ghosted_x_domain, x_pre_ghost, x_post_ghost] + = ddc::init_discrete_space(DDimX::init_ghosted( + x_domain_vect, + x_pre_ghost_vect, + x_post_ghost_vect)); + //! [build-domains] + + ddc::DiscreteDomain const + x_domain_begin(x_domain.front(), x_post_ghost.extents()); + ddc::DiscreteDomain const x_domain_end( + x_domain.back() - x_pre_ghost.extents() + 1, + x_pre_ghost.extents()); + + //! [Y-vectors] + std::vector y_domain_vect + = generate_random_vector(nb_y_points, y_start, y_end); + + std::size_t size_y = y_domain_vect.size(); + + //! [ghost_points_y] + std::vector y_pre_ghost_vect { + y_domain_vect.front() + - (y_domain_vect.back() - y_domain_vect[size_y - 2])}; + std::vector y_post_ghost_vect { + y_domain_vect.back() + + (y_domain_vect[1] - y_domain_vect.front())}; + //! [ghost_points_y] + + //! [Y-vectors] + + //! [build-Y-domain] + auto const [y_domain, ghosted_y_domain, y_pre_ghost, y_post_ghost] + = ddc::init_discrete_space(DDimY::init_ghosted( + y_domain_vect, + y_pre_ghost_vect, + y_post_ghost_vect)); + //! [build-Y-domain] + + ddc::DiscreteDomain const + y_domain_begin(y_domain.front(), y_post_ghost.extents()); + + ddc::DiscreteDomain const y_domain_end( + y_domain.back() - y_pre_ghost.extents() + 1, + y_pre_ghost.extents()); + + //! [CFL-condition] + + double const invdx2_max = ddc::transform_reduce( + x_domain, + 0., + ddc::reducer::max(), + [](ddc::DiscreteElement ix) { + return 1. + / (ddc::distance_at_left(ix) + * ddc::distance_at_right(ix)); + }); + + double const invdy2_max = ddc::transform_reduce( + y_domain, + 0., + ddc::reducer::max(), + [](ddc::DiscreteElement iy) { + return 1. + / (ddc::distance_at_left(iy) + * ddc::distance_at_right(iy)); + }); + + ddc::Coordinate const max_dt { + .5 / (kx * invdx2_max + ky * invdy2_max)}; + + //! [CFL-condition] + + //! [time-domain] + ddc::DiscreteVector const nb_time_steps( + std::ceil((end_time - start_time) / max_dt) + .2); + + ddc::DiscreteDomain const time_domain + = ddc::init_discrete_space(DDimT::init( + ddc::Coordinate(start_time), + ddc::Coordinate(end_time), + nb_time_steps + 1)); + + //! [time-domain] + + //! [data allocation] + ddc::Chunk ghosted_last_temp( + "ghosted_last_temp", + ddc::DiscreteDomain< + DDimX, + DDimY>(ghosted_x_domain, ghosted_y_domain), + ddc::DeviceAllocator()); + + ddc::Chunk ghosted_next_temp( + "ghosted_next_temp", + ddc::DiscreteDomain< + DDimX, + DDimY>(ghosted_x_domain, ghosted_y_domain), + ddc::DeviceAllocator()); + //! [data allocation] + + //! [initial-chunkspan] + ddc::ChunkSpan const ghosted_initial_temp + = ghosted_last_temp.span_view(); + //! [initial-chunkspan] + + //! [fill-initial-chunkspan] + ddc::parallel_for_each( + ddc::DiscreteDomain(x_domain, y_domain), + KOKKOS_LAMBDA(ddc::DiscreteElement const ixy) { + double const x = ddc::coordinate( + ddc::DiscreteElement(ixy)); + double const y = ddc::coordinate( + ddc::DiscreteElement(ixy)); + ghosted_initial_temp(ixy) + = 9.999 * ((x * x + y * y) < 0.25); + }); + //! [fill-initial-chunkspan] + + //! [host-chunk] + ddc::Chunk ghosted_temp( + "ghost_temp", + ddc::DiscreteDomain< + DDimX, + DDimY>(ghosted_x_domain, ghosted_y_domain), + ddc::HostAllocator()); + //! [host-chunk] + + //! [initial-deepcopy] + ddc::parallel_deepcopy(ghosted_temp, ghosted_last_temp); + //! [initial-deepcopy] + + //! [initial-display] + display(ddc::coordinate(time_domain.front()), + ghosted_temp[x_domain][y_domain]); + //! [initial-display] + + //! [last-output-iter] + ddc::DiscreteElement last_output_iter = time_domain.front(); + //! [last-output-iter] + + //! [time iteration] + for (ddc::DiscreteElement const iter : + time_domain.remove_first(ddc::DiscreteVector(1))) { + //! [time iteration] + + //! [boundary conditions] + ddc::parallel_deepcopy( + ghosted_last_temp[ddc::DiscreteDomain< + DDimX, + DDimY>(x_pre_ghost, y_domain)], + ghosted_last_temp[ddc::DiscreteDomain< + DDimX, + DDimY>(y_domain, x_domain_end)]); + ddc::parallel_deepcopy( + ghosted_last_temp[ddc::DiscreteDomain< + DDimX, + DDimY>(y_domain, x_post_ghost)], + ghosted_last_temp[ddc::DiscreteDomain< + DDimX, + DDimY>(y_domain, x_domain_begin)]); + ddc::parallel_deepcopy( + ghosted_last_temp[ddc::DiscreteDomain< + DDimX, + DDimY>(x_domain, y_pre_ghost)], + ghosted_last_temp[ddc::DiscreteDomain< + DDimX, + DDimY>(x_domain, y_domain_end)]); + ddc::parallel_deepcopy( + ghosted_last_temp[ddc::DiscreteDomain< + DDimX, + DDimY>(x_domain, y_post_ghost)], + ghosted_last_temp[ddc::DiscreteDomain< + DDimX, + DDimY>(x_domain, y_domain_begin)]); + //! [boundary conditions] + + //! [manipulated views] + ddc::ChunkSpan const next_temp( + ghosted_next_temp[ddc::DiscreteDomain< + DDimX, + DDimY>(x_domain, y_domain)]); + ddc::ChunkSpan const last_temp(ghosted_last_temp.span_cview()); + //! [manipulated views] + + //! [numerical scheme] + ddc::parallel_for_each( + next_temp.domain(), + KOKKOS_LAMBDA( + ddc::DiscreteElement const ixy) { + ddc::DiscreteElement const ix + = ddc::select(ixy); + ddc::DiscreteElement const iy + = ddc::select(ixy); + double const dx_l = ddc::distance_at_left(ix); + double const dx_r = ddc::distance_at_right(ix); + double const dx_m = 0.5 * (dx_l + dx_r); + double const dy_l = ddc::distance_at_left(iy); + double const dy_r = ddc::distance_at_right(iy); + double const dy_m = 0.5 * (dy_l + dy_r); + next_temp(ix, iy) = last_temp(ix, iy); + next_temp(ix, iy) + += kx * ddc::step() + * (dx_l * last_temp(ix + 1, iy) + - 2.0 * dx_m * last_temp(ix, iy) + + dx_r * last_temp(ix - 1, iy)) + / (dx_l * dx_m * dx_r); + next_temp(ix, iy) + += ky * ddc::step() + * (dy_l * last_temp(ix, iy + 1) + - 2.0 * dy_m * last_temp(ix, iy) + + dy_r * last_temp(ix, iy - 1)) + / (dy_l * dy_m * dy_r); + }); + //! [numerical scheme] + + //! [output] + if (iter - last_output_iter >= t_output_period) { + last_output_iter = iter; + ddc::parallel_deepcopy(ghosted_temp, ghosted_last_temp); + display(ddc::coordinate(iter), + ghosted_temp[ddc::DiscreteDomain< + DDimX, + DDimY>(x_domain, y_domain)]); + } + //! [output] + + //! [swap] + std::swap(ghosted_last_temp, ghosted_next_temp); + //! [swap] + } + + + //! [final output] + if (last_output_iter < time_domain.back()) { + ddc::parallel_deepcopy(ghosted_temp, ghosted_last_temp); + display(ddc::coordinate(time_domain.back()), + ghosted_temp[ddc::DiscreteDomain< + DDimX, + DDimY>(x_domain, y_domain)]); + } + //! [final output] + +#ifdef DDC_BUILD_PDI_WRAPPER + PDI_finalize(); + PC_tree_destroy(&pdi_conf); +#endif +} diff --git a/examples/uniform_heat_equation.cpp b/examples/uniform_heat_equation.cpp index 91ea4623e..53c20d649 100644 --- a/examples/uniform_heat_equation.cpp +++ b/examples/uniform_heat_equation.cpp @@ -36,13 +36,17 @@ struct DDimT : ddc::UniformPointSampling }; //! [time-space] +//! [display] + /** A function to pretty print the temperature * @tparam ChunkType The type of chunk span. This way the template parameters are avoided, * should be deduced by the compiler. * @param time The time at which the output is made. * @param temp The temperature at this time-step. */ + //! [display] + template void display(double time, ChunkType temp) { @@ -107,6 +111,7 @@ int main(int argc, char** argv) //! [main-start-t-parameters] std::ptrdiff_t const t_output_period = 10; //! [parameters] + //! [main-start] //! [X-parameters] @@ -149,12 +154,14 @@ int main(int argc, char** argv) //! [Y-domains] //! [CFL-condition] + double const dx = ddc::step(); double const dy = ddc::step(); double const invdx2 = 1. / (dx * dx); double const invdy2 = 1. / (dy * dy); ddc::Coordinate const dt(.5 / (kx * invdx2 + ky * invdy2)); + //! [CFL-condition] //! [time-domain] @@ -166,6 +173,7 @@ int main(int argc, char** argv) ddc::Coordinate(start_time), ddc::Coordinate(end_time), nb_time_steps + 1)); + //! [time-domain] //! [data allocation] diff --git a/include/ddc/kernels/splines/splines_linear_problem_sparse.hpp b/include/ddc/kernels/splines/splines_linear_problem_sparse.hpp index 209e3d247..38ba8d28c 100644 --- a/include/ddc/kernels/splines/splines_linear_problem_sparse.hpp +++ b/include/ddc/kernels/splines/splines_linear_problem_sparse.hpp @@ -192,7 +192,7 @@ class SplinesLinearProblemSparse : public SplinesLinearProblem * * Removes the zeros from the CSR object and instantiate a Ginkgo solver. It also constructs a transposed version of the solver. * - * The stopping criterion is a reduction factor ||x||/||b||<1e-15 with 1000 maximum iterations. + * The stopping criterion is a reduction factor ||Ax-b||/||b||<1e-15 with 1000 maximum iterations. */ void setup_solver() override { diff --git a/include/ddc/non_uniform_point_sampling.hpp b/include/ddc/non_uniform_point_sampling.hpp index 7bb957b4e..7a97d6708 100644 --- a/include/ddc/non_uniform_point_sampling.hpp +++ b/include/ddc/non_uniform_point_sampling.hpp @@ -113,6 +113,12 @@ class NonUniformPointSampling : detail::NonUniformPointSamplingBase return m_points.size(); } + /// @brief Lower bound index of the mesh + KOKKOS_FUNCTION discrete_element_type front() const noexcept + { + return discrete_element_type {0}; + } + /// @brief Convert a mesh index into a position in `CDim` KOKKOS_FUNCTION continuous_element_type coordinate(discrete_element_type const& icoord) const noexcept @@ -120,6 +126,78 @@ class NonUniformPointSampling : detail::NonUniformPointSamplingBase return m_points(icoord.uid()); } }; + + /** Construct an Impl and associated discrete_domain_type from an iterator + * containing the points coordinates along the `DDim` dimension. + * + * @param non_uniform_points a vector containing the coordinates of the points of the domain. + */ + template + static std::tuple, DiscreteDomain> + init(InputRange const non_uniform_points) + { + auto a = non_uniform_points.begin(); + auto b = non_uniform_points.end(); + auto n = std::distance(non_uniform_points.begin(), non_uniform_points.end()); + assert(a < b); + assert(n > 0); + typename DDim::template Impl disc(non_uniform_points); + DiscreteDomain domain {disc.front(), DiscreteVector {n}}; + return std::make_tuple(std::move(disc), std::move(domain)); + } + + /** Construct 4 non-uniform `DiscreteDomain` and an Impl from 3 iterators containing the points coordinates along the `DDim` dimension. + * + * @param domain_r an iterator containing the coordinates of the points of the main domain along the DDim position + * @param pre_ghost_r an iterator containing the positions of the ghost points before the main domain the DDim position + * @param post_ghost_r an iterator containing the positions of the ghost points after the main domain the DDim position + */ + template + static std::tuple< + typename DDim::template Impl, + DiscreteDomain, + DiscreteDomain, + DiscreteDomain, + DiscreteDomain> + init_ghosted( + InputRange const& domain_r, + InputRange const& pre_ghost_r, + InputRange const& post_ghost_r) + { + using discrete_domain_type = DiscreteDomain; + auto n = DiscreteVector {std::distance(domain_r.begin(), domain_r.end())}; + + assert(domain_r.begin() < domain_r.end()); + assert(n > 1); + + auto n_ghosts_before + = DiscreteVector {std::distance(pre_ghost_r.begin(), pre_ghost_r.end())}; + auto n_ghosts_after + = DiscreteVector {std::distance(post_ghost_r.begin(), post_ghost_r.end())}; + + std::vector full_domain; + + std::copy(pre_ghost_r.begin(), pre_ghost_r.end(), std::back_inserter(full_domain)); + std::copy(domain_r.begin(), domain_r.end(), std::back_inserter(full_domain)); + std::copy(post_ghost_r.begin(), post_ghost_r.end(), std::back_inserter(full_domain)); + + typename DDim::template Impl disc(full_domain); + + discrete_domain_type ghosted_domain + = discrete_domain_type(disc.front(), n + n_ghosts_before + n_ghosts_after); + discrete_domain_type pre_ghost + = discrete_domain_type(ghosted_domain.front(), n_ghosts_before); + discrete_domain_type main_domain + = discrete_domain_type(ghosted_domain.front() + n_ghosts_before, n); + discrete_domain_type post_ghost + = discrete_domain_type(main_domain.back() + 1, n_ghosts_after); + return std::make_tuple( + std::move(disc), + std::move(main_domain), + std::move(ghosted_domain), + std::move(pre_ghost), + std::move(post_ghost)); + } }; template