1"""Solve a 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
5 python run.py --problem cymbal --backend fmm3dbie --case small --selftest
7The problems are `room` (the default; a point source in a room with impedance
8walls) and `cymbal` (a listener near a rigid cymbal in free space). Results go
9to results/<case>-<backend>/ for the room and results/cymbal-<case>-<backend>/
10for the cymbal (summary.json, field.npz). The two backends solve the same
11discrete problem on the same surface nodes; compare.py puts their results side
12by side. `python run.py --problem cymbal --help` lists the cymbal's options.
13"""
15import argparse
16import os
17import sys
18import time
20import numpy as np
22HERE = os.path.dirname(os.path.abspath(__file__))
23sys.path.insert(0, HERE)
25import problems # noqa: E402
26from problems.common import impedance_data, machine_info, save, source_sum # noqa: E402
29class GPUBackend:
30 """The CuPy/CUDA solver in fmm3dbie_gpu."""
32 name = "gpu"
34 def __init__(self, pb, case, args):
35 import cupy as cp
36 from fmm3dbie_gpu import cuda_ops, solver
37 self.cp = self.xp = cp
38 t0 = time.time()
39 cuda_ops.module() # compile or load the CUDA kernels once, outside the timings
40 self.kernel_load_s = time.time() - t0
41 self.surf = pb.gpu_surface(case, args)
42 t0 = time.time()
43 self.op = solver.Operator(self.surf, case["k"], pb.ALPHA, case["lam"], ifinout=pb.IFINOUT,
44 far_dtype=np.float64 if args.far64 else np.float32, verbose=False)
45 cp.cuda.Device().synchronize()
46 self.setup_s = time.time() - t0
47 self.eps_gmres = args.eps_gmres
48 self.options = dict(far_field="float64" if args.far64 else "float32",
49 geometry="exact" if self.surf.exact is not None else "interpolated")
51 def solve(self, rhs):
52 sigma, phi1, niter, rres = self.op.solve(rhs, eps_gmres=self.eps_gmres, verbose=False)
53 return (sigma, phi1), niter, rres
55 def eval_points(self, state, pts):
56 return self.op.eval(pts.T, *state).get()
58 def eval_wall(self, state):
59 return self.op.eval(None, *state, on_surface=True).get()
61 def to_numpy(self, q):
62 return self.cp.asnumpy(q)
64 def sync(self):
65 self.cp.cuda.Device().synchronize()
67 def timed_solve(self, rhs):
68 return timed(self, self.solve, rhs)
71class Fmm3dbieBackend:
72 """fmm3dbie on the CPU (helm_rpcomb_imp_solver), through ctypes."""
74 name = "fmm3dbie"
76 def __init__(self, pb, case, args):
77 sys.path.insert(0, os.path.join(HERE, "reference"))
78 import fmm3dbie_ctypes as fb
79 self.fb = fb
80 self.xp = np
81 self.log = os.path.join(args.outdir, "fortran.log")
82 os.makedirs(args.outdir, exist_ok=True)
83 open(self.log, "w").close()
84 self.kernel_load_s = 0.0
85 with fb.fortran_output_to(self.log):
86 self.surf = pb.fmm3dbie_surface(case)
87 self.k, self.lam, self.alpha, self.ifinout = case["k"], case["lam"], pb.ALPHA, pb.IFINOUT
88 self.setup_s = 0.0 # the near quadrature is computed inside each solver call
89 self.eps, self.eps_gmres = args.eps, args.eps_gmres
90 self.options = dict(eps=args.eps, threads=os.environ.get("OMP_NUM_THREADS", "all"))
92 def solve(self, rhs):
93 with self.fb.fortran_output_to(self.log):
94 sigma, sik, niter, errs, rres = self.fb.helm_rpcomb_imp_solver(
95 self.surf, self.k, self.alpha, self.lam * np.ones(self.surf.npts), rhs,
96 ifinout=self.ifinout, eps=self.eps, eps_gmres=self.eps_gmres)
97 return (sigma, sik), int(niter), float(rres)
99 def eval_points(self, state, pts):
100 with self.fb.fortran_output_to(self.log):
101 return self.fb.helm_rpcomb_eval(self.surf, self.k, self.alpha, *state, pts, eps=self.eps)
103 def eval_wall(self, state):
104 with self.fb.fortran_output_to(self.log):
105 return self.fb.helm_rpcomb_eval_boundary(self.surf, self.k, self.alpha, *state,
106 eps=self.eps, ifinout=self.ifinout)
108 def to_numpy(self, q):
109 return q
111 def sync(self):
112 pass
114 def timed_solve(self, rhs):
115 return timed(self, self.solve, rhs)
118def timed(be, f, *a):
119 t0 = time.time()
120 out = f(*a)
121 be.sync()
122 return out, time.time() - t0
125def main():
126 pre = argparse.ArgumentParser(add_help=False)
127 pre.add_argument("--problem", choices=problems.PROBLEMS, default="room")
128 pb = problems.get(pre.parse_known_args()[0].problem)
130 ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
131 ap.add_argument("--problem", choices=problems.PROBLEMS, default="room")
132 ap.add_argument("--backend", choices=["gpu", "fmm3dbie"], required=True)
133 ap.add_argument("--case", choices=list(pb.CASES), default="small")
134 pb.add_arguments(ap)
135 ap.add_argument("--selftest", action="store_true",
136 help="also solve a problem with a known solution (a second solve)")
137 ap.add_argument("--eps-gmres", type=float, default=1e-6)
138 ap.add_argument("--eps", type=float, default=1e-4, help="fmm3dbie quadrature and FMM tolerance")
139 ap.add_argument("--far64", action="store_true", help="gpu: double precision far field")
140 ap.add_argument("--interp-geometry", action="store_true",
141 help="gpu: interpolated patch geometry inside the integrals, as fmm3dbie "
142 "(the cymbal always uses it)")
143 ap.add_argument("--ngrid", type=int, default=pb.NGRID,
144 help=f"grid size of the field slices (default {pb.NGRID}; 0 to skip)")
145 ap.add_argument("--outdir", help="default results/<case>-<backend> (room), "
146 "results/cymbal-<case>-<backend> (cymbal)")
147 args = ap.parse_args()
149 case = pb.make_case(args.case, args)
150 args.outdir = args.outdir or os.path.join(HERE, "results", pb.results_name(args.case, args.backend))
151 k, lam = case["k"], case["lam"]
153 be = (GPUBackend if args.backend == "gpu" else Fmm3dbieBackend)(pb, case, args)
154 surf = be.surf
155 xs, nrm = np.asarray(surf.srcvals)[0:3], np.asarray(surf.srcvals)[9:12]
156 print(f"[{args.backend}] {pb.NAME} {args.case}: {pb.describe(case)}; "
157 f"{surf.npatches} patches of order {case['norder']}, {surf.npts} points")
158 if args.backend == "gpu":
159 print(f" kernel load {be.kernel_load_s:.2f} s, operator setup {be.setup_s:.2f} s")
161 summary = dict(problem=pb.NAME, backend=args.backend, case=args.case,
162 **{key: val for key, val in case.items()}, npts=int(surf.npts),
163 eps_gmres=args.eps_gmres, options=be.options, machine=machine_info(),
164 kernel_load_s=be.kernel_load_s, setup_s=be.setup_s)
166 if args.selftest:
167 x0s, ws = pb.selftest_sources(case)
168 (state, niter, rres), t = be.timed_solve(source_sum(xs, x0s, ws, k, nrm=nrm, lam=lam, xp=be.xp))
169 xt = pb.selftest_points(case)
170 ut = be.eval_points(state, xt)
171 ub = be.eval_wall(state)
172 uex, ubex = source_sum(xt, x0s, ws, k, xp=be.xp), source_sum(xs, x0s, ws, k, xp=be.xp)
173 err = float(np.linalg.norm(ut - uex) / np.linalg.norm(uex))
174 err_wall = float(np.linalg.norm(ub - ubex) / np.linalg.norm(ubex))
175 print(f" self-test: {niter} iterations, {t:.2f} s; relative error {err:.2e} at "
176 f"{xt.shape[1]} points off the surface, {err_wall:.2e} on the surface")
177 summary.update(selftest_err=err, selftest_err_wall=err_wall, selftest_solve_s=t)
179 (state, niter, rres), t_solve = be.timed_solve(-impedance_data(xs, nrm, pb.source(case), k, lam))
180 print(f" main problem: {niter} GMRES iterations, residual {rres:.2e}, solve {t_solve:.2f} s")
181 summary.update(niter=niter, rres=rres, solve_s=t_solve, setup_plus_solve_s=be.setup_s + t_solve)
183 arrays = {}
184 pb.report(be, state, case, args, summary, arrays)
185 save(args.outdir, summary, **arrays)
186 print(f" saved {args.outdir}")
189if __name__ == "__main__":
190 main()