utilities.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan
// SPDX-License-Identifier: GNU Lesser General Public Version 2.1
 
#pragma once
 
#include <deal.II/base/config.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/mpi.h>
#include <deal.II/base/point.h>
#include <deal.II/base/tensor.h>
#include <deal.II/base/types.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/vectorization.h>
#include <deal.II/matrix_free/evaluation_flags.h>
 
#include <boost/range/algorithm.hpp>
#include <boost/range/algorithm_ext.hpp>
 
#include <prismspf/utilities/vectorized_operations.h>
 
#include <prismspf/config.h>
 
#include <string>
#include <string_view>
#include <unordered_map>
 
#ifdef PRISMS_PF_WITH_CALIPER
#  include <caliper/cali-manager.h>
#  include <caliper/cali.h>
#endif
 
PRISMS_PF_BEGIN_NAMESPACE

class MPIInitFinalize : public dealii::Utilities::MPI::MPI_InitFinalize
{
public:
  MPIInitFinalize(int &argc, char **&argv, unsigned int _max_n_threads = 0);
 
  ~MPIInitFinalize();
 
private:
#ifdef PRISMS_PF_THREADS
  static constexpr unsigned int max_n_threads = dealii::numbers::invalid_unsigned_int;
#else
  static constexpr unsigned int max_n_threads = 1;
#endif

#ifdef PRISMS_PF_WITH_CALIPER
  cali::ConfigManager mgr;
#endif
};
 
template <typename dealii_vector>
void
set_zero_entries_to_one(dealii_vector &vec)
{
  for (const auto index : vec.locally_owned_elements())
    {
      if (vec[index] == 0.0)
        {
          vec[index] = 1.0;
        }
    }
}

template <typename Type>
Type
string_to_type(const std::string                           &string,
               const std::unordered_map<std::string, Type> &table)
{
  auto iterator = table.find(string);
  AssertThrow(iterator != table.end(),
              dealii::ExcMessage("Unknown table entry: " + string));
  return iterator->second;
}

template <typename Type, typename OtherType>
std::pair<OtherType, Type>
string_to_type_pair_with_delimiters(
  const std::string                                &string,
  const std::unordered_map<std::string, Type>      &table,
  const std::unordered_map<std::string, OtherType> &other_table,
  const std::pair<char, char>                      &delimiters = {'(', ')'})
{
  // Check to see if the string has any delimiters
  const auto opening_delimiter_position = string.find(delimiters.first);
  const auto closing_delimiter_position = string.find(delimiters.second);
 
  bool has_opening_delimiter = opening_delimiter_position != std::string_view::npos;
  bool has_closing_delimiter = closing_delimiter_position != std::string_view::npos;
 
  // Some checks for malformed input
  AssertThrow(closing_delimiter_position >= opening_delimiter_position,
              dealii::ExcMessage(
                "Opening delimiter must precede the closing delimiter. You had " +
                string));
  AssertThrow(has_closing_delimiter == has_opening_delimiter,
              dealii::ExcMessage("You must either have a delimiter pair or not. You seem "
                                 "to only have one delimiter. You had " +
                                 string));
 
  // Case 2 - no delimiters
  if (!has_closing_delimiter && !has_opening_delimiter)
    {
      const auto second_type = string_to_type("", table);
      const auto first_type  = string_to_type(string, other_table);
      return {first_type, second_type};
    }
 
  // More checks for malformed input
  AssertThrow(closing_delimiter_position == string.size() - 1,
              dealii::ExcMessage(
                "Closing delimiter must be at the end of the string. You had " + string));
 
  // Case 1 - delimiters
  std::string_view input_view {string};
  std::string_view outer_string = input_view.substr(0, opening_delimiter_position);
  std::string_view inner_string =
    input_view.substr(opening_delimiter_position + 1,
                      closing_delimiter_position - opening_delimiter_position - 1);
 
  AssertThrow(!outer_string.empty() && !inner_string.empty(),
              dealii::ExcMessage(
                "Inner and outer map entries must not be empty. You had " + string));
  // NOTE: This will allocate the string views to strings on the heap. Hopefully this
  // isn't so bad in performance, but I didn't wanna support heterogeneous lookup tables.
  // I can't imagine this becoming a performance problem right now since we only deal with
  // strings in the initialization of simulations; however, if it does, it shouldn't be
  // two hard to support heterogeneous lookup or a table that uses std::string_view.
  const auto second_type = string_to_type(std::string(outer_string), table);
  const auto first_type  = string_to_type(std::string(inner_string), other_table);
  return {first_type, second_type};
}

