[cpp] Cleanup sim2d

tobanteAudio · Sep 28, 2024 · 976f5eb · 976f5eb
1 parent a6ff6a7
commit 976f5eb
Show file tree

Hide file tree

Showing 14 changed files with 93 additions and 88 deletions.
diff --git a/models/Modes2D/Modes2D.py b/models/Modes2D/Modes2D.py
@@ -19,13 +19,13 @@ def model(*, Lx=None, Ly=None, Nx=None, Ny=None, dx=None, X=None, Y=None, in_mas
     room=(3.65, 6),
     Tc=20,
     rh=50,
-    fmax=10000.0,
+    fmax=1000.0,
     ppw=10.5,
     duration=8.0,
     refl_coeff=0.99,
     model_factory=model,
     apply_loss=True,
-    diff=True,
+    diff=False,
     image=True,
     verbose=True,
 )
diff --git a/models/ProStudio/ProStudio.py b/models/ProStudio/ProStudio.py
@@ -38,10 +38,10 @@ class ProStudio(Setup3D):
     rh = 50
     fcc = False
     ppw = 7.7
-    fmax = 2000.0
+    fmax = 800.0
     save_folder = '../../sim_data/ProStudio/cpu'
     save_folder_gpu = '../../sim_data/ProStudio/gpu'
-    draw_vox = True
+    draw_vox = False
     draw_backend = 'polyscope'
     compress = 0
     rot_az_el = [0, 0]

diff --git a/models/ProStudio/model.json b/models/ProStudio/model.json
diff --git a/run_2d.sh b/run_2d.sh
@@ -7,7 +7,7 @@
 set -e
 
 build_dir=build
-build_dir=cmake-build-sycl
+# build_dir=cmake-build-sycl
 
 root_dir="$(cd "$(dirname "$0")" && pwd)"
 python_dir="$root_dir/src/python"
@@ -28,7 +28,11 @@ 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 -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 -s "$sim_dir" -e cpu
 # pffdtd sim2d run --sim_dir "$sim_dir" --video
 
 # Post-process

diff --git a/run_3d.sh b/run_3d.sh
@@ -7,17 +7,17 @@ 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="ProStudio"
+sim_name="Modes"
 sim_setup="${sim_name}.py"
 sim_dir="$root_dir/sim_data/$sim_name/gpu"
 
 model_dir="$root_dir/models/$sim_name"
 materials_dir="$root_dir/materials"
 
 fmin=20
-fmax=2000
+fmax=800
 smoothing=0
 
 # Delete old sim
@@ -28,7 +28,7 @@ cd "$model_dir"
 pffdtd sim3d setup "$sim_setup"
 
 # Run sim
-$pffdtd_engine sim3d -p "64" -s "$sim_dir"
+$pffdtd_engine sim3d -e cpu -p "32" -s "$sim_dir"
 # pffdtd sim3d engine --sim_dir="$sim_dir" --plot --draw_backend="mayavi" --json_model="${model_dir}/model.json"
 
 # Post-process

diff --git a/src/cpp/pffdtd/engine_cpu_2d.cpp b/src/cpp/pffdtd/engine_cpu_2d.cpp
@@ -3,6 +3,7 @@
 
 #include "engine_cpu_2d.hpp"
 
+#include "pffdtd/double.hpp"
 #include "pffdtd/progress.hpp"
 #include "pffdtd/time.hpp"
 
