1"""Solve the showcase problem with one backend and save timings and fields.
3 python run.py --backend gpu --case small --selftest
4 python run.py --backend fmm3dbie --case small --selftest
6Results go to results/<case>-<backend>/ (summary.json, field.npz). The two
7backends solve the same discrete problem on the same surface nodes; compare.py
8puts their results side by side.
9"""
11import argparse
12import os
13import sys
14import time
16import numpy as np
18HERE = os.path.dirname(os.path.abspath(__file__))
19sys.path.insert(0, HERE)
21import problem as pb # noqa: E402
24class GPUBackend:
25 """The CuPy/CUDA solver in fmm3dbie_gpu."""
27 name = "gpu"
29 def __init__(self, case, args):
30 import cupy as cp
31 from fmm3dbie_gpu import cuda_ops, solver, surface
32 self.cp = cp
33 t0 = time.time()
34 cuda_ops.module() # compile or load the CUDA kernels once, outside the timings
35 self.kernel_load_s = time.time() - t0
36 self.surf = surface.rounded_cube(a=pb.A, p=pb.P, m=case["m"], norder=case["norder"],
37 exact=not args.interp_geometry)
38 self.lam = -case["beta"]
39 t0 = time.time()
40 self.op = solver.Operator(self.surf, case["k"], pb.ALPHA, self.lam, ifinout=0,
41 far_dtype=np.float64 if args.far64 else np.float32, verbose=False)
42 cp.cuda.Device().synchronize()
43 self.setup_s = time.time() - t0
44 self.eps_gmres = args.eps_gmres
45 self.options = dict(far_field="float64" if args.far64 else "float32",
46 geometry="interpolated" if args.interp_geometry else "exact")
48 def solve(self, rhs):
49 sigma, phi1, niter, rres = self.op.solve(rhs, eps_gmres=self.eps_gmres, verbose=False)
50 return (sigma, phi1), niter, rres
52 def eval_points(self, state, pts):
53 return self.op.eval(pts.T, *state).get()
55 def eval_wall(self, state):
56 return self.op.eval(None, *state, on_surface=True).get()
58 def sync(self):
59 self.cp.cuda.Device().synchronize()
62class Fmm3dbieBackend:
63 """fmm3dbie on the CPU (helm_rpcomb_imp_solver), through ctypes."""
65 name = "fmm3dbie"
67 def __init__(self, case, args):
68 sys.path.insert(0, os.path.join(HERE, "reference"))
69 import fmm3dbie_ctypes as fb
70 import rounded_cube as rc
71 self.fb = fb
72 self.log = os.path.join(args.outdir, "fortran.log")
73 os.makedirs(args.outdir, exist_ok=True)
74 open(self.log, "w").close()
75 self.kernel_load_s = 0.0
76 with fb.fortran_output_to(self.log):
77 self.surf = rc.rounded_cube(a=pb.A, p=pb.P, m=case["m"], norder=case["norder"])
78 self.k, self.lam = case["k"], -case["beta"]
79 self.setup_s = 0.0 # the near quadrature is computed inside each solver call
80 self.eps, self.eps_gmres = args.eps, args.eps_gmres
81 self.options = dict(eps=args.eps, threads=os.environ.get("OMP_NUM_THREADS", "all"))
83 def solve(self, rhs):
84 with self.fb.fortran_output_to(self.log):
85 sigma, sik, niter, errs, rres = self.fb.helm_rpcomb_imp_solver(
86 self.surf, self.k, pb.ALPHA, self.lam * np.ones(self.surf.npts), rhs, ifinout=0,
87 eps=self.eps, eps_gmres=self.eps_gmres)
88 return (sigma, sik), int(niter), float(rres)
90 def eval_points(self, state, pts):
91 with self.fb.fortran_output_to(self.log):
92 return self.fb.helm_rpcomb_eval(self.surf, self.k, pb.ALPHA, *state, pts, eps=self.eps)
94 def eval_wall(self, state):
95 with self.fb.fortran_output_to(self.log):
96 return self.fb.helm_rpcomb_eval_boundary(self.surf, self.k, pb.ALPHA, *state, eps=self.eps)
98 def sync(self):
99 pass
102def timed(be, f, *a):
103 t0 = time.time()
104 out = f(*a)
105 be.sync()
106 return out, time.time() - t0
109def main():
110 ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
111 ap.add_argument("--backend", choices=["gpu", "fmm3dbie"], required=True)
112 ap.add_argument("--case", choices=list(pb.CASES), default="small")
113 ap.add_argument("--k", type=float, help="override the case's wavenumber")
114 ap.add_argument("--beta", type=float, help="override the case's wall admittance")
115 ap.add_argument("--m", type=int, help="override the patches per face edge")
116 ap.add_argument("--norder", type=int, help="override the order per patch")
117 ap.add_argument("--selftest", action="store_true",
118 help="also solve a problem with a known solution (a second solve)")
119 ap.add_argument("--eps-gmres", type=float, default=1e-6)
120 ap.add_argument("--eps", type=float, default=1e-4, help="fmm3dbie quadrature and FMM tolerance")
121 ap.add_argument("--far64", action="store_true", help="gpu: double precision far field")
122 ap.add_argument("--interp-geometry", action="store_true",
123 help="gpu: interpolated patch geometry inside the integrals, as fmm3dbie")
124 ap.add_argument("--ngrid", type=int, default=101, help="grid size of the field slices (0 to skip)")
125 ap.add_argument("--outdir", help="default results/<case>-<backend>")
126 args = ap.parse_args()
128 case = dict(pb.CASES[args.case])
129 for key in ["k", "beta", "m", "norder"]:
130 if getattr(args, key) is not None:
131 case[key] = getattr(args, key)
132 args.outdir = args.outdir or os.path.join(HERE, "results", f"{args.case}-{args.backend}")
133 k, lam = case["k"], -case["beta"]
135 be = (GPUBackend if args.backend == "gpu" else Fmm3dbieBackend)(case, args)
136 surf = be.surf
137 xs, nrm = np.asarray(surf.srcvals)[0:3], np.asarray(surf.srcvals)[9:12]
138 wts = np.asarray(surf.wts)
139 print(f"[{args.backend}] case {args.case}: k = {k}, beta = {case['beta']}, "
140 f"{surf.npatches} patches of order {case['norder']}, {surf.npts} points, "
141 f"room width {2 * pb.A * k / (2 * np.pi):.2f} wavelengths")
142 if args.backend == "gpu":
143 print(f" kernel load {be.kernel_load_s:.2f} s, operator setup {be.setup_s:.2f} s")
145 summary = dict(backend=args.backend, case=args.case, **case, npts=int(surf.npts),
146 eps_gmres=args.eps_gmres, options=be.options, machine=pb.machine_info(),
147 kernel_load_s=be.kernel_load_s, setup_s=be.setup_s)
149 if args.selftest:
150 (state, niter, rres), t = timed(be, be.solve, pb.selftest_rhs(xs, nrm, k, lam))
151 xt = pb.selftest_points()
152 ut = be.eval_points(state, xt)
153 ub = be.eval_wall(state)
154 uex, ubex = pb.green(xt, pb.EXTERIOR_SOURCE, k), pb.green(xs, pb.EXTERIOR_SOURCE, k)
155 err = float(np.linalg.norm(ut - uex) / np.linalg.norm(uex))
156 err_wall = float(np.linalg.norm(ub - ubex) / np.linalg.norm(ubex))
157 print(f" self-test: {niter} iterations, {t:.2f} s; relative error {err:.2e} at 50 interior "
158 f"points, {err_wall:.2e} on the walls")
159 summary.update(selftest_err=err, selftest_err_wall=err_wall, selftest_solve_s=t)
161 (state, niter, rres), t_solve = timed(be, be.solve, pb.room_rhs(xs, nrm, k, lam))
162 print(f" room problem: {niter} GMRES iterations, residual {rres:.2e}, solve {t_solve:.2f} s")
163 v0 = complex(be.eval_points(state, pb.SOURCE[:, None])[0])
164 u_wall = pb.green(xs, pb.SOURCE, k) + be.eval_wall(state)
165 p_src, p_wall = pb.energy_balance(k, lam, v0, wts, u_wall)
166 mismatch = abs(p_src - p_wall) / p_src
167 print(f" energy balance: source {p_src:.8f}, walls {p_wall:.8f}, relative mismatch {mismatch:.2e}")
168 summary.update(niter=niter, rres=rres, solve_s=t_solve, setup_plus_solve_s=be.setup_s + t_solve,
169 p_src=p_src, p_wall=p_wall, energy_mismatch=mismatch)
171 arrays = dict(x0=pb.SOURCE, k=k, lam=lam, a=pb.A, p=pb.P, xs=xs, wts=wts, u_wall=u_wall)
172 if args.ngrid > 0:
173 g, sl = pb.slice_points(args.ngrid)
174 t0 = time.time()
175 for name, (pts, mask, shape) in sl.items():
176 u = np.full(mask.size, np.nan + 1j * np.nan)
177 u[mask] = pb.green(pts[:, mask], pb.SOURCE, k) + be.eval_points(state, pts[:, mask])
178 arrays[f"u_{name}"] = u.reshape(shape)
179 be.sync()
180 summary["slices_s"] = time.time() - t0
181 arrays["grid"] = g
182 print(f" field on three {args.ngrid} x {args.ngrid} slices: {summary['slices_s']:.2f} s")
183 pb.save(args.outdir, summary, **arrays)
184 print(f" saved {args.outdir}")
187if __name__ == "__main__":
188 main()