[cpp] Make precision a runtime parameter for sim3d

run_2d.sh

@@ -28,7 +28,7 @@ cd "$model_dir"
 python "$sim_setup"
 
 # Run sim
-DPCPP_CPU_PLACES=cores DPCPP_CPU_CU_AFFINITY=spread DPCPP_CPU_NUM_CUS=$jobs OMP_NUM_THREADS=$jobs "$engine_exe" sim2d -p 32 -s "$sim_dir" -e sycl
+DPCPP_CPU_PLACES=cores DPCPP_CPU_CU_AFFINITY=spread DPCPP_CPU_NUM_CUS=$jobs OMP_NUM_THREADS=$jobs "$engine_exe" sim2d -p 64 -s "$sim_dir" -e sycl
 # pffdtd sim2d run --sim_dir "$sim_dir" --video
 
 # Post-process

run_3d.sh

@@ -7,7 +7,7 @@ set -e
 
 root_dir="$(cd "$(dirname "$0")" && pwd)"
 pffdtd_engine="$root_dir/build/src/cpp/pffdtd-engine"
-# pffdtd_engine="$root_dir/cmake-build-cuda/src/cpp/pffdtd-engine"
+pffdtd_engine="$root_dir/cmake-build-cuda/src/cpp/pffdtd-engine"
 
 sim_name="Modes"
 sim_setup="${sim_name}.py"
@@ -28,7 +28,7 @@ cd "$model_dir"
 pffdtd sim3d setup "$sim_setup"
 
 # Run sim
-$pffdtd_engine sim3d -e cpu -p "32" -s "$sim_dir"
+$pffdtd_engine sim3d -e cuda -p "64" -s "$sim_dir"
 # pffdtd sim3d engine --sim_dir="$sim_dir" --plot --draw_backend="mayavi" --json_model="${model_dir}/model.json"
 
 # Post-process

src/cpp/main.cpp

@@ -69,18 +69,15 @@ namespace {
   return engines;
 }
 
