Dual mesh

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 = symbols("mass")
+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 + 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 + linear_coupling_constant * phi
 
 # [dH/dphi, dH/dpi0, dH/dpi1, ...]
 hamiltonian_diff = Matrix(
@@ -69,8 +72,17 @@ def preprocess_cxx(expr: str) -> str:
     static constexpr std::size_t N = {N};
 
     const double mass;
-
-    constexpr FreeScalarFieldHamiltonian(double mass_) : mass(mass_) {{}}
+    const double quartic_coupling_constant;
+    const double linear_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<const double, N>& pi) const
     {{

src/similie/exterior/coboundary.hpp

@@ -12,6 +12,7 @@
 #include <similie/misc/portable_stl.hpp>
 #include <similie/misc/specialization.hpp>
 #include <similie/misc/type_seq_conversion.hpp>
+#include <similie/mesher/dualizer.hpp>
 #include <similie/tensor/antisymmetric_tensor.hpp>
 #include <similie/tensor/dummy_index.hpp>
 #include <similie/tensor/tensor_impl.hpp>
@@ -20,6 +21,7 @@
 
 #include "cochain.hpp"
 #include "cosimplex.hpp"
+#include "form.hpp"
 
 
 namespace sil {
@@ -334,6 +336,75 @@ coboundary_tensor_t<TagToAddToCochain, CochainTag, TensorType> deriv(
             TensorType>(exec_space, coboundary_tensor, tensor);
 }
 
+template <
+        class TagToAddToCochain,
+        tensor::TensorIndex CochainTag,
+        misc::Specialization<tensor::Tensor> OutTensorType,
+        misc::Specialization<tensor::Tensor> 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<typename OutTensorType::indices_domain_t>>;
+    using in_index_t = ddc::type_seq_element_t<
+            0,
+            ddc::to_type_seq_t<typename TensorType::indices_domain_t>>;
+    using in_d_dim_t = mesher::detail::
+            discrete_dimension_for_t<TagToAddToCochain, typename TensorType::non_indices_domain_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<in_d_dim_t>);
+
+    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<in_d_dim_t>(1);
+                double const dx = static_cast<double>(
+                        ddc::coordinate(ddc::DiscreteElement<in_d_dim_t>(in_right))
+                        - ddc::coordinate(ddc::DiscreteElement<in_d_dim_t>(in_left)));
+                out_tensor.mem(out_elem, ddc::DiscreteElement<out_index_t>(0))
+                        = (tensor.get(in_right, ddc::DiscreteElement<in_index_t>(0))
+                           - tensor.get(in_left, ddc::DiscreteElement<in_index_t>(0)))
+                          / dx;
+            });
+    return out_tensor;
+}
+
+template <
+        class SupportTag,
+        class... Components,
+        misc::Specialization<tensor::Tensor> TensorType,
+        class ExecSpace>
+TensorForm<SupportTag, Components...> deriv(
+        ExecSpace const& exec_space,
+        TensorForm<SupportTag, Components...> out_form,
+        TensorType tensor)
+{
+    using in_index_t = ddc::type_seq_element_t<
+            0,
+            ddc::to_type_seq_t<typename TensorType::indices_domain_t>>;
+    static_assert(in_index_t::rank() == 0);
+
+    (...,
+     sil::exterior::deriv<typename Components::tag, in_index_t>(
+             exec_space,
+             out_form.template component<typename Components::tag>(),
+             tensor,
+             mesher::HalfShiftDualizer<typename Components::tag> {}));
+    return out_form;
+}
+
 } // namespace exterior
 
 } // namespace sil

src/similie/exterior/codifferential.hpp

@@ -6,6 +6,7 @@
 #include <ddc/ddc.hpp>
 
 #include <similie/exterior/hodge_star.hpp>
+#include <similie/mesher/dualizer.hpp>
 #include <similie/misc/macros.hpp>
 #include <similie/misc/specialization.hpp>
 #include <similie/tensor/character.hpp>
@@ -15,6 +16,7 @@
 
 #include "cochain.hpp"
 #include "cosimplex.hpp"
+#include "form.hpp"
 
 
 namespace sil {
@@ -156,6 +158,125 @@ struct CodifferentialDummyIndexSeq<EndId, T>
     using type = typename CodifferentialDummyIndexSeq_<std::make_index_sequence<EndId>, T>::type;
 };
 