@@ -16,30 +17,24 @@
 
 namespace pffdtd {
 
-namespace {
-[[nodiscard]] constexpr auto
-to_ixy(std::integral auto x, std::integral auto y, std::integral auto /*Nx*/, std::integral auto Ny) -> std::integral
-    auto {
-  return x * Ny + y;
-}
-} // namespace
-
 auto EngineCPU2D::operator()(Simulation2D const& sim) const -> stdex::mdarray<double, stdex::dextents<size_t, 2>> {
+  // using Real = double;
+  using Real = Double<double>;
+
+  summary(sim);
 
   auto const Nx         = sim.Nx;
   auto const Ny         = sim.Ny;
+  auto const Npts       = sim.Nx * sim.Ny;
   auto const Nt         = sim.Nt;
   auto const Nb         = static_cast<int64_t>(sim.adj_bn.size());
-  auto const inx        = sim.inx;
-  auto const iny        = sim.iny;
+  auto const Ns         = static_cast<int64_t>(sim.in_ixy.size());
   auto const Nr         = static_cast<int64_t>(sim.out_ixy.size());
   auto const lossFactor = sim.loss_factor;
 
-  summary(sim);
-
-  auto u0_buf  = stdex::mdarray<double, stdex::dextents<size_t, 2>>(Nx, Ny);
-  auto u1_buf  = stdex::mdarray<double, stdex::dextents<size_t, 2>>(Nx, Ny);
-  auto u2_buf  = stdex::mdarray<double, stdex::dextents<size_t, 2>>(Nx, Ny);
+  auto u0_buf  = stdex::mdarray<Real, stdex::dextents<size_t, 1>>(Npts);
+  auto u1_buf  = stdex::mdarray<Real, stdex::dextents<size_t, 1>>(Npts);
+  auto u2_buf  = stdex::mdarray<Real, stdex::dextents<size_t, 1>>(Npts);
   auto out_buf = stdex::mdarray<double, stdex::dextents<size_t, 2>>(Nr, Nt);
 
   auto u0  = u0_buf.to_mdspan();
@@ -61,15 +56,15 @@ auto EngineCPU2D::operator()(Simulation2D const& sim) const -> stdex::mdarray<do
 #pragma omp parallel for
       for (int64_t x = 1; x < Nx - 1; ++x) {
         for (int64_t y = 1; y < Ny - 1; ++y) {
-          auto const idx    = to_ixy(x, y, 0, Ny);
-          auto const left   = u1.data_handle()[idx - 1];
-          auto const right  = u1.data_handle()[idx + 1];
-          auto const bottom = u1.data_handle()[idx - Ny];
-          auto const top    = u1.data_handle()[idx + Ny];
-          auto const last   = u2.data_handle()[idx];
-          auto const tmp    = 0.5 * (left + right + bottom + top) - last;
-
-          u0.data_handle()[idx] = sim.in_mask[idx] * tmp;
+          auto const idx    = x * Ny + y;
+          auto const left   = u1[idx - 1];
+          auto const right  = u1[idx + 1];
+          auto const bottom = u1[idx - Ny];
+          auto const top    = u1[idx + Ny];
+          auto const last   = u2[idx];
+          auto const tmp    = static_cast<Real>(0.5) * (left + right + bottom + top) - last;
+
+          u0[idx] = sim.in_mask[idx] * tmp;
         }
       }
     });
@@ -79,61 +74,63 @@ auto EngineCPU2D::operator()(Simulation2D const& sim) const -> stdex::mdarray<do
 #pragma omp parallel for
       for (int64_t i = 0; i < Nb; ++i) {
         auto const ib = sim.bn_ixy[i];
-        auto const K  = static_cast<double>(sim.adj_bn[i]);
+        auto const K  = static_cast<Real>(sim.adj_bn[i]);
 
-        auto const last1 = u1.data_handle()[ib];
-        auto const last2 = u2.data_handle()[ib];
+        auto const last1 = u1[ib];
+        auto const last2 = u2[ib];
 
-        auto const left      = u1.data_handle()[ib - 1];
-        auto const right     = u1.data_handle()[ib + 1];
-        auto const bottom    = u1.data_handle()[ib - Ny];
-        auto const top       = u1.data_handle()[ib + Ny];
+        auto const left      = u1[ib - 1];
+        auto const right     = u1[ib + 1];
+        auto const bottom    = u1[ib - Ny];
+        auto const top       = u1[ib + Ny];
         auto const neighbors = left + right + top + bottom;
 
-        u0.data_handle()[ib] = (2 - 0.5 * K) * last1 + 0.5 * neighbors - last2;
+        u0[ib] = (Real(2) - Real(0.5) * K) * last1 + Real(0.5) * neighbors - last2;
       }
 
 // Boundary Loss
 #pragma omp parallel for
       for (int64_t i = 0; i < Nb; ++i) {
         auto const ib = sim.bn_ixy[i];
         auto const K  = sim.adj_bn[i];
-        auto const K4 = static_cast<double>(4 - K);
+        auto const K4 = static_cast<Real>(4 - K);
         auto const lf = lossFactor;
 
-        auto const current = u0.data_handle()[ib];
-        auto const prev    = u2.data_handle()[ib];
+        auto const current = u0[ib];
+        auto const prev    = u2[ib];
 
-        u0.data_handle()[ib] = (current + lf * K4 * prev) / (1 + lf * K4);
+        u0[ib] = (current + lf * K4 * prev) / (Real(1) + lf * K4);
       }
     });
 
-    // Copy Input
-    u0(inx, iny) += sim.src_sig[n];
+    // Add sources
+    for (int64_t s = 0; s < Ns; ++s) {
+      u0[sim.in_ixy[s]] += sim.in_sigs[n];
+    }
 
