diff --git a/tools/build_mfem_bundle.py b/tools/build_mfem_bundle.py index 476a3764..24cd348d 100644 --- a/tools/build_mfem_bundle.py +++ b/tools/build_mfem_bundle.py @@ -87,6 +87,11 @@ def parse_args() -> argparse.Namespace: parser.add_argument("--cc-launcher", action="append", default=[]) parser.add_argument("--cxx-launcher", action="append", default=[]) parser.add_argument("--buildtype", required=True) + parser.add_argument("--optimization", choices=("plain", "0", "g", "1", "2", "3", "s"), required=True) + parser.add_argument("--debug", choices=("true", "false"), required=True) + parser.add_argument("--ndebug", choices=("true", "false"), required=True) + parser.add_argument("--compiler-syntax", choices=("gcc", "msvc"), required=True) + parser.add_argument("--compiler-id", default="") parser.add_argument("--jobs", type=int, required=True) parser.add_argument("--host-system", required=True) parser.add_argument("--cross-build", choices=("true", "false"), required=True) @@ -288,6 +293,32 @@ def compiler_flags(build_type: str) -> str: }[build_type] +def mfem_compile_flags(args: argparse.Namespace) -> str: + """Use Meson's effective options, not CMake's implicit -O3 -DNDEBUG. + + MFEM_DEBUG is separately selected by CMAKE_BUILD_TYPE. In particular, + debugoptimized retains symbols without enabling MFEM's expensive checks. + """ + if args.compiler_syntax == "msvc": + flags = { + "plain": [], "0": ["/Od"], "g": [], "1": ["/O1"], + "2": ["/O2"], "3": ["/O2", "/Gw"], "s": ["/O1", "/Gw"], + }[args.optimization].copy() + if args.debug == "true": + flags.append("/Zi") + if args.ndebug == "true": + flags.append("/DNDEBUG") + else: + flags = [] if args.optimization == "plain" else ["-O" + args.optimization] + if args.optimization == "s" and args.compiler_id in {"clang", "apple-clang", "emscripten"}: + flags = ["-Oz"] + if args.debug == "true": + flags.append("-g") + if args.ndebug == "true": + flags.append("-DNDEBUG") + return " ".join(flags) + + def joined_flags(values: Sequence[str]) -> str: return " ".join(shlex.quote(value) for value in values) @@ -1411,6 +1442,8 @@ def build_mfem( options = common_cmake_args( args, build_type, c_compiler=args.cc, cxx_compiler=args.cxx ) + [ + f"-DCMAKE_C_FLAGS_{build_type.upper()}={mfem_compile_flags(args)}", + f"-DCMAKE_CXX_FLAGS_{build_type.upper()}={mfem_compile_flags(args)}", "-DCMAKE_FIND_PACKAGE_PREFER_CONFIG=OFF", "-DMFEM_ENABLE_TESTING=OFF", "-DMFEM_ENABLE_EXAMPLES=OFF", @@ -1606,6 +1639,14 @@ def validate_mfem_feature_header( if not header.is_file(): raise RuntimeError("MFEM feature validation failed: _config.hpp is missing") text = header.read_text(encoding="utf-8") + debug_defined = re.search( + r"^\s*#\s*define\s+MFEM_DEBUG(?:\s|$)", text, flags=re.MULTILINE + ) is not None + if debug_defined != (cmake_build_type(args.buildtype) == "Debug"): + raise RuntimeError( + "MFEM feature validation failed: MFEM_DEBUG disagrees with " + f"Meson buildtype={args.buildtype}; rebuild the MFEM bundle" + ) for feature_name, macro in CMAKE_FEATURES.items(): is_defined = re.search( rf"^\s*#\s*define\s+{re.escape(macro)}(?:\s|$)", @@ -2432,6 +2473,7 @@ def main() -> int: if config_file.is_file(): old = json.loads(config_file.read_text(encoding="utf-8")) if old.get("fingerprint") == fingerprint and existing_libraries: + validate_mfem_feature_header(args, features) if args.stamp: write_text_atomic(args.stamp, fingerprint + "\n") print("[mfem-bundle] configuration unchanged; reusing private prefix", flush=True) @@ -2586,6 +2628,8 @@ def main() -> int: "features": features, "compiler": compiler_identity, "build_type": build_type, + "mfem_compile_flags": mfem_compile_flags(args), + "mfem_debug": build_type == "Debug", "prefix": os.fspath(args.prefix), } write_text_atomic(config_file, json.dumps(manifest, indent=2, sort_keys=True) + "\n") diff --git a/tools/read_serif_archive.py b/tools/read_serif_archive.py new file mode 100644 index 00000000..c32ff5a0 --- /dev/null +++ b/tools/read_serif_archive.py @@ -0,0 +1,171 @@ +#!/usr/bin/env python3 +"""Inspect SERiF archive data without compiling SERiF or installing dependencies. + +This is a data reader, not a physics/model factory. The C++ loader additionally +checks SERiF version tolerance, component schemas and discretization compatibility. +""" + +import argparse +import csv +import json +import shlex +import struct +import sys + + +def fingerprint(data): + """Match the archive's versioned corruption-detection fingerprint.""" + value = 14695981039346656037 + for byte in data: + value = ((value ^ byte) * 1099511628211) & ((1 << 64) - 1) + return f"fnv1a64-v1:{value:016x}" + + +def unquote(data): + """Decode std::quoted text, whose escapes differ from JSON escapes.""" + if len(data) < 2 or data[0] != 34 or data[-1] != 34: + raise ValueError("Invalid quoted string") + result = bytearray() + index = 1 + while index < len(data) - 1: + if data[index] == 92: + index += 1 + if index >= len(data) - 1: + raise ValueError("Unfinished string escape") + elif data[index] == 34: + raise ValueError("Unescaped quote inside string") + result.append(data[index]) + index += 1 + return bytes(result) + + +class Archive: + """Read the ledger once and seek directly to individual data sections.""" + + def __init__(self, path, maximum_bytes=1024 * 1024 * 1024): + self.file = open(path, "rb") + self.maximum_bytes = maximum_bytes + try: + self.file.seek(0, 2) + size = self.file.tell() + self.file.seek(0) + header = self.file.read(512).decode("ascii") + lines = header.split("END_HEADER\n", 1)[0].splitlines() + if not lines or lines[0] != "SERIF_ARCHIVE 1": + raise ValueError("Unsupported SERiF archive format") + self.header = dict(shlex.split(line) for line in lines[1:]) + offset = int(self.header["ledger_offset"]) + count = int(self.header["ledger_bytes"]) + if (offset != 512 or count < 0 or count > min(maximum_bytes, 64 * 1024 * 1024) + or offset + count > size or int(self.header["file_bytes"]) != size + or int(self.header["data_offset"]) != offset + count): + raise ValueError("Invalid archive section bounds") + ledger = self.file.read(count) + if fingerprint(ledger) != self.header["ledger_checksum"]: + raise ValueError("Damaged archive ledger") + self.entries = {} + previous_end = offset + count + previous_path = "" + for line in ledger.decode("utf-8").splitlines()[1:]: + path, schema, version, kind, extent, start, length, checksum = shlex.split(line) + entry = dict(schema=schema, version=int(version), kind=int(kind), count=int(extent), + offset=int(start), bytes=int(length), checksum=checksum) + if path <= previous_path or not 0 <= entry["kind"] <= 8 or entry["count"] < 0: + raise ValueError("Invalid or duplicate ledger entry") + if entry["kind"]: + if (entry["offset"] < previous_end or entry["bytes"] < 0 + or entry["bytes"] > maximum_bytes or entry["offset"] + entry["bytes"] > size): + raise ValueError("Invalid or overlapping payload bounds") + previous_end = entry["offset"] + entry["bytes"] + elif entry["offset"] or entry["bytes"] or entry["count"]: + raise ValueError("Directory contains a payload") + self.entries[path] = entry + previous_path = path + if len(self.entries) != int(self.header["entries"]): + raise ValueError("Incorrect ledger extent") + if self.header["format"] not in ("binary", "text"): + raise ValueError("Unknown numeric encoding") + except Exception: + self.file.close() + raise + + def __enter__(self): + """Allow a with-statement to own the open file.""" + return self + + def __exit__(self, *_): + """Close the reader after analysis, including when decoding fails.""" + self.file.close() + + def read(self, path): + """Return one string, number or numerical array without reading sibling fields.""" + entry = self.entries[path] + kind, count = entry["kind"], entry["count"] + self.file.seek(entry["offset"]) + data = self.file.read(entry["bytes"]) + if len(data) != entry["bytes"] or fingerprint(data) != entry["checksum"]: + raise ValueError(f"Damaged payload: {path}") + binary = self.header["format"] == "binary" + if kind == 0: + return None + if kind == 1: + allowed = (b"\x00", b"\x01") if binary else (b"0", b"1") + if count != 1 or data not in allowed: + raise ValueError("Invalid boolean") + return data == allowed[1] + if kind == 5: + text = data if binary else unquote(data) + if len(text) != count: + raise ValueError("Invalid string length") + return text.decode("utf-8", errors="surrogateescape") + if kind in (2, 3, 4) and count != 1: + raise ValueError("Invalid scalar extent") + if count > len(data) or count > self.maximum_bytes // 8: + raise ValueError("Invalid numerical extent") + code = {2: "q", 3: "Q", 4: "d", 6: "d", 7: "q", 8: "Q"}[kind] + if binary: + if len(data) != 8 * count: + raise ValueError("Invalid binary array length") + values = list(struct.unpack(f"<{count}{code}", data)) + else: + tokens = data.decode("ascii").split() + if len(tokens) != count: + raise ValueError("Invalid text array length") + values = [] + for token in tokens: + if code == "d": + value = struct.unpack(" int: - parsed = int(value) - if parsed < 1: - raise argparse.ArgumentTypeError("must be at least one") - return parsed - - -def integer_list(value: str) -> list[int]: - try: - parsed = [int(item.strip()) for item in value.split(",") if item.strip()] - except ValueError as error: - raise argparse.ArgumentTypeError("expected comma-separated integers") from error - if not parsed or any(item < 1 for item in parsed): - raise argparse.ArgumentTypeError("all entries must be positive") - return parsed - - -def command_output(command: list[str]) -> str | None: - try: - result = subprocess.run( - command, - check=False, - capture_output=True, - text=True, - timeout=10, - ) - except (OSError, subprocess.TimeoutExpired): - return None - if result.returncode != 0: - return None - output = result.stdout.strip() - return output if output else None - - -def resolve_program(value: str) -> str | None: - candidate = pathlib.Path(value) - if candidate.is_file(): - return str(candidate.resolve()) - discovered = shutil.which(value) - return str(pathlib.Path(discovered).resolve()) if discovered else None - - -def compact_cpu_list(cpus: list[int]) -> str: - if not cpus: - return "" - groups: list[str] = [] - start = previous = cpus[0] - for current in cpus[1:]: - if current == previous + 1: - previous = current - continue - groups.append(str(start) if start == previous else f"{start}-{previous}") - start = previous = current - groups.append(str(start) if start == previous else f"{start}-{previous}") - return ",".join(groups) - - -def affinity_prefix(case: Case, setting: str) -> tuple[list[str], str | None]: - if setting == "none" or shutil.which("taskset") is None: - return [], None - if setting != "auto": - return ["taskset", "-c", setting], setting - if not hasattr(os, "sched_getaffinity"): - return [], None - available = sorted(os.sched_getaffinity(0)) - count = min(len(available), max(case.ranks, case.threads)) - chosen = compact_cpu_list(available[:count]) - return (["taskset", "-c", chosen] if chosen else []), chosen or None - - -def expand_cases(args: argparse.Namespace) -> list[Case]: - requested = [item.strip() for item in args.cases.split(",") if item.strip()] - valid = {"cpu", "ceed-cpu", "omp", "mpi", "cuda", "ceed-cuda"} - unknown = sorted(set(requested) - valid) - if unknown: - raise ValueError(f"unknown benchmark case(s): {', '.join(unknown)}") - cases: list[Case] = [] - for requested_case in requested: - if requested_case == "cpu": - cases.append(Case("cpu-1r-1t", "cpu", 1, 1)) - elif requested_case == "ceed-cpu": - cases.append(Case("ceed-cpu-1r-1t", args.ceed_cpu_device, 1, 1)) - elif requested_case == "omp": - cases.extend( - Case(f"omp-1r-{threads}t", "omp", 1, threads) - for threads in args.omp_threads - ) - elif requested_case == "mpi": - cases.extend( - Case(f"mpi-{ranks}r-1t", "cpu", ranks, 1) - for ranks in args.mpi_ranks - ) - elif requested_case == "cuda": - cases.append(Case("cuda-1r-1t", "cuda", 1, 1)) - elif requested_case == "ceed-cuda": - cases.append(Case("ceed-cuda-1r-1t", args.ceed_cuda_device, 1, 1)) - deduplicated: list[Case] = [] - for case in cases: - if case not in deduplicated: - deduplicated.append(case) - return deduplicated - - -def run_case( - args: argparse.Namespace, case: Case, trial: int, warmup_trial: bool -) -> dict[str, Any]: - prefix, affinity = affinity_prefix(case, args.affinity) - command = prefix + [args.launcher, "-n", str(case.ranks)] - command.extend(args.launcher_arg) - command.extend( - [ - args.executable, - "--device", - case.device, - "--mesh-n", - str(args.mesh_n), - "--order", - str(args.order), - "--applications", - str(args.applications), - "--max-applications", - str(args.max_applications), - "--minimum-apply-seconds", - str(args.minimum_apply_seconds), - "--warmup-applications", - str(args.warmup_applications), - "--relative-tolerance", - str(args.relative_tolerance), - "--max-iterations", - str(args.max_iterations), - "--trial", - str(trial), - ] - ) - if args.solve: - command.append("--solve") - environment = os.environ.copy() - environment.update( - { - "OMP_NUM_THREADS": str(case.threads), - "OMP_DYNAMIC": "FALSE", - "OMP_PROC_BIND": "close", - "OMP_PLACES": "cores", - "FI_PROVIDER": environment.get("FI_PROVIDER", "sockets"), - } - ) - environment.setdefault("OMPI_ALLOW_RUN_AS_ROOT", "1") - environment.setdefault("OMPI_ALLOW_RUN_AS_ROOT_CONFIRM", "1") - selected_interface = args.network_interface - if selected_interface == "local": - selected_interface = "lo0" if sys.platform == "darwin" else "lo" - if selected_interface == "auto": - environment.pop("FI_SOCKETS_IFACE", None) - environment.pop("HYDRA_IFACE", None) - environment.pop("MPICH_INTERFACE_HOSTNAME", None) - else: - environment["FI_SOCKETS_IFACE"] = selected_interface - environment["HYDRA_IFACE"] = selected_interface - environment.pop("MPICH_INTERFACE_HOSTNAME", None) - if sys.platform == "darwin" and args.network_interface == "local": - environment["MPICH_INTERFACE_HOSTNAME"] = "127.0.0.1" - if case.device.startswith("cuda") or case.device.startswith("ceed-cuda"): - if args.gpu_aware_mpi: - environment["MFEM_GPU_AWARE_MPI"] = "1" - else: - environment.pop("MFEM_GPU_AWARE_MPI", None) - if args.cuda_visible_devices is not None: - environment["CUDA_VISIBLE_DEVICES"] = args.cuda_visible_devices - - started = time.perf_counter() - completed = subprocess.run( - command, - check=False, - capture_output=True, - text=True, - env=environment, - timeout=args.timeout, - ) - wall_seconds = time.perf_counter() - started - combined_output = "\n".join((completed.stdout, completed.stderr)) - records = [ - json.loads(line[len(JSON_PREFIX) :]) - for line in combined_output.splitlines() - if line.startswith(JSON_PREFIX) - ] - if completed.returncode != 0 or len(records) != 1: - tail = "\n".join(combined_output.splitlines()[-30:]) - raise RuntimeError( - f"{case.name} trial {trial} failed with exit code " - f"{completed.returncode}; JSON records={len(records)}\n{tail}" - ) - record = records[0] - record.update( - { - "case": case.name, - "requested_device": case.device, - "runner_trial": trial, - "discarded_warmup_trial": warmup_trial, - "wall_seconds": wall_seconds, - "cpu_affinity": affinity, - "network_interface": environment.get("HYDRA_IFACE"), - "cuda_visible_devices": environment.get("CUDA_VISIBLE_DEVICES"), - "gpu_aware_mpi": environment.get("MFEM_GPU_AWARE_MPI") == "1", - "command": command, - } - ) - if not record.get("valid", False): - raise RuntimeError(f"{case.name} trial {trial} failed numerical validation") - return record - - -def median_absolute_deviation(values: list[float]) -> float: - center = statistics.median(values) - return statistics.median(abs(value - center) for value in values) - - -def validate_cross_backend( - records: list[dict[str, Any]], relative_tolerance: float -) -> dict[str, float]: - measured = [row for row in records if not row["discarded_warmup_trial"]] - baseline = [row for row in measured if row["case"] == "cpu-1r-1t"] - if not baseline: - return {} - validation_fields = ["probe_norm"] - if measured and all(row.get("solve_ran", False) for row in measured): - validation_fields.append("solution_norm") - references = { - field: statistics.median(float(row[field]) for row in baseline) - for field in validation_fields - } - maximum_errors = {field: 0.0 for field in references} - for row in measured: - for field, reference in references.items(): - relative_error = abs(float(row[field]) - reference) / max( - abs(reference), 1.0e-300 - ) - maximum_errors[field] = max(maximum_errors[field], relative_error) - if relative_error > relative_tolerance: - raise RuntimeError( - f"{row['case']} disagrees with CPU for {field}: relative " - f"error {relative_error:.3e} > {relative_tolerance:.3e}" - ) - return maximum_errors - - -def checkpoint_results( - output_dir: pathlib.Path, - metadata: dict[str, Any], - records: list[dict[str, Any]], -) -> None: - output_dir.mkdir(parents=True, exist_ok=True) - destination = output_dir / "partial-results.json" - temporary = output_dir / "partial-results.json.tmp" - temporary.write_text( - json.dumps({"schema_version": 1, "metadata": metadata, "raw": records}, indent=2) - + "\n", - encoding="utf-8", - ) - temporary.replace(destination) - - -def summarize(cases: list[Case], records: list[dict[str, Any]]) -> list[dict[str, Any]]: - summaries: list[dict[str, Any]] = [] - for case in cases: - rows = [ - row - for row in records - if row["case"] == case.name and not row["discarded_warmup_trial"] - ] - if not rows: - continue - throughput = [float(row["apply_mdof_per_second"]) for row in rows] - apply_seconds = [float(row["apply_seconds"]) for row in rows] - solve_seconds = [float(row["solve_seconds"]) for row in rows] - total_seconds = [float(row["total_seconds"]) for row in rows] - timed_applications = [int(row["applications"]) for row in rows] - summaries.append( - { - "case": case.name, - "device": case.device, - "ranks": case.ranks, - "threads": case.threads, - "trials": len(rows), - "global_true_dofs": rows[0]["global_true_dofs"], - "requested_applications": rows[0]["requested_applications"], - "applications_min": min(timed_applications), - "applications_max": max(timed_applications), - "apply_mdof_per_second_median": statistics.median(throughput), - "apply_mdof_per_second_mad": median_absolute_deviation(throughput), - "apply_mdof_per_second_min": min(throughput), - "apply_mdof_per_second_max": max(throughput), - "apply_seconds_median": statistics.median(apply_seconds), - "solve_seconds_median": statistics.median(solve_seconds), - "total_seconds_median": statistics.median(total_seconds), - "cg_iterations": sorted({row["cg_iterations"] for row in rows}), - "solve_ran": all(bool(row["solve_ran"]) for row in rows), - "verified_relative_residual_max": max( - float(row["verified_relative_residual"]) for row in rows - ), - } - ) - baseline = next((row for row in summaries if row["case"] == "cpu-1r-1t"), None) - baseline_rate = ( - float(baseline["apply_mdof_per_second_median"]) if baseline else None - ) - for row in summaries: - row["speedup_vs_cpu"] = ( - float(row["apply_mdof_per_second_median"]) / baseline_rate - if baseline_rate - else None - ) - return summaries - - -def markdown_table(summaries: list[dict[str, Any]]) -> str: - lines = [ - "| Case | Device | Ranks | Threads | Timed apps | Apply window (s) | Apply MDoF/s, median [range] | MAD | Speedup | Solve (s) | CG iters |", - "|---|---|---:|---:|---:|---:|---:|---:|---:|---:|---:|", - ] - for row in summaries: - speedup = row["speedup_vs_cpu"] - speedup_text = f"{speedup:.2f}x" if speedup is not None else "n/a" - solve_text = f"{row['solve_seconds_median']:.4f}" if row["solve_ran"] else "—" - iteration_text = ( - ",".join(str(value) for value in row["cg_iterations"]) - if row["solve_ran"] - else "—" - ) - applications_text = ( - str(row["applications_min"]) - if row["applications_min"] == row["applications_max"] - else f"{row['applications_min']}-{row['applications_max']}" - ) - lines.append( - f"| {row['case']} | `{row['device']}` | {row['ranks']} | " - f"{row['threads']} | {applications_text} | " - f"{row['apply_seconds_median']:.3f} | " - f"{row['apply_mdof_per_second_median']:.2f} " - f"[{row['apply_mdof_per_second_min']:.2f}, " - f"{row['apply_mdof_per_second_max']:.2f}] | " - f"{row['apply_mdof_per_second_mad']:.2f} | " - f"{speedup_text} | {solve_text} | {iteration_text} |" - ) - return "\n".join(lines) + "\n" - - -def collect_metadata(args: argparse.Namespace) -> dict[str, Any]: - executable = pathlib.Path(args.executable).resolve() - source_root = pathlib.Path(__file__).resolve().parents[1] - return { - "timestamp_utc": dt.datetime.now(dt.timezone.utc).isoformat(), - "hostname": platform.node(), - "platform": platform.platform(), - "python": sys.version, - "cpu_count": os.cpu_count(), - "lscpu": command_output(["lscpu"]) if shutil.which("lscpu") else None, - "executable": str(executable), - "launcher": str(pathlib.Path(args.launcher).resolve()), - "launcher_version": command_output([args.launcher, "--version"]), - "compiler": command_output(["c++", "--version"]), - "nvidia_smi": command_output( - [ - "nvidia-smi", - "--query-gpu=name,driver_version,temperature.gpu,pstate,clocks.sm,power.draw", - "--format=csv,noheader", - ] - ) - if shutil.which("nvidia-smi") - else None, - "git_commit": command_output( - ["git", "-C", str(source_root), "rev-parse", "HEAD"] - ), - "git_dirty": bool( - command_output(["git", "-C", str(source_root), "status", "--porcelain"]) - ), - "parameters": { - "mesh_n": args.mesh_n, - "order": args.order, - "applications": args.applications, - "max_applications": args.max_applications, - "minimum_apply_seconds": args.minimum_apply_seconds, - "warmup_applications": args.warmup_applications, - "trials": args.trials, - "discarded_warmup_trials": args.warmup_trials, - "relative_tolerance": args.relative_tolerance, - "max_iterations": args.max_iterations, - "solve": args.solve, - "affinity": args.affinity, - "network_interface": args.network_interface, - "cuda_visible_devices": args.cuda_visible_devices, - "gpu_aware_mpi": args.gpu_aware_mpi, - "cross_backend_relative_tolerance": args.cross_backend_relative_tolerance, - }, - } - - -def write_results( - output_dir: pathlib.Path, - metadata: dict[str, Any], - records: list[dict[str, Any]], - summaries: list[dict[str, Any]], -) -> None: - output_dir.mkdir(parents=True, exist_ok=True) - (output_dir / "results.json").write_text( - json.dumps( - {"schema_version": 1, "metadata": metadata, "raw": records, "summary": summaries}, - indent=2, - ) - + "\n", - encoding="utf-8", - ) - fieldnames = [ - "case", - "device", - "ranks", - "threads", - "trials", - "global_true_dofs", - "requested_applications", - "applications_min", - "applications_max", - "apply_mdof_per_second_median", - "apply_mdof_per_second_mad", - "apply_mdof_per_second_min", - "apply_mdof_per_second_max", - "speedup_vs_cpu", - "apply_seconds_median", - "solve_seconds_median", - "solve_ran", - "total_seconds_median", - "cg_iterations", - "verified_relative_residual_max", - ] - with (output_dir / "results.csv").open("w", newline="", encoding="utf-8") as handle: - writer = csv.DictWriter(handle, fieldnames=fieldnames) - writer.writeheader() - for summary in summaries: - writer.writerow(summary) - first = records[0] - if first["minimum_apply_seconds"] > 0.0: - timing_description = ( - "Application counts were selected with a " - f"{first['minimum_apply_seconds']:.3g}-second calibration window, then " - "measured in a fresh timing window." - ) - else: - timing_description = ( - "The requested application count was measured directly in one timing " - "window." - ) - report = "# MFEM backend benchmark\n\n" - report += ( - f"Host: `{metadata['hostname']}` \n" - f"UTC timestamp: `{metadata['timestamp_utc']}` \n" - f"Workload: {first['global_elements']:,} 3D hex elements, H1 order " - f"{first['order']}, {first['global_true_dofs']:,} global true DoFs, " - f"at least {first['requested_applications']} timed operator applications " - f"per trial. {timing_description}\n\n" - ) - report += markdown_table(summaries) - report += ( - "\nThe table reports medians from fresh processes. Operator throughput uses " - "the same global partial-assembly diffusion operator and true-DOF field in " - "every case. MPI timings use the slowest rank.\n" - ) - (output_dir / "results.md").write_text(report, encoding="utf-8") - - -def parse_arguments() -> argparse.Namespace: - parser = argparse.ArgumentParser(description=__doc__) - parser.add_argument("--executable", required=True) - parser.add_argument("--launcher", required=True) - parser.add_argument( - "--cases", - default="cpu,ceed-cpu,omp,mpi,cuda,ceed-cuda", - help="Comma-separated subset of cpu, ceed-cpu, omp, mpi, cuda, ceed-cuda.", - ) - parser.add_argument("--omp-threads", type=integer_list, default=[8]) - parser.add_argument("--mpi-ranks", type=integer_list, default=[4]) - parser.add_argument("--trials", type=positive_int, default=5) - parser.add_argument("--warmup-trials", type=int, default=0) - parser.add_argument("--mesh-n", type=positive_int, default=48) - parser.add_argument("--order", type=positive_int, default=3) - parser.add_argument("--applications", type=positive_int, default=50) - parser.add_argument("--max-applications", type=positive_int, default=1_000_000) - parser.add_argument( - "--minimum-apply-seconds", - type=float, - default=1.0, - help=( - "Double the requested applications until a calibration window is this " - "long, then measure a fresh window." - ), - ) - parser.add_argument("--warmup-applications", type=int, default=5) - parser.add_argument("--relative-tolerance", type=float, default=1.0e-6) - parser.add_argument("--max-iterations", type=positive_int, default=1000) - parser.add_argument( - "--solve", - action="store_true", - help="Also run and time the Jacobi-preconditioned CG validation solve.", - ) - parser.add_argument("--timeout", type=positive_int, default=900) - parser.add_argument( - "--affinity", - default="auto", - help="CPU list for taskset, 'auto' for a compact set, or 'none'.", - ) - parser.add_argument( - "--network-interface", - default="local", - help="'local' for platform loopback, 'auto' for routing, or an interface name.", - ) - parser.add_argument( - "--cuda-visible-devices", - default=None, - help="Optional CUDA_VISIBLE_DEVICES value for GPU cases.", - ) - parser.add_argument( - "--gpu-aware-mpi", - action="store_true", - help="Opt in only when the selected MPI implementation is CUDA-aware.", - ) - parser.add_argument( - "--cross-backend-relative-tolerance", - type=float, - default=1.0e-5, - help="Maximum relative difference from CPU for probe and solution norms.", - ) - parser.add_argument( - "--ceed-cpu-device", - default="ceed-cpu", - help="Full libCEED CPU resource or portable ceed-cpu alias.", - ) - parser.add_argument( - "--ceed-cuda-device", - default="ceed-cuda", - help="Full libCEED CUDA resource or portable ceed-cuda alias.", - ) - parser.add_argument( - "--launcher-arg", - action="append", - default=[], - help="Additional MPI launcher argument; repeat as needed.", - ) - parser.add_argument( - "--output-dir", - type=pathlib.Path, - default=pathlib.Path("benchmark-results"), - ) - args = parser.parse_args() - if args.warmup_trials < 0 or args.warmup_applications < 0: - parser.error("warm-up counts may not be negative") - if args.relative_tolerance <= 0.0 or args.cross_backend_relative_tolerance <= 0.0: - parser.error("tolerances must be positive") - if args.minimum_apply_seconds < 0.0: - parser.error("minimum apply seconds may not be negative") - if args.max_applications < args.applications: - parser.error("max applications must be at least the requested applications") - for attribute in ("executable", "launcher"): - value = getattr(args, attribute) - resolved = resolve_program(value) - if resolved is None: - parser.error(f"program not found: {value}") - setattr(args, attribute, resolved) - return args - - -def main() -> int: - args = parse_arguments() - try: - cases = expand_cases(args) - except ValueError as error: - print(f"error: {error}", file=sys.stderr) - return 2 - metadata = collect_metadata(args) - records: list[dict[str, Any]] = [] - - for warmup in range(args.warmup_trials): - for case in cases: - print(f"warm-up {warmup + 1}/{args.warmup_trials}: {case.name}", flush=True) - records.append(run_case(args, case, -(warmup + 1), True)) - checkpoint_results(args.output_dir, metadata, records) - - schedule = [(trial, case) for trial in range(1, args.trials + 1) for case in cases] - random.Random(20260905).shuffle(schedule) - for index, (trial, case) in enumerate(schedule, start=1): - print( - f"run {index}/{len(schedule)}: {case.name}, trial {trial}/{args.trials}", - flush=True, - ) - records.append(run_case(args, case, trial, False)) - checkpoint_results(args.output_dir, metadata, records) - - summaries = summarize(cases, records) - workload_signatures = { - ( - row["global_elements"], - row["global_true_dofs"], - row["mesh_n"], - row["order"], - row["requested_applications"], - row["minimum_apply_seconds"], - ) - for row in records - if not row["discarded_warmup_trial"] - } - if len(workload_signatures) != 1: - raise RuntimeError( - f"cases did not use one fixed global workload: {workload_signatures}" - ) - metadata["cross_backend_maximum_relative_errors"] = validate_cross_backend( - records, args.cross_backend_relative_tolerance - ) - if shutil.which("nvidia-smi"): - metadata["nvidia_smi_after"] = command_output( - [ - "nvidia-smi", - "--query-gpu=name,driver_version,temperature.gpu,pstate,clocks.sm,power.draw", - "--format=csv,noheader", - ] - ) - write_results(args.output_dir, metadata, records, summaries) - table = markdown_table(summaries) - print("\n" + table) - print(f"JSON, CSV, and Markdown written to {args.output_dir.resolve()}") - return 0 - - -if __name__ == "__main__": - try: - raise SystemExit(main()) - except (OSError, subprocess.TimeoutExpired, RuntimeError) as error: - print(f"benchmark failed: {error}", file=sys.stderr) - raise SystemExit(1) diff --git a/tools/run_mpi_test.py b/tools/run_mpi_test.py deleted file mode 100644 index ab28156b..00000000 --- a/tools/run_mpi_test.py +++ /dev/null @@ -1,33 +0,0 @@ -#!/usr/bin/env python3 -"""Run one MPI test without baking a particular MPI launcher into the meson build tree.""" - -from __future__ import annotations - -import argparse -import os -import sys -import subprocess - - -def main() -> int: - parser = argparse.ArgumentParser() - parser.add_argument("--launcher", required=True) - parser.add_argument("--processes", type=int, default=2) - parser.add_argument("executable") - parser.add_argument("arguments", nargs="*") - args = parser.parse_args() - command = [args.launcher, "-n", str(args.processes), args.executable, *args.arguments] - environment = dict(os.environ) - environment.setdefault("OMPI_ALLOW_RUN_AS_ROOT", "1") - environment.setdefault("OMPI_ALLOW_RUN_AS_ROOT_CONFIRM", "1") - loopback = "lo0" if sys.platform == "darwin" else "lo" - environment.setdefault("FI_PROVIDER", "sockets") - environment.setdefault("FI_SOCKETS_IFACE", loopback) - environment.setdefault("HYDRA_IFACE", loopback) - if sys.platform == "darwin": - environment.setdefault("MPICH_INTERFACE_HOSTNAME", "127.0.0.1") - return subprocess.run(command, env=environment, check=False).returncode - - -if __name__ == "__main__": - raise SystemExit(main()) diff --git a/tools/test_mfem_build_policy.py b/tools/test_mfem_build_policy.py new file mode 100644 index 00000000..c437b79b --- /dev/null +++ b/tools/test_mfem_build_policy.py @@ -0,0 +1,211 @@ +"""MFEM build-policy regressions; no MFEM sources or downloaded dependencies needed. + +Run with ``python tools/test_mfem_build_policy.py``. The integration test compiles +a tiny C/C++ fixture with local CMake; the cache test never launches a compiler. +""" + +import contextlib +import io +import json +import os +from pathlib import Path +import platform +import shlex +import shutil +import subprocess +import sys +import tempfile +from types import SimpleNamespace +import unittest +from unittest.mock import patch + +import build_mfem_bundle as bundle + + +FEATURES = dict.fromkeys(bundle.CMAKE_FEATURES, False) +CMAKE = shutil.which("cmake") +CC = shutil.which("cc") +CXX = shutil.which("c++") + + +def arguments(root, mode="release", optimization="2", debug=False, ndebug=False): + return [ + "build_mfem_bundle.py", + "--source", str(root / "source"), + "--work-dir", str(root / "bundle" / "work"), + "--prefix", str(root / "prefix"), + "--cmake", CMAKE or "cmake", + "--cc", CC or "cc", "--cxx", CXX or "c++", + "--buildtype", mode, "--optimization", optimization, + "--debug", str(debug).lower(), "--ndebug", str(ndebug).lower(), + "--compiler-syntax", "gcc", "--compiler-id", "clang", + "--jobs", "1", "--host-system", platform.system().lower(), + "--cross-build", "false", "--library-kind", "static", + "--allow-system", "true", "--vendor-dependency-prefix", "false", + *[item for name in FEATURES for item in ("--feature", name + "=OFF")], + ] + + +def parse(root, **kwargs): + with patch.object(sys, "argv", arguments(root, **kwargs)): + return bundle.parse_args() + + +def public_header(prefix, debug=False): + header = prefix / "include" / "mfem" / "config" / "_config.hpp" + header.parent.mkdir(parents=True, exist_ok=True) + header.write_text("#define MFEM_USE_DOUBLE\n" + + ("#define MFEM_DEBUG\n" if debug else "")) + return header + + +class MfemBuildPolicyTests(unittest.TestCase): + def test_assertions_symbols_and_optimization_are_independent(self): + args = SimpleNamespace(compiler_syntax="gcc", compiler_id="clang", + optimization="2", debug="false", ndebug="false") + self.assertEqual(shlex.split(bundle.mfem_compile_flags(args)), ["-O2"]) + args.debug = "true" + self.assertNotIn("-DNDEBUG", shlex.split(bundle.mfem_compile_flags(args))) + self.assertIn("-g", shlex.split(bundle.mfem_compile_flags(args))) + args.optimization = "s" + self.assertIn("-Oz", shlex.split(bundle.mfem_compile_flags(args))) + args.compiler_id = "gcc" + self.assertIn("-Os", shlex.split(bundle.mfem_compile_flags(args))) + args.optimization = "plain" + args.debug = "false" + self.assertEqual(bundle.mfem_compile_flags(args), "") + args.compiler_syntax = "msvc" + args.optimization = "2" + args.ndebug = "true" + self.assertEqual(shlex.split(bundle.mfem_compile_flags(args)), + ["/O2", "/DNDEBUG"]) + + def test_public_header_checks_mfem_debug_in_both_directions(self): + with tempfile.TemporaryDirectory() as directory: + root = Path(directory) + for mode, debug in [("debug", True), ("release", False), + ("debugoptimized", False)]: + with self.subTest(mode=mode): + args = parse(root, mode=mode) + public_header(args.prefix, debug) + bundle.validate_mfem_feature_header(args, FEATURES) + public_header(args.prefix, not debug) + with self.assertRaisesRegex(RuntimeError, "MFEM_DEBUG"): + bundle.validate_mfem_feature_header(args, FEATURES) + + def test_cache_reconfiguration_and_cache_hit_header_validation(self): + with tempfile.TemporaryDirectory() as directory: + root = Path(directory) + source = root / "source" + source.mkdir() + (source / "CMakeLists.txt").write_text("# Stable source fingerprint\n") + + def fake_install(args, features, env, build_type, mpicc, mpicxx): + public_header(args.prefix, build_type == "Debug") + (args.prefix / "include" / "mfem.hpp").write_text("// fixture\n") + (args.prefix / "lib").mkdir(parents=True, exist_ok=True) + (args.prefix / "lib" / "libmfem.a").write_bytes(b"fixture archive") + + with patch.object(bundle, "output", return_value="fixture compiler"), \ + patch.object(bundle, "build_mfem", side_effect=fake_install) as build, \ + patch.object(bundle, "install_licenses"), \ + patch.object(bundle, "sanitize_installed_metadata"), \ + contextlib.redirect_stdout(io.StringIO()): + def invoke(**kwargs): + with patch.object(sys, "argv", arguments(root, **kwargs)): + return bundle.main() + + self.assertEqual(invoke(mode="debug", optimization="0", debug=True), 0) + self.assertEqual(invoke(), 0) # Same directory; release invalidates debug. + self.assertEqual(build.call_count, 2) + self.assertEqual(invoke(), 0) # Unchanged valid prefix is reusable. + self.assertEqual(build.call_count, 2) + header = public_header(root / "prefix", debug=True) + with self.assertRaisesRegex(RuntimeError, "MFEM_DEBUG"): + invoke() # Same fingerprint must not hide an incompatible header. + self.assertEqual(build.call_count, 2) + header.unlink() + with self.assertRaisesRegex(RuntimeError, "_config.hpp.*missing"): + invoke() + + @unittest.skipUnless(CMAKE and CC and CXX and os.name != "nt", + "requires local CMake and GCC-compatible C/C++ compilers") + def test_real_cmake_reconfigure_replaces_implicit_flags_for_both_languages(self): + with tempfile.TemporaryDirectory() as directory: + root = Path(directory) + source = root / "source" + source.mkdir() + (source / "CMakeLists.txt").write_text('''\ +cmake_minimum_required(VERSION 3.18) +project(mfem_policy_fixture LANGUAGES C CXX) +set(CMAKE_EXPORT_COMPILE_COMMANDS ON) +# Match the upstream MFEM build-type rule, including debugoptimized behavior. +if(CMAKE_BUILD_TYPE MATCHES "Debug|debug|DEBUG") + set(MFEM_DEBUG ON) +else() + set(MFEM_DEBUG OFF) +endif() +configure_file(config.hpp.in mfem_config.hpp) +add_executable(policy_probe probe.cpp probe.c) +target_include_directories(policy_probe PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}" "${CMAKE_CURRENT_BINARY_DIR}") +install(TARGETS policy_probe DESTINATION bin) +''') + (source / "config.hpp.in").write_text("#cmakedefine MFEM_DEBUG\n") + checks = '''\ +#include "expected.hpp" +#include "mfem_config.hpp" +#if defined(NDEBUG) != EXPECT_NDEBUG +#error CMake changed the requested assertion policy +#endif +#if defined(MFEM_DEBUG) != EXPECT_MFEM_DEBUG +#error Stale MFEM debug configuration +#endif +#if defined(__OPTIMIZE__) != EXPECT_OPTIMIZED +#error Requested optimization was not applied +#endif +''' + (source / "probe.c").write_text(checks + "int c_probe(void) { return 17; }\n") + (source / "probe.cpp").write_text(checks + '''\ +extern "C" int c_probe(void); +int main() { return c_probe() == 17 ? 0 : 1; } +''') + + def quiet_run(command, *, cwd=None, env=None): + result = subprocess.run([os.fspath(x) for x in command], cwd=cwd, + env=env, stdout=subprocess.PIPE, + stderr=subprocess.STDOUT, text=True, timeout=20) + self.assertEqual(result.returncode, 0, result.stdout) + + cases = [("debug", "0", True, False), + ("release", "3", False, False), + ("release", "2", False, False), + ("release", "2", False, True), + ("debugoptimized", "2", True, False), + ("debug", "0", True, False)] + with patch.object(bundle, "prepare_mfem_source", return_value=source), \ + patch.object(bundle, "run", side_effect=quiet_run): + for mode, optimization, debug, ndebug in cases: + with self.subTest(mode=mode, ndebug=ndebug): + args = parse(root, mode=mode, optimization=optimization, + debug=debug, ndebug=ndebug) + (source / "expected.hpp").write_text( + f"#define EXPECT_NDEBUG {int(ndebug)}\n" + f"#define EXPECT_MFEM_DEBUG {int(mode == 'debug')}\n" + f"#define EXPECT_OPTIMIZED {int(optimization != '0')}\n") + bundle.build_mfem(args, FEATURES, os.environ.copy(), + bundle.cmake_build_type(mode), "", "") + build = args.work_dir / "mfem" + commands = json.loads((build / "compile_commands.json").read_text()) + self.assertEqual({Path(x["file"]).suffix for x in commands}, {".c", ".cpp"}) + for command in commands: + tokens = shlex.split(command["command"]) + self.assertEqual([token for token in tokens + if token.startswith("-O")], + ["-O" + optimization]) + self.assertEqual("-DNDEBUG" in tokens, ndebug) + self.assertEqual("-g" in tokens, debug) + quiet_run([args.prefix / "bin" / "policy_probe"]) + + +if __name__ == "__main__": + unittest.main() diff --git a/tools/validate_ester.py b/tools/validate_ester.py new file mode 100644 index 00000000..0033ec65 --- /dev/null +++ b/tools/validate_ester.py @@ -0,0 +1,335 @@ +#!/usr/bin/env python3 +"""Run fresh SERiF solves and compare the unchanged physical requests with the archived ESTER study. + +Example: + python tools/validate_ester.py --report ../MeanField/reports/ester_validation_2026-09-21 \ + --executable build/tests/solver/ester_validation --output /tmp/serif-ester + +Omit --executable to analyze already generated SERiF output. No archived solver +field is used as a new measurement. Archive coefficients are warm initial guesses; +solve.json must demonstrate independent SERiF nonlinear convergence. +""" +import argparse +import csv +import hashlib +import importlib.util +import json +import math +from pathlib import Path +import subprocess +import sys + +import numpy as np + +CASES = ("n1_J0", "n1_J0p05", "n1_J0p1", "n1p5_J0", "n1p625_J0") +SCALARS = {"rho": "rho", "h": "h", "P": "pressure_h", "Phi": "phi", "bernoulli_defect": "bernoulli_defect"} + + +def read_rows(path): + """Retain archived point order and names, including failed point locations.""" + with path.open() as source: + return list(csv.DictReader(source)) + + +def column(rows, name): + """Read one named numerical measurement from every requested point.""" + return np.asarray([float(row[name]) for row in rows]) + + +def magnitude(values): + """Use absolute scalar values or Euclidean vector lengths for error norms.""" + return np.linalg.norm(values, axis=1) if values.ndim == 2 else np.abs(values) + + +def norms(values, mask, scale, weights=None): + """Equal-weight requested-grid norms, as specified by the original protocol.""" + selected = magnitude(values[mask]) / scale + if not len(selected) or not np.isfinite(selected).all(): + raise ValueError("A required physical comparison contains missing or nonfinite samples") + selected_weights = np.ones(len(selected)) if weights is None else weights[mask] + if not np.isfinite(selected_weights).all() or np.any(selected_weights <= 0): + raise ValueError("A required physical comparison has invalid quadrature weights") + return {"count": len(selected), "rms_scaled": float(np.sqrt(np.average(selected**2, weights=selected_weights))), + "max_scaled": float(np.max(selected))} + + +def digest(path): + """Record the exact bytes of a source or independently generated measurement.""" + return hashlib.sha256(path.read_bytes()).hexdigest() + + +def compare_surface(report, case, output, models): + """Repeat the original actual-surface quadrature and sampled-axis measurements.""" + old = read_rows(report / "meanfield_runs" / case / "surface.csv") + new = read_rows(output / "surface.csv") + identity = ("kind", "boundary_element", "sample", "reference_u", "reference_v") + if [[row[key] for key in identity[:3]] for row in old] != [[row[key] for key in identity[:3]] for row in new]: + raise ValueError(f"{case}: actual surface sample identities changed") + for key in identity[3:]: + if not np.array_equal(column(old, key), column(new, key)): + raise ValueError(f"{case}: actual surface reference requests changed") + base = models["base"] + scale = base.globals["R_eq"] + kind = np.asarray([row["kind"] for row in old]) + old_radius, radius = column(old, "radius"), column(new, "radius") + xyz = np.stack([column(new, c) for c in ("x", "y", "z")], axis=1) + if not np.isfinite(xyz).all() or not np.isfinite(radius).all() or np.any(radius <= 0): + raise ValueError(f"{case}: a required actual surface sample is nonfinite or nonpositive") + if not np.allclose(np.linalg.norm(xyz, axis=1), radius, rtol=32*np.finfo(float).eps, atol=1e-15): + raise ValueError(f"{case}: actual surface positions and radii disagree") + reference = {level: model.surface_radius(column(new, "theta")) for level, model in models.items()} + old_reference = base.surface_radius(column(old, "theta")) + archived = json.loads((report / "analysis" / case / "metrics.json").read_text())["surface"] + quadrature = {} + for rule in ("q14", "q18"): + mask = kind == rule + fresh = norms(radius - reference["base"], mask, scale, column(new, "area_weight")) + reconstructed = norms(old_radius - old_reference, mask, scale, column(old, "area_weight")) + recorded = {key: archived[rule][key] for key in ("rms_scaled", "max_scaled")} + if fresh["count"] != archived[rule]["eligible"] or any(abs(reconstructed[key] - recorded[key]) > 128*np.finfo(float).eps for key in recorded): + raise ValueError(f"{case}: archived {rule} surface metric did not reproduce") + sensitivity = {level: norms(reference[level] - reference["base"], mask, scale, column(new, "area_weight")) for level in ("low", "high")} + quadrature[rule] = {"fresh": fresh, "archived": recorded, "archived_reconstructed": reconstructed, + "fresh_minus_archived": norms(radius - old_radius, mask, scale, column(old, "area_weight")), + "area": float(column(new, "area_weight")[mask].sum()), "archived_area": archived[rule]["area"], + "ester_grid_change": sensitivity, + "rms_change": fresh["rms_scaled"] - recorded["rms_scaled"], "max_change": fresh["max_scaled"] - recorded["max_scaled"], + "within_archived_plus_high_grid_change": all(fresh[key] <= recorded[key] + sensitivity["high"][key] for key in recorded)} + # These are selected actual nodes, never extrema or ideal harmonic radii. + axis_tolerance = 2e-12 * scale + old_xyz = np.stack([column(old, c) for c in ("x", "y", "z")], axis=1) + axes = {} + # Freeze the original node visits and exact-coordinate deduplication. Repeating that + # deduplication on fresh floating-point positions can change weights through roundoff alone. + for name, select in {"equator": np.abs(old_xyz[:, 2]) <= axis_tolerance, + "north_pole": (np.hypot(old_xyz[:, 0], old_xyz[:, 1]) <= axis_tolerance) & (old_xyz[:, 2] > 0), + "south_pole": (np.hypot(old_xyz[:, 0], old_xyz[:, 1]) <= axis_tolerance) & (old_xyz[:, 2] < 0)}.items(): + indices = np.flatnonzero((kind == "node") & select) + if not len(indices): + raise ValueError(f"{case}: missing required {name} actual surface nodes") + _, unique = np.unique(old_xyz[indices], axis=0, return_index=True) + indices = indices[unique] + transverse = np.abs(xyz[indices, 2]) if name == "equator" else np.hypot(xyz[indices, 0], xyz[indices, 1]) + if np.any(transverse > axis_tolerance): + raise ValueError(f"{case}: a required {name} surface request moved off its physical axis") + values = radius[indices] + axes[name] = {"mean": float(values.mean()), "minimum": float(values.min()), "maximum": float(values.max()), + "archived_node_visits": int(((kind == "node") & select).sum()), "fixed_archived_unique_coordinate_samples": len(values), + "coordinate_tolerance": axis_tolerance} + return {"quadrature": quadrature, "axis_samples": axes} + + +def compare_globals(report, case, output, old, new, models, surface): + """Compare measured center traces, typed global coefficients and actual surface nodes.""" + center = [index for index, row in enumerate(old) if row["id"] == "center:00000"] + if len(center) != 1 or any(float(new[center[0]][c]) != 0 for c in ("x", "y", "z")): + raise ValueError(f"{case}: the explicit center trace is missing") + index = center[0] + summary = json.loads((output / "physical_summary.json").read_text()) + measured = {name: summary[name] for name in ("Omega", "C", "surface_volume_radius_q14", "surface_volume_radius_q18")} + measured.update(rho_c=float(new[index]["rho"]), h_c=float(new[index]["h"])) + bulk_rows = read_rows(output / "bulk.csv") + bulk = {(row["quadrature"], row["metric"]): float(row["value"]) for row in bulk_rows} + if len(bulk) != len(bulk_rows): + raise ValueError(f"{case}: duplicate physical bulk measurements") + for quadrature in ("14", "18"): + for quantity, metric in {"M": "mass", "I": "moment_of_inertia_z", "J": "angular_momentum_z", "R_volume": "volume_radius"}.items(): + measured[quantity + "_q" + quadrature] = bulk[(quadrature, metric)] + for axis, values in surface["axis_samples"].items(): + quantity = ("R_eq" if axis == "equator" else "R_p") + "_sampled_" + axis + "_mean" + measured[quantity] = values["mean"] + records = {(row["quantity"], row["quadrature"]): row for row in read_rows(report / "analysis" / "global_comparison.csv") if row["case"] == case} + result = {} + for name, value in measured.items(): + if name.startswith("surface_volume_radius_q"): + recorded = records[("R_volume", name.rsplit("q", 1)[1])] + quantity = "R_volume" + elif name.endswith(("_q14", "_q18")): + quantity, quadrature = name.rsplit("_q", 1) + recorded = records[(quantity, quadrature)] + else: + recorded = records[(name, "14" if name in ("Omega", "C") else "")] + quantity = name.split("_sampled_")[0] + scale = float(recorded["normalization"]) + if not math.isfinite(value): + raise ValueError(f"{case}: {name} is nonfinite") + reference = {level: model.globals[quantity] for level, model in models.items()} + previous = float(recorded["meanfield"]) + if any(abs(reference[level] - float(recorded["ester_" + level])) > 128*np.finfo(float).eps*max(1, abs(reference[level])) for level in reference): + raise ValueError(f"{case}: archived {name} ESTER global reference no longer reproduces") + fresh_error = (value - reference["base"]) / scale + old_error = (previous - reference["base"]) / scale + sensitivity = {level: (reference[level] - reference["base"]) / scale for level in ("low", "high")} + result[name] = {"fresh": value, "archived": previous, "ester": reference, "normalization": scale, + "fresh_scaled_difference": fresh_error, "archived_scaled_difference": old_error, + "fresh_minus_archived_scaled": (value - previous) / scale, "ester_grid_change_scaled": sensitivity, + "within_archived_plus_high_grid_change": bool(abs(fresh_error) <= abs(old_error) + abs(sensitivity["high"]))} + return {"quantities": result, "center_boundary_trace_possible": bool(float(new[index]["boundary_trace_possible"])), + "surface_volume_measurement": "Fresh divergence-theorem boundary quadrature compared with original volume-integral radii.", + "physical_summary": summary} + + +def compare(report, case, output, reference_class): + """Keep the original common-material masks fixed so missing samples cannot improve a score.""" + solve = json.loads((output / "solve.json").read_text()) + if solve["case"] != case or solve["converged"] is not True or not math.isfinite(solve["final_residual_norm"]) or not 0.0 <= solve["final_residual_norm"] <= 1e-10: + raise ValueError(f"{case}: fresh SERiF nonlinear qualification failed") + metadata = dict(line.split("=", 1) for line in (report / "inputs" / case / "metadata.txt").read_text().splitlines() if "=" in line) + if digest(report / "inputs" / case / "mesh.smesh") != metadata["mesh_sha256"]: + raise ValueError(f"{case}: the reference mesh no longer matches its recorded hash") + catalog = json.loads((report / "reference_catalog.json").read_text())[case] + old = read_rows(report / "meanfield_runs" / case / "samples.csv") + new = read_rows(output / "samples.csv") + if [row["id"] for row in old] != [row["id"] for row in new]: + raise ValueError(f"{case}: point identity or order changed") + xyz = np.stack([column(new, c) for c in ("x", "y", "z")], axis=1) + if not np.array_equal(xyz, np.stack([column(old, c) for c in ("x", "y", "z")], axis=1)): + raise ValueError(f"{case}: requested physical coordinates changed") + old_located = column(old, "located") == 1 + new_located = column(new, "located") == 1 + if np.any(old_located & ~new_located): + raise ValueError(f"{case}: {np.count_nonzero(old_located & ~new_located)} previously located samples were lost") + def physical_coordinates(rows, located): + """Validate actual located positions against their unchanged physical requests.""" + mapped = np.stack([column(rows, "mapped_" + c) for c in ("x", "y", "z")], axis=1) + error = np.linalg.norm(mapped - xyz, axis=1) + tolerance = 2e-11 * np.maximum(1.0, np.linalg.norm(xyz, axis=1)) + roundoff = 64 * np.finfo(float).eps * np.maximum(1.0, np.linalg.norm(xyz, axis=1)) + if not np.isfinite(mapped[located]).all() or np.any(error[located] > tolerance[located] + roundoff[located]): + raise ValueError(f"{case}: sampled physical location exceeds the original locator tolerance") + if np.any(np.abs(error[located] - column(rows, "location_error")[located]) > roundoff[located]): + raise ValueError(f"{case}: sampled location error is inconsistent with mapped coordinates") + coordinates = xyz.copy() + coordinates[located] = mapped[located] + return coordinates + old_xyz = physical_coordinates(old, old_located) + new_xyz = physical_coordinates(new, new_located) + models = {level: reference_class(report / "ester_runs" / Path(catalog[level]["directory"]).name, + n=float(metadata["n"]), K=float(metadata["K"]), M=float(metadata["M"]), G=float(metadata["G"])) + for level in ("base", "low", "high")} + if not all(model.globals["reference_converged"] for model in models.values()): + raise ValueError(f"{case}: an ESTER reference did not meet its recorded convergence criterion") + match_tolerance = json.loads((report / "matching_plan.json").read_text())["angular_momentum_absolute_tolerance"] + if any(abs(model.globals["J"] - float(metadata["J"])) > match_tolerance + 16 * np.finfo(float).eps for model in models.values()): + raise ValueError(f"{case}: an ESTER reference does not match the physical angular momentum") + reference = {level: model.evaluate(new_xyz) for level, model in models.items()} + old_reference = models["base"].evaluate(old_xyz) + base = models["base"] + radii = np.linalg.norm(old_xyz, axis=1) + theta = np.arccos(np.clip(np.divide(old_xyz[:, 2], radii, out=np.ones(len(radii)), where=radii > 0), -1, 1)) + common = old_located & (column(old, "stellar_material") == 1) & old_reference["inside_star"] + strict = common & (column(old, "boundary_trace_possible") == 0) & (radii / base.surface_radius(theta) < 1 - 1e-8) + if np.any(strict & (column(new, "stellar_material") != 1)): + raise ValueError(f"{case}: a required interior physical request changed material support") + if np.any(strict & (column(new, "boundary_trace_possible") != 0)): + raise ValueError(f"{case}: a required interior sample moved onto an ambiguous FE trace") + scales = {"rho": base.rho_c, "h": base.H, "P": base.rho_c * base.H / (base.n + 1), + "Phi": base.G * base.M / base.globals["R_eq"], "gradPhi": base.G * base.M / base.globals["R_eq"]**2, + "gradPhi_broken": base.G * base.M / base.globals["R_eq"]**2, "bernoulli_defect": base.H} + recorded_metrics = {(row["category"], row["region"], row["field"]): row + for row in read_rows(report / "analysis" / "metrics.csv") if row["case"] == case and row["group"] == "all"} + + def field_values(rows, name): + """Read the specified physical scalar or one of the two gravity gradients.""" + if name in ("gradPhi", "gradPhi_broken"): + prefix = "mixed_dphi_d" if name == "gradPhi" else "broken_phi_d" + return np.stack([column(rows, prefix + c) for c in ("x", "y", "z")], axis=1) + return column(rows, SCALARS[name]) + + def field_comparison(name, mask, category, region="all"): + """Reconstruct the exact archived metric before comparing a new independently sampled field.""" + archived_name = "gradPhi_mixed" if name == "gradPhi" else name + reference_name = "gradPhi" if name == "gradPhi_broken" else name + recorded = recorded_metrics[(category, region, archived_name)] + if int(recorded["eligible"]) != int(mask.sum()): + raise ValueError(f"{case}: {category}/{region} support no longer reproduces the original study") + old_field, new_field = field_values(old, name), field_values(new, name) + fresh = norms(new_field - reference["base"][reference_name], mask, scales[name]) + preserved = norms(new_field - old_field, mask, scales[name]) + sensitivity = {level: norms(reference[level][reference_name] - reference["base"][reference_name], mask, scales[name]) for level in ("low", "high")} + archived = {key: float(recorded[key]) for key in ("rms_scaled", "max_scaled")} + reconstructed = norms(old_field - old_reference[reference_name], mask, scales[name]) + if any(abs(reconstructed[key] - archived[key]) > 128 * np.finfo(float).eps * max(1.0, abs(archived[key])) + for key in ("rms_scaled", "max_scaled")): + raise ValueError(f"{case}: archived {category}/{region}/{name} metrics did not reproduce at their original mapped coordinates") + return {"fresh": fresh, "archived": archived, "archived_reconstructed": reconstructed, + "reference_relocation": norms(reference["base"][reference_name] - old_reference[reference_name], mask, scales[name]), + "fresh_minus_archived": preserved, "ester_grid_change": sensitivity, + "rms_change": fresh["rms_scaled"] - archived["rms_scaled"], + "max_change": fresh["max_scaled"] - archived["max_scaled"], + "within_archived_plus_high_grid_change": all(fresh[key] <= archived[key] + sensitivity["high"][key] + for key in ("rms_scaled", "max_scaled"))} + + fields = {name: field_comparison(name, strict, "common_material_strict") for name in (*SCALARS, "gradPhi", "gradPhi_broken")} + # Preserve the original radial-region definitions, including gravity in physical vacuum. + fraction = radii / base.surface_radius(theta) + regions = {"all": old_located, "core": old_located & (fraction < .5), + "envelope": old_located & (fraction >= .5) & (fraction < .9), + "surface_layer": old_located & (fraction >= .9) & (fraction <= 1), + "outside_reference": old_located & (fraction > 1)} + gravity = {region: {name: field_comparison(name, mask, "located_all_materials", region) for name in ("Phi", "gradPhi", "gradPhi_broken")} + for region, mask in regions.items() if mask.any()} + def material_diagnostics(rows, located): + """Expose negative material tails and undefined pressure without clipping them.""" + material = located & (column(rows, "stellar_material") == 1) + density, enthalpy, pressure = (column(rows, name) for name in ("rho", "h", "pressure_h")) + if not np.isfinite(density[material]).all() or not np.isfinite(enthalpy[material]).all(): + raise ValueError(f"{case}: material density or enthalpy contains a nonfinite value") + return {"material_points": int(material.sum()), "negative_density_points": int((material & (density < 0)).sum()), + "negative_enthalpy_points": int((material & (enthalpy < 0)).sum()), + "nonfinite_pressure_points": int((material & ~np.isfinite(pressure)).sum()), + "minimum_density": float(density[material].min()), "minimum_enthalpy": float(enthalpy[material].min())} + tails = {"archived": material_diagnostics(old, old_located), "fresh": material_diagnostics(new, new_located)} + surface = compare_surface(report, case, output, models) + globals_comparison = compare_globals(report, case, output, old, new, models, surface) + provenance = {str(path): digest(path) for path in [report / "inputs" / case / "accepted_state.txt", + report / "inputs" / case / "mesh.smesh", output / "accepted_state.txt", output / "samples.csv", report / "ester_reference.py", + report / "inputs" / case / "layout.csv", report / "inputs" / case / "metadata.txt", + output / "physical_summary.json", output / "surface.csv", output / "bulk.csv"]} + result = {"case": case, "fresh_solve": solve, "strict_common_material_points": int(strict.sum()), + "located_points": int(new_located.sum()), "fields": fields, "located_gravity_by_region": gravity, "surface": surface, "globals": globals_comparison, "material_tail_diagnostics": tails, "sha256": provenance, + "ester_reference_sha256": {level: model.hashes for level, model in models.items()}, + "qualification": "Fresh nonlinear atol=1e-10; fixed archived physical requests, strict material masks and all-located gravity masks; no universal physical percentage tolerance exists in the source protocol.", + "envelope": "Archived error plus the separately measured high-minus-base ESTER grid change is reported, not substituted for the recorded errors."} + qualified = list(fields.values()) + [item for region in gravity.values() for item in region.values()] + qualified += list(surface["quadrature"].values()) + list(globals_comparison["quantities"].values()) + result["within_archived_plus_high_grid_change"] = all(item["within_archived_plus_high_grid_change"] for item in qualified) + (output / "comparison.json").write_text(json.dumps(result, indent=2, allow_nan=False) + "\n") + return result + + +def main(): + """Optionally solve selected stars, then write their independently checked comparisons.""" + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--report", type=Path, required=True) + parser.add_argument("--output", type=Path, required=True) + parser.add_argument("--executable", type=Path) + parser.add_argument("--cases", nargs="+", choices=CASES, default=CASES) + args = parser.parse_args() + args.output.mkdir(parents=True, exist_ok=True) + sys.dont_write_bytecode = True + specification = importlib.util.spec_from_file_location("serif_ester_reference", args.report / "ester_reference.py") + reference = importlib.util.module_from_spec(specification) + specification.loader.exec_module(reference) + results = [] + for case in args.cases: + target = args.output / case + target.mkdir(parents=True, exist_ok=True) + if args.executable: + with (target / "solve.log").open("w") as log: + subprocess.run([str(args.executable.resolve()), str(args.report.resolve()), case, str(target.resolve())], + stdout=log, stderr=subprocess.STDOUT, check=True) + result = compare(args.report, case, target, reference.Model) + results.append(result) + print(case, "converged", result["fresh_solve"]["final_residual_norm"], "Newton steps", result["fresh_solve"]["accepted_iterations"], flush=True) + for field, data in result["fields"].items(): + print(f" {field}: RMS {data['fresh']['rms_scaled']:.9g}; archived {data['archived']['rms_scaled']:.9g}; max physical drift {data['fresh_minus_archived']['max_scaled']:.3g}", flush=True) + for rule, data in result["surface"]["quadrature"].items(): + print(f" surface {rule}: RMS {data['fresh']['rms_scaled']:.9g}; archived {data['archived']['rms_scaled']:.9g}; max radius drift {data['fresh_minus_archived']['max_scaled']:.3g}", flush=True) + print(" all reported errors within archived plus measured ESTER high-grid change:", + result["within_archived_plus_high_grid_change"], flush=True) + (args.output / "comparison.json").write_text(json.dumps(results, indent=2, allow_nan=False) + "\n") + + +if __name__ == "__main__": + main() diff --git a/tools/validate_matrix.sh b/tools/validate_matrix.sh index 000b363e..306df2e9 100755 --- a/tools/validate_matrix.sh +++ b/tools/validate_matrix.sh @@ -1,5 +1,5 @@ #!/bin/sh -# Build and run all six sandboxes with each Linux compiler/library/linker pair. +# Build and run all sandboxes with each Linux compiler/library/linker pair. set -eu project_root=$(CDPATH= cd -- "$(dirname -- "$0")/.." && pwd)