From f0be9bdb147ff2fa88f57e85179829a1bb0be544 Mon Sep 17 00:00:00 2001 From: ddundo Date: Tue, 13 Aug 2024 08:13:44 +0000 Subject: [PATCH 01/12] #55 #144: Add bubble shear demo --- animate/metric.py | 12 +- demos/bubble_shear.py | 430 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 436 insertions(+), 6 deletions(-) create mode 100644 demos/bubble_shear.py diff --git a/animate/metric.py b/animate/metric.py index cbcc4e7..141e9fa 100644 --- a/animate/metric.py +++ b/animate/metric.py @@ -84,12 +84,12 @@ def __init__(self, function_space, *args, **kwargs): if isinstance(fs.dof_count, Iterable): raise ValueError("Riemannian metric cannot be built in a mixed space.") self._check_space() - rank = len(fs.dof_dset.dim) - if rank != 2: - raise ValueError( - "Riemannian metric should be matrix-valued," - f" not rank-{rank} tensor-valued." - ) + # rank = len(fs.dof_dset.dim) + # if rank != 2: + # raise ValueError( + # "Riemannian metric should be matrix-valued," + # f" not rank-{rank} tensor-valued." + # ) # Stash mesh data plex = mesh.topology_dm.clone() diff --git a/demos/bubble_shear.py b/demos/bubble_shear.py new file mode 100644 index 0000000..7560e2a --- /dev/null +++ b/demos/bubble_shear.py @@ -0,0 +1,430 @@ +# On-the-fly time-dependent mesh adaptation +########################################### + +# In this demo we consider the 2-dimensional version of the mesh adaptation experiment +# presented in :cite:. The problem comprises a bubble of tracer +# concentration field advected by a time-varying flow. +# We will consider two different mesh adaptation strategies: the classical mesh +# adaptation algorithm, which adapts the mesh several times throughout the simulation, +# before solving the advection equation, and the metric advection algorithm, which +# advects the initial metric tensor along the flow in order to predict where to +# prescribe fine resolution in the future. +# +# We begin by defining the advection problem. We consider the advection equation +# +# .. math:: +# \begin{array}{rl} +# \frac{\partial c}{\partial t} + \nabla\cdot(\mathbf{u}c)=0& \text{in}\:\Omega,\\ +# c=0 & \text{on}\:\partial\Omega,\\ +# c=c_0(x,y) & \text{at}\:t=0, +# \end{array}, +# +# where :math:`c=c(x,y,t)` is the sought tracer concentration, +# :math:`\mathbf{u}=\mathbf{u}(x,y,t)` is the background velocity field, and +# :math:`\Omega=[0, 1]^2` is the spatial domain of interest. +# +# The background velocity field :math:`\mathbf{u}(x, y, t)` is chosen to be periodic in +# time, and is given by +# +# .. math:: +# \mathbf{u}(x, y, t) := \left(2\sin^2(\pi x)\sin(2\pi y), -\sin(2\pi x)\sin^2(\pi y) \right) \cos(2\pi t/T), +# +# where :math:`T` is the period. Note that the velocity field is solenoidal. At each +# timestep of the simulation we will update this field so we define a function that will +# return its vector expression. :: + +from firedrake import * + +T = 6.0 + + +def velocity_expression(mesh, t): + x, y = SpatialCoordinate(mesh) + u_expr = as_vector( + [ + 2 * sin(pi * x) ** 2 * sin(2 * pi * y) * cos(2 * pi * t / T), + -sin(2 * pi * x) * sin(pi * y) ** 2 * cos(2 * pi * t / T), + ] + ) + return u_expr + + +# We proceed similarly with prescribing initial conditions. At :math:`t=0`, we +# initialise the tracer concentration :math:`c_0 = c(x, y, 0)` to be :math:`1` inside +# a circular region of radius :math:`r_0=0.15` centred at :math:`(x_0, y_0)=(0.5, 0.65)` +# and :math:`0` elsewhere in the domain. Note that such a discontinuous function will +# not be represented well on a coarse uniform mesh. :: + + +def get_initial_condition(mesh): + x, y = SpatialCoordinate(mesh) + ball_r, ball_x0, ball_y0 = 0.15, 0.5, 0.65 + r = sqrt(pow(x - ball_x0, 2) + pow(y - ball_y0, 2)) + c0 = Function(FunctionSpace(mesh, "CG", 1)) + c0.interpolate(conditional(r < ball_r, 1.0, 0.0)) + # return conditional(r < ball_r0, 1.0, 0.0) + return c0 + + +# Now we are ready to solve the advection problem. Since we will solve the problem +# several times, we will wrap up the code in a function that we can easily call +# for each subinterval. The function takes the mesh over which to solve the problem, +# time interval :math:`[t_{\text{start}}, t_{\text{end}}]` over which to solve it, and +# the initial condition :math:`c_0 = c(x, y, t_{\text{start}})`. The function returns +# the solution :math:`c(x, y, t_{\text{end}})`. :: + + +def run_simulation(mesh, t_start, t_end, c0): + Q = FunctionSpace(mesh, "CG", 1) + V = VectorFunctionSpace(mesh, "CG", 1) + R = FunctionSpace(mesh, "R", 0) + + c_ = Function(Q).project(c0) # project initial condition onto the current mesh + c = Function(Q) # solution at current timestep + + u_ = Function(V).interpolate(velocity_expression(mesh, t_start)) # vel. at t_start + u = Function(V) # velocity field at current timestep + + # SUPG stabilisation parameters + D = Function(R).assign(0.1) # diffusivity coefficient + h = CellSize(mesh) # mesh cell size + U = sqrt(dot(u, u)) # velocity magnitude + tau = 0.5 * h / U + tau = min_value(tau, U * h / (6 * D)) + + # Apply SUPG stabilisation + phi = TestFunction(Q) + phi += tau * dot(u, grad(phi)) + + # Time-stepping parameters + dt = Function(R).assign(0.01) # timestep size + theta = Function(R).assign(0.5) # Crank-Nicolson implicitness + + # Variational form of the advection equation + trial = TrialFunction(Q) + a = inner(trial, phi) * dx + dt * theta * inner(dot(u, grad(trial)), phi) * dx + L = inner(c_, phi) * dx - dt * (1 - theta) * inner(dot(u_, grad(c_)), phi) * dx + + # Define variational problem + lvp = LinearVariationalProblem(a, L, c, bcs=DirichletBC(Q, 0.0, "on_boundary")) + lvs = LinearVariationalSolver(lvp) + + # Integrate from t_start to t_end + t = t_start + float(dt) + while t < t_end + 0.5 * float(dt): + # Update the background velocity field at the current timestep + u.interpolate(velocity_expression(mesh, t)) + + # Solve the advection equation + lvs.solve() + + # Update the solution at the previous timestep + c_.assign(c) + u_.assign(u) + t += float(dt) + + return c + + +# Finally, we are ready to run the simulation. We will first solve the entire +# problem over two uniform meshes: one with 32 elements in each direction and another +# with 128 elements in each direction. Since the flow reverts to its initial state at +# time :math:`t=T/2`, we run the simulations over the interval :math:`[0, T/2]`. In +# order to compare the efficacy of mesh adaptation methods, we will also keep track +# of the time it takes to run the simulation. :: + +import time + +simulation_end_time = T / 2.0 + +# mesh_coarse = UnitSquareMesh(32, 32) +mesh_coarse = UnitSquareMesh(256, 256) +c0_coarse = get_initial_condition(mesh_coarse) +cpu_time_coarse = time.time() +c_coarse_final = run_simulation(mesh_coarse, 0.0, simulation_end_time, c0_coarse) +cpu_time_coarse = time.time() - cpu_time_coarse + +mesh_fine = UnitSquareMesh(128, 128) +c0_fine = get_initial_condition(mesh_fine) +cpu_time_fine = time.time() +c_fine_final = run_simulation(mesh_fine, 0.0, simulation_end_time, c0_fine) +cpu_time_fine = time.time() - cpu_time_fine + +print(f"CPU time on the coarse mesh: {cpu_time_coarse:.2f} s") +print(f"CPU time on the fine mesh: {cpu_time_fine:.2f} s") + +# We can now visualise the final concentration fields on the two meshes. :: + +import matplotlib.pyplot as plt +from firedrake.pyplot import * + +fig, axes = plt.subplots(1, 2, figsize=(8, 5), sharey=True) +for ax, c, mesh in zip(axes, [c_coarse_final, c_fine_final], [mesh_coarse, mesh_fine]): + tripcolor(c, axes=ax) + ax.set_title(f"{mesh.num_cells()} mesh elements") + ax.set_xlim(0, 1) + ax.set_ylim(0, 1) + ax.set_aspect("equal") + ax.set_xlabel("x") + ax.set_ylabel("y") +fig.tight_layout() +fig.savefig("bubble_shear-uniform.jpg", dpi=300, bbox_inches="tight") + +# .. figure:: bubble_shear-uniform.jpg +# :figwidth: 80% +# :align: center +# +# We observe that the concentration fields have not returned to their initial value +# at the end of the simulation, as we would have expected. This is particularly obvious +# in the coarse resolution simulation. This is due to the addition of the SUPG +# stabilisation to the weak form of the advection equation, which adds numerical +# diffusion. Numerical diffusion is necessary for numerical stability and for preventing +# oscillations, but it also makes the solution irreversible. The amount of difussion +# added is related to the grid Péclet number :math:`Pe = U\,h/2D`: the coarser the mesh +# is, the more diffusion is added. We encourage the reader to verify this by running the +# simulation on a sequence of finer uniform meshes. +# +# In ... :: + + +def compute_rel_error(c_init, c_final): + init_l2_norm = norm(c_init, norm_type="L2") + abs_l2_error = errornorm(c_init, c_final, norm_type="L2") + return 100 * abs_l2_error / init_l2_norm + + +coarse_error = compute_rel_error(c0_coarse, c_coarse_final) +print(f"Relative L2 error on the coarse mesh: {coarse_error:.2f}%") +fine_error = compute_rel_error(c0_fine, c_fine_final) +print(f"Relative L2 error on the fine mesh: {fine_error:.2f}%") + +# .. code-block:: console +# +# Relative L2 error on the coarse mesh: 58.85% +# Relative L2 error on the fine mesh: 34.60% +# +# TODO Explain the error +# +# Instead of running the simulation on a very fine mesh, we can use mesh adaptation +# techniques to refine the mesh only in regions and at times where that is necessary. +# For the purposes of this demo, we are going to adapt the mesh 20 times throughout the +# simulation, at equal time intervals (i.e., every 0.15s of simulation time). :: + +num_adaptations = 20 +interval_length = simulation_end_time / num_adaptations + +# As mentioned in the introduction, we shall demonstrate two different mesh adaptation +# strategies. Let us begin with the classical mesh adaptation algorithm, which adapts +# each mesh before solving the advection equation over the corresponding subinterval. +# Here we will use the :class:`RiemannianMetric` class to define the metric, which we +# will compute based on the Hessian of the concentration field. Let us therefore define +# a function which takes the original mesh and the concentration field as arguments. We +# will also define parameters for computing the metric. :: + +from animate.adapt import adapt +from animate.metric import RiemannianMetric + +metric_params = { + "dm_plex_metric": { + "target_complexity": 2000.0, + "p": 2.0, + "h_min": 1e-04, # minimum edge length + "h_max": 1.0, # maximum edge length + } +} + + +def adapt_classical(mesh, c): + P1_ten = TensorFunctionSpace(mesh, "CG", 1) + metric = RiemannianMetric(P1_ten) + metric.set_parameters(metric_params) + metric.compute_hessian(c) + metric.normalise() + adapted_mesh = adapt(mesh, metric) + return adapted_mesh + + +# We can now run the simulation over each subinterval. :: + +mesh = mesh_coarse +c = get_initial_condition(mesh_coarse) +cpu_time_classical = time.time() +for i in range(num_adaptations): + t0 = i * interval_length # subinterval start time + t1 = (i + 1) * interval_length # subinterval end time + + # Adapt the mesh based on the concentration field at t0 + mesh = adapt_classical(mesh, c) + + # Solve the advection equation over the subinterval (t0, t1] + c = run_simulation(mesh, t0, t1, c) +cpu_time_classical = time.time() - cpu_time_classical +print(f"CPU time with classical adaptation algorithm: {cpu_time_classical:.2f} s") + +# Now let us plot the final adapted mesh and final concentration field computed on it. +# We will also compute the relative L2 error. :: + +fig, axes = plt.subplots(1, 2, figsize=(8, 5), sharey=True) +triplot(mesh, axes=axes[0]) +tripcolor(c, axes=axes[1]) + +for ax in axes.flat: + ax.set_xlim(0, 1) + ax.set_ylim(0, 1) + ax.set_aspect("equal") + ax.set_xlabel("x") + ax.set_ylabel("y") +fig.tight_layout() +fig.savefig("bubble_shear-classical.jpg", dpi=300, bbox_inches="tight") + +# Redefine the initial condition on the final adapted mesh +c0 = get_initial_condition(mesh) +classical_c0 = Function(FunctionSpace(mesh, "CG", 1)).interpolate(c0) +classical_error = compute_rel_error(classical_c0, c) + +print(f"Relative L2 error with classical adaptation algorithm: {classical_error:.2f}%") + +# .. code-block:: console +# +# Relative L2 error with classical adaptation algorithm: 26.54% +# +# .. figure:: bubble_shear-classical.jpg +# :figwidth: 80% +# :align: center +# +# TODO talk about this :: + + +def adapt_metric_advection(mesh, t_start, t_end, c): + P1_ten = TensorFunctionSpace(mesh, "CG", 1) + metric = RiemannianMetric(P1_ten) + metric_ = RiemannianMetric(P1_ten) + metric_intersect = RiemannianMetric(P1_ten) + + for mtrc in [metric, metric_, metric_intersect]: + mtrc.set_parameters(metric_params) + # metric_.set_parameters(mp) + metric_.compute_hessian(c) + metric_.normalise() + + Q = FunctionSpace(mesh, "CG", 1) + m_ = Function(Q) + m = Function(Q) + h = Function(Q) + + V = VectorFunctionSpace(mesh, "CG", 1) + u = Function(V) + u_ = Function(V).interpolate(velocity_expression(mesh, t_start)) + + R = FunctionSpace(mesh, "R", 0) + dt = Function(R).assign(0.01) # timestep size + theta = Function(R).assign(0.5) # Crank-Nicolson implicitness + + # Metric advection by component + trial = TrialFunction(Q) + phi = TestFunction(Q) + a = inner(trial, phi) * dx + dt * theta * inner(dot(u, grad(trial)), phi) * dx + L = inner(m_, phi) * dx - dt * (1 - theta) * inner(dot(u_, grad(m_)), phi) * dx + lvp = LinearVariationalProblem(a, L, m, bcs=DirichletBC(Q, h, "on_boundary")) + lvs = LinearVariationalSolver(lvp) + + t = t_start + float(dt) + while t < t_end + 0.5 * float(dt): + u.interpolate(velocity_expression(mesh, t)) + + # Advect each metric component + for i in range(2): + for j in range(2): + h_max = metric_params["dm_plex_metric"]["h_max"] + h.assign(1.0 / h_max**2 if i == j else 0.0) + m_.assign(metric_.sub(2 * i + j)) + lvs.solve() + metric.sub(2 * i + j).assign(m) + + # Ensure symmetry + m_.assign(0.5 * (metric.sub(1) + metric.sub(2))) + metric.sub(1).assign(m_) + metric.sub(2).assign(m_) + + # Intersect metrics at each timestep + metric.enforce_spd(restrict_sizes=True, restrict_anisotropy=True) + metric_intersect.intersect(metric) + + metric_.assign(metric) + u_.assign(u) + t += float(dt) + + metric_intersect.normalise() + amesh = adapt(mesh, metric_intersect) + return amesh + + +mesh = UnitSquareMesh(32, 32) +c = get_initial_condition(mesh) +cpu_time_metric_advection = time.time() +for i in range(num_adaptations): + t0 = i * interval_length # subinterval start time + t1 = (i + 1) * interval_length # subinterval end time + + # Adapt the mesh based on the concentration field at t0 + mesh = adapt_metric_advection(mesh, t0, t1, c) + + # Solve the advection equation over the subinterval (t0, t1] + c = run_simulation(mesh, t0, t1, c) +cpu_time_metric_advection = time.time() - cpu_time_metric_advection + +print(f"CPU time with metric advection: {cpu_time_metric_advection:.2f} s") + +fig, axes = plt.subplots(1, 2, figsize=(8, 5), sharey=True) +triplot(mesh, axes=axes[0]) +tripcolor(c, axes=axes[1]) + +for ax in axes.flat: + ax.set_xlim(0, 1) + ax.set_ylim(0, 1) + ax.set_aspect("equal") + ax.set_xlabel("x") + ax.set_ylabel("y") +fig.tight_layout() +fig.savefig("bubble_shear-metric_advection.jpg", dpi=300, bbox_inches="tight") + +# Redefine the initial condition on the final adapted mesh +c0 = get_initial_condition(mesh) +metric_adv_c0 = Function(FunctionSpace(mesh, "CG", 1)).interpolate(c0) +metric_adv_error = compute_rel_error(metric_adv_c0, c) + +print(f"Relative L2 error with metric advection: {metric_adv_error:.2f}%") + +# .. code-block:: console +# +# Relative L2 error with metric advection: 29.31% +# +# .. figure:: bubble_shear-metric_advection.jpg +# :figwidth: 80% +# :align: center + +fig, ax = plt.subplots() +errors = [coarse_error, fine_error, classical_error, metric_adv_error] +cpu_times = [ + cpu_time_coarse, + cpu_time_fine, + cpu_time_classical, + cpu_time_metric_advection, +] +ax.plot(cpu_times, errors, "o") +# label each point +for i, txt in enumerate(["Coarse mesh", "Fine mesh", "Classical", "Metric advection"]): + ax.annotate( + txt, + (cpu_times[i], errors[i]), + textcoords="offset points", + xytext=(0, 10), + ha="center", + ) +ax.set_xlabel("CPU time (s)") +ax.set_ylabel("Relative L2 error (%)") +fig.savefig("bubble_shear-time_error.jpg", dpi=300, bbox_inches="tight") + +# .. figure:: bubble_shear-time_error.jpg +# :figwidth: 80% +# :align: center From 3284d6e9f1c853cc9f1dd6da3aa2ad565158d7a2 Mon Sep 17 00:00:00 2001 From: ddundo Date: Thu, 12 Sep 2024 11:09:49 +0000 Subject: [PATCH 02/12] #55: Finalise demo --- demos/bubble_shear.py | 325 ++++++++++++++++++++++---------------- demos/demo_references.bib | 14 ++ 2 files changed, 207 insertions(+), 132 deletions(-) diff --git a/demos/bubble_shear.py b/demos/bubble_shear.py index 7560e2a..8ac8fd8 100644 --- a/demos/bubble_shear.py +++ b/demos/bubble_shear.py @@ -1,14 +1,17 @@ # On-the-fly time-dependent mesh adaptation -########################################### - +# ######################################### +# # In this demo we consider the 2-dimensional version of the mesh adaptation experiment # presented in :cite:. The problem comprises a bubble of tracer # concentration field advected by a time-varying flow. # We will consider two different mesh adaptation strategies: the classical mesh -# adaptation algorithm, which adapts the mesh several times throughout the simulation, -# before solving the advection equation, and the metric advection algorithm, which +# adaptation algorithm, which adapts the mesh several times throughout the simulation +# based on the solution at the current time, and the metric advection algorithm, which # advects the initial metric tensor along the flow in order to predict where to -# prescribe fine resolution in the future. +# prescribe fine resolution in the future. Both algorithms are an example of +# *on-the-fly* time-dependent mesh adaptation algorithms, where the mesh is adapted +# before each subinterval of the simulation, as opposed to fixed-point iteration +# algorithms, where the mesh is iteratively adapted at the end of the simulation. # # We begin by defining the advection problem. We consider the advection equation # @@ -21,7 +24,8 @@ # # where :math:`c=c(x,y,t)` is the sought tracer concentration, # :math:`\mathbf{u}=\mathbf{u}(x,y,t)` is the background velocity field, and -# :math:`\Omega=[0, 1]^2` is the spatial domain of interest. +# :math:`\Omega=[0, 1]^2` is the spatial domain of interest with boundary +# :math:`\partial\Omega`. # # The background velocity field :math:`\mathbf{u}(x, y, t)` is chosen to be periodic in # time, and is given by @@ -71,7 +75,14 @@ def get_initial_condition(mesh): # for each subinterval. The function takes the mesh over which to solve the problem, # time interval :math:`[t_{\text{start}}, t_{\text{end}}]` over which to solve it, and # the initial condition :math:`c_0 = c(x, y, t_{\text{start}})`. The function returns -# the solution :math:`c(x, y, t_{\text{end}})`. :: +# the solution :math:`c(x, y, t_{\text{end}})`. +# +# Note that we must include streamline upwind Petrov Galerkin (SUPG) stabilisation in +# the test function in order to ensure numerical stability. :: + +# Time-stepping parameters +# dt = 0.01 # timestep size +# theta = 0.5 # Crank-Nicolson implicitness def run_simulation(mesh, t_start, t_end, c0): @@ -85,14 +96,14 @@ def run_simulation(mesh, t_start, t_end, c0): u_ = Function(V).interpolate(velocity_expression(mesh, t_start)) # vel. at t_start u = Function(V) # velocity field at current timestep - # SUPG stabilisation parameters + # SUPG stabilisation D = Function(R).assign(0.1) # diffusivity coefficient h = CellSize(mesh) # mesh cell size U = sqrt(dot(u, u)) # velocity magnitude tau = 0.5 * h / U tau = min_value(tau, U * h / (6 * D)) - # Apply SUPG stabilisation + # Apply SUPG stabilisation to the test function phi = TestFunction(Q) phi += tau * dot(u, grad(phi)) @@ -129,29 +140,17 @@ def run_simulation(mesh, t_start, t_end, c0): # Finally, we are ready to run the simulation. We will first solve the entire # problem over two uniform meshes: one with 32 elements in each direction and another # with 128 elements in each direction. Since the flow reverts to its initial state at -# time :math:`t=T/2`, we run the simulations over the interval :math:`[0, T/2]`. In -# order to compare the efficacy of mesh adaptation methods, we will also keep track -# of the time it takes to run the simulation. :: - -import time +# time :math:`t=T/2`, we run the simulations over the interval :math:`[0, T/2]`. :: simulation_end_time = T / 2.0 -# mesh_coarse = UnitSquareMesh(32, 32) -mesh_coarse = UnitSquareMesh(256, 256) +mesh_coarse = UnitSquareMesh(32, 32) c0_coarse = get_initial_condition(mesh_coarse) -cpu_time_coarse = time.time() c_coarse_final = run_simulation(mesh_coarse, 0.0, simulation_end_time, c0_coarse) -cpu_time_coarse = time.time() - cpu_time_coarse mesh_fine = UnitSquareMesh(128, 128) c0_fine = get_initial_condition(mesh_fine) -cpu_time_fine = time.time() c_fine_final = run_simulation(mesh_fine, 0.0, simulation_end_time, c0_fine) -cpu_time_fine = time.time() - cpu_time_fine - -print(f"CPU time on the coarse mesh: {cpu_time_coarse:.2f} s") -print(f"CPU time on the fine mesh: {cpu_time_fine:.2f} s") # We can now visualise the final concentration fields on the two meshes. :: @@ -161,7 +160,7 @@ def run_simulation(mesh, t_start, t_end, c0): fig, axes = plt.subplots(1, 2, figsize=(8, 5), sharey=True) for ax, c, mesh in zip(axes, [c_coarse_final, c_fine_final], [mesh_coarse, mesh_fine]): tripcolor(c, axes=ax) - ax.set_title(f"{mesh.num_cells()} mesh elements") + ax.set_title(f"{mesh.num_vertices()} mesh vertices") ax.set_xlim(0, 1) ax.set_ylim(0, 1) ax.set_aspect("equal") @@ -184,7 +183,8 @@ def run_simulation(mesh, t_start, t_end, c0): # is, the more diffusion is added. We encourage the reader to verify this by running the # simulation on a sequence of finer uniform meshes. # -# In ... :: +# In order to quantify the above observation, we will compute the relative L2 error +# between the initial condition and the final concentration field on the final mesh. :: def compute_rel_error(c_init, c_final): @@ -196,40 +196,59 @@ def compute_rel_error(c_init, c_final): coarse_error = compute_rel_error(c0_coarse, c_coarse_final) print(f"Relative L2 error on the coarse mesh: {coarse_error:.2f}%") fine_error = compute_rel_error(c0_fine, c_fine_final) -print(f"Relative L2 error on the fine mesh: {fine_error:.2f}%") +print(f"Relative L2 error on the fine mesh: {fine_error:.2f}%.") # .. code-block:: console # -# Relative L2 error on the coarse mesh: 58.85% -# Relative L2 error on the fine mesh: 34.60% -# -# TODO Explain the error +# Relative L2 error on the coarse mesh: 56.52% +# Relative L2 error on the fine mesh: 32.29% # -# Instead of running the simulation on a very fine mesh, we can use mesh adaptation -# techniques to refine the mesh only in regions and at times where that is necessary. -# For the purposes of this demo, we are going to adapt the mesh 20 times throughout the -# simulation, at equal time intervals (i.e., every 0.15s of simulation time). :: +# Since accurate simulations require very fine resolution, which may be computationally +# too prohitibive, we will now demonstrate how to use mesh adaptation techniques to +# refine the mesh only in regions and at times where that is necessary. +# For the purposes of this demo, we are going to adapt the mesh 15 times throughout the +# simulation, at equal time intervals (i.e., every 0.2s of simulation time). :: -num_adaptations = 20 +num_adaptations = 15 interval_length = simulation_end_time / num_adaptations +# We will also define a function that will allow us to easily plot the adapted mesh, as +# well as the solution fields at the beginning and end of each subinterval. :: + + +def plot_mesh(mesh, c0, c1, i, method): + fig, ax = plt.subplots(1, 3, sharey=True, figsize=(12, 4)) + triplot(mesh, axes=ax[0]) + tripcolor(c0, axes=ax[1]) + tripcolor(c1, axes=ax[2]) + ax[0].set_title(f"Mesh {i} ({mesh.num_vertices()} vertices)") + ax[1].set_title(f"Solution at t={i*interval_length:.1f}s") + ax[2].set_title(f"Solution at t={(i+1)*interval_length:.1f}s") + for axes in ax: + axes.set_xlim(0, 1) + axes.set_ylim(0, 1) + axes.set_aspect("equal") + fig.savefig(f"bubble_shear-{method}_{i}.jpg", dpi=300, bbox_inches="tight") + plt.close(fig) + + # As mentioned in the introduction, we shall demonstrate two different mesh adaptation # strategies. Let us begin with the classical mesh adaptation algorithm, which adapts # each mesh before solving the advection equation over the corresponding subinterval. # Here we will use the :class:`RiemannianMetric` class to define the metric, which we # will compute based on the Hessian of the concentration field. Let us therefore define -# a function which takes the original mesh and the concentration field as arguments. We -# will also define parameters for computing the metric. :: +# a function which takes the original mesh and the concentration field as arguments, and +# returns the adapted mesh. We will also define parameters for computing the metric. :: from animate.adapt import adapt from animate.metric import RiemannianMetric metric_params = { "dm_plex_metric": { - "target_complexity": 2000.0, - "p": 2.0, - "h_min": 1e-04, # minimum edge length - "h_max": 1.0, # maximum edge length + "target_complexity": 1500.0, + "p": 2.0, # normalisation order + "h_min": 1e-04, # minimum allowed edge length + "h_max": 1.0, # maximum allowed edge length } } @@ -244,55 +263,78 @@ def adapt_classical(mesh, c): return adapted_mesh -# We can now run the simulation over each subinterval. :: +# We can now run the simulation, but we adapt the mesh before each subinterval. We begin +# with a coarse uniform mesh and track the number of vertices of each adapted mesh. :: -mesh = mesh_coarse -c = get_initial_condition(mesh_coarse) -cpu_time_classical = time.time() +mesh = UnitSquareMesh(32, 32) +mesh_numVertices = [] +c = get_initial_condition(mesh) for i in range(num_adaptations): t0 = i * interval_length # subinterval start time t1 = (i + 1) * interval_length # subinterval end time # Adapt the mesh based on the concentration field at t0 mesh = adapt_classical(mesh, c) + mesh_numVertices.append(mesh.num_vertices()) + + # Make a copy of the initial condition for plotting purposes + c0 = c.copy(deepcopy=True) # Solve the advection equation over the subinterval (t0, t1] c = run_simulation(mesh, t0, t1, c) -cpu_time_classical = time.time() - cpu_time_classical -print(f"CPU time with classical adaptation algorithm: {cpu_time_classical:.2f} s") -# Now let us plot the final adapted mesh and final concentration field computed on it. -# We will also compute the relative L2 error. :: - -fig, axes = plt.subplots(1, 2, figsize=(8, 5), sharey=True) -triplot(mesh, axes=axes[0]) -tripcolor(c, axes=axes[1]) + # Plot the adapted mesh and the concentration field at t0 and t1 + plot_mesh(mesh, c0, c, i, "classical") -for ax in axes.flat: - ax.set_xlim(0, 1) - ax.set_ylim(0, 1) - ax.set_aspect("equal") - ax.set_xlabel("x") - ax.set_ylabel("y") -fig.tight_layout() -fig.savefig("bubble_shear-classical.jpg", dpi=300, bbox_inches="tight") +# Now let us examine the final adapted mesh and final concentration field computed on it. +# We will also compute the relative L2 error. :: # Redefine the initial condition on the final adapted mesh c0 = get_initial_condition(mesh) classical_c0 = Function(FunctionSpace(mesh, "CG", 1)).interpolate(c0) classical_error = compute_rel_error(classical_c0, c) -print(f"Relative L2 error with classical adaptation algorithm: {classical_error:.2f}%") +print( + f"Classical mesh adaptation.\n" + f" Avg. number of vertices: {np.average(mesh_numVertices):.1f}\n" + f" Relative L2 error: {classical_error:.2f}%" +) # .. code-block:: console # -# Relative L2 error with classical adaptation algorithm: 26.54% +# Classical mesh adaptation. +# Avg. number of vertices: 2302.3 +# Relative L2 error: 31.76% # -# .. figure:: bubble_shear-classical.jpg +# .. figure:: bubble_shear-classical_14.jpg # :figwidth: 80% # :align: center # -# TODO talk about this :: +# As we can see, the relative L2 error is lower than the one obtained with the uniform +# fine mesh, even though the average number of vertices is more than 8 times smaller. +# This demonstrates the effectiveness of mesh adaptation, which we can also observe in +# the figure above. We see that final mesh has clearly been refined around the bubble +# at the beginning of the final subinterval and coarsened elsewhere. +# +# However, we also observe that the bubble has advected out of the fine resolution +# region by the end of the subinterval. This is a common occurence in time-dependent +# mesh adaptation, known as the *lagging mesh* problem, where the mesh is said to lag +# with respect to the solution. +# Ensuring that the bubble remains well-resolved throughout the simulation +# is not a trivial task, as it requires predicting where the bubble will be in the +# future. Earliest attempts at preventing the lagging mesh problem introduced a safety +# margin around the fine-resolution region. While potentially very effective, depending +# on the size of the margin, this approach is also likely to be inefficient as it may +# prescribe fine resolution in regions where it is not needed. +# +# For advection-dominated problems, as is the case here, a *metric advection* algorithm +# has been proposed in :cite:. The idea is to still compute the metric +# based on the solution at the current time, but then to advect the metric along the +# flow in order to predict where to prescribe fine resolution in the future. By +# combining the advected metrics in time, we obtain a final metric that is +# representative of the evolving solution throughout the subinterval. We achieve this in +# the following function, where we solve the before-seen advection equation for each +# component of the metric tensor. :: def adapt_metric_advection(mesh, t_start, t_end, c): @@ -301,16 +343,16 @@ def adapt_metric_advection(mesh, t_start, t_end, c): metric_ = RiemannianMetric(P1_ten) metric_intersect = RiemannianMetric(P1_ten) + # Compute the Hessian metric at t_start for mtrc in [metric, metric_, metric_intersect]: mtrc.set_parameters(metric_params) - # metric_.set_parameters(mp) metric_.compute_hessian(c) - metric_.normalise() + metric_.enforce_spd(restrict_sizes=False, restrict_anisotropy=False) Q = FunctionSpace(mesh, "CG", 1) - m_ = Function(Q) - m = Function(Q) - h = Function(Q) + h_bc = Function(Q) # used to set the boundary condition below + m_ = Function(Q) # metric component at previous timestep + m = Function(Q) # metric component at current timestep V = VectorFunctionSpace(mesh, "CG", 1) u = Function(V) @@ -320,36 +362,45 @@ def adapt_metric_advection(mesh, t_start, t_end, c): dt = Function(R).assign(0.01) # timestep size theta = Function(R).assign(0.5) # Crank-Nicolson implicitness - # Metric advection by component - trial = TrialFunction(Q) + # SUPG stabilisation + D = Function(R).assign(0.1) + h = CellSize(mesh) + U = sqrt(dot(u, u)) + tau = 0.5 * h / U + tau = min_value(tau, U * h / (6 * D)) + + # Apply SUPG stabilisation phi = TestFunction(Q) + phi += tau * dot(u, grad(phi)) + + # Variational form of the advection equation, as before + trial = TrialFunction(Q) a = inner(trial, phi) * dx + dt * theta * inner(dot(u, grad(trial)), phi) * dx L = inner(m_, phi) * dx - dt * (1 - theta) * inner(dot(u_, grad(m_)), phi) * dx - lvp = LinearVariationalProblem(a, L, m, bcs=DirichletBC(Q, h, "on_boundary")) + lvp = LinearVariationalProblem(a, L, m, bcs=DirichletBC(Q, h_bc, "on_boundary")) lvs = LinearVariationalSolver(lvp) + # Integrate from t_start to t_end t = t_start + float(dt) while t < t_end + 0.5 * float(dt): u.interpolate(velocity_expression(mesh, t)) - # Advect each metric component - for i in range(2): - for j in range(2): - h_max = metric_params["dm_plex_metric"]["h_max"] - h.assign(1.0 / h_max**2 if i == j else 0.0) - m_.assign(metric_.sub(2 * i + j)) - lvs.solve() - metric.sub(2 * i + j).assign(m) - - # Ensure symmetry - m_.assign(0.5 * (metric.sub(1) + metric.sub(2))) - metric.sub(1).assign(m_) - metric.sub(2).assign(m_) - - # Intersect metrics at each timestep + # Advect each metric component M_ij individually + for i, j in [(0, 0), (0, 1), (1, 1)]: + # print(i, j) + h_max = metric_params["dm_plex_metric"]["h_max"] + h_bc.assign(1.0 / h_max**2 if i == j else 0.0) + m_.assign(metric_.sub(2 * i + j)) + lvs.solve() + metric.sub(2 * i + j).assign(m) + # Metric is symmetric so we can copy the other off-diagonal component + metric.sub(2).assign(metric.sub(1)) + + # Intersect metrics at every timestep metric.enforce_spd(restrict_sizes=True, restrict_anisotropy=True) metric_intersect.intersect(metric) + # Update fields at the previous timestep metric_.assign(metric) u_.assign(u) t += float(dt) @@ -359,72 +410,82 @@ def adapt_metric_advection(mesh, t_start, t_end, c): return amesh +# We can now run the simulation over the entire time interval, but now we will adapt the +# mesh at the beginning of each subinterval using the metric advection algorithm. :: + mesh = UnitSquareMesh(32, 32) +mesh_numVertices = [] c = get_initial_condition(mesh) -cpu_time_metric_advection = time.time() for i in range(num_adaptations): t0 = i * interval_length # subinterval start time t1 = (i + 1) * interval_length # subinterval end time # Adapt the mesh based on the concentration field at t0 mesh = adapt_metric_advection(mesh, t0, t1, c) + mesh_numVertices.append(mesh.num_vertices()) + + c0 = c.copy(deepcopy=True) # Solve the advection equation over the subinterval (t0, t1] c = run_simulation(mesh, t0, t1, c) -cpu_time_metric_advection = time.time() - cpu_time_metric_advection -print(f"CPU time with metric advection: {cpu_time_metric_advection:.2f} s") - -fig, axes = plt.subplots(1, 2, figsize=(8, 5), sharey=True) -triplot(mesh, axes=axes[0]) -tripcolor(c, axes=axes[1]) - -for ax in axes.flat: - ax.set_xlim(0, 1) - ax.set_ylim(0, 1) - ax.set_aspect("equal") - ax.set_xlabel("x") - ax.set_ylabel("y") -fig.tight_layout() -fig.savefig("bubble_shear-metric_advection.jpg", dpi=300, bbox_inches="tight") + plot_mesh(mesh, c0, c, i, "metric_advection") -# Redefine the initial condition on the final adapted mesh c0 = get_initial_condition(mesh) metric_adv_c0 = Function(FunctionSpace(mesh, "CG", 1)).interpolate(c0) metric_adv_error = compute_rel_error(metric_adv_c0, c) -print(f"Relative L2 error with metric advection: {metric_adv_error:.2f}%") +print( + f"Metric advection mesh adaptation.\n" + f" Avg. number of vertices: {np.average(mesh_numVertices):.1f}\n" + f" Relative L2 error: {metric_adv_error:.2f}%" +) # .. code-block:: console # -# Relative L2 error with metric advection: 29.31% +# Metric advection mesh adaptation. +# Avg. number of vertices: 1932.5 +# Relative L2 error: 31.07% # -# .. figure:: bubble_shear-metric_advection.jpg +# .. figure:: bubble_shear-metric_advection_14.jpg # :figwidth: 80% # :align: center - -fig, ax = plt.subplots() -errors = [coarse_error, fine_error, classical_error, metric_adv_error] -cpu_times = [ - cpu_time_coarse, - cpu_time_fine, - cpu_time_classical, - cpu_time_metric_advection, -] -ax.plot(cpu_times, errors, "o") -# label each point -for i, txt in enumerate(["Coarse mesh", "Fine mesh", "Classical", "Metric advection"]): - ax.annotate( - txt, - (cpu_times[i], errors[i]), - textcoords="offset points", - xytext=(0, 10), - ha="center", - ) -ax.set_xlabel("CPU time (s)") -ax.set_ylabel("Relative L2 error (%)") -fig.savefig("bubble_shear-time_error.jpg", dpi=300, bbox_inches="tight") - -# .. figure:: bubble_shear-time_error.jpg +# +# The relative L2 error of 31.07% is slightly lower than the one obtained with the +# classical mesh adaptation algorithm, but note that the average number of vertices is +# about 20% smaller. Looking into the final adapted mesh and concentration fields in +# the above figure, we now observe that the mesh is indeed refined in a much wider +# region - ensuring that the bubble remains well-resolved throughout the subinterval. +# However, this also means that the available resolution is more widely distributed, +# leading to a coarser local resolution compared to classical mesh adaptation. +# +# Let us now consider the implications and limitations of each approach. Firstly, +# given a high enough adaptation frequency (i.e. number of adaptations/subintervals), +# the lagging mesh problem can be mitigated with the classical mesh adaptation +# algorithm. This may end up being more computationally efficient than the metric +# advection algorithm, which requires solving the advection equation four times in total +# at each timestep, as well as potentially expensive metric computations. +# However, while increasing the mesh adaptation frequency would alleviate the mesh lag, +# doing so may introduce large errors due to frequent solution transfers between meshes. +# +# This is where the advantage of the metric advection algorithm lies: it predicts where +# to prescribe fine resolution in the future, and thus avoids the need for frequent +# solution transfers. We can assure ourselves of that by repeating the above simulations +# with `num_adaptations = 5`, which yields relative L2 errors of 61.06% and 42.23% for +# the classical and metric advection algorithms, respectively. Furthermore, the problem +# considered in this example is relatively well-suited for classical mesh adaptation, +# as the bubble concentration field reverses and therefore often indeed remains in the +# finely-resolved region. We can observe that in the below figure, at the subinterval +# :math:`(1.4 s, 1.6 s]`. +# +# .. figure:: bubble_shear-classical_7.jpg # :figwidth: 80% # :align: center +# +# In conclusion, the choice of mesh adaptation algorithm depends on the specific problem +# at hand, as well as the computational resources available. We encourage the reader to +# experiment with different metric parameters, different adaptation frequencies, and +# even different velocity fields to further explore the capabilities and limitations of +# the above-presented mesh adaptation algorithms. +# +# This demo can also be accessed as a `Python script `__. diff --git a/demos/demo_references.bib b/demos/demo_references.bib index 3c4d091..be1e854 100644 --- a/demos/demo_references.bib +++ b/demos/demo_references.bib @@ -8,3 +8,17 @@ @article{FPP+:2009 year={2009}, publisher={Elsevier} } + +@article{Barral:2016, + title={Anisotropic mesh adaptation in Firedrake with PETSc DMPlex}, + author={Barral, Nicolas and Knepley, Matthew G and Lange, Michael and Piggott, Matthew D and Gorman, Gerard J}, + journal={arXiv preprint arXiv:1610.09874}, + year={2016} +} + +@phdthesis{Wilson:2010, + title={Modelling multiple-material flows on adaptive unstructured meshes}, + author={Wilson, Cian}, + year={2010}, + school={Imperial College London} +} \ No newline at end of file From 5510e9551d740eb3d8ccca6422841fa81a9d73eb Mon Sep 17 00:00:00 2001 From: ddundo Date: Thu, 12 Sep 2024 11:35:12 +0000 Subject: [PATCH 03/12] #55: Shorten run time of bubble shear demo --- test/test_demos.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/test/test_demos.py b/test/test_demos.py index 5375e1d..972fed2 100644 --- a/test/test_demos.py +++ b/test/test_demos.py @@ -21,6 +21,10 @@ # -- The value is the replacement string (use "" to remove the original code) modifications = { "ping_pong.py": {"niter = 50": "niter = 1"}, + "bubble_shear.py": { + "simulation_end_time = T / 2.0": "simulation_end_time = T / 10.0", + "num_adaptations = 15": "num_adaptations = 2", + }, } From e652186f3837642a84c77e4ae34e6789af98eaf6 Mon Sep 17 00:00:00 2001 From: ddundo Date: Sun, 15 Sep 2024 11:18:25 +0000 Subject: [PATCH 04/12] #144: Solve for the entire metric --- animate/metric.py | 12 ++++----- demos/bubble_shear.py | 58 +++++++++++++++++++------------------------ 2 files changed, 31 insertions(+), 39 deletions(-) diff --git a/animate/metric.py b/animate/metric.py index 141e9fa..cbcc4e7 100644 --- a/animate/metric.py +++ b/animate/metric.py @@ -84,12 +84,12 @@ def __init__(self, function_space, *args, **kwargs): if isinstance(fs.dof_count, Iterable): raise ValueError("Riemannian metric cannot be built in a mixed space.") self._check_space() - # rank = len(fs.dof_dset.dim) - # if rank != 2: - # raise ValueError( - # "Riemannian metric should be matrix-valued," - # f" not rank-{rank} tensor-valued." - # ) + rank = len(fs.dof_dset.dim) + if rank != 2: + raise ValueError( + "Riemannian metric should be matrix-valued," + f" not rank-{rank} tensor-valued." + ) # Stash mesh data plex = mesh.topology_dm.clone() diff --git a/demos/bubble_shear.py b/demos/bubble_shear.py index 8ac8fd8..df3dbb7 100644 --- a/demos/bubble_shear.py +++ b/demos/bubble_shear.py @@ -209,7 +209,7 @@ def compute_rel_error(c_init, c_final): # For the purposes of this demo, we are going to adapt the mesh 15 times throughout the # simulation, at equal time intervals (i.e., every 0.2s of simulation time). :: -num_adaptations = 15 +num_adaptations = 5 interval_length = simulation_end_time / num_adaptations # We will also define a function that will allow us to easily plot the adapted mesh, as @@ -339,20 +339,20 @@ def adapt_classical(mesh, c): def adapt_metric_advection(mesh, t_start, t_end, c): P1_ten = TensorFunctionSpace(mesh, "CG", 1) - metric = RiemannianMetric(P1_ten) - metric_ = RiemannianMetric(P1_ten) + m = RiemannianMetric(P1_ten) # metric at current timestep + m_ = RiemannianMetric(P1_ten) # metric at previous timestep metric_intersect = RiemannianMetric(P1_ten) # Compute the Hessian metric at t_start - for mtrc in [metric, metric_, metric_intersect]: + for mtrc in [m, m_, metric_intersect]: mtrc.set_parameters(metric_params) - metric_.compute_hessian(c) - metric_.enforce_spd(restrict_sizes=False, restrict_anisotropy=False) + m_.compute_hessian(c) + m_.enforce_spd(restrict_sizes=True, restrict_anisotropy=True) - Q = FunctionSpace(mesh, "CG", 1) - h_bc = Function(Q) # used to set the boundary condition below - m_ = Function(Q) # metric component at previous timestep - m = Function(Q) # metric component at current timestep + # Set the boundary condition for the metric tensor + h_bc = Function(P1_ten) + h_max = metric_params["dm_plex_metric"]["h_max"] + h_bc.interpolate(Constant([[1.0 / h_max**2, 0.0], [0.0, 1.0 / h_max**2]])) V = VectorFunctionSpace(mesh, "CG", 1) u = Function(V) @@ -370,14 +370,15 @@ def adapt_metric_advection(mesh, t_start, t_end, c): tau = min_value(tau, U * h / (6 * D)) # Apply SUPG stabilisation - phi = TestFunction(Q) + phi = TestFunction(P1_ten) phi += tau * dot(u, grad(phi)) - # Variational form of the advection equation, as before - trial = TrialFunction(Q) + # Variational form of the advection equation for the metric tensor + trial = TrialFunction(P1_ten) a = inner(trial, phi) * dx + dt * theta * inner(dot(u, grad(trial)), phi) * dx L = inner(m_, phi) * dx - dt * (1 - theta) * inner(dot(u_, grad(m_)), phi) * dx - lvp = LinearVariationalProblem(a, L, m, bcs=DirichletBC(Q, h_bc, "on_boundary")) + bcs = DirichletBC(P1_ten, h_bc, "on_boundary") + lvp = LinearVariationalProblem(a, L, m, bcs=bcs) lvs = LinearVariationalSolver(lvp) # Integrate from t_start to t_end @@ -385,23 +386,14 @@ def adapt_metric_advection(mesh, t_start, t_end, c): while t < t_end + 0.5 * float(dt): u.interpolate(velocity_expression(mesh, t)) - # Advect each metric component M_ij individually - for i, j in [(0, 0), (0, 1), (1, 1)]: - # print(i, j) - h_max = metric_params["dm_plex_metric"]["h_max"] - h_bc.assign(1.0 / h_max**2 if i == j else 0.0) - m_.assign(metric_.sub(2 * i + j)) - lvs.solve() - metric.sub(2 * i + j).assign(m) - # Metric is symmetric so we can copy the other off-diagonal component - metric.sub(2).assign(metric.sub(1)) + lvs.solve() # Intersect metrics at every timestep - metric.enforce_spd(restrict_sizes=True, restrict_anisotropy=True) - metric_intersect.intersect(metric) + m.enforce_spd(restrict_sizes=True, restrict_anisotropy=True) + metric_intersect.intersect(m) # Update fields at the previous timestep - metric_.assign(metric) + m_.assign(m) u_.assign(u) t += float(dt) @@ -444,16 +436,16 @@ def adapt_metric_advection(mesh, t_start, t_end, c): # .. code-block:: console # # Metric advection mesh adaptation. -# Avg. number of vertices: 1932.5 -# Relative L2 error: 31.07% +# Avg. number of vertices: 2032.7 +# Relative L2 error: 31.30% # # .. figure:: bubble_shear-metric_advection_14.jpg # :figwidth: 80% # :align: center # -# The relative L2 error of 31.07% is slightly lower than the one obtained with the +# The relative L2 error of 31.00% is slightly lower than the one obtained with the # classical mesh adaptation algorithm, but note that the average number of vertices is -# about 20% smaller. Looking into the final adapted mesh and concentration fields in +# about 15% smaller. Looking into the final adapted mesh and concentration fields in # the above figure, we now observe that the mesh is indeed refined in a much wider # region - ensuring that the bubble remains well-resolved throughout the subinterval. # However, this also means that the available resolution is more widely distributed, @@ -463,7 +455,7 @@ def adapt_metric_advection(mesh, t_start, t_end, c): # given a high enough adaptation frequency (i.e. number of adaptations/subintervals), # the lagging mesh problem can be mitigated with the classical mesh adaptation # algorithm. This may end up being more computationally efficient than the metric -# advection algorithm, which requires solving the advection equation four times in total +# advection algorithm, which requires solving the advection equation once more # at each timestep, as well as potentially expensive metric computations. # However, while increasing the mesh adaptation frequency would alleviate the mesh lag, # doing so may introduce large errors due to frequent solution transfers between meshes. @@ -471,7 +463,7 @@ def adapt_metric_advection(mesh, t_start, t_end, c): # This is where the advantage of the metric advection algorithm lies: it predicts where # to prescribe fine resolution in the future, and thus avoids the need for frequent # solution transfers. We can assure ourselves of that by repeating the above simulations -# with `num_adaptations = 5`, which yields relative L2 errors of 61.06% and 42.23% for +# with `num_adaptations = 5`, which yields relative L2 errors of 61.06% and 36.86% for # the classical and metric advection algorithms, respectively. Furthermore, the problem # considered in this example is relatively well-suited for classical mesh adaptation, # as the bubble concentration field reverses and therefore often indeed remains in the From 69d8404d389be8b1c0899631dffa1e2a25f5b776 Mon Sep 17 00:00:00 2001 From: ddundo Date: Sun, 15 Sep 2024 11:54:22 +0000 Subject: [PATCH 05/12] Fix changed-files pathing --- .github/workflows/test_suite.yml | 4 ++-- demos/bubble_shear.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/test_suite.yml b/.github/workflows/test_suite.yml index 5a3e636..f4e1ade 100644 --- a/.github/workflows/test_suite.yml +++ b/.github/workflows/test_suite.yml @@ -39,8 +39,8 @@ jobs: with: files: | .github/workflows/test_suite.yml - *.py - *.cxx + **/*.py + **/*.cxx - name: 'Cleanup' if: ${{ always() }} diff --git a/demos/bubble_shear.py b/demos/bubble_shear.py index df3dbb7..cb567f9 100644 --- a/demos/bubble_shear.py +++ b/demos/bubble_shear.py @@ -209,7 +209,7 @@ def compute_rel_error(c_init, c_final): # For the purposes of this demo, we are going to adapt the mesh 15 times throughout the # simulation, at equal time intervals (i.e., every 0.2s of simulation time). :: -num_adaptations = 5 +num_adaptations = 15 interval_length = simulation_end_time / num_adaptations # We will also define a function that will allow us to easily plot the adapted mesh, as From 118de05d0291868951e7469d53a715764e30dc5f Mon Sep 17 00:00:00 2001 From: ddundo Date: Mon, 16 Sep 2024 12:36:10 +0000 Subject: [PATCH 06/12] #55: Various minor fixes --- demos/bubble_shear.py | 55 ++++++++++++++++++++----------------------- 1 file changed, 25 insertions(+), 30 deletions(-) diff --git a/demos/bubble_shear.py b/demos/bubble_shear.py index cb567f9..7eb27ac 100644 --- a/demos/bubble_shear.py +++ b/demos/bubble_shear.py @@ -1,8 +1,8 @@ # On-the-fly time-dependent mesh adaptation # ######################################### # -# In this demo we consider the 2-dimensional version of the mesh adaptation experiment -# presented in :cite:. The problem comprises a bubble of tracer +# In this demo we consider the two-dimensional version of the mesh adaptation experiment +# presented in :cite:`Barral:2016``. The problem comprises a bubble of tracer # concentration field advected by a time-varying flow. # We will consider two different mesh adaptation strategies: the classical mesh # adaptation algorithm, which adapts the mesh several times throughout the simulation @@ -33,9 +33,8 @@ # .. math:: # \mathbf{u}(x, y, t) := \left(2\sin^2(\pi x)\sin(2\pi y), -\sin(2\pi x)\sin^2(\pi y) \right) \cos(2\pi t/T), # -# where :math:`T` is the period. Note that the velocity field is solenoidal. At each -# timestep of the simulation we will update this field so we define a function that will -# return its vector expression. :: +# where :math:`T` is the period. At each timestep of the simulation we will update this +# field so we define a function that will return its vector expression. :: from firedrake import * @@ -66,7 +65,6 @@ def get_initial_condition(mesh): r = sqrt(pow(x - ball_x0, 2) + pow(y - ball_y0, 2)) c0 = Function(FunctionSpace(mesh, "CG", 1)) c0.interpolate(conditional(r < ball_r, 1.0, 0.0)) - # return conditional(r < ball_r0, 1.0, 0.0) return c0 @@ -77,13 +75,9 @@ def get_initial_condition(mesh): # the initial condition :math:`c_0 = c(x, y, t_{\text{start}})`. The function returns # the solution :math:`c(x, y, t_{\text{end}})`. # -# Note that we must include streamline upwind Petrov Galerkin (SUPG) stabilisation in +# Note that we include streamline upwind Petrov Galerkin (SUPG) stabilisation in # the test function in order to ensure numerical stability. :: -# Time-stepping parameters -# dt = 0.01 # timestep size -# theta = 0.5 # Crank-Nicolson implicitness - def run_simulation(mesh, t_start, t_end, c0): Q = FunctionSpace(mesh, "CG", 1) @@ -183,8 +177,9 @@ def run_simulation(mesh, t_start, t_end, c0): # is, the more diffusion is added. We encourage the reader to verify this by running the # simulation on a sequence of finer uniform meshes. # -# In order to quantify the above observation, we will compute the relative L2 error -# between the initial condition and the final concentration field on the final mesh. :: +# In order to quantify the above observation, we will compute the relative :math:`L^2` +# error between the initial condition and the final concentration field on the final +# mesh. :: def compute_rel_error(c_init, c_final): @@ -235,7 +230,7 @@ def plot_mesh(mesh, c0, c1, i, method): # As mentioned in the introduction, we shall demonstrate two different mesh adaptation # strategies. Let us begin with the classical mesh adaptation algorithm, which adapts # each mesh before solving the advection equation over the corresponding subinterval. -# Here we will use the :class:`RiemannianMetric` class to define the metric, which we +# Here we will use the :class:`~.RiemannianMetric` class to define the metric, which we # will compute based on the Hessian of the concentration field. Let us therefore define # a function which takes the original mesh and the concentration field as arguments, and # returns the adapted mesh. We will also define parameters for computing the metric. :: @@ -287,7 +282,7 @@ def adapt_classical(mesh, c): plot_mesh(mesh, c0, c, i, "classical") # Now let us examine the final adapted mesh and final concentration field computed on it. -# We will also compute the relative L2 error. :: +# We will also compute the relative :math:`L^2` error. :: # Redefine the initial condition on the final adapted mesh c0 = get_initial_condition(mesh) @@ -310,11 +305,11 @@ def adapt_classical(mesh, c): # :figwidth: 80% # :align: center # -# As we can see, the relative L2 error is lower than the one obtained with the uniform -# fine mesh, even though the average number of vertices is more than 8 times smaller. -# This demonstrates the effectiveness of mesh adaptation, which we can also observe in -# the figure above. We see that final mesh has clearly been refined around the bubble -# at the beginning of the final subinterval and coarsened elsewhere. +# As we can see, the relative :math:`L^2` error is lower than the one obtained with the +# uniform fine mesh, even though the average number of vertices is more than 8 times +# smaller. This demonstrates the effectiveness of mesh adaptation, which we can also +# observe in the figure above. We see that final mesh has clearly been refined around +# the bubble at the beginning of the final subinterval and coarsened elsewhere. # # However, we also observe that the bubble has advected out of the fine resolution # region by the end of the subinterval. This is a common occurence in time-dependent @@ -328,7 +323,7 @@ def adapt_classical(mesh, c): # prescribe fine resolution in regions where it is not needed. # # For advection-dominated problems, as is the case here, a *metric advection* algorithm -# has been proposed in :cite:. The idea is to still compute the metric +# has been proposed in :cite:`Wilson:2010``. The idea is to still compute the metric # based on the solution at the current time, but then to advect the metric along the # flow in order to predict where to prescribe fine resolution in the future. By # combining the advected metrics in time, we obtain a final metric that is @@ -443,9 +438,9 @@ def adapt_metric_advection(mesh, t_start, t_end, c): # :figwidth: 80% # :align: center # -# The relative L2 error of 31.00% is slightly lower than the one obtained with the -# classical mesh adaptation algorithm, but note that the average number of vertices is -# about 15% smaller. Looking into the final adapted mesh and concentration fields in +# The relative :math:`L^2` error of 31.30% is slightly lower than the one obtained with +# the classical mesh adaptation algorithm, but note that the average number of vertices +# is about 15% smaller. Looking into the final adapted mesh and concentration fields in # the above figure, we now observe that the mesh is indeed refined in a much wider # region - ensuring that the bubble remains well-resolved throughout the subinterval. # However, this also means that the available resolution is more widely distributed, @@ -463,12 +458,12 @@ def adapt_metric_advection(mesh, t_start, t_end, c): # This is where the advantage of the metric advection algorithm lies: it predicts where # to prescribe fine resolution in the future, and thus avoids the need for frequent # solution transfers. We can assure ourselves of that by repeating the above simulations -# with `num_adaptations = 5`, which yields relative L2 errors of 61.06% and 36.86% for -# the classical and metric advection algorithms, respectively. Furthermore, the problem -# considered in this example is relatively well-suited for classical mesh adaptation, -# as the bubble concentration field reverses and therefore often indeed remains in the -# finely-resolved region. We can observe that in the below figure, at the subinterval -# :math:`(1.4 s, 1.6 s]`. +# with ``num_adaptations = 5``, which yields relative the errors of 61.06% and 36.86% +# for the classical and metric advection algorithms, respectively. Furthermore, the +# problem considered in this example is relatively well-suited for classical mesh +# adaptation, as the bubble concentration field reverses and therefore often indeed +# remains in the finely-resolved region. We can observe that in the below figure, at the +# subinterval :math:`(1.4 s, 1.6 s]`. # # .. figure:: bubble_shear-classical_7.jpg # :figwidth: 80% From 8af5e03cea6bafbdcfb62616f8ece41431c3a64f Mon Sep 17 00:00:00 2001 From: ddundo Date: Mon, 16 Sep 2024 12:46:50 +0000 Subject: [PATCH 07/12] #55: Expand intro on FPI --- demos/bubble_shear.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/demos/bubble_shear.py b/demos/bubble_shear.py index 7eb27ac..e589c66 100644 --- a/demos/bubble_shear.py +++ b/demos/bubble_shear.py @@ -11,7 +11,11 @@ # prescribe fine resolution in the future. Both algorithms are an example of # *on-the-fly* time-dependent mesh adaptation algorithms, where the mesh is adapted # before each subinterval of the simulation, as opposed to fixed-point iteration -# algorithms, where the mesh is iteratively adapted at the end of the simulation. +# algorithms. +# Fixed-point iteration algorithms solve the equation of interest at each iteration, +# before adapting the mesh based on computed solutions at the end of each iteration. +# This clearly makes them effective at predicting where to prescribe fine resolution +# throughout the simulation, but imply substantial computational cost. # # We begin by defining the advection problem. We consider the advection equation # From 414a6a90769add123a14510d67ae217efc7335b0 Mon Sep 17 00:00:00 2001 From: ddundo Date: Mon, 16 Sep 2024 12:50:52 +0000 Subject: [PATCH 08/12] #55: Add exercise rubric --- demos/bubble_shear.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/demos/bubble_shear.py b/demos/bubble_shear.py index e589c66..db4371c 100644 --- a/demos/bubble_shear.py +++ b/demos/bubble_shear.py @@ -178,8 +178,7 @@ def run_simulation(mesh, t_start, t_end, c0): # diffusion. Numerical diffusion is necessary for numerical stability and for preventing # oscillations, but it also makes the solution irreversible. The amount of difussion # added is related to the grid Péclet number :math:`Pe = U\,h/2D`: the coarser the mesh -# is, the more diffusion is added. We encourage the reader to verify this by running the -# simulation on a sequence of finer uniform meshes. +# is, the more diffusion is added. # # In order to quantify the above observation, we will compute the relative :math:`L^2` # error between the initial condition and the final concentration field on the final @@ -474,9 +473,12 @@ def adapt_metric_advection(mesh, t_start, t_end, c): # :align: center # # In conclusion, the choice of mesh adaptation algorithm depends on the specific problem -# at hand, as well as the computational resources available. We encourage the reader to -# experiment with different metric parameters, different adaptation frequencies, and -# even different velocity fields to further explore the capabilities and limitations of -# the above-presented mesh adaptation algorithms. +# at hand, as well as the computational resources available. +# +# .. rubric:: Exercise +# +# We encourage the reader to experiment with different metric parameters, different +# adaptation frequencies, and even different velocity fields to further explore the +# capabilities and limitations of the above-presented mesh adaptation algorithms. # # This demo can also be accessed as a `Python script `__. From 4196b34000f439387cbac45e60fb82103a98e2b0 Mon Sep 17 00:00:00 2001 From: ddundo Date: Mon, 16 Sep 2024 13:55:44 +0000 Subject: [PATCH 09/12] #55: Expand num adaptations discussion --- demos/bubble_shear.py | 31 ++++++++++++++++++++++++++----- 1 file changed, 26 insertions(+), 5 deletions(-) diff --git a/demos/bubble_shear.py b/demos/bubble_shear.py index db4371c..8cafa9e 100644 --- a/demos/bubble_shear.py +++ b/demos/bubble_shear.py @@ -461,12 +461,30 @@ def adapt_metric_advection(mesh, t_start, t_end, c): # This is where the advantage of the metric advection algorithm lies: it predicts where # to prescribe fine resolution in the future, and thus avoids the need for frequent # solution transfers. We can assure ourselves of that by repeating the above simulations -# with ``num_adaptations = 5``, which yields relative the errors of 61.06% and 36.86% -# for the classical and metric advection algorithms, respectively. Furthermore, the +# with ``num_adaptations = 5``, which yields relative errors of 61.06% and 36.86% +# for the classical and metric advection algorithms, respectively. Conversely, +# increasing the adaptation frequency to ``num_adaptations = 50`` yields again relative +# errors closer to one another. Note that the algorithms are identical if we adapt at +# every timestep. We summarise these results in the table below, noting also the average +# number of vertices, :math:`N_v`. +# +# .. table:: +# +# ======================= ============================== ===================================== +# Number of adaptations Classical (avg. :math:`N_v`) Metric advection (avg. :math:`N_v`) +# ======================= ============================== ===================================== +# 5 61.06% (2200.2) 36.86% (1922.2) +# 15 31.76% (2302.3) 31.30% (2032.7) +# 50 26.80% (2522.1) 27.99% (2166.3) +# ======================= ============================== ===================================== +# +# Furthermore, the # problem considered in this example is relatively well-suited for classical mesh # adaptation, as the bubble concentration field reverses and therefore often indeed # remains in the finely-resolved region. We can observe that in the below figure, at the -# subinterval :math:`(1.4 s, 1.6 s]`. +# subinterval :math:`(1.4 s, 1.6 s]`. This also means that the meshes adapted using +# classical and metric advection algorithms are qualitatively similar at this +# subinterval. # # .. figure:: bubble_shear-classical_7.jpg # :figwidth: 80% @@ -477,8 +495,11 @@ def adapt_metric_advection(mesh, t_start, t_end, c): # # .. rubric:: Exercise # -# We encourage the reader to experiment with different metric parameters, different +# We encourage the reader to repeat above experiments and investigate each of the +# adapted meshes. At what subintervals do the two algorithms produce most similar and +# most different meshes? +# We further encourage experimentation with different metric parameters, different # adaptation frequencies, and even different velocity fields to further explore the -# capabilities and limitations of the above-presented mesh adaptation algorithms. +# capabilities and limitations of the above-presented algorithms. # # This demo can also be accessed as a `Python script `__. From 789bb67c0947ede49253d35a74e3afa54f1bb038 Mon Sep 17 00:00:00 2001 From: ddundo Date: Mon, 16 Sep 2024 14:24:05 +0000 Subject: [PATCH 10/12] #55: Add IC to uniform fig --- demos/bubble_shear.py | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/demos/bubble_shear.py b/demos/bubble_shear.py index 8cafa9e..098b987 100644 --- a/demos/bubble_shear.py +++ b/demos/bubble_shear.py @@ -150,20 +150,24 @@ def run_simulation(mesh, t_start, t_end, c0): c0_fine = get_initial_condition(mesh_fine) c_fine_final = run_simulation(mesh_fine, 0.0, simulation_end_time, c0_fine) -# We can now visualise the final concentration fields on the two meshes. :: +# We can now compare computed final concentration fields on the two meshes to the +# initial condition. :: import matplotlib.pyplot as plt from firedrake.pyplot import * -fig, axes = plt.subplots(1, 2, figsize=(8, 5), sharey=True) -for ax, c, mesh in zip(axes, [c_coarse_final, c_fine_final], [mesh_coarse, mesh_fine]): - tripcolor(c, axes=ax) - ax.set_title(f"{mesh.num_vertices()} mesh vertices") - ax.set_xlim(0, 1) - ax.set_ylim(0, 1) - ax.set_aspect("equal") - ax.set_xlabel("x") - ax.set_ylabel("y") +fig, axes = plt.subplots(2, 2, figsize=(8, 8), sharex=True, sharey=True) +time_labels = [r"$t=0$", r"$t=T/2$"] +c_values = [[c0_coarse, c0_fine], [c_coarse_final, c_fine_final]] +meshes = [mesh_coarse, mesh_fine] +for i in range(2): + for ax, c, mesh in zip(axes[i], c_values[i], meshes): + im = tripcolor(c, axes=ax) + ax.set_title(f"{mesh.num_vertices()} mesh vertices") + ax.text(0.05, 0.05, time_labels[i], fontsize=12, color="white", ha="left") + ax.set_xlim(0, 1) + ax.set_ylim(0, 1) + ax.set_aspect("equal") fig.tight_layout() fig.savefig("bubble_shear-uniform.jpg", dpi=300, bbox_inches="tight") From 6111ca59ad7687dada97c42d8b839f6ad0ddccfe Mon Sep 17 00:00:00 2001 From: ddundo Date: Mon, 16 Sep 2024 14:25:57 +0000 Subject: [PATCH 11/12] #55: Fix citation --- demos/bubble_shear.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/demos/bubble_shear.py b/demos/bubble_shear.py index 098b987..c660b2a 100644 --- a/demos/bubble_shear.py +++ b/demos/bubble_shear.py @@ -2,7 +2,7 @@ # ######################################### # # In this demo we consider the two-dimensional version of the mesh adaptation experiment -# presented in :cite:`Barral:2016``. The problem comprises a bubble of tracer +# presented in :cite:`Barral:2016`. The problem comprises a bubble of tracer # concentration field advected by a time-varying flow. # We will consider two different mesh adaptation strategies: the classical mesh # adaptation algorithm, which adapts the mesh several times throughout the simulation @@ -330,7 +330,7 @@ def adapt_classical(mesh, c): # prescribe fine resolution in regions where it is not needed. # # For advection-dominated problems, as is the case here, a *metric advection* algorithm -# has been proposed in :cite:`Wilson:2010``. The idea is to still compute the metric +# has been proposed in :cite:`Wilson:2010`. The idea is to still compute the metric # based on the solution at the current time, but then to advect the metric along the # flow in order to predict where to prescribe fine resolution in the future. By # combining the advected metrics in time, we obtain a final metric that is From 6fff07a7f296fc5b1530a3ad84b68ce2bb793cf5 Mon Sep 17 00:00:00 2001 From: ddundo Date: Tue, 17 Sep 2024 07:26:28 +0000 Subject: [PATCH 12/12] #55: Minor comment fixes --- demos/bubble_shear.py | 33 +++++++++++++++------------------ 1 file changed, 15 insertions(+), 18 deletions(-) diff --git a/demos/bubble_shear.py b/demos/bubble_shear.py index c660b2a..76836c8 100644 --- a/demos/bubble_shear.py +++ b/demos/bubble_shear.py @@ -11,11 +11,10 @@ # prescribe fine resolution in the future. Both algorithms are an example of # *on-the-fly* time-dependent mesh adaptation algorithms, where the mesh is adapted # before each subinterval of the simulation, as opposed to fixed-point iteration -# algorithms. -# Fixed-point iteration algorithms solve the equation of interest at each iteration, -# before adapting the mesh based on computed solutions at the end of each iteration. -# This clearly makes them effective at predicting where to prescribe fine resolution -# throughout the simulation, but imply substantial computational cost. +# algorithms. Fixed-point iteration algorithms solve the equation of interest at each +# iteration, before adapting the mesh based on computed solutions at the end of each +# iteration. This clearly makes them effective at predicting where to prescribe fine +# resolution throughout the simulation, but implies substantial computational cost. # # We begin by defining the advection problem. We consider the advection equation # @@ -482,13 +481,12 @@ def adapt_metric_advection(mesh, t_start, t_end, c): # 50 26.80% (2522.1) 27.99% (2166.3) # ======================= ============================== ===================================== # -# Furthermore, the -# problem considered in this example is relatively well-suited for classical mesh -# adaptation, as the bubble concentration field reverses and therefore often indeed -# remains in the finely-resolved region. We can observe that in the below figure, at the -# subinterval :math:`(1.4 s, 1.6 s]`. This also means that the meshes adapted using -# classical and metric advection algorithms are qualitatively similar at this -# subinterval. +# Furthermore, the problem considered in this example is relatively well-suited for +# classical mesh adaptation, as the bubble concentration field reverses and therefore +# often indeed remains in the finely-resolved region. We can observe that in the below +# figure, at the subinterval :math:`(1.4 s, 1.6 s]`. This also means that the meshes +# adapted using classical and metric advection algorithms are qualitatively similar at +# this subinterval. # # .. figure:: bubble_shear-classical_7.jpg # :figwidth: 80% @@ -499,11 +497,10 @@ def adapt_metric_advection(mesh, t_start, t_end, c): # # .. rubric:: Exercise # -# We encourage the reader to repeat above experiments and investigate each of the -# adapted meshes. At what subintervals do the two algorithms produce most similar and -# most different meshes? -# We further encourage experimentation with different metric parameters, different -# adaptation frequencies, and even different velocity fields to further explore the -# capabilities and limitations of the above-presented algorithms. +# Repeat above experiments and investigate each of the adapted meshes. At what +# subintervals do the two algorithms produce most similar and most different meshes? +# Experiment with different metric parameters, different adaptation frequencies, and +# even different velocity fields to further explore the capabilities and limitations of +# the algorithms presented above. # # This demo can also be accessed as a `Python script `__.