-template<typename Real>
 [[nodiscard]] auto getEngines3D() {
   using namespace pffdtd;
-  auto engines   = std::map<std::string, std::function<void(Simulation3D<Real> const&)>>{};
+  auto engines   = std::map<std::string, std::function<void(Simulation3D const&)>>{};
   engines["cpu"] = EngineCPU3D{};
 #if defined(PFFDTD_HAS_CUDA)
   engines["cuda"] = EngineCUDA3D{};
 #endif
 #if defined(PFFDTD_HAS_METAL)
-  if constexpr (std::same_as<Real, float>) {
-    engines["metal"] = EngineMETAL3D{};
-  }
+  engines["metal"] = EngineMETAL3D{};
 #endif
 #if defined(PFFDTD_HAS_SYCL)
   engines["sycl"] = EngineSYCL3D{};
@@ -111,18 +108,17 @@ struct Arguments {
   Sim3D sim3d;
 };
 
-template<typename Real>
 auto run3D(Arguments::Sim3D const& args) {
   using namespace pffdtd;
   fmt::println("Running: {} on {} with precision {}", args.simDir, args.engine, toString(args.precision));
 
-  auto const engines = getEngines3D<Real>();
+  auto const engines = getEngines3D();
   auto const& engine = engines.at(args.engine);
 
   auto const simDir = std::filesystem::path{args.simDir};
   auto const start  = getTime();
 
-  auto sim = loadSimulation3D<Real>(simDir);
+  auto sim = loadSimulation3D(simDir, args.precision);
   scaleInput(sim);
   engine(sim);
   rescaleOutput(sim);
@@ -177,13 +173,7 @@ auto main(int argc, char** argv) -> int {
   }
 
   if (*sim3d) {
-    if (args.sim3d.precision == pffdtd::Precision::Float) {
-      run3D<float>(args.sim3d);
-    } else if (args.sim3d.precision == pffdtd::Precision::Double) {
-      run3D<double>(args.sim3d);
-    } else {
-      pffdtd::raisef<std::invalid_argument>("invalid precision '{}'", toString(args.sim3d.precision));
-    }
+    run3D(args.sim3d);
   }
 
   if (*test) {

src/cpp/pffdtd/engine_cpu_3d.cpp

@@ -4,6 +4,8 @@
 
 #include "engine_cpu_3d.hpp"
 
+#include "pffdtd/double.hpp"
+#include "pffdtd/exception.hpp"
 #include "pffdtd/progress.hpp"
 #include "pffdtd/time.hpp"
 #include "pffdtd/utility.hpp"
@@ -79,7 +81,11 @@ auto process_bnl_fd(
 }
 
 template<typename Real>
-auto run(Simulation3D<Real> const& sim) -> void {
+auto run(Simulation3D const& sim) -> void {
+  auto const ssaf_bnl_real  = convertTo<Real>(sim.ssaf_bnl);
+  auto const mat_beta_real  = convertTo<Real>(sim.mat_beta);
+  auto const mat_quads_real = convertTo<Real>(sim.mat_quads);
+
   // keep local ints, scalars
   int64_t const Ns   = sim.Ns;
   int64_t const Nr   = sim.Nr;
@@ -105,9 +111,9 @@ auto run(Simulation3D<Real> const& sim) -> void {
   int8_t const* mat_bnl          = sim.mat_bnl.data();
   int8_t const* Q_bna            = sim.Q_bna.data();
   double const* in_sigs          = sim.in_sigs.data();
-  Real const* ssaf_bnl           = sim.ssaf_bnl.data();
-  Real const* mat_beta           = sim.mat_beta.data();
-  MatQuad<Real> const* mat_quads = sim.mat_quads.data();
+  Real const* ssaf_bnl           = ssaf_bnl_real.data();
+  Real const* mat_beta           = mat_beta_real.data();
+  MatQuad<Real> const* mat_quads = mat_quads_real.data();
   double* u_out                  = sim.u_out.get();
 
   // allocate memory
@@ -130,11 +136,11 @@ auto run(Simulation3D<Real> const& sim) -> void {
   auto* gh1  = gh1_buf.data();
 
   // sim coefficients
-  auto const lo2 = sim.lo2;
-  auto const sl2 = sim.sl2;
-  auto const l   = sim.l;
-  auto const a1  = sim.a1;
-  auto const a2  = sim.a2;
+  auto const lo2 = static_cast<Real>(sim.lo2);
+  auto const sl2 = static_cast<Real>(sim.sl2);
+  auto const l   = static_cast<Real>(sim.l);
+  auto const a1  = static_cast<Real>(sim.a1);
+  auto const a2  = static_cast<Real>(sim.a2);
 
   // can control outside with OMP_NUM_THREADS env variable
   int const numWorkers = omp_get_max_threads();
@@ -373,8 +379,14 @@ auto run(Simulation3D<Real> const& sim) -> void {
 
 } // namespace
 
-auto EngineCPU3D::operator()(Simulation3D<float> const& sim) const -> void { run(sim); }
-
-auto EngineCPU3D::operator()(Simulation3D<double> const& sim) const -> void { run(sim); }
+auto EngineCPU3D::operator()(Simulation3D const& sim) const -> void {
+  switch (sim.precision) {
+    case Precision::Float: return run<float>(sim);
+    case Precision::Double: return run<double>(sim);
+    case Precision::DoubleFloat: return run<Double<float>>(sim);
+    case Precision::DoubleDouble: return run<Double<double>>(sim);
+    default: raisef<std::invalid_argument>("invalid precision {}", static_cast<int>(sim.precision));
+  }
+}
 
 } // namespace pffdtd

src/cpp/pffdtd/engine_cpu_3d.hpp

@@ -4,13 +4,13 @@
 
 #pragma once
 
+#include "pffdtd/precision.hpp"
 #include "pffdtd/simulation_3d.hpp"
 
 namespace pffdtd {
 
 struct EngineCPU3D {
-  auto operator()(Simulation3D<float> const& sim) const -> void;
-  auto operator()(Simulation3D<double> const& sim) const -> void;
+  auto operator()(Simulation3D const& sim) const -> void;
 };
 
 } // namespace pffdtd

src/cpp/pffdtd/engine_cuda_3d.cu

@@ -543,8 +543,7 @@ auto print_gpu_details(int i) -> uint64_t {
 }
 
 // input indices need to be sorted for multi-device allocation
-template<typename Real>
-void checkSorted(Simulation3D<Real> const& sim) {
+void checkSorted(Simulation3D const& sim) {
   int64_t const* bn_ixyz  = sim.bn_ixyz.data();
   int64_t const* bnl_ixyz = sim.bnl_ixyz.data();
   int64_t const* bna_ixyz = sim.bna_ixyz.data();
@@ -574,7 +573,7 @@ void checkSorted(Simulation3D<Real> const& sim) {
 
 // counts for splitting data across GPUs
 template<typename Real>
-void splitData(Simulation3D<Real> const& sim, std::span<HostData<Real>> ghds) {
+void splitData(Simulation3D const& sim, std::span<HostData<Real>> ghds) {
   auto const Nx    = sim.Nx;
   auto const Ny    = sim.Ny;
   auto const Nz    = sim.Nz;
@@ -724,7 +723,7 @@ void splitData(Simulation3D<Real> const& sim, std::span<HostData<Real>> ghds) {
 }
 
 template<typename Real>
-auto run(Simulation3D<Real> const& sim) -> void { // NOLINT(readability-function-cognitive-complexity)
+auto run(Simulation3D const& sim) -> void { // NOLINT(readability-function-cognitive-complexity)
   // if you want to test synchronous, env variable for that
   auto const* s = std::getenv("CUDA_LAUNCH_BLOCKING"); // NOLINT(concurrency-mt-unsafe)
   if (s != nullptr) {
@@ -745,6 +744,10 @@ auto run(Simulation3D<Real> const& sim) -> void { // NOLINT(readability-function
   auto ghds = std::vector<HostData<Real>>(static_cast<size_t>(ngpus));
   auto gds  = std::vector<DeviceData<Real>>(static_cast<size_t>(ngpus));
 
+  auto const ssaf_bnl_real  = convertTo<Real>(sim.ssaf_bnl);
+  auto const mat_beta_real  = convertTo<Real>(sim.mat_beta);
+  auto const mat_quads_real = convertTo<Real>(sim.mat_quads);
+
   if (ngpus > 1) {
     checkSorted(sim); // needs to be sorted for multi-GPU
   }
@@ -756,11 +759,11 @@ auto run(Simulation3D<Real> const& sim) -> void { // NOLINT(readability-function
     gds[gid].totalmembytes = print_gpu_details(gid);
   }
 
-  Real lo2 = sim.lo2;
-  Real a1  = sim.a1;
-  Real a2  = sim.a2;
-  Real l   = sim.l;
-  Real sl2 = sim.sl2;
+  auto lo2 = static_cast<Real>(sim.lo2);
+  auto a1  = static_cast<Real>(sim.a1);
+  auto a2  = static_cast<Real>(sim.a2);
+  auto l   = static_cast<Real>(sim.l);
+  auto sl2 = static_cast<Real>(sim.sl2);
 
   // timing stuff
   auto elapsed               = std::chrono::nanoseconds{0};
@@ -829,7 +832,7 @@ auto run(Simulation3D<Real> const& sim) -> void { // NOLINT(readability-function
 
     // aliased pointers (to memory already allocated)
     host.in_sigs   = sim.in_sigs.data() + Ns_read * sim.Nt;
-    host.ssaf_bnl  = sim.ssaf_bnl.data() + Nbl_read;
+    host.ssaf_bnl  = ssaf_bnl_real.data() + Nbl_read;
     host.adj_bn    = sim.adj_bn.data() + Nb_read;
     host.mat_bnl   = sim.mat_bnl.data() + Nbl_read;
     host.K_bn      = sim.K_bn.data() + Nb_read;
@@ -942,12 +945,12 @@ auto run(Simulation3D<Real> const& sim) -> void { // NOLINT(readability-function
     gpuErrchk(cudaMemcpy(gpu.mat_bnl, host.mat_bnl, (size_t)host.Nbl * sizeof(int8_t), cudaMemcpyHostToDevice));
 
     gpuErrchk(cudaMalloc(&(gpu.mat_beta), (size_t)sim.Nm * sizeof(Real)));
-    gpuErrchk(cudaMemcpy(gpu.mat_beta, sim.mat_beta.data(), (size_t)sim.Nm * sizeof(Real), cudaMemcpyHostToDevice));
+    gpuErrchk(cudaMemcpy(gpu.mat_beta, mat_beta_real.data(), (size_t)sim.Nm * sizeof(Real), cudaMemcpyHostToDevice));
 
     gpuErrchk(cudaMalloc(&(gpu.mat_quads), (size_t)sim.Nm * MMb * sizeof(MatQuad<Real>)));
     gpuErrchk(cudaMemcpy(
         gpu.mat_quads,
-        sim.mat_quads.data(),
+        mat_quads_real.data(),
         (size_t)sim.Nm * MMb * sizeof(MatQuad<Real>),
         cudaMemcpyHostToDevice
     ));
@@ -1407,8 +1410,12 @@ auto run(Simulation3D<Real> const& sim) -> void { // NOLINT(readability-function
 
 } // namespace
 
-auto EngineCUDA3D::operator()(Simulation3D<float> const& sim) const -> void { run(sim); }
-
-auto EngineCUDA3D::operator()(Simulation3D<double> const& sim) const -> void { run(sim); }
+auto EngineCUDA3D::operator()(Simulation3D const& sim) const -> void {
+  switch (sim.precision) {
+    case Precision::Float: return run<float>(sim);
+    case Precision::Double: return run<double>(sim);
+    default: throw std::invalid_argument("invalid precision " + std::to_string(static_cast<int>(sim.precision)));
+  }
+}
 
 } // namespace pffdtd

src/cpp/pffdtd/engine_cuda_3d.hpp

@@ -13,8 +13,7 @@
 namespace pffdtd {
 
 struct EngineCUDA3D {
-  auto operator()(Simulation3D<float> const& sim) const -> void;
-  auto operator()(Simulation3D<double> const& sim) const -> void;
+  auto operator()(Simulation3D const& sim) const -> void;
 };
 
 } // namespace pffdtd

src/cpp/pffdtd/engine_metal_2d.mm

@@ -19,7 +19,7 @@
 
 auto toFloat(std::vector<double> const& buf) {
   auto buf32 = std::vector<float>(buf.size());
-  std::ranges::transform(buf32, buf32.begin(), [](auto v) { return static_cast<float>(v); });
+  std::ranges::transform(buf, buf32.begin(), [](auto v) { return static_cast<float>(v); });
   return buf32;
 }
 

src/cpp/pffdtd/engine_metal_3d.hpp

@@ -11,7 +11,7 @@
 namespace pffdtd {
 
 struct EngineMETAL3D {
-  auto operator()(Simulation3D<float> const& sim) const -> void;
+  auto operator()(Simulation3D const& sim) const -> void;
 };
 
 } // namespace pffdtd

src/cpp/pffdtd/engine_metal_3d.mm

@@ -17,7 +17,7 @@
 namespace {
 
 template<typename Real>
-auto run(Simulation3D<Real> const& sim) {
+auto run(Simulation3D const& sim) {
   @autoreleasepool {
 
     // Device
@@ -286,6 +286,9 @@ auto run(Simulation3D<Real> const& sim) {
 
 } // namespace
 
-auto EngineMETAL3D::operator()(Simulation3D<float> const& sim) const -> void { run(sim); }
+auto EngineMETAL3D::operator()(Simulation3D const& sim) const -> void {
+  PFFDTD_ASSERT(precision == Precision::Float);
+  run<float>(sim);
+}
 
 } // namespace pffdtd

src/cpp/pffdtd/engine_sycl_3d.cpp

@@ -3,6 +3,9 @@
 #include "engine_sycl_3d.hpp"
 
 #include "pffdtd/assert.hpp"
+#include "pffdtd/double.hpp"
+#include "pffdtd/exception.hpp"
+#include "pffdtd/precision.hpp"
 #include "pffdtd/progress.hpp"
 #include "pffdtd/sycl.hpp"
 #include "pffdtd/time.hpp"
@@ -39,7 +42,7 @@ template<typename T>
 struct ReadOutput;
 
 template<typename Real>
-auto run(Simulation3D<Real> const& sim) -> void {
+auto run(Simulation3D const& sim) -> void {
   PFFDTD_ASSERT(sim.grid == Grid::CART);
 
   auto queue  = sycl::queue{sycl::property::queue::enable_profiling{}};
@@ -64,6 +67,10 @@ auto run(Simulation3D<Real> const& sim) -> void {
   auto const a1  = static_cast<Real>(sim.a1);
   auto const a2  = static_cast<Real>(sim.a2);
 
+  auto const ssaf_bnl_real  = convertTo<Real>(sim.ssaf_bnl);
+  auto const mat_beta_real  = convertTo<Real>(sim.mat_beta);
+  auto const mat_quads_real = convertTo<Real>(sim.mat_quads);
+
   auto Q_bna_buf     = sycl::buffer{sim.Q_bna};
   auto bn_mask_buf   = sycl::buffer{sim.bn_mask};
   auto adj_bn_buf    = sycl::buffer{sim.adj_bn};
@@ -73,11 +80,11 @@ auto run(Simulation3D<Real> const& sim) -> void {
   auto in_ixyz_buf   = sycl::buffer{sim.in_ixyz};
   auto out_ixyz_buf  = sycl::buffer{sim.out_ixyz};
   auto in_sigs_buf   = sycl::buffer{sim.in_sigs};
-  auto mat_beta_buf  = sycl::buffer{sim.mat_beta};
+  auto mat_beta_buf  = sycl::buffer{mat_beta_real};
   auto mat_bnl_buf   = sycl::buffer{sim.mat_bnl};
-  auto mat_quads_buf = sycl::buffer{sim.mat_quads};
+  auto mat_quads_buf = sycl::buffer{mat_quads_real};
   auto Mb_buf        = sycl::buffer{sim.Mb};
-  auto ssaf_bnl_buf  = sycl::buffer{sim.ssaf_bnl};
+  auto ssaf_bnl_buf  = sycl::buffer{ssaf_bnl_real};
 
   auto u0_buf    = sycl::buffer<Real>(static_cast<size_t>(Npts));
   auto u1_buf    = sycl::buffer<Real>(static_cast<size_t>(Npts));
@@ -343,14 +350,27 @@ auto run(Simulation3D<Real> const& sim) -> void {
   // Copy output to host
   auto host = sycl::host_accessor{u_out_buf, sycl::read_only};
   for (auto i{0UL}; i < static_cast<size_t>(Nr * Nt); ++i) {
-    sim.u_out[i] = host[i];
+    sim.u_out[i] = static_cast<double>(host[i]);
   }
 }
 
 } // namespace
 
-auto EngineSYCL3D::operator()(Simulation3D<float> const& sim) const -> void { run(sim); }
+auto EngineSYCL3D::operator()(Simulation3D const& sim) const -> void {
+  switch (sim.precision) {
+#if defined(__APPLE__) or defined(__clang__)
+    case Precision::Half: return run<_Float16>(sim);
+    case Precision::DoubleHalf: return run<Double<_Float16>>(sim);
+#endif
+
+    case Precision::Float: return run<float>(sim);
+    case Precision::DoubleFloat: return run<Double<float>>(sim);
+
+    case Precision::Double: return run<double>(sim);
+    case Precision::DoubleDouble: return run<Double<double>>(sim);
 
-auto EngineSYCL3D::operator()(Simulation3D<double> const& sim) const -> void { run(sim); }
+    default: raisef<std::invalid_argument>("invalid precision {}", static_cast<int>(sim.precision));
+  }
+}
 
 } // namespace pffdtd

src/cpp/pffdtd/engine_sycl_3d.hpp

@@ -11,8 +11,7 @@
 namespace pffdtd {
 
 struct EngineSYCL3D {
-  auto operator()(Simulation3D<float> const& sim) const -> void;
-  auto operator()(Simulation3D<double> const& sim) const -> void;
+  auto operator()(Simulation3D const& sim) const -> void;
 };
 
 } // namespace pffdtd