diff --git a/src/pcms/utility/CMakeLists.txt b/src/pcms/utility/CMakeLists.txt
index 041069a8..0016de28 100644
--- a/src/pcms/utility/CMakeLists.txt
+++ b/src/pcms/utility/CMakeLists.txt
@@ -5,6 +5,7 @@ set(
         bounding_box.h
         common.h
         entity_types.h
+        eqdsk.h
         mesh_geometry.h
         memory_spaces.h
         mpi_type.h
@@ -20,6 +21,7 @@ set(
 set(
         PCMS_UTILITY_SOURCES
         assert.cpp
+        eqdsk.cpp
         print.cpp
 )
 add_library(pcms_utility ${PCMS_UTILITY_SOURCES})
@@ -33,6 +35,9 @@ target_include_directories(
         "$<BUILD_INTERFACE:${CMAKE_CURRENT_BINARY_DIR}/../..>"
         "$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>")
 target_link_libraries(pcms_utility PUBLIC Kokkos::kokkos perfstubs)
+if(PCMS_ENABLE_OMEGA_H)
+    target_link_libraries(pcms_utility PUBLIC Omega_h::omega_h)
+endif()
 target_compile_features(pcms_utility PUBLIC cxx_std_20)
 
 set_target_properties(
diff --git a/src/pcms/utility/eqdsk.cpp b/src/pcms/utility/eqdsk.cpp
new file mode 100644
index 00000000..0737f3c4
--- /dev/null
+++ b/src/pcms/utility/eqdsk.cpp
@@ -0,0 +1,309 @@
+#include "pcms/utility/eqdsk.h"
+#include "pcms/utility/assert.h"
+#include <fstream>
+#include <sstream>
+#include <iomanip>
+#include <cmath>
+#include <algorithm>
+#include <cctype>
+
+namespace pcms
+{
+
+// ============================================================================
+// File readers for EQDSKData
+// ============================================================================
+
+void read_geqdsk(const std::string& filename, EQDSKData& eqdsk)
+{
+  constexpr double psi_factor = 1.0; // Default psi factor
+
+  std::ifstream fin(filename);
+
+  if (!fin) {
+    throw std::runtime_error("Cannot open G-EQDSK file: " + filename);
+  }
+
+  // ----------------------------------------
+  // Read header
+  // ----------------------------------------
+  std::string first_line;
+  std::getline(fin, first_line);
+
+  if (first_line.length() < 60) {
+    first_line.resize(60, ' ');
+  }
+
+  int mw = std::stoi(first_line.substr(52, 4));
+  int mh = std::stoi(first_line.substr(56, 4));
+
+  // ----------------------------------------
+  // Read geometry
+  // ----------------------------------------
+  double rdim, zdim, rzero, r_min, zmid;
+  fin >> rdim >> zdim >> rzero >> r_min >> zmid;
+
+  double ssimag, ssibry, dummy;
+  fin >> dummy >> dummy >> ssimag >> ssibry >> dummy;
+  fin >> dummy >> dummy >> dummy >> dummy >> dummy;
+  fin >> dummy >> dummy >> dummy >> dummy >> dummy;
+
+  // ----------------------------------------
+  // Read fpol (poloidal current function)
+  // ----------------------------------------
+  std::vector<double> fpol(mw);
+  for (int i = 0; i < mw; i++) {
+    fin >> fpol[i];
+  }
+
+  // Skip pres, ffprim, pprime arrays (3 * mw values)
+  for (int i = 0; i < 3 * mw; i++) {
+    fin >> dummy;
+  }
+
+  // ----------------------------------------
+  // Read psirz
+  // ----------------------------------------
+  std::vector<std::vector<double>> psirz(mw, std::vector<double>(mh));
+  for (int j = 0; j < mh; j++) {
+    for (int i = 0; i < mw; i++) {
+      fin >> psirz[i][j];
+    }
+  }
+
+  // Skip qpsi array
+  for (int i = 0; i < mw; i++) {
+    fin >> dummy;
+  }
+
+  // Skip boundary and limiter info
+  int nbdry, nlim;
+  fin >> nbdry >> nlim;
+  for (int i = 0; i < nbdry + nlim; i++) {
+    fin >> dummy >> dummy;
+  }
+
+  fin.close();
+
+  // ----------------------------------------
+  // Build EQDSKData
+  // ----------------------------------------
+  ssimag *= psi_factor;
+  ssibry *= psi_factor;
+  double z_min = zmid - zdim / 2.0;
+
+  // Initialize grid structure
+  eqdsk.grid.divisions[0] = static_cast<LO>(mw);
+  eqdsk.grid.divisions[1] = static_cast<LO>(mh);
+  eqdsk.grid.edge_length[0] = static_cast<Real>(rdim);
+  eqdsk.grid.edge_length[1] = static_cast<Real>(zdim);
+  eqdsk.grid.bot_left[0] = static_cast<Real>(r_min);
+  eqdsk.grid.bot_left[1] = static_cast<Real>(z_min);
+
+  // Allocate and copy PSIZR
+  eqdsk.PSIZR = Kokkos::View<Real*, DeviceMemorySpace>(
+    "PSIZR", static_cast<size_t>(mw) * mh);
+  auto psizr_host = Kokkos::create_mirror_view(eqdsk.PSIZR);
+  for (int j = 0; j < mh; ++j) {
+    for (int i = 0; i < mw; ++i) {
+      psizr_host(i + j * mw) = static_cast<Real>(psirz[i][j] * psi_factor);
+    }
+  }
+  Kokkos::deep_copy(eqdsk.PSIZR, psizr_host);
+
+  // Build psi_grid
+  eqdsk.psi_grid = Kokkos::View<Real*, DeviceMemorySpace>("psi_grid", mw);
+  auto psi_grid_host = Kokkos::create_mirror_view(eqdsk.psi_grid);
+  for (int i = 0; i < mw; ++i) {
+    psi_grid_host(i) = static_cast<Real>(
+      double(i) / double(mw - 1) * (ssibry - ssimag) + ssimag);
+  }
+  Kokkos::deep_copy(eqdsk.psi_grid, psi_grid_host);
+
+  // Allocate and copy poloidal current function
+  eqdsk.I_psi = Kokkos::View<Real*, DeviceMemorySpace>("I_psi", mw);
+  auto I_psi_host = Kokkos::create_mirror_view(eqdsk.I_psi);
+  for (int i = 0; i < mw; ++i) {
+    I_psi_host(i) = static_cast<Real>(fpol[i]);
+  }
+  Kokkos::deep_copy(eqdsk.I_psi, I_psi_host);
+}
+
+void read_eqd(const std::string& filename, EQDSKData& eqdsk)
+{
+  std::ifstream fin(filename);
+
+  if (!fin) {
+    throw std::runtime_error("Cannot open EQD file: " + filename);
+  }
+
+  // Read header
+  std::string eq_header;
+  std::getline(fin, eq_header);
+
+  // Read dimensions
+  int mw, mh, mpsi;
+  fin >> mw >> mh >> mpsi;
+
+  // Read geometry
+  double r_min, r_max, z_min, z_max;
+  fin >> r_min >> r_max >> z_min >> z_max;
+
+  double rmaxis, zmaxis, eq_axis_b;
+  fin >> rmaxis >> zmaxis >> eq_axis_b;
+
+  // Skip dummy variables
+  double dummy;
+  fin >> dummy >> dummy >> dummy;
+
+  // Read psi_grid
+  std::vector<double> psi_grid_vec(mpsi);
+  for (int i = 0; i < mpsi; i++) {
+    fin >> psi_grid_vec[i];
+  }
+
+  std::cout << "Axis (R,Z,B) = " << rmaxis << " " << zmaxis << " " << eq_axis_b
+            << "\\n";
+
+  // Read I(psi)
+  std::vector<double> eq_I_vec(mpsi);
+  for (int i = 0; i < mpsi; i++) {
+    fin >> eq_I_vec[i];
+  }
+
+  // Read psi(R,Z)
+  std::vector<std::vector<double>> psirz(mw, std::vector<double>(mh));
+  for (int j = 0; j < mh; j++) {
+    for (int i = 0; i < mw; i++) {
+      fin >> psirz[i][j];
+    }
+  }
+
+  // Verify end flag
+  int end_flag;
+  fin >> end_flag;
+  if (end_flag != -1) {
+    throw std::runtime_error("Wrong EQD file format");
+  }
+
+  fin.close();
+
+  // ----------------------------------------
+  // Build EQDSKData
+  // ----------------------------------------
+
+  // Initialize grid structure
+  eqdsk.grid.divisions[0] = static_cast<LO>(mw);
+  eqdsk.grid.divisions[1] = static_cast<LO>(mh);
+  eqdsk.grid.edge_length[0] = static_cast<Real>(r_max - r_min);
+  eqdsk.grid.edge_length[1] = static_cast<Real>(z_max - z_min);
+  eqdsk.grid.bot_left[0] = static_cast<Real>(r_min);
+  eqdsk.grid.bot_left[1] = static_cast<Real>(z_min);
+
+  // Allocate and copy PSIZR
+  eqdsk.PSIZR = Kokkos::View<Real*, DeviceMemorySpace>(
+    "PSIZR", static_cast<size_t>(mw) * mh);
+  auto psizr_host = Kokkos::create_mirror_view(eqdsk.PSIZR);
+  for (int j = 0; j < mh; ++j) {
+    for (int i = 0; i < mw; ++i) {
+      psizr_host(i + j * mw) = static_cast<Real>(psirz[i][j]);
+    }
+  }
+  Kokkos::deep_copy(eqdsk.PSIZR, psizr_host);
+
+  // Allocate and copy psi_grid
+  eqdsk.psi_grid = Kokkos::View<Real*, DeviceMemorySpace>("psi_grid", mpsi);
+  auto psi_grid_host = Kokkos::create_mirror_view(eqdsk.psi_grid);
+  for (int i = 0; i < mpsi; ++i) {
+    psi_grid_host(i) = static_cast<Real>(psi_grid_vec[i]);
+  }
+  Kokkos::deep_copy(eqdsk.psi_grid, psi_grid_host);
+
+  // Allocate and copy poloidal current function
+  eqdsk.I_psi = Kokkos::View<Real*, DeviceMemorySpace>("I_psi", mpsi);
+  auto I_psi_host = Kokkos::create_mirror_view(eqdsk.I_psi);
+  for (int i = 0; i < mpsi; ++i) {
+    I_psi_host(i) = static_cast<Real>(eq_I_vec[i]);
+  }
+  Kokkos::deep_copy(eqdsk.I_psi, I_psi_host);
+
+  std::cout << "EQD file loaded successfully\n";
+}
+
+// ============================================================================
+// EQDSKData Constructor
+// ============================================================================
+EQDSKData::EQDSKData(const std::string& filename)
+{
+  // Determine file type by extension
+  std::string extension;
+  size_t dot_pos = filename.find_last_of('.');
+  if (dot_pos != std::string::npos) {
+    extension = filename.substr(dot_pos);
+    // Convert to lowercase for case-insensitive comparison
+    std::transform(extension.begin(), extension.end(), extension.begin(),
+                   [](unsigned char c) { return std::tolower(c); });
+  }
+
+  // Check for EQD format extensions
+  if (extension == ".eqd") {
+    read_eqd(filename, *this);
+  }
+  // Default to G-EQDSK format for .geq, .geqdsk, .g, or any other extension
+  else {
+    read_geqdsk(filename, *this);
+  }
+}
+
+// Manual constructor for programmatic construction
+EQDSKData::EQDSKData(const Uniform2DGrid& grid_in,
+                     const std::vector<Real>& psirz,
+                     const std::vector<Real>& psi_grid_vals,
+                     const std::vector<Real>& I_psi_vals)
+  : grid(grid_in)
+{
+  const int mw = grid.divisions[0];
+  const int mh = grid.divisions[1];
+
+  // Validate input dimensions
+  if (psirz.size() != static_cast<size_t>(mw) * mh) {
+    throw std::runtime_error("psirz size mismatch: expected " +
+                             std::to_string(mw * mh) + " but got " +
+                             std::to_string(psirz.size()));
+  }
+
+  if (psi_grid_vals.size() != I_psi_vals.size()) {
+    throw std::runtime_error("psi_grid and I_psi size mismatch: " +
+                             std::to_string(psi_grid_vals.size()) + " vs " +
+                             std::to_string(I_psi_vals.size()));
+  }
+
+  const int npsi = static_cast<int>(psi_grid_vals.size());
+
+  // Allocate and copy PSIZR
+  PSIZR = Kokkos::View<Real*, DeviceMemorySpace>("PSIZR",
+                                                 static_cast<size_t>(mw) * mh);
+  auto psizr_host = Kokkos::create_mirror_view(PSIZR);
+  for (int i = 0; i < mw * mh; ++i) {
+    psizr_host(i) = psirz[i];
+  }
+  Kokkos::deep_copy(PSIZR, psizr_host);
+
+  // Allocate and copy psi_grid
+  psi_grid = Kokkos::View<Real*, DeviceMemorySpace>("psi_grid", npsi);
+  auto psi_grid_host = Kokkos::create_mirror_view(psi_grid);
+  for (int i = 0; i < npsi; ++i) {
+    psi_grid_host(i) = psi_grid_vals[i];
+  }
+  Kokkos::deep_copy(psi_grid, psi_grid_host);
+
+  // Allocate and copy I_psi
+  I_psi = Kokkos::View<Real*, DeviceMemorySpace>("I_psi", npsi);
+  auto I_psi_host = Kokkos::create_mirror_view(I_psi);
+  for (int i = 0; i < npsi; ++i) {
+    I_psi_host(i) = I_psi_vals[i];
+  }
+  Kokkos::deep_copy(I_psi, I_psi_host);
+}
+
+} // namespace pcms
diff --git a/src/pcms/utility/eqdsk.h b/src/pcms/utility/eqdsk.h
new file mode 100644
index 00000000..aec34843
--- /dev/null
+++ b/src/pcms/utility/eqdsk.h
@@ -0,0 +1,133 @@
+#ifndef PCMS_UTILITY_EQDSK_H
+#define PCMS_UTILITY_EQDSK_H
+
+#include "pcms/utility/types.h"
+#include "pcms/utility/uniform_grid.h"
+#include "pcms/utility/memory_spaces.h"
+#include <vector>
+#include <string>
+#include <array>
+
+namespace pcms
+{
+
+/**
+ * @brief Data structure for EQDSK equilibrium data (simplified for field
+ * evaluation)
+ *
+ * Contains the essential data needed for poloidal flux field evaluation:
+ * - Uniform 2D grid structure (R-Z computational domain)
+ * - PSIZR: Poloidal flux values on grid points
+ * - psi_grid: Normalized poloidal flux grid values
+ * - I_psi: Poloidal current function at psi_grid points (F(psi) for G-EQDSK,
+ *   I(psi) for EQD format)
+ *
+ * All grid data is stored in column-major order (R varies fastest).
+ */
+struct EQDSKData
+{
+  // Grid structure for the R-Z computational domain
+  // grid.divisions[0] = NW (number of R grid points)
+  // grid.divisions[1] = NH (number of Z grid points)
+  // grid.edge_length[0] = RDIM (horizontal dimension in meter)
+  // grid.edge_length[1] = ZDIM (vertical dimension in meter)
+  // grid.bot_left[0] = RLEFT (minimum R in meter)
+  // grid.bot_left[1] = Z_min (minimum Z in meter)
+  Uniform2DGrid grid;
+
+  // PSIZR: Poloidal flux in Weber/rad on the rectangular grid points
+  Kokkos::View<Real*, DeviceMemorySpace> PSIZR;
+
+  // Psi grid: Normalized poloidal flux grid values
+  Kokkos::View<Real*, DeviceMemorySpace> psi_grid;
+
+  // I_psi: Poloidal current function I(psi) at psi_grid points
+  Kokkos::View<Real*, DeviceMemorySpace> I_psi;
+
+  // Constructors
+  // Constructor from file
+  EQDSKData(const std::string& filename);
+
+  // Manual constructor for programmatic construction
+  // grid: Uniform2DGrid structure defining the R-Z computational domain
+  // psirz: 2D flux array (flattened in column-major order, size =
+  // grid.divisions[0]*grid.divisions[1]) psi_grid_vals: psi grid values (size =
+  // npsi) I_psi_vals: poloidal current function values (size = npsi)
+  EQDSKData(const Uniform2DGrid& grid, const std::vector<Real>& psirz,
+            const std::vector<Real>& psi_grid_vals,
+            const std::vector<Real>& I_psi_vals);
+
+  // Accessor methods for backward compatibility
+  [[nodiscard]] int GetNW() const noexcept { return grid.divisions[0]; }
+  [[nodiscard]] int GetNH() const noexcept { return grid.divisions[1]; }
+  [[nodiscard]] Real GetRDIM() const noexcept { return grid.edge_length[0]; }
+  [[nodiscard]] Real GetZDIM() const noexcept { return grid.edge_length[1]; }
+  [[nodiscard]] Real GetRLEFT() const noexcept { return grid.bot_left[0]; }
+  [[nodiscard]] Real GetZMID() const noexcept
+  {
+    return grid.bot_left[1] + grid.edge_length[1] / 2.0;
+  }
+
+  // Grid bounds
+  [[nodiscard]] Real GetZMin() const noexcept { return grid.bot_left[1]; }
+
+  [[nodiscard]] Real GetZMax() const noexcept
+  {
+    return grid.bot_left[1] + grid.edge_length[1];
+  }
+
+  [[nodiscard]] Real GetRMin() const noexcept { return grid.bot_left[0]; }
+
+  [[nodiscard]] Real GetRMax() const noexcept
+  {
+    return grid.bot_left[0] + grid.edge_length[0];
+  }
+
+  [[nodiscard]] Real GetDeltaR() const noexcept
+  {
+    const int NW = grid.divisions[0];
+    return (NW > 1) ? grid.edge_length[0] / static_cast<Real>(NW - 1) : 0.0;
+  }
+
+  [[nodiscard]] Real GetDeltaZ() const noexcept
+  {
+    const int NH = grid.divisions[1];
+    return (NH > 1) ? grid.edge_length[1] / static_cast<Real>(NH - 1) : 0.0;
+  }
+
+  [[nodiscard]] int GetPsiIndex(int i_R, int i_Z) const noexcept
+  {
+    return i_R + i_Z * grid.divisions[0];
+  }
+
+  [[nodiscard]] Real GetPsi(int i_R, int i_Z) const
+  {
+    return PSIZR[GetPsiIndex(i_R, i_Z)];
+  }
+
+  [[nodiscard]] Real GetR(int i_R) const noexcept
+  {
+    return grid.bot_left[0] + i_R * GetDeltaR();
+  }
+
+  [[nodiscard]] Real GetZ(int i_Z) const noexcept
+  {
+    return grid.bot_left[1] + i_Z * GetDeltaZ();
+  }
+
+  // Accessor for poloidal current data size
+  [[nodiscard]] int GetNPsi() const noexcept
+  {
+    return static_cast<int>(psi_grid.extent(0));
+  }
+
+  // Accessor for psi grid value
+  [[nodiscard]] Real GetPsiGrid(int i) const { return psi_grid[i]; }
+
+  // Accessor for poloidal current at psi grid point
+  [[nodiscard]] Real GetI(int i) const { return I_psi[i]; }
+};
+
+} // namespace pcms
+
+#endif // PCMS_UTILITY_EQDSK_H
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 8412bac2..1c603c45 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -380,6 +380,7 @@ if(Catch2_FOUND)
       APPEND
       PCMS_UNIT_TEST_SOURCES
       test_error_handling.cpp
+      test_eqdsk.cpp
       test_uniform_grid.cpp
       test_field_evaluation.cpp
       test_field_interpolation.cpp
diff --git a/test/test_eqdsk.cpp b/test/test_eqdsk.cpp
new file mode 100644
index 00000000..1e3d164a
--- /dev/null
+++ b/test/test_eqdsk.cpp
@@ -0,0 +1,260 @@
+#include <pcms/utility/eqdsk.h>
+#include <pcms/utility/arrays.h>
+#include <pcms/field/function_space/spline.h>
+#include <pcms/field/coordinate_system.h>
+#include <pcms/field/evaluation_request.h>
+#include <Kokkos_Core.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/catch_approx.hpp>
+#include <iostream>
+#include <cmath>
+
+using Catch::Approx;
+using pcms::CoordinateSystem;
+using pcms::EQDSKData;
+
+TEST_CASE("EQDSKData manual constructor")
+{
+  auto lib = Omega_h::Library{};
+
+  SECTION("Basic grid construction")
+  {
+    // Create a simple test grid
+    pcms::Uniform2DGrid test_grid;
+    test_grid.divisions[0] = 5;     // 5 R points
+    test_grid.divisions[1] = 4;     // 4 Z points
+    test_grid.bot_left[0] = 1.0;    // R_min
+    test_grid.bot_left[1] = -1.0;   // Z_min
+    test_grid.edge_length[0] = 2.0; // R_max - R_min = 3.0 - 1.0
+    test_grid.edge_length[1] = 2.0; // Z_max - Z_min = 1.0 - (-1.0)
+
+    // Create test PSIZR data (5x4 = 20 values)
+    std::vector<pcms::Real> test_psirz(20);
+    for (int i = 0; i < 20; ++i) {
+      test_psirz[i] = static_cast<pcms::Real>(i * 0.1);
+    }
+
+    // Create test psi_grid and I_psi data
+    const int test_npsi = 5;
+    std::vector<pcms::Real> test_psi_grid(test_npsi);
+    std::vector<pcms::Real> test_I_psi(test_npsi);
+    for (int i = 0; i < test_npsi; ++i) {
+      test_psi_grid[i] = static_cast<pcms::Real>(i * 0.2);
+      test_I_psi[i] = static_cast<pcms::Real>(1.0 + i * 0.5);
+    }
+
+    // Construct EQDSKData manually
+    EQDSKData eqdsk_data(test_grid, test_psirz, test_psi_grid, test_I_psi);
+
+    // Verify grid properties
+    REQUIRE(eqdsk_data.GetNW() == 5);
+    REQUIRE(eqdsk_data.GetNH() == 4);
+    REQUIRE(eqdsk_data.GetRMin() == Approx(1.0));
+    REQUIRE(eqdsk_data.GetRMax() == Approx(3.0));
+    REQUIRE(eqdsk_data.GetZMin() == Approx(-1.0));
+    REQUIRE(eqdsk_data.GetZMax() == Approx(1.0));
+    REQUIRE(eqdsk_data.GetRDIM() == Approx(2.0));
+    REQUIRE(eqdsk_data.GetZDIM() == Approx(2.0));
+
+    // Verify grid spacing
+    REQUIRE(eqdsk_data.GetDeltaR() == Approx(0.5));
+    REQUIRE(eqdsk_data.GetDeltaZ() == Approx(2.0 / 3.0));
+
+    // Verify PSIZR data
+    REQUIRE(eqdsk_data.PSIZR.extent(0) == 20);
+    auto psizr_host = Kokkos::create_mirror_view_and_copy(
+      pcms::HostMemorySpace(), eqdsk_data.PSIZR);
+    for (int i = 0; i < 20; ++i) {
+      REQUIRE(psizr_host(i) == Approx(test_psirz[i]));
+    }
+
+    // Verify psi_grid and I_psi data
+    REQUIRE(eqdsk_data.GetNPsi() == test_npsi);
+    auto psi_grid_host = Kokkos::create_mirror_view_and_copy(
+      pcms::HostMemorySpace(), eqdsk_data.psi_grid);
+    auto I_psi_host = Kokkos::create_mirror_view_and_copy(
+      pcms::HostMemorySpace(), eqdsk_data.I_psi);
+
+    for (int i = 0; i < test_npsi; ++i) {
+      REQUIRE(psi_grid_host(i) == Approx(test_psi_grid[i]));
+      REQUIRE(I_psi_host(i) == Approx(test_I_psi[i]));
+    }
+  }
+
+  SECTION("Grid accessor methods")
+  {
+    // Create a test grid
+    pcms::Uniform2DGrid test_grid;
+    test_grid.divisions[0] = 5;
+    test_grid.divisions[1] = 4;
+    test_grid.bot_left[0] = 1.0;
+    test_grid.bot_left[1] = -1.0;
+    test_grid.edge_length[0] = 2.0;
+    test_grid.edge_length[1] = 2.0;
+
+    std::vector<pcms::Real> test_psirz(20, 0.0);
+    std::vector<pcms::Real> test_psi_grid(5, 0.0);
+    std::vector<pcms::Real> test_I_psi(5, 0.0);
+
+    EQDSKData eqdsk_data(test_grid, test_psirz, test_psi_grid, test_I_psi);
+
+    // Verify index calculations
+    REQUIRE(eqdsk_data.GetPsiIndex(0, 0) == 0);
+    REQUIRE(eqdsk_data.GetPsiIndex(4, 0) == 4);
+    REQUIRE(eqdsk_data.GetPsiIndex(0, 3) == 15);
+    REQUIRE(eqdsk_data.GetPsiIndex(4, 3) == 19);
+
+    // Verify coordinate calculations
+    REQUIRE(eqdsk_data.GetR(0) == Approx(1.0));
+    REQUIRE(eqdsk_data.GetR(4) == Approx(3.0));
+    REQUIRE(eqdsk_data.GetZ(0) == Approx(-1.0));
+    REQUIRE(eqdsk_data.GetZ(3) == Approx(1.0));
+  }
+
+  SECTION("Input validation")
+  {
+    pcms::Uniform2DGrid test_grid;
+    test_grid.divisions[0] = 5;
+    test_grid.divisions[1] = 4;
+    test_grid.bot_left[0] = 1.0;
+    test_grid.bot_left[1] = -1.0;
+    test_grid.edge_length[0] = 2.0;
+    test_grid.edge_length[1] = 2.0;
+
+    std::vector<pcms::Real> test_psi_grid(5, 0.0);
+    std::vector<pcms::Real> test_I_psi(5, 0.0);
+
+    // Wrong PSIZR size (should be 20)
+    std::vector<pcms::Real> wrong_psirz(15, 0.0);
+    REQUIRE_THROWS_AS(
+      EQDSKData(test_grid, wrong_psirz, test_psi_grid, test_I_psi),
+      std::runtime_error);
+
+    // Mismatched psi_grid and I_psi sizes
+    std::vector<pcms::Real> correct_psirz(20, 0.0);
+    std::vector<pcms::Real> wrong_I_psi(3, 0.0);
+    REQUIRE_THROWS_AS(
+      EQDSKData(test_grid, correct_psirz, test_psi_grid, wrong_I_psi),
+      std::runtime_error);
+  }
+}
+
+TEST_CASE("EQDSKData with SplineFunctionSpace")
+{
+  auto lib = Omega_h::Library{};
+
+  // Create a test grid
+  pcms::Uniform2DGrid test_grid;
+  test_grid.divisions[0] = 10; // 10 R points
+  test_grid.divisions[1] = 8;  // 8 Z points
+  test_grid.bot_left[0] = 1.0;
+  test_grid.bot_left[1] = -1.0;
+  test_grid.edge_length[0] = 2.0;
+  test_grid.edge_length[1] = 2.0;
+
+  const int NW = 10;
+  const int NH = 8;
+
+  // Create test PSIZR data with a simple quadratic pattern
+  std::vector<pcms::Real> test_psirz(static_cast<size_t>(NW) * NH);
+  for (int j = 0; j < NH; ++j) {
+    for (int i = 0; i < NW; ++i) {
+      double r = 1.0 + i * 2.0 / (NW - 1);
+      double z = -1.0 + j * 2.0 / (NH - 1);
+      test_psirz[i + j * NW] = static_cast<pcms::Real>(r * r + z * z);
+    }
+  }
+
+  // Create test psi_grid and I_psi data
+  const int test_npsi = 10;
+  std::vector<pcms::Real> test_psi_grid(test_npsi);
+  std::vector<pcms::Real> test_I_psi(test_npsi);
+  for (int i = 0; i < test_npsi; ++i) {
+    test_psi_grid[i] = static_cast<pcms::Real>(i * 0.5);
+    test_I_psi[i] = static_cast<pcms::Real>(1.0 + i * 0.3);
+  }
+
+  EQDSKData eqdsk_data(test_grid, test_psirz, test_psi_grid, test_I_psi);
+
+  SECTION("Spline function space creation and evaluation")
+  {
+    // Create a modified grid for spline interpolation
+    // EQDSK grid.divisions = NW x NH (number of data points)
+    // SplineFunctionSpace with order=1 needs grid.divisions = (NW-1) x (NH-1)
+    // cells
+    pcms::UniformGrid<2> spline_grid;
+    spline_grid.bot_left = eqdsk_data.grid.bot_left;
+    spline_grid.edge_length = eqdsk_data.grid.edge_length;
+    spline_grid.divisions = {NW - 1, NH - 1}; // Number of cells, not points
+
+    auto spline_space = pcms::SplineFunctionSpace::FromUniformGrid(
+      spline_grid, CoordinateSystem::Cartesian);
+
+    auto psi_field = spline_space.CreateField<pcms::Real>();
+
+    auto layout = spline_space.GetLayout();
+    const int layout_size = layout->OwnedSize();
+    const int psizr_size = static_cast<int>(eqdsk_data.PSIZR.extent(0));
+
+    REQUIRE(layout_size == psizr_size);
+
+    psi_field.GetData().SetDOFHolderData(
+      pcms::Rank1View<const pcms::Real, pcms::DeviceMemorySpace>(
+        eqdsk_data.PSIZR.data(), eqdsk_data.PSIZR.extent(0)));
+
+    // Evaluate at test points
+    const pcms::Real R_mid =
+      (eqdsk_data.GetRMin() + eqdsk_data.GetRMax()) / 2.0;
+    const pcms::Real Z_mid =
+      (eqdsk_data.GetZMin() + eqdsk_data.GetZMax()) / 2.0;
+
+    std::vector<pcms::Real> eval_coords = {
+      R_mid,
+      Z_mid, // center
+      eqdsk_data.GetRMin() + 0.25 * eqdsk_data.GetRDIM(),
+      eqdsk_data.GetZMin() + 0.25 * eqdsk_data.GetZDIM(), // lower left
+      eqdsk_data.GetRMin() + 0.75 * eqdsk_data.GetRDIM(),
+      eqdsk_data.GetZMin() + 0.75 * eqdsk_data.GetZDIM() // upper right
+    };
+
+    const int num_eval_points = 3;
+
+    auto eval_coords_host = Kokkos::View<pcms::Real**, pcms::HostMemorySpace>(
+      "eval_coords_host", num_eval_points, 2);
+    for (size_t i = 0; i < num_eval_points; ++i) {
+      eval_coords_host(i, 0) = eval_coords[2 * i];     // R
+      eval_coords_host(i, 1) = eval_coords[2 * i + 1]; // Z
+    }
+
+    auto eval_coords_device =
+      Kokkos::View<pcms::Real**, pcms::DeviceMemorySpace>("eval_coords_device",
+                                                          num_eval_points, 2);
+    pcms::DeepCopyMismatchLayouts(eval_coords_device, eval_coords_host);
+
+    auto coords_view = pcms::MakeRank2View(eval_coords_device);
+    auto coord_view = pcms::CoordinateView<pcms::DeviceMemorySpace>{
+      CoordinateSystem::Cartesian, coords_view};
+    auto eval_request = pcms::EvaluationRequest::FromCoordinates(coord_view);
+
+    auto evaluator =
+      spline_space.CreatePointEvaluator<pcms::Real>(eval_request);
+
+    auto eval_results_1d = Kokkos::View<pcms::Real*, pcms::DeviceMemorySpace>(
+      "eval_results", num_eval_points);
+
+    using LayoutPolicy =
+      pcms::detail::default_layout_for_memory_space_t<pcms::DeviceMemorySpace>;
+    pcms::Rank2View<pcms::Real, pcms::DeviceMemorySpace, LayoutPolicy>
+      eval_results(eval_results_1d.data(), num_eval_points, 1);
+
+    evaluator->Evaluate(psi_field, eval_results);
+
+    auto results_host = Kokkos::create_mirror_view_and_copy(
+      pcms::HostMemorySpace(), eval_results_1d);
+
+    // Verify that all results are finite
+    for (int i = 0; i < num_eval_points; ++i) {
+      REQUIRE(std::isfinite(results_host(i)));
+    }
+  }
+}