2D Poisson example

examples/2d_poisson.cpp

@@ -0,0 +1,333 @@
+// SPDX-FileCopyrightText: 2024 Baptiste Legouix
+// SPDX-License-Identifier: GPL-3.0-or-later
+
+#include <ddc/ddc.hpp>
+#include <ddc/kernels/splines.hpp>
+#include <ddc/pdi.hpp>
+
+#include <similie/similie.hpp>
+
+// 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(<?xml version="1.0" ?>
+<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>
+<Xdmf Version="2.0">
+ <Domain>
+   <Grid Name="mesh1" GridType="Uniform">
+     <Topology TopologyType="2DSMesh" NumberOfElements="%i %i"/>
+     <Geometry GeometryType="XY">
+       <DataItem Dimensions="%i %i 2" NumberType="Float" Precision="8" Format="HDF">
+        2d_poisson.h5:/position
+       </DataItem>
+     </Geometry>
+     <Attribute Name="Source" AttributeType="Scalar" Center="Cell">
+       <DataItem Dimensions="%i %i" NumberType="Float" Precision="8" Format="HDF">
+        2d_poisson.h5:/source
+       </DataItem>
+     </Attribute>
+     <Attribute Name="Potential" AttributeType="Scalar" Center="Node">
+       <DataItem Dimensions="%i %i" NumberType="Float" Precision="8" Format="HDF">
+        2d_poisson.h5:/potential
+       </DataItem>
+     </Attribute>
+     <Attribute Name="Laplacian" AttributeType="Scalar" Center="Cell">
+       <DataItem Dimensions="%i %i" NumberType="Float" Precision="8" Format="HDF">
+        2d_poisson.h5:/laplacian
+       </DataItem>
+     </Attribute>
+   </Grid>
+ </Domain>
+</Xdmf>
+)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::MetricIndex1<X, Y>>,
+        sil::tensor::Covariant<sil::tensor::MetricIndex2<X, Y>>>;
+
+using MesherXY = sil::mesher::Mesher<s_degree, X, Y>;
+
+struct BSplinesX : MesherXY::template bsplines_type<X>
+{
+};
+
+struct DDimX : MesherXY::template discrete_dimension_type<X>
+{
+};
+
+struct BSplinesY : MesherXY::template bsplines_type<Y>
+{
+};
+
+struct DDimY : MesherXY::template discrete_dimension_type<Y>
+{
+};
+
+// Declare natural indices taking values in {X, Y}
+struct Nu : sil::tensor::TensorNaturalIndex<X, Y>
+{
+};
+
+// Declare indices
+using NuLow = sil::tensor::Covariant<Nu>;
+using NuUp = sil::tensor::Contravariant<Nu>;
+
+using DummyIndex = sil::tensor::Covariant<sil::tensor::ScalarIndex>;
+
+// 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<X, Y> lower_bounds(-5., -5.);
+    ddc::Coordinate<X, Y> upper_bounds(5., 5.);
+    ddc::DiscreteVector<DDimX, DDimY> nb_cells(1000, 1000);
+    ddc::DiscreteDomain<DDimX, DDimY> mesh_xy = mesher.template mesh<
+            ddc::detail::TypeSeq<DDimX, DDimY>,
+            ddc::detail::TypeSeq<BSplinesX, BSplinesY>>(lower_bounds, upper_bounds, nb_cells);
+    assert(static_cast<std::size_t>(mesh_xy.template extent<DDimX>())
+           == static_cast<std::size_t>(mesh_xy.template extent<DDimY>()));
+    ddc::expose_to_pdi("Nx", static_cast<int>(mesh_xy.template extent<DDimX>()));
+    ddc::expose_to_pdi("Ny", static_cast<int>(mesh_xy.template extent<DDimY>()));
+
+    // Allocate and instantiate a position field (used only to be exported).
+    [[maybe_unused]] sil::tensor::TensorAccessor<NuUp> position_accessor;
+    ddc::DiscreteDomain<DDimX, DDimY, NuUp> position_dom(mesh_xy, position_accessor.domain());
+    ddc::Chunk position_alloc(position_dom, ddc::HostAllocator<double>());
+    sil::tensor::Tensor position(position_alloc);
+    ddc::parallel_for_each(
+            Kokkos::DefaultHostExecutionSpace(),
+            mesh_xy,
+            KOKKOS_LAMBDA(ddc::DiscreteElement<DDimX, DDimY> elem) {
+                position(elem, position.accessor().access_element<X>())
+                        = static_cast<double>(ddc::coordinate(ddc::DiscreteElement<DDimX>(elem)));
+                position(elem, position.accessor().access_element<Y>())
+                        = static_cast<double>(ddc::coordinate(ddc::DiscreteElement<DDimY>(elem)));
+            });
+
+    // Allocate and instantiate a metric tensor field.
+    [[maybe_unused]] sil::tensor::TensorAccessor<MetricIndex> metric_accessor;
+    ddc::DiscreteDomain<DDimX, DDimY, MetricIndex> metric_dom(mesh_xy, metric_accessor.domain());
+    ddc::Chunk metric_alloc(metric_dom, ddc::DeviceAllocator<double>());
+    sil::tensor::Tensor metric(metric_alloc);
+    ddc::parallel_for_each(
+            Kokkos::DefaultExecutionSpace(),
+            mesh_xy,
+            KOKKOS_LAMBDA(ddc::DiscreteElement<DDimX, DDimY> elem) {
+                metric(elem, metric.accessor().access_element<X, X>()) = 1.;
+                metric(elem, metric.accessor().access_element<X, Y>()) = 0.;
+                metric(elem, metric.accessor().access_element<Y, Y>()) = 1.;
+            });
+
+    // Invert metric
+    [[maybe_unused]] sil::tensor::TensorAccessor<sil::tensor::upper_t<MetricIndex>>
+            inv_metric_accessor;
+    ddc::DiscreteDomain<DDimX, DDimY, sil::tensor::upper_t<MetricIndex>>
+            inv_metric_dom(mesh_xy, inv_metric_accessor.domain());
+    ddc::Chunk inv_metric_alloc(inv_metric_dom, ddc::DeviceAllocator<double>());
+    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<DummyIndex> source_accessor;
+    ddc::DiscreteDomain<DDimX, DDimY, DummyIndex> source_dom(
+            mesh_xy.remove_last(ddc::DiscreteVector<DDimX, DDimY> {1, 1}),
+            source_accessor.domain());
+    ddc::Chunk source_alloc(source_dom, ddc::DeviceAllocator<double>());
+    sil::tensor::Tensor source(source_alloc);
+
+    double const R = 2.;
+    ddc::parallel_for_each(
+            source.domain(),
+            KOKKOS_LAMBDA(ddc::DiscreteElement<DDimX, DDimY, DummyIndex> elem) {
+                double const r = Kokkos::sqrt(
+                        static_cast<double>(
+                                ddc::coordinate(ddc::DiscreteElement<DDimX>(elem))
+                                * ddc::coordinate(ddc::DiscreteElement<DDimX>(elem)))
+                        + static_cast<double>(
+                                ddc::coordinate(ddc::DiscreteElement<DDimY>(elem))
+                                * ddc::coordinate(ddc::DiscreteElement<DDimY>(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<DummyIndex> potential_accessor;
+    ddc::DiscreteDomain<DDimX, DDimY, DummyIndex>
+            potential_dom(metric.non_indices_domain(), potential_accessor.domain());
+    ddc::Chunk potential_alloc(potential_dom, ddc::DeviceAllocator<double>());
+    sil::tensor::Tensor potential(potential_alloc);
+
+    double const L = ddc::coordinate(ddc::DiscreteElement<DDimX>(potential.domain().back()))
+                     - ddc::coordinate(ddc::DiscreteElement<DDimX>(potential.domain().front()));
+    double const alpha = (static_cast<double>(nb_cells.template get<DDimX>())
+                          * static_cast<double>(nb_cells.template get<DDimY>()))
+                         / L / 2 / L / 2;
+    ddc::parallel_for_each(
+            potential.domain(),
+            KOKKOS_LAMBDA(ddc::DiscreteElement<DDimX, DDimY, DummyIndex> elem) {
+                double const r = Kokkos::sqrt(
+                        static_cast<double>(
+                                ddc::coordinate(ddc::DiscreteElement<DDimX>(elem))
+                                * ddc::coordinate(ddc::DiscreteElement<DDimX>(elem)))
+                        + static_cast<double>(
+                                ddc::coordinate(ddc::DiscreteElement<DDimY>(elem))
+                                * ddc::coordinate(ddc::DiscreteElement<DDimY>(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<DummyIndex> laplacian_accessor;
+    ddc::DiscreteDomain<DDimX, DDimY, DummyIndex> laplacian_dom(
+            mesh_xy.remove_last(ddc::DiscreteVector<DDimX, DDimY> {1, 1}),
+            laplacian_accessor.domain());
+    ddc::Chunk laplacian_alloc(laplacian_dom, ddc::DeviceAllocator<double>());
+    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<int>(mesh_xy.template extent<DDimX>()),
+            static_cast<int>(mesh_xy.template extent<DDimY>()));
+    std::cout << "XDMF model exported in 2d_poisson.xmf." << std::endl;
+
+    // Finalize PDI
+    PC_tree_destroy(&conf_pdi);
+    PDI_finalize();
+
+    return EXIT_SUCCESS;
+}

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)