"""Automated grader for the cubic-Hermite FEM (Euler-Bernoulli beam) benchmark. Usage: python grade.py [--budget SECONDS] [--no-exec] Required files in the output dir: fem_solver.py submitted solver (REQUIRED; imported and executed) reference_arrays.npz numerical arrays results.json summary sidecar Scoring: 1.0 Tier 1 + Tier 2 + Tier 3 all pass 0.5 Tier 1 + Tier 2 (and all provenance checks) pass, Tier 3 missing/fails 0.0 required files missing, schema/consistency violations, wrong mesh sequence, disallowed dependencies, a solver that cannot reproduce its own submitted numbers, or Tier 1 / Tier 2 fail Design notes ------------ The grader is self-contained: it knows the manufactured exact solution, recomputes it analytically, recomputes every order of convergence from the submitted arrays, and re-runs the submitted solver. *Method-specific thresholds are legitimate here.* Unlike the sibling ADI task, this specification fixes the method completely -- cubic-Hermite C^1 elements, a given manufactured load, a given mesh sequence, a given dense evaluation grid. There is therefore one correct answer, and a band around the gold reference's dense errors is a statement about the *problem*, not a preference for one of several admissible schemes. The bands are kept deliberately tight because the signature failure mode -- a dense reconstruction that omits the Hermite slope DOFs -- produces perfect nodal values and an O(h) dense error, which only the dense band and the dense order check can see. *Provenance.* Submitted arrays are worthless unless they came from the submitted code, so fem_solver.py is imported in an isolated subprocess and used to regenerate the study. ALL SEVEN meshes are mandatory within a configurable time budget: n = 32, 64, 128, 256 (n = 256 is the artifact mesh, so this also re-derives the Tier 3 field) and n = 512, 1024, 2048. Nothing is best-effort -- the reference solves the three largest in 0.11 s combined, and n = 2048 carries the Tier 2 round-off-floor requirement, so those rows cannot be accepted unverified. The runner also observes the required 2(n-1)-dimensional global solve and poisons solve results at n=32 to establish that the returned field causally depends on the solve rather than on a discarded dummy computation. """ import argparse import ast import json import os import subprocess import sys import tempfile import numpy as np # ---- problem ground truth (must match fem_solver.py) ------------------------ N_VALUES = [2**p for p in range(5, 12)] # 32 .. 2048 DENSE_FACTOR = 10 # dense grid is 10x finer CLEAN = 4 # n=32,64,128,256 are pre-round-off ARTIFACT_N = 256 ARTIFACT_POINTS = DENSE_FACTOR * ARTIFACT_N + 1 # 2561 ARTIFACT_INDEX = N_VALUES.index(ARTIFACT_N) # 3 # ---- acceptance thresholds -------------------------------------------------- # Gold reference dense errors in the clean (pre-round-off) regime, verified. REF_DENSE_CLEAN = [5.402330e-07, 3.459441e-08, 2.188513e-09, 1.376110e-10] REL_BAND = 0.25 # +/-25% band on the coarse dense errors ORDER_LO, ORDER_HI = 3.8, 4.2 NODAL_SUPER_ABS = 1e-8 # nodal e_n at n=32 must be ~round-off small NODAL_SUPER_RATIO = 100.0 # ... and << dense error there ROUNDOFF_FACTOR = 10.0 # finest dense error must exceed n=256 (round-off floor) ARTIFACT_MAX_ERR = 1.0e-9 # reference 1.38e-10 # The dense evaluation grid is prescribed exactly: x_dense[i] = i/(10n). # 1e-15 is ~4.5 ULP at magnitude 1, i.e. "the exact grid up to double-precision # representation" -- it admits any sane way of forming i/(10n) but rejects any # other point set (shifted, Chebyshev, 10n points instead of 10n+1, ...). GRID_ATOL = 1e-15 # Provenance tolerances. REPRO_RTOL = 1.0e-6 # regenerated error scalars vs submitted # Pointwise agreement between the submitted n=256 field and the one the solver # produces when re-run. Re-running the same code normally reproduces it bit for # bit; a materially different BLAS/LAPACK could shift it by round-off (~1e-12 at # this conditioning), so a deviation in the grey band is reported rather than # failed, and Tier 3's environment-independent cross-run check decides. Anything # above the artifact accuracy cap is not round-off and is failed outright. ARTIFACT_FIELD_ATOL = 1e-11 # silent pass below this ARTIFACT_FIELD_FAIL = 1e-9 # hard fail above this (= ARTIFACT_MAX_ERR) DEFAULT_BUDGET = 600.0 # seconds per solver-regeneration phase COARSE_N = N_VALUES[:CLEAN] # 32, 64, 128, 256 -- mandatory FINE_N = N_VALUES[CLEAN:] # 512, 1024, 2048 -- also mandatory ALLOWED_TOP_LEVEL = set(sys.stdlib_module_names) | {"numpy", "scipy"} # ---- independent implementation, used to verify the returned FEM field ---- # Reproducing one's own numbers proves only that the arrays came from the # submitted code -- not that the code does any finite element work. A solver # returning u_exact plus a synthetic O(h^4) error reproduces its own arrays # exactly, lands in the dense band if it knows the magnitudes, shows order 4, # and fakes the round-off floor. The only robust answer is to solve the problem # here, independently, and compare fields. For this particular problem the # Ritz solution can agree to round-off with the exact cubic-Hermite interpolant, # so field agreement is paired below with observed-system and causal solve audits. # # The method is fully specified (C1 cubic Hermite, given load, >= 3-point # Gauss), so a tight pointwise comparison is legitimate -- unlike the sibling ADI # task, which is scheme-neutral and needs a hidden instance instead. # # Calibration: at n = 32 the genuine FEM error is 5.40e-07, while 3-point and # 5-point Gauss implementations differ by only 4.51e-12 (5-point vs 8-point: # 7e-16). The tolerance below sits 200x above quadrature noise and 540x below # the FEM error itself. # EVERY mesh is checked, not a cheap subset: a solver doing genuine FEM only on # the checked meshes and fabricating the rest passes every other test (the dense # band, the orders, the superconvergence, the round-off floor and the artifact). # All seven cost 0.30 s combined in the grader's own implementation. FIDELITY_N = tuple(N_VALUES) # One tolerance suffices at every mesh: gold-vs-grader deviation is 4.5e-12 at # n = 32 and <= 1.1e-14 beyond, while the method's own error ranges 5.4e-07 down # to 1.4e-10 and back up to 1.2e-07 in the round-off regime -- four to seven # orders of separation throughout. FIDELITY_ATOL = 1.0e-9 # u_nodes and u_dense come from the SAME degrees of freedom, so at the coincident # points (every DENSE_FACTOR-th dense sample) they are the same evaluation. The # gold solver agrees to exactly 0.0 at all seven meshes; 1e-12 leaves room for a # different-but-valid evaluation order without admitting substituted values. # Without this, a solver can compute an honest dense field and still hand back # u_exact as its nodal values, faking perfect superconvergence. NODAL_DENSE_ATOL = 1.0e-12 def _shape(xi, h): """Cubic Hermite shapes on the reference element; slope DOFs are true derivatives (the slope functions carry the factor h).""" return np.array([2 * xi**3 - 3 * xi**2 + 1, h * (xi**3 - 2 * xi**2 + xi), -2 * xi**3 + 3 * xi**2, h * (xi**3 - xi**2)]) def grader_fem(n): """The grader's own cubic-Hermite solve, returning the dense field. Deliberately uses 5-point Gauss-Legendre where the bundled reference uses 3-point, so agreement is evidence about the METHOD rather than about a shared implementation detail. """ from numpy.polynomial.legendre import leggauss from scipy.sparse import lil_matrix, csr_matrix from scipy.sparse.linalg import spsolve xg, wg = leggauss(5) xi_q, w_q = 0.5 * (xg + 1.0), 0.5 * wg h = 1.0 / n ndof = 2 * (n + 1) Ke = (1.0 / h**3) * np.array([ [12, 6 * h, -12, 6 * h], [6 * h, 4 * h * h, -6 * h, 2 * h * h], [-12, -6 * h, 12, -6 * h], [6 * h, 2 * h * h, -6 * h, 4 * h * h]]) K = lil_matrix((ndof, ndof)) b = np.zeros(ndof) for e in range(n): xL = e * h d = [2 * e, 2 * e + 1, 2 * e + 2, 2 * e + 3] be = np.zeros(4) for q, w in zip(xi_q, w_q): be += w * f_load(xL + h * q) * _shape(q, h) be *= h for a in range(4): b[d[a]] += be[a] for c in range(4): K[d[a], d[c]] += Ke[a, c] free = np.arange(2, ndof - 2) # clamped at both ends U = np.zeros(ndof) U[free] = spsolve(csr_matrix(K)[free][:, free], b[free]) npts = DENSE_FACTOR * n + 1 xd = np.arange(npts) / float(DENSE_FACTOR * n) ud = np.empty(npts) for i, x in enumerate(xd): e = min(int(x / h), n - 1) ud[i] = _shape((x - e * h) / h, h) @ U[2 * e:2 * e + 4] return xd, ud def f_load(x): """The manufactured load, f = d4/dx4 [ e^x x^2 (1-x)^2 ].""" return np.exp(x) * (x**4 + 14 * x**3 + 49 * x**2 + 32 * x - 12) def u_exact(x): return np.exp(x) * (x * (1 - x))**2 def fail(msg): print(f"FAIL: {msg}") print("SCORE: 0.0") sys.exit(0) def req_float(value, what): """Coerce a submission-supplied value to a finite float, or fail cleanly. Every numeric field in results.json is agent-controlled, so a bare float() can raise on a string or None and crash the grader with no verdict at all. A malformed submission must score 0.0, not produce a traceback. """ # Reject non-numeric JSON types outright. bool is checked first because it # is a subclass of int in Python -- float(True) == 1.0 would otherwise let # "a": true satisfy a numeric field. Strings are rejected too: "1" is not a # JSON number, and the published schema asks for numbers. if isinstance(value, bool) or not isinstance( value, (int, float, np.integer, np.floating)): fail(f"{what} must be a JSON number; got {value!r}") try: fv = float(value) except (TypeError, ValueError, OverflowError): # A JSON integer may be arbitrarily long; float() raises OverflowError # on one too large to represent. fail(f"{what} must be a finite number; got {type(value).__name__} " f"out of floating-point range") if not np.isfinite(fv): fail(f"{what} must be finite; got {value!r}") return fv def npz_array(arr, key, what=None): """Fetch one NPZ member safely and require a real numeric dtype. np.load(..., allow_pickle=False) raises on an object-typed member, which would otherwise escape as an uncaught ValueError and exit with no verdict. String, byte-string, boolean and complex members are rejected too: they would be silently coerced by np.asarray(..., float) and pass as valid data. """ what = what or f"reference_arrays.npz['{key}']" try: a = arr[key] except Exception as e: fail(f"{what} could not be read ({type(e).__name__}: {e}); object-typed " f"or pickled members are not accepted") a = np.asanyarray(a) if a.dtype == object or a.dtype.kind in ("U", "S", "b", "V", "O"): fail(f"{what} has dtype {a.dtype!r}; a real numeric array is required " f"(string, byte, boolean and object arrays are rejected)") if not np.issubdtype(a.dtype, np.number): fail(f"{what} has non-numeric dtype {a.dtype!r}") if np.iscomplexobj(a): fail(f"{what} is complex; a real array is required") return a def req_int(value, what): """Coerce to an EXACT integer, or fail cleanly (int() would truncate).""" fv = req_float(value, what) if fv != np.floor(fv): fail(f"{what} must be an exact integer; got {value!r}") return int(fv) def note(msg): print(f"NOTE: {msg}") REVIEW = [] # machine-readable manual-review reasons def flag(reason): """Record a check that could not be settled automatically. Printed as a machine-readable REVIEW_REQUIRED line and, under --strict, score-affecting: an inconclusive check must not be silently certified. """ REVIEW.append(reason) print(f"REVIEW_REQUIRED: {reason}") # ============================================================================= # Subprocess runner: imports the submitted solver and regenerates the study. # ============================================================================= RUNNER = r''' import json, os, sys import numpy as np OUT_DIR, N_LIST_JSON, ARTIFACT_N = sys.argv[1], sys.argv[2], int(sys.argv[3]) CAPTURE_DIR = sys.argv[4] if len(sys.argv) > 4 else "." FIDELITY_N = set(json.loads(sys.argv[5])) if len(sys.argv) > 5 else set() POISON_SOLVES = json.loads(sys.argv[6]) if len(sys.argv) > 6 else False N_LIST = json.loads(N_LIST_JSON) DENSE_FACTOR = 10 GRID_ATOL = 1e-15 def u_exact(x): return np.exp(x) * (x * (1 - x))**2 STATS = {"n_solves": 0, "solve_dims": []} def _record_dim(A): try: STATS["n_solves"] += 1 STATS["solve_dims"].append(int(A.shape[-1])) except Exception: pass def _patch_solves(): """Observe global solves and optionally poison their numerical results.""" import numpy.linalg as nla import scipy.linalg as sla import scipy.sparse.linalg as ssl def poison(value): # splu/spilu return objects whose answer is produced later by .solve(). if hasattr(value, "solve") and callable(value.solve): class PoisonedFactor: def __init__(self, wrapped): self._wrapped = wrapped def solve(self, *a, **kw): return poison(self._wrapped.solve(*a, **kw)) def __getattr__(self, name): return getattr(self._wrapped, name) return PoisonedFactor(value) if callable(value): def poisoned_callable(*a, **kw): return poison(value(*a, **kw)) return poisoned_callable if isinstance(value, tuple): return tuple(poison(v) for v in value) if isinstance(value, list): return [poison(v) for v in value] if isinstance(value, np.ndarray) and value.dtype.kind in "fc": out = np.array(value, copy=True) out[...] = np.nan return out return value def wrap(mod, name, matrix_arg=0): orig = getattr(mod, name, None) if orig is None: return def wrapper(*a, **kw): if len(a) > matrix_arg: _record_dim(a[matrix_arg]) result = orig(*a, **kw) return poison(result) if POISON_SOLVES else result wrapper.__name__ = name setattr(mod, name, wrapper) for mod, names in ((ssl, ("spsolve", "spsolve_triangular", "splu", "spilu", "factorized", "cg", "gmres", "bicgstab", "minres", "lgmres")), (sla, ("solve", "lu_factor", "cho_factor", "cholesky", "solveh_banded")), (nla, ("solve",))): for name in names: wrap(mod, name) wrap(sla, "solve_banded", matrix_arg=1) _patch_solves() sys.path.insert(0, OUT_DIR) try: import fem_solver except Exception as e: print(json.dumps({"error": f"import failed: {type(e).__name__}: {e}"})) sys.exit(0) if not hasattr(fem_solver, "solve"): print(json.dumps({"error": "fem_solver.py does not define solve(n)"})) sys.exit(0) results = [] artifact = None for n in N_LIST: STATS["n_solves"] = 0 STATS["solve_dims"] = [] try: out = fem_solver.solve(n) except Exception as e: print(json.dumps({"error": f"solve({n}) raised {type(e).__name__}: {e}"})) sys.exit(0) try: xn, un, xd, ud = out except Exception: print(json.dumps({"error": f"solve({n}) must return a 4-tuple " f"(x_nodes, u_nodes, x_dense, u_dense)"})) sys.exit(0) xn = np.asarray(xn, float).ravel() un = np.asarray(un, float).ravel() xd = np.asarray(xd, float).ravel() ud = np.asarray(ud, float).ravel() nd = DENSE_FACTOR * n + 1 if xn.shape != (n + 1,) or un.shape != (n + 1,): print(json.dumps({"error": f"solve({n}) nodal arrays have shapes " f"{xn.shape}/{un.shape}, expected {(n + 1,)}"})) sys.exit(0) if xd.shape != (nd,) or ud.shape != (nd,): print(json.dumps({"error": f"solve({n}) dense arrays have shapes " f"{xd.shape}/{ud.shape}, expected {(nd,)}"})) sys.exit(0) nodes_ok = bool(np.max(np.abs(xn - np.arange(n + 1) / float(n))) <= GRID_ATOL) # u_nodes must be the dense field sampled at the nodes -- same DOFs, same # evaluation. Substituted nodal values (e.g. the exact solution) show up here. nodal_dense_dev = float(np.max(np.abs(un - ud[::DENSE_FACTOR]))) grid_ok = bool(np.max(np.abs(xd - np.arange(nd) / float(DENSE_FACTOR * n))) <= GRID_ATOL) finite = bool(np.all(np.isfinite(un)) and np.all(np.isfinite(ud))) row = { "n": int(n), "nodal_error": float(np.max(np.abs(un - u_exact(xn)))), "dense_error": float(np.max(np.abs(ud - u_exact(xd)))), "finite": finite, "nodes_ok": nodes_ok, "nodal_dense_dev": nodal_dense_dev, "grid_ok": grid_ok, "n_solves": STATS["n_solves"], "solve_dims": STATS["solve_dims"], } results.append(row) if int(n) == ARTIFACT_N: artifact = {"x": [float(v) for v in xd], "u": [float(v) for v in ud]} if int(n) in FIDELITY_N: np.save(os.path.join(CAPTURE_DIR, f"fid_{int(n)}.npy"), ud, allow_pickle=False) print(json.dumps({"results": results, "artifact": artifact})) ''' def run_solver(out_dir, n_list, budget, artifact_n=-1, poison_solves=False): """Execute the submitted solver in an isolated subprocess. Returns (payload_dict, error_string). Exactly one is None. """ with tempfile.TemporaryDirectory() as td: runner_path = os.path.join(td, "_repro_runner.py") with open(runner_path, "w") as fh: fh.write(RUNNER) env = dict(os.environ) env["PYTHONDONTWRITEBYTECODE"] = "1" try: proc = subprocess.run( [sys.executable, runner_path, os.path.abspath(out_dir), json.dumps([int(v) for v in n_list]), str(int(artifact_n)), td, json.dumps([int(v) for v in FIDELITY_N]), json.dumps(bool(poison_solves))], capture_output=True, text=True, timeout=budget, env=env, cwd=td, ) except subprocess.TimeoutExpired: return None, f"timed out after {budget:.0f}s" # Everything below must stay INSIDE the TemporaryDirectory: the captured # fidelity fields are files in `td`, which is removed on exit. if proc.returncode != 0: tail = (proc.stderr or "").strip().splitlines()[-3:] return None, f"runner exited {proc.returncode}: {' | '.join(tail)}" out = (proc.stdout or "").strip().splitlines() if not out: return None, "runner produced no output" try: payload = json.loads(out[-1]) except Exception as e: return None, f"unparseable runner output ({e})" if "error" in payload: return None, payload["error"] try: payload["fidelity_fields"] = { int(v): np.load(os.path.join(td, f"fid_{int(v)}.npy"), allow_pickle=False) for v in FIDELITY_N if os.path.exists(os.path.join(td, f"fid_{int(v)}.npy")) } except Exception as e: return None, f"could not load the fidelity capture ({e})" return payload, None def check_fidelity(payload, expected_n): """Compare the submitted dense fields against the grader's own solve. Returns the meshes actually verified. Called for both the mandatory and the large-mesh phase so that no mesh escapes the comparison. """ fid = payload.get("fidelity_fields") or {} want = [int(n) for n in expected_n if int(n) in FIDELITY_N] missing = [n for n in want if n not in fid] if missing: fail(f"the regeneration did not return the dense field at n = {missing}; " f"solve(n) must return the dense reconstruction on i/(10n)") for n in want: _, ud_ref = grader_fem(n) ud_sub = np.asarray(fid[n], float) if ud_sub.shape != ud_ref.shape: fail(f"solve({n}) dense field has shape {ud_sub.shape}, " f"expected {ud_ref.shape}") dev = float(np.max(np.abs(ud_sub - ud_ref))) if dev > FIDELITY_ATOL: fail(f"fem_solver.py does not implement the specified method: at " f"n = {n} its dense field differs from the grader's own " f"cubic-Hermite solution by {dev:.3e} (tolerance " f"{FIDELITY_ATOL:.0e}). Two legitimate implementations differing " f"only in quadrature rule agree to ~5e-12 here. A field using " f"the wrong element space or produced by genuine FEM on only " f"*some* meshes lands here.") if want: print(f"PASS scheme fidelity: the dense field matches the grader's own " f"cubic-Hermite solve at n = {', '.join(str(v) for v in want)} " f"to better than {FIDELITY_ATOL:.0e}") return set(want) def check_imports(solver_path): """Static allowlist over top-level imports. Returns a manual-review flag.""" try: with open(solver_path) as fh: src = fh.read() except Exception as e: fail(f"fem_solver.py unreadable ({e})") try: tree = ast.parse(src) except SyntaxError as e: fail(f"fem_solver.py does not parse ({e})") offenders, dynamic = set(), False for node in ast.walk(tree): if isinstance(node, ast.Import): for a in node.names: top = a.name.split(".")[0] if top not in ALLOWED_TOP_LEVEL: offenders.add(top) elif isinstance(node, ast.ImportFrom): if node.level == 0 and node.module: top = node.module.split(".")[0] if top not in ALLOWED_TOP_LEVEL: offenders.add(top) elif isinstance(node, ast.Name) and node.id == "__import__": dynamic = True elif isinstance(node, ast.Attribute) and node.attr in ( "import_module", "exec_module", "spec_from_file_location"): dynamic = True if offenders: fail("fem_solver.py imports disallowed module(s): " + ", ".join(sorted(offenders)) + " (Python + NumPy + SciPy only; no FEM/PDE library)") return dynamic def _rows(table, key, what): """Extract (n, h, err, beta) columns from a results.json study table.""" if not isinstance(table, list) or len(table) != len(N_VALUES): fail(f"results.json {what} table must be a list of {len(N_VALUES)} rows") try: # Validate the RAW 'n' before coercion -- int() truncates, so 32.25 would # otherwise pass as 32 (the same trap as n_values in the NPZ). raw_n = [r["n"] for r in table] for v in raw_n: fv = float(v) if not np.isfinite(fv) or fv != np.floor(fv): fail(f"results.json {what} table has a non-integer 'n' value " f"({v!r}); mesh sizes must be exact integers") srt = sorted(table, key=lambda r: int(r["n"])) n = np.array([req_int(r["n"], f"results.json {what} 'n'") for r in srt]) h = np.array([req_float(r["h"], f"results.json {what} 'h'") for r in srt]) e = np.array([req_float(r[key], f"results.json {what} '{key}'") for r in srt]) b = [r.get("beta", "__missing__") for r in srt] except Exception as exc: fail(f"results.json {what} table rows malformed ({exc}); each row needs " f"'n', 'h', '{key}', 'beta'") if any(v == "__missing__" for v in b): fail(f"results.json {what} table rows must each carry a 'beta' field") return n, h, e, b # ============================================================================= def main(out_dir, budget, do_exec, strict=False): solver_path = os.path.join(out_dir, "fem_solver.py") if not os.path.isfile(solver_path): fail("fem_solver.py missing (the solver is a required deliverable)") try: with open(os.path.join(out_dir, "results.json")) as fh: res = json.load(fh) except Exception as e: fail(f"results.json missing or invalid JSON ({e})") if not isinstance(res, dict): fail("results.json top level must be a JSON object") try: arr = np.load(os.path.join(out_dir, "reference_arrays.npz")) except Exception as e: fail(f"reference_arrays.npz missing or unreadable ({e})") artifact_provenance_bad = None # ---------------------------------------------------------------- schema study_keys = {"n_values": (7,), "nodal_errors": (7,), "dense_errors": (7,), "nodal_eoc": (7,), "dense_eoc": (7,)} artifact_keys = ("x_dense", "uh_dense", "uexact_dense", "error_dense") for key in study_keys: if key not in arr: fail(f"reference_arrays.npz missing required array '{key}'") # The Tier 3 artifact arrays are NOT mandatory: the rubric promises 0.5 for an # honest Tier 3 omission, so their absence must degrade to a Tier 3 failure # rather than a schema failure. missing_artifact = [k for k in artifact_keys if k not in arr] for key, shape in study_keys.items(): got = npz_array(arr, key).shape if got != shape: fail(f"'{key}' has shape {got}, expected {shape}") artifact_bad = "" if missing_artifact: artifact_bad = f"missing Tier 3 array(s): {', '.join(missing_artifact)}" else: for key in artifact_keys: got = npz_array(arr, key).shape if got != (ARTIFACT_POINTS,): artifact_bad = (f"'{key}' has shape {got}, expected " f"({ARTIFACT_POINTS},) (the n = {ARTIFACT_N} artifact)") break # Validate the RAW values before any integer coercion: np.asarray(x, int) # silently truncates, so 32.25 would otherwise pass as 32. n_raw = npz_array(arr, "n_values").ravel() if not np.all(np.isfinite(np.asarray(n_raw, float))): fail("n_values contains non-finite entries") if not np.all(np.asarray(n_raw, float) == np.floor(np.asarray(n_raw, float))): fail(f"n_values must be exact integers; got " f"{np.array2string(np.asarray(n_raw, float), precision=6)}") n_values = np.asarray(n_raw, int).ravel() nodal = np.asarray(npz_array(arr, "nodal_errors"), float).ravel() dense = np.asarray(npz_array(arr, "dense_errors"), float).ravel() nodal_eoc = np.asarray(npz_array(arr, "nodal_eoc"), float).ravel() dense_eoc = np.asarray(npz_array(arr, "dense_eoc"), float).ravel() if not np.array_equal(n_values, N_VALUES): fail(f"n_values must be 2^p for p = 5..11 (32..2048) in ascending order, " f"got {n_values.tolist()}") if not (np.all(np.isfinite(nodal)) and np.all(nodal > 0)): fail("nodal_errors must be finite and strictly positive") if not (np.all(np.isfinite(dense)) and np.all(dense > 0)): fail("dense_errors must be finite and strictly positive") # EOC arrays: entry 0 has no predecessor; entries 1..6 recomputed. nodal_rec = np.log2(nodal[:-1] / nodal[1:]) dense_rec = np.log2(dense[:-1] / dense[1:]) for name, sub, rec in (("nodal_eoc", nodal_eoc, nodal_rec), ("dense_eoc", dense_eoc, dense_rec)): if not np.isnan(sub[0]): fail(f"{name}[0] must be NaN (the coarsest mesh has no predecessor); " f"got {sub[0]!r}") if not np.allclose(sub[1:], rec, rtol=1e-6, atol=1e-12): fail(f"{name} inconsistent with the submitted errors: submitted " f"{np.array2string(sub[1:], precision=4)} vs recomputed " f"{np.array2string(rec, precision=4)}") if not artifact_bad: for key in artifact_keys: a = np.asarray(npz_array(arr, key), float) if not np.all(np.isfinite(a)): artifact_bad = f"'{key}' contains non-finite entries" break print(f"PASS schema: {9 - len(missing_artifact)}/9 arrays present, correctly " f"shaped, finite; mesh sequence correct; nodal_eoc and dense_eoc " f"consistent with the errors") # --------------------------------------------------- results.json sidecar try: nodal_study = res["nodal_study"]["table"] dense_study = res["dense_study"]["table"] except Exception as e: fail(f"results.json must contain nodal_study.table and dense_study.table ({e})") jn_n, jn_h, jn_e, jn_b = _rows(nodal_study, "e_n", "nodal_study") jd_n, jd_h, jd_e, jd_b = _rows(dense_study, "e_star_n", "dense_study") for tag, jn, jh in (("nodal_study", jn_n, jn_h), ("dense_study", jd_n, jd_h)): if not np.array_equal(jn, N_VALUES): fail(f"results.json {tag} n values do not match the mesh sequence") if not np.allclose(jh, 1.0 / np.asarray(N_VALUES, float), rtol=1e-12, atol=0.0): fail(f"results.json {tag} h values must equal 1/n exactly") if not np.allclose(jn_e, nodal, rtol=1e-6, atol=1e-18): fail("results.json nodal_study e_n values inconsistent with " "reference_arrays.npz nodal_errors") if not np.allclose(jd_e, dense, rtol=1e-6, atol=1e-18): fail("results.json dense_study e_star_n values inconsistent with " "reference_arrays.npz dense_errors") for tag, jb, rec in (("nodal_study", jn_b, nodal_rec), ("dense_study", jd_b, dense_rec)): if jb[0] is not None: fail(f"results.json {tag} beta for the coarsest mesh (n = 32) must be " f"null; got {jb[0]!r}") tail = np.array([req_float(v, f"results.json {tag} 'beta'") for v in jb[1:]]) if not np.allclose(tail, rec, rtol=1e-6, atol=1e-12): fail(f"results.json {tag} beta values inconsistent with the submitted " f"errors: submitted {np.array2string(tail, precision=4)} vs " f"recomputed {np.array2string(rec, precision=4)}") if not str(res.get("exact_solution", "")).strip(): fail("results.json must record the manufactured exact_solution") dofs = res.get("dofs_per_node") if (not isinstance(dofs, (list, tuple)) or len(dofs) != 2 or [str(v).strip().lower() for v in dofs] != ["value", "slope"]): fail(f"results.json dofs_per_node must be exactly ['value', 'slope'] in " f"that order -- the C^1 Hermite element carries a value DOF and a " f"slope DOF at each node; got {dofs!r}") # Descriptive metadata required by the published schema. Content is not # machine-verified (these are free-form descriptions), but they must be # present and non-empty rather than silently omitted. for key in ("problem", "exact_solution", "f", "notes"): if key not in res: fail(f"results.json is missing the required '{key}' field") if not isinstance(res[key], str) or not res[key].strip(): fail(f"results.json '{key}' must be a non-empty string; " f"got {res[key]!r}") if "asymptotic_order" not in res: fail("results.json is missing 'asymptotic_order'") ao = res["asymptotic_order"] if not np.isclose(req_float(ao, "results.json asymptotic_order"), 4.0, rtol=0.0, atol=0.05): fail(f"results.json asymptotic_order = {ao!r}; the cubic-Hermite method " f"is fourth order in the dense (field) L-infinity norm") # Claimed clean-regime order must be supported by the submitted errors. ls_clean = float(-np.polyfit(np.arange(CLEAN), np.log2(dense[:CLEAN]), 1)[0]) ooc = res.get("dense_study", {}).get("observed_order_clean_regime") if ooc is None: fail("results.json dense_study must report observed_order_clean_regime") oov = req_float(ooc, "results.json observed_order_clean_regime") if abs(oov - ls_clean) > 0.25: fail(f"results.json observed_order_clean_regime = {oov:.3f} is not " f"supported by the submitted dense errors (fitted order " f"{ls_clean:.3f} over n = 32..256)") art = res.get("artifact") if not isinstance(art, dict): fail("results.json missing the artifact section") for key in ("n", "points", "dense_max_error"): if key not in art: fail(f"results.json artifact section is missing '{key}'") if req_int(art["n"], "results.json artifact.n") != ARTIFACT_N: fail(f"results.json artifact.n must be {ARTIFACT_N}") if req_int(art["points"], "results.json artifact.points") != ARTIFACT_POINTS: fail(f"results.json artifact.points must be {ARTIFACT_POINTS}") if not np.isclose(req_float(art["dense_max_error"], "results.json artifact.dense_max_error"), dense[ARTIFACT_INDEX], rtol=1e-6, atol=0.0): fail(f"results.json artifact.dense_max_error = " f"{art.get('dense_max_error')!r} is inconsistent with the n = " f"{ARTIFACT_N} row of the dense study ({dense[ARTIFACT_INDEX]:.6e})") print("PASS sidecar: results.json fully consistent with the arrays (both " "study tables, h values, betas, orders, artifact metadata)") # === Tier 1: correct assembly + load + reconstruction (dense, clean) ===== for i in range(3): # n = 32, 64, 128 r = REF_DENSE_CLEAN[i] if not (r * (1 - REL_BAND) <= dense[i] <= r * (1 + REL_BAND)): fail(f"Tier 1: dense error at n={N_VALUES[i]} = {dense[i]:.3e} outside " f"+/-25% of {r:.3e} (slope DOFs omitted in reconstruction?)") if not all(ORDER_LO <= dense_rec[i] <= ORDER_HI for i in (0, 1)): # 64/32, 128/64 fail(f"Tier 1: coarse dense order {np.array2string(dense_rec[:2], precision=3)} " f"not ~4 (reconstruction not 4th order)") print(f"PASS Tier 1: dense errors in band; coarse order = " f"{np.array2string(dense_rec[:2], precision=3)}") # === Tier 2: two error measures, order 4, superconvergence, round-off ==== clean_eoc = dense_rec[:CLEAN - 1] # 3 orders across n=32..256 if sum(ORDER_LO <= b <= ORDER_HI for b in clean_eoc) < 2: fail(f"Tier 2: dense order in clean regime " f"{np.array2string(clean_eoc, precision=3)} not ~4") # nodal superconvergence: essentially exact at the nodes if not (nodal[0] < NODAL_SUPER_ABS and nodal[0] < dense[0] / NODAL_SUPER_RATIO): fail(f"Tier 2: nodal error at n=32 = {nodal[0]:.3e} not superconvergent " f"(should be << dense {dense[0]:.3e})") # round-off floor: refining past the clean regime stops helping if not (dense[-1] > ROUNDOFF_FACTOR * dense[CLEAN - 1]): fail(f"Tier 2: no round-off floor - dense error at n=2048 ({dense[-1]:.3e}) " f"should exceed that at n=256 ({dense[CLEAN-1]:.3e})") print(f"PASS Tier 2: dense order ~4 {np.array2string(clean_eoc, precision=3)}; " f"nodal superconvergent ({nodal[0]:.1e}); round-off floor present") # ======================= Provenance / solver enforcement ================= dynamic_imports = check_imports(solver_path) print("PASS dependencies: fem_solver.py imports only stdlib + NumPy + SciPy") if dynamic_imports: flag("fem_solver.py performs dynamic imports; the static dependency " "check is not exhaustive -> flag for manual review") if not do_exec: note("--no-exec given: solver execution and the reproducibility checks " "were SKIPPED; provenance is UNVERIFIED") else: payload, err = run_solver(out_dir, COARSE_N, budget, ARTIFACT_N) if err is not None: fail(f"could not regenerate results from fem_solver.py: {err}") rows = payload["results"] unobserved = [] def audit_global_solves(run_rows): """Check the published 2(n-1)-unknown assembled-system solve.""" for row in run_rows: n = int(row["n"]) if not row.get("n_solves", 0): unobserved.append(n) continue expected = 2 * (n - 1) dims = [int(v) for v in row.get("solve_dims", [])] if expected not in dims: fail(f"FEM assembly: solve({n}) made {row['n_solves']} observed " f"linear solve(s), with dimensions {dims}, but never " f"solved the required global system of {expected} free " "value/slope DOFs") for r in rows: if not r["finite"]: fail(f"regenerated solution at n={r['n']} is non-finite") if r.get("nodal_dense_dev", 0.0) > NODAL_DENSE_ATOL: fail(f"solve({r['n']}) returns nodal values that are not its own " f"dense field sampled at the nodes (max deviation " f"{r['nodal_dense_dev']:.3e} > {NODAL_DENSE_ATOL:.0e}). Both " f"come from the same degrees of freedom, so they must agree; " f"substituting the exact solution at the nodes to fake " f"superconvergence lands here.") if not r["nodes_ok"]: fail(f"solve({r['n']}) returned mesh nodes that are not x_j = j/n") if not r["grid_ok"]: fail(f"solve({r['n']}) returned a dense grid that is not " f"x_i = i/(10n)") audit_global_solves(rows) regen_nodal = np.array([r["nodal_error"] for r in rows]) regen_dense = np.array([r["dense_error"] for r in rows]) if not np.allclose(regen_nodal, nodal[:len(regen_nodal)], rtol=REPRO_RTOL, atol=0.0): fail("fem_solver.py does not reproduce the submitted nodal errors: " f"regenerated {np.array2string(regen_nodal, precision=6)} vs " f"submitted {np.array2string(nodal[:len(regen_nodal)], precision=6)}") if not np.allclose(regen_dense, dense[:len(regen_dense)], rtol=REPRO_RTOL, atol=0.0): fail("fem_solver.py does not reproduce the submitted dense errors: " f"regenerated {np.array2string(regen_dense, precision=6)} vs " f"submitted {np.array2string(dense[:len(regen_dense)], precision=6)}") print(f"PASS provenance: fem_solver.py regenerated both error measures at " f"n = {', '.join(str(v) for v in COARSE_N)} and matched the " f"submitted values to {REPRO_RTOL:.0e} relative") # ---- scheme fidelity: does the solver actually do finite elements? ---- # Reproducing its own numbers proves only that the arrays came from this # code. Here the returned field is compared against the grader's OWN # cubic-Hermite solve of the same problem. A solver that returns u_exact # plus a synthetic error -- which satisfies every other check, including # the dense band, the order, the superconvergence and the round-off # floor -- differs from a real FEM solution by the FEM error itself and # is rejected here by a factor of ~500. checked_fidelity = check_fidelity(payload, COARSE_N) # Field agreement alone is not a proof of assembly here: in the clean # regime the cubic-Hermite Ritz solution agrees to round-off with the # Hermite interpolant of the known exact solution. Poison every # intercepted solve result at n=32 and require that the returned field # change, become non-finite, or raise. A discarded dummy solve then # cannot certify an assembly-free analytic interpolant. A hand-written # solve is unobservable and is routed to manual review below. if rows[0].get("n_solves", 0): poison_n = int(COARSE_N[0]) poisoned, poison_err = run_solver( out_dir, [poison_n], budget, poison_solves=True) if poison_err is not None: if f"solve({poison_n}) raised" in poison_err: print("PASS causal assembly audit: poisoning intercepted " "solve results caused solve(n) to raise") else: flag("causal assembly audit could not complete after solve-" f"result poisoning ({poison_err})") else: poison_rows = poisoned.get("results") or [] poison_fields = poisoned.get("fidelity_fields") or {} normal_fields = payload.get("fidelity_fields") or {} if (len(poison_rows) != 1 or poison_n not in poison_fields or poison_n not in normal_fields): fail("causal assembly audit did not return the expected dense " f"field at n = {poison_n}") changed = (not poison_rows[0]["finite"] or not np.allclose(poison_fields[poison_n], normal_fields[poison_n], rtol=0.0, atol=1e-12, equal_nan=True)) if not changed: fail("FEM authenticity: poisoning every intercepted linear-" "solve result left the dense field unchanged. The " "observed solve calls do not produce the returned FEM " "solution (discarded/dummy solves are not assembly).") print("PASS causal assembly audit: poisoning intercepted solve " "results materially changed the dense field") # The n = 256 run is also the artifact run: compare the field pointwise. # Skipped when the artifact was honestly omitted -- there is nothing to # compare against, and that is a Tier 3 failure, not a provenance failure. # The runner has already verified that solve(256) returns the prescribed # grid, and Tier 3 verifies the saved x_dense against it, so only the # field itself is compared here. if artifact_bad: note(f"skipping the artifact provenance comparison ({artifact_bad})") else: af = payload.get("artifact") if not af: fail(f"the regeneration did not return the n = {ARTIFACT_N} artifact") regen_u = np.asarray(af["u"], float) sub_u = np.asarray(npz_array(arr, "uh_dense"), float) dev = float(np.max(np.abs(sub_u - regen_u))) if dev > ARTIFACT_FIELD_FAIL: fail(f"uh_dense was not produced by fem_solver.py: re-running " f"solve({ARTIFACT_N}) gives a field differing by {dev:.3e}, " f"which exceeds the artifact's entire accuracy budget " f"({ARTIFACT_FIELD_FAIL:.0e}) -- this is not round-off") if dev > ARTIFACT_FIELD_ATOL: # The grey band must NOT reach 1.0. Re-running the same code in # the same pinned environment reproduces the field bitwise (the # gold reference deviates by 0.0), so any deviation here is # unexplained. A warning that still awards full credit is exactly # the loophole that let a synthetic field -- one constructed to # match dense_errors[3] and satisfy every analytic check -- score # 1.0. Tier 3 now fails, which is also what EVALUATION.md claims. artifact_provenance_bad = ( f"the submitted n = {ARTIFACT_N} field differs from the one " f"regenerated by fem_solver.py by {dev:.3e}, above the " f"bit-reproducibility tolerance {ARTIFACT_FIELD_ATOL:.0e}. " f"Re-running the same solver in the same pinned environment " f"reproduces the field exactly, so this artifact is not the " f"output of the submitted code (spliced from another run, or " f"synthesised). If your grading environment genuinely uses a " f"different BLAS/LAPACK, adjudicate manually") else: print(f"PASS provenance: the n = {ARTIFACT_N} dense artifact is " f"reproduced pointwise by fem_solver.py (max deviation " f"{dev:.2e})") # Mandatory regeneration of the three largest meshes. These are NOT # expensive -- the reference solves n = 512, 1024 and 2048 in 0.11 s # combined -- and n = 2048 carries the Tier 2 round-off-floor check, so # leaving them best-effort let a solver stall only on those meshes while # fabricated fine-mesh values earned full credit. payload2, err2 = run_solver(out_dir, FINE_N, budget) if err2 is not None: fail(f"could not regenerate the large meshes (n = " f"{', '.join(str(v) for v in FINE_N)}) from fem_solver.py " f"({err2}). These rows carry the round-off-floor requirement, so " f"they cannot be accepted unverified; the reference solves all " f"three in well under a second. Raise --budget if the submission " f"is merely slow.") else: rows2 = payload2["results"] for r in rows2: if not r["finite"]: fail(f"regenerated solution at n={r['n']} is non-finite") if r.get("nodal_dense_dev", 0.0) > NODAL_DENSE_ATOL: fail(f"solve({r['n']}) returns nodal values that are not its " f"own dense field sampled at the nodes (max deviation " f"{r['nodal_dense_dev']:.3e} > {NODAL_DENSE_ATOL:.0e})") if not (r["nodes_ok"] and r["grid_ok"]): fail(f"solve({r['n']}) returned a nonconforming grid") audit_global_solves(rows2) regen2_nodal = np.array([r["nodal_error"] for r in rows2]) regen2_dense = np.array([r["dense_error"] for r in rows2]) if not np.allclose(regen2_nodal, nodal[CLEAN:], rtol=REPRO_RTOL, atol=0.0): fail("fem_solver.py does not reproduce the submitted nodal errors on " f"the large meshes: regenerated " f"{np.array2string(regen2_nodal, precision=6)} vs submitted " f"{np.array2string(nodal[CLEAN:], precision=6)}") if not np.allclose(regen2_dense, dense[CLEAN:], rtol=REPRO_RTOL, atol=0.0): fail("fem_solver.py does not reproduce the submitted dense errors on " f"the large meshes: regenerated " f"{np.array2string(regen2_dense, precision=6)} vs submitted " f"{np.array2string(dense[CLEAN:], precision=6)}") checked_fidelity |= check_fidelity(payload2, FINE_N) print("PASS provenance: all seven meshes regenerated from fem_solver.py") uncovered = sorted(set(FIDELITY_N) - checked_fidelity) if uncovered: fail(f"the dense field was never compared against the grader's own " f"solve at n = {uncovered}; every mesh must be verified, " f"otherwise a solver can do genuine FEM only where it is " f"checked and fabricate the rest") if unobserved: flag("no library global solve was observed at mesh size(s) " + ", ".join(str(n) for n in sorted(set(unobserved))) + ". A hand-written assembled-system solver is permitted, but " "the required FEM assembly and causal data path must be " "confirmed by reading fem_solver.py") # === Tier 3: dense reconstruction artifact at n = 256 ==================== # An omitted or malformed artifact is a Tier 3 failure (0.5), not a schema # failure (0.0) -- the rubric promises partial credit for honest omission. if artifact_bad: print(f"PARTIAL Tier 3 failed: {artifact_bad}") return emit_score(0.5, do_exec, strict) if artifact_provenance_bad: print(f"PARTIAL Tier 3 failed: {artifact_provenance_bad}") return emit_score(0.5, do_exec, strict) x = np.asarray(npz_array(arr, "x_dense"), float) uh = np.asarray(npz_array(arr, "uh_dense"), float) ue_sub = np.asarray(npz_array(arr, "uexact_dense"), float) esurf = np.asarray(npz_array(arr, "error_dense"), float) grid = np.arange(ARTIFACT_POINTS) / float(DENSE_FACTOR * ARTIFACT_N) ue = u_exact(grid) rec = float(np.max(np.abs(uh - ue))) ok, reason = True, "" if x[0] != 0.0 or x[-1] != 1.0: ok, reason = False, (f"x_dense must run from exactly 0.0 to exactly 1.0, " f"got [{x[0]!r}, {x[-1]!r}]") elif not np.all(np.diff(x) > 0): ok, reason = False, "x_dense is not strictly increasing" elif float(np.max(np.abs(x - grid))) > GRID_ATOL: ok, reason = False, ( f"x_dense is not the prescribed grid x_i = i/{DENSE_FACTOR * ARTIFACT_N} " f"(max deviation {float(np.max(np.abs(x - grid))):.3e} > " f"{GRID_ATOL:.0e})") elif not np.allclose(ue_sub, ue, rtol=1e-12, atol=1e-15): ok, reason = False, ( f"uexact_dense does not match the analytic exact solution " f"e^x x^2 (1-x)^2 on the artifact grid (max deviation " f"{float(np.max(np.abs(ue_sub - ue))):.3e})") elif rec >= ARTIFACT_MAX_ERR: ok, reason = False, f"artifact max error {rec:.3e} >= {ARTIFACT_MAX_ERR:.1e}" elif not np.allclose(esurf, uh - ue_sub, rtol=0.0, atol=1e-14): ok, reason = False, "error_dense != uh_dense - uexact_dense" elif not np.isclose(rec, dense[ARTIFACT_INDEX], rtol=1e-6, atol=0.0): ok, reason = False, ( f"dense_errors[{ARTIFACT_INDEX}] (n = {ARTIFACT_N}) = " f"{dense[ARTIFACT_INDEX]:.6e} but the maximum error recomputed from " f"uh_dense is {rec:.6e}: the convergence table and the Tier 3 artifact " f"did not come from the same run") if ok: print(f"PASS Tier 3: n={ARTIFACT_N} dense artifact on the exact " f"i/{DENSE_FACTOR * ARTIFACT_N} grid, exact solution and error array " f"consistent, and matching the n={ARTIFACT_N} study row; " f"max-norm {rec:.6e}") else: print(f"PARTIAL Tier 3 failed: {reason}") emit_score(1.0 if ok else 0.5, do_exec, strict) def emit_score(score, do_exec, strict=False): """Emit exactly one verdict line. Without solver execution nothing about provenance has been established, so no official score is produced -- a diagnostic run must not be mistaken for a graded one. Exit code 2 distinguishes 'ungraded' from a graded 0.0. """ if not do_exec: print("RESULT: UNGRADED - provenance was not verified (--no-exec). " "No official score is produced; re-run without --no-exec to grade.") sys.exit(2) if REVIEW: print(f"REVIEW: {len(REVIEW)} item(s) need human confirmation:") for r in REVIEW: print(f" - {r}") if strict: print("STRICT: an inconclusive check cannot be certified " "automatically -> score forced to 0.0") score = 0.0 else: print("REVIEW: none required - every check was conclusive") print(f"SCORE: {score}") if __name__ == "__main__": ap = argparse.ArgumentParser( description="Grade a cubic-Hermite FEM (Euler-Bernoulli beam) submission.") ap.add_argument("output_dir") ap.add_argument("--budget", type=float, default=DEFAULT_BUDGET, help="seconds allowed per solver-regeneration phase") ap.add_argument("--strict", action="store_true", help="treat any REVIEW_REQUIRED flag as a failure. Use when " "no human will read the manual-review items.") ap.add_argument("--no-exec", action="store_true", help="skip executing the submitted solver (diagnostics only; " "provenance is then unverified)") a = ap.parse_args() main(a.output_dir, a.budget, not a.no_exec, a.strict)