BasisFunctions.hpp

Go to the documentation of this file.

#pragma once
 
#include <stdexcept>
 
#ifdef __CUDACC__
#include <cuda_runtime.h>
#endif
 
#ifdef __HIPCC__
#include <hip/hip_runtime.h>
#endif
 
#if defined(__CUDACC__) || defined(__HIPCC__)
 
#include <ginkgo/extensions/kokkos.hpp>
#include <ginkgo/ginkgo.hpp>
 
#define PRECICE_HOST_DEVICE __host__ __device__
#define PRECICE_MEMORY_SPACE __device__
#define PRECICE_FMA Kokkos::fma
#define PRECICE_LOG Kokkos::log
 
#else
 
#define PRECICE_HOST_DEVICE
#define PRECICE_MEMORY_SPACE
#define PRECICE_FMA std::fma
#define PRECICE_LOG std::log
 
#endif
 
#define NUMERICAL_ZERO_DIFFERENCE_DEVICE 1.0e-14
 
#include "math/math.hpp"
 
namespace precice::mapping {

struct RadialBasisParameters {
  double parameter1{};
  double parameter2{};
  double parameter3{};
};

struct CompactSupportBase {
  static constexpr bool hasCompactSupport()
  {
    return true;
  }
};

struct NoCompactSupportBase {
  static constexpr bool hasCompactSupport()
  {
    return false;
  }
 
