"""A8 extension: independent FvCB/Ball--Berry roots and native leaf coupling.""" import math import json from pathlib import Path import numpy as np import less import appendix_energy_balance as cavity def bisect(function, lower, upper): left = function(lower) assert left*function(upper) <= 0, (lower, upper, left, function(upper)) for _ in range(80): middle = (lower+upper)/2 value = function(middle) if left*value > 0: lower, left = middle, value else: upper = middle return (lower+upper)/2 def saturation(temperature): c = temperature-273.15 return 610.78*math.exp(17.27*c/(c+237.3)) def biochemical(config, case, temperature, apar, ci): """Published rate equations, in ppm; no native solver or its iteration.""" p = config["physiology"]["parameters"] constants = config["physiology"]["constants"] t0, gas = 298.15, 8.314 def arr(activation): return math.exp(activation/gas*(1/t0-1/temperature)) def inhibition(entropy, deactivation): return ((1+math.exp((entropy*t0-deactivation)/(gas*t0)))/ (1+math.exp((entropy*temperature-deactivation)/(gas*temperature)))) scope = case["formulation"] == "scope" vcmax = p["Vcmax25"]*arr(65330)*inhibition(485 if scope else 490, 149250 if scope else 150000) rd = p["Vcmax25"]*p["Rd25_ratio"]*arr(46390) if scope: rd *= inhibition(490, 150650) pressure_bar = constants["pressure_Pa"]*1e-5 oxygen = constants["oxygen_mmol_mol"]*.001*pressure_bar gamma = .5*constants["oxygen_mmol_mol"]*.001/2444/1e-6*arr(37830) km = 405/pressure_bar*arr(79430)*(1+oxygen/(.279*arr(36380))) kd = max(.8738, .0301*(temperature-273.15)+.0773) electron = .5*4/(.05+kd+4)*apar ac = vcmax*(ci-gamma)/(ci+km) aj = .2*electron*(ci-gamma)/(ci+2*gamma) # Stable lower roots of the two nonrectangular hyperbolas. def smooth(a, b, theta, negative_branch=False): discriminant = math.sqrt(max(0, (a+b)**2-4*theta*a*b)) if negative_branch: return (a+b+discriminant)/(2*theta) return 2*a*b/(a+b+discriminant) gross = smooth(smooth(ac, aj, .8, ac < 0), .5*vcmax, .98) else: gamma = 42.7*arr(37830) km = 404.9*arr(79430)*(1+209/(278.4*arr(36380))) jmax = p["Jmax25"]*arr(43540)*inhibition(490, 150000) absorbed = .3*apar electron = (2*absorbed*jmax/(absorbed+jmax+ math.sqrt(max(0, (absorbed+jmax)**2-2.8*absorbed*jmax)))) ac = vcmax*(ci-gamma)/(ci+km) aj = .25*electron*(ci-gamma)/(ci+2*gamma) gross = min(ac, aj) return dict(An=gross-rd, Rd=rd, Ac=ac, Aj=aj) def stomatal(config, case, temperature, an): p = config["physiology"]["parameters"] humidity = case["humidity_percent"]/100 if case["formulation"] == "scope": c = temperature-273.15 es = 610.7*10**(7.5*c/(237.3+c)) humidity = min(1, saturation(case["air_C"]+273.15)*humidity/es) return max(p["BallBerryIntercept"], p["BallBerryIntercept"]+ p["BallBerrySlope"]*an*humidity/case["co2"]) def physiology(config, case, temperature, apar): def residual(ci): an = biochemical(config, case, temperature, apar, ci)["An"] gs = stomatal(config, case, temperature, an) diffusion = an if case["formulation"] == "scope" else max(an, 0) return (case["co2"]-ci)*gs/1.6-diffusion lower = .3*case["co2"] if case["formulation"] == "scope" else 10 ci = bisect(residual, lower, case["co2"]+1000) rates = biochemical(config, case, temperature, apar, ci) return dict(**rates, gs=stomatal(config, case, temperature, rates["An"]), Ci=ci) def absorbed_photons(config, case): k = config["physiology"]["constants"] return cavity.shortwave(config, case)*config["wavelength_nm"]*1e-9/( k["planck"]*k["speed_of_light"]*k["avogadro"])*1e6 def balance(config, case, temperature, gs, incoming): constants = config["physiology"]["constants"] ra = case["ra"] rho_cp = config["constants"]["rho_air"]*constants["cp_air"] ta = case["air_C"]+273.15 ea = case["humidity_percent"]/100*saturation(ta) rs = constants["pressure_Pa"]/(8.314*temperature*gs) sensible = rho_cp/ra*(temperature-ta) latent = rho_cp/constants["psychrometric_Pa_K"]*(saturation(temperature)-ea)/(ra+rs) emitted = case["sides"]*case["emissivity"]*config["constants"]["sigma"]*temperature**4 net = cavity.shortwave(config, case)+incoming-emitted return dict(S_abs=cavity.shortwave(config, case), L_abs=incoming, L_emit=emitted, R_n=net, H=sensible, LE=latent, G=0., residual=net-sensible-latent) def reference(config, case): adapted = dict(case, stage="coupled") incoming = cavity.reference(config, adapted) apar = absorbed_photons(config, case) result = [] for element in incoming["elements"]: lw = element["L_abs"] def residual(t): photo = physiology(config, case, t, apar) return balance(config, case, t, photo["gs"], lw)["residual"] t = (case["leaf_temperature_K"] if case["mode"] == "static" else bisect(residual, 280, 330)) photo = physiology(config, case, t, apar) result.append(dict(temperature_K=t, APAR=apar, **photo, fluxes=balance(config, case, t, photo["gs"], lw))) return result def execute(config, case, backend, samples, tolerance=None, initial=None): adapted = dict(case, stage="coupled") start = (case["leaf_temperature_K"] if case["mode"] == "static" else config["initial_temperature_K"] if initial is None else initial) scene = cavity.scene_for(config, adapted, backend, start) profile = less.SpectralProfile([less.SpectralBand(config["wavelength_nm"], config["band_width_nm"], "par")]) rad = scene.simulate(less.RadiationFieldProcess(spectral_profile=profile, n_samples=32, max_depth=8), seed=config["base_seed"]) tolerance = tolerance or config["physiology"]["iteration_tolerances"][-1] product = scene.simulate(less.PhotosynthesisProcess(rad=rad, temperature_mode=case["mode"], n_samples_per_face=samples, max_depth=8, max_iterations=80, convergence_tol=tolerance[0], energy_tolerance=tolerance[1], physiology_tolerance=tolerance[2], damping=.7), microclimate=less.Microclimate(air_temperature=case["air_C"], humidity=case["humidity_percent"], ca=case["co2"], pressure=config["physiology"]["constants"]["pressure_Pa"], oxygen=config["physiology"]["constants"]["oxygen_mmol_mol"], wind_speed=1, downwelling_longwave=config["constants"]["sigma"]* case["environment_K"]**4), seed=config["base_seed"]) mask = np.asarray(product.active_mask, bool) assert mask.sum() == 2 identity = product.metadata["execution"] workflow = "photosynthesis_process_static" if case["mode"] == "static" else "joint_energy_balance" executor = "photosynthesis" if case["mode"] == "static" else "energy_balance" assert identity["native"] and identity["backend"] == backend assert identity["workflow"] == workflow, identity assert identity["executor"] == f"{backend}.slang.{executor}", identity row = dict(backend=backend, case=case["name"], samples=samples, tolerances=tolerance, initial_temperature_K=start, iterations=int(product.iterations), converged=bool(product.converged), execution={k: identity[k] for k in ("backend", "workflow", "executor", "native")}) for name, field in (("temperature_K", "T_leaf"), ("An", "An"), ("gs", "gs"), ("APAR", "APAR")): row[name] = np.asarray(getattr(product, field), float)[mask].tolist() row["sunlit_fraction"] = np.asarray(product.sunlit_fraction, float)[mask].tolist() if case["mode"] == "dynamic": energy = product.metadata["energy_balance_result"] row["fluxes"] = {q: np.asarray(getattr(energy, q), float)[mask].tolist() for q in cavity.QUANTITIES} return row def assess(config, case, row, exact): limits = config["physiology"]["thresholds"] errors = {q: float(np.max(abs(np.array(row[q])-[e[q] for e in exact]))) for q in ("temperature_K", "An", "gs", "APAR")} carbon, stomata, constitutive = [], [], [] for j, (t, an, gs, apar) in enumerate(zip(row["temperature_K"], row["An"], row["gs"], row["APAR"])): transport = an if case["formulation"] == "scope" else max(an, 0) ci = case["co2"]-1.6*transport/gs carbon.append(abs(biochemical(config, case, t, apar, ci)["An"]-an)) stomata.append(abs(stomatal(config, case, t, an)-gs)) if case["mode"] == "dynamic": independent = balance(config, case, t, gs, row["fluxes"]["L_abs"][j]) constitutive.append(abs(independent["residual"])) row["errors"] = errors row["carbon_residual_umol_m2_s"] = max(carbon) row["stomatal_residual_mol_m2_s"] = max(stomata) t_limit = limits["hot_wall_temperature_K"] if "east_wall_K" in case else limits["temperature_K"] passed = (row["converged"] and errors["temperature_K"] <= t_limit and errors["An"] <= limits["An_umol_m2_s"] and errors["gs"] <= limits["gs_mol_m2_s"] and errors["APAR"] <= limits["APAR_umol_m2_s"] and max(carbon) <= limits["An_umol_m2_s"] and max(stomata) <= limits["gs_mol_m2_s"]) if constitutive: f = row["fluxes"] row["closure_W_m2"] = float(np.max(abs(np.array(f["R_n"])-f["H"]-np.array(f["LE"])))) row["constitutive_closure_W_m2"] = max(constitutive) row["flux_errors_W_m2"] = {q: float(np.max(abs(np.array(f[q])- [e["fluxes"][q] for e in exact]))) for q in cavity.QUANTITIES} flux_limit = limits["hot_wall_flux_W_m2"] if "east_wall_K" in case else limits["flux_W_m2"] passed &= (row["closure_W_m2"] <= limits["closure_W_m2"] and max(constitutive) <= limits["constitutive_closure_W_m2"] and max(row["flux_errors_W_m2"].values()) <= flux_limit) row["passed"] = bool(passed) return row def compact(row): return {k: v for k, v in row.items() if k not in ("execution", "fluxes", "sunlit_fraction", "passed")} def run(config, backends=cavity.BACKENDS, selected=None): settings = config["physiology"] cases = [c for c in settings["cases"] if selected is None or c["name"] == selected] refs = {c["name"]: reference(config, c) for c in cases} records, sampling, iteration = [], [], [] for backend in backends: for case in cases: budgets = (settings["hot_wall_samples"] if "east_wall_K" in case else settings["samples"] if case["mode"] == "dynamic" else [32]) levels = [] for samples in budgets: row = assess(config, case, execute(config, case, backend, samples), refs[case["name"]]) levels.append(compact(row)) print(backend, case["name"], samples, row["passed"], row["errors"], flush=True) if len(levels) > 1: change = max(abs(np.array(levels[-1]["temperature_K"])-levels[-2]["temperature_K"])) row["sampling_change_K"] = float(change) row["passed"] = bool(row["passed"] and change <= settings["thresholds"]["sampling_temperature_change_K"]) sampling.append(dict(backend=backend, case=case["name"], levels=levels)) records.append(row) if case["name"] in settings["iteration_cases"] and selected is None: controls = [] for tol in settings["iteration_tolerances"][:-1]: control = execute(config, case, backend, budgets[-1], tol) controls.append(compact(assess(config, case, control, refs[case["name"]]))) controls.append(compact(row)) for initial in config["initial_temperature_controls_K"]: control = execute(config, case, backend, budgets[-1], initial=initial) controls.append(compact(assess(config, case, control, refs[case["name"]]))) changes = {q: float(max(np.max(abs(np.array(c[q])-row[q])) for c in [controls[1]]+controls[3:])) for q in ("temperature_K", "An", "gs")} limits = settings["thresholds"] passed = (all(c["converged"] for c in controls) and changes["temperature_K"] <= limits["iteration_temperature_change_K"] and changes["An"] <= limits["iteration_An_change"] and changes["gs"] <= limits["iteration_gs_change"]) iteration.append(dict(backend=backend, case=case["name"], controls=controls, changes=changes, passed=bool(passed))) return dict(references=refs, records=records, sampling=sampling, iteration=iteration, passed=all(r["passed"] for r in records+iteration)) def figures(config, result, output): import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt colors = ("#2563a6", "#d47713", "#27976b") markers = ("o", "s", "^") plt.rcParams.update({"font.size": 10, "svg.fonttype": "none", "axes.spines.top": False, "axes.spines.right": False}) fig, axes = plt.subplots(2, 2, figsize=(12, 8), layout="constrained") cases = config["physiology"]["cases"] for formulation, style in (("scope", "-"), ("jmax", "--")): names = [c["name"] for c in cases if c["mode"] == "static" and c["formulation"] == formulation and c["co2"] == 400 and c["humidity_percent"] == 50 and c["leaf_temperature_K"] == 298.15] names.sort(key=lambda n: result["references"][n][0]["APAR"]) axes[0, 0].plot([result["references"][n][0]["APAR"] for n in names], [result["references"][n][0]["An"] for n in names], style, color="#334155", label=f"{formulation}: reference") for b, color, marker in zip(cavity.BACKENDS, colors, markers): rows = {r["case"]: r for r in result["records"] if r["backend"] == b} if all(n in rows for n in names): axes[0, 0].plot([np.mean(rows[n]["APAR"]) for n in names], [np.mean(rows[n]["An"]) for n in names], ls="none", marker=marker, markerfacecolor="none", color=color, ms=8, label=b if formulation == "scope" else None) axes[0, 0].set(xlabel="Absorbed PAR / μmol m$^{-2}$ s$^{-1}$", ylabel="Net assimilation / μmol CO$_2$ m$^{-2}$ s$^{-1}$", title="(a) Fixed-temperature light response") axes[0, 0].legend(frameon=False, fontsize=9) names = ["static_strong", "static_high_co2", "static_dry", "static_warm"] x = np.arange(len(names)) axes[0, 1].bar(x, [np.mean([e["gs"] for e in result["references"][n]]) for n in names], color="#9dc7d6", label="Reference") for b, color, marker in zip(cavity.BACKENDS, colors, markers): rows = {r["case"]: r for r in result["records"] if r["backend"] == b} if not rows: continue axes[0, 1].plot(x, [np.mean(rows[n]["gs"]) for n in names], ls="none", marker=marker, markerfacecolor="none", ms=8, color=color, label=b) dynamic = [c for c in cases if c["mode"] == "dynamic"] indices = np.arange(len(dynamic)) axes[1, 0].plot(indices, [np.mean(rows[c["name"]]["temperature_K"]) for c in dynamic], ls="none", marker=marker, markerfacecolor="none", ms=8, color=color) control = next(r for r in result["iteration"] if r["backend"] == b and r["case"] == "dynamic_scope") axes[1, 1].semilogx([r["tolerances"][0] for r in control["controls"][:3]], [r["constitutive_closure_W_m2"] for r in control["controls"][:3]], marker=marker, markerfacecolor="none", color=color, label=b) axes[0, 1].set(xticks=x, xticklabels=["Baseline", "800 ppm", "RH 30%", "35°C leaf"], ylabel="Stomatal conductance / mol H$_2$O m$^{-2}$ s$^{-1}$", title="(b) Fixed-temperature stomatal response") axes[0, 1].legend(frameon=False, fontsize=9) axes[1, 0].plot(indices, [np.mean([e["temperature_K"] for e in result["references"][c["name"]]]) for c in dynamic], color="#334155", zorder=1, label="Reference") axes[1, 0].set(xticks=indices, xticklabels=["Baseline", "RH 90%", "800 ppm", "Hot wall", "Jmax"], ylabel="Leaf temperature / K", title="(c) Coupled leaf temperature") axes[1, 1].invert_xaxis() axes[1, 1].set(xlabel="Temperature iteration tolerance / K (tighter →)", ylabel="Independent energy residual / W m$^{-2}$", title="(d) Coupled iteration convergence") for ax in axes.flat: ax.grid(alpha=.18, axis="y") fig.savefig(output/"physiology.svg") fig.savefig(output/"physiology.png", dpi=150) plt.close(fig) if __name__ == "__main__": import argparse parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--backend", choices=cavity.BACKENDS) parser.add_argument("--case") parser.add_argument("--output", type=Path, default=Path("build/leaf-physiology-validation")) args = parser.parse_args() config = json.loads((cavity.DATA/"input.json").read_text(encoding="utf-8")) result = run(config, [args.backend] if args.backend else cavity.BACKENDS, args.case) args.output.mkdir(parents=True, exist_ok=True) (args.output/"results.json").write_text(json.dumps(result, indent=2, allow_nan=False), encoding="utf-8") raise SystemExit(0 if result["passed"] else 1)