template <typename Number>
inline DEAL_II_ALWAYS_INLINE Number
pmod(const Number &value, const Number &modulus)
{
  return std::fmod(std::fmod(value, modulus) + modulus, modulus);
}

template <unsigned int dim, unsigned int degree, typename T>
inline DEAL_II_ALWAYS_INLINE T
stabilization_parameter(const T                         &timestep,
                        const T                         &element_volume,
                        const dealii::Tensor<1, dim, T> &velocity,
                        const T                         &kinematic_viscosity)
{
  using dealii::Utilities::fixed_power;
  using std::cbrt;
  using std::sqrt;
 
  constexpr auto degree_modifier  = T(1) / T(degree);
  constexpr auto inv_pi           = T(1) / T(M_PI);
  const auto     velocity_l2_norm = velocity.norm_square() + T(1e-12);
 
  auto element_size = T(0);
  if constexpr (dim == 1)
    {
      element_size = element_volume;
    }
  else if constexpr (dim == 2)
    {
      element_size = T(2) * sqrt(element_volume * inv_pi);
    }
  else if constexpr (dim == 3)
    {
      element_size = cbrt(T(6) * element_volume * inv_pi);
    }
  else
    {
      Assert(false, dealii::ExcMessage("Invalid dimension for the stabilization term."));
    }
  element_size *= degree_modifier;
  const auto element_size_2     = element_size * element_size;
  const auto element_size_2_inv = T(1) / element_size_2;
  const auto element_size_4_inv = element_size_2_inv * element_size_2_inv;
 
  const T term_time = fixed_power<2>(T(2) / timestep);
  const T term_adv  = T(4) * velocity_l2_norm * element_size_2_inv;
  const T term_diff = T(144) * fixed_power<2>(kinematic_viscosity) * element_size_4_inv;
 
  return T(1) / sqrt(term_time + term_adv + term_diff);
}

template <unsigned int dim, unsigned int degree, typename T>
inline DEAL_II_ALWAYS_INLINE T
stabilization_parameter(const T                         &element_volume,
                        const dealii::Tensor<1, dim, T> &velocity,
                        const T                         &kinematic_viscosity)
{
  using dealii::Utilities::fixed_power;
  using std::cbrt;
  using std::sqrt;
 
  constexpr auto degree_modifier  = T(1) / T(degree);
  constexpr auto inv_pi           = T(1) / T(M_PI);
  const auto     velocity_l2_norm = velocity.norm_square() + T(1e-12);
 
  auto element_size = T(0);
  if constexpr (dim == 1)
    {
      element_size = element_volume;
    }
  else if constexpr (dim == 2)
    {
      element_size = T(2) * sqrt(element_volume * inv_pi);
    }
  else if constexpr (dim == 3)
    {
      element_size = cbrt(T(6) * element_volume * inv_pi);
    }
  else
    {
      Assert(false, dealii::ExcMessage("Invalid dimension for the stabilization term."));
    }
  element_size *= degree_modifier;
  const auto element_size_2     = element_size * element_size;
  const auto element_size_2_inv = T(1) / element_size_2;
  const auto element_size_4_inv = element_size_2_inv * element_size_2_inv;
 
  const T term_adv  = T(4) * velocity_l2_norm * element_size_2_inv;
  const T term_diff = T(144) * fixed_power<2>(kinematic_viscosity) * element_size_4_inv;
 
  return T(1) / sqrt(term_adv + term_diff);
}

