Improve 2D simulation

data/models/Diffusor/Diffusor.py

@@ -0,0 +1,129 @@
+import numpy as np
+
+from pffdtd.common.myfuncs import to_ixy, point_on_circle
+from pffdtd.diffusor.design import diffusor_bandwidth
+from pffdtd.diffusor.qrd import quadratic_residue_diffuser
+from pffdtd.diffusor.prd import primitive_root_diffuser
+from pffdtd.sim2d.setup import sim_setup_2d
+
+
+def quantize_point(pos, dx, ret_err=True, scale=False):
+    def quantize(i):
+        low = np.floor(i/dx)
+        high = np.ceil(i/dx)
+        err_low = np.fabs(i-low*dx)
+        err_high = np.fabs(i-high*dx)
+        err = err_low if err_low < err_high else err_high
+        i_q = low if err_low < err_high else high
+        if scale:
+            return i_q*dx, err
+        return i_q, err
+
+    x, y = pos
+    xq, yq = quantize(x), quantize(y)
+    if ret_err:
+        return xq, yq
+    return xq[0], yq[0]
+
+
+def add_diffusor(prime, well_width, max_depth, room, in_mask, X, Y, dx, c, verbose=False):
+    print('--SIM-SETUP: Quantize diffusor')
+    dim = (well_width, max_depth)
+    width = well_width
+    depth = max_depth
+    fmin_t, fmax_t = diffusor_bandwidth(dim[0], dim[1], c=c)
+
+    (width_q, werr_q), (depth_q, derr_q) = quantize_point(dim, dx)
+    width_q *= dx
+    depth_q *= dx
+    print(dim, width_q, depth_q)
+    fmin_q, fmax_q = diffusor_bandwidth(width_q, depth_q, c=c)
+
+    radius = 5
+    total_width = 5
+    pos = (room[0]/2-total_width/2, room[1]/2-radius-depth)
+    pos_q = quantize_point(pos, dx)
+    n = int(total_width/width_q)
+
+    if verbose:
+        print(f"  {pos_q=}")
+        print(f"  {width=:.4f}")
+        print(f"  {depth=:.4f}")
+        print(f"  {fmin_t=:.4f}")
+        print(f"  {fmax_t=:.4f}")
+        print(f"  {width_q=:.4f}")
+        print(f"  {depth_q=:.4f}")
+        print(f"  {fmin_q=:.4f}")
+        print(f"  {fmax_q=:.4f}")
+        print(f"  error_w={werr_q/width*100:.2f}%")
+        print(f"  error_d={derr_q/depth*100:.2f}%")
+
+    print('--SIM-SETUP: Locate diffusor')
+    depths, _ = primitive_root_diffuser(prime, g=None, depth=depth_q)
+    depths = quadratic_residue_diffuser(prime, depth_q)
+    prime = depths.shape[0]
+    for w in range(n):
+        xs = (room[0]/2-total_width/2)+w*width_q
+        xe = xs+width_q
+        ys = room[1]/2-5-depth_q
+        ye = ys+depths[w % prime]+0.05
+        in_mask[(X >= xs) & (Y >= ys) & (X < xe) & (Y < ye)] = False
+
+    return in_mask
+
+
+def make_receiver_arc(count, center, radius, dx, Nx, Ny):
+    x, y = center
+    angles = np.linspace(0.0, 180.0, count, endpoint=True)
+    out_ixy = []
+    out_cart = []
+    for i in range(angles.shape[0]):
+        x, y = point_on_circle(center, radius, np.deg2rad(angles[i]))
+        xq, yq = quantize_point((x, y), dx, ret_err=False)
+        idx = to_ixy(xq, yq, Nx, Ny)
+        out_ixy.append(idx)
+        out_cart.append((xq, yq))
+    return out_ixy, out_cart
+
+
+def diffusor_model_factory(*, Lx=None, Ly=None, Nx=None, Ny=None, dx=None, X=None, Y=None, in_mask=None):
+    in_mask[(X >= 0) & (Y >= 0) & (X < Lx) & (Y < Ly)] = True
+
+    # source input position
+    x_in, y_in = Lx*0.5, Ly*0.5
+    inx = int(np.round(x_in / dx + 0.5) + 1)
+    iny = int(np.round(y_in / dx + 0.5) + 1)
+    assert in_mask[inx, iny]
+
+    # Diffusor
+    print('--SIM-SETUP: Generate diffusor')
+    prime = 17
+    depth = 0.35
+    width = 0.048
+    in_mask = add_diffusor(prime, width, depth, (Lx, Ly),
+                           in_mask, X, Y, dx, 343.0, True)
+
+    # Receiver linear index
+    print('--SIM-SETUP: Generate receivers')
+    arc_radius = 5.0
+    arc_center = (x_in, y_in-arc_radius)
+    out_ixy, _ = make_receiver_arc(180, arc_center, arc_radius, dx, Nx, Ny)
+
+    return in_mask, inx, iny, out_ixy
+
+
+sim_setup_2d(
+    sim_dir='../../sim_data/Diffusor/cpu',
+    room=(30, 30),
+    Tc=20,
+    rh=50,
+    fmax=4000,
+    ppw=10.5,
+    duration=0.050,
+    refl_coeff=0.99,
+    model_factory=diffusor_model_factory,
+    apply_loss=True,
+    diff=True,
+    image=True,
+    verbose=True,
+)

