"""Exact multiple-reflection benchmark using the public LESS API.""" from pathlib import Path import argparse import json import less import numpy as np DATA = Path(__file__).resolve().parents[1] / "data/multiple-scattering" def reference(config, depth=None): r, b, t, g = [np.asarray(config[key], dtype=float) for key in ("front_reflectance", "back_reflectance", "transmittance", "ground_reflectance")] if depth is None: denominator = 1 - g * b return np.array([r + t*t*g/denominator, 1-(r+t) + t*g*(1-(b+t))/denominator, t*(1-g)/denominator, np.zeros_like(r)]) ground_visits, underside_visits = depth // 2, (depth - 1) // 2 q = g*b ground_sum = (1-q**ground_visits)/(1-q) underside_sum = (1-q**underside_visits)/(1-q) remaining = t*q**underside_visits if depth % 2 else t*g*q**(ground_visits-1) return np.array([r+t*t*g*underside_sum, 1-(r+t)+t*g*(1-(b+t))*underside_sum, t*(1-g)*ground_sum, remaining]) def identity(product, backend, workflow): value = product.metadata["execution"] assert value["backend"] == backend and value["native"] is True assert value["workflow"] == workflow assert value["executor"] == f"{backend}.slang.{workflow}" return {key: value[key] for key in ("backend", "workflow", "executor", "native", "random_seed", "artifact_sha256")} def run(config): rows = [] expected = reference(config) for backend in ("optix", "vulkan", "embree"): 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.Lambertian( wavelengths=config["wavelengths_nm"], reflectance=config["ground_reflectance"])) sheet = less.Object("sheet", mesh=str(DATA / "sheet.obj")) sheet.set_property(less.Lambertian(wavelengths=config["wavelengths_nm"], reflectance=config["front_reflectance"], back_reflectance=config["back_reflectance"], transmittance=config["transmittance"])) scene.add(sheet, positions=[[0, 0, config["sheet_height_m"]]]) scene.illumination = less.Illumination( source=less.Sun(zenith=config["sun_zenith_degrees"], wavelengths=config["wavelengths_nm"], irradiance=config["normal_irradiance_W_m2_nm"]), atmosphere=less.NoAtmosphere()) previous = None plans = [(n, config["full_interaction_limit"], True) for n in config["photon_budgets"]] plans += [(config["photon_budgets"][-1], d, False) for d in config["interaction_limits"]] for photons, depth, full in plans: field = scene.simulate(less.ForwardRadiationFieldSensor( bands=config["wavelengths_nm"], mode="grid", grid=(1,1), layer_heights=[0,2], photons_per_m2=photons//4), scattering_order=depth) incident = field.incident_energy.astype(float) actual = np.array([field.escaped_energy, field.absorption_grid.sum(axis=(0,1,2), dtype=float), field.terrain_absorbed, field.truncated_energy]) / incident target = expected if full else reference(config, depth) error = float(np.max(abs(actual-target))) closure = float(np.max(abs(actual.sum(axis=0)-1))) row = dict(backend=backend, photons=photons, depth=depth, full=full, fractions=actual.tolist(), reference=target.tolist(), error=error, closure_error=closure, execution=identity(field,backend,"optical_forward")) passed = error <= config["thresholds"]["fraction_absolute_error"] passed &= closure <= config["thresholds"]["closure_absolute_error"] if full: brf = scene.simulate(less.BRFSensor(bands=config["wavelengths_nm"], directions=config["view_directions"], num_photons=photons), scattering_order=depth) row["brf"] = brf.brf.tolist() row["brf_execution"] = identity(brf, backend, "optical_brf") row["brf_error"] = float(np.max(abs(brf.brf-expected[0]))) counts = brf.metadata["transport_counts"] assert counts["emitted_photons"] == photons assert counts["physical_events"] > photons row["physical_events"] = counts["physical_events"] passed &= row["brf_error"] <= config["thresholds"]["brf_absolute_error"] passed &= np.max(actual[3]) <= config["thresholds"]["full_truncated_fraction"] current = np.r_[actual.ravel(), brf.brf.ravel()] if previous is not None: row["budget_change"] = float(np.max(abs(current-previous))) passed &= row["budget_change"] <= config["thresholds"]["budget_change_absolute"] previous = current row["passed"] = bool(passed) rows.append(row) print(backend, photons, depth, error, row["passed"], flush=True) adopted = [r for r in rows if not r["full"] or r["photons"] == config["photon_budgets"][-1]] return dict(input=config, less_version=less.__version__, reference=expected.tolist(), records=rows, passed=all(r["passed"] for r in adopted)) def figures(result, output): import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt plt.rcParams.update({"svg.fonttype":"none", "font.size":10}) config=result["input"] fig,axes=plt.subplots(1,2,figsize=(10,3.8),layout="constrained") expected=np.asarray(result["reference"]) for band,wavelength in enumerate(config["wavelengths_nm"]): depths=config["interaction_limits"] exact=[reference(config,d)[0,band] for d in depths] axes[0].plot(depths,exact,"-",label=f"{wavelength} nm analytic") axes[0].axhline(expected[0,band],ls=":",lw=.8,color="#777777") colors=dict(optix="#0072B2",vulkan="#D55E00",embree="#009E73") for backend,color in colors.items(): partial=[r for r in result["records"] if r["backend"]==backend and not r["full"]] axes[0].plot([r["depth"] for r in partial],np.array([r["fractions"][0] for r in partial]), ".",color=color,label=[backend,None,None]) full=[r for r in result["records"] if r["backend"]==backend and r["full"]] axes[1].plot([r["photons"] for r in full], [max(r["error"],r["brf_error"]) for r in full],"o-",color=color,label=backend) axes[0].set(xscale="log",xlabel="Maximum surface interactions",ylabel="Escaping fraction", title="Multiple-reflection series (dotted: infinite sum)") axes[1].axhline(.005,color="#444444",ls="--",label="Error limit") axes[1].set(xscale="log",xlabel="Photon budget",ylabel="Max absolute flux / BRF error",title="Sampling refinement") for ax in axes: ax.grid(alpha=.2); ax.legend(fontsize=7) fig.savefig(output/"results.svg");plt.close(fig) fig,ax=plt.subplots(figsize=(7,3.8),layout="constrained") ax.plot([0,2],[0,0],color="#8b6c3e",lw=6) ax.plot([0,2],[1,1],color="#4f7aa3",lw=4) for start,end in [((.2,1.6),(.4,1.02)),((.4,.98),(.7,.02)),((.7,.02),(1,.98)), ((1,.98),(1.3,.02)),((1.3,.02),(1.6,.98)),((1.6,1.02),(1.85,1.6))]: ax.annotate("",xy=end,xytext=start,arrowprops=dict(arrowstyle="->",color="#ce641c",lw=1.8)) ax.text(.02,1.12,"Lambertian sheet: top r, bottom b, transmission t") ax.text(.02,-.2,"Lambertian ground: reflectance g") ax.set(xlim=(-.1,2.1),ylim=(-.3,1.85),xlabel="X / m (periodic)",ylabel="Z / m", title="Infinite parallel-plane cavity; schematic photon path") fig.savefig(output/"scene.svg");plt.close(fig) def main(): parser=argparse.ArgumentParser(description=__doc__) parser.add_argument("--output",type=Path,default=Path("build/multiple-scattering")) args=parser.parse_args();args.output.mkdir(parents=True,exist_ok=True) result=run(json.loads((DATA/"input.json").read_text())) (args.output/"results.json").write_text(json.dumps(result,indent=2,allow_nan=False),encoding="utf-8") figures(result,args.output) raise SystemExit(0 if result["passed"] else 1) if __name__=="__main__": main()