"""Leaf-temperature validation against radiation equilibrium and independent balances.""" from pathlib import Path from decimal import Decimal, localcontext import argparse import json import less import numpy as np DATA = Path(__file__).resolve().parents[1]/"data/energy-balance" BACKENDS = ("optix", "vulkan", "embree") QUANTITIES = ("S_abs", "L_abs", "L_emit", "R_n", "H", "LE", "G") def leaf_vertices(config, case): angle = np.deg2rad(case.get("leaf_tilt_degrees", 0)) half = config["leaf_side_m"]/2 vertices = np.array([[-half, -half, 0], [half, -half, 0], [half, half, 0], [-half, half, 0]]) vertices[:, 2] = vertices[:, 0]*np.sin(angle) vertices[:, 0] *= np.cos(angle) return vertices+config["leaf_center_m"] def walls(config, lower=0, upper=None): x, y, z = config["box_size_m"] upper = z if upper is None else upper return dict(west=[[0, 0, lower], [0, y, lower], [0, y, upper], [0, 0, upper]], east=[[x, 0, lower], [x, y, lower], [x, y, upper], [x, 0, upper]], north=[[0, 0, lower], [x, 0, lower], [x, 0, upper], [0, 0, upper]], south=[[0, y, lower], [x, y, lower], [x, y, upper], [0, y, upper]]) def point_view_factor(points, polygon, normal): """Projected solid angle of a convex polygon wholly in one hemisphere.""" vectors = np.asarray(polygon)[None, :, :]-points[:, None, :] vectors /= np.linalg.norm(vectors, axis=2, keepdims=True) other = np.roll(vectors, -1, axis=1) cross = np.cross(vectors, other) norm = np.linalg.norm(cross, axis=2) angles = np.arctan2(norm, np.sum(vectors*other, axis=2)) projected = cross@np.asarray(normal)/norm return abs(np.sum(angles*projected, axis=1))/(2*np.pi) def view_factors(config, case, order): """Gauss-Duffy average on each leaf triangle of exact polygon factors.""" assert not case.get("leaf_tilt_degrees", 0) nodes, weights = np.polynomial.legendre.leggauss(order) u, v = np.meshgrid((nodes+1)/2, (nodes+1)/2, indexing="ij") weight = (weights[:, None]*weights[None, :]/4*2*(1-u)).ravel() z = config["leaf_center_m"][2] upper, lower = walls(config, z), walls(config, upper=z) x, y, _ = config["box_size_m"] floor = [[0, 0, 0], [x, 0, 0], [x, y, 0], [0, y, 0]] vertices = leaf_vertices(config, case) result = [] for indices in ((0, 1, 2), (0, 2, 3)): a, b, c = vertices[list(indices)] points = a+u.ravel()[:, None]*(b-a)+((1-u)*v).ravel()[:, None]*(c-a) top = {name: float(weight@point_view_factor(points, rect, [0, 0, 1])) for name, rect in upper.items()} bottom = {name: float(weight@point_view_factor(points, rect, [0, 0, -1])) for name, rect in lower.items()} bottom["floor"] = float(weight@point_view_factor(points, floor, [0, 0, -1])) assert abs(sum(bottom.values())-1) < 1e-12 combined = {name: top[name]+(bottom[name] if case["sides"] == 2 else 0) for name in top} combined["floor"] = bottom["floor"] if case["sides"] == 2 else 0 combined["sky"] = 1-sum(top.values()) assert min(combined.values()) >= 0 assert abs(sum(combined.values())-case["sides"]) < 1e-12 result.append(combined) return result def shortwave(config, case): angle = case.get("sun_zenith_degrees", 0)+case.get("leaf_tilt_degrees", 0) return ((1-config["reflectance"])*case["irradiance"]*config["band_width_nm"]* np.cos(np.deg2rad(angle))) def radiative_temperature(config, case): with localcontext() as context: context.prec = 50 value = (Decimal(str(case["environment_K"]))**4+ Decimal(str(shortwave(config, case)))/Decimal(str(case["sides"]* case["emissivity"]*config["constants"]["sigma"]))) return float(value.sqrt().sqrt()) def climate(config, case): air = (radiative_temperature(config, case)-273.15 if case["stage"] == "radiative" else case["air_C"]) return less.Microclimate(air_temperature=air, humidity=case.get("humidity_percent", 50), wind_speed=case.get("wind_m_s", 1), downwelling_longwave=config["constants"]["sigma"]*case["environment_K"]**4) def reference(config, case): q = shortwave(config, case) constants = config["constants"] sigma, epsilon = constants["sigma"], case["emissivity"] factors, refinement = None, 0.0 longwave = np.full(2, case["sides"]*epsilon*sigma*case["environment_K"]**4) if "east_wall_K" in case: coarse, factors = [view_factors(config, case, n) for n in config["view_factor_orders"]] refinement = max(abs(a[key]-b[key]) for a, b in zip(coarse, factors) for key in a) assert refinement < config["thresholds"]["view_factor_refinement"] longwave += [epsilon*sigma*f["east"]*(case["east_wall_K"]**4- case["environment_K"]**4) for f in factors] met = climate(config, case) ra = case["ra"] heat_sides = 1 if ra is not None else case["sides"] if ra is None: ra = 200*np.sqrt(case.get("leaf_width_m", .05)/max(case.get("wind_m_s", 1), .1)) ch = heat_sides*constants["rho_air"]*constants["cp_air"]/ra ns = dict(none=0, top=1, bottom=1, both=2)[case["stomata"]] ce = ns*constants["rho_air"]*constants["cp_air"]/( constants["psychrometric_Pa_K"]*(ra+case.get("rs", 150))) with localcontext() as context: context.prec = 50 dec = lambda x: Decimal(str(x)) def saturated(t): c = t-dec(273.15) return dec(610.78)*(dec(17.27)*c/(c+dec(237.3))).exp() ta = dec(met.air_temperature_K) ea = saturated(ta)*dec(met.humidity)/100 records = [] for incoming in longwave: def evaluate(t): emitted = dec(case["sides"])*dec(epsilon)*dec(sigma)*t**4 net = dec(q)+dec(incoming)-emitted sensible = dec(ch)*(t-ta) latent = dec(ce)*(saturated(t)-ea) return dict(S_abs=q, L_abs=float(incoming), L_emit=float(emitted), R_n=float(net), H=float(sensible), LE=float(latent), G=0.0) if case["stage"] == "radiative": # Closed form. The ambient temperature is prescribed to # this radiation-equilibrium point, so H and LE vanish there. temperature = dec(radiative_temperature(config, case)) else: lower, upper = dec(200), dec(350) for _ in range(100): middle = (lower+upper)/2 flux = evaluate(middle) residual = flux["R_n"]-flux["H"]-flux["LE"] if residual > 0: lower = middle else: upper = middle temperature = (lower+upper)/2 records.append(dict(temperature_K=float(temperature), **evaluate(temperature))) return dict(elements=records, air_temperature_C=met.air_temperature, ra_s_m=float(ra), heat_coefficient=float(ch), latent_coefficient=float(ce), view_factors=factors, view_factor_refinement=refinement) def scene_for(config, case, backend, initial=None): scene = less.Scene(backend=backend) scene.size = tuple(config["box_size_m"][:2]) scene.repetitive = False scene.terrain = less.Terrain(property=less.Lambertian(reflectance=0)) scene.terrain.set_property(less.ThermalProperty(emissivity=1, temperature=case["environment_K"])) scene.terrain.set_property(less.BiophysicalProperty(two_sided=False, solve=False, stomata_side="none")) enclosure = less.Object("enclosure") for name in walls(config): enclosure.add_component(name, less.Mesh(str(DATA/"walls.obj"), group=name)) enclosure.set_property(name, less.Lambertian(reflectance=0)) enclosure.set_property(name, less.ThermalProperty(emissivity=1, temperature=case.get("east_wall_K", case["environment_K"]) if name == "east" else case["environment_K"])) enclosure.set_property(name, less.BiophysicalProperty(solve=False, stomata_side="none")) scene.add(enclosure, positions=[(0, 0, 0)]) leaf = less.Object("leaf", mesh=str(DATA/("leaf-tilted.obj" if case.get("leaf_tilt_degrees", 0) else "leaf.obj"))) leaf.set_property(less.Lambertian(reflectance=config["reflectance"])) leaf.set_property(less.ThermalProperty(emissivity=case["emissivity"], temperature=config["initial_temperature_K"] if initial is None else initial)) leaf.set_property(less.BiophysicalProperty(two_sided=case["sides"] == 2, stomata_side=case["stomata"], stomatal_resistance=case.get("rs", 150), aerodynamic_resistance=case["ra"], leaf_width=case.get("leaf_width_m", .05))) if "formulation" in case: parameters = dict(config["physiology"]["parameters"]) if case["formulation"] == "scope": parameters.pop("Jmax25") leaf.set_property(less.Farquhar(formulation=case["formulation"], absorption_component="total", **parameters)) scene.add(leaf, positions=[tuple(config["leaf_center_m"])]) scene.illumination = less.Illumination(source=less.Sun( zenith=case.get("sun_zenith_degrees", 0), azimuth=config["sun_azimuth_degrees"], irradiance=case["irradiance"]), atmosphere=less.NoAtmosphere()) return scene def execute(config, case, backend, samples, tolerance=None, seed=None, initial=None): scene = scene_for(config, case, backend, initial) profile = less.SpectralProfile([less.SpectralBand(config["wavelength_nm"], config["band_width_nm"], "par")]) seed = config["base_seed"] if seed is None else seed radiation = scene.simulate(less.RadiationFieldProcess(spectral_profile=profile, n_samples=32, max_depth=8), seed=seed) tolerance = tolerance or config["iteration_tolerances"][-1] product = scene.simulate(less.EnergyBalanceProcess(rad=radiation, n_samples_per_face=samples, max_depth=8, max_iterations=config["max_iterations"], convergence_tol=tolerance[0], energy_tolerance=tolerance[1], damping=config["damping"]), microclimate=climate(config, case), seed=seed) result = product.result mask = np.asarray(result.active_mask, bool) assert mask.sum() == 2 execution = product.metadata["execution"] assert execution["backend"] == backend and execution["native"] is True assert execution["workflow"] == "joint_energy_balance" assert execution["executor"] == f"{backend}.slang.energy_balance" values = {name: np.asarray(getattr(result, name), float)[mask].tolist() for name in QUANTITIES} temperature = np.asarray(result.temperatures, float)[mask] closure = np.asarray(values["S_abs"])+values["L_abs"]-values["L_emit"] closure -= np.asarray(values["H"])+values["LE"]+values["G"] return dict(backend=backend, case=case["name"], samples=samples, seed=seed, initial_temperature_K=config["initial_temperature_K"] if initial is None else initial, tolerances=tolerance, temperature_K=temperature.tolist(), fluxes=values, area_m2=np.asarray(result.primitive_areas, float)[mask].tolist(), closure_W_m2=float(np.max(abs(closure))), reported_residual_W_m2=np.asarray(result.residual, float)[mask].tolist(), iterations=int(result.n_iterations), converged=bool(result.converged), final_temperature_change_K=float(result.final_max_dT), execution={k: execution[k] for k in ("backend", "workflow", "executor", "native")}) def errors(row, exact): elements = exact["elements"] temperature = float(np.max(abs(np.asarray(row["temperature_K"])- [e["temperature_K"] for e in elements]))) flux = {name: float(np.max(abs(np.asarray(row["fluxes"][name])- [e[name] for e in elements]))) for name in QUANTITIES} return dict(temperature_error_K=temperature, flux_errors_W_m2=flux) def scalar_record(row): record = {k: v for k, v in row.items() if k not in ("fluxes", "execution", "area_m2", "reported_residual_W_m2")} record["mean_fluxes_W_m2"] = {k: float(np.mean(v)) for k, v in row["fluxes"].items()} return record def run(config, backends=BACKENDS, selected=None): cases = [c for c in config["cases"] if selected is None or c["name"] == selected] refs = {c["name"]: reference(config, c) for c in cases} rows, convergence, iteration, stability = [], [], [], [] limits = config["thresholds"] for backend in backends: for case in cases: budgets = config["hot_wall_samples"] if "east_wall_K" in case else config["samples"] previous = None for samples in budgets: row = execute(config, case, backend, samples) row.update(errors(row, refs[case["name"]])) change = None if previous is None else float(np.max(abs( np.array(row["temperature_K"])-previous))) row["sampling_change_K"] = change convergence.append(scalar_record(row)) previous = np.array(row["temperature_K"]) print(json.dumps(scalar_record(row)), flush=True) stochastic = "east_wall_K" in case temperature_limit = limits["hot_wall_temperature_K" if stochastic else "temperature_K"] longwave_limit = limits["hot_wall_longwave_W_m2" if stochastic else "longwave_W_m2"] flux_limit = limits["hot_wall_flux_W_m2" if stochastic else "flux_W_m2"] flux = row["flux_errors_W_m2"] row["passed"] = bool(row["converged"] and row["temperature_error_K"] <= temperature_limit and flux["S_abs"] <= limits["shortwave_W_m2"] and flux["L_abs"] <= longwave_limit and max(flux[k] for k in ("L_emit", "R_n", "H", "LE", "G")) <= flux_limit and row["closure_W_m2"] <= limits["closure_W_m2"] and change <= limits["sampling_temperature_change_K"]) rows.append(row) if selected is not None: continue if case["name"] in config["iteration_cases"]: controls = [] for tolerance in config["iteration_tolerances"][:-1]: control = execute(config, case, backend, budgets[-1], tolerance) control.update(errors(control, refs[case["name"]])) controls.append(scalar_record(control)) controls.append(scalar_record(row)) for initial in config["initial_temperature_controls_K"]: control = execute(config, case, backend, budgets[-1], initial=initial) control.update(errors(control, refs[case["name"]])) control["final_difference_K"] = float(np.max(abs( np.array(control["temperature_K"])-row["temperature_K"]))) controls.append(scalar_record(control)) difference = float(np.max(abs(np.array(controls[1]["temperature_K"])- controls[2]["temperature_K"]))) flux_change = max(abs(controls[1]["mean_fluxes_W_m2"][k]- controls[2]["mean_fluxes_W_m2"][k]) for k in QUANTITIES) iteration.append(dict(backend=backend, case=case["name"], controls=controls, final_difference_K=difference, final_flux_change_W_m2=flux_change, passed=difference <= limits["iteration_temperature_change_K"] and flux_change <= limits["iteration_flux_change_W_m2"] and all(c["converged"] and c["closure_W_m2"] <= c["tolerances"][1] for c in controls) and all(c.get("final_difference_K", 0) <= limits["iteration_temperature_change_K"] for c in controls))) if stochastic: repeats = [scalar_record(row)] for seed in config["stability_seeds"][1:]: control = execute(config, case, backend, budgets[-1], seed=seed) control.update(errors(control, refs[case["name"]])) repeats.append(scalar_record(control)) temperatures = np.array([c["temperature_K"] for c in repeats]) mean, std = temperatures.mean(0), temperatures.std(0, ddof=1) expected = np.array([e["temperature_K"] for e in refs[case["name"]]["elements"]]) stability.append(dict(backend=backend, case=case["name"], repeats=repeats, mean_K=mean.tolist(), std_K=std.tolist(), mean_confidence_95_K=(4.302652729911275*std/np.sqrt(len(repeats))).tolist(), passed=bool(np.max(abs(mean-expected)) <= temperature_limit and all(c["temperature_error_K"] <= temperature_limit and c["converged"] and c["closure_W_m2"] <= limits["closure_W_m2"] for c in repeats)))) parity = [] for left_index, left in enumerate(stability): for right in stability[left_index+1:]: if left["case"] != right["case"]: continue difference = float(np.max(abs(np.array(left["mean_K"])-right["mean_K"]))) uncertainty = float(np.max(np.hypot(left["mean_confidence_95_K"], right["mean_confidence_95_K"]))) parity.append(dict(backends=[left["backend"], right["backend"]], case=left["case"], mean_difference_K=difference, uncertainty_K=uncertainty, passed=difference <= uncertainty+limits["temperature_K"])) return dict(input=config, references=refs, records=rows, convergence=convergence, iteration=iteration, stability=stability, statistical_parity=parity, passed=all(r["passed"] for r in rows+iteration+stability+parity)) def main(): parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--output", type=Path, default=Path("build/energy-balance-validation")) parser.add_argument("--backend", choices=BACKENDS) parser.add_argument("--case") parser.add_argument("--figures-only", action="store_true") parser.add_argument("--physiology-only", action="store_true") args = parser.parse_args() args.output.mkdir(parents=True, exist_ok=True) config = json.loads((DATA/"input.json").read_text(encoding="utf-8")) path = args.output/"results.json" if args.figures_only: result = json.loads(path.read_text(encoding="utf-8")) else: if args.physiology_only: result = json.loads((path if path.exists() else DATA/"results.json").read_text(encoding="utf-8")) previous = {k: v for k, v in result["input"].items() if k != "physiology"} assert previous == {k: v for k, v in config.items() if k != "physiology"} result["input"] = config else: result = run(config, [args.backend] if args.backend else BACKENDS, args.case) if args.case is None: from appendix_leaf_physiology import run as run_physiology result["physiology"] = run_physiology(config, [args.backend] if args.backend else BACKENDS) result["passed"] = all(r["passed"] for group in (result["records"], result["iteration"], result["stability"], result["statistical_parity"]) for r in group) and result["physiology"]["passed"] path.write_text(json.dumps(result, indent=2, allow_nan=False), encoding="utf-8") figures(result, args.output) if "physiology" in result: from appendix_leaf_physiology import figures as physiology_figures physiology_figures(result["input"], result["physiology"], args.output) raise SystemExit(0 if result["passed"] else 1) def figures(result, output): import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt plt.rcParams.update({"font.size": 10, "svg.fonttype": "none", "axes.spines.top": False, "axes.spines.right": False}) colors = ["#2563a6", "#d47713", "#27976b"] markers = ["o", "s", "^"] fig, axes = plt.subplots(2, 2, figsize=(12, 8.5), layout="constrained") names = [c["name"] for c in result["input"]["cases"]] reference_t = [np.mean([e["temperature_K"] for e in result["references"][name]["elements"]]) for name in names] axes[0, 0].plot(range(len(names)), reference_t, color="#263442", lw=1.2, label="Reference") for backend, color, marker in zip(BACKENDS, colors, markers): rows = [r for r in result["records"] if r["backend"] == backend] if not rows: continue axes[0, 0].plot([names.index(r["case"]) for r in rows], [np.mean(r["temperature_K"]) for r in rows], ls="none", marker=marker, markerfacecolor="none", markersize=8, color=color, label=backend) levels = [r for r in result["convergence"] if r["backend"] == backend and r["case"] == "hot_wall"] if levels: axes[1, 0].loglog([r["samples"] for r in levels], [r["temperature_error_K"] for r in levels], marker=marker, markerfacecolor="none", color=color, label=backend) controls = [r for r in result["iteration"] if r["backend"] == backend and r["case"] == "hot_wall"] if controls: controls = controls[0]["controls"][:3] axes[1, 1].semilogx([r["tolerances"][0] for r in controls], [r["closure_W_m2"] for r in controls], marker=marker, markerfacecolor="none", color=color, label=backend) axes[0, 0].set(xticks=range(len(names)), xticklabels=[str(i+1) for i in range(len(names))], xlabel="Condition (table order)", ylabel="Mean leaf temperature / K", title="(a) Leaf temperature") axes[0, 0].legend(ncols=2, frameon=False, fontsize=9) selected = [n for n in ("sensible_only", "transpiring", "humid", "two_stomata", "wind_exchange", "hot_wall") if n in names] x = np.arange(len(selected)) for index, (quantity, color) in enumerate(zip(("R_n", "H", "LE"), colors)): values = [np.mean([e[quantity] for e in result["references"][n]["elements"]]) for n in selected] axes[0, 1].bar(x+(index-1)*.24, values, width=.24, color=color, alpha=.75, label={"R_n": "Net radiation", "H": "Sensible", "LE": "Latent"}[quantity]) for backend, marker in zip(BACKENDS, markers): rows = {r["case"]: r for r in result["records"] if r["backend"] == backend} if all(n in rows for n in selected): axes[0, 1].plot(x+(index-1)*.24, [np.mean(rows[n]["fluxes"][quantity]) for n in selected], marker=marker, markerfacecolor="none", color="#202c39", ls="none", ms=5) axes[0, 1].set(xticks=x, xticklabels=[str(names.index(n)+1) for n in selected], xlabel="Condition (table order)", ylabel="Flux / W m$^{-2}$", title="(b) Heat partition: bars reference, symbols LESS") axes[0, 1].legend(frameon=False, fontsize=9) axes[1, 0].axhline(result["input"]["thresholds"]["hot_wall_temperature_K"], color="#6a7280", ls="--", label="Temperature limit") axes[1, 0].set(xlabel="Longwave samples per triangle per side", ylabel="Maximum leaf-temperature error / K", title="(c) Hot-wall sampling convergence") axes[1, 0].legend(frameon=False, fontsize=9) axes[1, 1].invert_xaxis() axes[1, 1].set(xlabel="Temperature iteration tolerance / K (tighter →)", ylabel="Maximum energy residual / W m$^{-2}$", title="(d) Hot-wall iteration convergence") for ax in axes.flat: ax.grid(alpha=.18, axis="y") fig.savefig(output/"results.svg") fig.savefig(output/"results.png", dpi=150) plt.close(fig) if __name__ == "__main__": main()