state_vector.h

Go to the documentation of this file.

//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
// Copyright (C) 2024 - 2026 by the ryujin authors
//
 
#pragma once
 
#include <compile_time_options.h>
 
#include "computing_timer.h"
#include "multicomponent_vector.h"
 
#include <algorithm>
#include <concepts>
#include <limits>
#include <vector>
 
namespace ryujin
{
#ifndef DOXYGEN
  /* Forward declaration */
  template <int dim, typename Number>
  class OfflineData;
#endif
 
  template <typename Number>
  using ScalarHostVector = dealii::LinearAlgebra::distributed::Vector<Number>;
 
 
  template <typename Number>
  using ScalarVector = MultiComponentVector<Number>;
 
 
  template <typename Number>
  struct StateVector {
    MultiComponentVector<Number> U;
    MultiComponentVector<Number> precomputed;
    std::vector<ScalarVector<Number>> V;
  };
 
 
  enum class StateVectorPart {
    hyperbolic,
    precomputed,
    parabolic,
  };
 
 
  template <typename MemorySpace,
            typename Number,
            std::same_as<StateVectorPart>... Parts>
  void copy_to_memory_space(const StateVector<Number> &state_vector
                            [[maybe_unused]],
                            Parts... parts [[maybe_unused]])
  {
    static_assert(sizeof...(Parts) > 0, "No state vector part selected");
 
    if constexpr (have_separate_memory_spaces) {
      ComputingTimer::Scope scope("time step [X] _ - memory space transfers");
 
      const auto &[U, precomputed, parabolic] = state_vector;
 
      const auto copy = [&](const StateVectorPart part) {
        switch (part) {
        case StateVectorPart::hyperbolic:
          U.template copy_to_memory_space<MemorySpace>();
          break;
        case StateVectorPart::precomputed:
          precomputed.template copy_to_memory_space<MemorySpace>();
          break;
        case StateVectorPart::parabolic:
          for (const auto &V : parabolic)
            V.template copy_to_memory_space<MemorySpace>();
          break;
        }
      };
 
      (copy(parts), ...);
    }
  }
 
 
  template <typename MemorySpace, typename Number>
  void copy_to_memory_space(const MultiComponentVector<Number> &vector
                            [[maybe_unused]])
  {
    if constexpr (have_separate_memory_spaces) {
      ComputingTimer::Scope scope("time step [X] _ - memory space transfers");
      vector.template copy_to_memory_space<MemorySpace>();
    }
  }
 
 
  template <typename MemorySpace,
            typename Number,
            std::same_as<StateVectorPart>... Parts>
  void move_to_memory_space(StateVector<Number> &state_vector [[maybe_unused]],
                            Parts... parts [[maybe_unused]])
  {
    static_assert(sizeof...(Parts) > 0, "No state vector part selected");
 
    if constexpr (have_separate_memory_spaces) {
      ComputingTimer::Scope scope("time step [X] _ - memory space transfers");
 
      auto &[U, precomputed, parabolic] = state_vector;
 
      const auto move = [&](const StateVectorPart part) {
        switch (part) {
        case StateVectorPart::hyperbolic:
          U.template move_to_memory_space<MemorySpace>();
          break;
        case StateVectorPart::precomputed:
          precomputed.template move_to_memory_space<MemorySpace>();
          break;
        case StateVectorPart::parabolic:
          for (auto &V : parabolic)
            V.template move_to_memory_space<MemorySpace>();
          break;
        }
      };
 
      (move(parts), ...);
    }
  }
 
 
  template <typename MemorySpace, typename Number>
  void move_to_memory_space(MultiComponentVector<Number> &vector
                            [[maybe_unused]])
  {
    if constexpr (have_separate_memory_spaces) {
      ComputingTimer::Scope scope("time step [X] _ - memory space transfers");
      vector.template move_to_memory_space<MemorySpace>();
    }
  }
 
 
  template <typename MemorySpace, typename Number>
  void sadd(StateVector<Number> &dst,
            const Number s,
            const Number b,
            const StateVector<Number> &src)
  {
    copy_to_memory_space<MemorySpace>(
        src, StateVectorPart::hyperbolic, StateVectorPart::parabolic);
    move_to_memory_space<MemorySpace>(
        dst, StateVectorPart::hyperbolic, StateVectorPart::parabolic);
 
    auto &dst_U = dst.U.template deal_ii_vector<MemorySpace>();
    const auto &src_U = src.U.template deal_ii_vector<MemorySpace>();
    dst_U.zero_out_ghost_values();
    dst_U.sadd(s, b, src_U);
 
    AssertDimension(dst.V.size(), src.V.size());
 
    for (std::size_t k = 0; k < dst.V.size(); ++k) {
      auto &dst_V_k = dst.V[k].template deal_ii_vector<MemorySpace>();
      const auto &src_V_k = src.V[k].template deal_ii_vector<MemorySpace>();
      dst_V_k.zero_out_ghost_values();
      dst_V_k.sadd(s, b, src_V_k);
    }
  }
 
 
  template <typename Number, typename OfflineData>
  void debug_poison_invalid_values(StateVector<Number> &state_vector
                                   [[maybe_unused]],
                                   const OfflineData &offline_data
                                   [[maybe_unused]])
  {
#ifdef DEBUG
 
    move_to_memory_space<dealii::MemorySpace::Host>(
        state_vector,
        StateVectorPart::hyperbolic,
        StateVectorPart::precomputed);
 
    auto &[U, prec, V] = state_vector;
 
    constexpr auto nan = std::numeric_limits<Number>::signaling_NaN();
 
    const unsigned int n_owned = offline_data.n_locally_owned();
    const unsigned int n_relevant = offline_data.n_locally_relevant();
    const auto &partitioner = offline_data.scalar_partitioner();
 
    /*
     * The number of components is a run time property: work on the
     * underlying deal.II vectors with n_comp() as stride.
     */
    const unsigned int prob_dim = U.n_comp();
    const unsigned int prec_dim = prec.n_comp();
    Number *const U_data = U.deal_ii_vector().begin();
    Number *const prec_data = prec.deal_ii_vector().begin();
 
    for (unsigned int i = 0; i < n_owned; ++i) {
      std::fill_n(prec_data + i * prec_dim, prec_dim, nan);
 
      if (!offline_data.affine_constraints().is_constrained(
              partitioner->local_to_global(i)))
        continue;
      std::fill_n(U_data + i * prob_dim, prob_dim, nan);
    }
 
    for (unsigned int i = n_owned; i < n_relevant; ++i) {
      std::fill_n(prec_data + i * prec_dim, prec_dim, nan);
 
      std::fill_n(U_data + i * prob_dim, prob_dim, nan);
    }
#endif
  }
} // namespace ryujin