From 2671f71fe731753c9b6afaeddf0ea9066d95c8ed Mon Sep 17 00:00:00 2001 From: blegouix Date: Sat, 28 Mar 2026 23:37:04 +0100 Subject: [PATCH 01/14] wip --- .../2d_free_scalar_field.cpp | 186 +++++++++++++++--- src/similie/exterior/exterior.hpp | 1 + src/similie/exterior/staggered_difference.hpp | 180 +++++++++++++++++ 3 files changed, 337 insertions(+), 30 deletions(-) create mode 100644 src/similie/exterior/staggered_difference.hpp diff --git a/examples/free_scalar_field/2d_free_scalar_field.cpp b/examples/free_scalar_field/2d_free_scalar_field.cpp index e3c1828f..96535f10 100644 --- a/examples/free_scalar_field/2d_free_scalar_field.cpp +++ b/examples/free_scalar_field/2d_free_scalar_field.cpp @@ -238,6 +238,14 @@ struct DDimY : ddc::UniformPointSampling { }; +struct DDimXHalf : ddc::UniformPointSampling +{ +}; + +struct DDimYHalf : ddc::UniformPointSampling +{ +}; + /* struct BSplinesX : MesherXY::template bsplines_type { @@ -305,7 +313,29 @@ int main(int argc, char** argv) ddc::select(lower_bounds), ddc::select(upper_bounds), ddc::select(nb_cells))); + auto const x_face_x_dom = ddc::init_discrete_space(DDimXHalf::init( + ddc::Coordinate( + ddc::select(lower_bounds) + + (ddc::select(upper_bounds) - ddc::select(lower_bounds)) + / (2 * (ddc::select(nb_cells).value() - 1))), + ddc::Coordinate( + ddc::select(upper_bounds) + - (ddc::select(upper_bounds) - ddc::select(lower_bounds)) + / (2 * (ddc::select(nb_cells).value() - 1))), + ddc::DiscreteVector(ddc::select(nb_cells).value() - 1))); + auto const y_face_y_dom = ddc::init_discrete_space(DDimYHalf::init( + ddc::Coordinate( + ddc::select(lower_bounds) + + (ddc::select(upper_bounds) - ddc::select(lower_bounds)) + / (2 * (ddc::select(nb_cells).value() - 1))), + ddc::Coordinate( + ddc::select(upper_bounds) + - (ddc::select(upper_bounds) - ddc::select(lower_bounds)) + / (2 * (ddc::select(nb_cells).value() - 1))), + ddc::DiscreteVector(ddc::select(nb_cells).value() - 1))); ddc::DiscreteDomain mesh_xy(x_dom, y_dom); + ddc::DiscreteDomain x_face_dom(x_face_x_dom, y_dom); + ddc::DiscreteDomain y_face_dom(x_dom, y_face_y_dom); assert(static_cast(mesh_xy.template extent()) == static_cast(mesh_xy.template extent())); @@ -395,17 +425,40 @@ int main(int argc, char** argv) ddc::Chunk half_step_potential_alloc(potential_dom, ddc::DeviceAllocator()); sil::tensor::Tensor half_step_potential(half_step_potential_alloc); - // Potential gradient - [[maybe_unused]] sil::tensor::TensorAccessor potential_grad_accessor; + // Staggered potential gradients + [[maybe_unused]] sil::tensor::TensorAccessor scalar_accessor; + ddc::DiscreteDomain + potential_grad_x_dom(x_face_dom, scalar_accessor.domain()); + ddc::Chunk potential_grad_x_alloc(potential_grad_x_dom, ddc::DeviceAllocator()); + sil::tensor::Tensor potential_grad_x(potential_grad_x_alloc); + + ddc::DiscreteDomain + potential_grad_y_dom(y_face_dom, scalar_accessor.domain()); + ddc::Chunk potential_grad_y_alloc(potential_grad_y_dom, ddc::DeviceAllocator()); + sil::tensor::Tensor potential_grad_y(potential_grad_y_alloc); + + // Staggered spatial moments + ddc::DiscreteDomain + spatial_moment_x_dom(x_face_dom, scalar_accessor.domain()); + ddc::Chunk spatial_moment_x_alloc(spatial_moment_x_dom, ddc::DeviceAllocator()); + sil::tensor::Tensor spatial_moment_x(spatial_moment_x_alloc); + + ddc::DiscreteDomain + spatial_moment_y_dom(y_face_dom, scalar_accessor.domain()); + ddc::Chunk spatial_moment_y_alloc(spatial_moment_y_dom, ddc::DeviceAllocator()); + sil::tensor::Tensor spatial_moment_y(spatial_moment_y_alloc); + + auto spatial_moment_x_host = ddc:: + create_mirror_view_and_copy(Kokkos::DefaultHostExecutionSpace(), spatial_moment_x); + auto spatial_moment_y_host = ddc:: + create_mirror_view_and_copy(Kokkos::DefaultHostExecutionSpace(), spatial_moment_y); + + // Reconstructed spatial moments on nodes, only for export and diagnostics + [[maybe_unused]] sil::tensor::TensorAccessor spatial_moments_accessor; ddc::DiscreteDomain - potential_grad_dom(mesh_xy, potential_grad_accessor.domain()); - ddc::Chunk potential_grad_alloc(potential_grad_dom, ddc::DeviceAllocator()); - sil::tensor::Tensor potential_grad(potential_grad_alloc); - - // Spatial moments - auto& spatial_moments - = potential_grad; // We can perform the computations inplace so spatial_moments is just an alias of potential_grad - + spatial_moments_dom(mesh_xy, spatial_moments_accessor.domain()); + ddc::Chunk spatial_moments_alloc(spatial_moments_dom, ddc::DeviceAllocator()); + sil::tensor::Tensor spatial_moments(spatial_moments_alloc); auto spatial_moments_host = ddc:: create_mirror_view_and_copy(Kokkos::DefaultHostExecutionSpace(), spatial_moments); @@ -495,32 +548,50 @@ int main(int argc, char** argv) / 2.; }); - // Compute the potential gradient - sil::exterior::deriv< - AlphaLow, - DummyIndex>(Kokkos::DefaultExecutionSpace(), potential_grad, half_step_potential); + // Compute the staggered potential gradients + sil::exterior::staggered_deriv_component( + Kokkos::DefaultExecutionSpace(), + potential_grad_x, + half_step_potential); + sil::exterior::staggered_deriv_component( + Kokkos::DefaultExecutionSpace(), + potential_grad_y, + half_step_potential); // Compute the spatial moments pi_\alpha by solving dphi/dx^\alpha = dH/dpi_\alpha ddc::parallel_for_each( Kokkos::DefaultExecutionSpace(), - mesh_xy, - KOKKOS_LAMBDA(ddc::DiscreteElement elem) { - spatial_moments(elem, ddc::DiscreteElement(0)) + x_face_dom, + KOKKOS_LAMBDA(ddc::DiscreteElement elem) { + spatial_moment_x(elem, ddc::DiscreteElement(0)) = FreeScalarFieldHamiltonian(mass).pi1( - potential_grad(elem, ddc::DiscreteElement(0))); - spatial_moments(elem, ddc::DiscreteElement(1)) + potential_grad_x(elem, ddc::DiscreteElement(0))); + }); + ddc::parallel_for_each( + Kokkos::DefaultExecutionSpace(), + y_face_dom, + KOKKOS_LAMBDA(ddc::DiscreteElement elem) { + spatial_moment_y(elem, ddc::DiscreteElement(0)) = FreeScalarFieldHamiltonian(mass).pi2( - potential_grad(elem, ddc::DiscreteElement(1))); + potential_grad_y(elem, ddc::DiscreteElement(0))); }); - if (i % nb_iter_between_exports == 0) { - ddc::parallel_deepcopy(spatial_moments_host, spatial_moments); - } // Compute the divergence dpi_\alpha/dx^\alpha of the spatial moments, which is the codifferential \delta pi of the spatial moments - sil::exterior::codifferential( + ddc::parallel_for_each( + Kokkos::DefaultExecutionSpace(), + spatial_moments_div.domain(), + KOKKOS_LAMBDA(ddc::DiscreteElement elem) { + spatial_moments_div(elem) = 0.; + }); + sil::exterior::add_staggered_codifferential_component( + Kokkos::DefaultExecutionSpace(), + spatial_moments_div, + spatial_moment_x, + inv_metric); + sil::exterior::add_staggered_codifferential_component( Kokkos::DefaultExecutionSpace(), spatial_moments_div, - spatial_moments, + spatial_moment_y, inv_metric); // Compute dpi_0/dx^0 by solving - dpi_0/dx^0 + dpi_\alpha/dx^\alpha = -dH/dphi and advect pi_0 by a time step dx^0. Also Then, perform the second phi half-advection by solving dphi/dx^0 = -dH/dpi_0 @@ -538,17 +609,17 @@ int main(int argc, char** argv) // Advect temporal moment by half-step, this is what is needed to perform the whole-step potential advection const double half_step_temporal_moment_ = temporal_moment_ - + (FreeScalarFieldHamiltonian(mass).dH_dphi(half_step_potential_) - + spatial_moments_div_) + + (spatial_moments_div_ + - FreeScalarFieldHamiltonian(mass).dH_dphi(half_step_potential_)) * dt / 2; // Whole-step advection of field state potential(elem) - += FreeScalarFieldHamiltonian(mass).dH_dpi0(half_step_temporal_moment_) + -= FreeScalarFieldHamiltonian(mass).dH_dpi0(half_step_temporal_moment_) * dt; temporal_moment(elem) - += (FreeScalarFieldHamiltonian(mass).dH_dphi(half_step_potential_) - + spatial_moments_div_) + += (spatial_moments_div_ + - FreeScalarFieldHamiltonian(mass).dH_dphi(half_step_potential_)) * dt; }); @@ -556,6 +627,61 @@ int main(int argc, char** argv) ddc::parallel_deepcopy(temporal_moment_host, temporal_moment); ddc::parallel_deepcopy(potential_host, potential); ddc::parallel_deepcopy(h_spatial_moments_div, spatial_moments_div); + ddc::parallel_deepcopy(spatial_moment_x_host, spatial_moment_x); + ddc::parallel_deepcopy(spatial_moment_y_host, spatial_moment_y); + + { + ddc::DiscreteElement const mesh_front = mesh_xy.front(); + for (std::size_t ix = 0; ix < mesh_xy.template extent(); ++ix) { + for (std::size_t iy = 0; iy < mesh_xy.template extent(); ++iy) { + ddc::DiscreteElement const elem + = mesh_front + ddc::DiscreteVector(ix, iy); + spatial_moments_host(elem, ddc::DiscreteElement(0)) = 0.; + spatial_moments_host(elem, ddc::DiscreteElement(1)) = 0.; + } + } + } + { + ddc::DiscreteElement const x_face_front = x_face_dom.front(); + for (std::size_t ix = 0; ix < x_face_dom.template extent(); ++ix) { + for (std::size_t iy = 0; iy < x_face_dom.template extent(); ++iy) { + ddc::DiscreteElement const elem + = x_face_front + ddc::DiscreteVector(ix, iy); + ddc::DiscreteElement const left_node( + ddc::DiscreteElement(ddc::uid(elem)), + ddc::DiscreteElement(elem)); + double const value + = spatial_moment_x_host(elem, ddc::DiscreteElement(0)); + spatial_moments_host(left_node, ddc::DiscreteElement(0)) + += value / 2.; + spatial_moments_host( + left_node + ddc::DiscreteVector(1, 0), + ddc::DiscreteElement(0)) + += value / 2.; + } + } + } + { + ddc::DiscreteElement const y_face_front = y_face_dom.front(); + for (std::size_t ix = 0; ix < y_face_dom.template extent(); ++ix) { + for (std::size_t iy = 0; iy < y_face_dom.template extent(); ++iy) { + ddc::DiscreteElement const elem + = y_face_front + ddc::DiscreteVector(ix, iy); + ddc::DiscreteElement const lower_node( + ddc::DiscreteElement(elem), + ddc::DiscreteElement(ddc::uid(elem))); + double const value + = spatial_moment_y_host(elem, ddc::DiscreteElement(0)); + spatial_moments_host(lower_node, ddc::DiscreteElement(1)) + += value / 2.; + spatial_moments_host( + lower_node + ddc::DiscreteVector(0, 1), + ddc::DiscreteElement(1)) + += value / 2.; + } + } + } + ddc::parallel_deepcopy(spatial_moments, spatial_moments_host); ddc::parallel_for_each( Kokkos::DefaultExecutionSpace(), diff --git a/src/similie/exterior/exterior.hpp b/src/similie/exterior/exterior.hpp index 6ebf206f..e918cd55 100644 --- a/src/similie/exterior/exterior.hpp +++ b/src/similie/exterior/exterior.hpp @@ -14,3 +14,4 @@ #include "laplacian.hpp" #include "local_chain.hpp" #include "simplex.hpp" +#include "staggered_difference.hpp" diff --git a/src/similie/exterior/staggered_difference.hpp b/src/similie/exterior/staggered_difference.hpp new file mode 100644 index 00000000..62828278 --- /dev/null +++ b/src/similie/exterior/staggered_difference.hpp @@ -0,0 +1,180 @@ +// SPDX-FileCopyrightText: 2026 Baptiste Legouix +// SPDX-License-Identifier: MIT + +#pragma once + +#include + +#include + +#include +#include +#include + +namespace sil { + +namespace exterior { + +namespace detail { + +template +struct DiscreteDimensionForTypeSeq; + +template +struct DiscreteDimensionForTypeSeq> +{ + using type = std::conditional_t< + std::is_same_v, + DDim, + typename DiscreteDimensionForTypeSeq>::type>; +}; + +template +struct DiscreteDimensionForTypeSeq> +{ + using type = void; +}; + +template +using discrete_dimension_for_t + = typename DiscreteDimensionForTypeSeq>::type; + +template +KOKKOS_FUNCTION ddc::DiscreteElement remap_dimension(SourceElem const& source_elem) +{ + using source_d_dim_t + = discrete_dimension_for_t; + static_assert(!std::is_void_v); + return ddc::DiscreteElement(ddc::uid(source_elem)); +} + +template +struct RemapElement; + +template +struct RemapElement, SourceElem> +{ + KOKKOS_FUNCTION static ddc::DiscreteElement apply(SourceElem const& source_elem) + { + return ddc::DiscreteElement( + remap_dimension(source_elem)...); + } +}; + +template +KOKKOS_FUNCTION TargetElem remap_element(SourceElem const& source_elem) +{ + return RemapElement::apply(source_elem); +} + +template +KOKKOS_FUNCTION ddc::DiscreteElement component_element() +{ + using natural_index_t = tensor::uncharacterize_t; + return ddc::DiscreteElement( + ddc::type_seq_rank_v); +} + +} // namespace detail + +template < + class CDim, + misc::Specialization OutTensorType, + misc::Specialization InTensorType, + class ExecSpace> +OutTensorType staggered_deriv_component( + ExecSpace const& exec_space, + OutTensorType out, + InTensorType in) +{ + using out_index_t = ddc::type_seq_element_t< + 0, + ddc::to_type_seq_t>; + using in_index_t = ddc::type_seq_element_t< + 0, + ddc::to_type_seq_t>; + using in_d_dim_t = detail::discrete_dimension_for_t; + using in_elem_t = typename InTensorType::non_indices_domain_t::discrete_element_type; + using out_elem_t = typename OutTensorType::non_indices_domain_t::discrete_element_type; + static_assert(out_index_t::rank() == 0); + static_assert(in_index_t::rank() == 0); + static_assert(!std::is_void_v); + + ddc::parallel_for_each( + "similie_staggered_deriv_component", + exec_space, + out.non_indices_domain(), + KOKKOS_LAMBDA(out_elem_t out_elem) { + in_elem_t const in_left = detail::remap_element(out_elem); + in_elem_t const in_right = in_left + ddc::DiscreteVector(1); + double const dx = static_cast( + ddc::coordinate(ddc::DiscreteElement(in_right)) + - ddc::coordinate(ddc::DiscreteElement(in_left))); + out.mem(out_elem, ddc::DiscreteElement(0)) + = (in.get(in_right, ddc::DiscreteElement(0)) + - in.get(in_left, ddc::DiscreteElement(0))) + / dx; + }); + return out; +} + +template < + class AlphaCDim, + class BetaCDim, + misc::Specialization OutTensorType, + misc::Specialization FluxTensorType, + misc::Specialization MetricType, + class ExecSpace> +OutTensorType add_staggered_codifferential_component( + ExecSpace const& exec_space, + OutTensorType out, + FluxTensorType flux, + MetricType inv_metric) +{ + using out_index_t = ddc::type_seq_element_t< + 0, + ddc::to_type_seq_t>; + using out_d_dim_t + = detail::discrete_dimension_for_t; + using flux_d_dim_t = detail:: + discrete_dimension_for_t; + using out_elem_t = typename OutTensorType::non_indices_domain_t::discrete_element_type; + using flux_elem_t = typename FluxTensorType::non_indices_domain_t::discrete_element_type; + using flux_index_t = ddc::type_seq_element_t< + 0, + ddc::to_type_seq_t>; + static_assert(out_index_t::rank() == 0); + static_assert(!std::is_void_v); + static_assert(!std::is_void_v); + + ddc::parallel_for_each( + "similie_add_staggered_codifferential_component", + exec_space, + out.non_indices_domain(), + KOKKOS_LAMBDA(out_elem_t out_elem) { + ddc::DiscreteDomain dim_dom(out.non_indices_domain()); + ddc::DiscreteElement dim_elem(out_elem); + if (dim_elem.uid() == dim_dom.front().uid() || dim_elem.uid() == dim_dom.back().uid()) { + return; + } + + flux_elem_t const right_face = detail::remap_element(out_elem); + flux_elem_t const left_face = right_face - ddc::DiscreteVector(1); + double const dx = static_cast( + ddc::coordinate(ddc::DiscreteElement(right_face)) + - ddc::coordinate(ddc::DiscreteElement(left_face))); + double const deriv = (flux.get(right_face, ddc::DiscreteElement(0)) + - flux.get(left_face, ddc::DiscreteElement(0))) + / dx; + out.mem(out_elem, ddc::DiscreteElement(0)) + += inv_metric.get( + out_elem, + inv_metric.accessor().template access_element()) + * deriv; + }); + return out; +} + +} // namespace exterior + +} // namespace sil From 4e8769955eb0b840981ef6c68c06bb5d4edbe81d Mon Sep 17 00:00:00 2001 From: blegouix Date: Sat, 28 Mar 2026 23:58:06 +0100 Subject: [PATCH 02/14] wip --- .../2d_free_scalar_field.cpp | 75 +++++++------------ src/similie/similie.hpp | 1 + 2 files changed, 26 insertions(+), 50 deletions(-) diff --git a/examples/free_scalar_field/2d_free_scalar_field.cpp b/examples/free_scalar_field/2d_free_scalar_field.cpp index 96535f10..582c7cd6 100644 --- a/examples/free_scalar_field/2d_free_scalar_field.cpp +++ b/examples/free_scalar_field/2d_free_scalar_field.cpp @@ -238,14 +238,6 @@ struct DDimY : ddc::UniformPointSampling { }; -struct DDimXHalf : ddc::UniformPointSampling -{ -}; - -struct DDimYHalf : ddc::UniformPointSampling -{ -}; - /* struct BSplinesX : MesherXY::template bsplines_type { @@ -281,6 +273,11 @@ using BetaLow = sil::tensor::Covariant; using DummyIndex = sil::tensor::Covariant; +using XDualizer = sil::mesher::HalfShiftDualizer; +using YDualizer = sil::mesher::HalfShiftDualizer; +using DDimXDual = sil::mesher::dual_discrete_dimension_t; +using DDimYDual = sil::mesher::dual_discrete_dimension_t; + int main(int argc, char** argv) { // Initialize PDI, Kokkos and DDC @@ -313,29 +310,11 @@ int main(int argc, char** argv) ddc::select(lower_bounds), ddc::select(upper_bounds), ddc::select(nb_cells))); - auto const x_face_x_dom = ddc::init_discrete_space(DDimXHalf::init( - ddc::Coordinate( - ddc::select(lower_bounds) - + (ddc::select(upper_bounds) - ddc::select(lower_bounds)) - / (2 * (ddc::select(nb_cells).value() - 1))), - ddc::Coordinate( - ddc::select(upper_bounds) - - (ddc::select(upper_bounds) - ddc::select(lower_bounds)) - / (2 * (ddc::select(nb_cells).value() - 1))), - ddc::DiscreteVector(ddc::select(nb_cells).value() - 1))); - auto const y_face_y_dom = ddc::init_discrete_space(DDimYHalf::init( - ddc::Coordinate( - ddc::select(lower_bounds) - + (ddc::select(upper_bounds) - ddc::select(lower_bounds)) - / (2 * (ddc::select(nb_cells).value() - 1))), - ddc::Coordinate( - ddc::select(upper_bounds) - - (ddc::select(upper_bounds) - ddc::select(lower_bounds)) - / (2 * (ddc::select(nb_cells).value() - 1))), - ddc::DiscreteVector(ddc::select(nb_cells).value() - 1))); ddc::DiscreteDomain mesh_xy(x_dom, y_dom); - ddc::DiscreteDomain x_face_dom(x_face_x_dom, y_dom); - ddc::DiscreteDomain y_face_dom(x_dom, y_face_y_dom); + XDualizer const x_dualizer; + YDualizer const y_dualizer; + ddc::DiscreteDomain x_face_dom = x_dualizer(mesh_xy); + ddc::DiscreteDomain y_face_dom = y_dualizer(mesh_xy); assert(static_cast(mesh_xy.template extent()) == static_cast(mesh_xy.template extent())); @@ -427,23 +406,23 @@ int main(int argc, char** argv) // Staggered potential gradients [[maybe_unused]] sil::tensor::TensorAccessor scalar_accessor; - ddc::DiscreteDomain + ddc::DiscreteDomain potential_grad_x_dom(x_face_dom, scalar_accessor.domain()); ddc::Chunk potential_grad_x_alloc(potential_grad_x_dom, ddc::DeviceAllocator()); sil::tensor::Tensor potential_grad_x(potential_grad_x_alloc); - ddc::DiscreteDomain + ddc::DiscreteDomain potential_grad_y_dom(y_face_dom, scalar_accessor.domain()); ddc::Chunk potential_grad_y_alloc(potential_grad_y_dom, ddc::DeviceAllocator()); sil::tensor::Tensor potential_grad_y(potential_grad_y_alloc); // Staggered spatial moments - ddc::DiscreteDomain + ddc::DiscreteDomain spatial_moment_x_dom(x_face_dom, scalar_accessor.domain()); ddc::Chunk spatial_moment_x_alloc(spatial_moment_x_dom, ddc::DeviceAllocator()); sil::tensor::Tensor spatial_moment_x(spatial_moment_x_alloc); - ddc::DiscreteDomain + ddc::DiscreteDomain spatial_moment_y_dom(y_face_dom, scalar_accessor.domain()); ddc::Chunk spatial_moment_y_alloc(spatial_moment_y_dom, ddc::DeviceAllocator()); sil::tensor::Tensor spatial_moment_y(spatial_moment_y_alloc); @@ -562,7 +541,7 @@ int main(int argc, char** argv) ddc::parallel_for_each( Kokkos::DefaultExecutionSpace(), x_face_dom, - KOKKOS_LAMBDA(ddc::DiscreteElement elem) { + KOKKOS_LAMBDA(ddc::DiscreteElement elem) { spatial_moment_x(elem, ddc::DiscreteElement(0)) = FreeScalarFieldHamiltonian(mass).pi1( potential_grad_x(elem, ddc::DiscreteElement(0))); @@ -570,7 +549,7 @@ int main(int argc, char** argv) ddc::parallel_for_each( Kokkos::DefaultExecutionSpace(), y_face_dom, - KOKKOS_LAMBDA(ddc::DiscreteElement elem) { + KOKKOS_LAMBDA(ddc::DiscreteElement elem) { spatial_moment_y(elem, ddc::DiscreteElement(0)) = FreeScalarFieldHamiltonian(mass).pi2( potential_grad_y(elem, ddc::DiscreteElement(0))); @@ -642,14 +621,12 @@ int main(int argc, char** argv) } } { - ddc::DiscreteElement const x_face_front = x_face_dom.front(); - for (std::size_t ix = 0; ix < x_face_dom.template extent(); ++ix) { + ddc::DiscreteElement const x_face_front = x_face_dom.front(); + for (std::size_t ix = 0; ix < x_face_dom.template extent(); ++ix) { for (std::size_t iy = 0; iy < x_face_dom.template extent(); ++iy) { - ddc::DiscreteElement const elem - = x_face_front + ddc::DiscreteVector(ix, iy); - ddc::DiscreteElement const left_node( - ddc::DiscreteElement(ddc::uid(elem)), - ddc::DiscreteElement(elem)); + ddc::DiscreteElement const elem + = x_face_front + ddc::DiscreteVector(ix, iy); + ddc::DiscreteElement const left_node = x_dualizer.primal(elem); double const value = spatial_moment_x_host(elem, ddc::DiscreteElement(0)); spatial_moments_host(left_node, ddc::DiscreteElement(0)) @@ -662,14 +639,12 @@ int main(int argc, char** argv) } } { - ddc::DiscreteElement const y_face_front = y_face_dom.front(); + ddc::DiscreteElement const y_face_front = y_face_dom.front(); for (std::size_t ix = 0; ix < y_face_dom.template extent(); ++ix) { - for (std::size_t iy = 0; iy < y_face_dom.template extent(); ++iy) { - ddc::DiscreteElement const elem - = y_face_front + ddc::DiscreteVector(ix, iy); - ddc::DiscreteElement const lower_node( - ddc::DiscreteElement(elem), - ddc::DiscreteElement(ddc::uid(elem))); + for (std::size_t iy = 0; iy < y_face_dom.template extent(); ++iy) { + ddc::DiscreteElement const elem + = y_face_front + ddc::DiscreteVector(ix, iy); + ddc::DiscreteElement const lower_node = y_dualizer.primal(elem); double const value = spatial_moment_y_host(elem, ddc::DiscreteElement(0)); spatial_moments_host(lower_node, ddc::DiscreteElement(1)) diff --git a/src/similie/similie.hpp b/src/similie/similie.hpp index 8af73fbb..69c3178e 100644 --- a/src/similie/similie.hpp +++ b/src/similie/similie.hpp @@ -10,6 +10,7 @@ namespace sil { #include "csr/csr.hpp" #include "exterior/exterior.hpp" +#include "mesher/dualizer.hpp" #include "mesher/mesher.hpp" #include "tensor/tensor.hpp" #include "young_tableau/young_tableau.hpp" From 04675f8669d5852403820ffff242e9eaa576261d Mon Sep 17 00:00:00 2001 From: blegouix Date: Sun, 29 Mar 2026 00:12:47 +0100 Subject: [PATCH 03/14] wip --- .../2d_free_scalar_field.cpp | 30 ++- src/similie/exterior/coboundary.hpp | 46 +++++ src/similie/exterior/codifferential.hpp | 63 ++++++ src/similie/exterior/exterior.hpp | 1 - src/similie/exterior/staggered_difference.hpp | 180 ------------------ 5 files changed, 131 insertions(+), 189 deletions(-) delete mode 100644 src/similie/exterior/staggered_difference.hpp diff --git a/examples/free_scalar_field/2d_free_scalar_field.cpp b/examples/free_scalar_field/2d_free_scalar_field.cpp index 582c7cd6..37600eab 100644 --- a/examples/free_scalar_field/2d_free_scalar_field.cpp +++ b/examples/free_scalar_field/2d_free_scalar_field.cpp @@ -528,14 +528,20 @@ int main(int argc, char** argv) }); // Compute the staggered potential gradients - sil::exterior::staggered_deriv_component( + sil::exterior::deriv< + X, + DummyIndex>( Kokkos::DefaultExecutionSpace(), potential_grad_x, - half_step_potential); - sil::exterior::staggered_deriv_component( + half_step_potential, + x_dualizer); + sil::exterior::deriv< + Y, + DummyIndex>( Kokkos::DefaultExecutionSpace(), potential_grad_y, - half_step_potential); + half_step_potential, + y_dualizer); // Compute the spatial moments pi_\alpha by solving dphi/dx^\alpha = dH/dpi_\alpha ddc::parallel_for_each( @@ -562,16 +568,24 @@ int main(int argc, char** argv) KOKKOS_LAMBDA(ddc::DiscreteElement elem) { spatial_moments_div(elem) = 0.; }); - sil::exterior::add_staggered_codifferential_component( + sil::exterior::codifferential< + MetricIndex, + X, + DummyIndex>( Kokkos::DefaultExecutionSpace(), spatial_moments_div, spatial_moment_x, - inv_metric); - sil::exterior::add_staggered_codifferential_component( + inv_metric, + x_dualizer); + sil::exterior::codifferential< + MetricIndex, + Y, + DummyIndex>( Kokkos::DefaultExecutionSpace(), spatial_moments_div, spatial_moment_y, - inv_metric); + inv_metric, + y_dualizer); // Compute dpi_0/dx^0 by solving - dpi_0/dx^0 + dpi_\alpha/dx^\alpha = -dH/dphi and advect pi_0 by a time step dx^0. Also Then, perform the second phi half-advection by solving dphi/dx^0 = -dH/dpi_0 double const dS = (ddc::get(upper_bounds) - ddc::get(lower_bounds)) diff --git a/src/similie/exterior/coboundary.hpp b/src/similie/exterior/coboundary.hpp index b817ec14..885ab211 100644 --- a/src/similie/exterior/coboundary.hpp +++ b/src/similie/exterior/coboundary.hpp @@ -12,6 +12,7 @@ #include #include #include +#include #include #include #include @@ -334,6 +335,51 @@ coboundary_tensor_t deriv( TensorType>(exec_space, coboundary_tensor, tensor); } +template < + class TagToAddToCochain, + tensor::TensorIndex CochainTag, + misc::Specialization OutTensorType, + misc::Specialization TensorType, + class Dualizer, + class ExecSpace> +OutTensorType deriv( + ExecSpace const& exec_space, + OutTensorType out_tensor, + TensorType tensor, + Dualizer const& dualizer) +{ + using out_index_t = ddc::type_seq_element_t< + 0, + ddc::to_type_seq_t>; + using in_index_t = ddc::type_seq_element_t< + 0, + ddc::to_type_seq_t>; + using in_d_dim_t = mesher::detail:: + discrete_dimension_for_t; + using out_elem_t = typename OutTensorType::non_indices_domain_t::discrete_element_type; + using in_elem_t = typename TensorType::non_indices_domain_t::discrete_element_type; + static_assert(out_index_t::rank() == 0); + static_assert(in_index_t::rank() == 0); + static_assert(!std::is_void_v); + + ddc::parallel_for_each( + "similie_compute_centered_dualized_deriv", + exec_space, + out_tensor.non_indices_domain(), + KOKKOS_LAMBDA(out_elem_t out_elem) { + in_elem_t const in_left = dualizer.primal(out_elem); + in_elem_t const in_right = in_left + ddc::DiscreteVector(1); + double const dx = static_cast( + ddc::coordinate(ddc::DiscreteElement(in_right)) + - ddc::coordinate(ddc::DiscreteElement(in_left))); + out_tensor.mem(out_elem, ddc::DiscreteElement(0)) + = (tensor.get(in_right, ddc::DiscreteElement(0)) + - tensor.get(in_left, ddc::DiscreteElement(0))) + / dx; + }); + return out_tensor; +} + } // namespace exterior } // namespace sil diff --git a/src/similie/exterior/codifferential.hpp b/src/similie/exterior/codifferential.hpp index 49a0a2ea..0826f0e3 100644 --- a/src/similie/exterior/codifferential.hpp +++ b/src/similie/exterior/codifferential.hpp @@ -6,6 +6,7 @@ #include #include +#include #include #include #include @@ -279,6 +280,68 @@ codifferential_tensor_t codiffer return codifferential_tensor; } +template < + tensor::TensorIndex MetricIndex, + class TagToRemoveFromCochain, + tensor::TensorIndex CochainTag, + misc::Specialization OutTensorType, + misc::Specialization TensorType, + misc::Specialization MetricType, + class Dualizer, + class ExecSpace> +OutTensorType codifferential( + ExecSpace const& exec_space, + OutTensorType out_tensor, + TensorType tensor, + MetricType inv_metric, + Dualizer const& dualizer) +{ + using out_index_t = ddc::type_seq_element_t< + 0, + ddc::to_type_seq_t>; + using out_d_dim_t = mesher::detail:: + discrete_dimension_for_t; + using flux_d_dim_t = mesher::detail:: + discrete_dimension_for_t; + using out_elem_t = typename OutTensorType::non_indices_domain_t::discrete_element_type; + using flux_elem_t = typename TensorType::non_indices_domain_t::discrete_element_type; + using flux_index_t = ddc::type_seq_element_t< + 0, + ddc::to_type_seq_t>; + static_assert(out_index_t::rank() == 0); + static_assert(!std::is_void_v); + static_assert(!std::is_void_v); + + ddc::parallel_for_each( + "similie_compute_centered_dualized_codifferential", + exec_space, + out_tensor.non_indices_domain(), + KOKKOS_LAMBDA(out_elem_t out_elem) { + ddc::DiscreteDomain dim_dom(out_tensor.non_indices_domain()); + ddc::DiscreteElement dim_elem(out_elem); + if (dim_elem.uid() == dim_dom.front().uid() || dim_elem.uid() == dim_dom.back().uid()) { + return; + } + + flux_elem_t const right_face = dualizer(out_elem); + flux_elem_t const left_face = right_face - ddc::DiscreteVector(1); + double const dx = static_cast( + ddc::coordinate(ddc::DiscreteElement(right_face)) + - ddc::coordinate(ddc::DiscreteElement(left_face))); + double const deriv = (tensor.get(right_face, ddc::DiscreteElement(0)) + - tensor.get(left_face, ddc::DiscreteElement(0))) + / dx; + out_tensor.mem(out_elem, ddc::DiscreteElement(0)) + += inv_metric.get( + out_elem, + inv_metric.accessor().template access_element< + TagToRemoveFromCochain, + TagToRemoveFromCochain>()) + * deriv; + }); + return out_tensor; +} + } // namespace exterior } // namespace sil diff --git a/src/similie/exterior/exterior.hpp b/src/similie/exterior/exterior.hpp index e918cd55..6ebf206f 100644 --- a/src/similie/exterior/exterior.hpp +++ b/src/similie/exterior/exterior.hpp @@ -14,4 +14,3 @@ #include "laplacian.hpp" #include "local_chain.hpp" #include "simplex.hpp" -#include "staggered_difference.hpp" diff --git a/src/similie/exterior/staggered_difference.hpp b/src/similie/exterior/staggered_difference.hpp deleted file mode 100644 index 62828278..00000000 --- a/src/similie/exterior/staggered_difference.hpp +++ /dev/null @@ -1,180 +0,0 @@ -// SPDX-FileCopyrightText: 2026 Baptiste Legouix -// SPDX-License-Identifier: MIT - -#pragma once - -#include - -#include - -#include -#include -#include - -namespace sil { - -namespace exterior { - -namespace detail { - -template -struct DiscreteDimensionForTypeSeq; - -template -struct DiscreteDimensionForTypeSeq> -{ - using type = std::conditional_t< - std::is_same_v, - DDim, - typename DiscreteDimensionForTypeSeq>::type>; -}; - -template -struct DiscreteDimensionForTypeSeq> -{ - using type = void; -}; - -template -using discrete_dimension_for_t - = typename DiscreteDimensionForTypeSeq>::type; - -template -KOKKOS_FUNCTION ddc::DiscreteElement remap_dimension(SourceElem const& source_elem) -{ - using source_d_dim_t - = discrete_dimension_for_t; - static_assert(!std::is_void_v); - return ddc::DiscreteElement(ddc::uid(source_elem)); -} - -template -struct RemapElement; - -template -struct RemapElement, SourceElem> -{ - KOKKOS_FUNCTION static ddc::DiscreteElement apply(SourceElem const& source_elem) - { - return ddc::DiscreteElement( - remap_dimension(source_elem)...); - } -}; - -template -KOKKOS_FUNCTION TargetElem remap_element(SourceElem const& source_elem) -{ - return RemapElement::apply(source_elem); -} - -template -KOKKOS_FUNCTION ddc::DiscreteElement component_element() -{ - using natural_index_t = tensor::uncharacterize_t; - return ddc::DiscreteElement( - ddc::type_seq_rank_v); -} - -} // namespace detail - -template < - class CDim, - misc::Specialization OutTensorType, - misc::Specialization InTensorType, - class ExecSpace> -OutTensorType staggered_deriv_component( - ExecSpace const& exec_space, - OutTensorType out, - InTensorType in) -{ - using out_index_t = ddc::type_seq_element_t< - 0, - ddc::to_type_seq_t>; - using in_index_t = ddc::type_seq_element_t< - 0, - ddc::to_type_seq_t>; - using in_d_dim_t = detail::discrete_dimension_for_t; - using in_elem_t = typename InTensorType::non_indices_domain_t::discrete_element_type; - using out_elem_t = typename OutTensorType::non_indices_domain_t::discrete_element_type; - static_assert(out_index_t::rank() == 0); - static_assert(in_index_t::rank() == 0); - static_assert(!std::is_void_v); - - ddc::parallel_for_each( - "similie_staggered_deriv_component", - exec_space, - out.non_indices_domain(), - KOKKOS_LAMBDA(out_elem_t out_elem) { - in_elem_t const in_left = detail::remap_element(out_elem); - in_elem_t const in_right = in_left + ddc::DiscreteVector(1); - double const dx = static_cast( - ddc::coordinate(ddc::DiscreteElement(in_right)) - - ddc::coordinate(ddc::DiscreteElement(in_left))); - out.mem(out_elem, ddc::DiscreteElement(0)) - = (in.get(in_right, ddc::DiscreteElement(0)) - - in.get(in_left, ddc::DiscreteElement(0))) - / dx; - }); - return out; -} - -template < - class AlphaCDim, - class BetaCDim, - misc::Specialization OutTensorType, - misc::Specialization FluxTensorType, - misc::Specialization MetricType, - class ExecSpace> -OutTensorType add_staggered_codifferential_component( - ExecSpace const& exec_space, - OutTensorType out, - FluxTensorType flux, - MetricType inv_metric) -{ - using out_index_t = ddc::type_seq_element_t< - 0, - ddc::to_type_seq_t>; - using out_d_dim_t - = detail::discrete_dimension_for_t; - using flux_d_dim_t = detail:: - discrete_dimension_for_t; - using out_elem_t = typename OutTensorType::non_indices_domain_t::discrete_element_type; - using flux_elem_t = typename FluxTensorType::non_indices_domain_t::discrete_element_type; - using flux_index_t = ddc::type_seq_element_t< - 0, - ddc::to_type_seq_t>; - static_assert(out_index_t::rank() == 0); - static_assert(!std::is_void_v); - static_assert(!std::is_void_v); - - ddc::parallel_for_each( - "similie_add_staggered_codifferential_component", - exec_space, - out.non_indices_domain(), - KOKKOS_LAMBDA(out_elem_t out_elem) { - ddc::DiscreteDomain dim_dom(out.non_indices_domain()); - ddc::DiscreteElement dim_elem(out_elem); - if (dim_elem.uid() == dim_dom.front().uid() || dim_elem.uid() == dim_dom.back().uid()) { - return; - } - - flux_elem_t const right_face = detail::remap_element(out_elem); - flux_elem_t const left_face = right_face - ddc::DiscreteVector(1); - double const dx = static_cast( - ddc::coordinate(ddc::DiscreteElement(right_face)) - - ddc::coordinate(ddc::DiscreteElement(left_face))); - double const deriv = (flux.get(right_face, ddc::DiscreteElement(0)) - - flux.get(left_face, ddc::DiscreteElement(0))) - / dx; - out.mem(out_elem, ddc::DiscreteElement(0)) - += inv_metric.get( - out_elem, - inv_metric.accessor().template access_element()) - * deriv; - }); - return out; -} - -} // namespace exterior - -} // namespace sil From 2d9f7cab2304ecfd16a4edff4c093de6c54d766a Mon Sep 17 00:00:00 2001 From: blegouix Date: Sun, 29 Mar 2026 00:35:45 +0100 Subject: [PATCH 04/14] wip --- src/similie/mesher/dualizer.hpp | 273 ++++++++++++++++++++++++++++++++ 1 file changed, 273 insertions(+) create mode 100644 src/similie/mesher/dualizer.hpp diff --git a/src/similie/mesher/dualizer.hpp b/src/similie/mesher/dualizer.hpp new file mode 100644 index 00000000..2d6c0019 --- /dev/null +++ b/src/similie/mesher/dualizer.hpp @@ -0,0 +1,273 @@ +// SPDX-FileCopyrightText: 2026 Baptiste Legouix +// SPDX-License-Identifier: MIT + +#pragma once + +#include + +#include + +namespace sil { + +namespace mesher { + +struct HalfShifted +{ +}; + +template +struct Dual; + +template +struct HalfShiftDualizer; + +template +struct Dual + : ddc::UniformPointSampling +{ + using primal_discrete_dimension_type = PrimalDDim; + using dual_policy = HalfShifted; +}; + +namespace detail { + +template +struct DiscreteDimensionForTypeSeq; + +template +struct DiscreteDimensionForTypeSeq> +{ + using type = std::conditional_t< + std::is_same_v, + DDim, + typename DiscreteDimensionForTypeSeq>::type>; +}; + +template +struct DiscreteDimensionForTypeSeq> +{ + using type = void; +}; + +template +using discrete_dimension_for_t + = typename DiscreteDimensionForTypeSeq>::type; + +template +struct DualDiscreteDimension; + +template +struct DualDiscreteDimension, DDim> +{ + using type = std::conditional_t< + std::is_same_v, + Dual, + DDim>; +}; + +template +struct PrimalDiscreteDimension +{ + using type = DDim; +}; + +template +struct PrimalDiscreteDimension> +{ + using type = typename DDim::primal_discrete_dimension_type; +}; + +template +using primal_discrete_dimension_t = typename PrimalDiscreteDimension::type; + +template +struct Dualized; + +template +struct Dualized> +{ + using type = ddc::DiscreteElement::type...>; +}; + +template +struct Dualized> +{ + using type = ddc::DiscreteVector::type...>; +}; + +template +struct Dualized> +{ + using type = ddc::DiscreteDomain::type...>; +}; + +template +using dualized_t = typename Dualized::type; + +template +struct Primalized; + +template +struct Primalized> +{ + using type = ddc::DiscreteElement...>; +}; + +template +struct Primalized> +{ + using type = ddc::DiscreteVector...>; +}; + +template +struct Primalized> +{ + using type = ddc::DiscreteDomain...>; +}; + +template +using primalized_t = typename Primalized::type; + +template +KOKKOS_FUNCTION ddc::DiscreteElement dualizer_remap_dimension(SourceElem const& source_elem) +{ + using source_d_dim_t + = discrete_dimension_for_t; + static_assert(!std::is_void_v); + return ddc::DiscreteElement(ddc::uid(source_elem)); +} + +template +struct DualizerRemapElement; + +template +struct DualizerRemapElement, SourceElem> +{ + KOKKOS_FUNCTION static ddc::DiscreteElement apply(SourceElem const& source_elem) + { + return ddc::DiscreteElement( + dualizer_remap_dimension(source_elem)...); + } +}; + +template +KOKKOS_FUNCTION TargetElem dualizer_remap_element(SourceElem const& source_elem) +{ + return DualizerRemapElement::apply(source_elem); +} + +template +auto dualize_subdomain(ddc::DiscreteDomain const& subdomain) +{ + using dual_d_dim_t = typename DualDiscreteDimension::type; + if constexpr (std::is_same_v) { + return subdomain; + } else { + static_assert(ddc::is_uniform_point_sampling_v); + auto const nb_points = subdomain.template extent().value(); + assert(nb_points > 1); + auto const lower = ddc::coordinate(subdomain.front()) + + ddc::Coordinate( + ddc::step() / 2.); + auto const upper = ddc::coordinate(subdomain.back()) + - ddc::Coordinate( + ddc::step() / 2.); + return ddc::init_discrete_space(dual_d_dim_t::template init( + lower, + upper, + ddc::DiscreteVector(nb_points - 1))); + } +} + +template +struct DualizerApplication; + +template +struct DualizerApplication> +{ + using input_type = ddc::DiscreteElement; + using output_type = dualized_t; + + KOKKOS_FUNCTION static output_type apply(input_type const& elem) + { + return dualizer_remap_element(elem); + } +}; + +template +struct DualizerApplication> +{ + using input_type = ddc::DiscreteVector; + using output_type = dualized_t; + + KOKKOS_FUNCTION static output_type apply(input_type const& vect) + { + return output_type( + ddc::select(vect).value()...); + } +}; + +template +struct DualizerApplication> +{ + using input_type = ddc::DiscreteDomain; + using output_type = dualized_t; + + static output_type apply(input_type const& dom) + { + return output_type(dualize_subdomain(ddc::DiscreteDomain(dom))...); + } +}; + +template +struct PrimalizerApplication; + +template +struct PrimalizerApplication> +{ + using input_type = ddc::DiscreteElement; + using output_type = primalized_t; + + KOKKOS_FUNCTION static output_type apply(input_type const& elem) + { + return dualizer_remap_element(elem); + } +}; + +template +struct PrimalizerApplication> +{ + using input_type = ddc::DiscreteVector; + using output_type = primalized_t; + + KOKKOS_FUNCTION static output_type apply(input_type const& vect) + { + return output_type( + ddc::select(vect).value()...); + } +}; + +} // namespace detail + +template +using dual_discrete_dimension_t = typename detail::DualDiscreteDimension::type; + +template +struct HalfShiftDualizer +{ + template + constexpr detail::dualized_t, T> operator()(T const& obj) const + { + return detail::DualizerApplication, T>::apply(obj); + } + + template + constexpr detail::primalized_t primal(T const& obj) const + { + return detail::PrimalizerApplication::apply(obj); + } +}; + +} // namespace mesher + +} // namespace sil From eeb21d70a842c83649308c2d315f69fe16792c16 Mon Sep 17 00:00:00 2001 From: blegouix Date: Sun, 29 Mar 2026 00:49:37 +0100 Subject: [PATCH 05/14] wip --- .../2d_free_scalar_field.cpp | 42 +++------ src/similie/exterior/coboundary.hpp | 24 ++++++ src/similie/exterior/codifferential.hpp | 50 +++++++++-- src/similie/exterior/form.hpp | 85 +++++++++++++++++++ tests/exterior/exterior.cpp | 62 ++++++++++++++ 5 files changed, 226 insertions(+), 37 deletions(-) diff --git a/examples/free_scalar_field/2d_free_scalar_field.cpp b/examples/free_scalar_field/2d_free_scalar_field.cpp index 37600eab..194c42fe 100644 --- a/examples/free_scalar_field/2d_free_scalar_field.cpp +++ b/examples/free_scalar_field/2d_free_scalar_field.cpp @@ -415,6 +415,9 @@ int main(int argc, char** argv) potential_grad_y_dom(y_face_dom, scalar_accessor.domain()); ddc::Chunk potential_grad_y_alloc(potential_grad_y_dom, ddc::DeviceAllocator()); sil::tensor::Tensor potential_grad_y(potential_grad_y_alloc); + auto potential_grad = sil::exterior::make_tensor_form( + sil::exterior::component(potential_grad_x), + sil::exterior::component(potential_grad_y)); // Staggered spatial moments ddc::DiscreteDomain @@ -426,6 +429,9 @@ int main(int argc, char** argv) spatial_moment_y_dom(y_face_dom, scalar_accessor.domain()); ddc::Chunk spatial_moment_y_alloc(spatial_moment_y_dom, ddc::DeviceAllocator()); sil::tensor::Tensor spatial_moment_y(spatial_moment_y_alloc); + auto spatial_moment = sil::exterior::make_tensor_form( + sil::exterior::component(spatial_moment_x), + sil::exterior::component(spatial_moment_y)); auto spatial_moment_x_host = ddc:: create_mirror_view_and_copy(Kokkos::DefaultHostExecutionSpace(), spatial_moment_x); @@ -528,20 +534,7 @@ int main(int argc, char** argv) }); // Compute the staggered potential gradients - sil::exterior::deriv< - X, - DummyIndex>( - Kokkos::DefaultExecutionSpace(), - potential_grad_x, - half_step_potential, - x_dualizer); - sil::exterior::deriv< - Y, - DummyIndex>( - Kokkos::DefaultExecutionSpace(), - potential_grad_y, - half_step_potential, - y_dualizer); + sil::exterior::deriv(Kokkos::DefaultExecutionSpace(), potential_grad, half_step_potential); // Compute the spatial moments pi_\alpha by solving dphi/dx^\alpha = dH/dpi_\alpha ddc::parallel_for_each( @@ -569,23 +562,10 @@ int main(int argc, char** argv) spatial_moments_div(elem) = 0.; }); sil::exterior::codifferential< - MetricIndex, - X, - DummyIndex>( - Kokkos::DefaultExecutionSpace(), - spatial_moments_div, - spatial_moment_x, - inv_metric, - x_dualizer); - sil::exterior::codifferential< - MetricIndex, - Y, - DummyIndex>( - Kokkos::DefaultExecutionSpace(), - spatial_moments_div, - spatial_moment_y, - inv_metric, - y_dualizer); + MetricIndex>(Kokkos::DefaultExecutionSpace(), + spatial_moments_div, + spatial_moment, + inv_metric); // Compute dpi_0/dx^0 by solving - dpi_0/dx^0 + dpi_\alpha/dx^\alpha = -dH/dphi and advect pi_0 by a time step dx^0. Also Then, perform the second phi half-advection by solving dphi/dx^0 = -dH/dpi_0 double const dS = (ddc::get(upper_bounds) - ddc::get(lower_bounds)) diff --git a/src/similie/exterior/coboundary.hpp b/src/similie/exterior/coboundary.hpp index 885ab211..e0ee1011 100644 --- a/src/similie/exterior/coboundary.hpp +++ b/src/similie/exterior/coboundary.hpp @@ -21,6 +21,7 @@ #include "cochain.hpp" #include "cosimplex.hpp" +#include "form.hpp" namespace sil { @@ -380,6 +381,29 @@ OutTensorType deriv( return out_tensor; } +template < + class... Components, + misc::Specialization TensorType, + class ExecSpace> +TensorForm deriv( + ExecSpace const& exec_space, + TensorForm out_form, + TensorType tensor) +{ + using in_index_t = ddc::type_seq_element_t< + 0, + ddc::to_type_seq_t>; + static_assert(in_index_t::rank() == 0); + + (..., + sil::exterior::deriv( + exec_space, + out_form.template component(), + tensor, + mesher::HalfShiftDualizer {})); + return out_form; +} + } // namespace exterior } // namespace sil diff --git a/src/similie/exterior/codifferential.hpp b/src/similie/exterior/codifferential.hpp index 0826f0e3..bd971c58 100644 --- a/src/similie/exterior/codifferential.hpp +++ b/src/similie/exterior/codifferential.hpp @@ -16,6 +16,7 @@ #include "cochain.hpp" #include "cosimplex.hpp" +#include "form.hpp" namespace sil { @@ -299,6 +300,9 @@ OutTensorType codifferential( using out_index_t = ddc::type_seq_element_t< 0, ddc::to_type_seq_t>; + using metric_component_index_t = ddc::type_seq_element_t< + 0, + ddc::to_type_seq_t>; using out_d_dim_t = mesher::detail:: discrete_dimension_for_t; using flux_d_dim_t = mesher::detail:: @@ -331,17 +335,51 @@ OutTensorType codifferential( double const deriv = (tensor.get(right_face, ddc::DiscreteElement(0)) - tensor.get(left_face, ddc::DiscreteElement(0))) / dx; + double const metric_factor = [&]() -> double { + if constexpr ( + misc::Specialization) { + return 1.; + } else { + return inv_metric( + out_elem, + inv_metric.accessor().template access_element< + TagToRemoveFromCochain, + TagToRemoveFromCochain>()); + } + }(); out_tensor.mem(out_elem, ddc::DiscreteElement(0)) - += inv_metric.get( - out_elem, - inv_metric.accessor().template access_element< - TagToRemoveFromCochain, - TagToRemoveFromCochain>()) - * deriv; + += metric_factor * deriv; }); return out_tensor; } +template < + tensor::TensorIndex MetricIndex, + misc::Specialization OutTensorType, + class... Components, + misc::Specialization MetricType, + class ExecSpace> +OutTensorType codifferential( + ExecSpace const& exec_space, + OutTensorType out_tensor, + TensorForm tensor_form, + MetricType inv_metric) +{ + (..., + sil::exterior::codifferential< + MetricIndex, + typename Components::tag, + ddc::type_seq_element_t< + 0, + ddc::to_type_seq_t>>( + exec_space, + out_tensor, + tensor_form.template component(), + inv_metric, + mesher::HalfShiftDualizer {})); + return out_tensor; +} + } // namespace exterior } // namespace sil diff --git a/src/similie/exterior/form.hpp b/src/similie/exterior/form.hpp index 99b2b384..625a0b39 100644 --- a/src/similie/exterior/form.hpp +++ b/src/similie/exterior/form.hpp @@ -3,8 +3,14 @@ #pragma once +#include +#include + #include +#include +#include + #include "cochain.hpp" #include "cosimplex.hpp" @@ -43,6 +49,85 @@ FormWrapper(Head, ElementType, Allocator) -> FormWrapper using Form = typename detail::FormWrapper::type; +template TensorType> +struct FormComponent +{ + using tag = Tag; + using tensor_type = TensorType; + + TensorType tensor; + + KOKKOS_FUNCTION explicit constexpr FormComponent(TensorType tensor_) noexcept : tensor(tensor_) {} +}; + +template TensorType> +KOKKOS_FUNCTION constexpr FormComponent component(TensorType tensor) noexcept +{ + return FormComponent(tensor); +} + +namespace detail { + +template +struct FormComponentByTag; + +template +struct FormComponentByTag +{ + using type = std::conditional_t< + std::is_same_v, + HeadComponent, + typename FormComponentByTag::type>; +}; + +template +struct FormComponentByTag +{ + using type = void; +}; + +template +constexpr bool has_unique_component_tags_v + = ((... && !std::is_same_v) || true); + +template +constexpr bool has_unique_component_tags_v + = ((... && !std::is_same_v)) + && has_unique_component_tags_v; + +} // namespace detail + +template +class TensorForm +{ + static_assert(sizeof...(Components) > 0); + static_assert(detail::has_unique_component_tags_v); + + std::tuple m_components; + +public: + using components_tuple_type = std::tuple; + + constexpr explicit TensorForm(Components... components) noexcept + : m_components(std::move(components)...) + { + } + + template + KOKKOS_FUNCTION constexpr auto component() const noexcept + { + using component_t = typename detail::FormComponentByTag::type; + static_assert(!std::is_void_v); + return std::get(m_components).tensor; + } +}; + +template +constexpr TensorForm make_tensor_form(Components... components) noexcept +{ + return TensorForm(components...); +} + } // namespace exterior } // namespace sil diff --git a/tests/exterior/exterior.cpp b/tests/exterior/exterior.cpp index c19afd24..c52b1a37 100644 --- a/tests/exterior/exterior.cpp +++ b/tests/exterior/exterior.cpp @@ -514,6 +514,68 @@ TEST(Form, Alias) 0.); } +TEST(Form, TensorFormDeriv) +{ + using DummyIndex = sil::tensor::Covariant; + using XDualizer = sil::mesher::HalfShiftDualizer; + using YDualizer = sil::mesher::HalfShiftDualizer; + using DDimXDual = sil::mesher::dual_discrete_dimension_t; + using DDimYDual = sil::mesher::dual_discrete_dimension_t; + + auto const x_dom = ddc::init_discrete_space(DDimX::init( + ddc::Coordinate(0.), + ddc::Coordinate(1.), + ddc::DiscreteVector(6))); + auto const y_dom = ddc::init_discrete_space(DDimY::init( + ddc::Coordinate(0.), + ddc::Coordinate(1.), + ddc::DiscreteVector(6))); + ddc::DiscreteDomain const mesh(x_dom, y_dom); + XDualizer const x_dualizer; + YDualizer const y_dualizer; + ddc::DiscreteDomain const x_face_dom = x_dualizer(mesh); + ddc::DiscreteDomain const y_face_dom = y_dualizer(mesh); + + [[maybe_unused]] sil::tensor::TensorAccessor scalar_accessor; + ddc::Chunk scalar_alloc( + ddc::DiscreteDomain(mesh, scalar_accessor.domain()), + ddc::HostAllocator()); + sil::tensor::Tensor scalar(scalar_alloc); + ddc::parallel_for_each( + Kokkos::DefaultHostExecutionSpace(), + mesh, + KOKKOS_LAMBDA(ddc::DiscreteElement elem) { + double const x = ddc::coordinate(ddc::DiscreteElement(elem)); + double const y = ddc::coordinate(ddc::DiscreteElement(elem)); + scalar(elem, ddc::DiscreteElement(0)) = x * x + y; + }); + + ddc::Chunk grad_x_alloc( + ddc::DiscreteDomain( + x_face_dom, + scalar_accessor.domain()), + ddc::HostAllocator()); + sil::tensor::Tensor grad_x(grad_x_alloc); + ddc::Chunk grad_y_alloc( + ddc::DiscreteDomain( + y_face_dom, + scalar_accessor.domain()), + ddc::HostAllocator()); + sil::tensor::Tensor grad_y(grad_y_alloc); + auto grad_form = sil::exterior::make_tensor_form( + sil::exterior::component(grad_x), + sil::exterior::component(grad_y)); + + sil::exterior::deriv(Kokkos::DefaultHostExecutionSpace(), grad_form, scalar); + + ddc::DiscreteElement const x_face_elem + = x_face_dom.front() + ddc::DiscreteVector(2, 2); + ddc::DiscreteElement const y_face_elem + = y_face_dom.front() + ddc::DiscreteVector(2, 2); + EXPECT_NEAR(grad_x(x_face_elem, ddc::DiscreteElement(0)), 1.0, 1e-12); + EXPECT_NEAR(grad_y(y_face_elem, ddc::DiscreteElement(0)), 1.0, 1e-12); +} + TEST(Cochain, Test) { sil::exterior::Chain From 471f60fc24e83aad446105054b8ba9c04e36df75 Mon Sep 17 00:00:00 2001 From: blegouix Date: Sun, 29 Mar 2026 01:02:33 +0100 Subject: [PATCH 06/14] wip --- src/similie/exterior/coboundary.hpp | 5 +- src/similie/exterior/codifferential.hpp | 3 +- src/similie/exterior/form.hpp | 35 +++++++-- src/similie/exterior/hodge_star.hpp | 66 +++++++++++++++++ tests/exterior/hodge_star.cpp | 95 +++++++++++++++++++++++++ 5 files changed, 197 insertions(+), 7 deletions(-) diff --git a/src/similie/exterior/coboundary.hpp b/src/similie/exterior/coboundary.hpp index e0ee1011..9a05e0bc 100644 --- a/src/similie/exterior/coboundary.hpp +++ b/src/similie/exterior/coboundary.hpp @@ -382,12 +382,13 @@ OutTensorType deriv( } template < + class SupportTag, class... Components, misc::Specialization TensorType, class ExecSpace> -TensorForm deriv( +TensorForm deriv( ExecSpace const& exec_space, - TensorForm out_form, + TensorForm out_form, TensorType tensor) { using in_index_t = ddc::type_seq_element_t< diff --git a/src/similie/exterior/codifferential.hpp b/src/similie/exterior/codifferential.hpp index bd971c58..56524c8f 100644 --- a/src/similie/exterior/codifferential.hpp +++ b/src/similie/exterior/codifferential.hpp @@ -356,13 +356,14 @@ OutTensorType codifferential( template < tensor::TensorIndex MetricIndex, misc::Specialization OutTensorType, + class SupportTag, class... Components, misc::Specialization MetricType, class ExecSpace> OutTensorType codifferential( ExecSpace const& exec_space, OutTensorType out_tensor, - TensorForm tensor_form, + TensorForm tensor_form, MetricType inv_metric) { (..., diff --git a/src/similie/exterior/form.hpp b/src/similie/exterior/form.hpp index 625a0b39..c190c1e7 100644 --- a/src/similie/exterior/form.hpp +++ b/src/similie/exterior/form.hpp @@ -18,6 +18,32 @@ namespace sil { namespace exterior { +struct PrimalSupport +{ +}; + +struct DualSupport +{ +}; + +template +struct HodgeDualSupport; + +template <> +struct HodgeDualSupport +{ + using type = DualSupport; +}; + +template <> +struct HodgeDualSupport +{ + using type = PrimalSupport; +}; + +template +using hodge_dual_support_t = typename HodgeDualSupport::type; + namespace detail { template @@ -97,7 +123,7 @@ constexpr bool has_unique_component_tags_v } // namespace detail -template +template class TensorForm { static_assert(sizeof...(Components) > 0); @@ -106,6 +132,7 @@ class TensorForm std::tuple m_components; public: + using support_tag = SupportTag; using components_tuple_type = std::tuple; constexpr explicit TensorForm(Components... components) noexcept @@ -122,10 +149,10 @@ class TensorForm } }; -template -constexpr TensorForm make_tensor_form(Components... components) noexcept +template +constexpr TensorForm make_tensor_form(Components... components) noexcept { - return TensorForm(components...); + return TensorForm(components...); } } // namespace exterior diff --git a/src/similie/exterior/hodge_star.hpp b/src/similie/exterior/hodge_star.hpp index c3693e1c..f5071990 100644 --- a/src/similie/exterior/hodge_star.hpp +++ b/src/similie/exterior/hodge_star.hpp @@ -17,6 +17,8 @@ #include #include +#include "form.hpp" + namespace sil { namespace exterior { @@ -161,6 +163,70 @@ HodgeStarType fill_hodge_star( return fill_hodge_star(exec_space, hodge_star, metric_det, metric_prod); } +template < + class SupportTag, + class... OutComponents, + class FirstComponent, + class SecondComponent, + class MetricType, + class ExecSpace> +auto hodge_star( + ExecSpace const& exec_space, + TensorForm, OutComponents...> out_form, + TensorForm in_form, + MetricType const&) +{ + using first_in_tensor_t = typename FirstComponent::tensor_type; + using second_in_tensor_t = typename SecondComponent::tensor_type; + using out_first_component_t = typename detail:: + FormComponentByTag::type; + using out_second_component_t = typename detail:: + FormComponentByTag::type; + using first_in_index_t = ddc::type_seq_element_t< + 0, + ddc::to_type_seq_t>; + using second_in_index_t = ddc::type_seq_element_t< + 0, + ddc::to_type_seq_t>; + using first_out_index_t = ddc::type_seq_element_t< + 0, + ddc::to_type_seq_t>; + using second_out_index_t = ddc::type_seq_element_t< + 0, + ddc::to_type_seq_t>; + static_assert(!std::is_void_v); + static_assert(!std::is_void_v); + static_assert(std::is_same_v); + static_assert(std::is_same_v); + static_assert(first_in_index_t::rank() == 0); + static_assert(second_in_index_t::rank() == 0); + static_assert(first_out_index_t::rank() == 0); + static_assert(second_out_index_t::rank() == 0); + + auto const in_first = in_form.template component(); + auto const in_second = in_form.template component(); + auto const out_first = out_form.template component(); + auto const out_second = out_form.template component(); + + ddc::parallel_for_each( + "similie_compute_tensor_form_hodge_star_first_component", + exec_space, + out_first.non_indices_domain(), + KOKKOS_LAMBDA(typename second_in_tensor_t::non_indices_domain_t::discrete_element_type elem) { + out_first.mem(elem, ddc::DiscreteElement(0)) + = -in_second.get(elem, ddc::DiscreteElement(0)); + }); + ddc::parallel_for_each( + "similie_compute_tensor_form_hodge_star_second_component", + exec_space, + out_second.non_indices_domain(), + KOKKOS_LAMBDA(typename first_in_tensor_t::non_indices_domain_t::discrete_element_type elem) { + out_second.mem(elem, ddc::DiscreteElement(0)) + = in_first.get(elem, ddc::DiscreteElement(0)); + }); + return out_form; +} + } // namespace exterior } // namespace sil diff --git a/tests/exterior/hodge_star.cpp b/tests/exterior/hodge_star.cpp index fc8b76f6..5a67b21a 100644 --- a/tests/exterior/hodge_star.cpp +++ b/tests/exterior/hodge_star.cpp @@ -7,6 +7,7 @@ #include #include +#include #include #include #include @@ -63,6 +64,8 @@ using HodgeStarDomain = sil::exterior:: using HodgeStarDomain2 = sil::exterior:: hodge_star_domain_t, ddc::detail::TypeSeq>; +using DummyIndex = sil::tensor::Covariant; + TEST(HodgeStar, Test) { ddc::DiscreteDomain @@ -137,3 +140,95 @@ TEST(HodgeStar, Test) EXPECT_DOUBLE_EQ(form(elem, form.accessor().access_element()), 3.); }); } + +TEST(HodgeStar, TensorForm1In2D) +{ + using XDualizer = sil::mesher::HalfShiftDualizer; + using YDualizer = sil::mesher::HalfShiftDualizer; + using DDimXDual = sil::mesher::dual_discrete_dimension_t; + using DDimYDual = sil::mesher::dual_discrete_dimension_t; + + auto const x_dom = ddc::init_discrete_space(DDimX::init( + ddc::Coordinate(0.), + ddc::Coordinate(1.), + ddc::DiscreteVector(4))); + auto const y_dom = ddc::init_discrete_space(DDimY::init( + ddc::Coordinate(0.), + ddc::Coordinate(1.), + ddc::DiscreteVector(4))); + ddc::DiscreteDomain const mesh(x_dom, y_dom); + XDualizer const x_dualizer; + YDualizer const y_dualizer; + ddc::DiscreteDomain const x_face_dom = x_dualizer(mesh); + ddc::DiscreteDomain const y_face_dom = y_dualizer(mesh); + + [[maybe_unused]] sil::tensor::TensorAccessor scalar_accessor; + + ddc::Chunk primal_x_alloc( + ddc::DiscreteDomain(x_face_dom, scalar_accessor.domain()), + ddc::HostAllocator()); + sil::tensor::Tensor primal_x(primal_x_alloc); + ddc::Chunk primal_y_alloc( + ddc::DiscreteDomain(y_face_dom, scalar_accessor.domain()), + ddc::HostAllocator()); + sil::tensor::Tensor primal_y(primal_y_alloc); + + ddc::parallel_for_each( + Kokkos::DefaultHostExecutionSpace(), + primal_x.domain(), + KOKKOS_LAMBDA(ddc::DiscreteElement elem) { + primal_x(elem) = 2.; + }); + ddc::parallel_for_each( + Kokkos::DefaultHostExecutionSpace(), + primal_y.domain(), + KOKKOS_LAMBDA(ddc::DiscreteElement elem) { + primal_y(elem) = 3.; + }); + + auto primal_form = sil::exterior::make_tensor_form( + sil::exterior::component(primal_x), + sil::exterior::component(primal_y)); + + ddc::Chunk dual_x_alloc( + ddc::DiscreteDomain(y_face_dom, scalar_accessor.domain()), + ddc::HostAllocator()); + sil::tensor::Tensor dual_x(dual_x_alloc); + ddc::Chunk dual_y_alloc( + ddc::DiscreteDomain(x_face_dom, scalar_accessor.domain()), + ddc::HostAllocator()); + sil::tensor::Tensor dual_y(dual_y_alloc); + auto dual_form = sil::exterior::make_tensor_form( + sil::exterior::component(dual_x), + sil::exterior::component(dual_y)); + + sil::exterior::hodge_star(Kokkos::DefaultHostExecutionSpace(), dual_form, primal_form, 0); + + EXPECT_DOUBLE_EQ( + dual_x(y_face_dom.front(), ddc::DiscreteElement(0)), + -3.); + EXPECT_DOUBLE_EQ( + dual_y(x_face_dom.front(), ddc::DiscreteElement(0)), + 2.); + + ddc::Chunk primal_back_x_alloc( + ddc::DiscreteDomain(x_face_dom, scalar_accessor.domain()), + ddc::HostAllocator()); + sil::tensor::Tensor primal_back_x(primal_back_x_alloc); + ddc::Chunk primal_back_y_alloc( + ddc::DiscreteDomain(y_face_dom, scalar_accessor.domain()), + ddc::HostAllocator()); + sil::tensor::Tensor primal_back_y(primal_back_y_alloc); + auto primal_back = sil::exterior::make_tensor_form( + sil::exterior::component(primal_back_x), + sil::exterior::component(primal_back_y)); + + sil::exterior::hodge_star(Kokkos::DefaultHostExecutionSpace(), primal_back, dual_form, 0); + + EXPECT_DOUBLE_EQ( + primal_back_x(x_face_dom.front(), ddc::DiscreteElement(0)), + -2.); + EXPECT_DOUBLE_EQ( + primal_back_y(y_face_dom.front(), ddc::DiscreteElement(0)), + -3.); +} From 0c67e9de85bb6aecfab7db8e5fd8820c5bbe1691 Mon Sep 17 00:00:00 2001 From: blegouix Date: Sun, 29 Mar 2026 01:20:20 +0100 Subject: [PATCH 07/14] wip --- src/similie/exterior/codifferential.hpp | 135 +++++++++++++++++++++--- tests/exterior/exterior.cpp | 92 ++++++++++++++++ 2 files changed, 212 insertions(+), 15 deletions(-) diff --git a/src/similie/exterior/codifferential.hpp b/src/similie/exterior/codifferential.hpp index 56524c8f..83bbeba6 100644 --- a/src/similie/exterior/codifferential.hpp +++ b/src/similie/exterior/codifferential.hpp @@ -158,6 +158,107 @@ struct CodifferentialDummyIndexSeq using type = typename CodifferentialDummyIndexSeq_, T>::type; }; +template +constexpr bool is_dual_discrete_dimension_v + = !std::is_same_v>; + +template < + class AxisTag, + misc::Specialization OutTensorType, + misc::Specialization TensorType> +KOKKOS_FUNCTION double centered_form_derivative( + typename OutTensorType::non_indices_domain_t::discrete_element_type out_elem, + TensorType tensor) +{ + using out_d_dim_t = mesher::detail:: + discrete_dimension_for_t; + using comp_d_dim_t = mesher::detail:: + discrete_dimension_for_t; + using comp_elem_t = typename TensorType::non_indices_domain_t::discrete_element_type; + using comp_index_t = ddc::type_seq_element_t< + 0, + ddc::to_type_seq_t>; + static_assert(!std::is_void_v); + static_assert(!std::is_void_v); + static_assert(comp_index_t::rank() == 0); + + comp_elem_t const aligned_elem = mesher::detail::dualizer_remap_element(out_elem); + if constexpr (is_dual_discrete_dimension_v) { + auto const left = aligned_elem; + auto const right = left + ddc::DiscreteVector(1); + double const dx = static_cast( + ddc::coordinate(ddc::DiscreteElement(right)) + - ddc::coordinate(ddc::DiscreteElement(left))); + return (tensor.get(right, ddc::DiscreteElement(0)) + - tensor.get(left, ddc::DiscreteElement(0))) + / dx; + } else { + auto const right = aligned_elem; + auto const left = right - ddc::DiscreteVector(1); + double const dx = static_cast( + ddc::coordinate(ddc::DiscreteElement(right)) + - ddc::coordinate(ddc::DiscreteElement(left))); + return (tensor.get(right, ddc::DiscreteElement(0)) + - tensor.get(left, ddc::DiscreteElement(0))) + / dx; + } +} + +template < + misc::Specialization OutTensorType, + class SupportTag, + class FirstComponent, + class SecondComponent, + class ExecSpace> +OutTensorType exterior_derivative_of_tensor_form_2d( + ExecSpace const& exec_space, + OutTensorType out_tensor, + TensorForm tensor_form) +{ + using out_index_t = ddc::type_seq_element_t< + 0, + ddc::to_type_seq_t>; + using first_out_d_dim_t = mesher::detail:: + discrete_dimension_for_t; + using second_out_d_dim_t = mesher::detail:: + discrete_dimension_for_t; + static_assert(out_index_t::rank() == 0); + + auto const first_tensor = tensor_form.template component(); + auto const second_tensor = tensor_form.template component(); + + ddc::parallel_for_each( + "similie_compute_tensor_form_exterior_derivative_2d", + exec_space, + out_tensor.non_indices_domain(), + KOKKOS_LAMBDA(typename OutTensorType::non_indices_domain_t::discrete_element_type elem) { + if constexpr (!is_dual_discrete_dimension_v) { + ddc::DiscreteDomain dim_dom(out_tensor.non_indices_domain()); + ddc::DiscreteElement dim_elem(elem); + if (dim_elem.uid() == dim_dom.front().uid() + || dim_elem.uid() == dim_dom.back().uid()) { + return; + } + } + if constexpr (!is_dual_discrete_dimension_v) { + ddc::DiscreteDomain dim_dom(out_tensor.non_indices_domain()); + ddc::DiscreteElement dim_elem(elem); + if (dim_elem.uid() == dim_dom.front().uid() + || dim_elem.uid() == dim_dom.back().uid()) { + return; + } + } + out_tensor.mem(elem, ddc::DiscreteElement(0)) + += centered_form_derivative( + elem, + second_tensor) + - centered_form_derivative( + elem, + first_tensor); + }); + return out_tensor; +} + } // namespace detail template < @@ -357,28 +458,32 @@ template < tensor::TensorIndex MetricIndex, misc::Specialization OutTensorType, class SupportTag, - class... Components, + class FirstComponent, + class SecondComponent, misc::Specialization MetricType, class ExecSpace> OutTensorType codifferential( ExecSpace const& exec_space, OutTensorType out_tensor, - TensorForm tensor_form, + TensorForm tensor_form, MetricType inv_metric) { - (..., - sil::exterior::codifferential< - MetricIndex, - typename Components::tag, - ddc::type_seq_element_t< - 0, - ddc::to_type_seq_t>>( - exec_space, - out_tensor, - tensor_form.template component(), - inv_metric, - mesher::HalfShiftDualizer {})); - return out_tensor; + ddc::Chunk dual_first_alloc( + typename SecondComponent::tensor_type::discrete_domain_type( + tensor_form.template component().domain()), + ddc::KokkosAllocator()); + sil::tensor::Tensor dual_first(dual_first_alloc); + ddc::Chunk dual_second_alloc( + typename FirstComponent::tensor_type::discrete_domain_type( + tensor_form.template component().domain()), + ddc::KokkosAllocator()); + sil::tensor::Tensor dual_second(dual_second_alloc); + auto dual_form = make_tensor_form>( + component(dual_first), + component(dual_second)); + + sil::exterior::hodge_star(exec_space, dual_form, tensor_form, inv_metric); + return detail::exterior_derivative_of_tensor_form_2d(exec_space, out_tensor, dual_form); } } // namespace exterior diff --git a/tests/exterior/exterior.cpp b/tests/exterior/exterior.cpp index c52b1a37..bc27f123 100644 --- a/tests/exterior/exterior.cpp +++ b/tests/exterior/exterior.cpp @@ -576,6 +576,98 @@ TEST(Form, TensorFormDeriv) EXPECT_NEAR(grad_y(y_face_elem, ddc::DiscreteElement(0)), 1.0, 1e-12); } +TEST(Form, TensorFormCodifferential) +{ + struct DDimX2 : ddc::UniformPointSampling + { + }; + struct DDimY2 : ddc::UniformPointSampling + { + }; + using MetricIndex = sil::tensor::TensorIdentityIndex< + sil::tensor::Covariant>, + sil::tensor::Covariant>>; + using InverseMetricIndex = sil::tensor::upper_t; + using DummyIndex = sil::tensor::Covariant; + using XDualizer = sil::mesher::HalfShiftDualizer; + using YDualizer = sil::mesher::HalfShiftDualizer; + using DDimXDual = sil::mesher::dual_discrete_dimension_t; + using DDimYDual = sil::mesher::dual_discrete_dimension_t; + + auto const x_dom = ddc::init_discrete_space(DDimX2::init( + ddc::Coordinate(0.), + ddc::Coordinate(1.), + ddc::DiscreteVector(6))); + auto const y_dom = ddc::init_discrete_space(DDimY2::init( + ddc::Coordinate(0.), + ddc::Coordinate(1.), + ddc::DiscreteVector(6))); + ddc::DiscreteDomain const mesh(x_dom, y_dom); + XDualizer const x_dualizer; + YDualizer const y_dualizer; + ddc::DiscreteDomain const x_face_dom = x_dualizer(mesh); + ddc::DiscreteDomain const y_face_dom = y_dualizer(mesh); + + [[maybe_unused]] sil::tensor::TensorAccessor scalar_accessor; + ddc::Chunk grad_x_alloc( + ddc::DiscreteDomain( + x_face_dom, + scalar_accessor.domain()), + ddc::HostAllocator()); + sil::tensor::Tensor grad_x(grad_x_alloc); + ddc::Chunk grad_y_alloc( + ddc::DiscreteDomain( + y_face_dom, + scalar_accessor.domain()), + ddc::HostAllocator()); + sil::tensor::Tensor grad_y(grad_y_alloc); + + ddc::parallel_for_each( + Kokkos::DefaultHostExecutionSpace(), + grad_x.domain(), + KOKKOS_LAMBDA(ddc::DiscreteElement elem) { + double const x = ddc::coordinate(ddc::DiscreteElement(elem)); + grad_x(elem) = 2. * x; + }); + ddc::parallel_for_each( + Kokkos::DefaultHostExecutionSpace(), + grad_y.domain(), + KOKKOS_LAMBDA(ddc::DiscreteElement elem) { + grad_y(elem) = 1.; + }); + + auto form = sil::exterior::make_tensor_form( + sil::exterior::component(grad_x), + sil::exterior::component(grad_y)); + + [[maybe_unused]] sil::tensor::TensorAccessor inv_metric_accessor; + ddc::Chunk inv_metric_alloc( + ddc::DiscreteDomain( + mesh, + inv_metric_accessor.domain()), + ddc::HostAllocator()); + sil::tensor::Tensor inv_metric(inv_metric_alloc); + + ddc::Chunk div_alloc( + ddc::DiscreteDomain(mesh, scalar_accessor.domain()), + ddc::HostAllocator()); + sil::tensor::Tensor div(div_alloc); + ddc::parallel_for_each( + Kokkos::DefaultHostExecutionSpace(), + div.domain(), + KOKKOS_LAMBDA(ddc::DiscreteElement elem) { div(elem) = 0.; }); + + sil::exterior::codifferential( + Kokkos::DefaultHostExecutionSpace(), + div, + form, + inv_metric); + + ddc::DiscreteElement const center + = mesh.front() + ddc::DiscreteVector(2, 2); + EXPECT_NEAR(div(center, ddc::DiscreteElement(0)), 2.0, 1e-12); +} + TEST(Cochain, Test) { sil::exterior::Chain From c51b4ef735526aa357b57a6848e5e410f39374ed Mon Sep 17 00:00:00 2001 From: blegouix Date: Sun, 29 Mar 2026 01:30:19 +0100 Subject: [PATCH 08/14] wip --- src/similie/exterior/hodge_star.hpp | 106 ++++++++++++++++++++++++++++ tests/exterior/hodge_star.cpp | 91 ++++++++++++++++++++++++ 2 files changed, 197 insertions(+) diff --git a/src/similie/exterior/hodge_star.hpp b/src/similie/exterior/hodge_star.hpp index f5071990..e84f5574 100644 --- a/src/similie/exterior/hodge_star.hpp +++ b/src/similie/exterior/hodge_star.hpp @@ -7,6 +7,7 @@ #include #include +#include #include #include #include @@ -23,6 +24,62 @@ namespace sil { namespace exterior { +namespace detail { + +template MetricType, class Elem> +KOKKOS_FUNCTION double inverse_metric_component(MetricType metric, Elem elem) +{ + using metric_component_index_t = ddc::type_seq_element_t< + 0, + ddc::to_type_seq_t>; + if constexpr (misc::Specialization) { + return std::is_same_v ? 1. : 0.; + } else if constexpr ( + misc::Specialization) { + if constexpr (std::is_same_v) { + return metric( + elem, + metric.accessor().template access_element()); + } else { + return 0.; + } + } else { + return metric( + elem, + metric.accessor().template access_element()); + } +} + +template MetricType, class Elem> +KOKKOS_FUNCTION double primal_volume_factor(MetricType metric, Elem elem) +{ + using metric_index_1_t = tensor::metric_index_1>; + using axis1_t = ddc::type_seq_element_t<0, typename metric_index_1_t::type_seq_dimensions>; + using axis2_t = ddc::type_seq_element_t<1, typename metric_index_1_t::type_seq_dimensions>; + double const gxx = inverse_metric_component(metric, elem); + double const gxy = inverse_metric_component(metric, elem); + double const gyx = inverse_metric_component(metric, elem); + double const gyy = inverse_metric_component(metric, elem); + double const det_inv_metric = gxx * gyy - gxy * gyx; + return Kokkos::sqrt(Kokkos::abs(1. / det_inv_metric)); +} + +template +struct IsFullyDualDomain; + +template +struct IsFullyDualDomain> +{ + static constexpr bool value + = (true + && ... && !std::is_same_v>); +}; + +template +constexpr bool is_fully_dual_domain_v = IsFullyDualDomain::value; + +} // namespace detail + template < misc::Specialization Indices1, misc::Specialization Indices2> @@ -227,6 +284,55 @@ auto hodge_star( return out_form; } +template < + tensor::TensorIndex MetricIndex, + misc::Specialization OutTensorType, + misc::Specialization InTensorType, + misc::Specialization MetricType, + class ExecSpace> + requires( + ddc::type_seq_element_t<0, ddc::to_type_seq_t>::rank() + == 0 + && ddc::type_seq_element_t< + 0, + ddc::to_type_seq_t>::rank() + == 0 + && OutTensorType::non_indices_domain_t::rank() == 2 + && InTensorType::non_indices_domain_t::rank() == 2) +OutTensorType hodge_star( + ExecSpace const& exec_space, + OutTensorType out_tensor, + InTensorType in_tensor, + MetricType metric) +{ + using out_index_t = ddc::type_seq_element_t< + 0, + ddc::to_type_seq_t>; + using in_index_t = ddc::type_seq_element_t< + 0, + ddc::to_type_seq_t>; + using out_elem_t = typename OutTensorType::non_indices_domain_t::discrete_element_type; + using in_elem_t = typename InTensorType::non_indices_domain_t::discrete_element_type; + using metric_elem_t = typename MetricType::non_indices_domain_t::discrete_element_type; + constexpr bool in_is_fully_dual = detail::is_fully_dual_domain_v; + + ddc::parallel_for_each( + "similie_compute_scalar_tensor_hodge_star_2d", + exec_space, + out_tensor.non_indices_domain(), + KOKKOS_LAMBDA(out_elem_t out_elem) { + in_elem_t const in_elem + = sil::mesher::detail::dualizer_remap_element(out_elem); + metric_elem_t const metric_elem + = sil::mesher::detail::dualizer_remap_element(out_elem); + double const volume_factor = detail::primal_volume_factor(metric, metric_elem); + double const scale = in_is_fully_dual ? 1. / volume_factor : volume_factor; + out_tensor.mem(out_elem, ddc::DiscreteElement(0)) + = scale * in_tensor.get(in_elem, ddc::DiscreteElement(0)); + }); + return out_tensor; +} + } // namespace exterior } // namespace sil diff --git a/tests/exterior/hodge_star.cpp b/tests/exterior/hodge_star.cpp index 5a67b21a..11caa530 100644 --- a/tests/exterior/hodge_star.cpp +++ b/tests/exterior/hodge_star.cpp @@ -232,3 +232,94 @@ TEST(HodgeStar, TensorForm1In2D) primal_back_y(y_face_dom.front(), ddc::DiscreteElement(0)), -3.); } + +TEST(HodgeStar, ScalarTensor0And2In2D) +{ + struct X2 + { + }; + struct Y2 + { + }; + struct DDimX2 : ddc::UniformPointSampling + { + }; + struct DDimY2 : ddc::UniformPointSampling + { + }; + using ScalarIndex = sil::tensor::Covariant; + using Metric2DIndex = sil::tensor::TensorDiagonalIndex< + sil::tensor::Contravariant>, + sil::tensor::Contravariant>>; + using XDualizer = sil::mesher::HalfShiftDualizer; + using YDualizer = sil::mesher::HalfShiftDualizer; + using DDimXDual = sil::mesher::dual_discrete_dimension_t; + using DDimYDual = sil::mesher::dual_discrete_dimension_t; + + auto const x_dom = ddc::init_discrete_space(DDimX2::init( + ddc::Coordinate(0.), + ddc::Coordinate(1.), + ddc::DiscreteVector(4))); + auto const y_dom = ddc::init_discrete_space(DDimY2::init( + ddc::Coordinate(0.), + ddc::Coordinate(1.), + ddc::DiscreteVector(4))); + ddc::DiscreteDomain const mesh(x_dom, y_dom); + XDualizer const x_dualizer; + YDualizer const y_dualizer; + auto const dual_mesh = y_dualizer(x_dualizer(mesh)); + + [[maybe_unused]] sil::tensor::TensorAccessor scalar_accessor; + ddc::Chunk scalar_alloc( + ddc::DiscreteDomain(mesh, scalar_accessor.domain()), + ddc::HostAllocator()); + sil::tensor::Tensor scalar(scalar_alloc); + ddc::parallel_for_each( + Kokkos::DefaultHostExecutionSpace(), + scalar.domain(), + KOKKOS_LAMBDA(ddc::DiscreteElement elem) { scalar(elem) = 2.; }); + + [[maybe_unused]] sil::tensor::TensorAccessor metric_accessor; + ddc::Chunk metric_alloc( + ddc::DiscreteDomain(mesh, metric_accessor.domain()), + ddc::HostAllocator()); + sil::tensor::Tensor metric(metric_alloc); + ddc::parallel_for_each( + Kokkos::DefaultHostExecutionSpace(), + mesh, + KOKKOS_LAMBDA(ddc::DiscreteElement elem) { + metric(elem, metric.accessor().access_element()) = 4.; + metric(elem, metric.accessor().access_element()) = 9.; + }); + + ddc::Chunk dual_scalar_alloc( + ddc::DiscreteDomain(dual_mesh, scalar_accessor.domain()), + ddc::HostAllocator()); + sil::tensor::Tensor dual_scalar(dual_scalar_alloc); + + sil::exterior::hodge_star( + Kokkos::DefaultHostExecutionSpace(), + dual_scalar, + scalar, + metric); + + EXPECT_DOUBLE_EQ(dual_scalar(dual_scalar.domain().front()), 1. / 3.); + + ddc::DiscreteDomain const primal_submesh( + mesh.front(), + ddc::DiscreteVector( + dual_mesh.template extent().value(), + dual_mesh.template extent().value())); + ddc::Chunk scalar_back_alloc( + ddc::DiscreteDomain(primal_submesh, scalar_accessor.domain()), + ddc::HostAllocator()); + sil::tensor::Tensor scalar_back(scalar_back_alloc); + ddc::parallel_fill(scalar_back, 0.); + sil::exterior::hodge_star( + Kokkos::DefaultHostExecutionSpace(), + scalar_back, + dual_scalar, + metric); + + EXPECT_DOUBLE_EQ(scalar_back(scalar_back.domain().front()), 2.); +} From 7afaf0021485a3bc4668b06e0ad3b109d9c68a92 Mon Sep 17 00:00:00 2001 From: blegouix Date: Sun, 29 Mar 2026 03:15:00 +0200 Subject: [PATCH 09/14] wip --- src/similie/exterior/codifferential.hpp | 63 ++++++----- src/similie/exterior/hodge_star.hpp | 142 ++++++++++++++++-------- tests/exterior/exterior.cpp | 54 ++++----- tests/exterior/hodge_star.cpp | 64 +++++------ 4 files changed, 181 insertions(+), 142 deletions(-) diff --git a/src/similie/exterior/codifferential.hpp b/src/similie/exterior/codifferential.hpp index 83bbeba6..3ce2990c 100644 --- a/src/similie/exterior/codifferential.hpp +++ b/src/similie/exterior/codifferential.hpp @@ -162,6 +162,36 @@ template constexpr bool is_dual_discrete_dimension_v = !std::is_same_v>; +template +KOKKOS_FUNCTION bool is_on_primal_boundary(Domain const& domain, Elem const& elem) +{ + if constexpr (is_dual_discrete_dimension_v) { + return false; + } else { + ddc::DiscreteDomain const dim_dom(domain); + ddc::DiscreteElement const dim_elem(elem); + return dim_elem.uid() == dim_dom.front().uid() || dim_elem.uid() == dim_dom.back().uid(); + } +} + +template < + class MetricComponentIndex, + class TagToRemoveFromCochain, + misc::Specialization MetricType, + class Elem> +KOKKOS_FUNCTION double codifferential_metric_factor(MetricType const& inv_metric, Elem const& elem) +{ + if constexpr (misc::Specialization) { + return 1.; + } else { + return inv_metric( + elem, + inv_metric.accessor().template access_element< + TagToRemoveFromCochain, + TagToRemoveFromCochain>()); + } +} + template < class AxisTag, misc::Specialization OutTensorType, @@ -232,21 +262,9 @@ OutTensorType exterior_derivative_of_tensor_form_2d( exec_space, out_tensor.non_indices_domain(), KOKKOS_LAMBDA(typename OutTensorType::non_indices_domain_t::discrete_element_type elem) { - if constexpr (!is_dual_discrete_dimension_v) { - ddc::DiscreteDomain dim_dom(out_tensor.non_indices_domain()); - ddc::DiscreteElement dim_elem(elem); - if (dim_elem.uid() == dim_dom.front().uid() - || dim_elem.uid() == dim_dom.back().uid()) { - return; - } - } - if constexpr (!is_dual_discrete_dimension_v) { - ddc::DiscreteDomain dim_dom(out_tensor.non_indices_domain()); - ddc::DiscreteElement dim_elem(elem); - if (dim_elem.uid() == dim_dom.front().uid() - || dim_elem.uid() == dim_dom.back().uid()) { - return; - } + if (is_on_primal_boundary(out_tensor.non_indices_domain(), elem) + || is_on_primal_boundary(out_tensor.non_indices_domain(), elem)) { + return; } out_tensor.mem(elem, ddc::DiscreteElement(0)) += centered_form_derivative( @@ -436,18 +454,9 @@ OutTensorType codifferential( double const deriv = (tensor.get(right_face, ddc::DiscreteElement(0)) - tensor.get(left_face, ddc::DiscreteElement(0))) / dx; - double const metric_factor = [&]() -> double { - if constexpr ( - misc::Specialization) { - return 1.; - } else { - return inv_metric( - out_elem, - inv_metric.accessor().template access_element< - TagToRemoveFromCochain, - TagToRemoveFromCochain>()); - } - }(); + double const metric_factor = detail::codifferential_metric_factor< + metric_component_index_t, + TagToRemoveFromCochain>(inv_metric, out_elem); out_tensor.mem(out_elem, ddc::DiscreteElement(0)) += metric_factor * deriv; }); diff --git a/src/similie/exterior/hodge_star.hpp b/src/similie/exterior/hodge_star.hpp index e84f5574..e60ad3ae 100644 --- a/src/similie/exterior/hodge_star.hpp +++ b/src/similie/exterior/hodge_star.hpp @@ -78,6 +78,95 @@ struct IsFullyDualDomain> template constexpr bool is_fully_dual_domain_v = IsFullyDualDomain::value; +template +struct TensorFormHodgeStarFirstComponentFunctor +{ + OutTensor out_tensor; + InTensor in_tensor; + + KOKKOS_FUNCTION void operator()(typename InTensor::non_indices_domain_t::discrete_element_type elem) const + { + out_tensor.mem(elem, ddc::DiscreteElement(0)) + = -in_tensor.get(elem, ddc::DiscreteElement(0)); + } +}; + +template +struct TensorFormHodgeStarSecondComponentFunctor +{ + OutTensor out_tensor; + InTensor in_tensor; + + KOKKOS_FUNCTION void operator()(typename InTensor::non_indices_domain_t::discrete_element_type elem) const + { + out_tensor.mem(elem, ddc::DiscreteElement(0)) + = in_tensor.get(elem, ddc::DiscreteElement(0)); + } +}; + +template < + class FirstTag, + class SecondTag, + class SupportTag, + class... OutComponents, + class FirstComponent, + class SecondComponent, + class ExecSpace> +void fill_tensor_form_hodge_star_2d( + ExecSpace const& exec_space, + TensorForm, OutComponents...> out_form, + TensorForm in_form) +{ + using first_in_tensor_t = typename FirstComponent::tensor_type; + using second_in_tensor_t = typename SecondComponent::tensor_type; + using out_first_component_t = typename FormComponentByTag::type; + using out_second_component_t = typename FormComponentByTag::type; + using first_in_index_t = ddc::type_seq_element_t< + 0, + ddc::to_type_seq_t>; + using second_in_index_t = ddc::type_seq_element_t< + 0, + ddc::to_type_seq_t>; + using first_out_index_t = ddc::type_seq_element_t< + 0, + ddc::to_type_seq_t>; + using second_out_index_t = ddc::type_seq_element_t< + 0, + ddc::to_type_seq_t>; + static_assert(!std::is_void_v); + static_assert(!std::is_void_v); + static_assert(std::is_same_v); + static_assert(std::is_same_v); + static_assert(first_in_index_t::rank() == 0); + static_assert(second_in_index_t::rank() == 0); + static_assert(first_out_index_t::rank() == 0); + static_assert(second_out_index_t::rank() == 0); + + auto const in_first = in_form.template component(); + auto const in_second = in_form.template component(); + auto const out_first = out_form.template component(); + auto const out_second = out_form.template component(); + + ddc::parallel_for_each( + "similie_compute_tensor_form_hodge_star_first_component", + exec_space, + out_first.non_indices_domain(), + TensorFormHodgeStarFirstComponentFunctor< + decltype(out_first), + decltype(in_second), + first_out_index_t, + second_in_index_t> {out_first, in_second}); + ddc::parallel_for_each( + "similie_compute_tensor_form_hodge_star_second_component", + exec_space, + out_second.non_indices_domain(), + TensorFormHodgeStarSecondComponentFunctor< + decltype(out_second), + decltype(in_first), + second_out_index_t, + first_in_index_t> {out_second, in_first}); +} + } // namespace detail template < @@ -231,56 +320,11 @@ auto hodge_star( ExecSpace const& exec_space, TensorForm, OutComponents...> out_form, TensorForm in_form, - MetricType const&) + MetricType const&) -> TensorForm, OutComponents...> { - using first_in_tensor_t = typename FirstComponent::tensor_type; - using second_in_tensor_t = typename SecondComponent::tensor_type; - using out_first_component_t = typename detail:: - FormComponentByTag::type; - using out_second_component_t = typename detail:: - FormComponentByTag::type; - using first_in_index_t = ddc::type_seq_element_t< - 0, - ddc::to_type_seq_t>; - using second_in_index_t = ddc::type_seq_element_t< - 0, - ddc::to_type_seq_t>; - using first_out_index_t = ddc::type_seq_element_t< - 0, - ddc::to_type_seq_t>; - using second_out_index_t = ddc::type_seq_element_t< - 0, - ddc::to_type_seq_t>; - static_assert(!std::is_void_v); - static_assert(!std::is_void_v); - static_assert(std::is_same_v); - static_assert(std::is_same_v); - static_assert(first_in_index_t::rank() == 0); - static_assert(second_in_index_t::rank() == 0); - static_assert(first_out_index_t::rank() == 0); - static_assert(second_out_index_t::rank() == 0); - - auto const in_first = in_form.template component(); - auto const in_second = in_form.template component(); - auto const out_first = out_form.template component(); - auto const out_second = out_form.template component(); - - ddc::parallel_for_each( - "similie_compute_tensor_form_hodge_star_first_component", - exec_space, - out_first.non_indices_domain(), - KOKKOS_LAMBDA(typename second_in_tensor_t::non_indices_domain_t::discrete_element_type elem) { - out_first.mem(elem, ddc::DiscreteElement(0)) - = -in_second.get(elem, ddc::DiscreteElement(0)); - }); - ddc::parallel_for_each( - "similie_compute_tensor_form_hodge_star_second_component", - exec_space, - out_second.non_indices_domain(), - KOKKOS_LAMBDA(typename first_in_tensor_t::non_indices_domain_t::discrete_element_type elem) { - out_second.mem(elem, ddc::DiscreteElement(0)) - = in_first.get(elem, ddc::DiscreteElement(0)); - }); + detail::fill_tensor_form_hodge_star_2d< + typename FirstComponent::tag, + typename SecondComponent::tag>(exec_space, out_form, in_form); return out_form; } diff --git a/tests/exterior/exterior.cpp b/tests/exterior/exterior.cpp index bc27f123..87a2578e 100644 --- a/tests/exterior/exterior.cpp +++ b/tests/exterior/exterior.cpp @@ -41,6 +41,14 @@ struct DDimZ : ddc::UniformPointSampling { }; +struct DDimX2 : ddc::UniformPointSampling +{ +}; + +struct DDimY2 : ddc::UniformPointSampling +{ +}; + TEST(Chain, Optimization) { sil::exterior::Chain chain = sil::exterior:: @@ -541,14 +549,11 @@ TEST(Form, TensorFormDeriv) ddc::DiscreteDomain(mesh, scalar_accessor.domain()), ddc::HostAllocator()); sil::tensor::Tensor scalar(scalar_alloc); - ddc::parallel_for_each( - Kokkos::DefaultHostExecutionSpace(), - mesh, - KOKKOS_LAMBDA(ddc::DiscreteElement elem) { - double const x = ddc::coordinate(ddc::DiscreteElement(elem)); - double const y = ddc::coordinate(ddc::DiscreteElement(elem)); - scalar(elem, ddc::DiscreteElement(0)) = x * x + y; - }); + ddc::host_for_each(mesh, [&](ddc::DiscreteElement elem) { + double const x = ddc::coordinate(ddc::DiscreteElement(elem)); + double const y = ddc::coordinate(ddc::DiscreteElement(elem)); + scalar(elem, ddc::DiscreteElement(0)) = x * x + y; + }); ddc::Chunk grad_x_alloc( ddc::DiscreteDomain( @@ -578,12 +583,6 @@ TEST(Form, TensorFormDeriv) TEST(Form, TensorFormCodifferential) { - struct DDimX2 : ddc::UniformPointSampling - { - }; - struct DDimY2 : ddc::UniformPointSampling - { - }; using MetricIndex = sil::tensor::TensorIdentityIndex< sil::tensor::Covariant>, sil::tensor::Covariant>>; @@ -622,19 +621,13 @@ TEST(Form, TensorFormCodifferential) ddc::HostAllocator()); sil::tensor::Tensor grad_y(grad_y_alloc); - ddc::parallel_for_each( - Kokkos::DefaultHostExecutionSpace(), - grad_x.domain(), - KOKKOS_LAMBDA(ddc::DiscreteElement elem) { - double const x = ddc::coordinate(ddc::DiscreteElement(elem)); - grad_x(elem) = 2. * x; - }); - ddc::parallel_for_each( - Kokkos::DefaultHostExecutionSpace(), - grad_y.domain(), - KOKKOS_LAMBDA(ddc::DiscreteElement elem) { - grad_y(elem) = 1.; - }); + ddc::host_for_each(grad_x.domain(), [&](ddc::DiscreteElement elem) { + double const x = ddc::coordinate(ddc::DiscreteElement(elem)); + grad_x(elem) = 2. * x; + }); + ddc::host_for_each(grad_y.domain(), [&](ddc::DiscreteElement elem) { + grad_y(elem) = 1.; + }); auto form = sil::exterior::make_tensor_form( sil::exterior::component(grad_x), @@ -652,10 +645,9 @@ TEST(Form, TensorFormCodifferential) ddc::DiscreteDomain(mesh, scalar_accessor.domain()), ddc::HostAllocator()); sil::tensor::Tensor div(div_alloc); - ddc::parallel_for_each( - Kokkos::DefaultHostExecutionSpace(), - div.domain(), - KOKKOS_LAMBDA(ddc::DiscreteElement elem) { div(elem) = 0.; }); + ddc::host_for_each(div.domain(), [&](ddc::DiscreteElement elem) { + div(elem) = 0.; + }); sil::exterior::codifferential( Kokkos::DefaultHostExecutionSpace(), diff --git a/tests/exterior/hodge_star.cpp b/tests/exterior/hodge_star.cpp index 11caa530..da345bd4 100644 --- a/tests/exterior/hodge_star.cpp +++ b/tests/exterior/hodge_star.cpp @@ -36,6 +36,22 @@ struct DDimZ : ddc::UniformPointSampling { }; +struct X2 +{ +}; + +struct Y2 +{ +}; + +struct DDimX2 : ddc::UniformPointSampling +{ +}; + +struct DDimY2 : ddc::UniformPointSampling +{ +}; + struct Mu : sil::tensor::TensorNaturalIndex { }; @@ -173,18 +189,12 @@ TEST(HodgeStar, TensorForm1In2D) ddc::HostAllocator()); sil::tensor::Tensor primal_y(primal_y_alloc); - ddc::parallel_for_each( - Kokkos::DefaultHostExecutionSpace(), - primal_x.domain(), - KOKKOS_LAMBDA(ddc::DiscreteElement elem) { - primal_x(elem) = 2.; - }); - ddc::parallel_for_each( - Kokkos::DefaultHostExecutionSpace(), - primal_y.domain(), - KOKKOS_LAMBDA(ddc::DiscreteElement elem) { - primal_y(elem) = 3.; - }); + ddc::host_for_each(primal_x.domain(), [&](ddc::DiscreteElement elem) { + primal_x(elem) = 2.; + }); + ddc::host_for_each(primal_y.domain(), [&](ddc::DiscreteElement elem) { + primal_y(elem) = 3.; + }); auto primal_form = sil::exterior::make_tensor_form( sil::exterior::component(primal_x), @@ -235,18 +245,6 @@ TEST(HodgeStar, TensorForm1In2D) TEST(HodgeStar, ScalarTensor0And2In2D) { - struct X2 - { - }; - struct Y2 - { - }; - struct DDimX2 : ddc::UniformPointSampling - { - }; - struct DDimY2 : ddc::UniformPointSampling - { - }; using ScalarIndex = sil::tensor::Covariant; using Metric2DIndex = sil::tensor::TensorDiagonalIndex< sil::tensor::Contravariant>, @@ -274,23 +272,19 @@ TEST(HodgeStar, ScalarTensor0And2In2D) ddc::DiscreteDomain(mesh, scalar_accessor.domain()), ddc::HostAllocator()); sil::tensor::Tensor scalar(scalar_alloc); - ddc::parallel_for_each( - Kokkos::DefaultHostExecutionSpace(), - scalar.domain(), - KOKKOS_LAMBDA(ddc::DiscreteElement elem) { scalar(elem) = 2.; }); + ddc::host_for_each(scalar.domain(), [&](ddc::DiscreteElement elem) { + scalar(elem) = 2.; + }); [[maybe_unused]] sil::tensor::TensorAccessor metric_accessor; ddc::Chunk metric_alloc( ddc::DiscreteDomain(mesh, metric_accessor.domain()), ddc::HostAllocator()); sil::tensor::Tensor metric(metric_alloc); - ddc::parallel_for_each( - Kokkos::DefaultHostExecutionSpace(), - mesh, - KOKKOS_LAMBDA(ddc::DiscreteElement elem) { - metric(elem, metric.accessor().access_element()) = 4.; - metric(elem, metric.accessor().access_element()) = 9.; - }); + ddc::host_for_each(mesh, [&](ddc::DiscreteElement elem) { + metric(elem, metric.accessor().access_element()) = 4.; + metric(elem, metric.accessor().access_element()) = 9.; + }); ddc::Chunk dual_scalar_alloc( ddc::DiscreteDomain(dual_mesh, scalar_accessor.domain()), From 4b836b8be81fd2f88fd2cdd2df5ae19401a3a116 Mon Sep 17 00:00:00 2001 From: blegouix Date: Sun, 29 Mar 2026 03:51:42 +0200 Subject: [PATCH 10/14] wip --- .../2d_free_scalar_field.cpp | 33 ++++++++++--------- 1 file changed, 17 insertions(+), 16 deletions(-) diff --git a/examples/free_scalar_field/2d_free_scalar_field.cpp b/examples/free_scalar_field/2d_free_scalar_field.cpp index 194c42fe..980eb069 100644 --- a/examples/free_scalar_field/2d_free_scalar_field.cpp +++ b/examples/free_scalar_field/2d_free_scalar_field.cpp @@ -378,9 +378,9 @@ int main(int argc, char** argv) double const y_1 = -0.3; double const sigma = .5; - double const v = 10.; - double const k = 1.; - double const mass = 1.; + double const v = 0.9; + double const mass = 10; + double const k = mass * v / std::sqrt(1. - v * v); ddc::parallel_for_each( potential.domain(), @@ -392,10 +392,9 @@ int main(int argc, char** argv) * std::exp( -((x - x_0) * (x - x_0) + (y - y_0) * (y - y_0)) / 2. / sigma / sigma) - + 30. * std::sin(k * (x - x_1)) - * std::exp( - -((x - x_1) * (x - x_1) + (y - y_1) * (y - y_1)) - / 2. / sigma / sigma); + + std::exp( + -((x - x_1) * (x - x_1) + (y - y_1) * (y - y_1)) / 2. + / sigma / sigma); }); auto potential_host = ddc::create_mirror_view_and_copy(Kokkos::DefaultHostExecutionSpace(), potential); @@ -468,8 +467,8 @@ int main(int argc, char** argv) KOKKOS_LAMBDA(ddc::DiscreteElement elem) { double const x = ddc::coordinate(ddc::DiscreteElement(elem)); double const y = ddc::coordinate(ddc::DiscreteElement(elem)) - y_0; - // v*dphi/dx of the left wave packet only to get a pure kick along x toward the immobile right one - temporal_moment(elem) = -v + // -vg * dphi/dx of the left wave packet only, to launch it toward +x + temporal_moment(elem) = -k / std::sqrt(k * k + mass * mass) * (-k * std::cos(k * (x - x_0)) + (x - x_0) / sigma / sigma * std::sin(k * (x - x_0))) * std::exp( @@ -561,11 +560,11 @@ int main(int argc, char** argv) KOKKOS_LAMBDA(ddc::DiscreteElement elem) { spatial_moments_div(elem) = 0.; }); - sil::exterior::codifferential< - MetricIndex>(Kokkos::DefaultExecutionSpace(), - spatial_moments_div, - spatial_moment, - inv_metric); + sil::exterior::codifferential( + Kokkos::DefaultExecutionSpace(), + spatial_moments_div, + spatial_moment, + inv_metric); // Compute dpi_0/dx^0 by solving - dpi_0/dx^0 + dpi_\alpha/dx^\alpha = -dH/dphi and advect pi_0 by a time step dx^0. Also Then, perform the second phi half-advection by solving dphi/dx^0 = -dH/dpi_0 double const dS = (ddc::get(upper_bounds) - ddc::get(lower_bounds)) @@ -620,7 +619,8 @@ int main(int argc, char** argv) for (std::size_t iy = 0; iy < x_face_dom.template extent(); ++iy) { ddc::DiscreteElement const elem = x_face_front + ddc::DiscreteVector(ix, iy); - ddc::DiscreteElement const left_node = x_dualizer.primal(elem); + ddc::DiscreteElement const left_node + = x_dualizer.primal(elem); double const value = spatial_moment_x_host(elem, ddc::DiscreteElement(0)); spatial_moments_host(left_node, ddc::DiscreteElement(0)) @@ -638,7 +638,8 @@ int main(int argc, char** argv) for (std::size_t iy = 0; iy < y_face_dom.template extent(); ++iy) { ddc::DiscreteElement const elem = y_face_front + ddc::DiscreteVector(ix, iy); - ddc::DiscreteElement const lower_node = y_dualizer.primal(elem); + ddc::DiscreteElement const lower_node + = y_dualizer.primal(elem); double const value = spatial_moment_y_host(elem, ddc::DiscreteElement(0)); spatial_moments_host(lower_node, ddc::DiscreteElement(1)) From df5157f7aaa861f84a12321e3eef669cd158c152 Mon Sep 17 00:00:00 2001 From: blegouix Date: Sun, 29 Mar 2026 03:55:58 +0200 Subject: [PATCH 11/14] wip --- examples/free_scalar_field/2d_free_scalar_field.cpp | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/examples/free_scalar_field/2d_free_scalar_field.cpp b/examples/free_scalar_field/2d_free_scalar_field.cpp index 980eb069..2199ffbe 100644 --- a/examples/free_scalar_field/2d_free_scalar_field.cpp +++ b/examples/free_scalar_field/2d_free_scalar_field.cpp @@ -462,15 +462,14 @@ int main(int argc, char** argv) ddc::Chunk temporal_moment_alloc(temporal_moment_dom, ddc::DeviceAllocator()); sil::tensor::Tensor temporal_moment(temporal_moment_alloc); + double const omega = std::sqrt(k * k + mass * mass); ddc::parallel_for_each( potential.domain(), KOKKOS_LAMBDA(ddc::DiscreteElement elem) { double const x = ddc::coordinate(ddc::DiscreteElement(elem)); double const y = ddc::coordinate(ddc::DiscreteElement(elem)) - y_0; - // -vg * dphi/dx of the left wave packet only, to launch it toward +x - temporal_moment(elem) = -k / std::sqrt(k * k + mass * mass) - * (-k * std::cos(k * (x - x_0)) - + (x - x_0) / sigma / sigma * std::sin(k * (x - x_0))) + temporal_moment(elem) = (-omega * std::cos(k * (x - x_0)) + + v * (x - x_0) / sigma / sigma * std::sin(k * (x - x_0))) * std::exp( -((x - x_0) * (x - x_0) + (y - y_0) * (y - y_0)) / 2. / sigma / sigma); From 2ea5b2ae3ebb531503306b2bcf95d32454c9bb82 Mon Sep 17 00:00:00 2001 From: blegouix Date: Sun, 29 Mar 2026 23:27:38 +0200 Subject: [PATCH 12/14] wip --- examples/free_scalar_field/2d_free_scalar_field.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/free_scalar_field/2d_free_scalar_field.cpp b/examples/free_scalar_field/2d_free_scalar_field.cpp index 2199ffbe..28bda0a3 100644 --- a/examples/free_scalar_field/2d_free_scalar_field.cpp +++ b/examples/free_scalar_field/2d_free_scalar_field.cpp @@ -378,8 +378,8 @@ int main(int argc, char** argv) double const y_1 = -0.3; double const sigma = .5; - double const v = 0.9; - double const mass = 10; + double const v = 0.5; + double const mass = 100; double const k = mass * v / std::sqrt(1. - v * v); ddc::parallel_for_each( From 36a9d6c07da41b6d051d907a930a141343520677 Mon Sep 17 00:00:00 2001 From: blegouix Date: Mon, 30 Mar 2026 00:04:40 +0200 Subject: [PATCH 13/14] wip --- .../2d_free_scalar_field.cpp | 32 ++++++++++++------- ...lar_field_hamiltonian_functor_generator.py | 12 ++++--- 2 files changed, 28 insertions(+), 16 deletions(-) diff --git a/examples/free_scalar_field/2d_free_scalar_field.cpp b/examples/free_scalar_field/2d_free_scalar_field.cpp index 28bda0a3..9a701323 100644 --- a/examples/free_scalar_field/2d_free_scalar_field.cpp +++ b/examples/free_scalar_field/2d_free_scalar_field.cpp @@ -380,6 +380,7 @@ int main(int argc, char** argv) double const v = 0.5; double const mass = 100; + double const quartic_coupling_constant = 0; double const k = mass * v / std::sqrt(1. - v * v); ddc::parallel_for_each( @@ -527,8 +528,9 @@ int main(int argc, char** argv) KOKKOS_LAMBDA(ddc::DiscreteElement elem) { half_step_potential(elem) = potential(elem) - - FreeScalarFieldHamiltonian(mass).dH_dpi0(temporal_moment(elem)) * dt - / 2.; + - FreeScalarFieldHamiltonian(mass, quartic_coupling_constant) + .dH_dpi0(temporal_moment(elem)) + * dt / 2.; }); // Compute the staggered potential gradients @@ -540,16 +542,20 @@ int main(int argc, char** argv) x_face_dom, KOKKOS_LAMBDA(ddc::DiscreteElement elem) { spatial_moment_x(elem, ddc::DiscreteElement(0)) - = FreeScalarFieldHamiltonian(mass).pi1( - potential_grad_x(elem, ddc::DiscreteElement(0))); + = FreeScalarFieldHamiltonian(mass, quartic_coupling_constant) + .pi1(potential_grad_x( + elem, + ddc::DiscreteElement(0))); }); ddc::parallel_for_each( Kokkos::DefaultExecutionSpace(), y_face_dom, KOKKOS_LAMBDA(ddc::DiscreteElement elem) { spatial_moment_y(elem, ddc::DiscreteElement(0)) - = FreeScalarFieldHamiltonian(mass).pi2( - potential_grad_y(elem, ddc::DiscreteElement(0))); + = FreeScalarFieldHamiltonian(mass, quartic_coupling_constant) + .pi2(potential_grad_y( + elem, + ddc::DiscreteElement(0))); }); // Compute the divergence dpi_\alpha/dx^\alpha of the spatial moments, which is the codifferential \delta pi of the spatial moments @@ -581,16 +587,18 @@ int main(int argc, char** argv) const double half_step_temporal_moment_ = temporal_moment_ + (spatial_moments_div_ - - FreeScalarFieldHamiltonian(mass).dH_dphi(half_step_potential_)) + - FreeScalarFieldHamiltonian(mass, quartic_coupling_constant) + .dH_dphi(half_step_potential_)) * dt / 2; // Whole-step advection of field state - potential(elem) - -= FreeScalarFieldHamiltonian(mass).dH_dpi0(half_step_temporal_moment_) - * dt; + potential(elem) -= FreeScalarFieldHamiltonian(mass, quartic_coupling_constant) + .dH_dpi0(half_step_temporal_moment_) + * dt; temporal_moment(elem) += (spatial_moments_div_ - - FreeScalarFieldHamiltonian(mass).dH_dphi(half_step_potential_)) + - FreeScalarFieldHamiltonian(mass, quartic_coupling_constant) + .dH_dphi(half_step_potential_)) * dt; }); @@ -661,7 +669,7 @@ int main(int argc, char** argv) spatial_moments(elem, ddc::DiscreteElement(0)), spatial_moments(elem, ddc::DiscreteElement(1))}; hamiltonian(elem) - = FreeScalarFieldHamiltonian(mass) + = FreeScalarFieldHamiltonian(mass, quartic_coupling_constant) .H(potential(elem, ddc::DiscreteElement()), pi); }); diff --git a/examples/free_scalar_field/free_scalar_field_hamiltonian_functor_generator.py b/examples/free_scalar_field/free_scalar_field_hamiltonian_functor_generator.py index 772e0e46..c7acf020 100644 --- a/examples/free_scalar_field/free_scalar_field_hamiltonian_functor_generator.py +++ b/examples/free_scalar_field/free_scalar_field_hamiltonian_functor_generator.py @@ -14,17 +14,17 @@ output_path = Path(sys.argv[2]) # Path of the output file output_path.parent.mkdir(parents=True, exist_ok=True) -mass = symbols("mass") +mass, quartic_coupling_constant = symbols("mass quartic_coupling_constant") phi = symbols("phi") pi = symbols(f"pi0:{N}") # Minkowski signature (+, -, -, ..., -) metric_sign = [-1] + [1] * (N - 1) -# H = 1/2 m^2 phi^2 + 1/2 (-p0^2 + p1^2 - ... - p{N-1}^2) +# H = 1/2 m^2 phi^2 + lambda/4 phi^4 + 1/2 (-p0^2 + p1^2 - ... - p{N-1}^2) hamiltonian = 0.5 * ( mass**2 * phi**2 + sum(metric_sign[i] * pi[i] ** 2 for i in range(N)) -) +) + quartic_coupling_constant * phi**4 / 4 # [dH/dphi, dH/dpi0, dH/dpi1, ...] hamiltonian_diff = Matrix( @@ -69,8 +69,12 @@ def preprocess_cxx(expr: str) -> str: static constexpr std::size_t N = {N}; const double mass; + const double quartic_coupling_constant; - constexpr FreeScalarFieldHamiltonian(double mass_) : mass(mass_) {{}} + constexpr FreeScalarFieldHamiltonian(double mass_, double quartic_coupling_constant_) + : mass(mass_) + , quartic_coupling_constant(quartic_coupling_constant_) + {{}} constexpr double H(double phi, const std::span& pi) const {{ From 71d6d20e54d94af8034ae01341492b9dcd55b2f0 Mon Sep 17 00:00:00 2001 From: blegouix Date: Tue, 31 Mar 2026 18:54:10 +0200 Subject: [PATCH 14/14] coupling --- .../2d_free_scalar_field.cpp | 39 +++++++++++++++---- ...lar_field_hamiltonian_functor_generator.py | 16 ++++++-- 2 files changed, 43 insertions(+), 12 deletions(-) diff --git a/examples/free_scalar_field/2d_free_scalar_field.cpp b/examples/free_scalar_field/2d_free_scalar_field.cpp index 9a701323..35521457 100644 --- a/examples/free_scalar_field/2d_free_scalar_field.cpp +++ b/examples/free_scalar_field/2d_free_scalar_field.cpp @@ -380,6 +380,7 @@ int main(int argc, char** argv) double const v = 0.5; double const mass = 100; + double const linear_coupling_constant = 0; double const quartic_coupling_constant = 0; double const k = mass * v / std::sqrt(1. - v * v); @@ -528,7 +529,10 @@ int main(int argc, char** argv) KOKKOS_LAMBDA(ddc::DiscreteElement elem) { half_step_potential(elem) = potential(elem) - - FreeScalarFieldHamiltonian(mass, quartic_coupling_constant) + - FreeScalarFieldHamiltonian( + mass, + quartic_coupling_constant, + linear_coupling_constant) .dH_dpi0(temporal_moment(elem)) * dt / 2.; }); @@ -542,7 +546,10 @@ int main(int argc, char** argv) x_face_dom, KOKKOS_LAMBDA(ddc::DiscreteElement elem) { spatial_moment_x(elem, ddc::DiscreteElement(0)) - = FreeScalarFieldHamiltonian(mass, quartic_coupling_constant) + = FreeScalarFieldHamiltonian( + mass, + quartic_coupling_constant, + linear_coupling_constant) .pi1(potential_grad_x( elem, ddc::DiscreteElement(0))); @@ -552,7 +559,10 @@ int main(int argc, char** argv) y_face_dom, KOKKOS_LAMBDA(ddc::DiscreteElement elem) { spatial_moment_y(elem, ddc::DiscreteElement(0)) - = FreeScalarFieldHamiltonian(mass, quartic_coupling_constant) + = FreeScalarFieldHamiltonian( + mass, + quartic_coupling_constant, + linear_coupling_constant) .pi2(potential_grad_y( elem, ddc::DiscreteElement(0))); @@ -587,17 +597,27 @@ int main(int argc, char** argv) const double half_step_temporal_moment_ = temporal_moment_ + (spatial_moments_div_ - - FreeScalarFieldHamiltonian(mass, quartic_coupling_constant) + - FreeScalarFieldHamiltonian( + mass, + quartic_coupling_constant, + linear_coupling_constant) .dH_dphi(half_step_potential_)) * dt / 2; // Whole-step advection of field state - potential(elem) -= FreeScalarFieldHamiltonian(mass, quartic_coupling_constant) - .dH_dpi0(half_step_temporal_moment_) + potential(elem) + -= FreeScalarFieldHamiltonian( + mass, + quartic_coupling_constant, + linear_coupling_constant) + .dH_dpi0(half_step_temporal_moment_) * dt; temporal_moment(elem) += (spatial_moments_div_ - - FreeScalarFieldHamiltonian(mass, quartic_coupling_constant) + - FreeScalarFieldHamiltonian( + mass, + quartic_coupling_constant, + linear_coupling_constant) .dH_dphi(half_step_potential_)) * dt; }); @@ -669,7 +689,10 @@ int main(int argc, char** argv) spatial_moments(elem, ddc::DiscreteElement(0)), spatial_moments(elem, ddc::DiscreteElement(1))}; hamiltonian(elem) - = FreeScalarFieldHamiltonian(mass, quartic_coupling_constant) + = FreeScalarFieldHamiltonian( + mass, + quartic_coupling_constant, + linear_coupling_constant) .H(potential(elem, ddc::DiscreteElement()), pi); }); diff --git a/examples/free_scalar_field/free_scalar_field_hamiltonian_functor_generator.py b/examples/free_scalar_field/free_scalar_field_hamiltonian_functor_generator.py index c7acf020..08771123 100644 --- a/examples/free_scalar_field/free_scalar_field_hamiltonian_functor_generator.py +++ b/examples/free_scalar_field/free_scalar_field_hamiltonian_functor_generator.py @@ -14,17 +14,20 @@ output_path = Path(sys.argv[2]) # Path of the output file output_path.parent.mkdir(parents=True, exist_ok=True) -mass, quartic_coupling_constant = symbols("mass quartic_coupling_constant") +mass, quartic_coupling_constant, linear_coupling_constant = symbols( + "mass quartic_coupling_constant linear_coupling_constant" +) phi = symbols("phi") pi = symbols(f"pi0:{N}") # Minkowski signature (+, -, -, ..., -) metric_sign = [-1] + [1] * (N - 1) -# H = 1/2 m^2 phi^2 + lambda/4 phi^4 + 1/2 (-p0^2 + p1^2 - ... - p{N-1}^2) +# H = linear_coupling_constant * phi + 1/2 m^2 phi^2 + lambda/4 phi^4 +# + 1/2 (-p0^2 + p1^2 - ... - p{N-1}^2) hamiltonian = 0.5 * ( mass**2 * phi**2 + sum(metric_sign[i] * pi[i] ** 2 for i in range(N)) -) + quartic_coupling_constant * phi**4 / 4 +) + quartic_coupling_constant * phi**4 / 4 + linear_coupling_constant * phi # [dH/dphi, dH/dpi0, dH/dpi1, ...] hamiltonian_diff = Matrix( @@ -70,10 +73,15 @@ def preprocess_cxx(expr: str) -> str: const double mass; const double quartic_coupling_constant; + const double linear_coupling_constant; - constexpr FreeScalarFieldHamiltonian(double mass_, double quartic_coupling_constant_) + constexpr FreeScalarFieldHamiltonian( + double mass_, + double quartic_coupling_constant_, + double linear_coupling_constant_) : mass(mass_) , quartic_coupling_constant(quartic_coupling_constant_) + , linear_coupling_constant(linear_coupling_constant_) {{}} constexpr double H(double phi, const std::span& pi) const