Skip to content

Commit

Permalink
Expand MA Helmholtz demo to interpolate soln field (#145)
Browse files Browse the repository at this point in the history
Closes #27.
  • Loading branch information
ddundo authored Jan 9, 2025
1 parent 1ffa472 commit f8ac7b9
Showing 1 changed file with 227 additions and 6 deletions.
233 changes: 227 additions & 6 deletions demos/monge_ampere_helmholtz.py
Original file line number Diff line number Diff line change
Expand Up @@ -112,7 +112,7 @@ def solve_helmholtz(mesh):
alpha = Constant(5.0)


def monitor(mesh):
def monitor_exact(mesh):
V = FunctionSpace(mesh, "CG", 1)
Hnorm = Function(V, name="Hnorm")
H = grad(grad(u_exact(mesh)))
Expand All @@ -122,11 +122,11 @@ def monitor(mesh):
return m


# Plot the monitor function on the original mesh
# Plot the monitor function on the original mesh:

fig, axes = plt.subplots()
m = Function(u_h, name="monitor")
m.interpolate(monitor(mesh))
m.interpolate(monitor_exact(mesh))
contours = tripcolor(m, axes=axes)
fig.colorbar(contours)
plt.savefig("monge_ampere_helmholtz-monitor.jpg")
Expand All @@ -135,10 +135,19 @@ def monitor(mesh):
# .. figure:: monge_ampere_helmholtz-monitor.jpg
# :figwidth: 60%
# :align: center
#
# Now we can construct a :class:`~movement.monge_ampere.MongeAmpereMover` instance to
# adapt the mesh based on this monitor function. We will also time how long the mesh
# movement step takes to complete [#f1]_.

import time

rtol = 1.0e-08
mover = MongeAmpereMover(mesh, monitor, method="quasi_newton", rtol=rtol)
mover = MongeAmpereMover(mesh, monitor_exact, method="quasi_newton", rtol=rtol)

t0 = time.time()
mover.move()
print(f"Time taken: {time.time() - t0:.2f} seconds")

# For every iteration the MongeAmpereMover prints the minimum to maximum ratio of
# the cell areas in the mesh, the residual in the Monge-Ampère equation, and the
Expand All @@ -165,6 +174,7 @@ def monitor(mesh):
# 16 Volume ratio 1.15 Variation (σ/μ) 2.57e-02 Residual 1.94e-07
# 17 Volume ratio 1.15 Variation (σ/μ) 2.57e-02 Residual 5.86e-08
# Solver converged in 17 iterations.
# Time taken: 2.52 seconds
#
# Plotting the resulting mesh:

Expand Down Expand Up @@ -210,7 +220,7 @@ def monitor(mesh):
from animate import RiemannianMetric


def monitor2(mesh):
def monitor_solve(mesh):
u_h = solve_helmholtz(mesh)
V = FunctionSpace(mesh, "CG", 1)
TV = TensorFunctionSpace(mesh, "CG", 1)
Expand All @@ -224,8 +234,45 @@ def monitor2(mesh):
return m


mover = MongeAmpereMover(mesh, monitor2, method="quasi_newton", rtol=rtol)
mover = MongeAmpereMover(mesh, monitor_solve, method="quasi_newton", rtol=rtol)

t0 = time.time()
mover.move()
print(f"Time taken: {time.time() - t0:.2f} seconds")

# .. code-block:: none
#
# 0 Volume ratio 5.04 Variation (σ/μ) 4.90e-01 Residual 4.93e-01
# 1 Volume ratio 2.60 Variation (σ/μ) 2.27e-01 Residual 2.24e-01
# 2 Volume ratio 1.63 Variation (σ/μ) 1.09e-01 Residual 1.03e-01
# 3 Volume ratio 1.31 Variation (σ/μ) 5.30e-02 Residual 4.33e-02
# 4 Volume ratio 1.20 Variation (σ/μ) 3.31e-02 Residual 1.96e-02
# 5 Volume ratio 1.16 Variation (σ/μ) 2.65e-02 Residual 8.77e-03
# 6 Volume ratio 1.14 Variation (σ/μ) 2.46e-02 Residual 3.83e-03
# 7 Volume ratio 1.14 Variation (σ/μ) 2.40e-02 Residual 1.63e-03
# 8 Volume ratio 1.14 Variation (σ/μ) 2.38e-02 Residual 6.77e-04
# 9 Volume ratio 1.14 Variation (σ/μ) 2.37e-02 Residual 2.72e-04
# 10 Volume ratio 1.14 Variation (σ/μ) 2.37e-02 Residual 1.06e-04
# 11 Volume ratio 1.14 Variation (σ/μ) 2.37e-02 Residual 3.97e-05
# 12 Volume ratio 1.14 Variation (σ/μ) 2.37e-02 Residual 1.44e-05
# 13 Volume ratio 1.14 Variation (σ/μ) 2.37e-02 Residual 5.30e-06
# 14 Volume ratio 1.14 Variation (σ/μ) 2.37e-02 Residual 2.20e-06
# 15 Volume ratio 1.14 Variation (σ/μ) 2.37e-02 Residual 1.11e-06
# 16 Volume ratio 1.14 Variation (σ/μ) 2.37e-02 Residual 6.13e-07
# 17 Volume ratio 1.14 Variation (σ/μ) 2.37e-02 Residual 3.38e-07
# 18 Volume ratio 1.14 Variation (σ/μ) 2.37e-02 Residual 1.80e-07
# 19 Volume ratio 1.14 Variation (σ/μ) 2.37e-02 Residual 9.21e-08
# Solver converged in 19 iterations.
# Time taken: 6.17 seconds

fig, axes = plt.subplots()
triplot(mover.mesh, axes=axes)
axes.set_aspect(1)
plt.savefig("monge_ampere_helmholtz-adapted_mesh2.jpg")

# .. figure:: monge_ampere_helmholtz-adapted_mesh2.jpg
# :figwidth: 60%
# :align: center

u_h = solve_helmholtz(mover.mesh)
error = u_h - u_exact(mover.mesh)
Expand All @@ -235,4 +282,178 @@ def monitor2(mesh):
#
# L2-norm error on moved mesh: 0.00630874419681285
#
# As expected, we do not observe significant changes in final adapted meshes between the
# two approaches. However, the second approach is almost twice longer, as it requires
# solving the PDE in every iteration. In practice it might therefore be sufficient to
# compute the solution on the initial mesh and interpolate it onto adapted meshes at
# every iteration. This is demonstrated in the following implementation:

t0 = time.time()
u_h = solve_helmholtz(mesh)


def monitor_interp_soln(mesh):
V = FunctionSpace(mesh, "CG", 1)
u_h_interp = Function(V).interpolate(u_h)
TV = TensorFunctionSpace(mesh, "CG", 1)
H = RiemannianMetric(TV)
H.compute_hessian(u_h_interp, method="L2")

Hnorm = Function(V, name="Hnorm")
Hnorm.interpolate(sqrt(inner(H, H)))
Hnorm_max = Hnorm.dat.data.max()
m = 1 + alpha * Hnorm / Hnorm_max
return m


mover = MongeAmpereMover(mesh, monitor_interp_soln, method="quasi_newton", rtol=rtol)
mover.move()
print(f"Time taken: {time.time() - t0:.2f} seconds")

# .. code-block:: none
#
# 0 Volume ratio 5.04 Variation (σ/μ) 4.90e-01 Residual 4.93e-01
# 1 Volume ratio 1.98 Variation (σ/μ) 8.07e-02 Residual 7.72e-02
# 2 Volume ratio 1.38 Variation (σ/μ) 4.87e-02 Residual 4.72e-02
# 3 Volume ratio 1.23 Variation (σ/μ) 3.02e-02 Residual 2.78e-02
# 4 Volume ratio 1.19 Variation (σ/μ) 1.99e-02 Residual 1.31e-02
# 5 Volume ratio 1.18 Variation (σ/μ) 1.76e-02 Residual 5.84e-03
# 6 Volume ratio 1.18 Variation (σ/μ) 1.76e-02 Residual 2.66e-03
# 7 Volume ratio 1.18 Variation (σ/μ) 1.77e-02 Residual 1.26e-03
# 8 Volume ratio 1.18 Variation (σ/μ) 1.78e-02 Residual 6.20e-04
# 9 Volume ratio 1.18 Variation (σ/μ) 1.79e-02 Residual 3.16e-04
# 10 Volume ratio 1.18 Variation (σ/μ) 1.79e-02 Residual 1.67e-04
# 11 Volume ratio 1.18 Variation (σ/μ) 1.80e-02 Residual 9.05e-05
# 12 Volume ratio 1.18 Variation (σ/μ) 1.80e-02 Residual 5.07e-05
# 13 Volume ratio 1.18 Variation (σ/μ) 1.80e-02 Residual 2.93e-05
# 14 Volume ratio 1.18 Variation (σ/μ) 1.80e-02 Residual 1.74e-05
# 15 Volume ratio 1.18 Variation (σ/μ) 1.80e-02 Residual 1.06e-05
# 16 Volume ratio 1.18 Variation (σ/μ) 1.80e-02 Residual 6.56e-06
# 17 Volume ratio 1.18 Variation (σ/μ) 1.80e-02 Residual 4.13e-06
# 18 Volume ratio 1.18 Variation (σ/μ) 1.80e-02 Residual 2.63e-06
# 19 Volume ratio 1.18 Variation (σ/μ) 1.80e-02 Residual 1.68e-06
# 20 Volume ratio 1.18 Variation (σ/μ) 1.80e-02 Residual 1.08e-06
# 21 Volume ratio 1.18 Variation (σ/μ) 1.80e-02 Residual 6.98e-07
# 22 Volume ratio 1.18 Variation (σ/μ) 1.80e-02 Residual 4.51e-07
# 23 Volume ratio 1.18 Variation (σ/μ) 1.80e-02 Residual 2.91e-07
# 24 Volume ratio 1.18 Variation (σ/μ) 1.80e-02 Residual 1.88e-07
# Solver converged in 24 iterations.
# Time taken: 6.19 seconds

fig, axes = plt.subplots()
triplot(mover.mesh, axes=axes)
axes.set_aspect(1)
plt.savefig("monge_ampere_helmholtz-adapted_mesh3.jpg")

# .. figure:: monge_ampere_helmholtz-adapted_mesh3.jpg
# :figwidth: 60%
# :align: center

u_h = solve_helmholtz(mover.mesh)
error = u_h - u_exact(mover.mesh)
print("L2-norm error on moved mesh:", sqrt(assemble(dot(error, error) * dx)))

# .. code-block:: none
#
# L2-norm error on moved mesh: 0.008712487462902048
#
# Even though the number of iterations has increased (24 compared to 17 and 19),
# each iteration is slightly cheaper since the PDE is not solved again every time.
# Despite that, the increased total number of iterations lead to a slightly longer total
# runtime. Depending on the PDE in question and the solver used to solve it,
# interpolating the solution may turn out to be slower than recomputing it.
# We also observe that interpolating the solution lead to a significantly larger error.
# In conclusion, in this particular example interpolating the solution field did not
# lead to neither a lower error nor a faster mesh movement step compared to previous
# approaches. Let us explore one more approach.
#
# We may further speed up the mesh movement step by also avoiding recomputing
# the Hessian at every iteration, which, in this particular example, is the most
# expensive part of the monitor function. Similarly to the above example, we compute the
# solution on the initial (uniform) mesh, but now we also compute its Hessian and
# interpolate the Hessian onto adapted meshes.

t0 = time.time()
u_h = solve_helmholtz(mesh)
TV = TensorFunctionSpace(mesh, "CG", 1)
H = RiemannianMetric(TV)
H.compute_hessian(u_h, method="L2")


def monitor_interp_Hessian(mesh):
V = FunctionSpace(mesh, "CG", 1)
TV = TensorFunctionSpace(mesh, "CG", 1)
H_interp = RiemannianMetric(TV).interpolate(H)

Hnorm = Function(V, name="Hnorm")
Hnorm.interpolate(sqrt(inner(H_interp, H_interp)))
Hnorm_max = Hnorm.dat.data.max()
m = 1 + alpha * Hnorm / Hnorm_max
return m


mover = MongeAmpereMover(mesh, monitor_interp_Hessian, method="quasi_newton", rtol=rtol)
mover.move()
print(f"Time taken: {time.time() - t0:.2f} seconds")

# .. code-block:: none
#
# 0 Volume ratio 5.04 Variation (σ/μ) 4.90e-01 Residual 4.93e-01
# 1 Volume ratio 1.92 Variation (σ/μ) 9.68e-02 Residual 9.18e-02
# 2 Volume ratio 1.21 Variation (σ/μ) 2.31e-02 Residual 6.69e-03
# 3 Volume ratio 1.15 Variation (σ/μ) 2.11e-02 Residual 4.14e-05
# 4 Volume ratio 1.15 Variation (σ/μ) 2.12e-02 Residual 4.80e-08
# Solver converged in 4 iterations.
# Time taken: 1.09 seconds

fig, axes = plt.subplots()
triplot(mover.mesh, axes=axes)
axes.set_aspect(1)
plt.savefig("monge_ampere_helmholtz-adapted_mesh4.jpg")

# .. figure:: monge_ampere_helmholtz-adapted_mesh4.jpg
# :figwidth: 60%
# :align: center

u_h = solve_helmholtz(mover.mesh)
error = u_h - u_exact(mover.mesh)
print("L2-norm error on moved mesh:", sqrt(assemble(dot(error, error) * dx)))

# .. code-block:: none
#
# L2-norm error on moved mesh: 0.008385305585746483
#
# The mesh movement step now only took 4 iterations, with a total runtime of only 1.09
# seconds, which is up to five times shorter than previous examples. The final error is
# again larger than the example where the solution is recomputed at every iteration,
# but is smaller than the example where we interpolated the solution field.
#
# We can summarise these results in the following table:
#
# .. table::
#
# ============================ ========= ==============
# Monitor function Error CPU time (s)
# ============================ ========= ==============
# ``monitor_exact`` 0.00596 2.52
# ``monitor_solve`` 0.00631 6.17
# ``monitor_interp_soln`` 0.00871 6.19
# ``monitor_interp_Hessian`` 0.00839 1.09
# ============================ ========= ==============
#
# In this demo we demonstrated several examples of monitor functions and briefly
# evaluated their performance. Each approach has inherent advantages and limitations
# which should be considered for every new problem of interest. PDEs that are highly
# sensitive to changes in local resolution may require recomputing the solution at
# every iteration. In other cases we may obtain adequate results by defining a monitor
# function that computes the solution less frequently, or even only once per iteration.
# Movement allows and encourages such experimentation with different monitor functions.
#
# This tutorial can be dowloaded as a `Python script <monge_ampere_helmholtz.py>`__.
#
# .. rubric:: Footnotes
#
# .. [#f1] Running this script as it is will yield a slightly longer runtime of the
# first approach (3.14 seconds on our machine) due to the additional overhead of the
# first-time JIT compilation. The reported runtime of 2.52 seconds is obtained by
# running the first example twice in a row and taking the second runtime.

0 comments on commit f8ac7b9

Please sign in to comment.