+template <class DDim>
+constexpr bool is_dual_discrete_dimension_v
+        = !std::is_same_v<DDim, mesher::detail::primal_discrete_dimension_t<DDim>>;
+
+template <class DDim, class Domain, class Elem>
+KOKKOS_FUNCTION bool is_on_primal_boundary(Domain const& domain, Elem const& elem)
+{
+    if constexpr (is_dual_discrete_dimension_v<DDim>) {
+        return false;
+    } else {
+        ddc::DiscreteDomain<DDim> const dim_dom(domain);
+        ddc::DiscreteElement<DDim> 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<tensor::Tensor> MetricType,
+        class Elem>
+KOKKOS_FUNCTION double codifferential_metric_factor(MetricType const& inv_metric, Elem const& elem)
+{
+    if constexpr (misc::Specialization<MetricComponentIndex, tensor::TensorIdentityIndex>) {
+        return 1.;
+    } else {
+        return inv_metric(
+                elem,
+                inv_metric.accessor().template access_element<
+                        TagToRemoveFromCochain,
+                        TagToRemoveFromCochain>());
+    }
+}
+
+template <
+        class AxisTag,
+        misc::Specialization<tensor::Tensor> OutTensorType,
+        misc::Specialization<tensor::Tensor> 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<AxisTag, typename OutTensorType::non_indices_domain_t>;
+    using comp_d_dim_t = mesher::detail::
+            discrete_dimension_for_t<AxisTag, typename TensorType::non_indices_domain_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<typename TensorType::indices_domain_t>>;
+    static_assert(!std::is_void_v<out_d_dim_t>);
+    static_assert(!std::is_void_v<comp_d_dim_t>);
+    static_assert(comp_index_t::rank() == 0);
+
+    comp_elem_t const aligned_elem = mesher::detail::dualizer_remap_element<comp_elem_t>(out_elem);
+    if constexpr (is_dual_discrete_dimension_v<out_d_dim_t>) {
+        auto const left = aligned_elem;
+        auto const right = left + ddc::DiscreteVector<comp_d_dim_t>(1);
+        double const dx = static_cast<double>(
+                ddc::coordinate(ddc::DiscreteElement<comp_d_dim_t>(right))
+                - ddc::coordinate(ddc::DiscreteElement<comp_d_dim_t>(left)));
+        return (tensor.get(right, ddc::DiscreteElement<comp_index_t>(0))
+                - tensor.get(left, ddc::DiscreteElement<comp_index_t>(0)))
+               / dx;
+    } else {
+        auto const right = aligned_elem;
+        auto const left = right - ddc::DiscreteVector<comp_d_dim_t>(1);
+        double const dx = static_cast<double>(
+                ddc::coordinate(ddc::DiscreteElement<comp_d_dim_t>(right))
+                - ddc::coordinate(ddc::DiscreteElement<comp_d_dim_t>(left)));
+        return (tensor.get(right, ddc::DiscreteElement<comp_index_t>(0))
+                - tensor.get(left, ddc::DiscreteElement<comp_index_t>(0)))
+               / dx;
+    }
+}
+
+template <
+        misc::Specialization<tensor::Tensor> OutTensorType,
+        class SupportTag,
+        class FirstComponent,
+        class SecondComponent,
+        class ExecSpace>
+OutTensorType exterior_derivative_of_tensor_form_2d(
+        ExecSpace const& exec_space,
+        OutTensorType out_tensor,
+        TensorForm<SupportTag, FirstComponent, SecondComponent> tensor_form)
+{
+    using out_index_t = ddc::type_seq_element_t<
+            0,
+            ddc::to_type_seq_t<typename OutTensorType::indices_domain_t>>;
+    using first_out_d_dim_t = mesher::detail::
+            discrete_dimension_for_t<typename FirstComponent::tag, typename OutTensorType::non_indices_domain_t>;
+    using second_out_d_dim_t = mesher::detail::
+            discrete_dimension_for_t<typename SecondComponent::tag, typename OutTensorType::non_indices_domain_t>;
+    static_assert(out_index_t::rank() == 0);
+
+    auto const first_tensor = tensor_form.template component<typename FirstComponent::tag>();
+    auto const second_tensor = tensor_form.template component<typename SecondComponent::tag>();
+
+    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 (is_on_primal_boundary<first_out_d_dim_t>(out_tensor.non_indices_domain(), elem)
+                    || is_on_primal_boundary<second_out_d_dim_t>(out_tensor.non_indices_domain(), elem)) {
+                    return;
+                }
+                out_tensor.mem(elem, ddc::DiscreteElement<out_index_t>(0))
+                        += centered_form_derivative<typename FirstComponent::tag, OutTensorType>(
+                                   elem,
+                                   second_tensor)
+                           - centered_form_derivative<typename SecondComponent::tag, OutTensorType>(
+                                   elem,
+                                   first_tensor);
+            });
+    return out_tensor;
+}
+
 } // namespace detail
 
 template <