data/models/Modes2D/Modes2D.py

@@ -0,0 +1,28 @@
+from pffdtd.common.myfuncs import to_ixy
+from pffdtd.sim2d.setup import sim_setup_2d
+
+
+def room_modes_model_factory(*, Lx=None, Ly=None, Nx=None, Ny=None, dx=None, X=None, Y=None, in_mask=None):
+    in_mask[(X >= 0) & (Y >= 0) & (X < Lx) & (Y < Ly)] = True
+    inx = 2
+    iny = 2
+    out_ixy = [to_ixy(Nx-4, Ny-4, Nx, Ny)]
+    assert in_mask[inx, iny]
+    return in_mask, inx, iny, out_ixy
+
+
+sim_setup_2d(
+    sim_dir='../../sim_data/Modes2D/cpu',
+    room=(3, 3),
+    Tc=20,
+    rh=50,
+    fmax=1000.0,
+    ppw=10.5,
+    duration=6.0,
+    refl_coeff=0.991,
+    model_factory=room_modes_model_factory,
+    apply_loss=True,
+    diff=True,
+    image=True,
+    verbose=True,
+)

run_2d.sh

@@ -9,20 +9,21 @@ root_dir="$(cd "$(dirname "$0")" && pwd)"
 python_dir="$root_dir/src/python"
 engine_exe="$root_dir/$build_dir/src/cpp/main_2d/pffdtd_2d"
 
-sim_name="Diffusor"
+sim_name="Modes2D"
 sim_dir="$root_dir/data/sim_data/$sim_name/cpu"
-
-fmax=4000
-duration=0.050
+sim_setup="${sim_name}.py"
+model_dir="$root_dir/data/models/$sim_name"
 
 # Delete old sim
 rm -rf "$sim_dir"
 
 # Generate model
-python -m pffdtd.sim2d.fdtd --verbose --save --data_dir="$sim_dir" --duration="$duration" --fmax="$fmax"
+cd "$model_dir"
+python "$sim_setup"
 
 # Run sim
 DPCPP_CPU_PLACES=cores DPCPP_CPU_CU_AFFINITY=spread DPCPP_CPU_NUM_CUS=16 "$engine_exe" -s "$sim_dir"
+# pffdtd sim2d run --sim_dir "$sim_dir"
 
 # Report
 python -m pffdtd.sim2d.report --sim_dir="$sim_dir" "$sim_dir/out.h5"

src/cpp/pffdtd/simulation_2d.cpp

