From ddca9ab9c5973bc5b2fb2ddb9af78e820afc24af Mon Sep 17 00:00:00 2001 From: blegouix Date: Wed, 10 Jul 2024 13:46:39 +0200 Subject: [PATCH 01/17] init --- include/ddc/kernels/fft.hpp | 116 +++++++++++++++++++++++++++++++++--- 1 file changed, 108 insertions(+), 8 deletions(-) diff --git a/include/ddc/kernels/fft.hpp b/include/ddc/kernels/fft.hpp index 74a01083d..a56cdb76c 100644 --- a/include/ddc/kernels/fft.hpp +++ b/include/ddc/kernels/fft.hpp @@ -67,8 +67,18 @@ namespace ddc { template struct Fourier; -// named arguments for FFT (and their default values) +/** + * @brief A named argument too choose the direction of the FFT. + * + * @see kwArgs_core + */ enum class FFT_Direction { FORWARD, BACKWARD }; + +/** + * @brief A named argument too choose the type of normalization of the FFT. + * + * @see kwArgs_core + */ enum class FFT_Normalization { OFF, FORWARD, BACKWARD, ORTHO, FULL }; } // namespace ddc @@ -115,7 +125,12 @@ KOKKOS_FUNCTION constexpr T LastSelector(const T a, const T b) return LastSelector(a, b); } -// transform_type : trait to determine the type of transformation (R2C, C2R, C2C...) <- no information about base type (float or double) +// transform_type : +/** + * @brief A trait to identify the type of transformation (R2C, C2R, C2C...). + * + * It does not contain the information about the base type (float or double). + */ enum class TransformType { R2R, R2C, C2R, C2C }; template @@ -142,6 +157,14 @@ struct transform_type, Kokkos::complex> static constexpr TransformType value = TransformType::C2C; }; +/** + * @brief A trait to get the TransformType for the input and output types. + * + * Rely on if T1 and T2 are Kokkos::complex or not. + * + * @tparam T1 The input type. + * @tparam T2 The output type. + */ template constexpr TransformType transform_type_v = transform_type::value; @@ -318,6 +341,11 @@ hipfftResult _hipfftExec([[maybe_unused]] LastArg lastArg, Args... args) } #endif +/* + * @brief A structure embedding the configuration of the core FFT function, direction and type of normalization. + * + * @see FFT_core + */ struct kwArgs_core { ddc::FFT_Direction @@ -326,15 +354,17 @@ struct kwArgs_core }; // N,a,b from x_mesh +/// @brief Mesh size in dimension of interest. template int N(ddc::DiscreteDomain x_mesh) { static_assert( (is_uniform_point_sampling_v && ...), "DDimX dimensions should derive from UniformPointSampling"); - return ddc::get(x_mesh.extents()); + return static_cast(x_mesh.template extent()); } +/// @brief Lower cell of the domain along requested dimension. template double a(ddc::DiscreteDomain x_mesh) { @@ -346,6 +376,7 @@ double a(ddc::DiscreteDomain x_mesh) / 2 / (N(x_mesh) - 1); } +/// @brief Upper cell of the domain along requested dimension. template double b(ddc::DiscreteDomain x_mesh) { @@ -357,7 +388,7 @@ double b(ddc::DiscreteDomain x_mesh) / 2 / (N(x_mesh) - 1); } -// core +/// @brief Core internal function to perform the FFT. template void core( ExecSpace const& execSpace, @@ -581,6 +612,19 @@ void core( namespace ddc { +/** + * @brief Initialize a Fourier space. + * + * Initialize the (1D) discrete space representing the Fourier discrete dimension associated + * to the (1D) spatial mesh passed as argument. + * + * @tparam DDimFx A PeriodicSampling representing the Fourier discrete dimension. + * @tparam DDimX The type of the spatial discrete dimension. + * + * @param x_mesh The DiscreteDomain representing the (1D) spatial mesh. + * + * @see PeriodicSampling + */ template typename DDimFx::template Impl init_fourier_space( ddc::DiscreteDomain x_mesh) @@ -590,7 +634,7 @@ typename DDimFx::template Impl init_fourier_space( "DDimX dimensions should derive from UniformPointSampling"); static_assert( is_periodic_sampling_v, - "DDimFx dimensions should derive from PeriodicPointSampling"); + "DDimFx dimensions should derive from PeriodicSampling"); auto [impl, ddom] = DDimFx::template init( ddc::Coordinate(0), ddc::Coordinate( @@ -602,7 +646,18 @@ typename DDimFx::template Impl init_fourier_space( return std::move(impl); } -// FourierMesh, first element corresponds to mode 0 +/** + * @brief Get the Fourier mesh. + * + * Compute the Fourier (or spectral) mesh on which the Discrete Fourier Transform of a + * spatial discrete function is defined. + * + * @param x_mesh The DiscreteDomain representing the spatial mesh. + * @param C2C A flag indicating if a complex-to-complex DFT is going to be performed. Indeed, + * in this case the spatial and spectral meshes have same number of points, whereas for real-to-complex + * or complex-to-real DFT, each complex value of the Fourier-transformed function contains twice more + * information, and thus only N/2+1 points are needed. + */ template ddc::DiscreteDomain FourierMesh(ddc::DiscreteDomain x_mesh, bool C2C) { @@ -621,12 +676,37 @@ ddc::DiscreteDomain FourierMesh(ddc::DiscreteDomain x_mesh, ddc::detail::fft::N(x_mesh)))))...); } +/* + * @brief A structure embedding the configuration of the exposed FFT function with the type of normalization. + * + * @see fft, ifft + */ struct kwArgs_fft { ddc::FFT_Normalization normalization; }; -// FFT +/** + * @brief Perform a direct Fast Fourier Transform. + * + * Compute the discrete Fourier transform of a spatial function using the specialized implementation for the Kokkos::ExecutionSpace + * of the FFT algorithm. + * + * @tparam Tin The type of the input elements (float, Kokkos::complex, double or Kokkos::complex). + * @tparam Tout The type of the output elements (Kokkos::complex or Kokkos::complex). + * @tparam DDimFx... The parameter pack of the Fourier discrete dimensions. + * @tparam DDimX... The parameter pack of the spatial discrete dimensions. + * @tparam ExecSpace The type of the Kokkos::ExecutionSpace on which the FFT is performed. It determines which specialized + * backend is used (ie. fftw, cuFFT...). + * @tparam MemorySpace The type of the Kokkos::MemorySpace on which are stored the input and output discrete functions. + * @tparam layout_in The layout of the Chunkspan representing the input discrete function. + * @tparam layout_out The layout of the Chunkspan representing the output discrete function. + * + * @param execSpace The Kokkos::ExecutionSpace on which the FFT is performed. + * @param out The output discrete function, represented as a ChunkSpan storing values on a spectral mesh. + * @param in The input discrete function, represented as a ChunkSpan storing values on a spatial mesh. + * @param kwarg The kwArgs_fft configuring the FFT. + */ template < typename Tin, typename Tout, @@ -663,7 +743,27 @@ void fft( {ddc::FFT_Direction::FORWARD, kwargs.normalization}); } -// iFFT (deduced from the fact that "in" is identified as a function on the Fourier space) +/** + * @brief Perform an inverse Fast Fourier Transform. + * + * Compute the inverse discrete Fourier transform of a spectral function using the specialized implementation for the Kokkos::ExecutionSpace + * of the iFFT algorithm. + * + * @tparam Tin The type of the input elements (Kokkos::complex or Kokkos::complex). + * @tparam Tout The type of the output elements (float, Kokkos::complex, double or Kokkos::complex). + * @tparam DDimX... The parameter pack of the spatial discrete dimensions. + * @tparam DDimFx... The parameter pack of the Fourier discrete dimensions. + * @tparam ExecSpace The type of the Kokkos::ExecutionSpace on which the iFFT is performed. It determines which specialized + * backend is used (ie. fftw, cuFFT...). + * @tparam MemorySpace The type of the Kokkos::MemorySpace on which are stored the input and output discrete functions. + * @tparam layout_in The layout of the Chunkspan representing the input discrete function. + * @tparam layout_out The layout of the Chunkspan representing the output discrete function. + * + * @param execSpace The Kokkos::ExecutionSpace on which the iFFT is performed. + * @param out The output discrete function, represented as a ChunkSpan storing values on a spatial mesh. + * @param in The input discrete function, represented as a ChunkSpan storing values on a spectral mesh. + * @param kwarg The kwArgs_fft configuring the iFFT. + */ template < typename Tin, typename Tout, From e8e3dc4c09ed6c53e87c52bcbd8ba3b50d870963 Mon Sep 17 00:00:00 2001 From: blegouix Date: Wed, 10 Jul 2024 19:59:38 +0200 Subject: [PATCH 02/17] wip --- examples/heat_equation_spectral.cpp | 2 +- include/ddc/kernels/fft.hpp | 66 ++++++++++++++++++++++++++--- 2 files changed, 60 insertions(+), 8 deletions(-) diff --git a/examples/heat_equation_spectral.cpp b/examples/heat_equation_spectral.cpp index e906656b6..8a58bc169 100644 --- a/examples/heat_equation_spectral.cpp +++ b/examples/heat_equation_spectral.cpp @@ -68,7 +68,7 @@ void display(double time, ChunkType temp) / temp.domain().size(); std::cout << std::fixed << std::setprecision(3); std::cout << "At t = " << time << ",\n"; - std::cout << " * mean temperature = " << mean_temp << "\n"; + printf(" * mean temperature = %.15f \n", mean_temp); // take a slice in the middle of the box ddc::ChunkSpan temp_slice = temp[ddc::get_domain(temp).front() diff --git a/include/ddc/kernels/fft.hpp b/include/ddc/kernels/fft.hpp index a56cdb76c..6382fb811 100644 --- a/include/ddc/kernels/fft.hpp +++ b/include/ddc/kernels/fft.hpp @@ -64,18 +64,23 @@ static_assert(alignof(hipfftDoubleComplex) <= alignof(Kokkos::complex)); namespace ddc { // TODO : maybe transfert this somewhere else because Fourier space is not specific to FFT +/** + * @brief A templated tag representing a continuous dimension in the Fourier space associated to a continuous spatial dimension. + * + * @tparam The tag representing the spatial dimensions. + */ template struct Fourier; /** - * @brief A named argument too choose the direction of the FFT. + * @brief A named argument to choose the direction of the FFT. * * @see kwArgs_core */ enum class FFT_Direction { FORWARD, BACKWARD }; /** - * @brief A named argument too choose the type of normalization of the FFT. + * @brief A named argument to choose the type of normalization of the FFT. * * @see kwArgs_core */ @@ -125,7 +130,7 @@ KOKKOS_FUNCTION constexpr T LastSelector(const T a, const T b) return LastSelector(a, b); } -// transform_type : +// transform_type : /** * @brief A trait to identify the type of transformation (R2C, C2R, C2C...). * @@ -160,7 +165,7 @@ struct transform_type, Kokkos::complex> /** * @brief A trait to get the TransformType for the input and output types. * - * Rely on if T1 and T2 are Kokkos::complex or not. + * Internally check if T1 and T2 are Kokkos::complex or not. * * @tparam T1 The input type. * @tparam T2 The output type. @@ -354,7 +359,14 @@ struct kwArgs_core }; // N,a,b from x_mesh -/// @brief Mesh size in dimension of interest. +/** + * @brief Get the mesh size along a given dimension. + * + * @tparam DDim The dimension along which the mesh size is returned. + * @param x_mesh The mesh. + * + * @return The mesh size along the required dimension. + */ template int N(ddc::DiscreteDomain x_mesh) { @@ -364,7 +376,25 @@ int N(ddc::DiscreteDomain x_mesh) return static_cast(x_mesh.template extent()); } -/// @brief Lower cell of the domain along requested dimension. +/** + * @brief Get the lower boundary coordinate along a given dimension. + * + * The lower boundary of the spatial domain (which appears in Nyquist-Shannon theorem) is not + * xmin=ddc::coordinate(x_mesh.front()). Indeed, this coordinate identifies the lower cell, but + * the lower boundary is the left side of this lowest cell, which is a = xmin - cell_size/2, with + * cell_size = (xmax-xmin)/N. It leads to a = xmin-(b-a)/2N. The same derivation for the + * upper boundary coordinate gives b = xmax+(b-a)/2N. Inverting this linear system leads to: + * + * a = ((2N-1)*xmin-xmax)/2/(N-1) + * b = ((2N-1)*xmax-xmin)/2/(N-1) + * + * The current function implements the first equation. + * + * @tparam DDim The dimension along which the lower cell coordinate of the Fourier mesh is returned. + * @param x_mesh The spatial mesh. + * + * @return The mesh size along the required dimension. + */ template double a(ddc::DiscreteDomain x_mesh) { @@ -376,7 +406,25 @@ double a(ddc::DiscreteDomain x_mesh) / 2 / (N(x_mesh) - 1); } -/// @brief Upper cell of the domain along requested dimension. +/** + * @brief Get the upper boundary coordinate along a given dimension. + * + * The upper boundary of the spatial domain (which appears in Nyquist-Shannon theorem) is not + * xmax=ddc::coordinate(x_mesh.back()). Indeed, this coordinate identifies the upper cell, but + * the upper boundary is the right side of this upper cell, which is b = xmax + cell_size/2, with + * cell_size = (xmax-xmin)/N. It leads to b = xmax+(b-a)/2N. The same derivation for the + * lower boundary coordinate gives a = xmin-(b-a)/2N. Inverting this linear system leads to: + * + * a = ((2N-1)*xmin-xmax)/2/(N-1) + * b = ((2N-1)*xmax-xmin)/2/(N-1) + * + * The current function implements the second equation. + * + * @tparam DDim The dimension along which the upper cell coordinate of the Fourier mesh is returned. + * @param x_mesh The spatial mesh. + * + * @return The mesh size along the required dimension. + */ template double b(ddc::DiscreteDomain x_mesh) { @@ -618,6 +666,10 @@ namespace ddc { * Initialize the (1D) discrete space representing the Fourier discrete dimension associated * to the (1D) spatial mesh passed as argument. * + * The formula comes from the Nyquist-Shannon theorem: the lower bound of the spectral domain is kmin = -pi/lambda + * = -pi*N/(xmax-xmin). This is independent on the number of points in the spectral domain, which depends only on But in a discrete representation of a continous function, each itself is discretized (because the spatial domain is bounded), thus + * the leftest cell center coordinate of the spectral domain is kmin + (kmax-kmin)/(2N) = pi*(1-N)/(xmax-xmin) . + * * @tparam DDimFx A PeriodicSampling representing the Fourier discrete dimension. * @tparam DDimX The type of the spatial discrete dimension. * From fed827d673e99a0f814d6e5dae34f3fbded38995 Mon Sep 17 00:00:00 2001 From: blegouix Date: Thu, 11 Jul 2024 09:41:12 +0200 Subject: [PATCH 03/17] rc --- include/ddc/kernels/fft.hpp | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/include/ddc/kernels/fft.hpp b/include/ddc/kernels/fft.hpp index 6382fb811..fd0ec832f 100644 --- a/include/ddc/kernels/fft.hpp +++ b/include/ddc/kernels/fft.hpp @@ -664,12 +664,13 @@ namespace ddc { * @brief Initialize a Fourier space. * * Initialize the (1D) discrete space representing the Fourier discrete dimension associated - * to the (1D) spatial mesh passed as argument. - * - * The formula comes from the Nyquist-Shannon theorem: the lower bound of the spectral domain is kmin = -pi/lambda - * = -pi*N/(xmax-xmin). This is independent on the number of points in the spectral domain, which depends only on But in a discrete representation of a continous function, each itself is discretized (because the spatial domain is bounded), thus - * the leftest cell center coordinate of the spectral domain is kmin + (kmax-kmin)/(2N) = pi*(1-N)/(xmax-xmin) . + * to the (1D) spatial mesh passed as argument. It is a N-periodic PeriodicSampling defined between + * ka=0 and kb=2*N/(b-a)*pi. * + * This value for kb comes from the Nyquist-Shannon theorem: the period of the spectral domain + * is kb-ka = 2*pi/cell_size = 2*pi*N/(b-a). The PeriodicSampling then contains cells between coordinates + * k=0 and k=2*pi*(N-1)/(b-a), because the cell at coordinate k=2*pi*N/(b-a) is a periodic point (f(ka)=f(kb)). + * * @tparam DDimFx A PeriodicSampling representing the Fourier discrete dimension. * @tparam DDimX The type of the spatial discrete dimension. * @@ -735,7 +736,8 @@ ddc::DiscreteDomain FourierMesh(ddc::DiscreteDomain x_mesh, */ struct kwArgs_fft { - ddc::FFT_Normalization normalization; + ddc::FFT_Normalization + normalization; ///< Enum member to identify the type of normalization performed }; /** From c3cab86a50d8ab9479b314c1bfeb0090d58fb0dc Mon Sep 17 00:00:00 2001 From: blegouix Date: Thu, 11 Jul 2024 09:47:08 +0200 Subject: [PATCH 04/17] restore --- examples/heat_equation_spectral.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/heat_equation_spectral.cpp b/examples/heat_equation_spectral.cpp index 8a58bc169..e906656b6 100644 --- a/examples/heat_equation_spectral.cpp +++ b/examples/heat_equation_spectral.cpp @@ -68,7 +68,7 @@ void display(double time, ChunkType temp) / temp.domain().size(); std::cout << std::fixed << std::setprecision(3); std::cout << "At t = " << time << ",\n"; - printf(" * mean temperature = %.15f \n", mean_temp); + std::cout << " * mean temperature = " << mean_temp << "\n"; // take a slice in the middle of the box ddc::ChunkSpan temp_slice = temp[ddc::get_domain(temp).front() From 647c08ed90ff64757016d777b6bc524e3f6e7171 Mon Sep 17 00:00:00 2001 From: blegouix Date: Thu, 11 Jul 2024 10:44:55 +0200 Subject: [PATCH 05/17] enums --- include/ddc/kernels/fft.hpp | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/include/ddc/kernels/fft.hpp b/include/ddc/kernels/fft.hpp index fd0ec832f..869930269 100644 --- a/include/ddc/kernels/fft.hpp +++ b/include/ddc/kernels/fft.hpp @@ -77,14 +77,27 @@ struct Fourier; * * @see kwArgs_core */ -enum class FFT_Direction { FORWARD, BACKWARD }; +enum class FFT_Direction { + FORWARD, ///< Forward, corresponds to direct FFT up to normalization + BACKWARD ///< Backward, corresponds to inverse FFT up to normalization +}; /** * @brief A named argument to choose the type of normalization of the FFT. * * @see kwArgs_core */ -enum class FFT_Normalization { OFF, FORWARD, BACKWARD, ORTHO, FULL }; +enum class FFT_Normalization { + OFF, ///< No normalization + FORWARD, ///< Multiply by 1/N for forward FFT, no normalization for backward FFT + BACKWARD, ///< No normalization for forward FFT, multiply by 1/N for backward FFT + ORTHO, ///< Multiply by 1/sqrt(N) + FULL /**< + * Multiply by (b-a)/N/sqrt(2*pi) for forward FFT and sqrt(2*pi)/(b-a) for forward + * FFT. It is aligned with the usual definition of the (continuous) Fourier transform + * 1/sqrt(2*pi)*int f(x)*e^-ikx*dx, and thus preserves the gaussian function exp(-x^2/2) numerically. + */ +}; } // namespace ddc namespace ddc::detail::fft { From b61226e835267c59702bd1e8298dadb8f0f510c6 Mon Sep 17 00:00:00 2001 From: blegouix Date: Thu, 11 Jul 2024 11:15:33 +0200 Subject: [PATCH 06/17] fixes --- include/ddc/kernels/fft.hpp | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/include/ddc/kernels/fft.hpp b/include/ddc/kernels/fft.hpp index 869930269..0f648cacd 100644 --- a/include/ddc/kernels/fft.hpp +++ b/include/ddc/kernels/fft.hpp @@ -674,7 +674,7 @@ void core( namespace ddc { /** - * @brief Initialize a Fourier space. + * @brief Initialize a discrete Fourier space. * * Initialize the (1D) discrete space representing the Fourier discrete dimension associated * to the (1D) spatial mesh passed as argument. It is a N-periodic PeriodicSampling defined between @@ -689,6 +689,8 @@ namespace ddc { * * @param x_mesh The DiscreteDomain representing the (1D) spatial mesh. * + * @return The initialized Impl representing the discrete Fourier space. + * * @see PeriodicSampling */ template @@ -723,6 +725,8 @@ typename DDimFx::template Impl init_fourier_space( * in this case the spatial and spectral meshes have same number of points, whereas for real-to-complex * or complex-to-real DFT, each complex value of the Fourier-transformed function contains twice more * information, and thus only N/2+1 points are needed. + * + * @return The domain representing the Fourier mesh. */ template ddc::DiscreteDomain FourierMesh(ddc::DiscreteDomain x_mesh, bool C2C) @@ -742,7 +746,7 @@ ddc::DiscreteDomain FourierMesh(ddc::DiscreteDomain x_mesh, ddc::detail::fft::N(x_mesh)))))...); } -/* +/** * @brief A structure embedding the configuration of the exposed FFT function with the type of normalization. * * @see fft, ifft @@ -772,7 +776,7 @@ struct kwArgs_fft * @param execSpace The Kokkos::ExecutionSpace on which the FFT is performed. * @param out The output discrete function, represented as a ChunkSpan storing values on a spectral mesh. * @param in The input discrete function, represented as a ChunkSpan storing values on a spatial mesh. - * @param kwarg The kwArgs_fft configuring the FFT. + * @param kwargs The kwArgs_fft configuring the FFT. */ template < typename Tin, @@ -829,7 +833,7 @@ void fft( * @param execSpace The Kokkos::ExecutionSpace on which the iFFT is performed. * @param out The output discrete function, represented as a ChunkSpan storing values on a spatial mesh. * @param in The input discrete function, represented as a ChunkSpan storing values on a spectral mesh. - * @param kwarg The kwArgs_fft configuring the iFFT. + * @param kwargs The kwArgs_fft configuring the iFFT. */ template < typename Tin, From 7d56842dd5fb4fc55261d2f6fd10aebf0f5e5aac Mon Sep 17 00:00:00 2001 From: Baptiste Legouix Date: Thu, 11 Jul 2024 13:30:14 +0200 Subject: [PATCH 07/17] autoreview --- include/ddc/kernels/fft.hpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/include/ddc/kernels/fft.hpp b/include/ddc/kernels/fft.hpp index 0f648cacd..1587a2706 100644 --- a/include/ddc/kernels/fft.hpp +++ b/include/ddc/kernels/fft.hpp @@ -75,7 +75,7 @@ struct Fourier; /** * @brief A named argument to choose the direction of the FFT. * - * @see kwArgs_core + * @see kwArgs_core, kwArgs_fft */ enum class FFT_Direction { FORWARD, ///< Forward, corresponds to direct FFT up to normalization @@ -85,7 +85,7 @@ enum class FFT_Direction { /** * @brief A named argument to choose the type of normalization of the FFT. * - * @see kwArgs_core + * @see kwArgs_core, kwArgs_fft */ enum class FFT_Normalization { OFF, ///< No normalization @@ -360,7 +360,7 @@ hipfftResult _hipfftExec([[maybe_unused]] LastArg lastArg, Args... args) #endif /* - * @brief A structure embedding the configuration of the core FFT function, direction and type of normalization. + * @brief A structure embedding the configuration of the core FFT function: direction and type of normalization. * * @see FFT_core */ @@ -395,7 +395,7 @@ int N(ddc::DiscreteDomain x_mesh) * The lower boundary of the spatial domain (which appears in Nyquist-Shannon theorem) is not * xmin=ddc::coordinate(x_mesh.front()). Indeed, this coordinate identifies the lower cell, but * the lower boundary is the left side of this lowest cell, which is a = xmin - cell_size/2, with - * cell_size = (xmax-xmin)/N. It leads to a = xmin-(b-a)/2N. The same derivation for the + * cell_size = (b-a)/N. It leads to a = xmin-(b-a)/2N. The same derivation for the * upper boundary coordinate gives b = xmax+(b-a)/2N. Inverting this linear system leads to: * * a = ((2N-1)*xmin-xmax)/2/(N-1) @@ -425,7 +425,7 @@ double a(ddc::DiscreteDomain x_mesh) * The upper boundary of the spatial domain (which appears in Nyquist-Shannon theorem) is not * xmax=ddc::coordinate(x_mesh.back()). Indeed, this coordinate identifies the upper cell, but * the upper boundary is the right side of this upper cell, which is b = xmax + cell_size/2, with - * cell_size = (xmax-xmin)/N. It leads to b = xmax+(b-a)/2N. The same derivation for the + * cell_size = (b-a)/N. It leads to b = xmax+(b-a)/2N. The same derivation for the * lower boundary coordinate gives a = xmin-(b-a)/2N. Inverting this linear system leads to: * * a = ((2N-1)*xmin-xmax)/2/(N-1) From 762f9013685fa95b47c2df4da657b9fec942626a Mon Sep 17 00:00:00 2001 From: blegouix Date: Thu, 11 Jul 2024 13:33:28 +0200 Subject: [PATCH 08/17] precision --- include/ddc/kernels/fft.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/include/ddc/kernels/fft.hpp b/include/ddc/kernels/fft.hpp index 1587a2706..9104f2172 100644 --- a/include/ddc/kernels/fft.hpp +++ b/include/ddc/kernels/fft.hpp @@ -720,6 +720,8 @@ typename DDimFx::template Impl init_fourier_space( * Compute the Fourier (or spectral) mesh on which the Discrete Fourier Transform of a * spatial discrete function is defined. * + * The uid identifies the mode (ie. ddc::DiscreteElement(0) corresponds to mode 0). + * * @param x_mesh The DiscreteDomain representing the spatial mesh. * @param C2C A flag indicating if a complex-to-complex DFT is going to be performed. Indeed, * in this case the spatial and spectral meshes have same number of points, whereas for real-to-complex From 9b4d39e5a679bdee20c5f166c6ef9dec15e23f67 Mon Sep 17 00:00:00 2001 From: Baptiste Legouix Date: Thu, 11 Jul 2024 13:35:18 +0200 Subject: [PATCH 09/17] Update include/ddc/kernels/fft.hpp --- include/ddc/kernels/fft.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/ddc/kernels/fft.hpp b/include/ddc/kernels/fft.hpp index 9104f2172..c4491e9c4 100644 --- a/include/ddc/kernels/fft.hpp +++ b/include/ddc/kernels/fft.hpp @@ -143,7 +143,7 @@ KOKKOS_FUNCTION constexpr T LastSelector(const T a, const T b) return LastSelector(a, b); } -// transform_type : +// transform_type: /** * @brief A trait to identify the type of transformation (R2C, C2R, C2C...). * From b6c164ce4d240a5904a0d9577cbfab0e8e5ac5bc Mon Sep 17 00:00:00 2001 From: blegouix Date: Thu, 11 Jul 2024 13:53:10 +0200 Subject: [PATCH 10/17] precisions --- include/ddc/kernels/fft.hpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/include/ddc/kernels/fft.hpp b/include/ddc/kernels/fft.hpp index c4491e9c4..e28ddb771 100644 --- a/include/ddc/kernels/fft.hpp +++ b/include/ddc/kernels/fft.hpp @@ -726,7 +726,8 @@ typename DDimFx::template Impl init_fourier_space( * @param C2C A flag indicating if a complex-to-complex DFT is going to be performed. Indeed, * in this case the spatial and spectral meshes have same number of points, whereas for real-to-complex * or complex-to-real DFT, each complex value of the Fourier-transformed function contains twice more - * information, and thus only N/2+1 points are needed. + * information, and thus only half (actually Nx*Ny*(Nz/2+1) for 3D R2C FFT to take in account mode 0) + * values are needed (cf. DFT conjugate symmetry property for more information about this). * * @return The domain representing the Fourier mesh. */ @@ -822,6 +823,8 @@ void fft( * Compute the inverse discrete Fourier transform of a spectral function using the specialized implementation for the Kokkos::ExecutionSpace * of the iFFT algorithm. * + * /!\ C2R iFFT does NOT preserve input ! + * * @tparam Tin The type of the input elements (Kokkos::complex or Kokkos::complex). * @tparam Tout The type of the output elements (float, Kokkos::complex, double or Kokkos::complex). * @tparam DDimX... The parameter pack of the spatial discrete dimensions. From 7e715d8bd0b79fdbb7c1a07af4384a11feb9425e Mon Sep 17 00:00:00 2001 From: Baptiste Legouix Date: Thu, 11 Jul 2024 13:55:00 +0200 Subject: [PATCH 11/17] Update include/ddc/kernels/fft.hpp --- include/ddc/kernels/fft.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/ddc/kernels/fft.hpp b/include/ddc/kernels/fft.hpp index e28ddb771..7ba100de2 100644 --- a/include/ddc/kernels/fft.hpp +++ b/include/ddc/kernels/fft.hpp @@ -394,7 +394,7 @@ int N(ddc::DiscreteDomain x_mesh) * * The lower boundary of the spatial domain (which appears in Nyquist-Shannon theorem) is not * xmin=ddc::coordinate(x_mesh.front()). Indeed, this coordinate identifies the lower cell, but - * the lower boundary is the left side of this lowest cell, which is a = xmin - cell_size/2, with + * the lower boundary is the left side of this lower cell, which is a = xmin - cell_size/2, with * cell_size = (b-a)/N. It leads to a = xmin-(b-a)/2N. The same derivation for the * upper boundary coordinate gives b = xmax+(b-a)/2N. Inverting this linear system leads to: * From ab4415279208aeabfb23366d9cb97a45384ec65a Mon Sep 17 00:00:00 2001 From: Baptiste Legouix Date: Thu, 11 Jul 2024 13:56:56 +0200 Subject: [PATCH 12/17] Apply suggestions from code review --- include/ddc/kernels/fft.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/ddc/kernels/fft.hpp b/include/ddc/kernels/fft.hpp index 7ba100de2..6df761f85 100644 --- a/include/ddc/kernels/fft.hpp +++ b/include/ddc/kernels/fft.hpp @@ -406,7 +406,7 @@ int N(ddc::DiscreteDomain x_mesh) * @tparam DDim The dimension along which the lower cell coordinate of the Fourier mesh is returned. * @param x_mesh The spatial mesh. * - * @return The mesh size along the required dimension. + * @return The lower boundary along the required dimension. */ template double a(ddc::DiscreteDomain x_mesh) @@ -436,7 +436,7 @@ double a(ddc::DiscreteDomain x_mesh) * @tparam DDim The dimension along which the upper cell coordinate of the Fourier mesh is returned. * @param x_mesh The spatial mesh. * - * @return The mesh size along the required dimension. + * @return The upper boundary along the required dimension. */ template double b(ddc::DiscreteDomain x_mesh) From 95f020cb70ff8c3a73a14ec80de32992f7bdc77c Mon Sep 17 00:00:00 2001 From: Baptiste Legouix Date: Fri, 12 Jul 2024 08:45:43 +0200 Subject: [PATCH 13/17] Apply suggestions from code review Co-authored-by: Thomas Padioleau --- include/ddc/kernels/fft.hpp | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/include/ddc/kernels/fft.hpp b/include/ddc/kernels/fft.hpp index 6df761f85..cac2967a9 100644 --- a/include/ddc/kernels/fft.hpp +++ b/include/ddc/kernels/fft.hpp @@ -63,11 +63,10 @@ static_assert(alignof(hipfftDoubleComplex) <= alignof(Kokkos::complex)); #endif namespace ddc { -// TODO : maybe transfert this somewhere else because Fourier space is not specific to FFT /** - * @brief A templated tag representing a continuous dimension in the Fourier space associated to a continuous spatial dimension. + * @brief A templated tag representing a continuous dimension in the Fourier space associated to the original continuous dimension. * - * @tparam The tag representing the spatial dimensions. + * @tparam The tag representing the original dimension. */ template struct Fourier; @@ -143,7 +142,6 @@ KOKKOS_FUNCTION constexpr T LastSelector(const T a, const T b) return LastSelector(a, b); } -// transform_type: /** * @brief A trait to identify the type of transformation (R2C, C2R, C2C...). * @@ -371,7 +369,6 @@ struct kwArgs_core ddc::FFT_Normalization normalization; }; -// N,a,b from x_mesh /** * @brief Get the mesh size along a given dimension. * @@ -674,7 +671,7 @@ void core( namespace ddc { /** - * @brief Initialize a discrete Fourier space. + * @brief Initialize a Fourier discrete dimension. * * Initialize the (1D) discrete space representing the Fourier discrete dimension associated * to the (1D) spatial mesh passed as argument. It is a N-periodic PeriodicSampling defined between From 66383203473e8713b54e61fbc25bb37d6ec51127 Mon Sep 17 00:00:00 2001 From: blegouix Date: Fri, 12 Jul 2024 12:26:45 +0200 Subject: [PATCH 14/17] remove "spatial" terminology --- include/ddc/kernels/fft.hpp | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/include/ddc/kernels/fft.hpp b/include/ddc/kernels/fft.hpp index cac2967a9..437dfc2a9 100644 --- a/include/ddc/kernels/fft.hpp +++ b/include/ddc/kernels/fft.hpp @@ -389,7 +389,7 @@ int N(ddc::DiscreteDomain x_mesh) /** * @brief Get the lower boundary coordinate along a given dimension. * - * The lower boundary of the spatial domain (which appears in Nyquist-Shannon theorem) is not + * The lower boundary of the domain (which appears in Nyquist-Shannon theorem) is not * xmin=ddc::coordinate(x_mesh.front()). Indeed, this coordinate identifies the lower cell, but * the lower boundary is the left side of this lower cell, which is a = xmin - cell_size/2, with * cell_size = (b-a)/N. It leads to a = xmin-(b-a)/2N. The same derivation for the @@ -401,7 +401,7 @@ int N(ddc::DiscreteDomain x_mesh) * The current function implements the first equation. * * @tparam DDim The dimension along which the lower cell coordinate of the Fourier mesh is returned. - * @param x_mesh The spatial mesh. + * @param x_mesh The mesh. * * @return The lower boundary along the required dimension. */ @@ -419,7 +419,7 @@ double a(ddc::DiscreteDomain x_mesh) /** * @brief Get the upper boundary coordinate along a given dimension. * - * The upper boundary of the spatial domain (which appears in Nyquist-Shannon theorem) is not + * The upper boundary of the domain (which appears in Nyquist-Shannon theorem) is not * xmax=ddc::coordinate(x_mesh.back()). Indeed, this coordinate identifies the upper cell, but * the upper boundary is the right side of this upper cell, which is b = xmax + cell_size/2, with * cell_size = (b-a)/N. It leads to b = xmax+(b-a)/2N. The same derivation for the @@ -431,7 +431,7 @@ double a(ddc::DiscreteDomain x_mesh) * The current function implements the second equation. * * @tparam DDim The dimension along which the upper cell coordinate of the Fourier mesh is returned. - * @param x_mesh The spatial mesh. + * @param x_mesh The mesh. * * @return The upper boundary along the required dimension. */ @@ -674,7 +674,7 @@ namespace ddc { * @brief Initialize a Fourier discrete dimension. * * Initialize the (1D) discrete space representing the Fourier discrete dimension associated - * to the (1D) spatial mesh passed as argument. It is a N-periodic PeriodicSampling defined between + * to the (1D) mesh passed as argument. It is a N-periodic PeriodicSampling defined between * ka=0 and kb=2*N/(b-a)*pi. * * This value for kb comes from the Nyquist-Shannon theorem: the period of the spectral domain @@ -682,9 +682,9 @@ namespace ddc { * k=0 and k=2*pi*(N-1)/(b-a), because the cell at coordinate k=2*pi*N/(b-a) is a periodic point (f(ka)=f(kb)). * * @tparam DDimFx A PeriodicSampling representing the Fourier discrete dimension. - * @tparam DDimX The type of the spatial discrete dimension. + * @tparam DDimX The type of the original discrete dimension. * - * @param x_mesh The DiscreteDomain representing the (1D) spatial mesh. + * @param x_mesh The DiscreteDomain representing the (1D) original mesh. * * @return The initialized Impl representing the discrete Fourier space. * @@ -715,13 +715,13 @@ typename DDimFx::template Impl init_fourier_space( * @brief Get the Fourier mesh. * * Compute the Fourier (or spectral) mesh on which the Discrete Fourier Transform of a - * spatial discrete function is defined. + * discrete function is defined. * * The uid identifies the mode (ie. ddc::DiscreteElement(0) corresponds to mode 0). * - * @param x_mesh The DiscreteDomain representing the spatial mesh. + * @param x_mesh The DiscreteDomain representing the original mesh. * @param C2C A flag indicating if a complex-to-complex DFT is going to be performed. Indeed, - * in this case the spatial and spectral meshes have same number of points, whereas for real-to-complex + * in this case the two meshes have same number of points, whereas for real-to-complex * or complex-to-real DFT, each complex value of the Fourier-transformed function contains twice more * information, and thus only half (actually Nx*Ny*(Nz/2+1) for 3D R2C FFT to take in account mode 0) * values are needed (cf. DFT conjugate symmetry property for more information about this). @@ -760,13 +760,13 @@ struct kwArgs_fft /** * @brief Perform a direct Fast Fourier Transform. * - * Compute the discrete Fourier transform of a spatial function using the specialized implementation for the Kokkos::ExecutionSpace + * Compute the discrete Fourier transform of a function using the specialized implementation for the Kokkos::ExecutionSpace * of the FFT algorithm. * * @tparam Tin The type of the input elements (float, Kokkos::complex, double or Kokkos::complex). * @tparam Tout The type of the output elements (Kokkos::complex or Kokkos::complex). * @tparam DDimFx... The parameter pack of the Fourier discrete dimensions. - * @tparam DDimX... The parameter pack of the spatial discrete dimensions. + * @tparam DDimX... The parameter pack of the original discrete dimensions. * @tparam ExecSpace The type of the Kokkos::ExecutionSpace on which the FFT is performed. It determines which specialized * backend is used (ie. fftw, cuFFT...). * @tparam MemorySpace The type of the Kokkos::MemorySpace on which are stored the input and output discrete functions. @@ -775,7 +775,7 @@ struct kwArgs_fft * * @param execSpace The Kokkos::ExecutionSpace on which the FFT is performed. * @param out The output discrete function, represented as a ChunkSpan storing values on a spectral mesh. - * @param in The input discrete function, represented as a ChunkSpan storing values on a spatial mesh. + * @param in The input discrete function, represented as a ChunkSpan storing values on a mesh. * @param kwargs The kwArgs_fft configuring the FFT. */ template < @@ -824,7 +824,7 @@ void fft( * * @tparam Tin The type of the input elements (Kokkos::complex or Kokkos::complex). * @tparam Tout The type of the output elements (float, Kokkos::complex, double or Kokkos::complex). - * @tparam DDimX... The parameter pack of the spatial discrete dimensions. + * @tparam DDimX... The parameter pack of the original discrete dimensions. * @tparam DDimFx... The parameter pack of the Fourier discrete dimensions. * @tparam ExecSpace The type of the Kokkos::ExecutionSpace on which the iFFT is performed. It determines which specialized * backend is used (ie. fftw, cuFFT...). @@ -833,7 +833,7 @@ void fft( * @tparam layout_out The layout of the Chunkspan representing the output discrete function. * * @param execSpace The Kokkos::ExecutionSpace on which the iFFT is performed. - * @param out The output discrete function, represented as a ChunkSpan storing values on a spatial mesh. + * @param out The output discrete function, represented as a ChunkSpan storing values on a mesh. * @param in The input discrete function, represented as a ChunkSpan storing values on a spectral mesh. * @param kwargs The kwArgs_fft configuring the iFFT. */ From b8ffaedff1321c0b58b19f997c214e96d7ce70b3 Mon Sep 17 00:00:00 2001 From: Baptiste Legouix Date: Tue, 16 Jul 2024 09:00:31 +0200 Subject: [PATCH 15/17] Update include/ddc/kernels/fft.hpp --- include/ddc/kernels/fft.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/ddc/kernels/fft.hpp b/include/ddc/kernels/fft.hpp index 437dfc2a9..15eac7053 100644 --- a/include/ddc/kernels/fft.hpp +++ b/include/ddc/kernels/fft.hpp @@ -92,7 +92,7 @@ enum class FFT_Normalization { BACKWARD, ///< No normalization for forward FFT, multiply by 1/N for backward FFT ORTHO, ///< Multiply by 1/sqrt(N) FULL /**< - * Multiply by (b-a)/N/sqrt(2*pi) for forward FFT and sqrt(2*pi)/(b-a) for forward + * Multiply by (b-a)/N/sqrt(2*pi) for forward FFT and (kb-ka)/N/sqrt(2*pi) for backward * FFT. It is aligned with the usual definition of the (continuous) Fourier transform * 1/sqrt(2*pi)*int f(x)*e^-ikx*dx, and thus preserves the gaussian function exp(-x^2/2) numerically. */ From 78a53410024d1e634895cc8377d412ea2aef46c0 Mon Sep 17 00:00:00 2001 From: Baptiste Legouix Date: Tue, 16 Jul 2024 09:00:45 +0200 Subject: [PATCH 16/17] Update include/ddc/kernels/fft.hpp --- include/ddc/kernels/fft.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/ddc/kernels/fft.hpp b/include/ddc/kernels/fft.hpp index 15eac7053..1926a72ea 100644 --- a/include/ddc/kernels/fft.hpp +++ b/include/ddc/kernels/fft.hpp @@ -92,7 +92,7 @@ enum class FFT_Normalization { BACKWARD, ///< No normalization for forward FFT, multiply by 1/N for backward FFT ORTHO, ///< Multiply by 1/sqrt(N) FULL /**< - * Multiply by (b-a)/N/sqrt(2*pi) for forward FFT and (kb-ka)/N/sqrt(2*pi) for backward + * Multiply by (b-a)/N/sqrt(2*pi) for forward FFT and (kb-ka)/N/sqrt(2*pi) for backward * FFT. It is aligned with the usual definition of the (continuous) Fourier transform * 1/sqrt(2*pi)*int f(x)*e^-ikx*dx, and thus preserves the gaussian function exp(-x^2/2) numerically. */ From e56516916e60fcfaf6fe2eefc014ee53a82d4ebe Mon Sep 17 00:00:00 2001 From: blegouix Date: Fri, 19 Jul 2024 11:30:43 +0000 Subject: [PATCH 17/17] core -> impl --- include/ddc/kernels/fft.hpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/include/ddc/kernels/fft.hpp b/include/ddc/kernels/fft.hpp index 4c89deccc..5566bc85a 100644 --- a/include/ddc/kernels/fft.hpp +++ b/include/ddc/kernels/fft.hpp @@ -74,7 +74,7 @@ struct Fourier; /** * @brief A named argument to choose the direction of the FFT. * - * @see kwArgs_core, kwArgs_fft + * @see kwArgs_impl, kwArgs_fft */ enum class FFT_Direction { FORWARD, ///< Forward, corresponds to direct FFT up to normalization @@ -84,7 +84,7 @@ enum class FFT_Direction { /** * @brief A named argument to choose the type of normalization of the FFT. * - * @see kwArgs_core, kwArgs_fft + * @see kwArgs_impl, kwArgs_fft */ enum class FFT_Normalization { OFF, ///< No normalization. Un-normalized FFT is sum_j f(x_j)*e^-ikx_j @@ -358,11 +358,11 @@ hipfftResult _hipfftExec([[maybe_unused]] LastArg lastArg, Args... args) #endif /* - * @brief A structure embedding the configuration of the core FFT function: direction and type of normalization. + * @brief A structure embedding the configuration of the impl FFT function: direction and type of normalization. * - * @see FFT_core + * @see FFT_impl */ -struct kwArgs_core +struct kwArgs_impl { ddc::FFT_Direction direction; // Only effective for C2C transform and for normalization BACKWARD and FORWARD @@ -388,12 +388,12 @@ int N(ddc::DiscreteDomain x_mesh) /// @brief Core internal function to perform the FFT. template -void core( +void impl( ExecSpace const& execSpace, Tout* out_data, Tin* in_data, ddc::DiscreteDomain mesh, - const kwArgs_core& kwargs) + const kwArgs_impl& kwargs) { static_assert( std::is_same_v, float> || std::is_same_v, double>, @@ -746,7 +746,7 @@ void fft( (is_periodic_sampling_v && ...), "DDimFx dimensions should derive from PeriodicPointSampling"); - ddc::detail::fft::core( + ddc::detail::fft::impl( execSpace, out.data_handle(), in.data_handle(), @@ -805,7 +805,7 @@ void ifft( (is_periodic_sampling_v && ...), "DDimFx dimensions should derive from PeriodicPointSampling"); - ddc::detail::fft::core( + ddc::detail::fft::impl( execSpace, out.data_handle(), in.data_handle(),