--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/experiments/laser_and_mirrors_aux.py Thu Feb 26 09:32:12 2026 -0500
@@ -0,0 +1,196 @@
+#
+# Laser and mirrors convection-diffusion experiment with unknown diffusivity and convectivity;
+# Includes `generic_setup` for the known case as well.
+#
+
+import numpy as np
+from dolfinx import fem
+from measures import DiscreteMeasure_2_f64
+from numpy.linalg import norm
+from pointsource_pde import Problem, RegTerm
+from pointsource_pde.compose import (
+ ComposeFnWithOperator,
+ SumOfSeparableFunctions,
+)
+from pointsource_pde.convection_diffusion import (
+ ConvectionDiffusion,
+ PlotFactory,
+ QuadraticRegularisation,
+ XBound,
+)
+from pointsource_pde.laser_sampling import LaserSampling
+from pointsource_pde.quadratic_dataterm import QuadraticDataTerm
+
+# Give name to the problem
+name = "laser_and_mirrors_aux"
+
+# Fine-tune algorithms
+alg_finetune = {
+ "tau0": 0.99,
+ "theta0": 0.01,
+ "inner_method": "PP",
+ "inner_pdps_τσ0": (19.8, 0.05),
+ "inner_pp_τ": [1000.0, 1000.0],
+ "merge": True,
+ "fitness_merging": True,
+}
+
+algorithm_overrides = {
+ "RadonFB": alg_finetune,
+ "RadonSlidingFB": alg_finetune,
+ "RadonPDPS": alg_finetune,
+ "RadonSlidingPDPS": alg_finetune,
+}
+
+
+# Setup shared between experiment with known/unknown conductivity and diffusivity
+def generic_setup(prefix, simple=True):
+
+ # Create true source
+ μ_true = DiscreteMeasure_2_f64(
+ [([0.1, 0.3], 0.08), ([0.4, 0.25], 0.05), ([0.25, 0.13], 0.06)]
+ )
+
+ # Create convection-diffusion PDE
+ pde = ConvectionDiffusion(
+ u0=lambda x: 0.0 * x[1],
+ g=lambda x: 0.0 * x[1],
+ domain_size=[0, 0.5, 0, 0.5],
+ num_steps=50,
+ )
+
+ # Set up lasers and mirrors
+ lasers = [
+ np.array([0.1, 0.1]),
+ np.array([0.1, 0.4]),
+ np.array([0.4, 0.1]),
+ np.array([0.4, 0.1]),
+ ]
+ mirrors_per_side = 10
+ x0, x1, y0, y1 = (
+ pde.domain_size[0],
+ pde.domain_size[1],
+ pde.domain_size[2],
+ pde.domain_size[3],
+ )
+ beams = [
+ (src, np.array(dst))
+ for src in lasers
+ for i in range(0, mirrors_per_side)
+ for z in [(i + 1) / (mirrors_per_side + 1)]
+ for dst in [
+ [z * (x1 - x0) + x0, y0],
+ [z * (x1 - x0) + x0, y1],
+ [x0, z * (y1 - y0) + y0],
+ [x1, z * (y1 - y0) + y0],
+ ]
+ ]
+ L_i = LaserSampling(pde.V, pde.domain, pde.nx, pde.ny, beams)
+
+ v0 = 1.0
+
+ if simple:
+ # Scalar convectivity and diffusivity
+ r = 0.5
+ θ = 30 * np.pi / 180
+ c1 = r * np.cos(θ)
+ c2 = r * np.sin(θ)
+ k = 0.01 # diffusive
+ # k = 0.001 # high wind
+ else:
+ # Convectivity and diffusivity fields
+ def velocity_x(x):
+ r = np.sqrt(x[0] ** 2 + x[1] ** 2 + 1e-10)
+ return v0 * x[0] / r
+
+ def velocity_y(x):
+ r = np.sqrt(x[0] ** 2 + x[1] ** 2 + 1e-10)
+ return v0 * x[1] / r
+
+ c1 = fem.Function(pde.V)
+ c1.interpolate(velocity_x)
+ c2 = fem.Function(pde.V)
+ c2.interpolate(velocity_y)
+ k = fem.Function(pde.V)
+ k0 = 1.0
+ k.interpolate(lambda x: k0 + 0.2 * np.sin(5 * x[0]) * np.cos(5 * x[1])) # k is
+
+ auxtrue = (k, (c1, c2))
+
+ # Simulate measurements
+ u_n_list = pde.apply((μ_true, (k, (c1, c2))))
+ b_true = [L_i.apply(u) for u in u_n_list]
+ b = [b_i + L_i.noise(0.01 * norm(b_i) / np.sqrt(np.size(b_i))) for b_i in b_true]
+
+ # Estimate bounds on domains and values
+ dat_bound = 10 * max([norm(b_i) for b_i in b])
+ μ_bound = 1.3 * sum([alpha for _point, alpha in μ_true.iter_padded()])
+ xbound = XBound(mu_dual=μ_bound, k_min=k, k_max=k, c_Linf=r)
+
+ # Create data term
+ dat = ComposeFnWithOperator(
+ SumOfSeparableFunctions(
+ [
+ QuadraticDataTerm(
+ L_i,
+ data,
+ xbound=dat_bound,
+ b_norm=norm(b),
+ λ=pde.dt / len(beams),
+ )
+ for data in b
+ ]
+ ),
+ pde,
+ xbound=[xbound],
+ xbound_pair=[xbound],
+ )
+
+ # Plot true and initial data
+ plot_factory = PlotFactory(pde)
+ plotter = plot_factory.produce(prefix)
+ ω, _ = dat.diff((DiscreteMeasure_2_f64([]), auxtrue))
+ ω_true, _ = dat.diff((μ_true, auxtrue))
+ np.savez_compressed(
+ "%s/omega_0.npz" % prefix,
+ true_mu=[[point[0], point[1], alpha] for point, alpha in μ_true.iter_padded()],
+ omega0=ω.x.array,
+ true_omega=ω_true.x.array,
+ true_u_n_list_array=[u_n_list[i].x.array for i in plotter.plot_frames],
+ )
+
+ return dat, auxtrue, μ_bound, μ_true, plot_factory
+
+
+def setup(prefix):
+ dat, auxtrue, _μ_bound, μ_true, plot_factory = generic_setup(prefix)
+
+ reg = RegTerm.NonnegRadon(10e-7)
+
+ μ0 = DiscreteMeasure_2_f64([])
+ aux0 = auxtrue
+ aux = QuadraticRegularisation(0.001, aux0)
+ ax = aux.apply(aux0)
+ inival = dat.apply((μ0, aux0))
+ print("Initial value:", inival + ax, ax)
+ print("Value at true μ:", dat.apply((μ_true, auxtrue)) + ax)
+
+ # auxinit = np.c_[
+ # np.zeros_like(k),
+ # np.zeros_like(c1),
+ # np.zeros_like(c2),
+ # ]
+ # print(np.shape(auxinit))
+
+ # We need Θ² such that ⟨F'(μ+Δ, aux-F'(μ, aux)|Δ⟩ ≤ Θ²|γ|(c_2), where Δ=(π_♯^1-π_♯^0)γ.
+ # We have
+ #
+ # ⟨F'(μ+Δ, aux)-F'(μ, aux)|Δ⟩ ≤ ‖⟨F'(μ+Δ, aux)-F'(μ, aux)‖ ‖Δ‖
+ # ≤ L_{F'(·, aux)} ‖Δ‖²,
+ #
+ # so Θ² ≥ L_{F'(·, aux)} ‖Δ‖² / |γ|(c_2). Thsi doesn't give anything useful if the
+ # transport is very small in distance.
+
+ dat.curvature_bound_components = lambda: (None, None) # 1.0
+
+ return Problem(dat, reg, aux, aux0, μ0, plot_factory=plot_factory)