From 5b643e7236af81186435cd7cad8ff50cbb595b87 Mon Sep 17 00:00:00 2001 From: Sichao25 Date: Thu, 14 May 2026 13:46:59 -0400 Subject: [PATCH 1/4] support eqdsk field --- src/pcms/utility/CMakeLists.txt | 5 + src/pcms/utility/eqdsk.cpp | 441 ++++++++++++++++++++++++++++++++ src/pcms/utility/eqdsk.h | 197 ++++++++++++++ test/CMakeLists.txt | 11 + test/test_eqdsk.cpp | 180 +++++++++++++ 5 files changed, 834 insertions(+) create mode 100644 src/pcms/utility/eqdsk.cpp create mode 100644 src/pcms/utility/eqdsk.h create mode 100644 test/test_eqdsk.cpp 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( "$" "$") 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..ad8b61ec --- /dev/null +++ b/src/pcms/utility/eqdsk.cpp @@ -0,0 +1,441 @@ +#include "pcms/utility/eqdsk.h" +#include "pcms/utility/assert.h" +#include +#include +#include +#include +#include +#include + +namespace pcms +{ + +void readgfile(const std::string& filename, GEQDData& eqd) +{ + constexpr double pi = 3.1415926535897932; + constexpr double tmu = 2.0e-07; + constexpr double twopi = 2.0 * pi; + + std::ifstream fin(filename); + + if (!fin) { + throw std::runtime_error("Cannot open file"); + } + + // ---------------------------------------- + // Read header + // ---------------------------------------- + // G-EQDSK format (Fortran format: 6a8,3i4): + // 6 strings of 8 characters each (characters 1-48) + // 3 integers of 4 characters each (characters 49-60) + // Example: " EFITD 05/12/98 # 96333 ,3337ms 3 65 65" + + std::string first_line; + std::getline(fin, first_line); + + // Ensure line is long enough (at least 60 characters for full format) + if (first_line.length() < 60) { + first_line.resize(60, ' '); + } + + // Read 6 strings of 8 characters each (positions 0-47) + for (int i = 0; i < 6; i++) { + eqd.eqd_case[i] = first_line.substr(i * 8, 8); + } + + // Read 3 integers of 4 characters each (positions 48-59) + std::string imfit_str = first_line.substr(48, 4); + std::string mw_str = first_line.substr(52, 4); + std::string mh_str = first_line.substr(56, 4); + + eqd.eqd_imfit = std::stoi(imfit_str); + eqd.eqd_mw = std::stoi(mw_str); + eqd.eqd_mh = std::stoi(mh_str); + + // ---------------------------------------- + // Allocate arrays + // ---------------------------------------- + + int mw = eqd.eqd_mw; + int mh = eqd.eqd_mh; + + eqd.eqd_fpol.resize(mw); + eqd.eqd_pres.resize(mw); + eqd.eqd_qpsi.resize(mw); + eqd.eqd_workk.resize(mw); + eqd.eqd_ffprim.resize(mw); + eqd.eqd_pprime.resize(mw); + eqd.eqd_xsi.resize(mw); + + eqd.eqd_rgrid.resize(mw); + eqd.eqd_zgrid.resize(mh); + eqd.eqd_psi_grid.resize(mw); + + eqd.eqd_psirz.resize(mw, std::vector(mh)); + + // ---------------------------------------- + // Read geometry + // ---------------------------------------- + + fin >> eqd.eqd_rdim >> eqd.eqd_zdim >> eqd.eqd_rzero >> eqd.eqd_r_min >> + eqd.eqd_zmid; + + fin >> eqd.eqd_rmaxis >> eqd.eqd_zmaxis >> eqd.eqd_ssimag >> eqd.eqd_ssibry >> + eqd.eqd_bzero; + + double xdum; + double sdum; + + fin >> eqd.eqd_cpasma >> xdum >> xdum >> eqd.eqd_rmaxis >> xdum; + + fin >> eqd.eqd_zmaxis >> xdum >> sdum >> xdum >> xdum; + + // ---------------------------------------- + // Read profiles + // ---------------------------------------- + + for (int i = 0; i < mw; i++) { + fin >> eqd.eqd_fpol[i]; + } + + for (int i = 0; i < mw; i++) { + fin >> eqd.eqd_pres[i]; + } + + for (int i = 0; i < mw; i++) { + fin >> eqd.eqd_workk[i]; + } + + eqd.eqd_drgrid = eqd.eqd_rdim / double(mw - 1); + eqd.eqd_dzgrid = eqd.eqd_zdim / double(mh - 1); + + eqd.eqd_z_min = eqd.eqd_zmid - eqd.eqd_zdim / 2.0; + + eqd.eqd_darea = eqd.eqd_drgrid * eqd.eqd_dzgrid; + + for (int i = 0; i < mw; i++) { + if (eqd.eqd_imfit >= 0) { + eqd.eqd_ffprim[i] = -eqd.eqd_workk[i] / (twopi * tmu); + } else { + eqd.eqd_ffprim[i] = -eqd.eqd_workk[i]; + } + } + + // ---------------------------------------- + // pprime + // ---------------------------------------- + + for (int i = 0; i < mw; i++) { + fin >> eqd.eqd_workk[i]; + } + + for (int i = 0; i < mw; i++) { + eqd.eqd_pprime[i] = -eqd.eqd_workk[i]; + } + + // ---------------------------------------- + // Read psirz + // ---------------------------------------- + + for (int j = 0; j < mh; j++) { + for (int i = 0; i < mw; i++) { + fin >> eqd.eqd_psirz[i][j]; + } + } + + // ---------------------------------------- + // qpsi + // ---------------------------------------- + + for (int i = 0; i < mw; i++) { + fin >> eqd.eqd_qpsi[i]; + } + + // ---------------------------------------- + // Boundary info + // ---------------------------------------- + + fin >> eqd.eqd_nbdry >> eqd.eqd_nlim; + + eqd.eqd_rbdry.resize(eqd.eqd_nbdry); + eqd.eqd_zbdry.resize(eqd.eqd_nbdry); + + eqd.eqd_rlim.resize(eqd.eqd_nlim); + eqd.eqd_zlim.resize(eqd.eqd_nlim); + + for (int i = 0; i < eqd.eqd_nbdry; i++) { + fin >> eqd.eqd_rbdry[i] >> eqd.eqd_zbdry[i]; + } + + for (int i = 0; i < eqd.eqd_nlim; i++) { + fin >> eqd.eqd_rlim[i] >> eqd.eqd_zlim[i]; + } + + fin.close(); + + // ---------------------------------------- + // Apply psi factor + // ---------------------------------------- + + for (int i = 0; i < mw; i++) { + for (int j = 0; j < mh; j++) { + eqd.eqd_psirz[i][j] *= eqd.eqd_psi_factor; + } + } + + eqd.eqd_ssimag *= eqd.eqd_psi_factor; + eqd.eqd_ssibry *= eqd.eqd_psi_factor; + + // ---------------------------------------- + // Build grids + // ---------------------------------------- + + eqd.eqd_r_max = eqd.eqd_r_min + eqd.eqd_rdim; + eqd.eqd_z_max = eqd.eqd_z_min + eqd.eqd_zdim; + + for (int i = 0; i < mw; i++) { + eqd.eqd_rgrid[i] = eqd.eqd_r_min + (eqd.eqd_r_max - eqd.eqd_r_min) * + double(i) / double(mw - 1); + } + + for (int i = 0; i < mh; i++) { + eqd.eqd_zgrid[i] = eqd.eqd_z_min + (eqd.eqd_z_max - eqd.eqd_z_min) * + double(i) / double(mh - 1); + } + + for (int i = 0; i < mw; i++) { + eqd.eqd_psi_grid[i] = + double(i) / double(mw - 1) * (eqd.eqd_ssibry - eqd.eqd_ssimag) + + eqd.eqd_ssimag; + } + + std::cout << "Finished reading G-EQDSK file\n"; +} + +void readeqdfile(const std::string& filename, EQDData& eqd) +{ + std::ifstream fin(filename); + + if (!fin) { + throw std::runtime_error("Cannot open eqd file"); + } + + // ---------------------------------------- + // Read header + // ---------------------------------------- + + std::string eq_header; + std::getline(fin, eq_header); + + // ---------------------------------------- + // Read dimensions + // ---------------------------------------- + + fin >> eqd.eqd_mw >> eqd.eqd_mh >> eqd.eqd_mpsi; + + int mw = eqd.eqd_mw; + int mh = eqd.eqd_mh; + int mpsi = eqd.eqd_mpsi; + + // ---------------------------------------- + // Allocate arrays + // ---------------------------------------- + + eqd.eqd_psi_grid.resize(mpsi); + + eqd.eq_I.resize(mpsi); + + eqd.eqd_rgrid.resize(mw); + eqd.eqd_zgrid.resize(mh); + + eqd.eqd_psirz.resize(mw, std::vector(mh)); + + // ---------------------------------------- + // Read geometry + // ---------------------------------------- + + fin >> eqd.eqd_r_min >> eqd.eqd_r_max >> eqd.eqd_z_min >> eqd.eqd_z_max; + + double eq_axis_b; + + fin >> eqd.eqd_rmaxis >> eqd.eqd_zmaxis >> eq_axis_b; + + // dummy variables + double eq_x_psi_loc; + double eq_x_r; + double eq_x_z; + + fin >> eq_x_psi_loc >> eq_x_r >> eq_x_z; + + // ---------------------------------------- + // Read psi grid + // ---------------------------------------- + + for (int i = 0; i < mpsi; i++) { + fin >> eqd.eqd_psi_grid[i]; + } + + std::cout << "Axis (R,Z,B) = " << eqd.eqd_rmaxis << " " << eqd.eqd_zmaxis + << " " << eq_axis_b << "\n"; + + // ---------------------------------------- + // Read I(psi) + // ---------------------------------------- + + for (int i = 0; i < mpsi; i++) { + fin >> eqd.eq_I[i]; + } + + // ---------------------------------------- + // Read psi(R,Z) + // ---------------------------------------- + + // Keep same indexing convention + // as original Fortran: + // + // eqd_psirz(i,j) + // + for (int j = 0; j < mh; j++) { + for (int i = 0; i < mw; i++) { + fin >> eqd.eqd_psirz[i][j]; + } + } + + // ---------------------------------------- + // End flag + // ---------------------------------------- + + int end_flag; + + fin >> end_flag; + + if (end_flag != -1) { + throw std::runtime_error("Wrong EQD file format"); + } + + fin.close(); + + // ---------------------------------------- + // Build R grid + // ---------------------------------------- + + for (int i = 0; i < mw; i++) { + + eqd.eqd_rgrid[i] = eqd.eqd_r_min + (eqd.eqd_r_max - eqd.eqd_r_min) * + double(i) / double(mw - 1); + } + + // ---------------------------------------- + // Build Z grid + // ---------------------------------------- + + for (int i = 0; i < mh; i++) { + + eqd.eqd_zgrid[i] = eqd.eqd_z_min + (eqd.eqd_z_max - eqd.eqd_z_min) * + double(i) / double(mh - 1); + } + + // ---------------------------------------- + // Construct rectangular limiter + // ---------------------------------------- + + eqd.eqd_nlim = 4; + + eqd.eqd_rlim.resize(4); + eqd.eqd_zlim.resize(4); + + eqd.eqd_rlim[0] = eqd.eqd_rgrid[0]; + eqd.eqd_zlim[0] = eqd.eqd_zgrid[0]; + + eqd.eqd_rlim[1] = eqd.eqd_rgrid[0]; + eqd.eqd_zlim[1] = eqd.eqd_zgrid[mh - 1]; + + eqd.eqd_rlim[2] = eqd.eqd_rgrid[mw - 1]; + eqd.eqd_zlim[2] = eqd.eqd_zgrid[mh - 1]; + + eqd.eqd_rlim[3] = eqd.eqd_rgrid[mw - 1]; + eqd.eqd_zlim[3] = eqd.eqd_zgrid[0]; + + std::cout << "EQD file loaded successfully\n"; +} + +// ============================================================================ +// EQDSKData Constructors +// ============================================================================ + +// Constructor from filename +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") { + EQDData eqd; + readeqdfile(filename, eqd); + *this = EQDSKData(eqd); + } + // Default to G-EQDSK format for .geq, .geqdsk, .g, or any other extension + else { + GEQDData geqd; + readgfile(filename, geqd); + *this = EQDSKData(geqd); + } +} + +// Constructor from GEQDData (G-EQDSK format) +EQDSKData::EQDSKData(const GEQDData& geqd) +{ + // Initialize grid structure + grid.divisions[0] = static_cast(geqd.eqd_mw); + grid.divisions[1] = static_cast(geqd.eqd_mh); + grid.edge_length[0] = static_cast(geqd.eqd_rdim); + grid.edge_length[1] = static_cast(geqd.eqd_zdim); + grid.bot_left[0] = static_cast(geqd.eqd_r_min); + grid.bot_left[1] = static_cast(geqd.eqd_z_min); + + // Allocate and copy 2D PSIZR array (column-major: R varies fastest) + const int nw = geqd.eqd_mw; + const int nh = geqd.eqd_mh; + PSIZR = Kokkos::View("PSIZR", nw * nh); + auto psizr_host = Kokkos::create_mirror_view(PSIZR); + for (int j = 0; j < nh; ++j) { + for (int i = 0; i < nw; ++i) { + psizr_host(i + j * nw) = static_cast(geqd.eqd_psirz[i][j]); + } + } + Kokkos::deep_copy(PSIZR, psizr_host); +} + +// Constructor from EQDData (EQD format) +EQDSKData::EQDSKData(const EQDData& eqd) +{ + // Initialize grid structure + grid.divisions[0] = static_cast(eqd.eqd_mw); + grid.divisions[1] = static_cast(eqd.eqd_mh); + grid.edge_length[0] = static_cast(eqd.eqd_r_max - eqd.eqd_r_min); + grid.edge_length[1] = static_cast(eqd.eqd_z_max - eqd.eqd_z_min); + grid.bot_left[0] = static_cast(eqd.eqd_r_min); + grid.bot_left[1] = static_cast(eqd.eqd_z_min); + + // Allocate and copy 2D PSIZR array (column-major: R varies fastest) + const int nw = eqd.eqd_mw; + const int nh = eqd.eqd_mh; + PSIZR = Kokkos::View("PSIZR", nw * nh); + auto psizr_host = Kokkos::create_mirror_view(PSIZR); + for (int j = 0; j < nh; ++j) { + for (int i = 0; i < nw; ++i) { + psizr_host(i + j * nw) = static_cast(eqd.eqd_psirz[i][j]); + } + } + Kokkos::deep_copy(PSIZR, psizr_host); +} + +} // namespace pcms diff --git a/src/pcms/utility/eqdsk.h b/src/pcms/utility/eqdsk.h new file mode 100644 index 00000000..a04a6b76 --- /dev/null +++ b/src/pcms/utility/eqdsk.h @@ -0,0 +1,197 @@ +#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 +#include +#include + +namespace pcms +{ + +struct GEQDData +{ + // Header + std::array eqd_case; + + int eqd_imfit; + int eqd_mw; + int eqd_mh; + + // Scalars + double eqd_rdim; + double eqd_zdim; + double eqd_rzero; + double eqd_r_min; + double eqd_zmid; + + double eqd_rmaxis; + double eqd_zmaxis; + double eqd_ssimag; + double eqd_ssibry; + double eqd_bzero; + + double eqd_cpasma; + + double eqd_drgrid; + double eqd_dzgrid; + double eqd_z_min; + double eqd_r_max; + double eqd_z_max; + double eqd_darea; + + int eqd_nbdry; + int eqd_nlim; + + double eqd_psi_factor = 1.0; + + // 1D arrays + std::vector eqd_fpol; + std::vector eqd_pres; + std::vector eqd_qpsi; + std::vector eqd_workk; + std::vector eqd_ffprim; + std::vector eqd_pprime; + std::vector eqd_xsi; + + std::vector eqd_rgrid; + std::vector eqd_zgrid; + std::vector eqd_psi_grid; + + std::vector eqd_rbdry; + std::vector eqd_zbdry; + + std::vector eqd_rlim; + std::vector eqd_zlim; + + // 2D array + std::vector> eqd_psirz; +}; + +struct EQDData +{ + int eqd_mw; + int eqd_mh; + int eqd_mpsi; + + int eqd_nlim; + + double eqd_r_min; + double eqd_r_max; + double eqd_z_min; + double eqd_z_max; + + double eqd_rmaxis; + double eqd_zmaxis; + + std::vector eqd_psi_grid; + std::vector eq_I; + + std::vector eqd_rgrid; + std::vector eqd_zgrid; + + std::vector eqd_rlim; + std::vector eqd_zlim; + + // 2D psi array + std::vector> eqd_psirz; +}; + +/** + * @brief Data structure for EQDSK equilibrium data (simplified for field + * evaluation) + * + * Contains only the essential data needed for poloidal flux field evaluation: + * - Uniform 2D grid structure (R-Z computational domain) + * - PSIZR: Poloidal flux values on grid points + * + * 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 PSIZR; + + // Constructors + // Constructor from file + EQDSKData(const std::string& filename); + + // Constructor from GEQDData (G-EQDSK format) + explicit EQDSKData(const GEQDData& geqd); + + // Constructor from EQDData (EQD format) + explicit EQDSKData(const EQDData& eqd); + + // 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(NW - 1) : 0.0; + } + + [[nodiscard]] Real GetDeltaZ() const noexcept + { + const int NH = grid.divisions[1]; + return (NH > 1) ? grid.edge_length[1] / static_cast(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(); + } +}; + +} // namespace pcms + +#endif // PCMS_UTILITY_EQDSK_H diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 8412bac2..4f8578dc 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -34,6 +34,17 @@ if(PCMS_ENABLE_OMEGA_H) ${notRendezvous} ${cyclone1p}) endif() + + # EQDSK test - takes filename as parameter + add_exe(test_eqdsk) + target_link_libraries(test_eqdsk pcms::utility) + # To run: ctest -R test_eqdsk_load -E none + # or: mpiexec -n 1 ./test_eqdsk /path/to/your/file.geqdsk + # Note: Test is not added by default - provide EQDSK file path to enable + if(DEFINED PCMS_EQDSK_TEST_FILE) + mpi_test(test_eqdsk_load 1 ./test_eqdsk ${PCMS_EQDSK_TEST_FILE}) + endif() + add_exe(test_ohClassPtn_appRibPtn) set(cyclone2p ${PCMS_TEST_DATA_DIR}/cyclone/23elements/2p.osh/) if(HOST_NPROC GREATER_EQUAL 4) diff --git a/test/test_eqdsk.cpp b/test/test_eqdsk.cpp new file mode 100644 index 00000000..cf759841 --- /dev/null +++ b/test/test_eqdsk.cpp @@ -0,0 +1,180 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using pcms::CoordinateSystem; +using pcms::EQDSKData; + +// This test only verifies that the EQDSK data can be loaded and that the +// SplineFunctionSpace can be created and evaluated. The actual values of the +// loaded data and the spline evaluation are not checked. +int main(int argc, char** argv) +{ + auto lib = Omega_h::Library{}; + + if (argc != 2) { + std::cerr << "Usage: " << argv[0] << " \n"; + return 1; + } + + const std::string filename(argv[1]); + + std::cout << "Loading EQDSK file: " << filename << "\n"; + + try { + // Test constructor from file + EQDSKData eqdsk_data(filename); + + std::cout << "Grid dimensions: NW = " << eqdsk_data.GetNW() + << ", NH = " << eqdsk_data.GetNH() << "\n"; + std::cout << "R range: [" << eqdsk_data.GetRMin() << ", " + << eqdsk_data.GetRMax() << "]\n"; + std::cout << "Z range: [" << eqdsk_data.GetZMin() << ", " + << eqdsk_data.GetZMax() << "]\n"; + std::cout << "Grid spacing: dR = " << eqdsk_data.GetDeltaR() + << ", dZ = " << eqdsk_data.GetDeltaZ() << "\n"; + + // Verify grid dimensions are positive + assert(eqdsk_data.GetNW() > 0); + assert(eqdsk_data.GetNH() > 0); + + // Verify grid bounds are valid + assert(eqdsk_data.GetRMin() < eqdsk_data.GetRMax()); + assert(eqdsk_data.GetZMin() < eqdsk_data.GetZMax()); + + // Verify grid dimensions match + const double r_tol = 1e-10; + const double z_tol = 1e-10; + assert(std::abs(eqdsk_data.GetRDIM() - + (eqdsk_data.GetRMax() - eqdsk_data.GetRMin())) < r_tol); + assert(std::abs(eqdsk_data.GetZDIM() - + (eqdsk_data.GetZMax() - eqdsk_data.GetZMin())) < z_tol); + + // Verify PSIZR array size + const int expected_size = eqdsk_data.GetNW() * eqdsk_data.GetNH(); + assert(eqdsk_data.PSIZR.extent(0) == static_cast(expected_size)); + + assert(eqdsk_data.GetDeltaR() > 0.0); + assert(eqdsk_data.GetDeltaZ() > 0.0); + + const int NW = eqdsk_data.GetNW(); + const int NH = eqdsk_data.GetNH(); + + assert(eqdsk_data.GetPsiIndex(0, 0) == 0); + assert(eqdsk_data.GetPsiIndex(NW - 1, 0) == NW - 1); + assert(eqdsk_data.GetPsiIndex(0, NH - 1) == (NH - 1) * NW); + assert(eqdsk_data.GetPsiIndex(NW - 1, NH - 1) == NW * NH - 1); + + const double coord_tol = 1e-10; + assert(std::abs(eqdsk_data.GetR(0) - eqdsk_data.GetRMin()) < coord_tol); + assert(std::abs(eqdsk_data.GetR(NW - 1) - eqdsk_data.GetRMax()) < + coord_tol); + assert(std::abs(eqdsk_data.GetZ(0) - eqdsk_data.GetZMin()) < coord_tol); + assert(std::abs(eqdsk_data.GetZ(NH - 1) - eqdsk_data.GetZMax()) < + coord_tol); + + std::cout << "\nTesting SplineFunctionSpace evaluation...\n"; + + // 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(); + + auto layout = spline_space.GetLayout(); + const int layout_size = layout->OwnedSize(); + const int psizr_size = static_cast(eqdsk_data.PSIZR.extent(0)); + std::cout << "Layout OwnedSize: " << layout_size << "\n"; + std::cout << "PSIZR size: " << psizr_size << " (NW=" << NW << " * NH=" << NH + << ")\n"; + + assert(layout_size == psizr_size); + + psi_field.GetData().SetDOFHolderData( + pcms::Rank1View( + eqdsk_data.PSIZR.data(), eqdsk_data.PSIZR.extent(0))); + + 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 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( + "eval_coords_host", num_eval_points, 2); + for (int 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("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{ + CoordinateSystem::Cartesian, coords_view}; + auto eval_request = pcms::EvaluationRequest::FromCoordinates(coord_view); + + auto evaluator = + spline_space.CreatePointEvaluator(eval_request); + + auto eval_results_1d = Kokkos::View( + "eval_results", num_eval_points); + + using LayoutPolicy = + pcms::detail::default_layout_for_memory_space_t; + pcms::Rank2View + 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); + + std::cout << "Spline evaluation results:\n"; + for (int i = 0; i < num_eval_points; ++i) { + std::cout << " Point " << i << ": (R=" << eval_coords[2 * i] + << ", Z=" << eval_coords[2 * i + 1] + << ") -> Psi = " << results_host(i) << "\n"; + } + + for (int i = 0; i < num_eval_points; ++i) { + assert(std::isfinite(results_host(i))); + } + + std::cout << "SplineFunctionSpace evaluation test passed!\n"; + + std::cout << "\nAll tests passed!\n"; + return 0; + + } catch (const std::exception& e) { + std::cerr << "Error: " << e.what() << "\n"; + return 1; + } +} From 9b4590bf77539869f2270df542a28ac212985c31 Mon Sep 17 00:00:00 2001 From: Sichao25 Date: Fri, 15 May 2026 22:21:07 -0400 Subject: [PATCH 2/4] simplify eqdsk input --- src/pcms/utility/eqdsk.cpp | 453 +++++++++++++------------------------ src/pcms/utility/eqdsk.h | 124 +++------- test/CMakeLists.txt | 23 +- test/test_eqdsk.cpp | 236 ++++++++++++------- 4 files changed, 360 insertions(+), 476 deletions(-) diff --git a/src/pcms/utility/eqdsk.cpp b/src/pcms/utility/eqdsk.cpp index ad8b61ec..2ff787d6 100644 --- a/src/pcms/utility/eqdsk.cpp +++ b/src/pcms/utility/eqdsk.cpp @@ -10,305 +10,177 @@ namespace pcms { -void readgfile(const std::string& filename, GEQDData& eqd) +// ============================================================================ +// File readers for EQDSKData +// ============================================================================ + +void read_geqdsk(const std::string& filename, EQDSKData& eqdsk) { - constexpr double pi = 3.1415926535897932; - constexpr double tmu = 2.0e-07; - constexpr double twopi = 2.0 * pi; + constexpr double psi_factor = 1.0; // Default psi factor std::ifstream fin(filename); if (!fin) { - throw std::runtime_error("Cannot open file"); + throw std::runtime_error("Cannot open G-EQDSK file: " + filename); } // ---------------------------------------- // Read header // ---------------------------------------- - // G-EQDSK format (Fortran format: 6a8,3i4): - // 6 strings of 8 characters each (characters 1-48) - // 3 integers of 4 characters each (characters 49-60) - // Example: " EFITD 05/12/98 # 96333 ,3337ms 3 65 65" - std::string first_line; std::getline(fin, first_line); - // Ensure line is long enough (at least 60 characters for full format) if (first_line.length() < 60) { first_line.resize(60, ' '); } - // Read 6 strings of 8 characters each (positions 0-47) - for (int i = 0; i < 6; i++) { - eqd.eqd_case[i] = first_line.substr(i * 8, 8); - } - - // Read 3 integers of 4 characters each (positions 48-59) - std::string imfit_str = first_line.substr(48, 4); - std::string mw_str = first_line.substr(52, 4); - std::string mh_str = first_line.substr(56, 4); - - eqd.eqd_imfit = std::stoi(imfit_str); - eqd.eqd_mw = std::stoi(mw_str); - eqd.eqd_mh = std::stoi(mh_str); - - // ---------------------------------------- - // Allocate arrays - // ---------------------------------------- - - int mw = eqd.eqd_mw; - int mh = eqd.eqd_mh; - - eqd.eqd_fpol.resize(mw); - eqd.eqd_pres.resize(mw); - eqd.eqd_qpsi.resize(mw); - eqd.eqd_workk.resize(mw); - eqd.eqd_ffprim.resize(mw); - eqd.eqd_pprime.resize(mw); - eqd.eqd_xsi.resize(mw); - - eqd.eqd_rgrid.resize(mw); - eqd.eqd_zgrid.resize(mh); - eqd.eqd_psi_grid.resize(mw); - - eqd.eqd_psirz.resize(mw, std::vector(mh)); + 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; - fin >> eqd.eqd_rdim >> eqd.eqd_zdim >> eqd.eqd_rzero >> eqd.eqd_r_min >> - eqd.eqd_zmid; - - fin >> eqd.eqd_rmaxis >> eqd.eqd_zmaxis >> eqd.eqd_ssimag >> eqd.eqd_ssibry >> - eqd.eqd_bzero; - - double xdum; - double sdum; - - fin >> eqd.eqd_cpasma >> xdum >> xdum >> eqd.eqd_rmaxis >> xdum; - - fin >> eqd.eqd_zmaxis >> xdum >> sdum >> xdum >> xdum; + double ssimag, ssibry, dummy; + fin >> dummy >> dummy >> ssimag >> ssibry >> dummy; + fin >> dummy >> dummy >> dummy >> dummy >> dummy; + fin >> dummy >> dummy >> dummy >> dummy >> dummy; // ---------------------------------------- - // Read profiles + // Read fpol (poloidal current function) // ---------------------------------------- - - for (int i = 0; i < mw; i++) { - fin >> eqd.eqd_fpol[i]; - } - - for (int i = 0; i < mw; i++) { - fin >> eqd.eqd_pres[i]; - } - + std::vector fpol(mw); for (int i = 0; i < mw; i++) { - fin >> eqd.eqd_workk[i]; + fin >> fpol[i]; } - eqd.eqd_drgrid = eqd.eqd_rdim / double(mw - 1); - eqd.eqd_dzgrid = eqd.eqd_zdim / double(mh - 1); - - eqd.eqd_z_min = eqd.eqd_zmid - eqd.eqd_zdim / 2.0; - - eqd.eqd_darea = eqd.eqd_drgrid * eqd.eqd_dzgrid; - - for (int i = 0; i < mw; i++) { - if (eqd.eqd_imfit >= 0) { - eqd.eqd_ffprim[i] = -eqd.eqd_workk[i] / (twopi * tmu); - } else { - eqd.eqd_ffprim[i] = -eqd.eqd_workk[i]; - } - } - - // ---------------------------------------- - // pprime - // ---------------------------------------- - - for (int i = 0; i < mw; i++) { - fin >> eqd.eqd_workk[i]; - } - - for (int i = 0; i < mw; i++) { - eqd.eqd_pprime[i] = -eqd.eqd_workk[i]; + // Skip pres, ffprim, pprime arrays (3 * mw values) + for (int i = 0; i < 3 * mw; i++) { + fin >> dummy; } // ---------------------------------------- // Read psirz // ---------------------------------------- - + std::vector> psirz(mw, std::vector(mh)); for (int j = 0; j < mh; j++) { for (int i = 0; i < mw; i++) { - fin >> eqd.eqd_psirz[i][j]; + fin >> psirz[i][j]; } } - // ---------------------------------------- - // qpsi - // ---------------------------------------- - + // Skip qpsi array for (int i = 0; i < mw; i++) { - fin >> eqd.eqd_qpsi[i]; - } - - // ---------------------------------------- - // Boundary info - // ---------------------------------------- - - fin >> eqd.eqd_nbdry >> eqd.eqd_nlim; - - eqd.eqd_rbdry.resize(eqd.eqd_nbdry); - eqd.eqd_zbdry.resize(eqd.eqd_nbdry); - - eqd.eqd_rlim.resize(eqd.eqd_nlim); - eqd.eqd_zlim.resize(eqd.eqd_nlim); - - for (int i = 0; i < eqd.eqd_nbdry; i++) { - fin >> eqd.eqd_rbdry[i] >> eqd.eqd_zbdry[i]; + fin >> dummy; } - for (int i = 0; i < eqd.eqd_nlim; i++) { - fin >> eqd.eqd_rlim[i] >> eqd.eqd_zlim[i]; + // Skip boundary and limiter info + int nbdry, nlim; + fin >> nbdry >> nlim; + for (int i = 0; i < nbdry + nlim; i++) { + fin >> dummy >> dummy; } fin.close(); // ---------------------------------------- - // Apply psi factor + // Build EQDSKData // ---------------------------------------- + ssimag *= psi_factor; + ssibry *= psi_factor; + double z_min = zmid - zdim / 2.0; - for (int i = 0; i < mw; i++) { - for (int j = 0; j < mh; j++) { - eqd.eqd_psirz[i][j] *= eqd.eqd_psi_factor; + // Initialize grid structure + eqdsk.grid.divisions[0] = static_cast(mw); + eqdsk.grid.divisions[1] = static_cast(mh); + eqdsk.grid.edge_length[0] = static_cast(rdim); + eqdsk.grid.edge_length[1] = static_cast(zdim); + eqdsk.grid.bot_left[0] = static_cast(r_min); + eqdsk.grid.bot_left[1] = static_cast(z_min); + + // Allocate and copy PSIZR + eqdsk.PSIZR = Kokkos::View("PSIZR", 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(psirz[i][j] * psi_factor); } } - - eqd.eqd_ssimag *= eqd.eqd_psi_factor; - eqd.eqd_ssibry *= eqd.eqd_psi_factor; - - // ---------------------------------------- - // Build grids - // ---------------------------------------- - - eqd.eqd_r_max = eqd.eqd_r_min + eqd.eqd_rdim; - eqd.eqd_z_max = eqd.eqd_z_min + eqd.eqd_zdim; - - for (int i = 0; i < mw; i++) { - eqd.eqd_rgrid[i] = eqd.eqd_r_min + (eqd.eqd_r_max - eqd.eqd_r_min) * - double(i) / double(mw - 1); - } - - for (int i = 0; i < mh; i++) { - eqd.eqd_zgrid[i] = eqd.eqd_z_min + (eqd.eqd_z_max - eqd.eqd_z_min) * - double(i) / double(mh - 1); + Kokkos::deep_copy(eqdsk.PSIZR, psizr_host); + + // Build psi_grid + eqdsk.psi_grid = Kokkos::View("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( + double(i) / double(mw - 1) * (ssibry - ssimag) + ssimag); } + Kokkos::deep_copy(eqdsk.psi_grid, psi_grid_host); - for (int i = 0; i < mw; i++) { - eqd.eqd_psi_grid[i] = - double(i) / double(mw - 1) * (eqd.eqd_ssibry - eqd.eqd_ssimag) + - eqd.eqd_ssimag; + // Allocate and copy poloidal current function + eqdsk.I_psi = Kokkos::View("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(fpol[i]); } - - std::cout << "Finished reading G-EQDSK file\n"; + Kokkos::deep_copy(eqdsk.I_psi, I_psi_host); } -void readeqdfile(const std::string& filename, EQDData& eqd) +void read_eqd(const std::string& filename, EQDSKData& eqdsk) { std::ifstream fin(filename); if (!fin) { - throw std::runtime_error("Cannot open eqd file"); + throw std::runtime_error("Cannot open EQD file: " + filename); } - // ---------------------------------------- // Read header - // ---------------------------------------- - std::string eq_header; std::getline(fin, eq_header); - // ---------------------------------------- // Read dimensions - // ---------------------------------------- - - fin >> eqd.eqd_mw >> eqd.eqd_mh >> eqd.eqd_mpsi; - - int mw = eqd.eqd_mw; - int mh = eqd.eqd_mh; - int mpsi = eqd.eqd_mpsi; - - // ---------------------------------------- - // Allocate arrays - // ---------------------------------------- + int mw, mh, mpsi; + fin >> mw >> mh >> mpsi; - eqd.eqd_psi_grid.resize(mpsi); - - eqd.eq_I.resize(mpsi); - - eqd.eqd_rgrid.resize(mw); - eqd.eqd_zgrid.resize(mh); - - eqd.eqd_psirz.resize(mw, std::vector(mh)); - - // ---------------------------------------- // Read geometry - // ---------------------------------------- - - fin >> eqd.eqd_r_min >> eqd.eqd_r_max >> eqd.eqd_z_min >> eqd.eqd_z_max; - - double eq_axis_b; + double r_min, r_max, z_min, z_max; + fin >> r_min >> r_max >> z_min >> z_max; - fin >> eqd.eqd_rmaxis >> eqd.eqd_zmaxis >> eq_axis_b; + double rmaxis, zmaxis, eq_axis_b; + fin >> rmaxis >> zmaxis >> eq_axis_b; - // dummy variables - double eq_x_psi_loc; - double eq_x_r; - double eq_x_z; - - fin >> eq_x_psi_loc >> eq_x_r >> eq_x_z; - - // ---------------------------------------- - // Read psi grid - // ---------------------------------------- + // Skip dummy variables + double dummy; + fin >> dummy >> dummy >> dummy; + // Read psi_grid + std::vector psi_grid_vec(mpsi); for (int i = 0; i < mpsi; i++) { - fin >> eqd.eqd_psi_grid[i]; + fin >> psi_grid_vec[i]; } - std::cout << "Axis (R,Z,B) = " << eqd.eqd_rmaxis << " " << eqd.eqd_zmaxis - << " " << eq_axis_b << "\n"; + std::cout << "Axis (R,Z,B) = " << rmaxis << " " << zmaxis << " " << eq_axis_b + << "\\n"; - // ---------------------------------------- // Read I(psi) - // ---------------------------------------- - + std::vector eq_I_vec(mpsi); for (int i = 0; i < mpsi; i++) { - fin >> eqd.eq_I[i]; + fin >> eq_I_vec[i]; } - // ---------------------------------------- // Read psi(R,Z) - // ---------------------------------------- - - // Keep same indexing convention - // as original Fortran: - // - // eqd_psirz(i,j) - // + std::vector> psirz(mw, std::vector(mh)); for (int j = 0; j < mh; j++) { for (int i = 0; i < mw; i++) { - fin >> eqd.eqd_psirz[i][j]; + fin >> psirz[i][j]; } } - // ---------------------------------------- - // End flag - // ---------------------------------------- - + // Verify end flag int end_flag; - fin >> end_flag; - if (end_flag != -1) { throw std::runtime_error("Wrong EQD file format"); } @@ -316,54 +188,49 @@ void readeqdfile(const std::string& filename, EQDData& eqd) fin.close(); // ---------------------------------------- - // Build R grid + // Build EQDSKData // ---------------------------------------- - for (int i = 0; i < mw; i++) { - - eqd.eqd_rgrid[i] = eqd.eqd_r_min + (eqd.eqd_r_max - eqd.eqd_r_min) * - double(i) / double(mw - 1); + // Initialize grid structure + eqdsk.grid.divisions[0] = static_cast(mw); + eqdsk.grid.divisions[1] = static_cast(mh); + eqdsk.grid.edge_length[0] = static_cast(r_max - r_min); + eqdsk.grid.edge_length[1] = static_cast(z_max - z_min); + eqdsk.grid.bot_left[0] = static_cast(r_min); + eqdsk.grid.bot_left[1] = static_cast(z_min); + + // Allocate and copy PSIZR + eqdsk.PSIZR = Kokkos::View("PSIZR", 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(psirz[i][j]); + } } + Kokkos::deep_copy(eqdsk.PSIZR, psizr_host); - // ---------------------------------------- - // Build Z grid - // ---------------------------------------- - - for (int i = 0; i < mh; i++) { - - eqd.eqd_zgrid[i] = eqd.eqd_z_min + (eqd.eqd_z_max - eqd.eqd_z_min) * - double(i) / double(mh - 1); + // Allocate and copy psi_grid + eqdsk.psi_grid = Kokkos::View("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(psi_grid_vec[i]); } + Kokkos::deep_copy(eqdsk.psi_grid, psi_grid_host); - // ---------------------------------------- - // Construct rectangular limiter - // ---------------------------------------- - - eqd.eqd_nlim = 4; - - eqd.eqd_rlim.resize(4); - eqd.eqd_zlim.resize(4); - - eqd.eqd_rlim[0] = eqd.eqd_rgrid[0]; - eqd.eqd_zlim[0] = eqd.eqd_zgrid[0]; - - eqd.eqd_rlim[1] = eqd.eqd_rgrid[0]; - eqd.eqd_zlim[1] = eqd.eqd_zgrid[mh - 1]; - - eqd.eqd_rlim[2] = eqd.eqd_rgrid[mw - 1]; - eqd.eqd_zlim[2] = eqd.eqd_zgrid[mh - 1]; - - eqd.eqd_rlim[3] = eqd.eqd_rgrid[mw - 1]; - eqd.eqd_zlim[3] = eqd.eqd_zgrid[0]; + // Allocate and copy poloidal current function + eqdsk.I_psi = Kokkos::View("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(eq_I_vec[i]); + } + Kokkos::deep_copy(eqdsk.I_psi, I_psi_host); std::cout << "EQD file loaded successfully\n"; } // ============================================================================ -// EQDSKData Constructors +// EQDSKData Constructor // ============================================================================ - -// Constructor from filename EQDSKData::EQDSKData(const std::string& filename) { // Determine file type by extension @@ -378,64 +245,62 @@ EQDSKData::EQDSKData(const std::string& filename) // Check for EQD format extensions if (extension == ".eqd") { - EQDData eqd; - readeqdfile(filename, eqd); - *this = EQDSKData(eqd); + read_eqd(filename, *this); } // Default to G-EQDSK format for .geq, .geqdsk, .g, or any other extension else { - GEQDData geqd; - readgfile(filename, geqd); - *this = EQDSKData(geqd); + read_geqdsk(filename, *this); } } -// Constructor from GEQDData (G-EQDSK format) -EQDSKData::EQDSKData(const GEQDData& geqd) +// Manual constructor for programmatic construction +EQDSKData::EQDSKData(const Uniform2DGrid& grid_in, + const std::vector& psirz, + const std::vector& psi_grid_vals, + const std::vector& I_psi_vals) + : grid(grid_in) { - // Initialize grid structure - grid.divisions[0] = static_cast(geqd.eqd_mw); - grid.divisions[1] = static_cast(geqd.eqd_mh); - grid.edge_length[0] = static_cast(geqd.eqd_rdim); - grid.edge_length[1] = static_cast(geqd.eqd_zdim); - grid.bot_left[0] = static_cast(geqd.eqd_r_min); - grid.bot_left[1] = static_cast(geqd.eqd_z_min); - - // Allocate and copy 2D PSIZR array (column-major: R varies fastest) - const int nw = geqd.eqd_mw; - const int nh = geqd.eqd_mh; - PSIZR = Kokkos::View("PSIZR", nw * nh); - auto psizr_host = Kokkos::create_mirror_view(PSIZR); - for (int j = 0; j < nh; ++j) { - for (int i = 0; i < nw; ++i) { - psizr_host(i + j * nw) = static_cast(geqd.eqd_psirz[i][j]); - } + const int mw = grid.divisions[0]; + const int mh = grid.divisions[1]; + + // Validate input dimensions + if (psirz.size() != static_cast(mw * mh)) { + throw std::runtime_error("psirz size mismatch: expected " + + std::to_string(mw * mh) + " but got " + + std::to_string(psirz.size())); } - Kokkos::deep_copy(PSIZR, psizr_host); -} -// Constructor from EQDData (EQD format) -EQDSKData::EQDSKData(const EQDData& eqd) -{ - // Initialize grid structure - grid.divisions[0] = static_cast(eqd.eqd_mw); - grid.divisions[1] = static_cast(eqd.eqd_mh); - grid.edge_length[0] = static_cast(eqd.eqd_r_max - eqd.eqd_r_min); - grid.edge_length[1] = static_cast(eqd.eqd_z_max - eqd.eqd_z_min); - grid.bot_left[0] = static_cast(eqd.eqd_r_min); - grid.bot_left[1] = static_cast(eqd.eqd_z_min); - - // Allocate and copy 2D PSIZR array (column-major: R varies fastest) - const int nw = eqd.eqd_mw; - const int nh = eqd.eqd_mh; - PSIZR = Kokkos::View("PSIZR", nw * nh); + 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(psi_grid_vals.size()); + + // Allocate and copy PSIZR + PSIZR = Kokkos::View("PSIZR", mw * mh); auto psizr_host = Kokkos::create_mirror_view(PSIZR); - for (int j = 0; j < nh; ++j) { - for (int i = 0; i < nw; ++i) { - psizr_host(i + j * nw) = static_cast(eqd.eqd_psirz[i][j]); - } + 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("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("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 index a04a6b76..aec34843 100644 --- a/src/pcms/utility/eqdsk.h +++ b/src/pcms/utility/eqdsk.h @@ -11,101 +11,16 @@ namespace pcms { -struct GEQDData -{ - // Header - std::array eqd_case; - - int eqd_imfit; - int eqd_mw; - int eqd_mh; - - // Scalars - double eqd_rdim; - double eqd_zdim; - double eqd_rzero; - double eqd_r_min; - double eqd_zmid; - - double eqd_rmaxis; - double eqd_zmaxis; - double eqd_ssimag; - double eqd_ssibry; - double eqd_bzero; - - double eqd_cpasma; - - double eqd_drgrid; - double eqd_dzgrid; - double eqd_z_min; - double eqd_r_max; - double eqd_z_max; - double eqd_darea; - - int eqd_nbdry; - int eqd_nlim; - - double eqd_psi_factor = 1.0; - - // 1D arrays - std::vector eqd_fpol; - std::vector eqd_pres; - std::vector eqd_qpsi; - std::vector eqd_workk; - std::vector eqd_ffprim; - std::vector eqd_pprime; - std::vector eqd_xsi; - - std::vector eqd_rgrid; - std::vector eqd_zgrid; - std::vector eqd_psi_grid; - - std::vector eqd_rbdry; - std::vector eqd_zbdry; - - std::vector eqd_rlim; - std::vector eqd_zlim; - - // 2D array - std::vector> eqd_psirz; -}; - -struct EQDData -{ - int eqd_mw; - int eqd_mh; - int eqd_mpsi; - - int eqd_nlim; - - double eqd_r_min; - double eqd_r_max; - double eqd_z_min; - double eqd_z_max; - - double eqd_rmaxis; - double eqd_zmaxis; - - std::vector eqd_psi_grid; - std::vector eq_I; - - std::vector eqd_rgrid; - std::vector eqd_zgrid; - - std::vector eqd_rlim; - std::vector eqd_zlim; - - // 2D psi array - std::vector> eqd_psirz; -}; - /** * @brief Data structure for EQDSK equilibrium data (simplified for field * evaluation) * - * Contains only the essential data needed for poloidal flux 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). */ @@ -123,15 +38,24 @@ struct EQDSKData // PSIZR: Poloidal flux in Weber/rad on the rectangular grid points Kokkos::View PSIZR; + // Psi grid: Normalized poloidal flux grid values + Kokkos::View psi_grid; + + // I_psi: Poloidal current function I(psi) at psi_grid points + Kokkos::View I_psi; + // Constructors // Constructor from file EQDSKData(const std::string& filename); - // Constructor from GEQDData (G-EQDSK format) - explicit EQDSKData(const GEQDData& geqd); - - // Constructor from EQDData (EQD format) - explicit EQDSKData(const EQDData& eqd); + // 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& psirz, + const std::vector& psi_grid_vals, + const std::vector& I_psi_vals); // Accessor methods for backward compatibility [[nodiscard]] int GetNW() const noexcept { return grid.divisions[0]; } @@ -190,6 +114,18 @@ struct EQDSKData { return grid.bot_left[1] + i_Z * GetDeltaZ(); } + + // Accessor for poloidal current data size + [[nodiscard]] int GetNPsi() const noexcept + { + return static_cast(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 diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 4f8578dc..74e072d8 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -35,16 +35,6 @@ if(PCMS_ENABLE_OMEGA_H) ${cyclone1p}) endif() - # EQDSK test - takes filename as parameter - add_exe(test_eqdsk) - target_link_libraries(test_eqdsk pcms::utility) - # To run: ctest -R test_eqdsk_load -E none - # or: mpiexec -n 1 ./test_eqdsk /path/to/your/file.geqdsk - # Note: Test is not added by default - provide EQDSK file path to enable - if(DEFINED PCMS_EQDSK_TEST_FILE) - mpi_test(test_eqdsk_load 1 ./test_eqdsk ${PCMS_EQDSK_TEST_FILE}) - endif() - add_exe(test_ohClassPtn_appRibPtn) set(cyclone2p ${PCMS_TEST_DATA_DIR}/cyclone/23elements/2p.osh/) if(HOST_NPROC GREATER_EQUAL 4) @@ -391,6 +381,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 @@ -457,6 +448,18 @@ if(Catch2_FOUND) --data_root ${PCMS_TEST_DATA_DIR}/ltx-interpolation-case/ ) + # EQDSK test with Catch2 + #if(PCMS_ENABLE_OMEGA_H) + #add_executable(test_eqdsk test_eqdsk.cpp) + #target_link_libraries(test_eqdsk PUBLIC + # Catch2::Catch2WithMain + # pcms::utility + # pcms::core + # ) + # target_include_directories(test_eqdsk PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) + # Catch_discover_tests(test_eqdsk) + #endif() + include(Catch) Catch_discover_tests(unit_tests) else() diff --git a/test/test_eqdsk.cpp b/test/test_eqdsk.cpp index cf759841..0c9c3874 100644 --- a/test/test_eqdsk.cpp +++ b/test/test_eqdsk.cpp @@ -4,83 +4,180 @@ #include #include #include +#include +#include #include #include -#include +using Catch::Approx; using pcms::CoordinateSystem; using pcms::EQDSKData; -// This test only verifies that the EQDSK data can be loaded and that the -// SplineFunctionSpace can be created and evaluated. The actual values of the -// loaded data and the spline evaluation are not checked. -int main(int argc, char** argv) +TEST_CASE("EQDSKData manual constructor") { auto lib = Omega_h::Library{}; - if (argc != 2) { - std::cerr << "Usage: " << argv[0] << " \n"; - return 1; - } - - const std::string filename(argv[1]); - - std::cout << "Loading EQDSK file: " << filename << "\n"; - - try { - // Test constructor from file - EQDSKData eqdsk_data(filename); + 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 test_psirz(20); + for (int i = 0; i < 20; ++i) { + test_psirz[i] = static_cast(i * 0.1); + } - std::cout << "Grid dimensions: NW = " << eqdsk_data.GetNW() - << ", NH = " << eqdsk_data.GetNH() << "\n"; - std::cout << "R range: [" << eqdsk_data.GetRMin() << ", " - << eqdsk_data.GetRMax() << "]\n"; - std::cout << "Z range: [" << eqdsk_data.GetZMin() << ", " - << eqdsk_data.GetZMax() << "]\n"; - std::cout << "Grid spacing: dR = " << eqdsk_data.GetDeltaR() - << ", dZ = " << eqdsk_data.GetDeltaZ() << "\n"; + // Create test psi_grid and I_psi data + const int test_npsi = 5; + std::vector test_psi_grid(test_npsi); + std::vector test_I_psi(test_npsi); + for (int i = 0; i < test_npsi; ++i) { + test_psi_grid[i] = static_cast(i * 0.2); + test_I_psi[i] = static_cast(1.0 + i * 0.5); + } - // Verify grid dimensions are positive - assert(eqdsk_data.GetNW() > 0); - assert(eqdsk_data.GetNH() > 0); + // 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 grid bounds are valid - assert(eqdsk_data.GetRMin() < eqdsk_data.GetRMax()); - assert(eqdsk_data.GetZMin() < eqdsk_data.GetZMax()); + // 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); - // Verify grid dimensions match - const double r_tol = 1e-10; - const double z_tol = 1e-10; - assert(std::abs(eqdsk_data.GetRDIM() - - (eqdsk_data.GetRMax() - eqdsk_data.GetRMin())) < r_tol); - assert(std::abs(eqdsk_data.GetZDIM() - - (eqdsk_data.GetZMax() - eqdsk_data.GetZMin())) < z_tol); + 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])); + } + } - // Verify PSIZR array size - const int expected_size = eqdsk_data.GetNW() * eqdsk_data.GetNH(); - assert(eqdsk_data.PSIZR.extent(0) == static_cast(expected_size)); + 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 test_psirz(20, 0.0); + std::vector test_psi_grid(5, 0.0); + std::vector 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)); + } - assert(eqdsk_data.GetDeltaR() > 0.0); - assert(eqdsk_data.GetDeltaZ() > 0.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 test_psi_grid(5, 0.0); + std::vector test_I_psi(5, 0.0); + + // Wrong PSIZR size (should be 20) + std::vector 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 correct_psirz(20, 0.0); + std::vector wrong_I_psi(3, 0.0); + REQUIRE_THROWS_AS( + EQDSKData(test_grid, correct_psirz, test_psi_grid, wrong_I_psi), + std::runtime_error); + } +} - const int NW = eqdsk_data.GetNW(); - const int NH = eqdsk_data.GetNH(); +TEST_CASE("EQDSKData with SplineFunctionSpace") +{ + auto lib = Omega_h::Library{}; - assert(eqdsk_data.GetPsiIndex(0, 0) == 0); - assert(eqdsk_data.GetPsiIndex(NW - 1, 0) == NW - 1); - assert(eqdsk_data.GetPsiIndex(0, NH - 1) == (NH - 1) * NW); - assert(eqdsk_data.GetPsiIndex(NW - 1, NH - 1) == NW * NH - 1); + // 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 test_psirz(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(r * r + z * z); + } + } - const double coord_tol = 1e-10; - assert(std::abs(eqdsk_data.GetR(0) - eqdsk_data.GetRMin()) < coord_tol); - assert(std::abs(eqdsk_data.GetR(NW - 1) - eqdsk_data.GetRMax()) < - coord_tol); - assert(std::abs(eqdsk_data.GetZ(0) - eqdsk_data.GetZMin()) < coord_tol); - assert(std::abs(eqdsk_data.GetZ(NH - 1) - eqdsk_data.GetZMax()) < - coord_tol); + // Create test psi_grid and I_psi data + const int test_npsi = 10; + std::vector test_psi_grid(test_npsi); + std::vector test_I_psi(test_npsi); + for (int i = 0; i < test_npsi; ++i) { + test_psi_grid[i] = static_cast(i * 0.5); + test_I_psi[i] = static_cast(1.0 + i * 0.3); + } - std::cout << "\nTesting SplineFunctionSpace evaluation...\n"; + 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) @@ -98,16 +195,14 @@ int main(int argc, char** argv) auto layout = spline_space.GetLayout(); const int layout_size = layout->OwnedSize(); const int psizr_size = static_cast(eqdsk_data.PSIZR.extent(0)); - std::cout << "Layout OwnedSize: " << layout_size << "\n"; - std::cout << "PSIZR size: " << psizr_size << " (NW=" << NW << " * NH=" << NH - << ")\n"; - assert(layout_size == psizr_size); + REQUIRE(layout_size == psizr_size); psi_field.GetData().SetDOFHolderData( pcms::Rank1View( 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 = @@ -157,24 +252,9 @@ int main(int argc, char** argv) auto results_host = Kokkos::create_mirror_view_and_copy( pcms::HostMemorySpace(), eval_results_1d); - std::cout << "Spline evaluation results:\n"; + // Verify that all results are finite for (int i = 0; i < num_eval_points; ++i) { - std::cout << " Point " << i << ": (R=" << eval_coords[2 * i] - << ", Z=" << eval_coords[2 * i + 1] - << ") -> Psi = " << results_host(i) << "\n"; + REQUIRE(std::isfinite(results_host(i))); } - - for (int i = 0; i < num_eval_points; ++i) { - assert(std::isfinite(results_host(i))); - } - - std::cout << "SplineFunctionSpace evaluation test passed!\n"; - - std::cout << "\nAll tests passed!\n"; - return 0; - - } catch (const std::exception& e) { - std::cerr << "Error: " << e.what() << "\n"; - return 1; } } From bc8603b9379422ea69713ea4896944c4a9546f8c Mon Sep 17 00:00:00 2001 From: Sichao25 Date: Sat, 16 May 2026 00:48:15 -0400 Subject: [PATCH 3/4] debug clang-tidy issues --- src/pcms/utility/eqdsk.cpp | 11 +++++++---- test/test_eqdsk.cpp | 4 ++-- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/src/pcms/utility/eqdsk.cpp b/src/pcms/utility/eqdsk.cpp index 2ff787d6..0737f3c4 100644 --- a/src/pcms/utility/eqdsk.cpp +++ b/src/pcms/utility/eqdsk.cpp @@ -101,7 +101,8 @@ void read_geqdsk(const std::string& filename, EQDSKData& eqdsk) eqdsk.grid.bot_left[1] = static_cast(z_min); // Allocate and copy PSIZR - eqdsk.PSIZR = Kokkos::View("PSIZR", mw * mh); + eqdsk.PSIZR = Kokkos::View( + "PSIZR", static_cast(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) { @@ -200,7 +201,8 @@ void read_eqd(const std::string& filename, EQDSKData& eqdsk) eqdsk.grid.bot_left[1] = static_cast(z_min); // Allocate and copy PSIZR - eqdsk.PSIZR = Kokkos::View("PSIZR", mw * mh); + eqdsk.PSIZR = Kokkos::View( + "PSIZR", static_cast(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) { @@ -264,7 +266,7 @@ EQDSKData::EQDSKData(const Uniform2DGrid& grid_in, const int mh = grid.divisions[1]; // Validate input dimensions - if (psirz.size() != static_cast(mw * mh)) { + if (psirz.size() != static_cast(mw) * mh) { throw std::runtime_error("psirz size mismatch: expected " + std::to_string(mw * mh) + " but got " + std::to_string(psirz.size())); @@ -279,7 +281,8 @@ EQDSKData::EQDSKData(const Uniform2DGrid& grid_in, const int npsi = static_cast(psi_grid_vals.size()); // Allocate and copy PSIZR - PSIZR = Kokkos::View("PSIZR", mw * mh); + PSIZR = Kokkos::View("PSIZR", + static_cast(mw) * mh); auto psizr_host = Kokkos::create_mirror_view(PSIZR); for (int i = 0; i < mw * mh; ++i) { psizr_host(i) = psirz[i]; diff --git a/test/test_eqdsk.cpp b/test/test_eqdsk.cpp index 0c9c3874..1e3d164a 100644 --- a/test/test_eqdsk.cpp +++ b/test/test_eqdsk.cpp @@ -156,7 +156,7 @@ TEST_CASE("EQDSKData with SplineFunctionSpace") const int NH = 8; // Create test PSIZR data with a simple quadratic pattern - std::vector test_psirz(NW * NH); + std::vector test_psirz(static_cast(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); @@ -221,7 +221,7 @@ TEST_CASE("EQDSKData with SplineFunctionSpace") auto eval_coords_host = Kokkos::View( "eval_coords_host", num_eval_points, 2); - for (int i = 0; i < num_eval_points; ++i) { + 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 } From bb5e0cde9a31e760d7a761ee3b96826d042f3795 Mon Sep 17 00:00:00 2001 From: Sichao25 Date: Mon, 18 May 2026 12:37:55 -0400 Subject: [PATCH 4/4] remove debug cmake comment --- test/CMakeLists.txt | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 74e072d8..1c603c45 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -34,7 +34,6 @@ if(PCMS_ENABLE_OMEGA_H) ${notRendezvous} ${cyclone1p}) endif() - add_exe(test_ohClassPtn_appRibPtn) set(cyclone2p ${PCMS_TEST_DATA_DIR}/cyclone/23elements/2p.osh/) if(HOST_NPROC GREATER_EQUAL 4) @@ -448,18 +447,6 @@ if(Catch2_FOUND) --data_root ${PCMS_TEST_DATA_DIR}/ltx-interpolation-case/ ) - # EQDSK test with Catch2 - #if(PCMS_ENABLE_OMEGA_H) - #add_executable(test_eqdsk test_eqdsk.cpp) - #target_link_libraries(test_eqdsk PUBLIC - # Catch2::Catch2WithMain - # pcms::utility - # pcms::core - # ) - # target_include_directories(test_eqdsk PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) - # Catch_discover_tests(test_eqdsk) - #endif() - include(Catch) Catch_discover_tests(unit_tests) else()