@@ -279,6 +400,101 @@ codifferential_tensor_t<TagToRemoveFromCochain, CochainTag, TensorType> codiffer
     return codifferential_tensor;
 }
 
+template <
+        tensor::TensorIndex MetricIndex,
+        class TagToRemoveFromCochain,
+        tensor::TensorIndex CochainTag,
+        misc::Specialization<tensor::Tensor> OutTensorType,
+        misc::Specialization<tensor::Tensor> TensorType,
+        misc::Specialization<tensor::Tensor> 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<typename OutTensorType::indices_domain_t>>;
+    using metric_component_index_t = ddc::type_seq_element_t<
+            0,
+            ddc::to_type_seq_t<typename MetricType::indices_domain_t>>;
+    using out_d_dim_t = mesher::detail::
+            discrete_dimension_for_t<TagToRemoveFromCochain, typename OutTensorType::non_indices_domain_t>;
+    using flux_d_dim_t = mesher::detail::
+            discrete_dimension_for_t<TagToRemoveFromCochain, typename TensorType::non_indices_domain_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<typename TensorType::indices_domain_t>>;
+    static_assert(out_index_t::rank() == 0);
+    static_assert(!std::is_void_v<out_d_dim_t>);
+    static_assert(!std::is_void_v<flux_d_dim_t>);
+
+    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<out_d_dim_t> dim_dom(out_tensor.non_indices_domain());
+                ddc::DiscreteElement<out_d_dim_t> 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<flux_d_dim_t>(1);
+                double const dx = static_cast<double>(
+                        ddc::coordinate(ddc::DiscreteElement<flux_d_dim_t>(right_face))
+                        - ddc::coordinate(ddc::DiscreteElement<flux_d_dim_t>(left_face)));
+                double const deriv = (tensor.get(right_face, ddc::DiscreteElement<flux_index_t>(0))
+                                      - tensor.get(left_face, ddc::DiscreteElement<flux_index_t>(0)))
+                                     / dx;
+                double const metric_factor = detail::codifferential_metric_factor<
+                        metric_component_index_t,
+                        TagToRemoveFromCochain>(inv_metric, out_elem);
+                out_tensor.mem(out_elem, ddc::DiscreteElement<out_index_t>(0))
+                        += metric_factor * deriv;
+            });
+    return out_tensor;
+}
+
+template <
+        tensor::TensorIndex MetricIndex,
+        misc::Specialization<tensor::Tensor> OutTensorType,
+        class SupportTag,
+        class FirstComponent,
+        class SecondComponent,
+        misc::Specialization<tensor::Tensor> MetricType,
+        class ExecSpace>
+OutTensorType codifferential(
+        ExecSpace const& exec_space,
+        OutTensorType out_tensor,
+        TensorForm<SupportTag, FirstComponent, SecondComponent> tensor_form,
+        MetricType inv_metric)
+{
+    ddc::Chunk dual_first_alloc(
+            typename SecondComponent::tensor_type::discrete_domain_type(
+                    tensor_form.template component<typename SecondComponent::tag>().domain()),
+            ddc::KokkosAllocator<double, typename ExecSpace::memory_space>());
+    sil::tensor::Tensor dual_first(dual_first_alloc);
+    ddc::Chunk dual_second_alloc(
+            typename FirstComponent::tensor_type::discrete_domain_type(
+                    tensor_form.template component<typename FirstComponent::tag>().domain()),
+            ddc::KokkosAllocator<double, typename ExecSpace::memory_space>());
+    sil::tensor::Tensor dual_second(dual_second_alloc);
+    auto dual_form = make_tensor_form<hodge_dual_support_t<SupportTag>>(
+            component<typename FirstComponent::tag>(dual_first),
+            component<typename SecondComponent::tag>(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
 
 } // namespace sil