Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Move linear interpolation functions to ablastr #5714

Merged
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 3 additions & 3 deletions Source/Initialization/WarpXInitData.cpp
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
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
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_
Loading