1- Hamiltonian-structure- preserving implementation of the Benjamin-Bona-Mahoney equation
2- =====================================================================================
1+ Mixed Hamiltonian-preserving formulation of the Benjamin-Bona-Mahoney equation
2+ ==============================================================================
33
44This demo solves the Benjamin-Bona-Mahony equation:
55
@@ -20,18 +20,18 @@ BBM is known to have a Hamiltonian structure, and there are several canonical po
2020 I_3 & = \int \frac {u^2 }{2 } + \frac {u^3 }{6 } \, \mathrm {d}x
2121
2222 The BBM invariants are the total momentum :math: `I_1 `, the :math: `H^1 `-energy
23- norm :math: `I_2 `, and the Hamiltonian :math: `I_3 `.
24- The Hamiltonian variational formulation reads
23+ :math: `I_2 `, and the Hamiltonian :math: `I_3 `.
24+ The Hamiltonian mixed variational formulation reads
2525
2626.. math ::
2727
28- (\partial _t u + \partial _x \tilde {w}_H, v )_{H^ 1 } & = 0
28+ (\partial _t u, v)_{H^ 1 } - ( \tilde {w}_H, \partial _x v_H )_{L^ 2 } & = 0
2929
30- (\tilde {w}_H, v_H)_{H^ 1 } & = \langle \frac {\delta I_3 }{\delta u}, v_H \rangle
30+ (\tilde {w}_H, v_H)_{L^ 2 } & = \langle \frac {\delta I_3 }{\delta u}, v_H \rangle
3131
3232 For all test functions :math: `v, v_H` in a suitable function space.
3333The numerical scheme in this demo introduces
34- the :math: `H^ 1 `-Riesz representative of the Fréchet derivative of the
34+ the :math: `L^ 2 `-Riesz representative of the Fréchet derivative of the
3535Hamiltonian :math: `\frac {\delta I_3 }{\delta u}`
3636as the auxiliary variable :math: `\tilde {w}_H`.
3737
@@ -52,8 +52,8 @@ Firedrake, Irksome, and other imports::
5252
5353 from firedrake import (Constant, Function, FunctionSpace,
5454 PeriodicIntervalMesh, SpatialCoordinate, TestFunction, TrialFunction,
55- assemble, derivative, dx, errornorm, exp, grad, inner,
56- interpolate, norm, plot, project, replace, solve, split
55+ assemble, derivative, dx, exp, grad, inner,
56+ norm, plot, project, replace, solve, split
5757 )
5858
5959 from irksome import Dt, GalerkinTimeStepper, TimeQuadratureLabel
@@ -88,11 +88,11 @@ Next, we define the domain and the exact solution ::
8888This sets up the function space for the unknown :math: `u` and
8989auxiliary variable :math: `\tilde {w}_H`::
9090
91- space_deg = 3
91+ space_deg = 1
9292 time_deg = 1
9393
94- V = FunctionSpace(msh, "Hermite ", space_deg)
95- Z = V * V
94+ V = FunctionSpace(msh, "CG ", space_deg)
95+ Z = V* V
9696
9797We next define the BBM invariants. Again, the discrete formulation preserves
9898:math: `I_1 ` and :math: `I_3 ` up to solver tolerances and roundoff errors,
@@ -119,27 +119,28 @@ initial condition for the auxiliary variable. We need to find :math:`\tilde{w}_
119119
120120 ::
121121
122- uwH = Function(Z)
123- u0, wH0 = uwH .subfunctions
122+ uw = Function(Z)
123+ u0, w0 = uw .subfunctions
124124
125125 v = TestFunction(V)
126126 w = TrialFunction(V)
127- a = h1inner(w, v) * dx
127+ a0 = inner(w, v) * dx
128+ a1 = h1inner(w, v) * dx
128129 dHdu = derivative(I3(u0), u0, v)
129130
130- solve(a == h1inner(uexact, v)*dx, u0)
131- solve(a == dHdu, wH0 )
131+ solve(a1 == h1inner(uexact, v)*dx, u0)
132+ solve(a0 == dHdu, w0 )
132133
133134Visualize the initial condition::
134135
135136 fig, axes = plt.subplots(1)
136- plot(Function(FunctionSpace(msh, "CG", 1)).interpolate(u0) , axes=axes)
137+ plot(u0 , axes=axes)
137138 axes.set_title("Initial condition")
138139 axes.set_xlabel("x")
139140 axes.set_ylabel("u")
140- plt.savefig("bbm_init .png")
141+ plt.savefig("bbm_aux_init .png")
141142
142- .. figure :: bbm_init .png
143+ .. figure :: bbm_aux_init .png
143144 :align: center
144145
145146Create time quadrature labels::
@@ -153,10 +154,9 @@ Create time quadrature labels::
153154This tags several of the terms with a low-order time integration scheme,
154155but forces a higher-order method on the nonlinear term::
155156
156- u, wH = split(uwH )
157+ u, w = split(uw )
157158 v, vH = split(TestFunction(Z))
158-
159- Flow = h1inner(Dt(u) + wH.dx(0), v) * dx + h1inner(wH, vH) * dx
159+ Flow = h1inner(Dt(u), v) * dx + inner(w.dx(0), v)*dx + inner(w, vH)*dx
160160 Fhigh = replace(dHdu, {u0: u})
161161
162162 F = Llow(Flow) - Lhigh(Fhigh(vH))
@@ -165,8 +165,7 @@ This sets up the cPG time stepper. There are two fields in the unknown, we
165165indicate the second one is an auxiliary and hence to be discretized in the DG
166166test space instead by passing the `aux_indices ` keyword::
167167
168- stepper = GalerkinTimeStepper(
169- F, time_deg, t, dt, uwH, aux_indices=[1])
168+ stepper = GalerkinTimeStepper(F, time_deg, t, dt, uw, aux_indices=[1])
170169
171170UFL expressions for the invariants, which we are going to track as we go
172171through time steps::
@@ -201,7 +200,7 @@ Visualize invariant preservation::
201200 axes.set_xlabel("Time")
202201 axes.set_ylabel("I(t)")
203202 axes.legend()
204- plt.savefig("invariants .png")
203+ plt.savefig("bbm_aux_invariants .png")
205204 axes.clear()
206205
207206 for i in (0, 1, 2):
@@ -210,22 +209,22 @@ Visualize invariant preservation::
210209 axes.set_xlabel("Time")
211210 axes.set_ylabel("|1-I/I(0)|")
212211 axes.legend()
213- plt.savefig("invariant_errors .png")
212+ plt.savefig("bbm_aux_errors .png")
214213
215- .. figure :: invariants .png
214+ .. figure :: bbm_aux_invariants .png
216215 :align: center
217216
218- .. figure :: invariant_errors .png
217+ .. figure :: bbm_aux_errors .png
219218 :align: center
220219
221220Visualize the solution at final time step::
222221
223222 axes.clear()
224- plot(Function(FunctionSpace(msh, "CG", 1)).interpolate(u0) , axes=axes)
223+ plot(u0 , axes=axes)
225224 axes.set_title(f"Solution at time {tfinal}")
226225 axes.set_xlabel("x")
227226 axes.set_ylabel("u")
228- plt.savefig("bbm_final .png")
227+ plt.savefig("bbm_aux_final .png")
229228
230- .. figure :: bbm_final .png
229+ .. figure :: bbm_aux_final .png
231230 :align: center
0 commit comments