Mercurial > repos > pointsource_pde / file comparison

comparison: experiments/laser_and_mirrors_aux.py

experiments/laser_and_mirrors_aux.py

changeset 1
a4137aedcb3a
child 3
c3a4f4bb87f7
equal deleted inserted replaced
0:7ec1cfe19a24 1:a4137aedcb3a
1 #
2 # Laser and mirrors convection-diffusion experiment with unknown diffusivity and convectivity;
3 # Includes `generic_setup` for the known case as well.
4 #
5
6 import numpy as np
7 from dolfinx import fem
8 from measures import DiscreteMeasure_2_f64
9 from numpy.linalg import norm
10 from pointsource_pde import Problem, RegTerm
11 from pointsource_pde.compose import (
12 ComposeFnWithOperator,
13 SumOfSeparableFunctions,
14 )
15 from pointsource_pde.convection_diffusion import (
16 ConvectionDiffusion,
17 PlotFactory,
18 QuadraticRegularisation,
19 XBound,
20 )
21 from pointsource_pde.laser_sampling import LaserSampling
22 from pointsource_pde.quadratic_dataterm import QuadraticDataTerm
23
24 # Give name to the problem
25 name = "laser_and_mirrors_aux"
26
27 # Fine-tune algorithms
28 alg_finetune = {
29 "tau0": 0.99,
30 "theta0": 0.01,
31 "inner_method": "PP",
32 "inner_pdps_τσ0": (19.8, 0.05),
33 "inner_pp_τ": [1000.0, 1000.0],
34 "merge": True,
35 "fitness_merging": True,
36 }
37
38 algorithm_overrides = {
39 "RadonFB": alg_finetune,
40 "RadonSlidingFB": alg_finetune,
41 "RadonPDPS": alg_finetune,
42 "RadonSlidingPDPS": alg_finetune,
43 }
44
45
46 # Setup shared between experiment with known/unknown conductivity and diffusivity
47 def generic_setup(prefix, simple=True):
48
49 # Create true source
50 μ_true = DiscreteMeasure_2_f64(
51 [([0.1, 0.3], 0.08), ([0.4, 0.25], 0.05), ([0.25, 0.13], 0.06)]
52 )
53
54 # Create convection-diffusion PDE
55 pde = ConvectionDiffusion(
56 u0=lambda x: 0.0 * x[1],
57 g=lambda x: 0.0 * x[1],
58 domain_size=[0, 0.5, 0, 0.5],
59 num_steps=50,
60 )
61
62 # Set up lasers and mirrors
63 lasers = [
64 np.array([0.1, 0.1]),
65 np.array([0.1, 0.4]),
66 np.array([0.4, 0.1]),
67 np.array([0.4, 0.1]),
68 ]
69 mirrors_per_side = 10
70 x0, x1, y0, y1 = (
71 pde.domain_size[0],
72 pde.domain_size[1],
73 pde.domain_size[2],
74 pde.domain_size[3],
75 )
76 beams = [
77 (src, np.array(dst))
78 for src in lasers
79 for i in range(0, mirrors_per_side)
80 for z in [(i + 1) / (mirrors_per_side + 1)]
81 for dst in [
82 [z * (x1 - x0) + x0, y0],
83 [z * (x1 - x0) + x0, y1],
84 [x0, z * (y1 - y0) + y0],
85 [x1, z * (y1 - y0) + y0],
86 ]
87 ]
88 L_i = LaserSampling(pde.V, pde.domain, pde.nx, pde.ny, beams)
89
90 v0 = 1.0
91
92 if simple:
93 # Scalar convectivity and diffusivity
94 r = 0.5
95 θ = 30 * np.pi / 180
96 c1 = r * np.cos(θ)
97 c2 = r * np.sin(θ)
98 k = 0.01 # diffusive
99 # k = 0.001 # high wind
100 else:
101 # Convectivity and diffusivity fields
102 def velocity_x(x):
103 r = np.sqrt(x[0] ** 2 + x[1] ** 2 + 1e-10)
104 return v0 * x[0] / r
105
106 def velocity_y(x):
107 r = np.sqrt(x[0] ** 2 + x[1] ** 2 + 1e-10)
108 return v0 * x[1] / r
109
110 c1 = fem.Function(pde.V)
111 c1.interpolate(velocity_x)
112 c2 = fem.Function(pde.V)
113 c2.interpolate(velocity_y)
114 k = fem.Function(pde.V)
115 k0 = 1.0
116 k.interpolate(lambda x: k0 + 0.2 * np.sin(5 * x[0]) * np.cos(5 * x[1])) # k is
117
118 auxtrue = (k, (c1, c2))
119
120 # Simulate measurements
121 u_n_list = pde.apply((μ_true, (k, (c1, c2))))
122 b_true = [L_i.apply(u) for u in u_n_list]
123 b = [b_i + L_i.noise(0.01 * norm(b_i) / np.sqrt(np.size(b_i))) for b_i in b_true]
124
125 # Estimate bounds on domains and values
126 dat_bound = 10 * max([norm(b_i) for b_i in b])
127 μ_bound = 1.3 * sum([alpha for _point, alpha in μ_true.iter_padded()])
128 xbound = XBound(mu_dual=μ_bound, k_min=k, k_max=k, c_Linf=r)
129
130 # Create data term
131 dat = ComposeFnWithOperator(
132 SumOfSeparableFunctions(
133 [
134 QuadraticDataTerm(
135 L_i,
136 data,
137 xbound=dat_bound,
138 b_norm=norm(b),
139 λ=pde.dt / len(beams),
140 )
141 for data in b
142 ]
143 ),
144 pde,
145 xbound=[xbound],
146 xbound_pair=[xbound],
147 )
148
149 # Plot true and initial data
150 plot_factory = PlotFactory(pde)
151 plotter = plot_factory.produce(prefix)
152 ω, _ = dat.diff((DiscreteMeasure_2_f64([]), auxtrue))
153 ω_true, _ = dat.diff((μ_true, auxtrue))
154 np.savez_compressed(
155 "%s/omega_0.npz" % prefix,
156 true_mu=[[point[0], point[1], alpha] for point, alpha in μ_true.iter_padded()],
157 omega0=ω.x.array,
158 true_omega=ω_true.x.array,
159 true_u_n_list_array=[u_n_list[i].x.array for i in plotter.plot_frames],
160 )
161
162 return dat, auxtrue, μ_bound, μ_true, plot_factory
163
164
165 def setup(prefix):
166 dat, auxtrue, _μ_bound, μ_true, plot_factory = generic_setup(prefix)
167
168 reg = RegTerm.NonnegRadon(10e-7)
169
170 μ0 = DiscreteMeasure_2_f64([])
171 aux0 = auxtrue
172 aux = QuadraticRegularisation(0.001, aux0)
173 ax = aux.apply(aux0)
174 inival = dat.apply((μ0, aux0))
175 print("Initial value:", inival + ax, ax)
176 print("Value at true μ:", dat.apply((μ_true, auxtrue)) + ax)
177
178 # auxinit = np.c_[
179 # np.zeros_like(k),
180 # np.zeros_like(c1),
181 # np.zeros_like(c2),
182 # ]
183 # print(np.shape(auxinit))
184
185 # We need Θ² such that ⟨F'(μ+Δ, aux-F'(μ, aux)|Δ⟩ ≤ Θ²|γ|(c_2), where Δ=(π_♯^1-π_♯^0)γ.
186 # We have
187 #
188 # ⟨F'(μ+Δ, aux)-F'(μ, aux)|Δ⟩ ≤ ‖⟨F'(μ+Δ, aux)-F'(μ, aux)‖ ‖Δ‖
189 # ≤ L_{F'(·, aux)} ‖Δ‖²,
190 #
191 # so Θ² ≥ L_{F'(·, aux)} ‖Δ‖² / |γ|(c_2). Thsi doesn't give anything useful if the
192 # transport is very small in distance.
193
194 dat.curvature_bound_components = lambda: (None, None) # 1.0
195
196 return Problem(dat, reg, aux, aux0, μ0, plot_factory=plot_factory)

Mercurial Repository: pointsource_pde

mercurial