From 288719e1e649c11c6037e20df46260220ef9dbc2 Mon Sep 17 00:00:00 2001 From: Clawrence Date: Tue, 14 Jul 2026 10:37:06 -0400 Subject: [PATCH] bench: add remaining test-wave probe scripts (PML sweep, solver frontier, wave driver) Completes reproducibility for TESTWAVE-2026-07-14.md: pml_bench.py (Mie-gated PML thickness), probe4.py (pc/solver frontier probe), wave.sh (phase-split + batching + basix-count + PML driver). Co-Authored-By: Claude Fable 5 --- bench/pml_bench.py | 43 ++++++++++++++++++++++++++++++++++ bench/probe4.py | 50 +++++++++++++++++++++++++++++++++++++++ bench/wave.sh | 58 ++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 151 insertions(+) create mode 100644 bench/pml_bench.py create mode 100644 bench/probe4.py create mode 100644 bench/wave.sh diff --git a/bench/pml_bench.py b/bench/pml_bench.py new file mode 100644 index 0000000..3c33ba0 --- /dev/null +++ b/bench/pml_bench.py @@ -0,0 +1,43 @@ +# PML-thickness economy probe: friend's testSphereScattering geometry (Mie ground +# truth), deg 2, sweep PML_thickness; report ndofs, factor memory, far-field error. +# Usage: python3 pml_bench.py +import os +os.environ["OMP_NUM_THREADS"] = "1" +import matplotlib +matplotlib.use("Agg") +import sys +import numpy as np +import dolfinx.fem.petsc +from mpi4py import MPI +from timeit import default_timer as timer + +sys.path.insert(0, os.environ.get("SCATT3D_SRC", "/work/Scatt3D")) +import meshMaker +import scatteringProblem + +comm = MPI.COMM_WORLD +T = float(sys.argv[1]) +h = 1 / 12 +deg = 2 + +os.makedirs("data3D", exist_ok=True) +t0 = timer() +refMesh = meshMaker.MeshInfo( + comm, reference=True, viewGMSH=False, verbosity=2, N_antennas=0, + object_radius=.33, domain_radius=.9, PML_thickness=T, h=h, + domain_geom="sphere", object_geom="sphere", FF_surface=True, order=deg) +freqs = np.linspace(10e9, 12e9, 1) +prob = scatteringProblem.Scatt3DProblem( + comm, refMesh, verbosity=2, name=f"pml{T}", dataFolder="data3D/", MPInum=comm.size, + makeOptVects=True, excitation="planewave", freqs=freqs, + material_epsrs=[2.0 * (1 - 0.01j)], fem_degree=deg, + solver_settings={"symmetric": True}) +t1 = timer() +# his own Mie-comparison path prints 'Forward-scattering intensity relative error: ...' +prob.calcFarField(reference=True, compareToMie=True, showPlots=False, + returnConvergenceVals=False) +nd = getattr(prob, "ndofs", None) or getattr(getattr(prob, "FEMmesh_ref", None), "ndofs", "?") +fmem = getattr(prob, "lastFactorMemMB", None) +if comm.rank == 0: + print(f"PML_RESULT T={T} ndofs={nd} solve_s={t1-t0:.1f} " + f"mumps_mem_mb={list(fmem) if fmem else None}") diff --git a/bench/probe4.py b/bench/probe4.py new file mode 100644 index 0000000..494bf51 --- /dev/null +++ b/bench/probe4.py @@ -0,0 +1,50 @@ +# Direct-solver frontier probe: factor the dumped 545k coax matrix with a chosen +# PC type + factor solver package; report factor time, peak RSS, MUMPS INFOG when +# applicable, and forward-solve accuracy. Usage: probe4.py [petsc opts] +import sys +import resource +import petsc4py +petsc4py.init(sys.argv) +from petsc4py import PETSc +import numpy as np +from timeit import default_timer as timer + +pctype, solver = sys.argv[1], sys.argv[2] +comm = PETSc.COMM_WORLD +viewer = PETSc.Viewer().createBinary("/work/bench/cableport/A545k.bin", "r", comm=comm) +A = PETSc.Mat().load(viewer) +viewer.destroy() +A.assemble() +n = A.getSize()[0] +A.setOption(PETSc.Mat.Option.SYMMETRIC, True) +A.setOption(PETSc.Mat.Option.SYMMETRY_ETERNAL, True) +rss0 = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / 1024 + +rng = np.random.default_rng(0) +x_ex = rng.standard_normal(n) + 1j * rng.standard_normal(n) +xv = PETSc.Vec().createWithArray(x_ex, comm=comm) +b = A.createVecLeft() +A.mult(xv, b) +x = A.createVecRight() + +ksp = PETSc.KSP().create(comm) +ksp.setOperators(A) +ksp.setType("preonly") +pc = ksp.getPC() +pc.setType(pctype) +pc.setFactorSolverType(solver) +ksp.setFromOptions() +t0 = timer() +ksp.setUp() +tf = timer() - t0 +rss1 = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / 1024 +infog = "" +try: + F = pc.getFactorMatrix() + infog = (f"INFOG21={F.getMumpsInfog(21)} INFOG22={F.getMumpsInfog(22)}") +except Exception: + pass +ksp.solve(b, x) +err = (x - xv).norm() / xv.norm() +print(f"PROBE4_RESULT pc={pctype} solver={solver} n={n} factor_s={tf:.1f} " + f"rss_before_mb={rss0:.0f} rss_peak_mb={rss1:.0f} {infog} rel_err={err:.2e}") diff --git a/bench/wave.sh b/bench/wave.sh new file mode 100644 index 0000000..b9b0873 --- /dev/null +++ b/bench/wave.sh @@ -0,0 +1,58 @@ +#!/bin/bash +# Wave-2 fast items: (a) log_view phase split / symbolic-reuse check, +# (b) matrix dump + MatMatSolve batching probe, (c) basix condensable-DOF count, +# (d) PML-thickness economy sweep with Mie accuracy gate. +exec > /work/bench/wave.log 2>&1 +source /usr/local/bin/dolfinx-complex-mode 2>/dev/null +cd /work/bench +mkdir -p wave_work && cd wave_work +cp ../bench_driver.py . + +echo "=== (a) LOGVIEW sym h=0.42 Nf=6 np=8 ($(date -u +%H:%M:%SZ)) ===" +PETSC_OPTIONS="-log_view :/work/bench/wave_work/logview_sym.txt" \ + SCATT3D_SRC=/work/Scatt3D timeout 7200 mpirun -np 8 python3 bench_driver.py 0.42 3 6 '{"symmetric": true}' PHASES > case_phases.log 2>&1 +echo "EXIT $?" +grep -E 'BENCH_RESULT|MUMPS factor memory' case_phases.log | head -3 +grep -E 'MatCholFctrSym|MatCholFctrNum|MatSolve |MatAssemblyEnd|SNESSolve|KSPSolve ' /work/bench/wave_work/logview_sym.txt | head -8 + +echo "=== (b) matrix dump for MatMatSolve probe ($(date -u +%H:%M:%SZ)) ===" +mkdir -p /work/bench/cableport +if [ ! -f /work/bench/cableport/A545k.bin ]; then + PETSC_OPTIONS="-theScatteringProblem_ksp_view_mat binary:/work/bench/cableport/A545k.bin" \ + SCATT3D_SRC=/work/Scatt3D timeout 3600 python3 ../cableport_validate.py 3 '{}' dump 0.2857142857142857 3 1 > case_dump.log 2>&1 + echo "dump EXIT $?" +fi +ls -la /work/bench/cableport/ +echo "=== (b) MATMAT probe ($(date -u +%H:%M:%SZ)) ===" +timeout 3600 python3 ../matmat_bench.py > case_matmat.log 2>&1 +echo "EXIT $?" +grep MATMAT case_matmat.log + +echo "=== (c) basix deg-3 N1curl condensable-DOF count ===" +python3 - <<'PY' +import basix +import numpy as np +import dolfinx.mesh as dm +from mpi4py import MPI +e = basix.create_element(basix.ElementFamily.N1E, basix.CellType.tetrahedron, 3) +ned = e.num_entity_dofs +per_edge = ned[1][0]; per_face = ned[2][0]; per_cell = ned[3][0] +print("BASIX per-entity dofs deg3 N1E: edge", per_edge, "face", per_face, "cell-interior", per_cell) +msh = dm.create_unit_cube(MPI.COMM_WORLD, 20, 20, 20, dm.CellType.tetrahedron) +msh.topology.create_entities(1); msh.topology.create_entities(2) +nE = msh.topology.index_map(1).size_global; nF = msh.topology.index_map(2).size_global +nC = msh.topology.index_map(3).size_global +tot = nE*per_edge + nF*per_face + nC*per_cell +print(f"BASIX mesh E={nE} F={nF} C={nC} -> dofs total={tot}, " + f"cell-interior {nC*per_cell} ({100*nC*per_cell/tot:.1f}%), " + f"face {nF*per_face} ({100*nF*per_face/tot:.1f}%), edge {nE*per_edge} ({100*nE*per_edge/tot:.1f}%)") +PY + +echo "=== (d) PML thickness sweep, Mie gate ($(date -u +%H:%M:%SZ)) ===" +for T in 0.5 0.35 0.25 0.15 0.1; do + SCATT3D_SRC=/work/Scatt3D timeout 7200 mpirun -np 8 python3 ../pml_bench.py $T > case_pml$T.log 2>&1 + echo "PML $T EXIT $?" + grep -E 'PML_RESULT|Traceback' case_pml$T.log | head -3 + rm -rf data3D/pml* +done +echo WAVE_DONE