combined_coordinates.hh

Go to the documentation of this file.

#pragma once
 
// standard library
#include <cmath>
#include <tuple>
 
// DiFfRG
#include <DiFfRG/common/cuda_prefix.hh>
#include <DiFfRG/common/math.hh>
#include <DiFfRG/discretization/grid/coordinates.hh>
#include <DiFfRG/discretization/grid/stack_coordinates.hh>
 
namespace DiFfRG
{
  template <typename Idx, typename NT> class BosonicCoordinates1DFiniteT
  {
  public:
    using ctype = NT;
    static constexpr uint dim = 2;
 
    BosonicCoordinates1DFiniteT(Idx m_start, Idx m_stop, NT m_T, uint grid_extent, NT p_start, NT p_stop, NT p_bias)
        : m_start(m_start), m_stop(m_stop), m_T(m_T), p_start(p_start), p_stop(p_stop), p_bias(p_bias),
          grid_extent(grid_extent), m_size(m_stop - m_start), matsubara_values(m_start, m_stop, m_T),
          logarithmic_coordinates(grid_extent, p_start, p_stop, p_bias)
    {
    }
 
    template <typename Idx2, typename NT2>
    BosonicCoordinates1DFiniteT(const BosonicCoordinates1DFiniteT<Idx2, NT2> &other)
        : BosonicCoordinates1DFiniteT(other.m_start, other.m_stop, other.m_T, other.grid_extent, other.p_start,
                                      other.p_stop, other.p_bias)
    {
    }
 
    template <typename NT2>
    BosonicCoordinates1DFiniteT(const LogarithmicCoordinates1D<NT2> &logarithmic_coordinates, Idx m_start, Idx m_stop,
                                NT m_T)
        : m_start(m_start), m_stop(m_stop), m_T(m_T), p_start(logarithmic_coordinates.start),
          p_stop(logarithmic_coordinates.stop), p_bias(logarithmic_coordinates.bias),
          grid_extent(logarithmic_coordinates.size()), m_size(m_stop - m_start), matsubara_values(m_start, m_stop, m_T),
          logarithmic_coordinates(logarithmic_coordinates)
    {
    }
 
    std::array<NT, 2> __forceinline__ __device__ __host__ forward(const uint m, const uint p) const
    {
      return {matsubara_values.forward(m), logarithmic_coordinates.forward(p)};
    }
 
    std::tuple<Idx, NT> __forceinline__ __device__ __host__ backward(const NT m, const NT p) const
    {
      Idx m_idx = matsubara_values.backward(m);
      NT p_idx = 0;
      if (m_idx >= int(m_size)) {
        m_idx = m_size - 1;
        const auto new_p = std::sqrt(powr<2>(p) + powr<2>(m - matsubara_values.forward(m_idx)));
        p_idx = logarithmic_coordinates.backward(new_p);
      } else if (m_idx < 0) {
        m_idx = 0;
        const auto new_p = std::sqrt(powr<2>(p) + powr<2>(m - matsubara_values.forward(m_idx)));
        p_idx = logarithmic_coordinates.backward(new_p);
      } else
        p_idx = logarithmic_coordinates.backward(p);
      return {m_idx, p_idx};
    }
 
    uint size() const { return m_size * grid_extent; }
    std::array<uint, 2> sizes() const { return {{m_size, grid_extent}}; }
 
    const NT m_start, m_stop, m_T, p_start, p_stop, p_bias;
 
  private:
    const uint grid_extent, m_size;
 
    BosonicMatsubaraValues<Idx, NT> matsubara_values;
    LogarithmicCoordinates1D<NT> logarithmic_coordinates;
  };
 
  template <typename Idx, typename NT> class FermionicCoordinates1DFiniteT
  {
  public:
    using ctype = NT;
    static constexpr uint dim = 2;
 
    FermionicCoordinates1DFiniteT(Idx m_start, Idx m_stop, NT m_T, uint grid_extent, NT p_start, NT p_stop, NT p_bias)
        : m_start(m_start), m_stop(m_stop), m_T(m_T), p_start(p_start), p_stop(p_stop), p_bias(p_bias),
          grid_extent(grid_extent), m_size(m_stop - m_start), matsubara_values(m_start, m_stop, m_T),
          logarithmic_coordinates(grid_extent, p_start, p_stop, p_bias)
    {
    }
 
    template <typename Idx2, typename NT2>
    FermionicCoordinates1DFiniteT(const FermionicCoordinates1DFiniteT<Idx2, NT2> &other)
        : FermionicCoordinates1DFiniteT(other.m_start, other.m_stop, other.m_T, other.grid_extent, other.p_start,
                                        other.p_stop, other.p_bias)
    {
    }
 
    template <typename NT2>
    FermionicCoordinates1DFiniteT(const LogarithmicCoordinates1D<NT2> &logarithmic_coordinates, Idx m_start, Idx m_stop,
                                  NT m_T)
        : m_start(m_start), m_stop(m_stop), m_T(m_T), p_start(logarithmic_coordinates.start),
          p_stop(logarithmic_coordinates.stop), p_bias(logarithmic_coordinates.bias),
          grid_extent(logarithmic_coordinates.size()), m_size(m_stop - m_start), matsubara_values(m_start, m_stop, m_T),
          logarithmic_coordinates(logarithmic_coordinates)
    {
    }
 
    std::array<NT, 2> __forceinline__ __device__ __host__ forward(const uint m, const uint p) const
    {
      return {matsubara_values.forward(m), logarithmic_coordinates.forward(p)};
    }
 
    std::tuple<Idx, NT> __forceinline__ __device__ __host__ backward(const NT m, const NT p) const
    {
      Idx m_idx = matsubara_values.backward(m);
      NT p_idx = 0;
      if (m_idx >= int(m_size)) {
        m_idx = m_size - 1;
        const auto new_p = std::sqrt(powr<2>(p) + powr<2>(m - matsubara_values.forward(m_idx)));
        p_idx = logarithmic_coordinates.backward(new_p);
      } else if (m_idx < 0) {
        m_idx = 0;
        const auto new_p = std::sqrt(powr<2>(p) + powr<2>(m - matsubara_values.forward(m_idx)));
        p_idx = logarithmic_coordinates.backward(new_p);
      } else
        p_idx = logarithmic_coordinates.backward(p);
      return {m_idx, p_idx};
    }
 
    uint size() const { return m_size * grid_extent; }
    std::array<uint, 2> sizes() const { return {{m_size, grid_extent}}; }
 
    const NT m_start, m_stop, m_T, p_start, p_stop, p_bias;
 
  private:
    const uint grid_extent, m_size;
 
    FermionicMatsubaraValues<Idx, NT> matsubara_values;
    LogarithmicCoordinates1D<NT> logarithmic_coordinates;
  };
} // namespace DiFfRG