"""Solve a 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 python run.py --problem cymbal --backend fmm3dbie --case small --selftest The problems are `room` (the default; a point source in a room with impedance walls) and `cymbal` (a listener near a rigid cymbal in free space). Results go to results/-/ for the room and results/cymbal--/ for the cymbal (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. `python run.py --problem cymbal --help` lists the cymbal's options. """ 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 problems # noqa: E402 from problems.common import impedance_data, machine_info, save, source_sum # noqa: E402 class GPUBackend: """The CuPy/CUDA solver in fmm3dbie_gpu.""" name = "gpu" def __init__(self, pb, case, args): import cupy as cp from fmm3dbie_gpu import cuda_ops, solver self.cp = self.xp = 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 = pb.gpu_surface(case, args) t0 = time.time() self.op = solver.Operator(self.surf, case["k"], pb.ALPHA, case["lam"], ifinout=pb.IFINOUT, 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="exact" if self.surf.exact is not None else "interpolated") 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 to_numpy(self, q): return self.cp.asnumpy(q) def sync(self): self.cp.cuda.Device().synchronize() def timed_solve(self, rhs): return timed(self, self.solve, rhs) class Fmm3dbieBackend: """fmm3dbie on the CPU (helm_rpcomb_imp_solver), through ctypes.""" name = "fmm3dbie" def __init__(self, pb, case, args): sys.path.insert(0, os.path.join(HERE, "reference")) import fmm3dbie_ctypes as fb self.fb = fb self.xp = np 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 = pb.fmm3dbie_surface(case) self.k, self.lam, self.alpha, self.ifinout = case["k"], case["lam"], pb.ALPHA, pb.IFINOUT 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, self.alpha, self.lam * np.ones(self.surf.npts), rhs, ifinout=self.ifinout, 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, self.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, self.alpha, *state, eps=self.eps, ifinout=self.ifinout) def to_numpy(self, q): return q def sync(self): pass def timed_solve(self, rhs): return timed(self, self.solve, rhs) def timed(be, f, *a): t0 = time.time() out = f(*a) be.sync() return out, time.time() - t0 def main(): pre = argparse.ArgumentParser(add_help=False) pre.add_argument("--problem", choices=problems.PROBLEMS, default="room") pb = problems.get(pre.parse_known_args()[0].problem) ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) ap.add_argument("--problem", choices=problems.PROBLEMS, default="room") ap.add_argument("--backend", choices=["gpu", "fmm3dbie"], required=True) ap.add_argument("--case", choices=list(pb.CASES), default="small") pb.add_arguments(ap) 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 " "(the cymbal always uses it)") ap.add_argument("--ngrid", type=int, default=pb.NGRID, help=f"grid size of the field slices (default {pb.NGRID}; 0 to skip)") ap.add_argument("--outdir", help="default results/- (room), " "results/cymbal-- (cymbal)") args = ap.parse_args() case = pb.make_case(args.case, args) args.outdir = args.outdir or os.path.join(HERE, "results", pb.results_name(args.case, args.backend)) k, lam = case["k"], case["lam"] be = (GPUBackend if args.backend == "gpu" else Fmm3dbieBackend)(pb, case, args) surf = be.surf xs, nrm = np.asarray(surf.srcvals)[0:3], np.asarray(surf.srcvals)[9:12] print(f"[{args.backend}] {pb.NAME} {args.case}: {pb.describe(case)}; " f"{surf.npatches} patches of order {case['norder']}, {surf.npts} points") if args.backend == "gpu": print(f" kernel load {be.kernel_load_s:.2f} s, operator setup {be.setup_s:.2f} s") summary = dict(problem=pb.NAME, backend=args.backend, case=args.case, **{key: val for key, val in case.items()}, npts=int(surf.npts), eps_gmres=args.eps_gmres, options=be.options, machine=machine_info(), kernel_load_s=be.kernel_load_s, setup_s=be.setup_s) if args.selftest: x0s, ws = pb.selftest_sources(case) (state, niter, rres), t = be.timed_solve(source_sum(xs, x0s, ws, k, nrm=nrm, lam=lam, xp=be.xp)) xt = pb.selftest_points(case) ut = be.eval_points(state, xt) ub = be.eval_wall(state) uex, ubex = source_sum(xt, x0s, ws, k, xp=be.xp), source_sum(xs, x0s, ws, k, xp=be.xp) 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 " f"{xt.shape[1]} points off the surface, {err_wall:.2e} on the surface") summary.update(selftest_err=err, selftest_err_wall=err_wall, selftest_solve_s=t) (state, niter, rres), t_solve = be.timed_solve(-impedance_data(xs, nrm, pb.source(case), k, lam)) print(f" main problem: {niter} GMRES iterations, residual {rres:.2e}, solve {t_solve:.2f} s") summary.update(niter=niter, rres=rres, solve_s=t_solve, setup_plus_solve_s=be.setup_s + t_solve) arrays = {} pb.report(be, state, case, args, summary, arrays) save(args.outdir, summary, **arrays) print(f" saved {args.outdir}") if __name__ == "__main__": main()