diff --git a/examples/2d_poisson.cpp b/examples/2d_poisson.cpp new file mode 100644 index 00000000..3462f4c5 --- /dev/null +++ b/examples/2d_poisson.cpp @@ -0,0 +1,333 @@ +// SPDX-FileCopyrightText: 2024 Baptiste Legouix +// SPDX-License-Identifier: GPL-3.0-or-later + +#include +#include +#include + +#include + +// PDI config +constexpr char const* const PDI_CFG = R"PDI_CFG( +metadata: + Nx : int + Ny : int + +data: + position: + type: array + subtype: double + size: [ '$Nx', '$Ny', 2] + source: + type: array + subtype: double + size: [ '$Nx-1', '$Ny-1' ] + potential: + type: array + subtype: double + size: [ '$Nx', '$Ny' ] + laplacian: + type: array + subtype: double + size: [ '$Nx-1', '$Ny-1' ] + +plugins: + decl_hdf5: + - file: '2d_poisson.h5' + on_event: [export] + collision_policy: replace_and_warn + write: [Nx, Ny, position, source, potential, laplacian] + #trace: ~ +)PDI_CFG"; + +// XDMF +int write_xdmf(int Nx, int Ny) +{ + constexpr char const* const xdmf = R"XDMF( + + + + + + + + 2d_poisson.h5:/position + + + + + 2d_poisson.h5:/source + + + + + 2d_poisson.h5:/potential + + + + + 2d_poisson.h5:/laplacian + + + + + +)XDMF"; + + FILE* file = fopen("2d_poisson.xmf", "w"); + fprintf(file, xdmf, Nx, Ny, Nx, Ny, Nx - 1, Ny - 1, Nx, Ny, Nx - 1, Ny - 1); + fclose(file); + + return 1; +} + +static constexpr std::size_t s_degree = 3; + +// Labelize the dimensions of space +struct X +{ + static constexpr bool PERIODIC = false; +}; + +struct Y +{ + static constexpr bool PERIODIC = false; +}; + +// Declare a metric +using MetricIndex = sil::tensor::TensorSymmetricIndex< + sil::tensor::Covariant>, + sil::tensor::Covariant>>; + +using MesherXY = sil::mesher::Mesher; + +struct BSplinesX : MesherXY::template bsplines_type +{ +}; + +struct DDimX : MesherXY::template discrete_dimension_type +{ +}; + +struct BSplinesY : MesherXY::template bsplines_type +{ +}; + +struct DDimY : MesherXY::template discrete_dimension_type +{ +}; + +// Declare natural indices taking values in {X, Y} +struct Nu : sil::tensor::TensorNaturalIndex +{ +}; + +// Declare indices +using NuLow = sil::tensor::Covariant; +using NuUp = sil::tensor::Contravariant; + +using DummyIndex = sil::tensor::Covariant; + +// Declare boundary conditions (Dirichlet) +static constexpr ddc::BoundCond BoundCond = ddc::BoundCond::GREVILLE; +using ExtrapolationRule = ddc::NullExtrapolationRule; + +int main(int argc, char** argv) +{ + // Initialize PDI, Kokkos and DDC + PC_tree_t conf_pdi = PC_parse_string(PDI_CFG); + PC_errhandler(PC_NULL_HANDLER); + PDI_init(conf_pdi); + + Kokkos::ScopeGuard const kokkos_scope(argc, argv); + ddc::ScopeGuard const ddc_scope(argc, argv); + + // Produce mesh + MesherXY mesher; + ddc::Coordinate lower_bounds(-5., -5.); + ddc::Coordinate upper_bounds(5., 5.); + ddc::DiscreteVector nb_cells(1000, 1000); + ddc::DiscreteDomain mesh_xy = mesher.template mesh< + ddc::detail::TypeSeq, + ddc::detail::TypeSeq>(lower_bounds, upper_bounds, nb_cells); + assert(static_cast(mesh_xy.template extent()) + == static_cast(mesh_xy.template extent())); + ddc::expose_to_pdi("Nx", static_cast(mesh_xy.template extent())); + ddc::expose_to_pdi("Ny", static_cast(mesh_xy.template extent())); + + // Allocate and instantiate a position field (used only to be exported). + [[maybe_unused]] sil::tensor::TensorAccessor position_accessor; + ddc::DiscreteDomain position_dom(mesh_xy, position_accessor.domain()); + ddc::Chunk position_alloc(position_dom, ddc::HostAllocator()); + sil::tensor::Tensor position(position_alloc); + ddc::parallel_for_each( + Kokkos::DefaultHostExecutionSpace(), + mesh_xy, + KOKKOS_LAMBDA(ddc::DiscreteElement elem) { + position(elem, position.accessor().access_element()) + = static_cast(ddc::coordinate(ddc::DiscreteElement(elem))); + position(elem, position.accessor().access_element()) + = static_cast(ddc::coordinate(ddc::DiscreteElement(elem))); + }); + + // Allocate and instantiate a metric tensor field. + [[maybe_unused]] sil::tensor::TensorAccessor metric_accessor; + ddc::DiscreteDomain metric_dom(mesh_xy, metric_accessor.domain()); + ddc::Chunk metric_alloc(metric_dom, ddc::DeviceAllocator()); + sil::tensor::Tensor metric(metric_alloc); + ddc::parallel_for_each( + Kokkos::DefaultExecutionSpace(), + mesh_xy, + KOKKOS_LAMBDA(ddc::DiscreteElement elem) { + metric(elem, metric.accessor().access_element()) = 1.; + metric(elem, metric.accessor().access_element()) = 0.; + metric(elem, metric.accessor().access_element()) = 1.; + }); + + // Invert metric + [[maybe_unused]] sil::tensor::TensorAccessor> + inv_metric_accessor; + ddc::DiscreteDomain> + inv_metric_dom(mesh_xy, inv_metric_accessor.domain()); + ddc::Chunk inv_metric_alloc(inv_metric_dom, ddc::DeviceAllocator()); + sil::tensor::Tensor inv_metric(inv_metric_alloc); + sil::tensor::fill_inverse_metric< + MetricIndex>(Kokkos::DefaultExecutionSpace(), inv_metric, metric); + + // Source + [[maybe_unused]] sil::tensor::TensorAccessor source_accessor; + ddc::DiscreteDomain source_dom( + mesh_xy.remove_last(ddc::DiscreteVector {1, 1}), + source_accessor.domain()); + ddc::Chunk source_alloc(source_dom, ddc::DeviceAllocator()); + sil::tensor::Tensor source(source_alloc); + + double const R = 2.; + ddc::parallel_for_each( + source.domain(), + KOKKOS_LAMBDA(ddc::DiscreteElement elem) { + double const r = Kokkos::sqrt( + static_cast( + ddc::coordinate(ddc::DiscreteElement(elem)) + * ddc::coordinate(ddc::DiscreteElement(elem))) + + static_cast( + ddc::coordinate(ddc::DiscreteElement(elem)) + * ddc::coordinate(ddc::DiscreteElement(elem)))); + if (r <= R) { + source.mem(elem) = 1.; + } else { + source.mem(elem) = 0.; + } + }); + + auto source_host + = ddc::create_mirror_view_and_copy(Kokkos::DefaultHostExecutionSpace(), source); + + // Spline builder & evaluator + ddc::SplineBuilder< + Kokkos::DefaultExecutionSpace, + Kokkos::DefaultExecutionSpace::memory_space, + BSplinesX, + DDimX, + BoundCond, + BoundCond, + ddc::SplineSolver::LAPACK, + DDimX, + DDimY> const spline_builder_x(mesh_xy); + ddc::SplineBuilder< + Kokkos::DefaultExecutionSpace, + Kokkos::DefaultExecutionSpace::memory_space, + BSplinesY, + DDimY, + BoundCond, + BoundCond, + ddc::SplineSolver::LAPACK, + DDimX, + DDimY> const spline_builder_y(mesh_xy); + ExtrapolationRule const extrapolation_rule; + ddc::SplineEvaluator< + Kokkos::DefaultExecutionSpace, + Kokkos::DefaultExecutionSpace::memory_space, + BSplinesX, + DDimX, + ExtrapolationRule, + ExtrapolationRule, + DDimX, + DDimY> const spline_evaluator_x(extrapolation_rule, extrapolation_rule); + ddc::SplineEvaluator< + Kokkos::DefaultExecutionSpace, + Kokkos::DefaultExecutionSpace::memory_space, + BSplinesY, + DDimY, + ExtrapolationRule, + ExtrapolationRule, + DDimX, + DDimY> const spline_evaluator_y(extrapolation_rule, extrapolation_rule); + + // Potential + [[maybe_unused]] sil::tensor::TensorAccessor potential_accessor; + ddc::DiscreteDomain + potential_dom(metric.non_indices_domain(), potential_accessor.domain()); + ddc::Chunk potential_alloc(potential_dom, ddc::DeviceAllocator()); + sil::tensor::Tensor potential(potential_alloc); + + double const L = ddc::coordinate(ddc::DiscreteElement(potential.domain().back())) + - ddc::coordinate(ddc::DiscreteElement(potential.domain().front())); + double const alpha = (static_cast(nb_cells.template get()) + * static_cast(nb_cells.template get())) + / L / 2 / L / 2; + ddc::parallel_for_each( + potential.domain(), + KOKKOS_LAMBDA(ddc::DiscreteElement elem) { + double const r = Kokkos::sqrt( + static_cast( + ddc::coordinate(ddc::DiscreteElement(elem)) + * ddc::coordinate(ddc::DiscreteElement(elem))) + + static_cast( + ddc::coordinate(ddc::DiscreteElement(elem)) + * ddc::coordinate(ddc::DiscreteElement(elem)))); + if (r <= R) { + potential.mem(elem) = -alpha * r * r; + } else { + potential.mem(elem) = alpha * R * R * (2 * Kokkos::log(R / r) - 1); + } + }); + auto potential_host + = ddc::create_mirror_view_and_copy(Kokkos::DefaultHostExecutionSpace(), potential); + + // Laplacian + [[maybe_unused]] sil::tensor::TensorAccessor laplacian_accessor; + ddc::DiscreteDomain laplacian_dom( + mesh_xy.remove_last(ddc::DiscreteVector {1, 1}), + laplacian_accessor.domain()); + ddc::Chunk laplacian_alloc(laplacian_dom, ddc::DeviceAllocator()); + sil::tensor::Tensor laplacian(laplacian_alloc); + + sil::exterior::laplacian< + MetricIndex, + NuLow, + DummyIndex>(Kokkos::DefaultExecutionSpace(), laplacian, potential, inv_metric); + Kokkos::fence(); + + auto laplacian_host + = ddc::create_mirror_view_and_copy(Kokkos::DefaultHostExecutionSpace(), laplacian); + + // Export HDF5 and XDMF + ddc::PdiEvent("export") + .with("position", position) + .and_with("source", source_host) + .and_with("potential", potential_host) + .and_with("laplacian", laplacian_host); + std::cout << "Computation result exported in 2d_poisson.h5." << std::endl; + + write_xdmf( + static_cast(mesh_xy.template extent()), + static_cast(mesh_xy.template extent())); + std::cout << "XDMF model exported in 2d_poisson.xmf." << std::endl; + + // Finalize PDI + PC_tree_destroy(&conf_pdi); + PDI_finalize(); + + return EXIT_SUCCESS; +} diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 812c458c..fc475336 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -17,6 +17,14 @@ target_link_libraries(2d_vector_laplacian sil::sil ) +add_executable(2d_poisson 2d_poisson.cpp) + +target_link_libraries(2d_poisson + PUBLIC + DDC::DDC + sil::sil +) + if("${SIMILIE_BUILD_YOUNG_TABLEAU}") add_executable(kretschmann_scalar kretschmann_scalar.cpp)