From f8ac7b9b7a632f85c57a607e97f12c3c7c5da7c1 Mon Sep 17 00:00:00 2001 From: Davor Dundovic <33790330+ddundo@users.noreply.github.com> Date: Thu, 9 Jan 2025 19:04:56 +0100 Subject: [PATCH] Expand MA Helmholtz demo to interpolate soln field (#145) Closes #27. --- demos/monge_ampere_helmholtz.py | 233 +++++++++++++++++++++++++++++++- 1 file changed, 227 insertions(+), 6 deletions(-) diff --git a/demos/monge_ampere_helmholtz.py b/demos/monge_ampere_helmholtz.py index c6bb099..3751860 100644 --- a/demos/monge_ampere_helmholtz.py +++ b/demos/monge_ampere_helmholtz.py @@ -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))) @@ -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") @@ -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 @@ -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: @@ -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) @@ -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) @@ -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 `__. +# +# .. 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.