"""Analytic thermal-radiation references and public-API simulations.""" from pathlib import Path import argparse import json import less import numpy as np DATA = Path(__file__).resolve().parents[1] / "data/thermal-cavity" def planck(wavelengths_nm, temperature_K): wavelength_m = np.asarray(wavelengths_nm, dtype=float) * 1e-9 h, c, k = 6.62607015e-34, 299792458.0, 1.380649e-23 return 2*h*c*c / (wavelength_m**5 * np.expm1(h*c/(wavelength_m*k*temperature_K))) * 1e-9 def brightness_temperature(wavelengths_nm, radiance_per_nm): wavelength_m = np.asarray(wavelengths_nm, dtype=float) * 1e-9 h, c, k = 6.62607015e-34, 299792458.0, 1.380649e-23 return h*c / (wavelength_m*k*np.log1p(2*h*c*c/(wavelength_m**5 * np.asarray(radiance_per_nm)*1e9))) def cavity_reference(wavelengths_nm, ground_K, sheet_K, sky_K, ground_emissivity, sheet_emissivity, sheet_transmittance=0): bg, bs, sky = [planck(wavelengths_nm, temperature) for temperature in (ground_K, sheet_K, sky_K)] eg, es, t = [np.asarray(v, dtype=float) for v in (ground_emissivity, sheet_emissivity, sheet_transmittance)] rg, rs = 1-eg, 1-es-t upward = (eg*bg + rg*es*bs + rg*t*sky) / (1-rg*rs) downward = es*bs + rs*upward + t*sky external = es*bs + rs*sky + t*upward return dict(inside_upward=upward, inside_downward=downward, external_upward=external) def run_blackbody(config): rows = [] for backend in ("optix", "vulkan", "embree"): for temperature in config["blackbody_temperatures_K"]: scene = less.Scene(backend=backend) scene.sampler = less.Sampler(config.get("sampling_method", "independent")) scene.size = tuple(config["size_m"]) scene.terrain = less.Terrain(property=less.ThermalProperty( temperature=temperature, emissivity=1)) scene.illumination = less.Illumination( source=less.Sun(irradiance=0), atmosphere=less.NoAtmosphere()) target = planck(config["wavelengths_nm"], temperature) for quality in config["qualities"]: product = scene.simulate(less.ThermalImager( less.Orthographic(image_size=config["image_size"], extent=1), bands=config["wavelengths_nm"], quality=quality, mode="property")) execution = product.metadata["execution"] assert execution["backend"] == backend and execution["native"] is True assert execution["executor"] == f"{backend}.slang.optical_backward" actual = product.data.mean(axis=(0, 1), dtype=float) brightness = product.to_brightness_temperature(unit="K").data.mean(axis=(0,1), dtype=float) error = float(np.max(np.abs(actual-target)/target)) temperature_error = float(np.max(np.abs(brightness-temperature))) rows.append(dict(backend=backend, temperature_K=temperature, quality=quality, radiance=actual.tolist(), reference=target.tolist(), brightness_K=brightness.tolist(), relative_error=error, temperature_error_K=temperature_error, execution={key:execution[key] for key in ("backend", "workflow", "executor", "native", "artifact_sha256")}, passed=bool(error <= config["thresholds"]["radiance_relative_error"] and temperature_error <= config["thresholds"]["brightness_temperature_error_K"]))) print(backend, temperature, error, temperature_error, flush=True) return rows def run_cavity(config): rows = [] bands = config["wavelengths_nm"] for backend in ("optix", "vulkan", "embree"): for temperatures in config["temperature_sets"]: scene = less.Scene(backend=backend) scene.sampler = less.Sampler(config.get("sampling_method", "independent")) scene.size = tuple(config["size_m"]) scene.repetitive = config["periodic_wrap_limit"] scene.terrain = less.Terrain(property=less.ThermalProperty( temperature=temperatures["ground_K"], wavelengths=bands, emissivity=config["ground_emissivity"])) sheet = less.Object("sheet", mesh=str(DATA / "sheet.obj")) sheet.set_property(less.ThermalProperty( temperature=temperatures["sheet_K"], wavelengths=bands, emissivity=config["sheet_emissivity"], transmittance=config["sheet_transmittance"])) scene.add(sheet, positions=[(0, 0, config["sheet_height_m"])]) scene.illumination = less.Illumination( source=less.Sun(irradiance=0), atmosphere=less.PrescribedAtmosphere( sky_temperature=temperatures["sky_K"])) references = cavity_reference(bands, **{k: temperatures[k] for k in ("ground_K", "sheet_K", "sky_K")}, **{k: config[k] for k in ("ground_emissivity", "sheet_emissivity", "sheet_transmittance")}) for quantity, height, sign in (("inside_upward", .5, -1), ("inside_downward", .5, 1), ("external_upward", 1.5, -1)): target = references[quantity] for angle in config["view_angles_degrees"]: radians = np.deg2rad(angle) projection = less.Perspective(image_size=config["image_size"], fov=config["fov_degrees"], position=(1, 1, height), target=(1+np.sin(radians), 1, height+sign*np.cos(radians))) previous = None for quality in config["qualities"]: product = scene.simulate(less.ThermalImager( projection, bands=bands, quality=quality, mode="property")) execution = product.metadata["execution"] assert execution["native"] is True and execution["backend"] == backend assert execution["executor"] == f"{backend}.slang.optical_backward" actual = product.data.mean(axis=(0, 1), dtype=float) error = float(np.max(np.abs(actual-target)/target)) bt_error = float(np.max(np.abs(brightness_temperature(bands, actual) - brightness_temperature(bands, target)))) change = (None if previous is None else float(np.max(np.abs(actual-previous)/target))) rows.append(dict(backend=backend, temperatures=temperatures["name"], quantity=quantity, angle_degrees=angle, quality=quality, radiance=actual.tolist(), reference=target.tolist(), relative_error=error, temperature_error_K=bt_error, budget_change=change, execution={key:execution[key] for key in ("backend", "workflow", "executor", "native", "artifact_sha256")}, passed=bool(error <= config["thresholds"]["radiance_relative_error"] and bt_error <= config["thresholds"]["brightness_temperature_error_K"] and (quality != config["qualities"][-1] or change <= config["thresholds"]["final_budget_relative_change"])))) previous = actual print(backend, temperatures["name"], quantity, angle, error, bt_error, change, flush=True) return rows def plot_results(result, output): import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt plt.rcParams.update({"font.size": 10, "svg.fonttype": "none"}) config = result["input"] wavelengths = np.asarray(config["wavelengths_nm"]) / 1000 colors = dict(optix="#0072b2", vulkan="#d55e00", embree="#009e73") markers = dict(optix="o", vulkan="x", embree="+") final = [r for r in result["records"] if r["quality"] == config["qualities"][-1]] fig, axes = plt.subplots(1, 3, figsize=(12, 3.6), layout="constrained") for ax, quantity, title in zip(axes, ("inside_upward", "inside_downward", "external_upward"), ("Inside: upward U", "Inside: downward D", "Outside: upward O")): rows = [r for r in final if r["temperatures"] == "unequal" and r["quantity"] == quantity and r["angle_degrees"] == 0] ax.plot(wavelengths, np.asarray(rows[0]["reference"])*1000, "k-", label="Analytic") for row in rows: b = row["backend"] ax.plot(wavelengths, np.asarray(row["radiance"])*1000, color=colors[b], marker=markers[b], ls="none", label=b) ax.set(title=title, xlabel="Wavelength (µm)", ylabel="Radiance (W m⁻² sr⁻¹ µm⁻¹)") ax.grid(alpha=.2) axes[0].legend(frameon=False) fig.savefig(output / "radiance.svg") fig.savefig(output / "radiance.png", dpi=130) plt.close(fig) fig, axes = plt.subplots(1, 3, figsize=(12, 3.7), layout="constrained") for temperature in config["blackbody_temperatures_K"]: r = next(r for r in result["blackbody"] if r["temperature_K"] == temperature) axes[0].plot(wavelengths, np.asarray(r["reference"])*1000, label=f"{temperature} K") for backend, color in colors.items(): controls = [r for r in result["blackbody"] if r["backend"] == backend and r["quality"] == config["qualities"][-1]] for r in controls: axes[0].plot(wavelengths, np.asarray(r["radiance"])*1000, markers[backend], color=color, ms=5) equilibrium = [r for r in final if r["backend"] == backend and r["temperatures"] == "equilibrium"] errors = np.max([np.abs(brightness_temperature(config["wavelengths_nm"], r["radiance"])-300) for r in equilibrium], axis=0) axes[1].plot(wavelengths, errors, marker=markers[backend], color=color, label=backend) changes = [max(r["relative_error"] for r in result["records"] if r["backend"] == backend and r["quality"] == q) for q in config["qualities"]] axes[2].plot(config["qualities"], np.asarray(changes)*100, marker=markers[backend], color=color, label=backend) axes[0].set(title="Blackbody plane control", xlabel="Wavelength (µm)", ylabel="Radiance (W m⁻² sr⁻¹ µm⁻¹)") axes[0].legend(frameon=False) axes[1].axhline(config["thresholds"]["brightness_temperature_error_K"], color="k", ls="--") axes[1].set(title="300 K equilibrium: worst of six views", xlabel="Wavelength (µm)", ylabel="Brightness-temperature error (K)") axes[1].legend(frameon=False) axes[2].axhline(config["thresholds"]["radiance_relative_error"]*100, color="k", ls="--") axes[2].set(title="Cavity: worst of all views and bands", xlabel="Samples per pixel", ylabel="Maximum relative radiance error (%)", xscale="log") axes[2].legend(frameon=False) for ax in axes: ax.grid(alpha=.2) fig.savefig(output / "validation.svg") fig.savefig(output / "validation.png", dpi=130) plt.close(fig) def main(): parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--output", type=Path, default=Path("build/thermal-cavity")) parser.add_argument("--plot-only", type=Path, help="Plot an existing results.json without simulation") args = parser.parse_args() args.output.mkdir(parents=True, exist_ok=True) if args.plot_only: plot_results(json.loads(args.plot_only.read_text(encoding="utf-8")), args.output) return config = json.loads((DATA/"input.json").read_text()) controls = run_blackbody(config) rows = run_cavity(config) final = [row for row in controls+rows if row["quality"] == config["qualities"][-1]] result = dict(input=config, blackbody=controls, records=rows, passed=all(row["passed"] for row in final)) (args.output/"results.json").write_text(json.dumps(result, indent=2), encoding="utf-8") raise SystemExit(0 if result["passed"] else 1) if __name__ == "__main__": main()