bench: add remaining test-wave probe scripts

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 <PML_thickness>
+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}")

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 <pctype> <solver> [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}")

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