Helmholtz equation

Let us consider the Helmholtz equation in two space dimensions

\[\begin{aligned} &\Delta u(x, y)+k^{2} u(x, y)=q(x, y), \quad(x, y) \in \Omega:=(-1,1)^2 \\ &u(x, y)=0, \quad(x, y) \in \partial \Omega \end{aligned}\]

where

\[q(x, y)=-\left(a_{1} \pi\right)^{2} \sin \left(a_{1} \pi x\right) \sin \left(a_{2} \pi y\right)-\left(a_{2} \pi\right)^{2} \sin \left(a_{1} \pi x\right) \sin \left(a_{2} \pi y\right)+k^{2} \sin \left(a_{1} \pi x\right) \sin \left(a_{2} \pi y\right).\]

The exact solution is $u(x,y)=\sin{a_1\pi x}\sin{a_2\pi y}$. We chose $k=1, a_1 = 1$ and $a_2 = 4$.

using ModelingToolkit, IntervalSets, Sophon, Lux, Zygote
using Optimization, OptimizationOptimJL

@parameters x,y
@variables u(..)
Dxx = Differential(x)^2
Dyy = Differential(y)^2

a1 = 1
a2 = 4
k = 1

q(x,y) = -(a1*π)^2 * sin(a1*π*x) * sin(a2*π*y) - (a2*π)^2 * sin(a1*π*x) * sin(a2*π*y) + k^2 * sin(a1*π*x) * sin(a2*π*y)
eq = Dxx(u(x,y)) + Dyy(u(x,y)) + k^2 * u(x,y) ~ q(x,y)
domains = [x ∈ Interval(-1,1), y ∈ Interval(-1,1)]
bcs = [u(-1,y) ~ 0, u(1,y) ~ 0, u(x, -1) ~ 0, u(x, 1) ~ 0]

@named helmholtz = PDESystem(eq, bcs, domains, [x,y], [u(x,y)])

\[ \begin{align} u\left( x, y \right) + \frac{\mathrm{d}}{\mathrm{d}y} \frac{\mathrm{d}}{\mathrm{d}y} u\left( x, y \right) + \frac{\mathrm{d}}{\mathrm{d}x} \frac{\mathrm{d}}{\mathrm{d}x} u\left( x, y \right) =& - 166.78 \sin\left( 12.566 y \right) \sin\left( 3.1416 x \right) \end{align} \]

Note that the boundary conditions are compatible with periocity, which allows us to apply BACON.

chain = BACON(2, 1, 5, 2; hidden_dims = 32, num_layers=5)
pinn = PINN(chain)
sampler = QuasiRandomSampler(300, 100)
strategy = NonAdaptiveTraining()

prob = Sophon.discretize(helmholtz, pinn, sampler, strategy)

@showprogress res = Optimization.solve(prob, BFGS(); maxiters=1000)

retcode: Success
u: ComponentVector{Float64}(filters = (filter_1 = (bias = [0.028183566753781385; -0.8192489481759255; … ; -0.8494417979599271; -1.1061428442021761;;]), filter_2 = (bias = [0.6736395572998042; 0.5613293022149991; … ; 0.4807571326668061; -0.9425557192172175;;]), filter_3 = (bias = [0.6791908548814871; 0.5826621158507546; … ; -0.0918765815912823; 0.767372992458497;;]), filter_4 = (bias = [-0.4535561162129571; -0.44739000853926464; … ; -0.7630260759887113; -0.3723271287313381;;]), filter_5 = (bias = [-0.5031143011140751; -0.5135322833506123; … ; -0.9181285523506575; 0.7564907939645017;;])), linear_layers = (layer_1 = (weight = [0.48861023840784784 -0.5380161894803988 … 0.2567403138111084 0.20953304253609267; -0.22737835944731544 -0.4650604735393933 … 0.0968382229526869 0.04151039487364857; … ; -0.3749527208560011 -0.21263689375197928 … -0.11736037418556353 0.3454034458681847; 0.3564695195010524 0.39708084042797537 … -0.20870191297960036 0.3934771801075347], bias = [0.025445261911640284; -0.018267029279022025; … ; -0.03000331261632699; 0.06112721146255124;;]), layer_2 = (weight = [-0.2532777274807004 -0.009778659996892786 … -0.16349484002276718 -0.2014155811970701; -0.48381217855786285 -0.34547435178515234 … -0.19263709207856886 0.08308867125353095; … ; 0.40298662318180145 0.3430510760363498 … -0.40962355980324094 0.14592932519091273; -0.5041481369078836 0.28636204592353437 … 0.059877887505800434 -0.2803438351996603], bias = [-0.03082431736895129; 0.03733137050579801; … ; 0.004474729545184752; 0.02988366860538561;;]), layer_3 = (weight = [-0.11738853071409813 0.19632401206799394 … -0.42065785911431625 -0.14949864869050938; 0.1630065656446501 -0.36188399865106524 … -0.4300766018399455 -0.11648240729664733; … ; -0.056388865391109454 -0.10490497752644401 … -0.20612849959928387 -0.24921880210508537; 0.017579642074646223 -0.32159389283947476 … -0.14384152680458262 0.17194359571875803], bias = [-0.02303272483436167; 0.002113800949979219; … ; -0.03873002027727485; -0.0020053347473351704;;]), layer_4 = (weight = [0.3368519589282945 -0.12049395416113537 … 0.2623515068413288 -0.05871443936889818; 0.12506121900367265 -0.02776589425119601 … -0.04769962029885092 -0.08269816109785966; … ; -0.2447386493940465 0.09645377564675675 … 0.19722380297040173 -0.019946467261616998; 0.036475774004251164 -0.38427377827698744 … -0.37528686709234343 -0.0014429434416557496], bias = [-0.006146627276299785; 0.007873676124658937; … ; 0.030227268246602122; 0.0004091599826215488;;])), output_layer = (weight = [-0.1263910112572273 0.2284525003085034 … 0.5634497792091683 -0.10733235685639825], bias = [0.031948711990690895;;]))

Let's plot the result.

phi = pinn.phi
ps = res.u
xs, ys= [infimum(d.domain):0.01:supremum(d.domain) for d in domains]
u_analytic(x,y) = sinpi(a1*x)*sinpi(a2*y)
u_real = [u_analytic(x,y) for x in xs, y in ys]

u_pred = [sum(phi(([x,y]), ps)) for x in xs, y in ys]

using CairoMakie
axis = (xlabel="x", ylabel="y", title="Analytical Solution")
fig, ax1, hm1 = heatmap(xs, ys, u_real, axis=axis)
Colorbar(fig[:, end+1], hm1)
ax2, hm2= heatmap(fig[1, end+1], xs, ys, u_pred, axis= merge(axis, (;title = "Prediction")))
Colorbar(fig[:, end+1], hm2)
ax3, hm3 = heatmap(fig[1, end+1], xs, ys, abs.(u_pred-u_real), axis= merge(axis, (;title = "Absolute Error")))
Colorbar(fig[:, end+1], hm3)
fig