-    // Copy Output
-    for (int64_t i = 0; i < Nr; ++i) {
-      auto r_ixy = sim.out_ixy[i];
-      out(i, n)  = u0.data_handle()[r_ixy];
+    // Read outputs
+    for (int64_t r = 0; r < Nr; ++r) {
+      out(r, n) = static_cast<double>(u0[sim.out_ixy[r]]);
     }
 
+    // Rotate buffers
     auto tmp = u2;
     u2       = u1;
     u1       = u0;
     u0       = tmp;
 
-    auto const now = getTime();
-
+    auto const now           = getTime();
     auto const elapsed       = now - start;
     auto const elapsedSample = now - sampleStart;
+
     elapsedAir += elapsedAirSample;
     elapsedBoundary += elapsedBoundarySample;
 
     print(ProgressReport{
         .n                     = n,
         .Nt                    = Nt,
         .Npts                  = Nx * Ny,
-        .Nb                    = static_cast<int64_t>(Nb),
+        .Nb                    = Nb,
         .elapsed               = elapsed,
         .elapsedSample         = elapsedSample,
         .elapsedAir            = elapsedAir,

diff --git a/src/cpp/pffdtd/engine_metal.metal b/src/cpp/pffdtd/engine_metal.metal
@@ -78,12 +78,12 @@ namespace sim2d {
 
 [[kernel]] void addSource(
     device float* u0 [[buffer(0)]],
-    device float const* src_sig [[buffer(1)]],
+    device float const* in_sigs [[buffer(1)]],
     constant Constants2D<float>& constants [[buffer(2)]],
     constant int64_t& timestep [[buffer(3)]],
     uint id [[thread_position_in_grid]]
 ) {
-  u0[constants.in_ixy] += src_sig[timestep];
+  u0[constants.in_ixy] += in_sigs[timestep];
 }
 
 [[kernel]] void readOutput(

diff --git a/src/cpp/pffdtd/engine_metal_2d.mm b/src/cpp/pffdtd/engine_metal_2d.mm
@@ -32,17 +32,15 @@ auto run(Simulation2D const& sim) {
     auto const Npts        = Nx * Ny;
     auto const Nt          = sim.Nt;
     auto const Nb          = static_cast<int64_t>(sim.adj_bn.size());
-    auto const inx         = sim.inx;
-    auto const iny         = sim.iny;
     auto const Nr          = static_cast<int64_t>(sim.out_ixy.size());
     auto const loss_factor = sim.loss_factor;
-    auto const src_sig_f32 = toFloat(sim.src_sig);
+    auto const in_sigs_f32 = toFloat(sim.in_sigs);
 
     auto const c = Constants2D<float>{
         .Nx         = Nx,
         .Ny         = Ny,
         .Nt         = Nt,
-        .in_ixy     = inx * Ny + iny,
+        .in_ixy     = sim.in_ixy.at(0),
         .lossFactor = static_cast<float>(loss_factor),
     };
 
@@ -67,7 +65,7 @@ auto run(Simulation2D const& sim) {
     id<MTLBuffer> bn_ixy    = makeBuffer(device, sim.bn_ixy, MTLResourceStorageModeShared);
     id<MTLBuffer> adj_bn    = makeBuffer(device, sim.adj_bn, MTLResourceStorageModeShared);
     id<MTLBuffer> out_ixy   = makeBuffer(device, sim.out_ixy, MTLResourceStorageModeShared);
-    id<MTLBuffer> src_sig   = makeBuffer(device, src_sig_f32, MTLResourceStorageModeShared);
+    id<MTLBuffer> in_sigs   = makeBuffer(device, in_sigs_f32, MTLResourceStorageModeShared);
     id<MTLBuffer> constants = [device newBufferWithBytes:&c
                                                   length:sizeof(Constants2D<float>)
                                                  options:MTLResourceStorageModeShared];
@@ -132,7 +130,7 @@ auto run(Simulation2D const& sim) {
       id<MTLComputeCommandEncoder> addSource = [commandBuffer computeCommandEncoder];
       [addSource setComputePipelineState:addSourceKernel];
       [addSource setBuffer:u0 offset:0 atIndex:0];
-      [addSource setBuffer:src_sig offset:0 atIndex:1];
+      [addSource setBuffer:in_sigs offset:0 atIndex:1];
       [addSource setBuffer:constants offset:0 atIndex:2];
       [addSource setBuffer:timestep offset:0 atIndex:3];
       [addSource dispatchThreads:MTLSizeMake(1, 1, 1) threadsPerThreadgroup:MTLSizeMake(1, 1, 1)];

diff --git a/src/cpp/pffdtd/engine_sycl_2d.cpp b/src/cpp/pffdtd/engine_sycl_2d.cpp
@@ -106,7 +106,7 @@ auto EngineSYCL2D::operator()(Simulation2D const& sim) const -> stdex::mdarray<d
   auto bn_ixy_buf  = sycl::buffer<int64_t, 1>{sim.bn_ixy};
   auto adj_bn_buf  = sycl::buffer<int64_t, 1>{sim.adj_bn};
   auto out_ixy_buf = sycl::buffer<int64_t, 1>{sim.out_ixy};
-  auto src_sig_buf = sycl::buffer<double, 1>{sim.src_sig};
+  auto src_sig_buf = sycl::buffer<double, 1>{sim.in_sigs};
 
   auto elapsedAir      = std::chrono::nanoseconds{0};
   auto elapsedBoundary = std::chrono::nanoseconds{0};
@@ -151,9 +151,9 @@ auto EngineSYCL2D::operator()(Simulation2D const& sim) const -> stdex::mdarray<d
 
     queue.submit([&](sycl::handler& cgh) {
       auto u0      = sycl::accessor{u0_buf, cgh, sycl::read_write};
-      auto src_sig = sycl::accessor{src_sig_buf, cgh, sycl::read_only};
+      auto in_sigs = sycl::accessor{src_sig_buf, cgh, sycl::read_only};
 
-      cgh.parallel_for<struct CopyInput>(1, [=](sycl::id<1>) { u0[inx][iny] += src_sig[n]; });
+      cgh.parallel_for<struct CopyInput>(1, [=](sycl::id<1>) { u0[inx][iny] += in_sigs[n]; });
     });
 
     queue.submit([&](sycl::handler& cgh) {

diff --git a/src/cpp/pffdtd/simulation_2d.cpp b/src/cpp/pffdtd/simulation_2d.cpp
@@ -12,8 +12,13 @@ namespace pffdtd {
 auto loadSimulation2D(std::filesystem::path const& dir) -> Simulation2D {
   auto sim = HDF5Reader{dir / "sim.h5"};
 
-  auto const Nx = sim.read<int64_t>("Nx");
-  auto const Ny = sim.read<int64_t>("Ny");
+  auto const Nx  = sim.read<int64_t>("Nx");
+  auto const Ny  = sim.read<int64_t>("Ny");
+  auto const inx = sim.read<int64_t>("inx");
+  auto const iny = sim.read<int64_t>("iny");
+
+  auto in_ixy = std::vector<int64_t>{};
+  in_ixy.push_back(inx * Ny + iny);
 
   return Simulation2D{
       .dir = dir,
@@ -27,10 +32,11 @@ auto loadSimulation2D(std::filesystem::path const& dir) -> Simulation2D {
       .bn_ixy      = sim.read<std::vector<int64_t>>("bn_ixy"),
       .loss_factor = sim.read<double>("loss_factor"),
 
-      .inx     = sim.read<int64_t>("inx"),
-      .iny     = sim.read<int64_t>("iny"),
-      .src_sig = sim.read<std::vector<double>>("src_sig"),
+      .inx = sim.read<int64_t>("inx"),
+      .iny = sim.read<int64_t>("iny"),
 
+      .in_sigs = sim.read<std::vector<double>>("in_sigs"),
+      .in_ixy  = std::move(in_ixy),
       .out_ixy = sim.read<std::vector<int64_t>>("out_ixy"),
   };
 }
@@ -40,13 +46,12 @@ auto summary(Simulation2D const& sim) -> void {
   fmt::println("Nx: {}", sim.Nx);
   fmt::println("Ny: {}", sim.Ny);
   fmt::println("N: {}", sim.Nx * sim.Ny);
-  fmt::println("inx: {}", sim.inx);
-  fmt::println("iny: {}", sim.iny);
   fmt::println("in_mask: {}", sim.in_mask.size());
   fmt::println("bn_ixy: {}", sim.bn_ixy.size());
   fmt::println("adj_bn: {}", sim.adj_bn.size());
+  fmt::println("in_sigs: {}", sim.in_sigs.size());
+  fmt::println("in_ixy: {}", sim.in_ixy.size());
   fmt::println("out_ixy: {}", sim.out_ixy.size());
-  fmt::println("src_sig: {}", sim.src_sig.size());
   fmt::println("loss_factor: {}", sim.loss_factor);
 }
 

diff --git a/src/cpp/pffdtd/simulation_2d.hpp b/src/cpp/pffdtd/simulation_2d.hpp
@@ -21,10 +21,11 @@ struct Simulation2D {
   std::vector<int64_t> bn_ixy;  // Boundary nodes
   double loss_factor;           // Boundary loss
 
-  int64_t inx;                 // Source position x
-  int64_t iny;                 // Source position y
-  std::vector<double> src_sig; // Source signal
+  int64_t inx; // Source position x
+  int64_t iny; // Source position y
 
+  std::vector<double> in_sigs;  // Source signal
+  std::vector<int64_t> in_ixy;  // Source nodes
   std::vector<int64_t> out_ixy; // Receiver nodes
 };