@@ -19,7 +19,7 @@ auto loadSimulation2D(std::filesystem::path const& dir, bool video)
   auto const videoRatio   = static_cast<double>(Ny) / static_cast<double>(Nx);
   auto const videoWidth   = std::min<size_t>(2000, static_cast<size_t>(Nx));
   auto const videoOptions = VideoWriter::Options{
-      .file      = dir.parent_path() / "out.avi",
+      .file      = dir / "out.avi",
       .width     = videoWidth,
       .height    = static_cast<size_t>(videoWidth * videoRatio),
       .fps       = sim.read<double>("video_fps"),

src/python/pffdtd/cli.py

@@ -2,6 +2,7 @@
 
 from pffdtd.analysis.cli import analysis
 from pffdtd.materials.cli import materials
+from pffdtd.sim2d.cli import sim2d
 
 
 @click.group()
@@ -14,3 +15,4 @@ def main(ctx, verbose):
 
 main.add_command(analysis)
 main.add_command(materials)
+main.add_command(sim2d)

src/python/pffdtd/common/myfuncs.py

@@ -20,17 +20,18 @@
 EPS = np.finfo(np.float64).eps
 
 
-def rotmatrix_ax_ang(Rax:Any, Rang:float):
-    assert isinstance(Rax,np.ndarray)
+def rotmatrix_ax_ang(Rax: Any, Rang: float):
+    assert isinstance(Rax, np.ndarray)
     assert Rax.shape == (3,)
     assert type(Rang) is float
 
-    Rax = Rax/npl.norm(Rax) #to be sure
+    Rax = Rax/npl.norm(Rax)  # to be sure
 
-    theta = Rang/180.0*pi #in rad
-    #see https://en.wikipedia.org/wiki/Rotation_matrix
-    R = np.array([[cos(theta) + Rax[0]*Rax[0]*(1-cos(theta)), Rax[0]*Rax[1]*(1-cos(theta)) - Rax[2]*sin(theta), Rax[0]*Rax[2]*(1-cos(theta)) + Rax[1]*sin(theta)],\
-                  [Rax[1]*Rax[0]*(1-cos(theta)) + Rax[2]*sin(theta), cos(theta) + Rax[1]*Rax[1]*(1-cos(theta)), Rax[1]*Rax[2]*(1-cos(theta)) - Rax[0]*sin(theta)],\
+    theta = Rang/180.0*pi  # in rad
+    # see https://en.wikipedia.org/wiki/Rotation_matrix
+    R = np.array([[cos(theta) + Rax[0]*Rax[0]*(1-cos(theta)), Rax[0]*Rax[1]*(1-cos(theta)) - Rax[2]*sin(theta), Rax[0]*Rax[2]*(1-cos(theta)) + Rax[1]*sin(theta)],
+                  [Rax[1]*Rax[0]*(1-cos(theta)) + Rax[2]*sin(theta), cos(theta) + Rax[1]*Rax[1]*(
+                      1-cos(theta)), Rax[1]*Rax[2]*(1-cos(theta)) - Rax[0]*sin(theta)],
                   [Rax[2]*Rax[0]*(1-cos(theta)) - Rax[1]*sin(theta), Rax[2]*Rax[1]*(1-cos(theta)) + Rax[0]*sin(theta), cos(theta) + Rax[2]*Rax[2]*(1-cos(theta))]])
     assert npl.norm(npl.inv(R)-R.T) < 1e-8
     return R
@@ -45,16 +46,16 @@ def rotate_xyz_deg(thx_d, thy_d, thz_d):
     thz = np.deg2rad(thz_d)
 
     Rx = np.array([[1, 0, 0],
-                    [0, np.cos(thx), -np.sin(thx)],
-                    [0, np.sin(thx), np.cos(thx)]])
+                   [0, np.cos(thx), -np.sin(thx)],
+                   [0, np.sin(thx), np.cos(thx)]])
 
     Ry = np.array([[np.cos(thy), 0, np.sin(thy)],
-                    [0, 1, 0],
-                    [-np.sin(thy), 0, np.cos(thy)]])
+                   [0, 1, 0],
+                   [-np.sin(thy), 0, np.cos(thy)]])
 
     Rz = np.array([[np.cos(thz), -np.sin(thz), 0],
-                    [np.sin(thz), np.cos(thz), 0],
-                    [0, 0, 1]])
+                   [np.sin(thz), np.cos(thz), 0],
+                   [0, 0, 1]])
 
     R = Rx @ Ry @ Rz
 
@@ -295,3 +296,27 @@ def wavwrite(fname, fs, data):
     # reads in (Nsamples,Nchannels)
     scipy.io.wavfile.write(fname, int(fs), np.float32(data.T))
     print(f'wrote {fname} at fs={fs/1000:.2f} kHz')
+
+
+def to_ixy(x, y, Nx, Ny, order="row"):
+    if order == "row":
+        return x*Ny+y
+    return y*Nx+x
+
+
+def point_on_circle(center, radius: float, angle: float):
+    """
+    Calculate the coordinates of a point on a circle arc.
+
+    Parameters:
+    center (tuple): (x, y) coordinates of the center of the circle.
+    radius: Radius of the circle.
+    angle: Angle in radians.
+
+    Returns:
+    tuple: (p_x, p_y) coordinates of the point on the circle arc.
+    """
+    x, y = center
+    p_x = x + radius * np.cos(angle)
+    p_y = y + radius * np.sin(angle)
+    return (p_x, p_y)

src/python/pffdtd/sim2d/cli.py

@@ -0,0 +1,17 @@
+import click
+
+from pffdtd.sim2d.engine import Engine2D
+
+
+@click.group(help="2D wave-equation.")
+def sim2d():
+    pass
+
+
+@sim2d.command(help="Run simulation.")
+@click.option('--sim_dir', type=click.Path(exists=True))
+@click.option('--video', is_flag=True)
+def run(sim_dir, video):
+    engine = Engine2D(sim_dir=sim_dir, video=video)
+    engine.run()
+    engine.save_output()