  static constexpr double getSupportRadius()
  {
    return std::numeric_limits<double>::max();
  }
};

template <bool isDefinite>
struct DefiniteFunction {
  static constexpr bool isStrictlyPositiveDefinite()
  {
    return isDefinite;
  }
};

class ThinPlateSplines : public NoCompactSupportBase,
                         public DefiniteFunction<false> {
public:
  double evaluate(double radius) const
  {
    return operator()(radius, _params);
  }
 
  PRECICE_HOST_DEVICE inline double operator()(const double radius, const RadialBasisParameters params) const
  {
    // We don't need to read any values from params since there is no need here
    return PRECICE_LOG(std::max(radius, NUMERICAL_ZERO_DIFFERENCE_DEVICE)) * math::pow_int<2>(radius);
  }
 
  RadialBasisParameters getFunctionParameters()
  {
    return _params;
  }
 
private:
  RadialBasisParameters _params;
};

class Multiquadrics : public NoCompactSupportBase,
                      public DefiniteFunction<false> {
public:
  explicit Multiquadrics(double c)
      : _cPow2(math::pow_int<2>(c))
  {
    _params.parameter1 = _cPow2;
  }
 
  double evaluate(double radius) const
  {
    return operator()(radius, _params);
  }
 
  PRECICE_HOST_DEVICE inline double operator()(const double radius, const RadialBasisParameters params) const
  {
    double cPow2 = params.parameter1;
    return std::sqrt(cPow2 + math::pow_int<2>(radius));
  }
 
  RadialBasisParameters getFunctionParameters()
  {
    return _params;
  }
 
private:
  double                _cPow2;
  RadialBasisParameters _params;
};

class InverseMultiquadrics : public NoCompactSupportBase,
                             public DefiniteFunction<true> {
public:
  explicit InverseMultiquadrics(double c)
      : _cPow2(math::pow_int<2>(c))
  {
    if (!math::greater(c, 0.0)) {
      throw std::invalid_argument{"Shape parameter for radial-basis-function inverse multiquadric has to be larger than zero. Please update the \"shape-parameter\" attribute."};
    }
    _params.parameter1 = _cPow2;
  }
 
  double evaluate(double radius) const
  {
    return operator()(radius, _params);
  }
 
  PRECICE_HOST_DEVICE inline double operator()(const double radius, const RadialBasisParameters params) const
  {
    double cPow2 = params.parameter1;
    return 1.0 / std::sqrt(cPow2 + math::pow_int<2>(radius));
  }
 
  RadialBasisParameters getFunctionParameters()
  {
    return _params;
  }
 
private:
  double const _cPow2;
 
  RadialBasisParameters _params;
};

class VolumeSplines : public NoCompactSupportBase,
                      public DefiniteFunction<false> {
public:
  double evaluate(double radius) const
  {
    return operator()(radius, _params);
  }
 
  PRECICE_HOST_DEVICE inline double operator()(const double radius, const RadialBasisParameters params) const
  {
    return std::abs(radius);
  }
 
  RadialBasisParameters getFunctionParameters()
  {
    return _params;
  }
 
private:
  RadialBasisParameters _params;
};

class Gaussian : public CompactSupportBase,
                 public DefiniteFunction<true> {
public:
  Gaussian(const double shape, const double supportRadius = std::numeric_limits<double>::infinity())
      : _shape(shape),
        _supportRadius(supportRadius)
  {
    if (math::greater(_shape, 0.0)) {
      std::invalid_argument{
          "Shape parameter for radial-basis-function gaussian has to be larger than zero. Please update the \"shape-parameter\" attribute."};
    }
    if (math::greater(_supportRadius, 0.0)) {
      std::invalid_argument{
          "Support radius for radial-basis-function gaussian has to be larger than zero. Please update the \"support-radius\" attribute."};
    }
 
    double threshold   = std::sqrt(-std::log(cutoffThreshold)) / shape;
    _supportRadius     = std::min(supportRadius, threshold);
    _params.parameter1 = _shape;
    _params.parameter2 = _supportRadius;
    _params.parameter3 = _deltaY;
    if (supportRadius < std::numeric_limits<double>::infinity()) {
      _deltaY            = evaluate(supportRadius);
      _params.parameter3 = _deltaY;
    }
  }
 
  double getSupportRadius() const
  {
    return _supportRadius;
  }
 
  double evaluate(const double radius) const
  {
    return operator()(radius, _params);
  }
 
  PRECICE_HOST_DEVICE inline double operator()(const double radius, const RadialBasisParameters params) const
  {
    double shape         = params.parameter1;
    double supportRadius = params.parameter2;
    double deltaY        = params.parameter3;
 
    if (radius > supportRadius)
      return 0.0;
    return std::exp(-math::pow_int<2>(shape * radius)) - deltaY;
  }
 
  RadialBasisParameters getFunctionParameters()
  {
    return _params;
  };
 
public:
  static constexpr double cutoffThreshold = 1e-9;
 
private:
  double const _shape;

  double _supportRadius;
 
  double _deltaY = 0;
 
  RadialBasisParameters _params;
};

class CompactThinPlateSplinesC2 : public CompactSupportBase,
                                  public DefiniteFunction<true> {
public:
  explicit CompactThinPlateSplinesC2(double supportRadius)
  {
    if (math::greater(supportRadius, 0.0)) {
      std::invalid_argument{
          "Support radius for radial-basis-function compact thin-plate-splines c2 has to be larger than zero. Please update the \"support-radius\" attribute."};
    }
    _r_inv             = 1. / supportRadius;
    _params.parameter1 = _r_inv;
  }
 
  double getSupportRadius() const
  {
    return 1. / _r_inv;
  }
 
  double evaluate(double radius) const
  {
    return operator()(radius, _params);
  }
 
  PRECICE_HOST_DEVICE inline double operator()(const double radius, const RadialBasisParameters params) const
  {
    double       r_inv = params.parameter1;
    const double p     = radius * r_inv;
    if (p >= 1)
      return 0.0;
    return 1.0 - 30.0 * math::pow_int<2>(p) - 10.0 * math::pow_int<3>(p) + 45.0 * math::pow_int<4>(p) - 6.0 * math::pow_int<5>(p) - math::pow_int<3>(p) * 60.0 * PRECICE_LOG(std::max(p, NUMERICAL_ZERO_DIFFERENCE_DEVICE));
  }
 
  RadialBasisParameters getFunctionParameters()
  {
    return _params;
  }
 
private:
  double                _r_inv;
  RadialBasisParameters _params;
};

class CompactPolynomialC0 : public CompactSupportBase,
                            public DefiniteFunction<true> {
public:
  explicit CompactPolynomialC0(double supportRadius)
  {
    if (math::greater(supportRadius, 0.0)) {
      std::invalid_argument{
          "Support radius for radial-basis-function compact polynomial c0 has to be larger than zero. Please update the \"support-radius\" attribute."};
    }
    _r_inv             = 1. / supportRadius;
    _params.parameter1 = _r_inv;
  }
 
  double getSupportRadius() const
  {
    return 1. / _r_inv;
  }
 
  double evaluate(double radius) const
  {
    return operator()(radius, _params);
  }
 
  PRECICE_HOST_DEVICE inline double operator()(const double radius, const RadialBasisParameters params) const
  {
    double       r_inv = params.parameter1;
    const double p     = radius * r_inv;
    if (p >= 1)
      return 0.0;
    return math::pow_int<2>(1.0 - p);
  }
 
  RadialBasisParameters getFunctionParameters()
  {
    return _params;
  }
 
private:
  double                _r_inv;
  RadialBasisParameters _params;
};

class CompactPolynomialC2 : public CompactSupportBase,
                            public DefiniteFunction<true> {
public:
  explicit CompactPolynomialC2(double supportRadius)
  {
    if (math::greater(supportRadius, 0.0)) {
      std::invalid_argument{
          "Support radius for radial-basis-function compact polynomial c2 has to be larger than zero. Please update the \"support-radius\" attribute."};
    }
 
    _r_inv             = 1. / supportRadius;
    _params.parameter1 = _r_inv;
  }
 
  double getSupportRadius() const
  {
    return 1. / _r_inv;
  }
 
  double evaluate(double radius) const
  {
    return operator()(radius, _params);
  }
 
  PRECICE_HOST_DEVICE inline double operator()(const double radius, const RadialBasisParameters params) const
  {
    double       r_inv = params.parameter1;
    const double p     = radius * r_inv;
    if (p >= 1)
      return 0.0;
    return math::pow_int<4>(1.0 - p) * PRECICE_FMA(4, p, 1);
  }
 
  RadialBasisParameters getFunctionParameters()
  {
    return _params;
  }
 
private:
  double                _r_inv;
  RadialBasisParameters _params;
};

class CompactPolynomialC4 : public CompactSupportBase,
                            public DefiniteFunction<true> {
public:
  explicit CompactPolynomialC4(double supportRadius)
  {
    if (math::greater(supportRadius, 0.0)) {
      std::invalid_argument{
          "Support radius for radial-basis-function compact polynomial c4 has to be larger than zero. Please update the \"support-radius\" attribute."};
    }
 
    _r_inv             = 1. / supportRadius;
    _params.parameter1 = _r_inv;
  }
 
  double getSupportRadius() const
  {
    return 1. / _r_inv;
  }
 
  double evaluate(double radius) const
  {
    return operator()(radius, _params);
  }
 
  PRECICE_HOST_DEVICE inline double operator()(const double radius, const RadialBasisParameters params) const
  {
    double       r_inv = params.parameter1;
    const double p     = radius * r_inv;
    if (p >= 1)
      return 0.0;
    return math::pow_int<6>(1.0 - p) * (35 * math::pow_int<2>(p) + PRECICE_FMA(18, p, 3));
  }
 
  RadialBasisParameters getFunctionParameters()
  {
    return _params;
  }
 
private:
  double                _r_inv;
  RadialBasisParameters _params;
};

class CompactPolynomialC6 : public CompactSupportBase,
                            public DefiniteFunction<true> {
public:
  explicit CompactPolynomialC6(double supportRadius)
  {
    if (math::greater(supportRadius, 0.0)) {
      std::invalid_argument{
          "Support radius for radial-basis-function compact polynomial c6 has to be larger than zero. Please update the \"support-radius\" attribute."};
    }
    _r_inv             = 1. / supportRadius;
    _params.parameter1 = _r_inv;
  }
 
  double getSupportRadius() const
  {
    return 1. / _r_inv;
  }
 
  double evaluate(double radius) const
  {
    return operator()(radius, _params);
  }
 
  PRECICE_HOST_DEVICE inline double operator()(const double radius, const RadialBasisParameters params) const
  {
    double       r_inv = params.parameter1;
    const double p     = radius * r_inv;
    if (p >= 1)
      return 0.0;
    return math::pow_int<8>(1.0 - p) * (32.0 * math::pow_int<3>(p) + 25.0 * math::pow_int<2>(p) + PRECICE_FMA(8.0, p, 1.0));
  };
 
  const RadialBasisParameters getFunctionParameters()
  {
    return _params;
  }
 
private:
  double                _r_inv;
  RadialBasisParameters _params;
};

class CompactPolynomialC8 : public CompactSupportBase,
                            public DefiniteFunction<true> {
public:
  explicit CompactPolynomialC8(double supportRadius)
  {
    if (math::greater(supportRadius, 0.0)) {
      std::invalid_argument{
          "Support radius for radial-basis-function compact polynomial c6 has to be larger than zero. Please update the \"support-radius\" attribute."};
    }
    _r_inv             = 1. / supportRadius;
    _params.parameter1 = _r_inv;
  }
 
  double getSupportRadius() const
  {
    return 1. / _r_inv;
  }
 
  double evaluate(double radius) const
  {
    return operator()(radius, _params);
  }
 
  PRECICE_HOST_DEVICE inline double operator()(const double radius, const RadialBasisParameters params) const
  {
    double       r_inv = params.parameter1;
    const double p     = radius * r_inv;
    if (p >= 1)
      return 0.0;
    return math::pow_int<10>(1.0 - p) * (1287.0 * math::pow_int<4>(p) + 1350.0 * math::pow_int<3>(p) + 630.0 * math::pow_int<2>(p) + 150.0 * p + 15);
  };
 
  const RadialBasisParameters getFunctionParameters()
  {
    return _params;
  }
 
private:
  double                _r_inv;
  RadialBasisParameters _params;
};
 
#undef PRECICE_FMA
#undef PRECICE_LOG
#undef PRECICE_MEMORY_SPACE
#undef PRECICE_HOST_DEVICE
#undef NUMERICAL_ZERO_DIFFERENCE_DEVICE
 
} // namespace precice::mapping