template <unsigned int dim, unsigned int degree, typename T>
inline DEAL_II_ALWAYS_INLINE T
stabilization_parameter(const T                         &timestep,
                        const T                         &element_volume,
                        const dealii::Tensor<1, dim, T> &velocity)
{
  using dealii::Utilities::fixed_power;
  using std::cbrt;
  using std::sqrt;
 
  constexpr auto degree_modifier  = T(1) / T(degree);
  constexpr auto inv_pi           = T(1) / T(M_PI);
  const auto     velocity_l2_norm = velocity.norm_square() + T(1e-12);
 
  auto element_size = T(0);
  if constexpr (dim == 1)
    {
      element_size = element_volume;
    }
  else if constexpr (dim == 2)
    {
      element_size = T(2) * sqrt(element_volume * inv_pi);
    }
  else if constexpr (dim == 3)
    {
      element_size = cbrt(T(6) * element_volume * inv_pi);
    }
  else
    {
      Assert(false, dealii::ExcMessage("Invalid dimension for the stabilization term."));
    }
  element_size *= degree_modifier;
  const auto element_size_2     = element_size * element_size;
  const auto element_size_2_inv = T(1) / element_size_2;
  const auto element_size_4_inv = element_size_2_inv * element_size_2_inv;
 
  const T term_time = fixed_power<2>(T(2) / timestep);
  const T term_adv  = T(4) * velocity_l2_norm * element_size_2_inv;
 
  return T(1) / sqrt(term_time + term_adv);
}

template <unsigned int dim, unsigned int degree, typename T>
inline DEAL_II_ALWAYS_INLINE T
stabilization_parameter(const T                         &element_volume,
                        const dealii::Tensor<1, dim, T> &velocity)
{
  using dealii::Utilities::fixed_power;
  using std::cbrt;
  using std::sqrt;
 
  constexpr auto degree_modifier  = T(1) / T(degree);
  constexpr auto inv_pi           = T(1) / T(M_PI);
  const auto     velocity_l2_norm = velocity.norm_square() + T(1e-12);
 
  auto element_size = T(0);
  if constexpr (dim == 1)
    {
      element_size = element_volume;
    }
  else if constexpr (dim == 2)
    {
      element_size = T(2) * sqrt(element_volume * inv_pi);
    }
  else if constexpr (dim == 3)
    {
      element_size = cbrt(T(6) * element_volume * inv_pi);
    }
  else
    {
      Assert(false, dealii::ExcMessage("Invalid dimension for the stabilization term."));
    }
  element_size *= degree_modifier;
  const auto element_size_2     = element_size * element_size;
  const auto element_size_2_inv = T(1) / element_size_2;
  const auto element_size_4_inv = element_size_2_inv * element_size_2_inv;
 
  const T term_adv = T(4) * velocity_l2_norm * element_size_2_inv;
  // NOTE: I could simplify the math here, but I won't bother for now since I imagine the
  // use case is pretty limited.
 
  return T(1) / sqrt(term_adv);
}

inline std::string
strip_whitespace(const std::string &_text)
{
  std::string text = _text;
  boost::range::remove_erase_if(text, ::isspace);
  return text;
}

inline const char *
bool_to_string(bool boolean)
{
  return boolean ? "true" : "false";
}

inline std::string
eval_flags_to_string(dealii::EvaluationFlags::EvaluationFlags flag)
{
  std::string result;
 
  if ((flag & dealii::EvaluationFlags::EvaluationFlags::values) != 0U)
    {
      result += "values";
    }
  if ((flag & dealii::EvaluationFlags::EvaluationFlags::gradients) != 0U)
    {
      if (!result.empty())
        {
          result += " | ";
        }
      result += "gradients";
    }
  if ((flag & dealii::EvaluationFlags::EvaluationFlags::hessians) != 0U)
    {
      if (!result.empty())
        {
          result += " | ";
        }
      result += "hessians";
    }
 
  if (result.empty())
    {
      return "nothing";
    }
 
  return result;
}
 
template <unsigned int dim, typename number>
inline DEAL_II_ALWAYS_INLINE std::vector<number>
dealii_point_to_vector(const dealii::Point<dim, number> &point)
{
  static_assert(dim < 4, "We only allow 3 space dimensions");
 
  std::vector<number> vec(3, 0.0);
 
  // NOLINTBEGIN(cppcoreguidelines-pro-bounds-constant-array-index)
 
  for (unsigned int i = 0; i < dim; ++i)
    {
      vec[i] = point[i];
    }
 
  // NOLINTEND(cppcoreguidelines-pro-bounds-constant-array-index)
 
  return vec;
}
 
PRISMS_PF_END_NAMESPACE
 
#include <prismspf/utilities/vectorized_operations.h>