Skip to content

Commit

Permalink
Move linear interpolation functions to ablastr (#5714)
Browse files Browse the repository at this point in the history 
Pure functions like `linear_interp`, `bilinear_interp`, and
`trilinear_interp` are very general. Therefore, we can consider to move
them to `ablastr::math`. Besides, I will need some of these functions to
move `picsar_qed` inside `ablastr`
(#5677)
  • Loading branch information
@lucafedeli88
lucafedeli88 authored Feb 28, 2025
1 parent cfd9d1d commit 2386c18
Show file tree
Hide file tree
Showing 3 changed files with 44 additions and 39 deletions.
6 changes: 3 additions & 3 deletions Source/Initialization/WarpXInitData.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,6 @@
#include "Initialization/ExternalField.H"
#include "Initialization/DivCleaner/ProjectionDivCleaner.H"
#include "Particles/MultiParticleContainer.H"
#include "Utils/Algorithms/LinearInterpolation.H"
#include "Utils/Logo/GetLogo.H"
#include "Utils/Parser/ParserUtils.H"
#include "Utils/TextMsg.H"
Expand All @@ -40,6 +39,7 @@
#include "Python/callbacks.H"

#include <ablastr/fields/MultiFabRegister.H>
#include <ablastr/math/LinearInterpolation.H>
#include <ablastr/parallelization/MPIInitHelpers.H>
#include <ablastr/utils/Communication.H>
#include <ablastr/utils/UsedInputsFile.H>
Expand Down Expand Up @@ -1698,7 +1698,7 @@ WarpX::ReadExternalFieldFromFile (
f01 = fc_array(0, iz , ir+1),
f10 = fc_array(0, iz+1, ir ),
f11 = fc_array(0, iz+1, ir+1);
mffab(i,j,k) = static_cast<amrex::Real>(utils::algorithms::bilinear_interp<double>
mffab(i,j,k) = static_cast<amrex::Real>(ablastr::math::bilinear_interp<double>
(xx0, xx0+file_dr, xx1, xx1+file_dz,
f00, f01, f10, f11,
x0, x1));
Expand All @@ -1713,7 +1713,7 @@ WarpX::ReadExternalFieldFromFile (
f101 = fc_array(iz+1, iy , ix+1),
f110 = fc_array(iz , iy+1, ix+1),
f111 = fc_array(iz+1, iy+1, ix+1);
mffab(i,j,k) = static_cast<amrex::Real>(utils::algorithms::trilinear_interp<double>
mffab(i,j,k) = static_cast<amrex::Real>(ablastr::math::trilinear_interp<double>
(xx0, xx0+file_dx, xx1, xx1+file_dy, xx2, xx2+file_dz,
f000, f001, f010, f011, f100, f101, f110, f111,
x0, x1, x2));
Expand Down
58 changes: 30 additions & 28 deletions Source/Laser/LaserProfilesImpl/LaserProfileFromFile.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,12 +6,12 @@
*/
#include "Laser/LaserProfiles.H"

#include "Utils/Algorithms/LinearInterpolation.H"
#include "Utils/Parser/ParserUtils.H"
#include "Utils/TextMsg.H"
#include "Utils/WarpX_Complex.H"
#include "Utils/WarpXConst.H"

#include <ablastr/math/LinearInterpolation.H>
#include <ablastr/warn_manager/WarnManager.H>

#include <AMReX.H>
Expand Down Expand Up @@ -489,7 +489,7 @@ WarpXLaserProfiles::FromFileLaserProfile::internal_fill_amplitude_uniform_cartes
(i_interp-tmp_idx_first_time)*tmp_nx*tmp_ny+
j_interp*tmp_nx + k_interp;
};
const Complex val = utils::algorithms::trilinear_interp(
const Complex val = ablastr::math::trilinear_interp(
t_left, t_right,
x_0, x_1,
y_0, y_1,
Expand Down Expand Up @@ -574,33 +574,35 @@ WarpXLaserProfiles::FromFileLaserProfile::internal_fill_amplitude_uniform_cylind
Complex fact = Complex{costheta, sintheta};

// azimuthal mode 0
val += utils::algorithms::bilinear_interp(
t_left, t_right,
r_0, r_1,
p_E_lasy_data[idx(0, idx_t_left, idx_r_left)],
p_E_lasy_data[idx(0, idx_t_left, idx_r_right)],
p_E_lasy_data[idx(0, idx_t_right, idx_r_left)],
p_E_lasy_data[idx(0, idx_t_right, idx_r_right)],
t, Rp_i);
val +=
ablastr::math::bilinear_interp(
t_left, t_right,
r_0, r_1,
p_E_lasy_data[idx(0, idx_t_left, idx_r_left)],
p_E_lasy_data[idx(0, idx_t_left, idx_r_right)],
p_E_lasy_data[idx(0, idx_t_right, idx_r_left)],
p_E_lasy_data[idx(0, idx_t_right, idx_r_right)],
t, Rp_i);

// higher modes
for (int m=1 ; m <= tmp_n_rz_azimuthal_components/2; m++) {
val += utils::algorithms::bilinear_interp(
t_left, t_right,
r_0, r_1,
p_E_lasy_data[idx(2*m-1, idx_t_left, idx_r_left)],
p_E_lasy_data[idx(2*m-1, idx_t_left, idx_r_right)],
p_E_lasy_data[idx(2*m-1, idx_t_right, idx_r_left)],
p_E_lasy_data[idx(2*m-1, idx_t_right, idx_r_right)],
t, Rp_i)*(fact.real()) +
utils::algorithms::bilinear_interp(
t_left, t_right,
r_0, r_1,
p_E_lasy_data[idx(2*m, idx_t_left, idx_r_left)],
p_E_lasy_data[idx(2*m, idx_t_left, idx_r_right)],
p_E_lasy_data[idx(2*m, idx_t_right, idx_r_left)],
p_E_lasy_data[idx(2*m, idx_t_right, idx_r_right)],
t, Rp_i)*(fact.imag()) ;
val +=
ablastr::math::bilinear_interp(
t_left, t_right,
r_0, r_1,
p_E_lasy_data[idx(2*m-1, idx_t_left, idx_r_left)],
p_E_lasy_data[idx(2*m-1, idx_t_left, idx_r_right)],
p_E_lasy_data[idx(2*m-1, idx_t_right, idx_r_left)],
p_E_lasy_data[idx(2*m-1, idx_t_right, idx_r_right)],
t, Rp_i)*(fact.real()) +
ablastr::math::bilinear_interp(
t_left, t_right,
r_0, r_1,
p_E_lasy_data[idx(2*m, idx_t_left, idx_r_left)],
p_E_lasy_data[idx(2*m, idx_t_left, idx_r_right)],
p_E_lasy_data[idx(2*m, idx_t_right, idx_r_left)],
p_E_lasy_data[idx(2*m, idx_t_right, idx_r_right)],
t, Rp_i)*(fact.imag()) ;
fact = fact*Complex{costheta, sintheta};
}
amplitude[i] = (val*exp_omega_t).real();
Expand Down Expand Up @@ -683,7 +685,7 @@ WarpXLaserProfiles::FromFileLaserProfile::internal_fill_amplitude_uniform_binary
(i_interp-tmp_idx_first_time)*tmp_nx*tmp_ny+
j_interp*tmp_ny + k_interp;
};
amplitude[i] = utils::algorithms::trilinear_interp(
amplitude[i] = ablastr::math::trilinear_interp(
t_left, t_right,
x_0, x_1,
y_0, y_1,
Expand All @@ -702,7 +704,7 @@ WarpXLaserProfiles::FromFileLaserProfile::internal_fill_amplitude_uniform_binary
const auto idx = [=](int i_interp, int j_interp){
return (i_interp-tmp_idx_first_time) * tmp_nx + j_interp;
};
amplitude[i] = utils::algorithms::bilinear_interp(
amplitude[i] = ablastr::math::bilinear_interp(
t_left, t_right,
x_0, x_1,
p_E_binary_data[idx(idx_t_left, idx_x_left)],
Expand Down
19 changes: 11 additions & 8 deletions ...ce/Utils/Algorithms/LinearInterpolation.H → Source/ablastr/math/LinearInterpolation.H
View file
Original file line number Diff line number Diff line change
@@ -1,24 +1,25 @@
/* Copyright 2022 Luca Fedeli
/* Copyright 2022-2025 Luca Fedeli
*
* This file is part of WarpX.
*
* License: BSD-3-Clause-LBNL
*/

#ifndef WARPX_UTILS_ALGORITHMS_LINEAR_INTERPOLATION_H_
#define WARPX_UTILS_ALGORITHMS_LINEAR_INTERPOLATION_H_
#ifndef ABLASTR_MATH_LINEAR_INTERPOLATION_H_
#define ABLASTR_MATH_LINEAR_INTERPOLATION_H_

#include <AMReX_Extension.H>
#include <AMReX_GpuQualifiers.H>

namespace utils::algorithms
namespace ablastr::math
{
/** \brief Performs a linear interpolation
*
* Performs a linear interpolation at x given the 2 points
* (x0, f0) and (x1, f1)
*/
template<typename TCoord, typename TVal> AMREX_GPU_DEVICE AMREX_FORCE_INLINE
template<typename TCoord, typename TVal>
AMREX_GPU_DEVICE AMREX_FORCE_INLINE
constexpr auto linear_interp(
TCoord x0, TCoord x1,
TVal f0, TVal f1,
Expand All @@ -32,7 +33,8 @@ namespace utils::algorithms
* Performs a bilinear interpolation at (x,y) given the 4 points
* (x0, y0, f00), (x0, y1, f01), (x1, y0, f10), (x1, y1, f11).
*/
template<typename TCoord, typename TVal> AMREX_GPU_DEVICE AMREX_FORCE_INLINE
template<typename TCoord, typename TVal>
AMREX_GPU_DEVICE AMREX_FORCE_INLINE
constexpr auto bilinear_interp(
TCoord x0, TCoord x1, TCoord y0, TCoord y1,
TVal f00, TVal f01, TVal f10, TVal f11,
Expand All @@ -49,7 +51,8 @@ namespace utils::algorithms
* (x0, y0, z0, f000), (x0, y0, z1, f001), (x0, y1, z0, f010), (x0, y1, z1, f011),
* (x1, y0, z0, f100), (x1, y0, z1, f101), (x1, y1, z0, f110), (x1, y1, z1, f111)
*/
template<typename TCoord, typename TVal> AMREX_GPU_DEVICE AMREX_FORCE_INLINE
template<typename TCoord, typename TVal>
AMREX_GPU_DEVICE AMREX_FORCE_INLINE
constexpr auto trilinear_interp(
TCoord x0, TCoord x1, TCoord y0, TCoord y1, TCoord z0, TCoord z1,
TVal f000, TVal f001, TVal f010, TVal f011, TVal f100, TVal f101, TVal f110, TVal f111,
Expand All @@ -63,4 +66,4 @@ namespace utils::algorithms
}
}

#endif //WARPX_UTILS_ALGORITHMS_LINEAR_INTERPOLATION_H_
#endif //ABLASTR_MATH_LINEAR_INTERPOLATION_H_

0 comments on commit 2386c18

Please sign in to comment.