diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp
index c70188f07bc..3089fd7304c 100644
--- a/Source/Initialization/WarpXInitData.cpp
+++ b/Source/Initialization/WarpXInitData.cpp
@@ -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"
@@ -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>
@@ -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));
@@ -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));
diff --git a/Source/Laser/LaserProfilesImpl/LaserProfileFromFile.cpp b/Source/Laser/LaserProfilesImpl/LaserProfileFromFile.cpp
index 934a537be5a..0ce5e71337f 100644
--- a/Source/Laser/LaserProfilesImpl/LaserProfileFromFile.cpp
+++ b/Source/Laser/LaserProfilesImpl/LaserProfileFromFile.cpp
@@ -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>
@@ -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,
@@ -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();
@@ -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,
@@ -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)],
diff --git a/Source/Utils/Algorithms/LinearInterpolation.H b/Source/ablastr/math/LinearInterpolation.H
similarity index 79%
rename from Source/Utils/Algorithms/LinearInterpolation.H
rename to Source/ablastr/math/LinearInterpolation.H
index 716c2d829b6..502ef68c080 100644
--- a/Source/Utils/Algorithms/LinearInterpolation.H
+++ b/Source/ablastr/math/LinearInterpolation.H
@@ -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,
@@ -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,
@@ -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,
@@ -63,4 +66,4 @@ namespace utils::algorithms
     }
 }
 
-#endif //WARPX_UTILS_ALGORITHMS_LINEAR_INTERPOLATION_H_
+#endif //ABLASTR_MATH_LINEAR_INTERPOLATION_H_