"""Solve the showcase problem with one backend and save timings and fields. python run.py --backend gpu --case small --selftest python run.py --backend fmm3dbie --case small --selftest Results go to results/-/ (summary.json, field.npz). The two backends solve the same discrete problem on the same surface nodes; compare.py puts their results side by side. """ import argparse import os import sys import time import numpy as np HERE = os.path.dirname(os.path.abspath(__file__)) sys.path.insert(0, HERE) import problem as pb # noqa: E402 class GPUBackend: """The CuPy/CUDA solver in fmm3dbie_gpu.""" name = "gpu" def __init__(self, case, args): import cupy as cp from fmm3dbie_gpu import cuda_ops, solver, surface self.cp = cp t0 = time.time() cuda_ops.module() # compile or load the CUDA kernels once, outside the timings self.kernel_load_s = time.time() - t0 self.surf = surface.rounded_cube(a=pb.A, p=pb.P, m=case["m"], norder=case["norder"], exact=not args.interp_geometry) self.lam = -case["beta"] t0 = time.time() self.op = solver.Operator(self.surf, case["k"], pb.ALPHA, self.lam, ifinout=0, far_dtype=np.float64 if args.far64 else np.float32, verbose=False) cp.cuda.Device().synchronize() self.setup_s = time.time() - t0 self.eps_gmres = args.eps_gmres self.options = dict(far_field="float64" if args.far64 else "float32", geometry="interpolated" if args.interp_geometry else "exact") def solve(self, rhs): sigma, phi1, niter, rres = self.op.solve(rhs, eps_gmres=self.eps_gmres, verbose=False) return (sigma, phi1), niter, rres def eval_points(self, state, pts): return self.op.eval(pts.T, *state).get() def eval_wall(self, state): return self.op.eval(None, *state, on_surface=True).get() def sync(self): self.cp.cuda.Device().synchronize() class Fmm3dbieBackend: """fmm3dbie on the CPU (helm_rpcomb_imp_solver), through ctypes.""" name = "fmm3dbie" def __init__(self, case, args): sys.path.insert(0, os.path.join(HERE, "reference")) import fmm3dbie_ctypes as fb import rounded_cube as rc self.fb = fb self.log = os.path.join(args.outdir, "fortran.log") os.makedirs(args.outdir, exist_ok=True) open(self.log, "w").close() self.kernel_load_s = 0.0 with fb.fortran_output_to(self.log): self.surf = rc.rounded_cube(a=pb.A, p=pb.P, m=case["m"], norder=case["norder"]) self.k, self.lam = case["k"], -case["beta"] self.setup_s = 0.0 # the near quadrature is computed inside each solver call self.eps, self.eps_gmres = args.eps, args.eps_gmres self.options = dict(eps=args.eps, threads=os.environ.get("OMP_NUM_THREADS", "all")) def solve(self, rhs): with self.fb.fortran_output_to(self.log): sigma, sik, niter, errs, rres = self.fb.helm_rpcomb_imp_solver( self.surf, self.k, pb.ALPHA, self.lam * np.ones(self.surf.npts), rhs, ifinout=0, eps=self.eps, eps_gmres=self.eps_gmres) return (sigma, sik), int(niter), float(rres) def eval_points(self, state, pts): with self.fb.fortran_output_to(self.log): return self.fb.helm_rpcomb_eval(self.surf, self.k, pb.ALPHA, *state, pts, eps=self.eps) def eval_wall(self, state): with self.fb.fortran_output_to(self.log): return self.fb.helm_rpcomb_eval_boundary(self.surf, self.k, pb.ALPHA, *state, eps=self.eps) def sync(self): pass def timed(be, f, *a): t0 = time.time() out = f(*a) be.sync() return out, time.time() - t0 def main(): ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) ap.add_argument("--backend", choices=["gpu", "fmm3dbie"], required=True) ap.add_argument("--case", choices=list(pb.CASES), default="small") ap.add_argument("--k", type=float, help="override the case's wavenumber") ap.add_argument("--beta", type=float, help="override the case's wall admittance") ap.add_argument("--m", type=int, help="override the patches per face edge") ap.add_argument("--norder", type=int, help="override the order per patch") ap.add_argument("--selftest", action="store_true", help="also solve a problem with a known solution (a second solve)") ap.add_argument("--eps-gmres", type=float, default=1e-6) ap.add_argument("--eps", type=float, default=1e-4, help="fmm3dbie quadrature and FMM tolerance") ap.add_argument("--far64", action="store_true", help="gpu: double precision far field") ap.add_argument("--interp-geometry", action="store_true", help="gpu: interpolated patch geometry inside the integrals, as fmm3dbie") ap.add_argument("--ngrid", type=int, default=101, help="grid size of the field slices (0 to skip)") ap.add_argument("--outdir", help="default results/-") args = ap.parse_args() case = dict(pb.CASES[args.case]) for key in ["k", "beta", "m", "norder"]: if getattr(args, key) is not None: case[key] = getattr(args, key) args.outdir = args.outdir or os.path.join(HERE, "results", f"{args.case}-{args.backend}") k, lam = case["k"], -case["beta"] be = (GPUBackend if args.backend == "gpu" else Fmm3dbieBackend)(case, args) surf = be.surf xs, nrm = np.asarray(surf.srcvals)[0:3], np.asarray(surf.srcvals)[9:12] wts = np.asarray(surf.wts) print(f"[{args.backend}] case {args.case}: k = {k}, beta = {case['beta']}, " f"{surf.npatches} patches of order {case['norder']}, {surf.npts} points, " f"room width {2 * pb.A * k / (2 * np.pi):.2f} wavelengths") if args.backend == "gpu": print(f" kernel load {be.kernel_load_s:.2f} s, operator setup {be.setup_s:.2f} s") summary = dict(backend=args.backend, case=args.case, **case, npts=int(surf.npts), eps_gmres=args.eps_gmres, options=be.options, machine=pb.machine_info(), kernel_load_s=be.kernel_load_s, setup_s=be.setup_s) if args.selftest: (state, niter, rres), t = timed(be, be.solve, pb.selftest_rhs(xs, nrm, k, lam)) xt = pb.selftest_points() ut = be.eval_points(state, xt) ub = be.eval_wall(state) uex, ubex = pb.green(xt, pb.EXTERIOR_SOURCE, k), pb.green(xs, pb.EXTERIOR_SOURCE, k) err = float(np.linalg.norm(ut - uex) / np.linalg.norm(uex)) err_wall = float(np.linalg.norm(ub - ubex) / np.linalg.norm(ubex)) print(f" self-test: {niter} iterations, {t:.2f} s; relative error {err:.2e} at 50 interior " f"points, {err_wall:.2e} on the walls") summary.update(selftest_err=err, selftest_err_wall=err_wall, selftest_solve_s=t) (state, niter, rres), t_solve = timed(be, be.solve, pb.room_rhs(xs, nrm, k, lam)) print(f" room problem: {niter} GMRES iterations, residual {rres:.2e}, solve {t_solve:.2f} s") v0 = complex(be.eval_points(state, pb.SOURCE[:, None])[0]) u_wall = pb.green(xs, pb.SOURCE, k) + be.eval_wall(state) p_src, p_wall = pb.energy_balance(k, lam, v0, wts, u_wall) mismatch = abs(p_src - p_wall) / p_src print(f" energy balance: source {p_src:.8f}, walls {p_wall:.8f}, relative mismatch {mismatch:.2e}") summary.update(niter=niter, rres=rres, solve_s=t_solve, setup_plus_solve_s=be.setup_s + t_solve, p_src=p_src, p_wall=p_wall, energy_mismatch=mismatch) arrays = dict(x0=pb.SOURCE, k=k, lam=lam, a=pb.A, p=pb.P, xs=xs, wts=wts, u_wall=u_wall) if args.ngrid > 0: g, sl = pb.slice_points(args.ngrid) t0 = time.time() for name, (pts, mask, shape) in sl.items(): u = np.full(mask.size, np.nan + 1j * np.nan) u[mask] = pb.green(pts[:, mask], pb.SOURCE, k) + be.eval_points(state, pts[:, mask]) arrays[f"u_{name}"] = u.reshape(shape) be.sync() summary["slices_s"] = time.time() - t0 arrays["grid"] = g print(f" field on three {args.ngrid} x {args.ngrid} slices: {summary['slices_s']:.2f} s") pb.save(args.outdir, summary, **arrays) print(f" saved {args.outdir}") if __name__ == "__main__": main()