From 7ee83312c69cef099fe857dea1d77e139dd47c16 Mon Sep 17 00:00:00 2001 From: ddundo Date: Tue, 12 Mar 2024 00:04:33 +0000 Subject: [PATCH 01/10] #73: Bubble shear demo --- demos/bubble_shear.py | 232 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 232 insertions(+) create mode 100644 demos/bubble_shear.py diff --git a/demos/bubble_shear.py b/demos/bubble_shear.py new file mode 100644 index 00000000..b5f5197c --- /dev/null +++ b/demos/bubble_shear.py @@ -0,0 +1,232 @@ +# Advection of a tracer bubble in a non-uniform shear flow +# =================================================== + +# In this demo, we consider the advection of a scalar concentration field, similarly to +# our previous solid body rotation demo. However, we now consider a background velocity +# field that is not uniform in time. We further prescribe homogeneous Dirichlet boundary +# conditions for the concentration, so in this demo we solve +# +# .. 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. +# +# First, we import Firedrake and Goalie. :: + +from firedrake import * +from goalie import * + +# We begin by definining the background velocity field :math:`\mathbf{u}(x, y, t)`. +# Specifically, we choose the velocity field that is periodic in time, such as +# in TODO: cite, and is given by +# +# .. math:: +# \mathbf{u}(x, y, t) := \cos(2\pi t/T)\left(2\sin^2(\pi x)\sin(2\pi y), -\sin(2\pi x)\sin^2(\pi y) \right), +# +# where :math:`T` is the period. At each timestep of the simulation we will have to +# update this field so we define a function that will return this vector expression. :: + +period = 6.0 + + +def velocity_expression(x, y, t): + u_expr = as_vector( + [ + 2 * sin(pi * x) ** 2 * sin(2 * pi * y) * cos(2 * pi * t / period), + -sin(2 * pi * x) * sin(pi * y) ** 2 * cos(2 * pi * t / period), + ] + ) + return u_expr + + +# However, since the background velocity :math:`\mathbf{u}` is already known, we do not +# solve for it and pass it as a field to :class:`MeshSeq`. Instead, we pass only the +# tracer concentration :math:`c` and, for now, avoid specifying how to handle the +# velocity field. Therefore, at this stage we only define how to build the +# :class:`FunctionSpace` for the tracer concentration. :: + +fields = ["c"] + + +def get_function_spaces(mesh): + return {"c": FunctionSpace(mesh, "CG", 1)} + + +# We proceed similarly with prescribing initial and boundary 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. :: + + +def get_bcs(mesh_seq): + def bcs(index): + return [DirichletBC(mesh_seq.function_spaces["c"][index], 0.0, "on_boundary")] + + return bcs + + +def ball_initial_condition(x, y): + ball_r0, ball_x0, ball_y0 = 0.15, 0.5, 0.65 + r = sqrt(pow(x - ball_x0, 2) + pow(y - ball_y0, 2)) + return conditional(r < ball_r0, 1.0, 0.0) + + +def get_initial_condition(mesh_seq): + fs = mesh_seq.function_spaces["c"][0] + x, y = SpatialCoordinate(mesh_seq[0]) + ball = ball_initial_condition(x, y) + c0 = assemble(interpolate(ball, fs)) + return {"c": c0} + + +# To compile the weak form, we require both the concentration field :math:`c` and the +# background velocity field :math:`\mathbf{u}` in :meth:`get_form`. In this demo, we +# will always pass both these fields from :meth:`get_solver`. Note that we still do not +# define the velocity field. +# As in the point discharge with diffusion demo, we include additional `streamline +# upwind Petrov Galerkin (SUPG)` stabilisation by modifying the test function :math:`\psi`. +# To advance in time, we use Crank-Nicolson discretisation. :: + + +def get_form(mesh_seq): + def form(index, form_fields): + Q = mesh_seq.function_spaces["c"][index] + c, c_ = form_fields["c"] + u, u_ = form_fields["u"] + + R = FunctionSpace(mesh_seq[index], "R", 0) + dt = Function(R).assign(mesh_seq.time_partition.timesteps[index]) + theta = Function(R).assign(0.5) # Crank-Nicolson theta + + # SUPG stabilisation parameter + D = Function(R).assign(0.1) # diffusivity coefficient + h = CellSize(mesh_seq[index]) # mesh cell size + U = sqrt(dot(u, u)) # velocity magnitude + tau = 0.5 * h / U + tau = min_value(tau, U * h / (6 * D)) + + phi = TestFunction(Q) + phi += tau * dot(u, grad(phi)) + + a = c * phi * dx + dt * theta * dot(u, grad(c)) * phi * dx + L = c_ * phi * dx - dt * (1 - theta) * dot(u_, grad(c_)) * phi * dx + F = a - L + + return {"c": F} + + return form + + +# Finally, it remains to define the solver. Since :class:`MeshSeq` contains information +# on how to build the function spaces, initial and boundary conditions for the tracer +# concentration, we can simply define them as in the previous demos. +# On the other hand, so far we have not specified how to handle the background velocity +# field. We will do this here. :: + + +def get_solver(mesh_seq): + def solver(index, ic): + tp = mesh_seq.time_partition + t_start, t_end = tp.subintervals[index] + dt = tp.timesteps[index] + + # Initialise the concentration fields + Q = mesh_seq.function_spaces["c"][index] + c = Function(Q, name="c") + c_ = Function(Q, name="c_old").assign(ic["c"]) + + # Specify the velocity function space + V = VectorFunctionSpace(mesh_seq[index], "CG", 1) + u = Function(V) + u_ = Function(V) + + # Compute the velocity field at t_start and assign it to u_ + x, y = SpatialCoordinate(mesh_seq[index]) + u_.interpolate(velocity_expression(x, y, t_start)) + + # We pass both the concentration and velocity Functions to get_form + form_fields = {"c": (c, c_), "u": (u, u_)} + F = mesh_seq.form(index, form_fields)["c"] + nlvp = NonlinearVariationalProblem(F, c, bcs=mesh_seq.bcs(index)) + nlvs = NonlinearVariationalSolver(nlvp, ad_block_tag="c") + + # Time integrate from t_start to t_end + t = t_start + while t < t_end - 0.5 * dt: + # update the background velocity field at the current timestep + u.interpolate(velocity_expression(x, y, t)) + + # solve the advection equation + nlvs.solve() + + # update the 'lagged' concentration and velocity field + c_.assign(c) + u_.assign(u) + t += dt + + return {"c": c} + + return solver + + +# Now that we have defined the solver and specified how to build and update the +# background velocity, we are ready to solve the problem. We run the simulation +# until :math:`t=3` on a sequence of two coarse meshes. :: + +n = 64 +meshes = [UnitSquareMesh(n, n), UnitSquareMesh(n, n)] +end_time = period / 2 +dt = 0.01 +time_partition = TimePartition( + end_time, + len(meshes), + dt, + fields, + num_timesteps_per_export=50, +) + +msq = MeshSeq( + time_partition, + meshes, + get_function_spaces=get_function_spaces, + get_initial_condition=get_initial_condition, + get_bcs=get_bcs, + get_form=get_form, + get_solver=get_solver, +) +solutions = msq.solve_forward() + +# Let us plot the solutions :math:`c` at exported timesteps. :: + +plot_kwargs = {"cmap": "coolwarm", "levels": 100} +fig, _, _ = plot_snapshots(solutions, time_partition, "c", "forward", **plot_kwargs) +fig.savefig("bubble_shear.jpg") + +# .. figure:: bubble_shear.jpg +# :figwidth: 80% +# :align: center +# +# We observe that the initial tracer bubble has been deformed by the strong shear forces. +# However, we notice that the deformation appears to be reverting. This is due to the +# periodicity of the velocity field :math:`\mathbf{u}`, which we chose to be an odd +# function of time. This means that the velocity field reverses direction at :math:`t=T/2`, +# and so does the deformation of the tracer bubble, which returns to its initial shape. +# +# However, while the tracer bubble at :math:`t=T/2` does resemble a very blurry circular +# shape, it is far from matching the initial condition. This "bluriness" is the result of +# adding the SUPG stabilisation to the weak form, which adds numerical diffusion. This +# 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 Peclet 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, and to visit the `bubble shear mesh adaptation demo <./bubble_shear-goal.py.html>`__, +# where we use a goal-oriented approach to identify areas of the domain that require +# refinement to reduce the numerical diffusion. +# +# This tutorial can be dowloaded as a `Python script `__. From 5da1a0465037ef9d085719622880dcfa5851a802 Mon Sep 17 00:00:00 2001 From: ddundo Date: Tue, 12 Mar 2024 00:18:11 +0000 Subject: [PATCH 02/10] #28: Add reference --- demos/demo_references.bib | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/demos/demo_references.bib b/demos/demo_references.bib index 31fe2797..094f853d 100644 --- a/demos/demo_references.bib +++ b/demos/demo_references.bib @@ -16,3 +16,10 @@ @Book{Hundsdorfer:2003 year={2003}, publisher={Springer} } + +@article{barral2016anisotropic, + 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} +} From 92b5dd3a3c8a3dd384fb6bd502a9525385097b30 Mon Sep 17 00:00:00 2001 From: ddundo Date: Tue, 12 Mar 2024 02:57:49 +0000 Subject: [PATCH 03/10] #28: Bubble shear GO adaptation demo --- demos/bubble_shear-goal.py | 296 +++++++++++++++++++++++++++++++++++++ demos/bubble_shear.py | 2 +- 2 files changed, 297 insertions(+), 1 deletion(-) create mode 100644 demos/bubble_shear-goal.py diff --git a/demos/bubble_shear-goal.py b/demos/bubble_shear-goal.py new file mode 100644 index 00000000..c4172689 --- /dev/null +++ b/demos/bubble_shear-goal.py @@ -0,0 +1,296 @@ +# Goal-oriented mesh adaptation for a 2D time-dependent problem +# =================================================== + +# In a `previous demo <./bubble_shear.py.html>`__, we considered the advection of a scalar +# concentration field in a non-uniform background velocity field. The demo concluded with +# labelling the problem as a good candidate for a goal-oriented mesh adaptation approach. +# This is what we set out to do in this demo. +# +# We will use the same setup as in the previous demo, but a small modifiction to the +# :meth:`get_form` function is necessary in order to take into account a prescribed +# time-dependent background velocity field which is not passed to :class:`GoalOrientedMeshSeq`. :: + +from firedrake import * +from goalie_adjoint import * +from animate.metric import RiemannianMetric +from animate.adapt import adapt + +period = 6.0 + + +def velocity_expression(x, y, t): + u_expr = as_vector( + [ + 2 * sin(pi * x) ** 2 * sin(2 * pi * y) * cos(2 * pi * t / period), + -sin(2 * pi * x) * sin(pi * y) ** 2 * cos(2 * pi * t / period), + ] + ) + return u_expr + + +fields = ["c"] + + +def get_function_spaces(mesh): + return {"c": FunctionSpace(mesh, "CG", 1)} + + +def get_bcs(mesh_seq): + def bcs(index): + return [DirichletBC(mesh_seq.function_spaces["c"][index], 0.0, "on_boundary")] + + return bcs + + +def ball_initial_condition(x, y): + ball_r0, ball_x0, ball_y0 = 0.15, 0.5, 0.65 + r = sqrt(pow(x - ball_x0, 2) + pow(y - ball_y0, 2)) + return conditional(r < ball_r0, 1.0, 0.0) + + +def get_initial_condition(mesh_seq): + fs = mesh_seq.function_spaces["c"][0] + x, y = SpatialCoordinate(mesh_seq[0]) + ball = ball_initial_condition(x, y) + c0 = assemble(interpolate(ball, fs)) + return {"c": c0} + + +def get_solver(mesh_seq): + def solver(index, ic): + tp = mesh_seq.time_partition + t_start, t_end = tp.subintervals[index] + dt = tp.timesteps[index] + + # Initialise the concentration fields + Q = mesh_seq.function_spaces["c"][index] + c = Function(Q, name="c") + c_ = Function(Q, name="c_old").assign(ic["c"]) + + # Specify the velocity function space + V = VectorFunctionSpace(mesh_seq[index], "CG", 1) + u = Function(V) + u_ = Function(V) + + # Compute the velocity field at t_start and assign it to u_ + x, y = SpatialCoordinate(mesh_seq[index]) + u_.interpolate(velocity_expression(x, y, t_start)) + + # We pass both the concentration and velocity Functions to get_form + form_fields = {"c": (c, c_), "u": (u, u_)} + F = mesh_seq.form(index, form_fields)["c"] + nlvp = NonlinearVariationalProblem(F, c, bcs=mesh_seq.bcs(index)) + nlvs = NonlinearVariationalSolver(nlvp, ad_block_tag="c") + + # Time integrate from t_start to t_end + t = t_start + while t < t_end - 0.5 * dt: + # update the background velocity field at the current timestep + u.interpolate(velocity_expression(x, y, t)) + + # solve the advection equation + nlvs.solve() + + # update the 'lagged' concentration and velocity field + c_.assign(c) + u_.assign(u) + t += dt + + return {"c": c} + + return solver + + +# In the previous demo, the :meth:`get_form` method was only called from within the +# :meth:`get_solver`, so we were able to pass the velocity field as an argument when +# calling :meth:`get_form`. However, in the context of goal-oriented mesh adaptation, +# the :meth:`get_form` method is called from within :class:`GoalOrientedMeshSeq` +# while computing error indicators. There, only those fields that are passed to +# :class:`GoalOrientedMeshSeq` are passed to :meth:`get_form`. This means that the velocity +# field would not be updated in time. In such a case, we will manually compute the +# velocity field using the current simulation time, which will be passed automatically +# as an :arg:`err_ind_time` kwarg in :meth:`GoalOrientedMeshSeq.indicate_errors`. :: + + +def get_form(mesh_seq): + def form(index, form_fields, err_ind_time=None): + Q = mesh_seq.function_spaces["c"][index] + c, c_ = form_fields["c"] + + if "u" in form_fields: + u, u_ = form_fields["u"] + else: + x, y = SpatialCoordinate(mesh_seq[index]) + V = VectorFunctionSpace(mesh_seq[index], "CG", 1) + u = Function(V).interpolate(velocity_expression(x, y, err_ind_time)) + u_ = Function(V).interpolate(velocity_expression(x, y, err_ind_time)) + + # the rest remains unchanged + + R = FunctionSpace(mesh_seq[index], "R", 0) + dt = Function(R).assign(mesh_seq.time_partition.timesteps[index]) + theta = Function(R).assign(0.5) # Crank-Nicolson theta + + # SUPG stabilisation parameter + D = Function(R).assign(0.1) # diffusivity coefficient + h = CellSize(mesh_seq[index]) # mesh cell size + U = sqrt(dot(u, u)) # velocity magnitude + tau = 0.5 * h / U + tau = min_value(tau, U * h / (6 * D)) + + phi = TestFunction(Q) + phi += tau * dot(u, grad(phi)) + + a = c * phi * dx + dt * theta * dot(u, grad(c)) * phi * dx + L = c_ * phi * dx - dt * (1 - theta) * dot(u_, grad(c_)) * phi * dx + F = a - L + + return {"c": F} + + return form + + +# To compute the adjoint, we must also define the goal functional. As motivated in +# the previous demo, we will use the L2 norm of the difference between the +# concentration field at :math:`t=0` and :math:`t=T/2`, i.e. the simulation end time. +# +# .. math:: +# J(c) := \int_{\Omega} (c(x, y, T/2) - c_0(x, y))^2 \, dx \, dy. :: + + +def get_qoi(self, sols, index): + def qoi(): + c0 = self.get_initial_condition()["c"] + c0_proj = project(c0, self.function_spaces["c"][index]) + c = sols["c"] + + J = (c - c0_proj) * (c - c0_proj) * dx + return J + + return qoi + + +# We must also define the adaptor function to drive the mesh adaptation process. Here +# we compute the anisotropic DWR metric which requires us to compute the hessian of the +# tracer concentration field. We combine each metric in a subinterval by intersection. +# We will only run two iterations of the fixed point iteration algorithm so we do not +# ramp up the complexity. :: + + +def adaptor(mesh_seq, solutions, indicators): + iteration = mesh_seq.fp_iteration + mp = { + "dm_plex_metric": { + "target_complexity": 1000, + "p": 1.0, + "h_min": 1e-04, + "h_max": 1.0, + } + } + + metrics = [] + for i, mesh in enumerate(mesh_seq): + c_solutions = solutions["c"]["forward"][i] + + P1_ten = TensorFunctionSpace(mesh, "CG", 1) + subinterval_metric = RiemannianMetric(P1_ten) + + for j, c_sol in enumerate(c_solutions): + # Recover the Hessian of the forward solution + hessian = RiemannianMetric(P1_ten) + hessian.compute_hessian(c_sol) + + metric = RiemannianMetric(P1_ten) + metric.set_parameters(mp) + + # Compute an anisotropic metric from the error indicator field and Hessian + metric.compute_anisotropic_dwr_metric(indicators["c"][i][j], hessian) + + # Combine the metrics from each exported timestep in the subinterval + subinterval_metric.intersect(metric) + metrics.append(metric) + + # Normalise the metrics in space and time + space_time_normalise(metrics, mesh_seq.time_partition, mp) + + complexities = np.zeros(len(mesh_seq)) + for i, metric in enumerate(metrics): + complexities[i] = metric.complexity() + mesh_seq[i] = adapt(mesh_seq[i], metric) + + num_dofs = mesh_seq.count_vertices() + num_elem = mesh_seq.count_elements() + pyrint(f"Fixed point iteration {iteration}:") + for i, (complexity, ndofs, nelem) in enumerate( + zip(complexities, num_dofs, num_elem) + ): + pyrint( + f" subinterval {i}, complexity: {complexity:4.0f}" + f", dofs: {ndofs:4d}, elements: {nelem:4d}" + ) + + return True + + +# Finally, we can define the :class:`GoalOrientedMeshSeq` and run the fixed point iteration. +# For the purposes of this demo, we divide the time interval into 25 subintervals and only +# run two iterations of the fixed point iteration, which is not enough to reach convergence. :: + +num_subintervals = 25 +n = 50 +meshes = [UnitSquareMesh(n, n) for _ in range(num_subintervals)] +end_time = period / 2 +dt = 0.01 +time_partition = TimePartition( + end_time, + len(meshes), + dt, + fields, + num_timesteps_per_export=6, +) + +parameters = GoalOrientedMetricParameters({"maxiter": 2}) +msq = GoalOrientedMeshSeq( + time_partition, + meshes, + get_function_spaces=get_function_spaces, + get_initial_condition=get_initial_condition, + get_bcs=get_bcs, + get_form=get_form, + get_solver=get_solver, + get_qoi=get_qoi, + qoi_type="end_time", + parameters=parameters, +) +solutions, indicators = msq.fixed_point_iteration(adaptor) + +# Let us plot the intermediate and final solutions, as well as the final adapted meshes. + +import matplotlib.pyplot as plt +from firedrake.pyplot import tripcolor, triplot + +fig, axes = plt.subplots(3, 3, figsize=(8, 8), sharex=True, sharey=True) +for i, ax in enumerate(axes.flat): + ax.set_aspect("equal") + ax.set_xlim(0, 1) + ax.set_ylim(0, 1) + tripcolor(solutions["c"]["forward"][int(i * 3)][-1], axes=ax, cmap="coolwarm") + triplot(msq[int(i * 3)], axes=ax, interior_kw={"linewidth": 0.08}) + time = time_partition.subintervals[int(i * 3)][-1] + ax.annotate(f"t={time:.2f}", (0.05, 0.05), color="white") +fig.tight_layout(pad=0.4) +fig.savefig("bubble_shear-goal_final_meshes.jpg", dpi=300, bbox_inches="tight") + +# .. figure:: bubble_shear-goal_final_meshes.jpg +# :figwidth: 80% +# :align: center +# +# Compared to the solution obtained on a fixed uniform mesh in the previous demo, we +# observe that the final solution at :math:`t=T/2` better matches the initial tracer +# bubble :math:`c_0` and is less diffused. Despite employing only two fixed +# point iterations, the goal-oriented mesh adaptation process was still able to +# significantly improve the accuracy of the solution while reducing the number of +# degrees of freedom by half. We encourage further experimentation with the number of +# subintervals, adaptor functions and fixed point iterations necessary to reach convergence. +# +# This tutorial can be dowloaded as a `Python script `__. diff --git a/demos/bubble_shear.py b/demos/bubble_shear.py index b5f5197c..2b490df2 100644 --- a/demos/bubble_shear.py +++ b/demos/bubble_shear.py @@ -179,7 +179,7 @@ def solver(index, ic): # background velocity, we are ready to solve the problem. We run the simulation # until :math:`t=3` on a sequence of two coarse meshes. :: -n = 64 +n = 50 meshes = [UnitSquareMesh(n, n), UnitSquareMesh(n, n)] end_time = period / 2 dt = 0.01 From 06403d58534494aaf194c1d2489dca430a3fd2cd Mon Sep 17 00:00:00 2001 From: ddundo Date: Tue, 12 Mar 2024 03:02:35 +0000 Subject: [PATCH 04/10] #73: Accounting for time-dep constants in indicate_errors --- goalie/go_mesh_seq.py | 25 +++++++++++++++++++------ 1 file changed, 19 insertions(+), 6 deletions(-) diff --git a/goalie/go_mesh_seq.py b/goalie/go_mesh_seq.py index 1f67d86b..ab537ff1 100644 --- a/goalie/go_mesh_seq.py +++ b/goalie/go_mesh_seq.py @@ -13,6 +13,7 @@ import numpy as np from typing import Tuple import ufl +from inspect import signature __all__ = ["GoalOrientedMeshSeq"] @@ -130,6 +131,10 @@ def indicate_errors( self.solve_adjoint(**adj_kwargs) enriched_mesh_seq.solve_adjoint(**adj_kwargs) + tp = self.time_partition + # check whether get_form contains a kwarg indicating time-dependent constants + timedep_const = "err_ind_time" in signature(enriched_mesh_seq.form).parameters + FWD, ADJ = "forward", "adjoint" FWD_OLD = "forward" if self.steady else "forward_old" ADJ_NEXT = "adjoint" if self.steady else "adjoint_next" @@ -150,18 +155,26 @@ def indicate_errors( u_star_e[f] = Function(fs_e) # Get forms for each equation in enriched space - forms = enriched_mesh_seq.form(i, mapping) - if not isinstance(forms, dict): - raise TypeError( - "The function defined by get_form should return a dictionary" - f", not type '{type(forms)}'." - ) + if not timedep_const: + forms = enriched_mesh_seq.form(i, mapping) + if not isinstance(forms, dict): + raise TypeError( + "The function defined by get_form should return a dictionary" + f", not type '{type(forms)}'." + ) # Loop over each strongly coupled field for f in self.fields: # Loop over each timestep for j in range(len(self.solutions[f]["forward"][i])): # Update fields + if timedep_const: + # recompile the form with updated time-dependent constants + time = ( + tp.subintervals[i][0] + + (j + 1) * tp.timesteps[i] * tp.num_timesteps_per_export[i] + ) + forms = enriched_mesh_seq.form(i, mapping, err_ind_time=time) transfer(self.solutions[f][FWD][i][j], u[f]) transfer(self.solutions[f][FWD_OLD][i][j], u_[f]) transfer(self.solutions[f][ADJ][i][j], u_star[f]) From 57045a86c5cd45391474268d262ca0205f71b9be Mon Sep 17 00:00:00 2001 From: ddundo Date: Tue, 12 Mar 2024 13:47:52 +0000 Subject: [PATCH 05/10] #28: Reduce the bubble shear demos cost during testing --- demos/bubble_shear-goal.py | 25 ++++++++++++++++++------- demos/bubble_shear.py | 14 +++++++++----- 2 files changed, 27 insertions(+), 12 deletions(-) diff --git a/demos/bubble_shear-goal.py b/demos/bubble_shear-goal.py index c4172689..49a40253 100644 --- a/demos/bubble_shear-goal.py +++ b/demos/bubble_shear-goal.py @@ -15,6 +15,8 @@ from animate.metric import RiemannianMetric from animate.adapt import adapt +set_log_level(DEBUG) + period = 6.0 @@ -83,8 +85,9 @@ def solver(index, ic): nlvs = NonlinearVariationalSolver(nlvp, ad_block_tag="c") # Time integrate from t_start to t_end - t = t_start - while t < t_end - 0.5 * dt: + t = t_start + dt + while t < t_end + 0.5 * dt: + print(t) # update the background velocity field at the current timestep u.interpolate(velocity_expression(x, y, t)) @@ -236,11 +239,15 @@ def adaptor(mesh_seq, solutions, indicators): # For the purposes of this demo, we divide the time interval into 25 subintervals and only # run two iterations of the fixed point iteration, which is not enough to reach convergence. :: +# Reduce the cost of the demo during testing +test = os.environ.get("GOALIE_REGRESSION_TEST") is not None +n = 50 if not test else 5 +dt = 0.01 if not test else 0.02 +maxiter = 2 if not test else 1 # maximum number of fixed point iterations + num_subintervals = 25 -n = 50 meshes = [UnitSquareMesh(n, n) for _ in range(num_subintervals)] end_time = period / 2 -dt = 0.01 time_partition = TimePartition( end_time, len(meshes), @@ -249,7 +256,7 @@ def adaptor(mesh_seq, solutions, indicators): num_timesteps_per_export=6, ) -parameters = GoalOrientedMetricParameters({"maxiter": 2}) +parameters = GoalOrientedMetricParameters({"maxiter": maxiter}) msq = GoalOrientedMeshSeq( time_partition, meshes, @@ -290,7 +297,11 @@ def adaptor(mesh_seq, solutions, indicators): # bubble :math:`c_0` and is less diffused. Despite employing only two fixed # point iterations, the goal-oriented mesh adaptation process was still able to # significantly improve the accuracy of the solution while reducing the number of -# degrees of freedom by half. We encourage further experimentation with the number of -# subintervals, adaptor functions and fixed point iterations necessary to reach convergence. +# degrees of freedom by half. +# +# We encourage users to experiment with different numbers of subintervals, adaptor +# functions and metric parameters to explore the convergence of the fixed point iteration. +# Note that you will likely also need to decrease the timestep size to maintain stability +# as the mesh is refined. # # This tutorial can be dowloaded as a `Python script `__. diff --git a/demos/bubble_shear.py b/demos/bubble_shear.py index 2b490df2..2b2658a4 100644 --- a/demos/bubble_shear.py +++ b/demos/bubble_shear.py @@ -61,7 +61,8 @@ def get_function_spaces(mesh): # We proceed similarly with prescribing initial and boundary 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. :: +# and :math:`0` elsewhere in the domain. Note that this is a discontinuous function which +# will not be represented well on a coarse uniform mesh. :: def get_bcs(mesh_seq): @@ -157,8 +158,8 @@ def solver(index, ic): nlvs = NonlinearVariationalSolver(nlvp, ad_block_tag="c") # Time integrate from t_start to t_end - t = t_start - while t < t_end - 0.5 * dt: + t = t_start + dt + while t < t_end + 0.5 * dt: # update the background velocity field at the current timestep u.interpolate(velocity_expression(x, y, t)) @@ -177,9 +178,12 @@ def solver(index, ic): # Now that we have defined the solver and specified how to build and update the # background velocity, we are ready to solve the problem. We run the simulation -# until :math:`t=3` on a sequence of two coarse meshes. :: +# until :math:`t=T/2` on a sequence of two coarse meshes. :: + +# Reduce the cost of the demo during testing. +test = os.environ.get("GOALIE_REGRESSION_TEST") is not None +n = 50 if not test else 5 -n = 50 meshes = [UnitSquareMesh(n, n), UnitSquareMesh(n, n)] end_time = period / 2 dt = 0.01 From a7d2a6401fde6011ab0b62840d00292fbcdfef4c Mon Sep 17 00:00:00 2001 From: ddundo Date: Tue, 12 Mar 2024 16:29:45 +0000 Subject: [PATCH 06/10] #28: Further reduce cost of bubble shear demos during testing --- demos/bubble_shear-goal.py | 42 ++++++++++++++++++++++---------------- demos/bubble_shear.py | 15 +++++++------- 2 files changed, 32 insertions(+), 25 deletions(-) diff --git a/demos/bubble_shear-goal.py b/demos/bubble_shear-goal.py index 49a40253..b401ed54 100644 --- a/demos/bubble_shear-goal.py +++ b/demos/bubble_shear-goal.py @@ -15,7 +15,6 @@ from animate.metric import RiemannianMetric from animate.adapt import adapt -set_log_level(DEBUG) period = 6.0 @@ -87,7 +86,6 @@ def solver(index, ic): # Time integrate from t_start to t_end t = t_start + dt while t < t_end + 0.5 * dt: - print(t) # update the background velocity field at the current timestep u.interpolate(velocity_expression(x, y, t)) @@ -132,7 +130,7 @@ def form(index, form_fields, err_ind_time=None): R = FunctionSpace(mesh_seq[index], "R", 0) dt = Function(R).assign(mesh_seq.time_partition.timesteps[index]) - theta = Function(R).assign(0.5) # Crank-Nicolson theta + theta = Function(R).assign(0.5) # Crank-Nicolson implicitness # SUPG stabilisation parameter D = Function(R).assign(0.1) # diffusivity coefficient @@ -236,16 +234,16 @@ def adaptor(mesh_seq, solutions, indicators): # Finally, we can define the :class:`GoalOrientedMeshSeq` and run the fixed point iteration. -# For the purposes of this demo, we divide the time interval into 25 subintervals and only +# For the purposes of this demo, we divide the time interval into 6 subintervals and only # run two iterations of the fixed point iteration, which is not enough to reach convergence. :: # Reduce the cost of the demo during testing test = os.environ.get("GOALIE_REGRESSION_TEST") is not None n = 50 if not test else 5 -dt = 0.01 if not test else 0.02 +dt = 0.01 if not test else 0.05 maxiter = 2 if not test else 1 # maximum number of fixed point iterations -num_subintervals = 25 +num_subintervals = 6 meshes = [UnitSquareMesh(n, n) for _ in range(num_subintervals)] end_time = period / 2 time_partition = TimePartition( @@ -253,7 +251,7 @@ def adaptor(mesh_seq, solutions, indicators): len(meshes), dt, fields, - num_timesteps_per_export=6, + num_timesteps_per_export=5, ) parameters = GoalOrientedMetricParameters({"maxiter": maxiter}) @@ -271,21 +269,26 @@ def adaptor(mesh_seq, solutions, indicators): ) solutions, indicators = msq.fixed_point_iteration(adaptor) -# Let us plot the intermediate and final solutions, as well as the final adapted meshes. +# Let us plot the intermediate and final solutions, as well as the final adapted meshes. :: import matplotlib.pyplot as plt from firedrake.pyplot import tripcolor, triplot -fig, axes = plt.subplots(3, 3, figsize=(8, 8), sharex=True, sharey=True) -for i, ax in enumerate(axes.flat): +fig, axes = plt.subplots(2, 3, figsize=(7.5, 5), sharex=True, sharey=True) + +for i, ax in enumerate(axes.flatten()): ax.set_aspect("equal") ax.set_xlim(0, 1) ax.set_ylim(0, 1) - tripcolor(solutions["c"]["forward"][int(i * 3)][-1], axes=ax, cmap="coolwarm") - triplot(msq[int(i * 3)], axes=ax, interior_kw={"linewidth": 0.08}) - time = time_partition.subintervals[int(i * 3)][-1] + # ax.tick_params(axis='both', which='major', labelsize=6) + + # plot the solution at the final subinterval timestep + time = time_partition.subintervals[i][-1] + tripcolor(solutions["c"]["forward"][i][-1], axes=ax, cmap="coolwarm") + triplot(msq[i], axes=ax, interior_kw={"linewidth": 0.08}) ax.annotate(f"t={time:.2f}", (0.05, 0.05), color="white") -fig.tight_layout(pad=0.4) + +fig.tight_layout(pad=0.3) fig.savefig("bubble_shear-goal_final_meshes.jpg", dpi=300, bbox_inches="tight") # .. figure:: bubble_shear-goal_final_meshes.jpg @@ -299,9 +302,12 @@ def adaptor(mesh_seq, solutions, indicators): # significantly improve the accuracy of the solution while reducing the number of # degrees of freedom by half. # -# We encourage users to experiment with different numbers of subintervals, adaptor -# functions and metric parameters to explore the convergence of the fixed point iteration. -# Note that you will likely also need to decrease the timestep size to maintain stability -# as the mesh is refined. +# As shown in the above figure, the bubble experiences extreme deformation and requires +# frequent adaptation to be resolved accurately (surely more than 5, as in this demo!). +# We encourage users to experiment with different numbers of subintervals, number of +# exports per subinterval, adaptor functions, and other parameters to explore the +# convergence of the fixed point iteration. Another interesting experiment would be to +# compare the impact of switching to an explicit time integration and using a smaller +# timestep to maintain numerical stability (look at CFL condition). # # This tutorial can be dowloaded as a `Python script `__. diff --git a/demos/bubble_shear.py b/demos/bubble_shear.py index 2b2658a4..40ed77ee 100644 --- a/demos/bubble_shear.py +++ b/demos/bubble_shear.py @@ -92,7 +92,7 @@ def get_initial_condition(mesh_seq): # define the velocity field. # As in the point discharge with diffusion demo, we include additional `streamline # upwind Petrov Galerkin (SUPG)` stabilisation by modifying the test function :math:`\psi`. -# To advance in time, we use Crank-Nicolson discretisation. :: +# To advance in time, we use Crank-Nicolson timestepping. :: def get_form(mesh_seq): @@ -103,7 +103,7 @@ def form(index, form_fields): R = FunctionSpace(mesh_seq[index], "R", 0) dt = Function(R).assign(mesh_seq.time_partition.timesteps[index]) - theta = Function(R).assign(0.5) # Crank-Nicolson theta + theta = Function(R).assign(0.5) # Crank-Nicolson implicitness # SUPG stabilisation parameter D = Function(R).assign(0.1) # diffusivity coefficient @@ -182,17 +182,18 @@ def solver(index, ic): # Reduce the cost of the demo during testing. test = os.environ.get("GOALIE_REGRESSION_TEST") is not None -n = 50 if not test else 5 +n = 50 if not test else 10 +dt = 0.01 if not test else 0.025 -meshes = [UnitSquareMesh(n, n), UnitSquareMesh(n, n)] +num_subintervals = 2 +meshes = [UnitSquareMesh(n, n) for _ in range(num_subintervals)] end_time = period / 2 -dt = 0.01 time_partition = TimePartition( end_time, len(meshes), dt, fields, - num_timesteps_per_export=50, + num_timesteps_per_export=30, ) msq = MeshSeq( @@ -206,7 +207,7 @@ def solver(index, ic): ) solutions = msq.solve_forward() -# Let us plot the solutions :math:`c` at exported timesteps. :: +# Let us plot the solution :math:`c` at exported timesteps. :: plot_kwargs = {"cmap": "coolwarm", "levels": 100} fig, _, _ = plot_snapshots(solutions, time_partition, "c", "forward", **plot_kwargs) From 262c78a51a99714a4794b55259edf6fc5aea3567 Mon Sep 17 00:00:00 2001 From: ddundo Date: Fri, 15 Mar 2024 10:10:11 +0000 Subject: [PATCH 07/10] #28: Minor changes --- ...-goal.py => bubble_shear-goal_oriented.py} | 53 ++++++++------- demos/bubble_shear.py | 65 ++++++++++--------- demos/demo_references.bib | 2 +- 3 files changed, 63 insertions(+), 57 deletions(-) rename demos/{bubble_shear-goal.py => bubble_shear-goal_oriented.py} (84%) diff --git a/demos/bubble_shear-goal.py b/demos/bubble_shear-goal_oriented.py similarity index 84% rename from demos/bubble_shear-goal.py rename to demos/bubble_shear-goal_oriented.py index b401ed54..a559a131 100644 --- a/demos/bubble_shear-goal.py +++ b/demos/bubble_shear-goal_oriented.py @@ -1,14 +1,15 @@ # Goal-oriented mesh adaptation for a 2D time-dependent problem -# =================================================== +# ============================================================= -# In a `previous demo <./bubble_shear.py.html>`__, we considered the advection of a scalar -# concentration field in a non-uniform background velocity field. The demo concluded with -# labelling the problem as a good candidate for a goal-oriented mesh adaptation approach. -# This is what we set out to do in this demo. +# In a `previous demo <./bubble_shear.py.html>`__, we considered the advection of a +# scalar concentration field in a non-uniform background velocity field. The demo +# concluded with labelling the problem as a good candidate for a goal-oriented mesh +# adaptation approach. This is what we set out to do in this demo. # -# We will use the same setup as in the previous demo, but a small modifiction to the +# We will use the same setup as in the original demo, but a small modifiction to the # :meth:`get_form` function is necessary in order to take into account a prescribed -# time-dependent background velocity field which is not passed to :class:`GoalOrientedMeshSeq`. :: +# time-dependent background velocity field which is not passed to +# :class:`GoalOrientedMeshSeq`. :: from firedrake import * from goalie_adjoint import * @@ -102,15 +103,16 @@ def solver(index, ic): return solver -# In the previous demo, the :meth:`get_form` method was only called from within the -# :meth:`get_solver`, so we were able to pass the velocity field as an argument when -# calling :meth:`get_form`. However, in the context of goal-oriented mesh adaptation, -# the :meth:`get_form` method is called from within :class:`GoalOrientedMeshSeq` -# while computing error indicators. There, only those fields that are passed to -# :class:`GoalOrientedMeshSeq` are passed to :meth:`get_form`. This means that the velocity -# field would not be updated in time. In such a case, we will manually compute the -# velocity field using the current simulation time, which will be passed automatically -# as an :arg:`err_ind_time` kwarg in :meth:`GoalOrientedMeshSeq.indicate_errors`. :: +# In the `first demo `__ where we considered this problem, the +# :meth:`get_form` method was only called from within the :meth:`get_solver`, so we were +# able to pass the velocity field as an argument when calling :meth:`get_form`. However, +# in the context of goal-oriented mesh adaptation, the :meth:`get_form` method is +# called from within :class:`GoalOrientedMeshSeq` while computing error indicators. +# There, only those fields that are passed to :class:`GoalOrientedMeshSeq` are passed to +# :meth:`get_form`. This means that the velocity field would not be updated in time. In +# such a case, we will manually compute the velocity field using the current simulation +# time, which will be passed automatically as an :arg:`err_ind_time` kwarg in +# :meth:`GoalOrientedMeshSeq.indicate_errors`. :: def get_form(mesh_seq): @@ -126,7 +128,7 @@ def form(index, form_fields, err_ind_time=None): u = Function(V).interpolate(velocity_expression(x, y, err_ind_time)) u_ = Function(V).interpolate(velocity_expression(x, y, err_ind_time)) - # the rest remains unchanged + # The rest remains unchanged R = FunctionSpace(mesh_seq[index], "R", 0) dt = Function(R).assign(mesh_seq.time_partition.timesteps[index]) @@ -233,9 +235,10 @@ def adaptor(mesh_seq, solutions, indicators): return True -# Finally, we can define the :class:`GoalOrientedMeshSeq` and run the fixed point iteration. -# For the purposes of this demo, we divide the time interval into 6 subintervals and only -# run two iterations of the fixed point iteration, which is not enough to reach convergence. :: +# Finally, we can define the :class:`GoalOrientedMeshSeq` and run the fixed point +# iteration. For the purposes of this demo, we divide the time interval into 6 +# subintervals and only run two iterations of the fixed point iteration, which is not +# enough to reach convergence. :: # Reduce the cost of the demo during testing test = os.environ.get("GOALIE_REGRESSION_TEST") is not None @@ -248,7 +251,7 @@ def adaptor(mesh_seq, solutions, indicators): end_time = period / 2 time_partition = TimePartition( end_time, - len(meshes), + num_subintervals, dt, fields, num_timesteps_per_export=5, @@ -282,16 +285,16 @@ def adaptor(mesh_seq, solutions, indicators): ax.set_ylim(0, 1) # ax.tick_params(axis='both', which='major', labelsize=6) - # plot the solution at the final subinterval timestep + # Plot the solution at the final subinterval timestep time = time_partition.subintervals[i][-1] tripcolor(solutions["c"]["forward"][i][-1], axes=ax, cmap="coolwarm") triplot(msq[i], axes=ax, interior_kw={"linewidth": 0.08}) ax.annotate(f"t={time:.2f}", (0.05, 0.05), color="white") fig.tight_layout(pad=0.3) -fig.savefig("bubble_shear-goal_final_meshes.jpg", dpi=300, bbox_inches="tight") +fig.savefig("bubble_shear-goal_oriented.jpg", dpi=300, bbox_inches="tight") -# .. figure:: bubble_shear-goal_final_meshes.jpg +# .. figure:: bubble_shear-goal_oriented.jpg # :figwidth: 80% # :align: center # @@ -310,4 +313,4 @@ def adaptor(mesh_seq, solutions, indicators): # compare the impact of switching to an explicit time integration and using a smaller # timestep to maintain numerical stability (look at CFL condition). # -# This tutorial can be dowloaded as a `Python script `__. +# This tutorial can be dowloaded as a `Python script `__. diff --git a/demos/bubble_shear.py b/demos/bubble_shear.py index 40ed77ee..b95e12e7 100644 --- a/demos/bubble_shear.py +++ b/demos/bubble_shear.py @@ -1,10 +1,11 @@ # Advection of a tracer bubble in a non-uniform shear flow -# =================================================== +# ======================================================== # In this demo, we consider the advection of a scalar concentration field, similarly to -# our previous solid body rotation demo. However, we now consider a background velocity -# field that is not uniform in time. We further prescribe homogeneous Dirichlet boundary -# conditions for the concentration, so in this demo we solve +# the `solid body rotation demo <./solid_body_rotation.py.html>`__. However, we now +# consider a background velocity field that is not uniform in time. We further prescribe +# homogeneous Dirichlet boundary conditions for the concentration, so in this demo we +# solve the following advection problem: # # .. math:: # \begin{array}{rl} @@ -13,9 +14,9 @@ # 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. +# 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. # # First, we import Firedrake and Goalie. :: @@ -24,7 +25,7 @@ # We begin by definining the background velocity field :math:`\mathbf{u}(x, y, t)`. # Specifically, we choose the velocity field that is periodic in time, such as -# in TODO: cite, and is given by +# in :cite:, and is given by # # .. math:: # \mathbf{u}(x, y, t) := \cos(2\pi t/T)\left(2\sin^2(\pi x)\sin(2\pi y), -\sin(2\pi x)\sin^2(\pi y) \right), @@ -59,10 +60,10 @@ def get_function_spaces(mesh): # We proceed similarly with prescribing initial and boundary 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 this is a discontinuous function which -# will not be represented well on a coarse uniform mesh. :: +# 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 this is a discontinuous function +# which will not be represented well on a coarse uniform mesh. :: def get_bcs(mesh_seq): @@ -91,8 +92,8 @@ def get_initial_condition(mesh_seq): # will always pass both these fields from :meth:`get_solver`. Note that we still do not # define the velocity field. # As in the point discharge with diffusion demo, we include additional `streamline -# upwind Petrov Galerkin (SUPG)` stabilisation by modifying the test function :math:`\psi`. -# To advance in time, we use Crank-Nicolson timestepping. :: +# upwind Petrov Galerkin (SUPG)` stabilisation by modifying the test function +# :math:`\psi`. To advance in time, we use Crank-Nicolson timestepping. :: def get_form(mesh_seq): @@ -177,8 +178,8 @@ def solver(index, ic): # Now that we have defined the solver and specified how to build and update the -# background velocity, we are ready to solve the problem. We run the simulation -# until :math:`t=T/2` on a sequence of two coarse meshes. :: +# background velocity, we are ready to solve the problem. We run the simulation until +# :math:`t=T/2` on a sequence of two coarse meshes. :: # Reduce the cost of the demo during testing. test = os.environ.get("GOALIE_REGRESSION_TEST") is not None @@ -190,7 +191,7 @@ def solver(index, ic): end_time = period / 2 time_partition = TimePartition( end_time, - len(meshes), + num_subintervals, dt, fields, num_timesteps_per_export=30, @@ -217,21 +218,23 @@ def solver(index, ic): # :figwidth: 80% # :align: center # -# We observe that the initial tracer bubble has been deformed by the strong shear forces. -# However, we notice that the deformation appears to be reverting. This is due to the -# periodicity of the velocity field :math:`\mathbf{u}`, which we chose to be an odd -# function of time. This means that the velocity field reverses direction at :math:`t=T/2`, -# and so does the deformation of the tracer bubble, which returns to its initial shape. +# In the first part of the simulation, we observe that the initial tracer bubble becomes +# increasingly deformed by the strong shear forces. However, we notice that the +# deformation appears to be reverting. This is due to the periodicity of the velocity +# field :math:`\mathbf{u}`, which we chose to be an odd function of time. This means +# that the velocity field reverses direction at :math:`t=T/2`, and so does the +# deformation of the tracer bubble, which returns to its initial shape. # # However, while the tracer bubble at :math:`t=T/2` does resemble a very blurry circular -# shape, it is far from matching the initial condition. This "bluriness" is the result of -# adding the SUPG stabilisation to the weak form, which adds numerical diffusion. This -# 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 Peclet 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, and to visit the `bubble shear mesh adaptation demo <./bubble_shear-goal.py.html>`__, -# where we use a goal-oriented approach to identify areas of the domain that require -# refinement to reduce the numerical diffusion. +# shape, it is far from matching the initial condition. This "bluriness" is the result +# of adding the SUPG stabilisation to the weak form, which adds numerical diffusion. +# Adding 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, and to visit the +# `bubble shear mesh adaptation demo <./bubble_shear-goal.py.html>`__, where we use a +# goal-oriented approach to identify areas of the domain that require refinement to +# reduce the numerical diffusion. # # This tutorial can be dowloaded as a `Python script `__. diff --git a/demos/demo_references.bib b/demos/demo_references.bib index 094f853d..d8f10407 100644 --- a/demos/demo_references.bib +++ b/demos/demo_references.bib @@ -17,7 +17,7 @@ @Book{Hundsdorfer:2003 publisher={Springer} } -@article{barral2016anisotropic, +@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}, From e1558f544dcadf4e19eb26c9c2f3d311f3a59458 Mon Sep 17 00:00:00 2001 From: Joe Wallwork <22053413+jwallwork23@users.noreply.github.com> Date: Sun, 21 Apr 2024 06:31:31 +0100 Subject: [PATCH 08/10] Introduce R-space time `Function` (#172) Introduce R-space time `Function`; pull in recent changes from `main`. --------- Co-authored-by: Davor Dundovic <33790330+ddundo@users.noreply.github.com> Co-authored-by: acse-ej321 <89605848+acse-ej321@users.noreply.github.com> Co-authored-by: acse-ej321 Co-authored-by: ddundo --- .flake8 | 18 - .github/workflows/test_suite.yml | 28 +- .gitignore | 4 +- .pre-commit-config.yaml | 13 +- Makefile | 5 +- README.md | 54 +- demos/bubble_shear-goal_oriented.py | 50 +- demos/bubble_shear.py | 17 +- demos/burgers-hessian.py | 12 +- demos/burgers.py | 1 + demos/burgers1.py | 1 + demos/burgers2.py | 2 +- demos/burgers_ee.py | 1 + demos/burgers_oo.py | 2 +- demos/burgers_time_integrated.py | 1 + demos/gray_scott.py | 1 + demos/gray_scott_split.py | 1 + demos/mesh_seq.py | 5 +- demos/ode.py | 5 +- demos/point_discharge2d-goal_oriented.py | 8 +- demos/point_discharge2d-hessian.py | 4 +- demos/point_discharge2d.py | 1 + demos/solid_body_rotation.py | 4 +- demos/solid_body_rotation_split.py | 1 - goalie/__init__.py | 5 +- goalie/adjoint.py | 213 +++++--- goalie/error_estimation.py | 30 +- goalie/function_data.py | 185 +++++-- goalie/go_mesh_seq.py | 185 ++++--- goalie/interpolation.py | 150 ------ goalie/log.py | 17 +- goalie/math.py | 1 - goalie/mesh_seq.py | 507 +++++++++++++----- goalie/metric.py | 174 +++--- goalie/options.py | 54 +- goalie/point_seq.py | 3 +- goalie/time_partition.py | 87 +-- goalie/utility.py | 280 +--------- goalie_adjoint/__init__.py | 4 +- pyproject.toml | 65 ++- requirements.txt | 5 - setup.py | 12 - test/conftest.py | 4 +- test/sensors.py | 2 +- test/test_error_estimation.py | 6 +- test/test_function_data.py | 217 ++++++++ test/test_interpolation.py | 170 ------ test/test_math.py | 8 +- test/test_mesh_seq.py | 17 +- ...est_metric_normalise.py => test_metric.py} | 106 +++- test/test_metric_utils.py | 46 -- test/test_options.py | 8 +- test/test_time_partition.py | 4 +- test/test_utility.py | 230 +------- test/utility.py | 4 +- test_adjoint/conftest.py | 7 +- test_adjoint/examples/burgers.py | 1 - test_adjoint/examples/point_discharge2d.py | 5 +- test_adjoint/examples/point_discharge3d.py | 5 +- test_adjoint/examples/steady_flow_past_cyl.py | 4 +- test_adjoint/test_adjoint.py | 10 +- test_adjoint/test_demos.py | 5 +- test_adjoint/test_fp_iteration.py | 8 +- test_adjoint/test_mesh_seq.py | 18 +- test_adjoint/test_utils.py | 8 +- 65 files changed, 1526 insertions(+), 1583 deletions(-) delete mode 100644 .flake8 delete mode 100644 goalie/interpolation.py delete mode 100644 requirements.txt delete mode 100644 setup.py create mode 100644 test/test_function_data.py delete mode 100644 test/test_interpolation.py rename test/{test_metric_normalise.py => test_metric.py} (63%) delete mode 100644 test/test_metric_utils.py diff --git a/.flake8 b/.flake8 deleted file mode 100644 index 2ea09c6e..00000000 --- a/.flake8 +++ /dev/null @@ -1,18 +0,0 @@ -[flake8] -ignore = - # Line too long - E501, - # Missing whitespace around arithmetic operator - E226, - # Module level import not at top of file - E402, - # Do not use variables named "I", "O", or "l" - E741, - # Unable to detect undefined names - F403, - # Name may be undefined, or defined from star imports - F405, - # Line break occurred before binary operator - W503 -max-line-length = 88 -select = B,C,E,F,W,T4,B9 diff --git a/.github/workflows/test_suite.yml b/.github/workflows/test_suite.yml index c8c17297..8e37390b 100644 --- a/.github/workflows/test_suite.yml +++ b/.github/workflows/test_suite.yml @@ -1,4 +1,4 @@ -name: Build Goalie +name: 'Run Goalie test suite' on: # Run test suite whenever main is updated @@ -15,7 +15,7 @@ on: jobs: build: - name: "Build Goalie" + name: 'Test suite' # The type of runner that the job will run on runs-on: ubuntu-latest # The docker container to use. @@ -24,30 +24,36 @@ jobs: options: --user root steps: - uses: actions/checkout@v2 - - name: Cleanup + + - name: 'Cleanup' if: ${{ always() }} run: | cd .. rm -rf build - - name: Setup python + + - name: 'Setup Python' uses: actions/setup-python@v2 with: python-version: 3.8 - - name: Install Goalie + + - name: 'Install Goalie' run: | . /home/firedrake/firedrake/bin/activate python -m pip uninstall -y goalie python -m pip install -e . - - name: Test Goalie + + - name: 'Lint' + if: ${{ always() }} + run: | + . /home/firedrake/firedrake/bin/activate + make lint + + - name: 'Test Goalie' run: | . /home/firedrake/firedrake/bin/activate python $(which firedrake-clean) + export GITHUB_ACTIONS_TEST_RUN=1 python -m coverage erase python -m coverage run -a --source=goalie -m pytest -v --durations=20 test python -m coverage run -a --source=goalie -m pytest -v --durations=10 test_adjoint python -m coverage report - - name: Lint - if: ${{ always() }} - run: | - . /home/firedrake/firedrake/bin/activate - make lint diff --git a/.gitignore b/.gitignore index ffe94357..afe9f8ee 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,6 @@ -# Directories +# .gitignore + +# Subdirectories htmlcov* outputs/ plots/ diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 75ded481..b1ce0bc6 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,10 +1,7 @@ repos: -- repo: https://github.com/psf/black - rev: stable +- repo: https://github.com/astral-sh/ruff-pre-commit + rev: v0.3.7 hooks: - - id: black - language_version: python3.10 -- repo: https://github.com/pycqa/flake8 - rev: 3.7.9 - hooks: - - id: flake8 + - id: ruff # run the linter + args: [ --fix ] + - id: ruff-format # run the formatter \ No newline at end of file diff --git a/Makefile b/Makefile index 6da7a7a8..2cf64768 100644 --- a/Makefile +++ b/Makefile @@ -3,9 +3,6 @@ all: install .PHONY: demos test install: - @echo "Installing dependencies..." - @python3 -m pip install -r requirements.txt - @echo "Done." @echo "Installing Goalie..." @python3 -m pip install -e . @echo "Done." @@ -15,7 +12,7 @@ install: lint: @echo "Checking lint..." - @flake8 + @ruff check @echo "PASS" test: lint diff --git a/README.md b/README.md index 1089018c..d4a5e074 100644 --- a/README.md +++ b/README.md @@ -1,8 +1,7 @@ ## Goalie Goal-Oriented Mesh Adaptation Toolkit -![GitHub top language](https://img.shields.io/github/languages/top/pyroteus/goalie) -![GitHub repo size](https://img.shields.io/github/repo-size/pyroteus/goalie) -[![Code style: black](https://img.shields.io/badge/code%20style-black-000000.svg)](https://github.com/psf/black) +![GitHub top language](https://img.shields.io/github/languages/top/mesh-adaptation/goalie) +![GitHub repo size](https://img.shields.io/github/repo-size/mesh-adaptation/goalie) Goalie provides goal-oriented mesh adaptation pipelines for solving partial differential equations (PDEs) using adapted meshes built on the Python-based finite element library [Firedrake](http://www.firedrakeproject.org/). It runs a fixed point iteration loop for progressively solving time-dependent PDEs and their adjoints on sequences of meshes, performing goal-oriented error estimation, and adapting the meshes in sequence with a user-provided adaptor function until defined convergence criteria have been met. It is recommended that users are familiar with adjoint methods, mesh adaptation and the goal-oriented framework before starting with Goalie. @@ -10,51 +9,8 @@ For more information on Firedrake, please see: [Firedrake documentation](https: For more information on the implementation of the adjoint method, please see: [dolfin-adjoint documentation](http://www.dolfin-adjoint.org/en/latest/documentation/maths/index.html).  -For more information on the goal-oriented mesh adaptation, please see: [Goalie documentation](https://pyroteus.github.io/goalie/index.html) +For more information on the goal-oriented mesh adaptation, please see: [Goalie documentation](https://mesh-adaptation.github.io/goalie/index.html) -## Linux Installation +## Installation -The following installation instructions assume a Linux or WSL operating system. The options below include installation from a custom shell script which also installs the custom setup for Firedrake and PETSc. - -To use Goalie you will need a bespoke Firedrake installation. - -If Firedrake is already installed, please see instructions "To install Goalie via `git clone`". - -Additionally, although the Animate anisotropic mesh adaptation package is not a technical dependency for Goalie, as any bespoke adaptation method can be applied, the associated demos and tests do rely on Animate. For detail on installing Animate see: [Animate](https://github.com/pyroteus/animate) - -### To install Firedrake via shell script - -Firedrake, along with PETSc, is required by the Goalie package and is available for installation via a shell script. - -Instructions: -- Download installation files either: - - manually from: - - `install/install_firedrake_custom_mpi.sh` - - `install/petsc_options.txt` - - via curl: - - `curl -O https://raw.githubusercontent.com/pyroteus/animate/main/install/install_firedrake_custom_mpi.sh` - - `curl -O https://raw.githubusercontent.com/pyroteus/animate/main/install/petsc_options.txt` -- Install firedrake and associated dependencies to a local environment via `source install_firedrake_custom_mpi.sh` -- Continue to follow the instructions below in "To install Goalie via `git clone`" to complete the installation of Goalie. - -### To install Firedrake via Docker image - -A Firedrake docker image exists and can alternatively be downloaded and installed before installing Goalie. - -To install the Firedrake docker image: -- Pull the docker image: `docker pull jwallwork/firedrake-parmmg` -- Run the docker image: `docker run --rm -it -v ${HOME}:${HOME} jwallwork/firedrake-parmmg` - -Please note, that by installing via a Docker image with `${HOME}` you are giving Docker access to your home space. - -Continue to follow the instructions below in "To install Goalie via `git clone`" to complete the installation of Goalie. - -### To install Goalie via `git clone` - -Installing Goalie via cloning the GitHub repository assumes prior installation of Firedrake and its dependencies. For separate instructions for installing Firedrake please see: [Firedrake download instructions](https://www.firedrakeproject.org/download.html). - -To install Goalie via `git clone`: -- Activate your local virtual environment containing the bespoke Firedrake installation and navigate to the `src` folder. -- from the main Goalie GitHub, `git clone` the repository using HTTPS or SSH into the `src` folder -- `cd goalie` and run `make install` to install Goalie -- Execute the test suite to confirm installation was successful via `make test` +For installation instructions, we refer to the [Wiki page](https://github.com/mesh-adaptation/mesh-adaptation-docs/wiki/Installation-Instructions). diff --git a/demos/bubble_shear-goal_oriented.py b/demos/bubble_shear-goal_oriented.py index a559a131..adeddde0 100644 --- a/demos/bubble_shear-goal_oriented.py +++ b/demos/bubble_shear-goal_oriented.py @@ -11,11 +11,11 @@ # time-dependent background velocity field which is not passed to # :class:`GoalOrientedMeshSeq`. :: -from firedrake import * -from goalie_adjoint import * -from animate.metric import RiemannianMetric from animate.adapt import adapt +from animate.metric import RiemannianMetric +from firedrake import * +from goalie_adjoint import * period = 6.0 @@ -63,6 +63,7 @@ def solver(index, ic): tp = mesh_seq.time_partition t_start, t_end = tp.subintervals[index] dt = tp.timesteps[index] + time = mesh_seq.get_time(index) # Initialise the concentration fields Q = mesh_seq.function_spaces["c"][index] @@ -76,7 +77,7 @@ def solver(index, ic): # Compute the velocity field at t_start and assign it to u_ x, y = SpatialCoordinate(mesh_seq[index]) - u_.interpolate(velocity_expression(x, y, t_start)) + u_.interpolate(velocity_expression(x, y, time)) # We pass both the concentration and velocity Functions to get_form form_fields = {"c": (c, c_), "u": (u, u_)} @@ -85,18 +86,17 @@ def solver(index, ic): nlvs = NonlinearVariationalSolver(nlvp, ad_block_tag="c") # Time integrate from t_start to t_end - t = t_start + dt - while t < t_end + 0.5 * dt: - # update the background velocity field at the current timestep - u.interpolate(velocity_expression(x, y, t)) + while float(time) < t_end + 0.5 * dt: + # Update the background velocity field at the current timestep + u.interpolate(velocity_expression(x, y, time)) - # solve the advection equation + # Solve the advection equation nlvs.solve() - # update the 'lagged' concentration and velocity field + # Update the 'lagged' concentration and velocity field c_.assign(c) u_.assign(u) - t += dt + time += dt return {"c": c} @@ -109,24 +109,25 @@ def solver(index, ic): # in the context of goal-oriented mesh adaptation, the :meth:`get_form` method is # called from within :class:`GoalOrientedMeshSeq` while computing error indicators. # There, only those fields that are passed to :class:`GoalOrientedMeshSeq` are passed to -# :meth:`get_form`. This means that the velocity field would not be updated in time. In -# such a case, we will manually compute the velocity field using the current simulation -# time, which will be passed automatically as an :arg:`err_ind_time` kwarg in -# :meth:`GoalOrientedMeshSeq.indicate_errors`. :: +# :meth:`get_form`. Currently, Goalie does not consider passing non-prognostic fields +# to this method during the error-estimation step, so the velocity would not be updated +# in time. To account for this, we need to add some code for updating the velocity +# field when it is not present in the dictionary of fields passed. :: def get_form(mesh_seq): - def form(index, form_fields, err_ind_time=None): + def form(index, form_fields): Q = mesh_seq.function_spaces["c"][index] c, c_ = form_fields["c"] + time = mesh_seq.get_time(index) if "u" in form_fields: u, u_ = form_fields["u"] else: x, y = SpatialCoordinate(mesh_seq[index]) V = VectorFunctionSpace(mesh_seq[index], "CG", 1) - u = Function(V).interpolate(velocity_expression(x, y, err_ind_time)) - u_ = Function(V).interpolate(velocity_expression(x, y, err_ind_time)) + u = Function(V).interpolate(velocity_expression(x, y, time)) + u_ = Function(V).interpolate(velocity_expression(x, y, time)) # The rest remains unchanged @@ -216,19 +217,14 @@ def adaptor(mesh_seq, solutions, indicators): # Normalise the metrics in space and time space_time_normalise(metrics, mesh_seq.time_partition, mp) - complexities = np.zeros(len(mesh_seq)) - for i, metric in enumerate(metrics): - complexities[i] = metric.complexity() - mesh_seq[i] = adapt(mesh_seq[i], metric) - + # Apply mesh adaptation + mesh_seq.set_meshes(map(adapt, mesh_seq, metrics)) num_dofs = mesh_seq.count_vertices() num_elem = mesh_seq.count_elements() pyrint(f"Fixed point iteration {iteration}:") - for i, (complexity, ndofs, nelem) in enumerate( - zip(complexities, num_dofs, num_elem) - ): + for i, (ndofs, nelem, metric) in enumerate(zip(num_dofs, num_elem, metrics)): pyrint( - f" subinterval {i}, complexity: {complexity:4.0f}" + f" subinterval {i}, complexity: {metric.complexity():4.0f}" f", dofs: {ndofs:4d}, elements: {nelem:4d}" ) diff --git a/demos/bubble_shear.py b/demos/bubble_shear.py index b95e12e7..0d0c0bac 100644 --- a/demos/bubble_shear.py +++ b/demos/bubble_shear.py @@ -21,6 +21,7 @@ # First, we import Firedrake and Goalie. :: from firedrake import * + from goalie import * # We begin by definining the background velocity field :math:`\mathbf{u}(x, y, t)`. @@ -137,6 +138,7 @@ def solver(index, ic): tp = mesh_seq.time_partition t_start, t_end = tp.subintervals[index] dt = tp.timesteps[index] + time = mesh_seq.get_time(index) # Initialise the concentration fields Q = mesh_seq.function_spaces["c"][index] @@ -150,7 +152,7 @@ def solver(index, ic): # Compute the velocity field at t_start and assign it to u_ x, y = SpatialCoordinate(mesh_seq[index]) - u_.interpolate(velocity_expression(x, y, t_start)) + u_.interpolate(velocity_expression(x, y, time)) # We pass both the concentration and velocity Functions to get_form form_fields = {"c": (c, c_), "u": (u, u_)} @@ -159,18 +161,17 @@ def solver(index, ic): nlvs = NonlinearVariationalSolver(nlvp, ad_block_tag="c") # Time integrate from t_start to t_end - t = t_start + dt - while t < t_end + 0.5 * dt: - # update the background velocity field at the current timestep - u.interpolate(velocity_expression(x, y, t)) + while float(time) < t_end + 0.5 * dt: + # Update the background velocity field at the current timestep + u.interpolate(velocity_expression(x, y, time)) - # solve the advection equation + # Solve the advection equation nlvs.solve() - # update the 'lagged' concentration and velocity field + # Update the 'lagged' concentration and velocity field c_.assign(c) u_.assign(u) - t += dt + time += dt return {"c": c} diff --git a/demos/burgers-hessian.py b/demos/burgers-hessian.py index 98ff7311..009e6856 100644 --- a/demos/burgers-hessian.py +++ b/demos/burgers-hessian.py @@ -9,12 +9,12 @@ # # As before, we copy over what is now effectively boiler plate to set up our problem. :: -from firedrake import * +import matplotlib.pyplot as plt from animate.adapt import adapt from animate.metric import RiemannianMetric -from goalie import * -import matplotlib.pyplot as plt +from firedrake import * +from goalie import * fields = ["u"] @@ -164,10 +164,8 @@ def adaptor(mesh_seq, solutions): space_time_normalise(metrics, mesh_seq.time_partition, mp) # Adapt each mesh w.r.t. the corresponding metric, provided it hasn't converged - for i, metric in enumerate(metrics): - if not mesh_seq.converged[i]: - mesh_seq[i] = adapt(mesh_seq[i], metric) - complexities.append(metric.complexity()) + complexities = [metric.complexity() for metric in metrics] + mesh_seq.set_meshes(map(adapt, mesh_seq, metrics)) num_dofs = mesh_seq.count_vertices() num_elem = mesh_seq.count_elements() pyrint(f"fixed point iteration {iteration + 1}:") diff --git a/demos/burgers.py b/demos/burgers.py index 5dda6a64..99610dc2 100644 --- a/demos/burgers.py +++ b/demos/burgers.py @@ -21,6 +21,7 @@ # side. See the Firedrake demo for details on the discretisation used. from firedrake import * + from goalie import * # In this problem, we have a single prognostic variable, diff --git a/demos/burgers1.py b/demos/burgers1.py index 170b42be..c819541f 100644 --- a/demos/burgers1.py +++ b/demos/burgers1.py @@ -11,6 +11,7 @@ # :: from firedrake import * + from goalie_adjoint import * # For ease, the field list and functions for obtaining the diff --git a/demos/burgers2.py b/demos/burgers2.py index ada6c7fc..01f7c282 100644 --- a/demos/burgers2.py +++ b/demos/burgers2.py @@ -9,8 +9,8 @@ # Again, begin by importing Goalie with adjoint mode activated. :: from firedrake import * -from goalie_adjoint import * +from goalie_adjoint import * set_log_level(DEBUG) diff --git a/demos/burgers_ee.py b/demos/burgers_ee.py index bc41c8bc..7e6f9655 100644 --- a/demos/burgers_ee.py +++ b/demos/burgers_ee.py @@ -25,6 +25,7 @@ # Goalie. :: from firedrake import * + from goalie_adjoint import * set_log_level(DEBUG) diff --git a/demos/burgers_oo.py b/demos/burgers_oo.py index b9b1a046..cb3c5ac0 100644 --- a/demos/burgers_oo.py +++ b/demos/burgers_oo.py @@ -19,8 +19,8 @@ # :: from firedrake import * -from goalie_adjoint import * +from goalie_adjoint import * set_log_level(DEBUG) diff --git a/demos/burgers_time_integrated.py b/demos/burgers_time_integrated.py index a08e2f50..93b50520 100644 --- a/demos/burgers_time_integrated.py +++ b/demos/burgers_time_integrated.py @@ -9,6 +9,7 @@ # Begin by importing from Goalie and the first Burgers demo. :: from firedrake import * + from goalie_adjoint import * # Redefine the ``get_initial_condition``, ``get_function_spaces``, diff --git a/demos/gray_scott.py b/demos/gray_scott.py index 8aaf6eee..5bad7847 100644 --- a/demos/gray_scott.py +++ b/demos/gray_scott.py @@ -10,6 +10,7 @@ # diffusivities and react with one another nonlinearly. :: from firedrake import * + from goalie_adjoint import * # The problem is defined on a doubly periodic mesh of squares. :: diff --git a/demos/gray_scott_split.py b/demos/gray_scott_split.py index c735e97b..5225d81a 100644 --- a/demos/gray_scott_split.py +++ b/demos/gray_scott_split.py @@ -9,6 +9,7 @@ # equations differ in both the diffusion and reaction terms. :: from firedrake import * + from goalie_adjoint import * # This time, we have two fields instead of one, as well as two function spaces. :: diff --git a/demos/mesh_seq.py b/demos/mesh_seq.py index c0ff4dca..667f5b92 100644 --- a/demos/mesh_seq.py +++ b/demos/mesh_seq.py @@ -12,6 +12,7 @@ # of both Firedrake and Goalie. :: from firedrake import * + from goalie import * # Again, turn debugging mode on to get verbose output. :: @@ -57,7 +58,7 @@ # :figwidth: 90% # :align: center -# In the `next demo <./burgers.py.html>`__, we actually solve a -# PDE using a :class:`MeshSeq`. +# In the `next demo <./ode.py.html>`__, we solve an ordinary differential equation (ODE) +# using a special kind of :class:`MeshSeq`. # # This demo can also be accessed as a `Python script `__. diff --git a/demos/ode.py b/demos/ode.py index b05baca8..9222252e 100644 --- a/demos/ode.py +++ b/demos/ode.py @@ -42,8 +42,8 @@ # First, import from the namespaces of Firedrake and Goalie. :: from firedrake import * -from goalie import * +from goalie import * # Next, create a simple :class:`~.TimeInterval` object to hold information related to # the time discretisation. This is a simplified version of :class:`~.TimePartition`, @@ -331,4 +331,7 @@ def form(index, solutions): # Backward Euler: 3.051757812 # Crank-Nicolson: 2.727412827 # +# In the `next demo <./burgers.py.html>`__, we move on to solve a *partial* differential +# equation (PDE) using a :class:`MeshSeq`. +# # This demo can also be accessed as a `Python script `__. diff --git a/demos/point_discharge2d-goal_oriented.py b/demos/point_discharge2d-goal_oriented.py index f398c110..72c87c48 100644 --- a/demos/point_discharge2d-goal_oriented.py +++ b/demos/point_discharge2d-goal_oriented.py @@ -10,14 +10,14 @@ # We copy over the setup as before. The only difference is that we import from # `goalie_adjoint` rather than `goalie`. :: -from firedrake import * +import matplotlib.colors as mcolors +import matplotlib.pyplot as plt from animate.adapt import adapt from animate.metric import RiemannianMetric -from goalie_adjoint import * -import matplotlib.pyplot as plt -import matplotlib.colors as mcolors +from firedrake import * from matplotlib import ticker +from goalie_adjoint import * fields = ["c"] diff --git a/demos/point_discharge2d-hessian.py b/demos/point_discharge2d-hessian.py index ee0c0b9e..2e0a7242 100644 --- a/demos/point_discharge2d-hessian.py +++ b/demos/point_discharge2d-hessian.py @@ -13,11 +13,11 @@ # adaptation functionality from Firedrake, which can be found in ``firedrake.meshadapt``. # :: -from firedrake import * from animate.adapt import adapt from animate.metric import RiemannianMetric -from goalie import * +from firedrake import * +from goalie import * # We again consider the "point discharge with diffusion" test case from the # `previous demo <./point_discharge2d.py.html>`__, approximating the tracer concentration diff --git a/demos/point_discharge2d.py b/demos/point_discharge2d.py index 32e85378..bc4187b1 100644 --- a/demos/point_discharge2d.py +++ b/demos/point_discharge2d.py @@ -26,6 +26,7 @@ # As always, start by importing Firedrake and Goalie. :: from firedrake import * + from goalie_adjoint import * # We solve the advection-diffusion problem in :math:`\mathbb P1` space. :: diff --git a/demos/solid_body_rotation.py b/demos/solid_body_rotation.py index 889ea49c..c09c98fb 100644 --- a/demos/solid_body_rotation.py +++ b/demos/solid_body_rotation.py @@ -37,6 +37,7 @@ # adjoint mode activated. :: from firedrake import * + from goalie_adjoint import * # For simplicity, we use a :math:`\mathbb P1` space for each @@ -113,9 +114,8 @@ def get_initial_condition(mesh_seq, field="c"): # only the :meth:`get_function_spaces` and :meth:`get_initial_condition` # methods implemented. :: -from firedrake.pyplot import tricontourf import matplotlib.pyplot as plt - +from firedrake.pyplot import tricontourf mesh_seq = MeshSeq( time_partition, diff --git a/demos/solid_body_rotation_split.py b/demos/solid_body_rotation_split.py index b4872439..586c9ea6 100644 --- a/demos/solid_body_rotation_split.py +++ b/demos/solid_body_rotation_split.py @@ -18,7 +18,6 @@ from solid_body_rotation import * - fields = ["bell", "cone", "slot_cyl"] diff --git a/goalie/__init__.py b/goalie/__init__.py index d943c9c4..afa8221a 100644 --- a/goalie/__init__.py +++ b/goalie/__init__.py @@ -7,8 +7,11 @@ from goalie.mesh_seq import * # noqa from goalie.options import * # noqa from goalie.point_seq import * # noqa -from goalie.interpolation import * # noqa +from goalie.function_data import * # noqa from goalie.error_estimation import * # noqa + +from animate.utility import Mesh, VTKFile # noqa + from firedrake.__future__ import interpolate # noqa import numpy as np # noqa diff --git a/goalie/adjoint.py b/goalie/adjoint.py index 6211a04e..710ad471 100644 --- a/goalie/adjoint.py +++ b/goalie/adjoint.py @@ -2,29 +2,31 @@ Drivers for solving adjoint problems on sequences of meshes. """ +from functools import wraps + import firedrake -from firedrake.petsc import PETSc +import numpy as np +from animate.utility import norm from firedrake.adjoint import pyadjoint +from firedrake.petsc import PETSc + from .function_data import AdjointSolutionData -from .interpolation import project +from .log import pyrint from .mesh_seq import MeshSeq from .options import GoalOrientedParameters -from .time_partition import TimePartition -from .utility import AttrDict, norm -from .log import pyrint -from collections.abc import Callable -from functools import wraps -import numpy as np - +from .utility import AttrDict __all__ = ["AdjointMeshSeq", "annotate_qoi"] -def annotate_qoi(get_qoi: Callable) -> Callable: +def annotate_qoi(get_qoi): """ Decorator that ensures QoIs are annotated properly. - Should be applied to the :meth:`~.AdjointMeshSeq.get_qoi` method. + To be applied to the :meth:`~.AdjointMeshSeq.get_qoi` method. + + :arg get_qoi: a function mapping a dictionary of solution data and an integer index + to a QoI function """ @wraps(get_qoi) @@ -72,12 +74,23 @@ class AdjointMeshSeq(MeshSeq): :attr:`~AdjointMeshSeq.J`, which holds the QoI value. """ - def __init__(self, time_partition: TimePartition, initial_meshes: list, **kwargs): - """ - :kwarg get_qoi: a function, with two arguments, a :class:`~.AdjointMeshSeq`, - which returns a function of either one or two variables, corresponding to - either an end time or time integrated quantity of interest, respectively, - as well as an index for the :class:`~.MeshSeq` + def __init__(self, time_partition, initial_meshes, **kwargs): + r""" + :arg time_partition: a partition of the temporal domain + :type time_partition: :class:`~.TimePartition` + :arg initial_meshes: a list of meshes corresponding to the subinterval of the + time partition, or a single mesh to use for all subintervals + :type initial_meshes: :class:`list` or :class:`~.MeshGeometry` + :kwarg get_function_spaces: a function as described in + :meth:`~.MeshSeq.get_function_spaces` + :kwarg get_initial_condition: a function as described in + :meth:`~.MeshSeq.get_initial_condition` + :kwarg get_form: a function as described in :meth:`~.MeshSeq.get_form` + :kwarg get_solver: a function as described in :meth:`~.MeshSeq.get_solver` + :kwarg get_bcs: a function as described in :meth:`~.MeshSeq.get_bcs` + :kwarg parameters: parameters to apply to the mesh adaptation process + :type parameters: :class:`~.AdaptParameters` + :kwarg get_qoi: a function as described in :meth:`~.AdjointMeshSeq.get_qoi` """ if kwargs.get("parameters") is None: kwargs["parameters"] = GoalOrientedParameters() @@ -98,7 +111,7 @@ def __init__(self, time_partition: TimePartition, initial_meshes: list, **kwargs ) self._get_qoi = kwargs.get("get_qoi") self.J = 0 - self.controls = None + self._controls = None self.qoi_values = [] @property @@ -107,26 +120,50 @@ def initial_condition(self): return super().initial_condition @annotate_qoi - def get_qoi(self, solution_map: dict, i: int) -> Callable: + def get_qoi(self, solution_map, subinterval): + """ + Get the function for evaluating the QoI, which has either zero or one arguments, + corresponding to either an end time or time integrated quantity of interest, + respectively. If the QoI has an argument then it is for the current time. + + Signature for the function to be returned: + ``` + :arg t: the current time (for time-integrated QoIs) + :type t: :class:`float` + :return: the QoI as a 0-form + :rtype: :class:`ufl.form.Form` + ``` + + :arg solution_map: a dictionary whose keys are the solution field names and whose + values are the corresponding solutions + :type solution_map: :class:`dict` with :class:`str` keys and values and + :class:`firedrake.function.Function` values + :arg subinterval: the subinterval index + :type subinterval: :class:`int` + :returns: the function for obtaining the QoI + :rtype: see docstring above + """ if self._get_qoi is None: raise NotImplementedError("'get_qoi' is not implemented.") - return self._get_qoi(self, solution_map, i) + return self._get_qoi(self, solution_map, subinterval) @pyadjoint.no_annotations @PETSc.Log.EventDecorator() - def get_checkpoints( - self, solver_kwargs: dict = {}, run_final_subinterval: bool = False - ) -> list: - """ + def get_checkpoints(self, solver_kwargs={}, run_final_subinterval=False): + r""" Solve forward on the sequence of meshes, extracting checkpoints corresponding to the starting fields on each subinterval. The QoI is also evaluated. - :kwarg solver_kwargs: additional keyword arguments which will be passed to the - solver - :kwarg run_final_subinterval: toggle whether to solve the PDE on the final + :kwarg solver_kwargs: additional keyword arguments to be passed to the solver + :type solver_kwargs: :class:`dict` with :class:`str` keys and values which may + take various types + :kwarg run_final_subinterval: if ``True``, the solver is run on the final subinterval + :type run_final_subinterval: :class:`bool` + :returns: checkpoints for each subinterval + :rtype: :class:`list` of :class:`firedrake.function.Function`\s """ # In some cases we run over all subintervals to check the QoI that is computed @@ -145,17 +182,20 @@ def get_checkpoints( return checkpoints @PETSc.Log.EventDecorator() - def get_solve_blocks( - self, field: str, subinterval: int, has_adj_sol: bool = True - ) -> list: - """ + def get_solve_blocks(self, field, subinterval, has_adj_sol=True): + r""" Get all blocks of the tape corresponding to solve steps for prognostic solution field on a given subinterval. :arg field: name of the prognostic solution field + :type field: :class:`str` :arg subinterval: subinterval index + :type subinterval: :class:`int` :kwarg has_adj_sol: if ``True``, only blocks with ``adj_sol`` attributes will be considered + :type has_adj_sol: :class:`bool` + :returns: list of solve blocks + :rtype: :class:`list` of :class:`pyadjoint.block.Block`\s """ solve_blocks = super().get_solve_blocks(field, subinterval) if not has_adj_sol: @@ -176,49 +216,51 @@ def get_solve_blocks( return solve_blocks def _create_solutions(self): + """ + Create the :class:`~.FunctionData` instance for holding solution data. + """ self._solutions = AdjointSolutionData(self.time_partition, self.function_spaces) @PETSc.Log.EventDecorator() def solve_adjoint( self, - solver_kwargs: dict = {}, - adj_solver_kwargs: dict = {}, - get_adj_values: bool = False, - test_checkpoint_qoi: bool = False, - ) -> dict: + solver_kwargs={}, + adj_solver_kwargs={}, + get_adj_values=False, + test_checkpoint_qoi=False, + ): """ Solve an adjoint problem on a sequence of subintervals. - As well as the quantity of interest value, a dictionary of solution fields is - computed, the contents of which give values at all exported timesteps, indexed - first by the field label and then by type. The contents of these nested - dictionaries are lists which are indexed first by subinterval and then by - export. For a given exported timestep, the field types are: - - * ``'forward'``: the forward solution after taking the timestep; - * ``'forward_old'``: the forward solution before taking the timestep (provided - the problem is not steady-state) - * ``'adjoint'``: the adjoint solution after taking the timestep; - * ``'adjoint_next'``: the adjoint solution before taking the timestep - backwards (provided the problem is not steady-state). - - :kwarg solver_kwargs: a dictionary providing parameters to the solver. Any - keyword arguments for the QoI should be included as a subdict with label + As well as the quantity of interest value, solution fields are computed - see + :class:`~.AdjointSolutionData` for more information. + + :kwarg solver_kwargs: parameters for the forward solver, as well as any + parameters for the QoI, which should be included as a sub-dictionary with key 'qoi_kwargs' - :kwarg adj_solver_kwargs: a dictionary providing parameters to the adjoint - solver. - :kwarg get_adj_values: additionally output adjoint actions at exported timesteps + :type solver_kwargs: :class:`dict` with :class:`str` keys and values which may + take various types + :kwarg adj_solver_kwargs: parameters for the adjoint solver + :type adj_solver_kwargs: :class:`dict` with :class:`str` keys and values which + may take various types + :kwarg get_adj_values: if ``True``, adjoint actions are also returned at exported + timesteps + :type get_adj_values: :class:`bool` :kwarg test_checkpoint_qoi: solve over the final subinterval when checkpointing so that the QoI value can be checked across runs - - :return solution: an :class:`~.AttrDict` containing solution fields and their - lagged versions. This can also be accessed as :meth:`solutions`. + :returns: the solution data of the forward and adjoint solves + :rtype: :class:`~.AdjointSolutionData` """ + # TODO #125: Support get_adj_values in AdjointSolutionData + # TODO #126: Separate out qoi_kwargs P = self.time_partition num_subintervals = len(self) solver = self.solver qoi_kwargs = solver_kwargs.get("qoi_kwargs", {}) + # Reinitialise the solution data object + self._create_solutions() + # Solve forward to get checkpoints and evaluate QoI checkpoints = self.get_checkpoints( solver_kwargs=solver_kwargs, @@ -233,9 +275,9 @@ def solve_adjoint( if get_adj_values: for field in self.fields: - self.solutions[field]["adj_value"] = [] + self.solutions.extract(layout="field")[field]["adj_value"] = [] for i, fs in enumerate(self.function_spaces[field]): - self.solutions[field]["adj_value"].append( + self.solutions.extract(layout="field")[field]["adj_value"].append( [ firedrake.Cofunction(fs.dual(), name=f"{field}_adj_value") for j in range(P.num_exports_per_subinterval[i] - 1) @@ -244,20 +286,27 @@ def solve_adjoint( @PETSc.Log.EventDecorator("goalie.AdjointMeshSeq.solve_adjoint.evaluate_fwd") @wraps(solver) - def wrapped_solver(subinterval, ic, **kwargs): + def wrapped_solver(subinterval, initial_condition_map, **kwargs): """ Decorator to allow the solver to stash its initial conditions as controls. :arg subinterval: the subinterval index - :arg ic: the dictionary of initial condition :class:`~.Functions` + :type subinterval: :class:`int` + :arg initial_condition_map: a dictionary of initial conditions, keyed by + field name + :type initial_condition_map: :class:`dict` with :class:`str` keys and + :class:`firedrake.function.Function` values All keyword arguments are passed to the solver. """ - init = AttrDict( - {field: ic[field].copy(deepcopy=True) for field in self.fields} + copy_map = AttrDict( + { + field: initial_condition.copy(deepcopy=True) + for field, initial_condition in initial_condition_map.items() + } ) - self.controls = [pyadjoint.Control(init[field]) for field in self.fields] - return solver(subinterval, init, **kwargs) + self._controls = list(map(pyadjoint.Control, copy_map.values())) + return solver(subinterval, copy_map, **kwargs) # Loop over subintervals in reverse seeds = {} @@ -287,7 +336,7 @@ def wrapped_solver(subinterval, ic, **kwargs): pyadjoint.pause_annotation() else: for field, fs in self.function_spaces.items(): - checkpoint[field].block_variable.adj_value = project( + checkpoint[field].block_variable.adj_value = self._transfer( seeds[field], fs[i] ) @@ -299,7 +348,7 @@ def wrapped_solver(subinterval, ic, **kwargs): # Solve adjoint problem tape = pyadjoint.get_working_tape() with PETSc.Log.Event("goalie.AdjointMeshSeq.solve_adjoint.evaluate_adj"): - m = pyadjoint.enlisting.Enlist(self.controls) + m = pyadjoint.enlisting.Enlist(self._controls) with pyadjoint.stop_annotating(): with tape.marked_nodes(m): tape.evaluate_adj(markings=True) @@ -335,32 +384,32 @@ def wrapped_solver(subinterval, ic, **kwargs): # Update forward and adjoint solution data based on block dependencies # and outputs - sols = self.solutions[field] + solutions = self.solutions.extract(layout="field")[field] for j, block in enumerate(reversed(solve_blocks[::-stride])): # Current forward solution is determined from outputs out = self._output(field, i, block) if out is not None: - sols.forward[i][j].assign(out.saved_output) + solutions.forward[i][j].assign(out.saved_output) # Current adjoint solution is determined from the adj_sol attribute if block.adj_sol is not None: - sols.adjoint[i][j].assign(block.adj_sol) + solutions.adjoint[i][j].assign(block.adj_sol) # Lagged forward solution comes from dependencies dep = self._dependency(field, i, block) if not self.steady and dep is not None: - sols.forward_old[i][j].assign(dep.saved_output) + solutions.forward_old[i][j].assign(dep.saved_output) # Adjoint action also comes from dependencies if get_adj_values and dep is not None: - sols.adj_value[i][j].assign(dep.adj_value) + solutions.adj_value[i][j].assign(dep.adj_value) # The adjoint solution at the 'next' timestep is determined from the # adj_sol attribute of the next solve block if not steady: if (j + 1) * stride < num_solve_blocks: if solve_blocks[(j + 1) * stride].adj_sol is not None: - sols.adjoint_next[i][j].assign( + solutions.adjoint_next[i][j].assign( solve_blocks[(j + 1) * stride].adj_sol ) elif (j + 1) * stride > num_solve_blocks: @@ -372,15 +421,17 @@ def wrapped_solver(subinterval, ic, **kwargs): # The initial timestep of the current subinterval is the 'next' timestep # after the final timestep of the previous subinterval if i > 0 and solve_blocks[0].adj_sol is not None: - project(solve_blocks[0].adj_sol, sols.adjoint_next[i - 1][-1]) + self._transfer( + solve_blocks[0].adj_sol, solutions.adjoint_next[i - 1][-1] + ) # Check non-zero adjoint solution/value - if np.isclose(norm(self.solutions[field].adjoint[i][0]), 0.0): + if np.isclose(norm(solutions.adjoint[i][0]), 0.0): self.warning( f"Adjoint solution for field '{field}' on {self.th(i)}" " subinterval is zero." ) - if get_adj_values and np.isclose(norm(sols.adj_value[i][0]), 0.0): + if get_adj_values and np.isclose(norm(solutions.adj_value[i][0]), 0.0): self.warning( f"Adjoint action for field '{field}' on {self.th(i)}" " subinterval is zero." @@ -388,7 +439,7 @@ def wrapped_solver(subinterval, ic, **kwargs): # Get adjoint action on each subinterval with pyadjoint.stop_annotating(): - for field, control in zip(self.fields, self.controls): + for field, control in zip(self.fields, self._controls): seeds[field] = firedrake.Cofunction( self.function_spaces[field][i].dual() ) @@ -415,9 +466,14 @@ def wrapped_solver(subinterval, ic, **kwargs): return self.solutions @staticmethod - def th(num: int) -> str: + def th(num): """ Convert from cardinal to ordinal. + + :arg num: the cardinal number to convert + :type num: :class:`int` + :returns: the corresponding ordinal number + :rtype: :class:`str` """ end = int(str(num)[-1]) try: @@ -436,6 +492,7 @@ def check_qoi_convergence(self): difference in QoI value being smaller than the specified tolerance. :return: ``True`` if QoI convergence is detected, else ``False`` + :rtype: :class:`bool` """ if not self.check_convergence.any(): self.info( diff --git a/goalie/error_estimation.py b/goalie/error_estimation.py index f8472a9a..3f83b6eb 100644 --- a/goalie/error_estimation.py +++ b/goalie/error_estimation.py @@ -1,19 +1,18 @@ """ Tools to automate goal-oriented error estimation. """ + import firedrake +import ufl from firedrake import Function, FunctionSpace from firedrake.functionspaceimpl import WithGeometry from firedrake.petsc import PETSc -import ufl -from typing import Dict, Optional, Union - __all__ = ["get_dwr_indicator"] @PETSc.Log.EventDecorator() -def form2indicator(F: ufl.form.Form) -> Function: +def form2indicator(F): r""" Given a 0-form, multiply the integrand of each of its integrals by a :math:`\mathbb P0` test function and reassemble to give an element-wise error @@ -23,7 +22,9 @@ def form2indicator(F: ufl.form.Form) -> Function: or :class:`firedrake.ufl_expr.TrialFunction`\s. :arg F: the 0-form + :type F: :class:`ufl.form.Form` :return: the corresponding error indicator field + :rtype: `firedrake.function.Function` """ if not isinstance(F, ufl.form.Form): raise TypeError(f"Expected 'F' to be a Form, not '{type(F)}'.") @@ -59,9 +60,7 @@ def form2indicator(F: ufl.form.Form) -> Function: @PETSc.Log.EventDecorator() -def get_dwr_indicator( - F, adjoint_error: Function, test_space: Optional[Union[WithGeometry, Dict]] = None -) -> Function: +def get_dwr_indicator(F, adjoint_error, test_space=None): r""" Given a 1-form and an approximation of the error in the adjoint solution, compute a dual weighted residual (DWR) error indicator. @@ -73,11 +72,18 @@ def get_dwr_indicator( :class:`firedrake.ufl_expr.TrialFunction`\s. :arg F: the form - :arg adjoint_error: the approximation to the adjoint error, either as a single - :class:`firedrake.function.Function`, or in a dictionary - :kwarg test_space: the - :class:`firedrake.functionspaceimpl.WithGeometry` that the test function lives - in, or an appropriate dictionary + :type F: :class:`ufl.form.Form` + :arg adjoint_error: a dictionary whose keys are field names and whose values are the + approximations to the corresponding components of the adjoint error, or a single + such component + :type adjoint_error: :class:`firedrake.function.Function` or :class:`dict` with + :class:`str` keys and :class:`firedrake.function.Function` values + :kwarg test_space: a dictionary whose keys are field names and whose values are the + test spaces for the corresponding fields, or a single such test space (or + ``None`` to determine the test space(s) automatically) + :type test_space: :class:`firedrake.functionspaceimpl.WithGeometry` + :returns: the DWR indicator + :rtype: :class:`firedrake.function.Function` """ mapping = {} if not isinstance(F, ufl.form.Form): diff --git a/goalie/function_data.py b/goalie/function_data.py index 2b674d61..ae511ab7 100644 --- a/goalie/function_data.py +++ b/goalie/function_data.py @@ -1,10 +1,13 @@ r""" Nested dictionaries of solution data :class:`~.Function`\s. """ + +import abc + import firedrake.function as ffunc import firedrake.functionspace as ffs + from .utility import AttrDict -import abc __all__ = [ "ForwardSolutionData", @@ -18,8 +21,6 @@ class FunctionData(abc.ABC): Abstract base class for classes holding field data. """ - labels = {} - def __init__(self, time_partition, function_spaces): r""" :arg time_partition: the :class:`~.TimePartition` used to discretise the problem @@ -30,10 +31,13 @@ def __init__(self, time_partition, function_spaces): self.time_partition = time_partition self.function_spaces = function_spaces self._data = None + self.labels = self._label_dict[ + "steady" if time_partition.steady else "unsteady" + ] def _create_data(self): - assert self.labels - P = self.time_partition + assert self._label_dict + tp = self.time_partition self._data = AttrDict( { field: AttrDict( @@ -41,57 +45,143 @@ def _create_data(self): label: [ [ ffunc.Function(fs, name=f"{field}_{label}") - for j in range(P.num_exports_per_subinterval[i] - 1) + for j in range(tp.num_exports_per_subinterval[i] - 1) ] for i, fs in enumerate(self.function_spaces[field]) ] - for label in self.labels[field_type] + for label in self._label_dict[field_type] } ) - for field, field_type in zip(P.fields, P.field_types) + for field, field_type in zip(tp.fields, tp.field_types) } ) @property - def data(self): + def _data_by_field(self): + """ + Extract field data array in the default layout: as a doubly-nested dictionary + whose first key is the field name and second key is the field label. Entries + of the doubly-nested dictionary are doubly-nested lists, indexed first by + subinterval and then by export. + """ if self._data is None: self._create_data() return self._data def __getitem__(self, key): - return self.data[key] + return self._data_by_field[key] def items(self): - return self.data.items() - + return self._data_by_field.items() -class SolutionData(FunctionData, abc.ABC): - """ - Abstract base class that defines the API for solution data classes. - """ + @property + def _data_by_label(self): + """ + Extract field data array in an alternative layout: as a doubly-nested dictionary + whose first key is the field label and second key is the field name. Entries + of the doubly-nested dictionary are doubly-nested lists, which retain the default + layout: indexed first by subinterval and then by export. + """ + tp = self.time_partition + return AttrDict( + { + label: AttrDict( + {field: self._data_by_field[field][label] for field in tp.fields} + ) + for label in self.labels + } + ) @property - def solutions(self): - return self.data + def _data_by_subinterval(self): + """ + Extract field data array in an alternative format: as a list indexed by + subinterval. Entries of the list are doubly-nested dictionaries, which retain + the default layout: with the first key being field name and the second key being + the field label. Entries of the doubly-nested dictionaries are lists of field + data, indexed by export. + """ + tp = self.time_partition + return [ + AttrDict( + { + field: AttrDict( + { + label: self._data_by_field[field][label][subinterval] + for label in self.labels + } + ) + for field in tp.fields + } + ) + for subinterval in range(tp.num_subintervals) + ] + + def extract(self, layout="field"): + """ + Extract field data array in a specified layout. + + The default layout is a doubly-nested dictionary whose first key is the field + name and second key is the field label. Entries of the doubly-nested dictionary + are doubly-nested lists, indexed first by subinterval and then by export. That + is: ``data[field][label][subinterval][export]``. + + Choosing a different layout simply promotes the specified variable to first + access: + * ``layout == "label"`` implies ``data[label][field][subinterval][export]`` + * ``layout == "subinterval"`` implies ``data[subinterval][field][label][export]`` + The export index is not promoted because the number of exports may differ across + subintervals. -class ForwardSolutionData(SolutionData): + :kwarg layout: the layout to promote, as described above + :type layout: :class:`str` + """ + if layout == "field": + return self._data_by_field + elif layout == "label": + return self._data_by_label + elif layout == "subinterval": + return self._data_by_subinterval + else: + raise ValueError(f"Layout type '{layout}' not recognised.") + + +class ForwardSolutionData(FunctionData): """ Class representing solution data for general forward problems. + + For a given exported timestep, the field types are: + + * ``'forward'``: the forward solution after taking the timestep; + * ``'forward_old'``: the forward solution before taking the timestep (provided + the problem is not steady-state). """ def __init__(self, *args, **kwargs): - self.labels = {"steady": ("forward",), "unsteady": ("forward", "forward_old")} + self._label_dict = { + "steady": ("forward",), + "unsteady": ("forward", "forward_old"), + } super().__init__(*args, **kwargs) -class AdjointSolutionData(SolutionData): +class AdjointSolutionData(FunctionData): """ Class representing solution data for general adjoint problems. + + For a given exported timestep, the field types are: + + * ``'forward'``: the forward solution after taking the timestep; + * ``'forward_old'``: the forward solution before taking the timestep (provided + the problem is not steady-state) + * ``'adjoint'``: the adjoint solution after taking the timestep; + * ``'adjoint_next'``: the adjoint solution before taking the timestep + backwards (provided the problem is not steady-state). """ def __init__(self, *args, **kwargs): - self.labels = { + self._label_dict = { "steady": ("forward", "adjoint"), "unsteady": ("forward", "forward_old", "adjoint", "adjoint_next"), } @@ -112,25 +202,52 @@ def __init__(self, time_partition, meshes): in time :arg meshes: the list of meshes used to discretise the problem in space """ - self.labels = { + self._label_dict = { field_type: ("error_indicator",) for field_type in ("steady", "unsteady") } - P0_spaces = [ffs.FunctionSpace(mesh, "DG", 0) for mesh in meshes] super().__init__( - time_partition, {key: P0_spaces for key in time_partition.fields} + time_partition, + { + key: [ffs.FunctionSpace(mesh, "DG", 0) for mesh in meshes] + for key in time_partition.fields + }, ) - def _create_data(self): - assert all(len(labels) == 1 for labels in self.labels.values()) - super()._create_data() - P = self.time_partition - self._data = AttrDict( + @property + def _data_by_field(self): + """ + Extract indicator data array in the default layout: as a dictionary keyed with + the field name. Entries of the dictionary are doubly-nested lists, indexed first + by subinterval and then by export. + """ + if self._data is None: + self._create_data() + return AttrDict( { - field: self.data[field][self.labels[field_type][0]] - for field, field_type in zip(P.fields, P.field_types) + field: self._data[field]["error_indicator"] + for field in self.time_partition.fields } ) @property - def indicators(self): - return self.data + def _data_by_label(self): + """ + For indicator data there is only one field label (``"error_indicator"``), so + this method just delegates to :meth:`~._data_by_field`. + """ + return self._data_by_field + + @property + def _data_by_subinterval(self): + """ + Extract indicator data array in an alternative format: as a list indexed by + subinterval. Entries of the list are dictionaries, keyed by field label. + Entries of the dictionaries are lists of field data, indexed by export. + """ + tp = self.time_partition + return [ + AttrDict( + {field: self._data_by_field[field][subinterval] for field in tp.fields} + ) + for subinterval in range(tp.num_subintervals) + ] diff --git a/goalie/go_mesh_seq.py b/goalie/go_mesh_seq.py index ab537ff1..3fe72346 100644 --- a/goalie/go_mesh_seq.py +++ b/goalie/go_mesh_seq.py @@ -2,19 +2,18 @@ Drivers for goal-oriented error estimation on sequences of meshes. """ +from collections.abc import Iterable +from inspect import signature + +import numpy as np +import ufl +from firedrake import Function, FunctionSpace, MeshHierarchy, TransferManager +from firedrake.petsc import PETSc + from .adjoint import AdjointMeshSeq from .error_estimation import get_dwr_indicator from .function_data import IndicatorData from .log import pyrint -from .utility import AttrDict -from firedrake import Function, FunctionSpace, MeshHierarchy, TransferManager, project -from firedrake.petsc import PETSc -from collections.abc import Callable, Iterable -import numpy as np -from typing import Tuple -import ufl -from inspect import signature - __all__ = ["GoalOrientedMeshSeq"] @@ -29,18 +28,22 @@ def __init__(self, *args, **kwargs): self.estimator_values = [] @PETSc.Log.EventDecorator() - def get_enriched_mesh_seq( - self, enrichment_method: str = "p", num_enrichments: int = 1 - ) -> AdjointMeshSeq: + def get_enriched_mesh_seq(self, enrichment_method="p", num_enrichments=1): """ Construct a sequence of globally enriched spaces. - Currently, global enrichment may be achieved using one of: - * h-refinement (``enrichment_method = 'h'``); - * p-refinement (``enrichment_method = 'p'``). - - The number of refinements may be controlled by the keyword argument - ``num_enrichments``. + The following global enrichment methods are supported: + * h-refinement (``enrichment_method='h'``) - refine each mesh element + uniformly in each direction; + * p-refinement (``enrichment_method='p'``) - increase the function space + polynomial order by one globally. + + :kwarg enrichment_method: the method for enriching the mesh sequence + :type enrichment_method: :class:`str` + :kwarg num_enrichments: the number of enrichments to apply + :type num_enrichments: :class:`int` + :returns: the enriched mesh sequence + :type: the type is inherited from the parent mesh sequence """ if enrichment_method not in ("h", "p"): raise ValueError(f"Enrichment method '{enrichment_method}' not supported.") @@ -58,7 +61,7 @@ def get_enriched_mesh_seq( meshes = self.meshes # Construct object to hold enriched spaces - enriched_mesh_seq = self.__class__( + enriched_mesh_seq = type(self)( self.time_partition, meshes, get_function_spaces=self._get_function_spaces, @@ -70,7 +73,7 @@ def get_enriched_mesh_seq( qoi_type=self.qoi_type, parameters=self.params, ) - enriched_mesh_seq.update_function_spaces() + enriched_mesh_seq._update_function_spaces() # Apply p-refinement if enrichment_method == "p": @@ -88,18 +91,32 @@ def get_enriched_mesh_seq( @staticmethod def _get_transfer_function(enrichment_method): + """ + Get the function for transferring function data between a mesh sequence and its + enriched counterpart. + + :arg enrichment_method: the enrichment method used to generate the counterpart + - see :meth:`~.GoalOrientedMeshSeq.get_enriched_mesh_seq` for the supported + enrichment methods + :type enrichment_method: :class:`str` + :returns: the function for mapping function data between mesh sequences + """ if enrichment_method == "h": return TransferManager().prolong else: return lambda source, target: target.interpolate(source) def _create_indicators(self): + """ + Create the :class:`~.FunctionData` instance for holding error indicator data. + """ self._indicators = IndicatorData(self.time_partition, self.meshes) @property def indicators(self): """ - Arrays holding exported error indicators. + :returns: the error indicator data object + :rtype: :class:`~.IndicatorData` """ if not hasattr(self, "_indicators"): self._create_indicators() @@ -107,29 +124,44 @@ def indicators(self): @PETSc.Log.EventDecorator() def indicate_errors( - self, - enrichment_kwargs: dict = {}, - adj_kwargs: dict = {}, - indicator_fn: Callable = get_dwr_indicator, - ) -> Tuple[dict, AttrDict]: + self, enrichment_kwargs={}, solver_kwargs={}, indicator_fn=get_dwr_indicator + ): """ Compute goal-oriented error indicators for each subinterval based on solving the adjoint problem in a globally enriched space. :kwarg enrichment_kwargs: keyword arguments to pass to the global enrichment - method - :kwarg adj_kwargs: keyword arguments to pass to the adjoint solver - :kwarg indicator_fn: function for error indication, which takes the form, - adjoint error and enriched space(s) as arguments + method - see :meth:`~.GoalOrientedMeshSeq.get_enriched_mesh_seq` for the + supported enrichment methods and options + :type enrichment_kwargs: :class:`dict` with :class:`str` keys and values which + may take various types + :kwarg solver_kwargs: parameters for the forward solver, as well as any + parameters for the QoI, which should be included as a sub-dictionary with key + 'qoi_kwargs' + :type solver_kwargs: :class:`dict` with :class:`str` keys and values which may + take various types + :kwarg indicator_fn: function which maps the form, adjoint error and enriched + space(s) as arguments to the error indicator + :class:`firedrake.function.Function` + :returns: solution and indicator data objects + :rtype1: :class:`~.AdjointSolutionData + :rtype2: :class:`~.IndicatorData """ enrichment_kwargs.setdefault("enrichment_method", "p") enrichment_kwargs.setdefault("num_enrichments", 1) enriched_mesh_seq = self.get_enriched_mesh_seq(**enrichment_kwargs) transfer = self._get_transfer_function(enrichment_kwargs["enrichment_method"]) + # Reinitialise the error indicator data object + self._create_indicators() + # Solve the forward and adjoint problems on the MeshSeq and its enriched version - self.solve_adjoint(**adj_kwargs) - enriched_mesh_seq.solve_adjoint(**adj_kwargs) + self.solve_adjoint(**solver_kwargs) + enriched_mesh_seq.solve_adjoint(**solver_kwargs) + + tp = self.time_partition + # check whether get_form contains a kwarg indicating time-dependent constants + timedep_const = "err_ind_time" in signature(enriched_mesh_seq.form).parameters tp = self.time_partition # check whether get_form contains a kwarg indicating time-dependent constants @@ -140,6 +172,8 @@ def indicate_errors( ADJ_NEXT = "adjoint" if self.steady else "adjoint_next" P0_spaces = [FunctionSpace(mesh, "DG", 0) for mesh in self] for i, mesh in enumerate(self): + time = self.get_time(i) + # Get Functions u, u_, u_star, u_star_next, u_star_e = {}, {}, {}, {}, {} enriched_spaces = { @@ -166,15 +200,17 @@ def indicate_errors( # Loop over each strongly coupled field for f in self.fields: # Loop over each timestep - for j in range(len(self.solutions[f]["forward"][i])): + solutions = self.solutions.extract(layout="field") + indicators = self.indicators.extract(layout="field") + for j in range(len(solutions[f]["forward"][i])): # Update fields if timedep_const: # recompile the form with updated time-dependent constants - time = ( + time.assign( tp.subintervals[i][0] + (j + 1) * tp.timesteps[i] * tp.num_timesteps_per_export[i] ) - forms = enriched_mesh_seq.form(i, mapping, err_ind_time=time) + forms = enriched_mesh_seq.form(i, mapping) transfer(self.solutions[f][FWD][i][j], u[f]) transfer(self.solutions[f][FWD_OLD][i][j], u_[f]) transfer(self.solutions[f][ADJ][i][j], u_star[f]) @@ -194,20 +230,23 @@ def indicate_errors( # Evaluate error indicator indi_e = indicator_fn(forms[f], u_star_e[f]) - # Project back to the base space - indi = project(indi_e, P0_spaces[i]) + # Transfer back to the base space + indi = self._transfer(indi_e, P0_spaces[i]) indi.interpolate(abs(indi)) - self.indicators[f][i][j].interpolate(ufl.max_value(indi, 1.0e-16)) + indicators[f][i][j].interpolate(ufl.max_value(indi, 1.0e-16)) return self.solutions, self.indicators @PETSc.Log.EventDecorator() - def error_estimate(self, absolute_value: bool = False) -> float: + def error_estimate(self, absolute_value=False): r""" - Deduce the error estimator value associated with error indicator fields defined over - a :class:`~.MeshSeq`. + Deduce the error estimator value associated with error indicator fields defined + over the mesh sequence. - :kwarg absolute_value: toggle whether to take the modulus on each element + :kwarg absolute_value: if ``True``, the modulus is taken on each element + :type absolute_value: :class:`bool` + :returns: the error estimator value + :rtype: :class:`float` """ assert isinstance(self.indicators, IndicatorData) if not isinstance(absolute_value, bool): @@ -237,6 +276,7 @@ def check_estimator_convergence(self): difference in error estimator value being smaller than the specified tolerance. :return: ``True`` if estimator convergence is detected, else ``False`` + :rtype: :class:`bool` """ if not self.check_convergence.any(): self.info( @@ -257,32 +297,42 @@ def check_estimator_convergence(self): @PETSc.Log.EventDecorator() def fixed_point_iteration( self, - adaptor: Callable, - enrichment_kwargs: dict = {}, - adaptor_kwargs: dict = {}, - adj_kwargs: dict = {}, - indicator_fn: Callable = get_dwr_indicator, - **kwargs, + adaptor, + update_params=None, + enrichment_kwargs={}, + adaptor_kwargs={}, + solver_kwargs={}, + indicator_fn=get_dwr_indicator, ): r""" - Apply goal-oriented mesh adaptation using a fixed point iteration loop. - - :arg adaptor: function for adapting the mesh sequence. Its arguments are the - :class:`~.MeshSeq` instance, the dictionary of solution - :class:`firedrake.function.Function`\s and the list of error indicators. - It should return ``True`` if the convergence criteria checks are to be - skipped for this iteration. Otherwise, it should return ``False``. - :kwarg update_params: function for updating :attr:`~.GoalOrientedMeshSeq.params` - at each iteration. Its arguments are the parameter class and the fixed point - iteration number + Apply goal-oriented mesh adaptation using a fixed point iteration loop approach. + + :arg adaptor: function for adapting the mesh sequence. Its arguments are the mesh + sequence and the solution and indicator data objects. It should return + ``True`` if the convergence criteria checks are to be skipped for this + iteration. Otherwise, it should return ``False``. + :kwarg update_params: function for updating :attr:`~.MeshSeq.params` at each + iteration. Its arguments are the parameter class and the fixed point + iteration :kwarg enrichment_kwargs: keyword arguments to pass to the global enrichment method - :kwarg adaptor_kwargs: a dictionary providing parameters to the adaptor - :kwarg adj_kwargs: keyword arguments to pass to the adjoint solver - :kwarg indicator_fn: function for error indication, which takes the form, adjoint - error and enriched space(s) as arguments + :type enrichment_kwargs: :class:`dict` with :class:`str` keys and values which + may take various types + :kwarg solver_kwargs: parameters to pass to the solver + :type solver_kwargs: :class:`dict` with :class:`str` keys and values which may + take various types + :kwarg adaptor_kwargs: parameters to pass to the adaptor + :type adaptor_kwargs: :class:`dict` with :class:`str` keys and values which may + take various types + :kwarg indicator_fn: function which maps the form, adjoint error and enriched + space(s) as arguments to the error indicator + :class:`firedrake.function.Function` + :returns: solution and indicator data objects + :rtype1: :class:`~.AdjointSolutionData + :rtype2: :class:`~.IndicatorData """ - update_params = kwargs.get("update_params") + # TODO #124: adaptor no longer needs solution and indicator data to be passed + # explicitly self.element_counts = [self.count_elements()] self.vertex_counts = [self.count_vertices()] self.qoi_values = [] @@ -295,18 +345,15 @@ def fixed_point_iteration( update_params(self.params, self.fp_iteration) # Indicate errors over all meshes - self._create_solutions() - self._create_indicators() self.indicate_errors( enrichment_kwargs=enrichment_kwargs, - adj_kwargs=adj_kwargs, + solver_kwargs=solver_kwargs, indicator_fn=indicator_fn, ) # Check for QoI convergence - # TODO: Put this check inside the adjoint solve as - # an optional return condition so that we - # can avoid unnecessary extra solves + # TODO #23: Put this check inside the adjoint solve as an optional return + # condition so that we can avoid unnecessary extra solves self.qoi_values.append(self.J) qoi_converged = self.check_qoi_convergence() if self.params.convergence_criteria == "any" and qoi_converged: diff --git a/goalie/interpolation.py b/goalie/interpolation.py deleted file mode 100644 index 46c65bc2..00000000 --- a/goalie/interpolation.py +++ /dev/null @@ -1,150 +0,0 @@ -""" -Driver functions for mesh-to-mesh data transfer. -""" -from .utility import assemble_mass_matrix, cofunction2function, function2cofunction -import firedrake -from firedrake.functionspaceimpl import WithGeometry -from firedrake.petsc import PETSc -from petsc4py import PETSc as petsc4py - - -__all__ = ["project"] - - -def project(source, target_space, **kwargs): - r""" - Overload :func:`firedrake.projection.project` to account for the case of two mixed - function spaces defined on different meshes and for the adjoint projection operator - when applied to :class:`firedrake.cofunction.Cofunction`\s. - - Extra keyword arguments are passed to :func:`firedrake.projection.project`. - - :arg source: the :class:`firedrake.function.Function` or - :class:`firedrake.cofunction.Cofunction` to be projected - :arg target_space: the :class:`firedrake.functionspaceimpl.FunctionSpace` which we - seek to project into, or the :class:`firedrake.function.Function` or - :class:`firedrake.cofunction.Cofunction` to use as the target - """ - if not isinstance(source, (firedrake.Function, firedrake.Cofunction)): - raise NotImplementedError( - "Can only currently project Functions and Cofunctions." - ) - if isinstance(target_space, WithGeometry): - target = firedrake.Function(target_space) - elif isinstance(target_space, (firedrake.Cofunction, firedrake.Function)): - target = target_space - else: - raise TypeError( - "Second argument must be a FunctionSpace, Function, or Cofunction." - ) - if isinstance(source, firedrake.Cofunction): - return _project_adjoint(source, target, **kwargs) - elif source.function_space() == target.function_space(): - return target.assign(source) - else: - return _project(source, target, **kwargs) - - -@PETSc.Log.EventDecorator("goalie.interpolation.project") -def _project(source, target, **kwargs): - """ - Apply a mesh-to-mesh conservative projection to some source - :class:`firedrake.function.Function`, mapping to a target - :class:`firedrake.function.Function`. - - This function extends to the case of mixed spaces. - - Extra keyword arguments are passed to Firedrake's - :func:`firedrake.projection.project`` function. - - :arg source: the `Function` to be projected - :arg target: the `Function` which we seek to project onto - """ - Vs = source.function_space() - Vt = target.function_space() - if hasattr(Vs, "num_sub_spaces"): - if not hasattr(Vt, "num_sub_spaces"): - raise ValueError( - "Source space has multiple components but target space does not." - ) - if Vs.num_sub_spaces() != Vt.num_sub_spaces(): - raise ValueError( - "Inconsistent numbers of components in source and target spaces:" - f" {Vs.num_sub_spaces()} vs. {Vt.num_sub_spaces()}." - ) - elif hasattr(Vt, "num_sub_spaces"): - raise ValueError( - "Target space has multiple components but source space does not." - ) - assert isinstance(target, firedrake.Function) - if hasattr(Vt, "num_sub_spaces"): - for s, t in zip(source.subfunctions, target.subfunctions): - t.project(s, **kwargs) - else: - target.project(source, **kwargs) - return target - - -@PETSc.Log.EventDecorator("goalie.interpolation.project_adjoint") -def _project_adjoint(target_b, source_b, **kwargs): - """ - Apply the adjoint of a mesh-to-mesh conservative projection to some seed - :class:`firedrake.cofunction.Cofunction`, mapping to an output - :class:`firedrake.cofunction.Cofunction`. - - The notation used here is in terms of the adjoint of standard projection. - However, this function may also be interpreted as a projector in its own right, - mapping ``target_b`` to ``source_b``. - - Extra keyword arguments are passed to :func:`firedrake.projection.project`. - - :arg target_b: seed :class:`firedrake.cofunction.Cofunction` from the target space - of the forward projection - :arg source_b: the :class:`firedrake.cofunction.Cofunction` from the source space - of the forward projection - """ - from firedrake.supermeshing import assemble_mixed_mass_matrix - - # Map to Functions to apply the adjoint projection - if not isinstance(target_b, firedrake.Function): - target_b = cofunction2function(target_b) - if not isinstance(source_b, firedrake.Function): - source_b = cofunction2function(source_b) - - Vt = target_b.function_space() - Vs = source_b.function_space() - if hasattr(Vs, "num_sub_spaces"): - if not hasattr(Vt, "num_sub_spaces"): - raise ValueError( - "Source space has multiple components but target space does not." - ) - if Vs.num_sub_spaces() != Vt.num_sub_spaces(): - raise ValueError( - "Inconsistent numbers of components in target and source spaces:" - f" {Vs.num_sub_spaces()} vs. {Vt.num_sub_spaces()}." - ) - target_b_split = target_b.subfunctions - source_b_split = source_b.subfunctions - elif hasattr(Vt, "num_sub_spaces"): - raise ValueError( - "Target space has multiple components but source space does not." - ) - else: - target_b_split = [target_b] - source_b_split = [source_b] - - # Apply adjoint projection operator to each component - if Vs == Vt: - source_b.assign(target_b) - else: - for i, (t_b, s_b) in enumerate(zip(target_b_split, source_b_split)): - ksp = petsc4py.KSP().create() - ksp.setOperators(assemble_mass_matrix(t_b.function_space())) - mixed_mass = assemble_mixed_mass_matrix(Vt[i], Vs[i]) - with t_b.dat.vec_ro as tb, s_b.dat.vec_wo as sb: - residual = tb.copy() - ksp.solveTranspose(tb, residual) - mixed_mass.mult(residual, sb) # NOTE: already transposed above - - # Map back to a Cofunction - return function2cofunction(source_b) diff --git a/goalie/log.py b/goalie/log.py index 6a9420d5..6ff63445 100644 --- a/goalie/log.py +++ b/goalie/log.py @@ -4,10 +4,11 @@ Code mostly copied from `the Thetis project `__. """ -import firedrake + import logging -from logging import DEBUG, WARNING, ERROR +from logging import DEBUG, ERROR, WARNING +import firedrake __all__ = [ "logger", @@ -25,9 +26,15 @@ ] -def get_new_logger( - name: str, fmt: str = "%(levelname)s %(message)s" -) -> logging.StreamHandler: +def get_new_logger(name, fmt="%(levelname)s %(message)s"): + """ + :arg name: the name for the logger + :type name: :class:`str` + :kwarg fmt: format string to use + :type fmt: :class:`str` + :returns: logger instance + :rtype: :class:`logging.StreamHandler` + """ logger = logging.getLogger(name) for handler in logger.handlers: logger.removeHandler(handler) diff --git a/goalie/math.py b/goalie/math.py index c06214a2..bd8971dd 100644 --- a/goalie/math.py +++ b/goalie/math.py @@ -2,7 +2,6 @@ import ufl import ufl.core.expr - __all__ = ["bessi0", "bessk0"] diff --git a/goalie/mesh_seq.py b/goalie/mesh_seq.py index 7df3e374..0837caab 100644 --- a/goalie/mesh_seq.py +++ b/goalie/mesh_seq.py @@ -2,24 +2,24 @@ Sequences of meshes corresponding to a :class:`~.TimePartition`. """ +from collections.abc import Iterable + import firedrake +import firedrake.function as ffunc +import firedrake.functionspace as ffs +import numpy as np +from animate.interpolation import transfer +from animate.quality import QualityMeasure +from animate.utility import Mesh from firedrake.adjoint import pyadjoint from firedrake.adjoint_utils.solving import get_solve_blocks from firedrake.petsc import PETSc from firedrake.pyplot import triplot + from .function_data import ForwardSolutionData -from .interpolation import project -from .log import pyrint, debug, warning, info, logger, DEBUG +from .log import DEBUG, debug, info, logger, pyrint, warning from .options import AdaptParameters -from animate.quality import QualityMeasure -from .time_partition import TimePartition -from .utility import AttrDict, Mesh -from collections import OrderedDict -from collections.abc import Callable, Iterable -import matplotlib -import numpy as np -from typing import Tuple - +from .utility import AttrDict __all__ = ["MeshSeq"] @@ -31,24 +31,25 @@ class MeshSeq: """ @PETSc.Log.EventDecorator() - def __init__(self, time_partition: TimePartition, initial_meshes: list, **kwargs): + def __init__(self, time_partition, initial_meshes, **kwargs): r""" - :arg time_partition: the :class:`~.TimePartition` which partitions the temporal - domain - :arg initial_meshes: list of meshes corresponding to the subintervals of the - :class:`~.TimePartition`, or a single mesh to use for all subintervals - :kwarg get_function_spaces: a function, whose only argument is a - :class:`~.MeshSeq`, which constructs prognostic - :class:`firedrake.functionspaceimpl.FunctionSpace`\s for each subinterval - :kwarg get_initial_condition: a function, whose only argument is a - :class:`~.MeshSeq`, which specifies initial conditions on the first mesh - :kwarg get_form: a function, whose only argument is a :class:`~.MeshSeq`, which - returns a function that generates the PDE weak form - :kwarg get_solver: a function, whose only argument is a :class:`~.MeshSeq`, - which returns a function that integrates initial data over a subinterval - :kwarg get_bcs: a function, whose only argument is a :class:`~.MeshSeq`, which - returns a function that determines any Dirichlet boundary conditions - :kwarg parameters: :class:`~.AdaptParameters` instance + :arg time_partition: a partition of the temporal domain + :type time_partition: :class:`~.TimePartition` + :arg initial_meshes: a list of meshes corresponding to the subinterval of the + time partition, or a single mesh to use for all subintervals + :type initial_meshes: :class:`list` or :class:`~.MeshGeometry` + :kwarg get_function_spaces: a function as described in + :meth:`~.MeshSeq.get_function_spaces` + :kwarg get_initial_condition: a function as described in + :meth:`~.MeshSeq.get_initial_condition` + :kwarg get_form: a function as described in :meth:`~.MeshSeq.get_form` + :kwarg get_solver: a function as described in :meth:`~.MeshSeq.get_solver` + :kwarg get_bcs: a function as described in :meth:`~.MeshSeq.get_bcs` + :kwarg transfer_method: the method to use for transferring fields between + meshes. Options are "interpolate" (default) and "project" + :type transfer_method: :class:`str` + :kwarg parameters: parameters to apply to the mesh adaptation process + :type parameters: :class:`~.AdaptParameters` """ self.time_partition = time_partition self.fields = time_partition.fields @@ -65,6 +66,7 @@ def __init__(self, time_partition: TimePartition, initial_meshes: list, **kwargs self._get_form = kwargs.get("get_form") self._get_solver = kwargs.get("get_solver") self._get_bcs = kwargs.get("get_bcs") + self._transfer_method = kwargs.get("transfer_method", "interpolate") self.params = kwargs.get("parameters") self.steady = time_partition.steady self.check_convergence = np.array([True] * len(self), dtype=bool) @@ -74,11 +76,11 @@ def __init__(self, time_partition: TimePartition, initial_meshes: list, **kwargs self.params = AdaptParameters() self.sections = [{} for mesh in self] - def __str__(self) -> str: + def __str__(self): return f"{[str(mesh) for mesh in self.meshes]}" - def __repr__(self) -> str: - name = self.__class__.__name__ + def __repr__(self): + name = type(self).__name__ if len(self) == 1: return f"{name}([{repr(self.meshes[0])}])" elif len(self) == 2: @@ -86,51 +88,99 @@ def __repr__(self) -> str: else: return f"{name}([{repr(self.meshes[0])}, ..., {repr(self.meshes[-1])}])" - def debug(self, msg: str): - debug(f"{self.__class__.__name__}: {msg}") + def debug(self, msg): + """ + Print a ``debug`` message. - def warning(self, msg: str): - warning(f"{self.__class__.__name__}: {msg}") + :arg msg: the message to print + :type msg: :class:`str` + """ + debug(f"{type(self).__name__}: {msg}") - def info(self, msg: str): - info(f"{self.__class__.__name__}: {msg}") + def warning(self, msg): + """ + Print a ``warning`` message. - def __len__(self) -> int: + :arg msg: the message to print + :type msg: :class:`str` + """ + warning(f"{type(self).__name__}: {msg}") + + def info(self, msg): + """ + Print an ``info`` level message. + + :arg msg: the message to print + :type msg: :class:`str` + """ + info(f"{type(self).__name__}: {msg}") + + def __len__(self): return len(self.meshes) - def __getitem__(self, i: int) -> firedrake.MeshGeometry: - return self.meshes[i] + def __getitem__(self, subinterval): + """ + :arg subinterval: a subinterval index + :type subinterval: :class:`int` + :returns: the corresponding mesh + :rtype: :class:`firedrake.MeshGeometry` + """ + return self.meshes[subinterval] - def __setitem__( - self, i: int, mesh: firedrake.MeshGeometry - ) -> firedrake.MeshGeometry: - self.meshes[i] = mesh + def __setitem__(self, subinterval, mesh): + """ + :arg subinterval: a subinterval index + :type subinterval: :class:`int` + :arg mesh: the mesh to use for that subinterval + :type subinterval: :class:`firedrake.MeshGeometry` + """ + self.meshes[subinterval] = mesh - def count_elements(self) -> list: - return [mesh.num_cells() for mesh in self] # TODO: make parallel safe + def count_elements(self): + r""" + Count the number of elements in each mesh in the sequence. - def count_vertices(self) -> list: - return [mesh.num_vertices() for mesh in self] # TODO: make parallel safe + :returns: list of element counts + :rtype: :class:`list` of :class:`int`\s + """ + return [mesh.num_cells() for mesh in self] # TODO #123: make parallel safe + + def count_vertices(self): + r""" + Count the number of vertices in each mesh in the sequence. + + :returns: list of vertex counts + :rtype: :class:`list` of :class:`int`\s + """ + return [mesh.num_vertices() for mesh in self] # TODO #123: make parallel safe - def reset_counts(self): + def _reset_counts(self): + """ + Reset the lists of element and vertex counts. + """ self.element_counts = [self.count_elements()] self.vertex_counts = [self.count_vertices()] def set_meshes(self, meshes): - """ - Update the meshes associated with the :class:`~.MeshSeq`, as well as the - associated attributes. + r""" + Set all meshes in the sequence and deduce various properties. - :arg meshes: mesh or list of meshes to use in the sequence + :arg meshes: list of meshes to use in the sequence, or a single mesh to use for + all subintervals + :type meshes: :class:`list` of :class:`firedrake.MeshGeometry`\s or + :class:`firedrake.MeshGeometry` """ - if not isinstance(meshes, Iterable): + # TODO #122: Refactor to use the set method + if isinstance(meshes, map): + meshes = list(meshes) + elif not isinstance(meshes, Iterable): meshes = [Mesh(meshes) for subinterval in self.subintervals] self.meshes = meshes dim = np.array([mesh.topological_dimension() for mesh in meshes]) if dim.min() != dim.max(): raise ValueError("Meshes must all have the same topological dimension.") self.dim = dim.min() - self.reset_counts() + self._reset_counts() if logger.level == DEBUG: for i, mesh in enumerate(meshes): nc = mesh.num_cells() @@ -142,18 +192,38 @@ def set_meshes(self, meshes): f"{i}: {nc:7d} cells, {nv:7d} vertices, max aspect ratio {mar:.2f}" ) debug(100 * "-") + self._time = [ + ffunc.Function(ffs.FunctionSpace(mesh, "R", 0)) for mesh in meshes + ] - def plot( - self, **kwargs - ) -> Tuple[matplotlib.figure.Figure, matplotlib.axes._axes.Axes]: + def get_time(self, subinterval): + """ + Get the time :class:`~.Function` associated with a given subinterval, + initialised to the value at the start of the subinterval, plus one timestep. + + :arg subinterval: the subinterval index + :type subinterval: :class:`int` + :return: the associated $R$-space time Function + :rtype: :class:`~.Function` + """ + start_time = self.time_partition[subinterval].start_time + dt = self.time_partition.timesteps[subinterval] + self._time[subinterval].assign(start_time + dt) + return self._time[subinterval] + + def plot(self, fig=None, axes=None, **kwargs): """ Plot the meshes comprising a 2D :class:`~.MeshSeq`. - :kwarg fig: matplotlib figure - :kwarg axes: matplotlib axes - :kwargs: parameters to pass to :func:`firedrake.pyplot.triplot` - function - :return: matplotlib figure and axes for the plots + :kwarg fig: matplotlib figure to use + :type fig: :class:`matplotlib.figure.Figure` + :kwarg axes: matplotlib axes to use + :type axes: :class:`matplotlib.axes._axes.Axes` + :returns: matplotlib figure and axes for the plots + :rtype1: :class:`matplotlib.figure.Figure` + :rtype2: :class:`matplotlib.axes._axes.Axes` + + All keyword arguments are passed to :func:`firedrake.pyplot.triplot`. """ from matplotlib.pyplot import subplots @@ -161,8 +231,6 @@ def plot( raise ValueError("MeshSeq plotting only supported in 2D") # Process kwargs - fig = kwargs.pop("fig", None) - axes = kwargs.pop("axes", None) interior_kw = {"edgecolor": "k"} interior_kw.update(kwargs.pop("interior_kw", {})) boundary_kw = {"edgecolor": "k"} @@ -191,12 +259,30 @@ def plot( axes = axes[0] return fig, axes - def get_function_spaces(self, mesh: firedrake.MeshGeometry) -> Callable: + def get_function_spaces(self, mesh): + """ + Construct the function spaces corresponding to each field, for a given mesh. + + :arg mesh: the mesh to base the function spaces on + :type mesh: :class:`firedrake.mesh.MeshGeometry` + :returns: a dictionary whose keys are field names and whose values are the + corresponding function spaces + :rtype: :class:`dict` with :class:`str` keys and + :class:`firedrake.functionspaceimpl.FunctionSpace` values + """ if self._get_function_spaces is None: raise NotImplementedError("'get_function_spaces' needs implementing.") return self._get_function_spaces(mesh) - def get_initial_condition(self) -> dict: + def get_initial_condition(self): + r""" + Get the initial conditions applied on the first mesh in the sequence. + + :returns: the dictionary, whose keys are field names and whose values are the + corresponding initial conditions applied + :rtype: :class:`dict` with :class:`str` keys and + :class:`firedrake.function.Function` values + """ if self._get_initial_condition is not None: return self._get_initial_condition(self) return { @@ -204,22 +290,101 @@ def get_initial_condition(self) -> dict: for field, fs in self.function_spaces.items() } - def get_form(self) -> Callable: + def get_form(self): + """ + Get the function mapping a subinterval index and a solution dictionary to a + dictionary containing parts of the PDE weak form corresponding to each solution + component. + + Signature for the function to be returned: + ``` + :arg index: the subinterval index + :type index: :class:`int` + :arg solutions: map from fields to tuples of current and previous solutions + :type solutions: :class:`dict` with :class:`str` keys and :class:`tuple` values + :return: map from fields to the corresponding forms + :rtype: :class:`dict` with :class:`str` keys and :class:`ufl.form.Form` values + ``` + + :returns: the function for obtaining the form + :rtype: see docstring above + """ if self._get_form is None: raise NotImplementedError("'get_form' needs implementing.") return self._get_form(self) - def get_solver(self) -> Callable: + def get_solver(self): + """ + Get the function mapping a subinterval index and an initial condition dictionary + to a dictionary of solutions for the corresponding solver step. + + Signature for the function to be returned: + ``` + :arg index: the subinterval index + :type index: :class:`int` + :arg ic: map from fields to the corresponding initial condition components + :type ic: :class:`dict` with :class:`str` keys and + :class:`firedrake.function.Function` values + :return: map from fields to the corresponding solutions + :rtype: :class:`dict` with :class:`str` keys and + :class:`firedrake.function.Function` values + ``` + + :returns: the function for obtaining the solver + :rtype: see docstring above + """ if self._get_solver is None: raise NotImplementedError("'get_solver' needs implementing.") return self._get_solver(self) - def get_bcs(self) -> Callable: + def get_bcs(self): + """ + Get the function mapping a subinterval index to a set of Dirichlet boundary + conditions. + + Signature for the function to be returned: + ``` + :arg index: the subinterval index + :type index: :class:`int` + :return: boundary conditions + :rtype: :class:`~.DirichletBC` or :class:`list` thereof + :rtype: see docstring above + ``` + + :returns: the function for obtaining the boundary conditions + """ if self._get_bcs is not None: return self._get_bcs(self) - @property - def _function_spaces_consistent(self) -> bool: + def _transfer(self, source, target_space, **kwargs): + """ + Transfer a field between meshes using the specified transfer method. + + :arg source: the function to be transferred + :type source: :class:`firedrake.function.Function` or + :class:`firedrake.cofunction.Cofunction` + :arg target_space: the function space which we seek to transfer onto, or the + function or cofunction to use as the target + :type target_space: :class:`firedrake.functionspaceimpl.FunctionSpace`, + :class:`firedrake.function.Function` + or :class:`firedrake.cofunction.Cofunction` + :returns: the transferred function + :rtype: :class:`firedrake.function.Function` or + :class:`firedrake.cofunction.Cofunction` + + Extra keyword arguments are passed to :func:`goalie.interpolation.transfer`. + """ + return transfer(source, target_space, self._transfer_method, **kwargs) + + def _function_spaces_consistent(self): + """ + Determine whether the mesh sequence's function spaces are consistent with its + meshes. + + :returns: ``True`` if the meshes and function spaces are consistent, otherwise + ``False`` + :rtype: `:class:`bool` + """ consistent = len(self.time_partition) == len(self) consistent &= all(len(self) == len(self._fs[field]) for field in self.fields) for field in self.fields: @@ -232,65 +397,93 @@ def _function_spaces_consistent(self) -> bool: ) return consistent - def update_function_spaces(self) -> list: + def _update_function_spaces(self): """ - Update the function space dictionary associated with the :class:`MeshSeq`. + Update the function space dictionary associated with the mesh sequence. """ - if self._fs is None or not self._function_spaces_consistent: - self._fs = [self.get_function_spaces(mesh) for mesh in self.meshes] + if self._fs is None or not self._function_spaces_consistent(): self._fs = AttrDict( { - field: [self._fs[i][field] for i in range(len(self))] + field: [self.get_function_spaces(mesh)[field] for mesh in self] for field in self.fields } ) assert ( - self._function_spaces_consistent + self._function_spaces_consistent() ), "Meshes and function spaces are inconsistent" - return self._fs @property def function_spaces(self): - return self.update_function_spaces() + """ + Get the function spaces associated with the mesh sequence. + + :returns: a dictionary whose keys are field names and whose values are the + corresponding function spaces + :rtype: :class:`~.AttrDict` with :class:`str` keys and + :class:`firedrake.functionspaceimpl.FunctionSpace` values + """ + self._update_function_spaces() + return self._fs @property - def initial_condition(self) -> AttrDict: - ic = OrderedDict(self.get_initial_condition()) - assert issubclass( - ic.__class__, dict + def initial_condition(self): + """ + Get the initial conditions associated with the first subinterval. + + :returns: a dictionary whose keys are field names and whose values are the + corresponding initial conditions applied on the first subinterval + :rtype: :class:`~.AttrDict` with :class:`str` keys and + :class:`firedrake.function.Function` values + """ + initial_condition_map = self.get_initial_condition() + assert isinstance( + initial_condition_map, dict ), "`get_initial_condition` should return a dict" - assert set(self.fields).issubset( - set(ic.keys()) + mesh_seq_fields = set(self.fields) + initial_condition_fields = set(initial_condition_map.keys()) + assert mesh_seq_fields.issubset( + initial_condition_fields ), "missing fields in initial condition" - assert set(ic.keys()).issubset( - set(self.fields) + assert initial_condition_fields.issubset( + mesh_seq_fields ), "more initial conditions than fields" - return AttrDict(ic) + return AttrDict(initial_condition_map) @property - def form(self) -> Callable: + def form(self): + """ + See :meth:`~.MeshSeq.get_form`. + """ return self.get_form() @property - def solver(self) -> Callable: + def solver(self): + """ + See :meth:`~.MeshSeq.get_solver`. + """ return self.get_solver() @property - def bcs(self) -> Callable: + def bcs(self): + """ + See :meth:`~.MeshSeq.get_bcs`. + """ return self.get_bcs() @PETSc.Log.EventDecorator() - def get_checkpoints( - self, solver_kwargs: dict = {}, run_final_subinterval: bool = False - ) -> list: - """ + def get_checkpoints(self, solver_kwargs={}, run_final_subinterval=False): + r""" Solve forward on the sequence of meshes, extracting checkpoints corresponding to the starting fields on each subinterval. - :kwarg solver_kwargs: additional keyword arguments which will be passed to the - solver - :kwarg run_final_subinterval: toggle whether to solve the PDE on the final + :kwarg solver_kwargs: parameters for the forward solver + :type solver_kwargs: :class:`dict` with :class:`str` keys and values which may + take various types + :kwarg run_final_subinterval: if ``True``, the solver is run on the final subinterval + :type run_final_subinterval: :class:`bool` + :returns: checkpoints for each subinterval + :rtype: :class:`list` of :class:`firedrake.function.Function`\s """ N = len(self) @@ -304,14 +497,14 @@ def get_checkpoints( # Otherwise, solve each subsequent subinterval, in each case making use of the # previous checkpoint for i in range(N if run_final_subinterval else N - 1): - sols = self.solver(i, checkpoints[i], **solver_kwargs) - if not isinstance(sols, dict): + solutions = self.solver(i, checkpoints[i], **solver_kwargs) + if not isinstance(solutions, dict): raise TypeError( - f"Solver should return a dictionary, not '{type(sols)}'." + f"Solver should return a dictionary, not '{type(solutions)}'." ) # Check that the output of the solver is as expected - fields = set(sols.keys()) + fields = set(solutions.keys()) if not set(self.fields).issubset(fields): diff = set(self.fields).difference(fields) raise ValueError(f"Fields are missing from the solver: {diff}.") @@ -319,12 +512,12 @@ def get_checkpoints( diff = fields.difference(set(self.fields)) raise ValueError(f"Unexpected solver outputs: {diff}.") - # Transfer between meshes using conservative projection + # Transfer between meshes if i < N - 1: checkpoints.append( AttrDict( { - field: project(sols[field], fs[i + 1]) + field: self._transfer(solutions[field], fs[i + 1]) for field, fs in self._fs.items() } ) @@ -333,13 +526,17 @@ def get_checkpoints( return checkpoints @PETSc.Log.EventDecorator() - def get_solve_blocks(self, field: str, subinterval: int) -> list: - """ + def get_solve_blocks(self, field, subinterval): + r""" Get all blocks of the tape corresponding to solve steps for prognostic solution field on a given subinterval. :arg field: name of the prognostic solution field + :type field: :class:`str` :arg subinterval: subinterval index + :type subinterval: :class:`int` + :returns: list of solve blocks + :rtype: :class:`list` of :class:`pyadjoint.block.Block`\s """ tape = pyadjoint.get_working_tape() if tape is None: @@ -403,13 +600,19 @@ def get_solve_blocks(self, field: str, subinterval: int) -> list: def _output(self, field, subinterval, solve_block): """ - For a given solve block and solution field, get the block's outputs which - corresponds to the solution from the current timestep. + For a given solve block and solution field, get the block's outputs corresponding + to the solution from the current timestep. :arg field: field of interest + :type field: :class:`str` :arg subinterval: subinterval index - :arg solve_block: taped :class:`firedrake.adjoint.blocks.GenericSolveBlock` + :type subinterval: :class:`int` + :arg solve_block: taped solve block + :type solve_block: :class:`firedrake.adjoint.blocks.GenericSolveBlock` + :returns: the output + :rtype: :class:`firedrake.function.Function` """ + # TODO #93: Inconsistent return value - can be None fs = self.function_spaces[field][subinterval] # Loop through the solve block's outputs @@ -450,9 +653,15 @@ def _dependency(self, field, subinterval, solve_block): corresponds to the solution from the previous timestep. :arg field: field of interest + :type field: :class:`str` :arg subinterval: subinterval index - :arg solve_block: taped :class:`firedrake.adjoint.blocks.GenericSolveBlock` + :type subinterval: :class:`int` + :arg solve_block: taped solve block + :type solve_block: :class:`firedrake.adjoint.blocks.GenericSolveBlock` + :returns: the dependency + :rtype: :class:`firedrake.function.Function` """ + # TODO #93: Inconsistent return value - can be None if self.field_types[field] == "steady": return fs = self.function_spaces[field][subinterval] @@ -490,43 +699,42 @@ def _dependency(self, field, subinterval, solve_block): ) def _create_solutions(self): + """ + Create the :class:`~.FunctionData` instance for holding solution data. + """ self._solutions = ForwardSolutionData(self.time_partition, self.function_spaces) @property def solutions(self): """ - Arrays holding exported solution fields and their lagged counterparts. + :returns: the solution data object + :rtype: :class:`~.FunctionData` """ if not hasattr(self, "_solutions"): self._create_solutions() return self._solutions @PETSc.Log.EventDecorator() - def solve_forward(self, solver_kwargs: dict = {}) -> AttrDict: - """ + def solve_forward(self, solver_kwargs={}): + r""" Solve a forward problem on a sequence of subintervals. - A dictionary of solution fields is computed, the contents of which give values - at all exported timesteps, indexed first by the field label and then by type. - The contents of these nested dictionaries are nested lists which are indexed - first by subinterval and then by export. For a given exported timestep, the - field types are: - - * ``'forward'``: the forward solution after taking the timestep; - * ``'forward_old'``: the forward solution before taking the timestep (provided - the problem is not steady-state). + A dictionary of solution fields is computed - see :class:`~.ForwardSolutionData` + for more details. - :kwarg solver_kwargs: a dictionary providing parameters to the solver. Any - keyword arguments for the QoI should be included as a subdict with label - 'qoi_kwargs' - - :return solution: an :class:`~.AttrDict` containing solution fields and their - lagged versions. This can also be accessed as :meth:`solutions`. + :kwarg solver_kwargs: parameters for the forward solver + :type solver_kwargs: :class:`dict` whose keys are :class:`str`\s and whose values + may take various types + :returns: the solution data of the forward solves + :rtype: :class:`~.ForwardSolutionData` """ num_subintervals = len(self) P = self.time_partition solver = self.solver + # Reinitialise the solution data object + self._create_solutions() + # Start annotating if pyadjoint.annotate_tape(): tape = pyadjoint.get_working_tape() @@ -569,23 +777,23 @@ def solve_forward(self, solver_kwargs: dict = {}) -> AttrDict: ) # Update solution data based on block dependencies and outputs - sols = self.solutions[field] + solutions = self.solutions.extract(layout="field")[field] for j, block in enumerate(reversed(solve_blocks[::-stride])): # Current solution is determined from outputs out = self._output(field, i, block) if out is not None: - sols.forward[i][j].assign(out.saved_output) + solutions.forward[i][j].assign(out.saved_output) # Lagged solution comes from dependencies dep = self._dependency(field, i, block) if not self.steady and dep is not None: - sols.forward_old[i][j].assign(dep.saved_output) + solutions.forward_old[i][j].assign(dep.saved_output) # Transfer the checkpoint between subintervals if i < num_subintervals - 1: checkpoint = AttrDict( { - field: project(checkpoint[field], fs[i + 1]) + field: self._transfer(checkpoint[field], fs[i + 1]) for field, fs in self._fs.items() } ) @@ -596,12 +804,13 @@ def solve_forward(self, solver_kwargs: dict = {}) -> AttrDict: return self.solutions def check_element_count_convergence(self): - """ + r""" Check for convergence of the fixed point iteration due to the relative difference in element count being smaller than the specified tolerance. - :return: Boolean array with ``True`` in the appropriate entry if convergence is - detected on a subinterval. + :return: an array, whose entries are ``True`` if convergence is detected on the + corresponding subinterval + :rtype: :class:`list` of :class:`bool`\s """ if self.params.drop_out_converged: converged = self.converged @@ -639,27 +848,28 @@ def check_element_count_convergence(self): @PETSc.Log.EventDecorator() def fixed_point_iteration( - self, - adaptor: Callable, - solver_kwargs: dict = {}, - adaptor_kwargs: dict = {}, - **kwargs, + self, adaptor, update_params=None, solver_kwargs={}, adaptor_kwargs={} ): r""" - Apply goal-oriented mesh adaptation using a fixed point iteration loop. + Apply mesh adaptation using a fixed point iteration loop approach. - :arg adaptor: function for adapting the mesh sequence. Its arguments are the - :class:`~.MeshSeq` instance and the dictionary of solution - :class:`firedrake.function.Function`\s. It should return ``True`` if the + :arg adaptor: function for adapting the mesh sequence. Its arguments are the mesh + sequence and the solution data object. It should return ``True`` if the convergence criteria checks are to be skipped for this iteration. Otherwise, it should return ``False``. :kwarg update_params: function for updating :attr:`~.MeshSeq.params` at each iteration. Its arguments are the parameter class and the fixed point iteration - :kwarg solver_kwargs: a dictionary providing parameters to the solver - :kwarg adaptor_kwargs: a dictionary providing parameters to the adaptor + :kwarg solver_kwargs: parameters to pass to the solver + :type solver_kwargs: :class:`dict` with :class:`str` keys and values which may + take various types + :kwarg adaptor_kwargs: parameters to pass to the adaptor + :type adaptor_kwargs: :class:`dict` with :class:`str` keys and values which may + take various types + :returns: solution data object + :rtype: :class:`~.ForwardSolutionData` """ - update_params = kwargs.get("update_params") + # TODO #124: adaptor no longer needs solution data to be passed explicitly self.element_counts = [self.count_elements()] self.vertex_counts = [self.count_vertices()] self.converged[:] = False @@ -670,7 +880,6 @@ def fixed_point_iteration( update_params(self.params, self.fp_iteration) # Solve the forward problem over all meshes - self._create_solutions() self.solve_forward(solver_kwargs=solver_kwargs) # Adapt meshes, logging element and vertex counts diff --git a/goalie/metric.py b/goalie/metric.py index 2d6c12e6..167e1b76 100644 --- a/goalie/metric.py +++ b/goalie/metric.py @@ -1,54 +1,46 @@ """ Driver functions for metric-based mesh adaptation. """ -from animate.interpolation import clement_interpolant -import animate.metric -import animate.recovery -from .log import debug -from .time_partition import TimePartition -from animate.metric import RiemannianMetric + +from collections.abc import Iterable + import firedrake -from firedrake.petsc import PETSc import numpy as np -from pyop2 import op2 -from typing import List, Optional, Union import ufl +from animate.metric import RiemannianMetric +from firedrake.petsc import PETSc +from .log import debug -__all__ = ["enforce_element_constraints", "space_time_normalise", "ramp_complexity"] +__all__ = ["enforce_variable_constraints", "space_time_normalise", "ramp_complexity"] -# TODO: Implement this on the PETSc level and port it through to Firedrake @PETSc.Log.EventDecorator() -def enforce_element_constraints( - metrics: List[RiemannianMetric], - h_min: List[firedrake.Function], - h_max: List[firedrake.Function], - a_max: List[firedrake.Function], - boundary_tag: Optional[Union[str, int]] = None, - optimise: bool = False, -) -> List[firedrake.Function]: - """ - Post-process a list of metrics to enforce minimum and - maximum element sizes, as well as maximum anisotropy. +def enforce_variable_constraints( + metrics, + h_min=1.0e-30, + h_max=1.0e30, + a_max=1.0e5, + boundary_tag=None, +): + r""" + Post-process a list of metrics to enforce minimum and maximum element sizes, as well + as maximum anisotropy. :arg metrics: the metrics - :arg h_min: minimum tolerated element size, - which could be a :class:`firedrake.function.Function` - or a number. - :arg h_max: maximum tolerated element size, - which could be a :class:`firedrake.function.Function` - or a number. - :arg a_max: maximum tolerated element anisotropy, - which could be a :class:`firedrake.function.Function` - or a number. + :type metrics: :class:`list` of :class:`~.RiemannianMetric`\s + :kwarg h_min: minimum tolerated element size + :type h_min: :class:`firedrake.function.Function`, :class:`float`, or :class:`int` + :kwarg h_max: maximum tolerated element size + :type h_max: :class:`firedrake.function.Function`, :class:`float`, or :class:`int` + :kwarg a_max: maximum tolerated element anisotropy + :type a_max: :class:`firedrake.function.Function`, :class:`float`, or :class:`int` :kwarg boundary_tag: optional tag to enforce sizes on. - :kwarg optimise: is this a timed run? + :type boundary_tag: :class:`str` or :class:`int` """ - from collections.abc import Iterable - if isinstance(metrics, RiemannianMetric): metrics = [metrics] + assert isinstance(metrics, Iterable) if not isinstance(h_min, Iterable): h_min = [h_min] * len(metrics) if not isinstance(h_max, Iterable): @@ -56,85 +48,53 @@ def enforce_element_constraints( if not isinstance(a_max, Iterable): a_max = [a_max] * len(metrics) for metric, hmin, hmax, amax in zip(metrics, h_min, h_max, a_max): - fs = metric.function_space() - mesh = fs.mesh() - P1 = firedrake.FunctionSpace(mesh, "CG", 1) - - def interp(f): - if isinstance(f, firedrake.Function): - return clement_interpolant(f) - else: - return firedrake.Function(P1).assign(f) - - # Interpolate hmin, hmax and amax into P1 - hmin = interp(hmin) - hmax = interp(hmax) - amax = interp(amax) - - # Check the values are okay - if not optimise: - _hmin = hmin.vector().gather().min() - if _hmin <= 0.0: - raise ValueError( - f"Encountered negative non-positive hmin value: {_hmin}." - ) - _hmax = hmax.vector().gather().min() - if _hmax < _hmin: - raise ValueError( - f"Encountered hmax value smaller than hmin: {_hmax} vs. {_hmin}." - ) - dx = ufl.dx(domain=mesh) - integral = firedrake.assemble(ufl.conditional(hmax < hmin, 1, 0) * dx) - if not np.isclose(integral, 0.0): - raise ValueError( - f"Encountered regions where hmax < hmin: volume {integral}." - ) - _amax = amax.vector().gather().min() - if _amax < 1.0: - raise ValueError(f"Encountered amax value smaller than unity: {_amax}.") - - # Enforce constraints - dim = fs.mesh().topological_dimension() - if boundary_tag is None: - node_set = fs.node_set - else: - node_set = firedrake.DirichletBC(fs, 0, boundary_tag).node_set - op2.par_loop( - animate.recovery.get_metric_kernel("postproc_metric", dim), - node_set, - metric.dat(op2.RW), - hmin.dat(op2.READ), - hmax.dat(op2.READ), - amax.dat(op2.READ), + metric.set_parameters( + { + "dm_plex_metric_h_min": hmin, + "dm_plex_metric_h_max": hmax, + "dm_plex_metric_a_max": amax, + "dm_plex_metric_boundary_tag": boundary_tag, + } ) + metric.enforce_spd() return metrics @PETSc.Log.EventDecorator() def space_time_normalise( - metrics: List[RiemannianMetric], - time_partition: TimePartition, - metric_parameters: Union[dict, list], - global_factor: Optional[float] = None, - boundary: bool = False, - restrict_sizes: bool = True, - restrict_anisotropy: bool = True, -) -> List[RiemannianMetric]: + metrics, + time_partition, + metric_parameters, + global_factor=None, + boundary=False, + restrict_sizes=True, + restrict_anisotropy=True, +): r""" Apply :math:`L^p` normalisation in both space and time. Based on Equation (1) in :cite:`Barral:2016`. - :arg metrics: list of :class:`RiemannianMetric`\s corresponding - to the metric associated with each subinterval - :arg time_partition: :class:`TimePartition` for the problem at hand + :arg metrics: the metrics associated with each subinterval + :type metrics: :class:`list` of :class:`~.RiemannianMetric`\s + :arg time_partition: temporal discretisation for the problem at hand + :type time_partition: :class:`TimePartition` :arg metric_parameters: dictionary containing the target *space-time* metric complexity under `dm_plex_metric_target_complexity` and the normalisation order under `dm_plex_metric_p`, or a list thereof + :type metric_parameters: :class:`list` of :class:`dict`\s or a single :class:`dict` + to use for all subintervals :kwarg global_factor: pre-computed global normalisation factor - :kwarg boundary: is the normalisation to be done over the boundary? - :kwarg restrict_sizes: should minimum and maximum metric magnitudes be enforced? - :kwarg restrict_anisotropy: should maximum anisotropy be enforced? + :type global_factor: :class:`float` + :kwarg boundary: if ``True``, the normalisation to be performed over the boundary + :type boundary: :class:`bool` + :kwarg restrict_sizes: if ``True``, minimum and maximum metric magnitudes are + enforced + :type restrict_sizes: :class:`bool` + :kwarg restrict_anisotropy: if ``True``, maximum anisotropy is enforced + :type restrict_anisotropy: :class:`bool` + :returns: the space-time normalised metrics + :rtype: :class:`list` of :class:`~.RiemannianMetric`\s """ if isinstance(metric_parameters, dict): metric_parameters = [metric_parameters for _ in range(len(time_partition))] @@ -217,26 +177,30 @@ def space_time_normalise( return metrics -def ramp_complexity( - base: float, target: float, i: int, num_iterations: int = 3 -) -> float: +def ramp_complexity(base, target, iteration, num_iterations=3): """ Ramp up the target complexity over the first few iterations. :arg base: the base complexity to start from + :type base: :class:`float` :arg target: the desired complexity - :arg i: the current iteration + :type target: :class:`float` + :arg iteration: the current iteration + :type iteration: :class:`int` :kwarg num_iterations: how many iterations to ramp over? + :type num_iterations: :class:`int` + :returns: the ramped target complexity + :rtype: :class:`float` """ if base <= 0.0: raise ValueError(f"Base complexity must be positive, not {base}.") if target <= 0.0: raise ValueError(f"Target complexity must be positive, not {target}.") - if i < 0: - raise ValueError(f"Current iteration must be non-negative, not {i}.") + if iteration < 0: + raise ValueError(f"Current iteration must be non-negative, not {iteration}.") if num_iterations < 0: raise ValueError( f"Number of iterations must be non-negative, not {num_iterations}." ) - alpha = 1 if num_iterations == 0 else min(i / num_iterations, 1) + alpha = 1 if num_iterations == 0 else min(iteration / num_iterations, 1) return alpha * target + (1 - alpha) * base diff --git a/goalie/options.py b/goalie/options.py index 1450d546..b11f9fb1 100644 --- a/goalie/options.py +++ b/goalie/options.py @@ -1,6 +1,6 @@ -from .utility import AttrDict from animate.adapt import RiemannianMetric +from .utility import AttrDict __all__ = [ "AdaptParameters", @@ -16,9 +16,11 @@ class AdaptParameters(AttrDict): loops. """ - def __init__(self, parameters: dict = {}): + def __init__(self, parameters={}): """ - :arg parameters: dictionary of parameters to set + :arg parameters: parameters to set + :type parameters: :class:`dict` with :class:`str` keys and values which may take + various types """ self["miniter"] = 3 # Minimum iteration count self["maxiter"] = 35 # Maximum iteration count @@ -42,6 +44,14 @@ def __init__(self, parameters: dict = {}): self._check_type("drop_out_converged", bool) def _check_type(self, key, expected): + """ + Check that a given parameter is of the expected type. + + :arg key: the parameter label + :type key: :class:`str` + :arg expected: the expected type + :type expected: :class:`type` + """ if not isinstance(self[key], expected): if isinstance(expected, tuple): name = "' or '".join([e.__name__ for e in expected]) @@ -53,18 +63,26 @@ def _check_type(self, key, expected): ) def _check_value(self, key, possibilities): + """ + Check that a given parameter takes one of the possible values. + + :arg key: the parameter label + :type key: :class:`str` + :arg possibilities: all possible values for the parameter + :type possibilities: :class:`list` + """ value = self[key] if value not in possibilities: raise ValueError( f"Unsupported value '{value}' for '{key}'. Choose from {possibilities}." ) - def __str__(self) -> str: + def __str__(self): return str({key: value for key, value in self.items()}) - def __repr__(self) -> str: + def __repr__(self): d = ", ".join([f"{key}={value}" for key, value in self.items()]) - return f"{self.__class__.__name__}({d})" + return f"{type(self).__name__}({d})" class MetricParameters(AdaptParameters): @@ -73,9 +91,11 @@ class MetricParameters(AdaptParameters): metric-based adaptive mesh fixed point iteration loops. """ - def __init__(self, parameters: dict = {}): + def __init__(self, parameters={}): """ - :arg parameters: dictionary of parameters to set + :arg parameters: parameters to set + :type parameters: :class:`dict` with :class:`str` keys and values which may take + various types """ self["num_ramp_iterations"] = 3 # Number of iterations to ramp over self["verbosity"] = -1 # -1 = silent, 10 = maximum @@ -124,7 +144,8 @@ def export(self, metric): """ Set parameters appropriate to a given :class:`RiemannianMetric`. - :arg metric: the :class:`RiemannianMetric` to apply parameters to + :arg metric: the metric to apply parameters to + :type metric: :class:`animate.metric.RiemannianMetric` """ if not isinstance(metric, RiemannianMetric): raise TypeError( @@ -146,11 +167,9 @@ def export(self, metric): "hausdorff_number", "gradation_factor", ) - metric_parameters = { - "dm_plex_metric": {key: self[key] for key in petsc_specific} - } - metric_parameters["num_iterations"] = self["num_parmmg_iterations"] - metric.set_parameters(metric_parameters) + metric_parameters_sub = {"num_iterations": self["num_parmmg_iterations"]} + metric_parameters_sub.update({key: self[key] for key in petsc_specific}) + metric.set_parameters({"dm_plex_metric": metric_parameters_sub}) class GoalOrientedParameters(AdaptParameters): @@ -160,7 +179,12 @@ class GoalOrientedParameters(AdaptParameters): loops. """ - def __init__(self, parameters: dict = {}): + def __init__(self, parameters={}): + """ + :arg parameters: parameters to set + :type parameters: :class:`dict` with :class:`str` keys and values which may take + various types + """ self["qoi_rtol"] = 0.001 # Relative tolerance for QoI self["estimator_rtol"] = 0.001 # Relative tolerance for estimator self["convergence_criteria"] = "any" # Mode for convergence checking diff --git a/goalie/point_seq.py b/goalie/point_seq.py index 389e2f43..f63d3028 100644 --- a/goalie/point_seq.py +++ b/goalie/point_seq.py @@ -1,5 +1,6 @@ import firedrake import firedrake.mesh as fmesh + from .mesh_seq import MeshSeq __all__ = ["PointSeq"] @@ -43,4 +44,4 @@ def set_meshes(self, mesh): self.meshes = [mesh for _ in self.subintervals] self.dim = mesh.topological_dimension() assert self.dim == 0 - self.reset_counts() + self._reset_counts() diff --git a/goalie/time_partition.py b/goalie/time_partition.py index 1a59aef0..f56df865 100644 --- a/goalie/time_partition.py +++ b/goalie/time_partition.py @@ -1,11 +1,12 @@ """ Partitioning for the temporal domain. """ -from .log import debug + from collections.abc import Iterable + import numpy as np -from typing import List, Optional, Union +from .log import debug __all__ = ["TimePartition", "TimeInterval", "TimeInstant"] @@ -20,27 +21,37 @@ class TimePartition: def __init__( self, - end_time: float, - num_subintervals: int, - timesteps: Union[List[float], float], - fields: Union[List[str], str], - num_timesteps_per_export: int = 1, - start_time: float = 0.0, - subintervals: Optional[List[float]] = None, - field_types: Optional[Union[List[str], str]] = None, + end_time, + num_subintervals, + timesteps, + fields, + num_timesteps_per_export=1, + start_time=0.0, + subintervals=None, + field_types=None, ): - """ + r""" :arg end_time: end time of the interval of interest + :type end_time: :class:`float` or :class:`int` :arg num_subintervals: number of subintervals in the partition - :arg timesteps: (list of values for the) timestep used on each subinterval - :arg fields: (list of) field names - :kwarg num_timesteps_per_export: (list of) number of timesteps per export + :type num_subintervals: :class:`int` + :arg timesteps: a list timesteps to be used on each subinterval, or a single + timestep to use for all subintervals + :type timesteps: :class:`list` of :class:`float`\s or :class:`float` + :arg fields: the list of field names to consider + :type fields: :class:`list` of :class:`str`\s or :class:`str` + :kwarg num_timesteps_per_export: a list of numbers of timesteps per export for + each subinterval, or a single number to use for all subintervals + :type num_timesteps_per_export: :class:`list` of :class`int`\s or :class:`int` :kwarg start_time: start time of the interval of interest - :kwarg subinterals: user-provided sequence of subintervals, which need not be of - uniform length - :kwarg field_types: (list of) strings indicating whether each field is + :type start_time: :class:`float` or :class:`int` + :kwarg subinterals: sequence of subintervals (which need not be of uniform + length), or ``None`` to use uniform subintervals (the default) + :type subintervals: :class:`list` of :class:`tuple`\s + :kwarg field_types: a list of strings indicating whether each field is 'unsteady' or 'steady', i.e., does the corresponding equation involve time derivatives or not? + :type field_types: :class:`list` of :class:`str`\s or :class:`str` """ debug(100 * "-") if isinstance(fields, str): @@ -120,9 +131,11 @@ def __init__( debug("field_types") debug(100 * "-") - def debug(self, attr: str): + def debug(self, attr): """ Print attribute 'msg' for debugging purposes. + + :arg attr: the attribute to display debugging information for """ try: val = self.__getattribute__(attr) @@ -133,10 +146,10 @@ def debug(self, attr: str): label = " ".join(attr.split("_")) debug(f"TimePartition: {label:25s} {val}") - def __str__(self) -> str: + def __str__(self): return f"{self.subintervals}" - def __repr__(self) -> str: + def __repr__(self): timesteps = ", ".join([str(dt) for dt in self.timesteps]) fields = ", ".join([f"'{field}'" for field in self.fields]) return ( @@ -147,14 +160,17 @@ def __repr__(self) -> str: f"fields=[{fields}])" ) - def __len__(self) -> int: + def __len__(self): return self.num_subintervals - def __getitem__(self, sl: Union[int, slice]) -> dict: + def __getitem__(self, index_or_slice): """ - :arg sl: index or slice - :return: subinterval bounds and timestep associated with that index + :arg index_or_slice: an index or slice to generate a sub-time partition for + :type index_or_slice: :class:`int` or :class:`slice` + :returns: a time partition for the given index or slice + :rtype: :class:`~.TimePartition` """ + sl = index_or_slice if not isinstance(sl, slice): sl = slice(sl, sl + 1, 1) step = sl.step or 1 @@ -175,6 +191,10 @@ def __getitem__(self, sl: Union[int, slice]) -> dict: @property def num_timesteps(self): + """ + :returns the total number of timesteps + :rtype: :class:`int` + """ return sum(self.num_timesteps_per_subinterval) def _check_subintervals(self): @@ -274,7 +294,7 @@ def __ne__(self, other): class TimeInterval(TimePartition): """ - A trivial :class:`~.TimePartition`. + A trivial :class:`~.TimePartition` with a single subinterval. """ def __init__(self, *args, **kwargs): @@ -288,7 +308,7 @@ def __init__(self, *args, **kwargs): fields = args[2] super().__init__(end_time, 1, timestep, fields, **kwargs) - def __repr__(self) -> str: + def __repr__(self): return ( f"TimeInterval(" f"end_time={self.end_time}, " @@ -297,7 +317,11 @@ def __repr__(self) -> str: ) @property - def timestep(self) -> float: + def timestep(self): + """ + :returns: the timestep used on the single interval + :rtype: :class:`float` + """ return self.timesteps[0] @@ -305,11 +329,10 @@ class TimeInstant(TimeInterval): """ A :class:`~.TimePartition` for steady-state problems. - Under the hood this means dividing :math:`[0,1)` into - a single timestep. + Under the hood this means dividing :math:`[0,1)` into a single timestep. """ - def __init__(self, fields: Union[List[str], str], **kwargs): + def __init__(self, fields, **kwargs): if "end_time" in kwargs: if "time" in kwargs: raise ValueError("Both 'time' and 'end_time' are set.") @@ -319,8 +342,8 @@ def __init__(self, fields: Union[List[str], str], **kwargs): timestep = time super().__init__(time, timestep, fields, **kwargs) - def __str__(self) -> str: + def __str__(self): return f"({self.end_time})" - def __repr__(self) -> str: + def __repr__(self): return f"TimeInstant(" f"time={self.end_time}, " f"fields={self.fields})" diff --git a/goalie/utility.py b/goalie/utility.py index 089c8969..4971f8a6 100644 --- a/goalie/utility.py +++ b/goalie/utility.py @@ -2,257 +2,15 @@ Utility functions and classes for mesh adaptation. """ -from collections import OrderedDict -import firedrake -import firedrake.mesh as fmesh -from firedrake.petsc import PETSc -from firedrake.__future__ import interpolate -import mpi4py -import numpy as np import os -import petsc4py -import ufl - - -@PETSc.Log.EventDecorator("goalie.Mesh") -def Mesh(arg, **kwargs) -> firedrake.mesh.MeshGeometry: - """ - Overload :func:`firedrake.mesh.Mesh` to - endow the output mesh with useful quantities. - - The following quantities are computed by default: - * cell size; - * facet area. - - The argument and keyword arguments are passed to - Firedrake's ``Mesh`` constructor, modified so - that the argument could also be a mesh. - """ - try: - mesh = firedrake.Mesh(arg, **kwargs) - except TypeError: - mesh = firedrake.Mesh(arg.coordinates, **kwargs) - if isinstance(mesh._topology, fmesh.VertexOnlyMeshTopology): - return mesh - P0 = firedrake.FunctionSpace(mesh, "DG", 0) - P1 = firedrake.FunctionSpace(mesh, "CG", 1) - dim = mesh.topological_dimension() - - # Facet area - boundary_markers = sorted(mesh.exterior_facets.unique_markers) - one = firedrake.Function(P1).assign(1.0) - bnd_len = OrderedDict( - {i: firedrake.assemble(one * ufl.ds(int(i))) for i in boundary_markers} - ) - if dim == 2: - mesh.boundary_len = bnd_len - else: - mesh.boundary_area = bnd_len - - # Cell size - if dim == 2 and mesh.coordinates.ufl_element().cell == ufl.triangle: - mesh.delta_x = firedrake.assemble(interpolate(ufl.CellDiameter(mesh), P0)) - - return mesh - - -class VTKFile(firedrake.output.VTKFile): - """ - Overload :class:`firedrake.output.VTKFile` so that it uses ``adaptive`` mode by - default. Whilst this means that the mesh topology is recomputed at every export, it - removes any need for the user to reset it manually. - """ - - def __init__(self, *args, **kwargs): - kwargs.setdefault("adaptive", True) - super().__init__(*args, **kwargs) - - def _write_vtu(self, *functions): - """ - Overload the Firedrake functionality - under the blind assumption that the - same list of functions are outputted - each time (albeit on different meshes). - """ - if self._fnames is not None: - if len(self._fnames) != len(functions): - raise ValueError( - "Writing different number of functions: expected" - f" {len(self._fnames)}, got {len(functions)}." - ) - for name, f in zip(self._fnames, functions): - if f.name() != name: - f.rename(name) - return super()._write_vtu(*functions) - - -@PETSc.Log.EventDecorator("goalie.assemble_mass_matrix") -def assemble_mass_matrix( - space: firedrake.FunctionSpace, norm_type: str = "L2" -) -> petsc4py.PETSc.Mat: - """ - Assemble the ``norm_type`` mass matrix - associated with some finite element ``space``. - """ - trial = firedrake.TrialFunction(space) - test = firedrake.TestFunction(space) - if norm_type == "L2": - lhs = ufl.inner(trial, test) * ufl.dx - elif norm_type == "H1": - lhs = ( - ufl.inner(trial, test) * ufl.dx - + ufl.inner(ufl.grad(trial), ufl.grad(test)) * ufl.dx - ) - else: - raise ValueError(f"Norm type '{norm_type}' not recognised.") - return firedrake.assemble(lhs).petscmat - - -def cofunction2function(c): - """ - Map a :class:`Cofunction` to a :class:`Function`. - """ - f = firedrake.Function(c.function_space().dual()) - if isinstance(f.dat.data_with_halos, tuple): - for i, arr in enumerate(f.dat.data_with_halos): - arr[:] = c.dat.data_with_halos[i] - else: - f.dat.data_with_halos[:] = c.dat.data_with_halos - return f - - -def function2cofunction(f): - """ - Map a :class:`Function` to a :class:`Cofunction`. - """ - c = firedrake.Cofunction(f.function_space().dual()) - if isinstance(c.dat.data_with_halos, tuple): - for i, arr in enumerate(c.dat.data_with_halos): - arr[:] = f.dat.data_with_halos[i] - else: - c.dat.data_with_halos[:] = f.dat.data_with_halos - return c - - -@PETSc.Log.EventDecorator() -def norm(v, norm_type="L2", **kwargs): - r""" - Overload :func:`firedrake.norms.norm` to allow for :math:`\ell^p` norms. - - Note that this version is case sensitive, i.e. ``'l2'`` and ``'L2'`` will give - different results in general. - - :arg v: the :class:`firedrake.function.Function` or - :class:`firedrake.cofunction.Cofunction` to take the norm of - :kwarg norm_type: choose from ``'l1'``, ``'l2'``, ``'linf'``, ``'L2'``, ``'Linf'``, - ``'H1'``, ``'Hdiv'``, ``'Hcurl'``, or any ``'Lp'`` with :math:`p >= 1`. - :kwarg condition: a UFL condition for specifying a subdomain to compute the norm - over - :kwarg boundary: should the norm be computed over the domain boundary? - """ - if isinstance(v, firedrake.Cofunction): - v = cofunction2function(v) - boundary = kwargs.get("boundary", False) - condition = kwargs.get("condition", firedrake.Constant(1.0)) - norm_codes = {"l1": 0, "l2": 2, "linf": 3} - p = 2 - if norm_type in norm_codes or norm_type == "Linf": - if boundary: - raise NotImplementedError("lp errors on the boundary not yet implemented.") - v.interpolate(condition * v) - if norm_type == "Linf": - with v.dat.vec_ro as vv: - return vv.max()[1] - else: - with v.dat.vec_ro as vv: - return vv.norm(norm_codes[norm_type]) - elif norm_type[0] == "l": - raise NotImplementedError( - "lp norm of order {:s} not supported.".format(norm_type[1:]) - ) - else: - dX = ufl.ds if boundary else ufl.dx - if norm_type.startswith("L"): - try: - p = int(norm_type[1:]) - except Exception: - raise ValueError(f"Don't know how to interpret '{norm_type}' norm.") - if p < 1: - raise ValueError(f"'{norm_type}' norm does not make sense.") - integrand = ufl.inner(v, v) - elif norm_type.lower() == "h1": - integrand = ufl.inner(v, v) + ufl.inner(ufl.grad(v), ufl.grad(v)) - elif norm_type.lower() == "hdiv": - integrand = ufl.inner(v, v) + ufl.div(v) * ufl.div(v) - elif norm_type.lower() == "hcurl": - integrand = ufl.inner(v, v) + ufl.inner(ufl.curl(v), ufl.curl(v)) - else: - raise ValueError(f"Unknown norm type '{norm_type}'.") - return firedrake.assemble(condition * integrand ** (p / 2) * dX) ** (1 / p) - -@PETSc.Log.EventDecorator() -def errornorm(u, uh, norm_type="L2", **kwargs): - r""" - Overload :func:`firedrake.norms.errornorm` to allow for :math:`\ell^p` norms. - - Note that this version is case sensitive, i.e. ``'l2'`` and ``'L2'`` will give - different results in general. - - :arg u: the 'true' value - :arg uh: the approximation of the 'truth' - :kwarg norm_type: choose from ``'l1'``, ``'l2'``, ``'linf'``, ``'L2'``, ``'Linf'``, - ``'H1'``, ``'Hdiv'``, ``'Hcurl'``, or any ``'Lp'`` with :math:`p >= 1`. - :kwarg boundary: should the norm be computed over the domain boundary? - """ - if isinstance(u, firedrake.Cofunction): - u = cofunction2function(u) - if isinstance(uh, firedrake.Cofunction): - uh = cofunction2function(uh) - if not isinstance(uh, firedrake.Function): - raise TypeError(f"uh should be a Function, is a '{type(uh)}'.") - if norm_type[0] == "l": - if not isinstance(u, firedrake.Function): - raise TypeError(f"u should be a Function, is a '{type(u)}'.") - - if len(u.ufl_shape) != len(uh.ufl_shape): - raise RuntimeError("Mismatching rank between u and uh.") - - if isinstance(u, firedrake.Function): - degree_u = u.function_space().ufl_element().degree() - degree_uh = uh.function_space().ufl_element().degree() - if degree_uh > degree_u: - firedrake.logging.warning( - "Degree of exact solution less than approximation degree" - ) - - # Case 1: point-wise norms - if norm_type[0] == "l": - v = u - v -= uh - - # Case 2: UFL norms for mixed function spaces - elif hasattr(uh.function_space(), "num_sub_spaces"): - if norm_type == "L2": - vv = [uu - uuh for uu, uuh in zip(u.subfunctions, uh.subfunctions)] - dX = ufl.ds if kwargs.get("boundary", False) else ufl.dx - return ufl.sqrt(firedrake.assemble(sum([ufl.inner(v, v) for v in vv]) * dX)) - else: - raise NotImplementedError( - f"Norm type '{norm_type}' not supported for mixed spaces." - ) - - # Case 3: UFL norms for non-mixed spaces - else: - v = u - uh - - return norm(v, norm_type=norm_type, **kwargs) +import firedrake +import numpy as np class AttrDict(dict): """ - Dictionary that provides both ``self[key]`` - and ``self.key`` access to members. + Dictionary that provides both ``self[key]`` and ``self.key`` access to members. **Disclaimer**: Copied from `stackoverflow `__. @@ -263,19 +21,20 @@ def __init__(self, *args, **kwargs): self.__dict__ = self -def effectivity_index(error_indicator: firedrake.Function, Je: float) -> float: +def effectivity_index(error_indicator, Je): r""" - Overestimation factor of some error estimator - for the QoI error. + Compute the overestimation factor of some error estimator for the QoI error. - Note that this is only typically used for simple - steady-state problems with analytical solutions. + Note that this is only typically used for simple steady-state problems with + analytical solutions. - :arg error_indicator: a :math:`\mathbb P0` - :class:`firedrake.function.Function` which - localises contributions to an error estimator - to individual elements - :arg Je: error in quantity of interest + :arg error_indicator: a :math:`\mathbb P0` error indicator which localises + contributions to an error estimator to individual elements + :type error_indicator: :class:`firedrake.function.Function` + :arg Je: the error in the quantity of interest + :type Je: :class:`float` + :returns: the effectivity index + :rtype: :class:`float` """ if not isinstance(error_indicator, firedrake.Function): raise ValueError("Error indicator must return a Function.") @@ -286,17 +45,18 @@ def effectivity_index(error_indicator: firedrake.Function, Je: float) -> float: return np.abs(eta / Je) -def create_directory( - path: str, comm: mpi4py.MPI.Intracomm = firedrake.COMM_WORLD -) -> str: +def create_directory(path, comm=firedrake.COMM_WORLD): """ Create a directory on disk. - Code copied from `Thetis - `__. + **Disclaimer**: Code copied from `Thetis `__. :arg path: path to the directory + :type path: :class:`str` :kwarg comm: MPI communicator + :type comm: :class:`mpi4py.MPI.Intracomm` + :returns: the path in absolute form + :rtype path: :class:`str` """ if comm.rank == 0: if not os.path.exists(path): diff --git a/goalie_adjoint/__init__.py b/goalie_adjoint/__init__.py index ac3567a2..bae10d55 100644 --- a/goalie_adjoint/__init__.py +++ b/goalie_adjoint/__init__.py @@ -1,5 +1,7 @@ from __future__ import absolute_import + +from pyadjoint import no_annotations # noqa + from goalie import * # noqa from goalie.adjoint import * # noqa from goalie.go_mesh_seq import * # noqa -from pyadjoint import no_annotations # noqa diff --git a/pyproject.toml b/pyproject.toml index 620169a9..94d4ab62 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,11 +1,62 @@ -[tool.black] +[build-system] +requires = ["setuptools"] + +[project] +name = "goalie" +version = "0.2" +dependencies = [ + "coverage", + "parameterized", + "pre-commit", + "ruff", +] +authors = [ + {name = "Joseph G. Wallwork", email = "joe.wallwork@outlook.com"}, + {name = "Davor Dundovic"}, + {name = "Eleda Johnson"}, + {name = "Stephan C. Kramer"}, +] +maintainers = [ + {name = "Joseph G. Wallwork", email = "joe.wallwork@outlook.com"}, + {name = "Davor Dundovic"}, + {name = "Eleda Johnson"}, + {name = "Stephan C. Kramer"}, +] +description = "Goal-oriented error estimation and mesh adaptation for finite element problems solved using Firedrake" +readme = "README.md" +license = {file = "LICENSE"} +classifiers = [ + "Development Status :: 4 - Beta", + "Programming Language :: Python", +] + +[project.urls] +Homepage = "https://mesh-adaptation.github.io" +Documentation = "https://mesh-adaptation.github.io/goalie/index.html" +Repository = "https://github.com/mesh-adaptation/goalie" + +[tool.setuptools] +packages = ["goalie", "goalie_adjoint"] + +[tool.ruff] line-length = 88 -include = '\.pyi?$' -exclude = ''' -/( - \.git -)/ -''' + +[tool.ruff.lint] +select = [ +# "B", # flake8-bugbear TODO: enable this (#165) +# "C", # mccabe complexity TODO: enable this (#165) + "E", "W", # Pycodestyle + "F", # Pyflakes + "I", # isort +] +ignore = [ + "E501", # line too long + "E226", # missing whitespace around arithmetic operator + "E402", # module level import not at top of file + "E741", # ambiguous variable name + "F403", # unable to detect undefined names + "F405", # name may be undefined, or defined from star imports +] [tool.pytest.ini_options] filterwarnings = [ diff --git a/requirements.txt b/requirements.txt deleted file mode 100644 index 2a14c02d..00000000 --- a/requirements.txt +++ /dev/null @@ -1,5 +0,0 @@ -black -coverage -flake8 -parameterized -pre-commit diff --git a/setup.py b/setup.py deleted file mode 100644 index f5fe27a8..00000000 --- a/setup.py +++ /dev/null @@ -1,12 +0,0 @@ -from setuptools import setup - - -setup( - name="goalie", - description="Goalie: Goal-Oriented Mesh Adaptation Toolkit", - author="Joseph Wallwork", - author_email="joe.wallwork@outlook.com", - url="https://github.com/pyroteus/goalie", - version="0.1", - packages=["goalie", "goalie_adjoint"], -) diff --git a/test/conftest.py b/test/conftest.py index 4a4820bd..b4293220 100644 --- a/test/conftest.py +++ b/test/conftest.py @@ -3,6 +3,7 @@ **Disclaimer: some functions copied from firedrake/src/tests/conftest.py """ + from subprocess import check_call @@ -81,9 +82,10 @@ def pytest_runtest_setup(item): if MPI.COMM_WORLD.size > 1: # Turn on source hash checking - from firedrake import parameters from functools import partial + from firedrake import parameters + def _reset(check): parameters["pyop2_options"]["check_src_hashes"] = check diff --git a/test/sensors.py b/test/sensors.py index 4de5f3d3..2d8330ce 100644 --- a/test/sensors.py +++ b/test/sensors.py @@ -5,11 +5,11 @@ adaptation for unsteady CFD simulations involving moving geometries. Diss. 2011. """ + import firedrake from ufl import * from utility import uniform_mesh - __all__ = ["bowl", "hyperbolic", "multiscale", "interweaved", "mesh_for_sensors"] diff --git a/test/test_error_estimation.py b/test/test_error_estimation.py index 3b54cf8a..bc74fe61 100644 --- a/test/test_error_estimation.py +++ b/test/test_error_estimation.py @@ -1,4 +1,8 @@ +import unittest + from firedrake import * +from parameterized import parameterized + from goalie.error_estimation import ( form2indicator, get_dwr_indicator, @@ -6,8 +10,6 @@ from goalie.function_data import IndicatorData from goalie.go_mesh_seq import GoalOrientedMeshSeq from goalie.time_partition import TimeInstant, TimePartition -from parameterized import parameterized -import unittest class ErrorEstimationTestCase(unittest.TestCase): diff --git a/test/test_function_data.py b/test/test_function_data.py new file mode 100644 index 00000000..aa93c96c --- /dev/null +++ b/test/test_function_data.py @@ -0,0 +1,217 @@ +""" +Unit tests for :class:`~.FunctionData` and its subclasses. +""" + +import abc +import unittest + +from firedrake import * + +from goalie import * + + +class BaseTestCases: + """ + Class containing abstract base classes for unit testing subclasses of + :class:`~.FunctionData`. + """ + + class TestFunctionData(unittest.TestCase, abc.ABC): + """ + Base class for unit testing subclasses of :class:`~.FunctionData`. + """ + + def setUpUnsteady(self): + end_time = 1.0 + self.num_subintervals = 2 + timesteps = [0.5, 0.25] + self.field = "field" + self.num_exports = [1, 2] + self.mesh = UnitTriangleMesh() + self.time_partition = TimePartition( + end_time, self.num_subintervals, timesteps, self.field + ) + self.function_spaces = { + self.field: [ + FunctionSpace(self.mesh, "DG", 0) + for _ in range(self.num_subintervals) + ] + } + self._create_function_data() + + def setUpSteady(self): + end_time = 1.0 + self.num_subintervals = 1 + timesteps = [1.0] + self.field = "field" + self.num_exports = [1] + self.mesh = UnitTriangleMesh() + self.time_partition = TimePartition( + end_time, self.num_subintervals, timesteps, self.field + ) + self.function_spaces = { + self.field: [ + FunctionSpace(self.mesh, "DG", 0) + for _ in range(self.num_subintervals) + ] + } + self._create_function_data() + + @abc.abstractmethod + def _create_function_data(self): + pass + + def test_extract_by_field(self): + data = self.solution_data.extract(layout="field") + self.assertTrue(isinstance(data, AttrDict)) + self.assertTrue(self.field in data) + for label in self.labels: + self.assertTrue(isinstance(data[self.field], AttrDict)) + self.assertTrue(label in data[self.field]) + self.assertTrue(isinstance(data[self.field][label], list)) + self.assertEqual(len(data[self.field][label]), self.num_subintervals) + for i, num_exports in enumerate(self.num_exports): + self.assertTrue(isinstance(data[self.field][label][i], list)) + self.assertEqual(len(data[self.field][label][i]), num_exports) + for f in data[self.field][label][i]: + self.assertTrue(isinstance(f, Function)) + + def test_extract_by_label(self): + data = self.solution_data.extract(layout="label") + self.assertTrue(isinstance(data, AttrDict)) + for label in self.labels: + self.assertTrue(label in data) + self.assertTrue(isinstance(data[label], AttrDict)) + self.assertTrue(self.field in data[label]) + self.assertTrue(isinstance(data[label][self.field], list)) + self.assertEqual(len(data[label][self.field]), self.num_subintervals) + for i, num_exports in enumerate(self.num_exports): + self.assertTrue(isinstance(data[label][self.field][i], list)) + self.assertEqual(len(data[label][self.field][i]), num_exports) + for f in data[label][self.field][i]: + self.assertTrue(isinstance(f, Function)) + + def test_extract_by_subinterval(self): + data = self.solution_data.extract(layout="subinterval") + self.assertTrue(isinstance(data, list)) + self.assertEqual(len(data), self.num_subintervals) + for i, sub_data in enumerate(data): + self.assertTrue(isinstance(sub_data, AttrDict)) + self.assertTrue(self.field in sub_data) + self.assertTrue(isinstance(sub_data[self.field], AttrDict)) + for label in self.labels: + self.assertTrue(label in sub_data[self.field]) + self.assertTrue(isinstance(sub_data[self.field][label], list)) + self.assertEqual( + len(sub_data[self.field][label]), self.num_exports[i] + ) + for f in sub_data[self.field][label]: + self.assertTrue(isinstance(f, Function)) + + +class TestSteadyForwardSolutionData(BaseTestCases.TestFunctionData): + """ + Unit tests for :class:`~.ForwardSolutionData`. + """ + + def setUp(self): + super().setUpSteady() + self.labels = ("forward",) + + def _create_function_data(self): + self.solution_data = ForwardSolutionData( + self.time_partition, self.function_spaces + ) + + +class TestUnsteadyForwardSolutionData(BaseTestCases.TestFunctionData): + """ + Unit tests for :class:`~.ForwardSolutionData`. + """ + + def setUp(self): + super().setUpUnsteady() + self.labels = ("forward", "forward_old") + + def _create_function_data(self): + self.solution_data = ForwardSolutionData( + self.time_partition, self.function_spaces + ) + + +class TestSteadyAdjointSolutionData(BaseTestCases.TestFunctionData): + """ + Unit tests for :class:`~.AdjointSolutionData`. + """ + + def setUp(self): + super().setUpSteady() + self.labels = ("forward", "adjoint") + + def _create_function_data(self): + self.solution_data = AdjointSolutionData( + self.time_partition, self.function_spaces + ) + + +class TestUnsteadyAdjointSolutionData(BaseTestCases.TestFunctionData): + """ + Unit tests for :class:`~.AdjointSolutionData`. + """ + + def setUp(self): + super().setUpUnsteady() + self.labels = ("forward", "forward_old", "adjoint", "adjoint_next") + + def _create_function_data(self): + self.solution_data = AdjointSolutionData( + self.time_partition, self.function_spaces + ) + + +class TestIndicatorData(BaseTestCases.TestFunctionData): + """ + Unit tests for :class:`~.Indicatordata`. + """ + + def setUp(self): + super().setUpUnsteady() + self.labels = ("error_indicator",) + + def _create_function_data(self): + self.solution_data = IndicatorData( + self.time_partition, [self.mesh for _ in range(self.num_subintervals)] + ) + + def _test_extract_by_field_or_label(self, data): + self.assertTrue(isinstance(data, AttrDict)) + self.assertTrue(self.field in data) + self.assertEqual(len(data[self.field]), self.num_subintervals) + for i, num_exports in enumerate(self.num_exports): + self.assertTrue(isinstance(data[self.field][i], list)) + self.assertEqual(len(data[self.field][i]), num_exports) + for f in data[self.field][i]: + self.assertTrue(isinstance(f, Function)) + + def test_extract_by_field(self): + data = self.solution_data.extract(layout="field") + self._test_extract_by_field_or_label(data) + + def test_extract_by_label(self): + data = self.solution_data.extract(layout="label") + self._test_extract_by_field_or_label(data) + + def test_extract_by_subinterval(self): + data = self.solution_data.extract(layout="subinterval") + self.assertTrue(isinstance(data, list)) + self.assertEqual(len(data), self.num_subintervals) + for i, sub_data in enumerate(data): + self.assertTrue(isinstance(sub_data, AttrDict)) + self.assertTrue(self.field in sub_data) + self.assertTrue(isinstance(sub_data[self.field], list)) + for f in sub_data[self.field]: + self.assertTrue(isinstance(f, Function)) + + +if __name__ == "__main__": + unittest.main() diff --git a/test/test_interpolation.py b/test/test_interpolation.py deleted file mode 100644 index 99d449b1..00000000 --- a/test/test_interpolation.py +++ /dev/null @@ -1,170 +0,0 @@ -""" -Test interpolation schemes. -""" - -from firedrake import * -from goalie import * -from goalie.utility import function2cofunction -import unittest - - -class TestProject(unittest.TestCase): - """ - Unit tests for mesh-to-mesh projection. - """ - - def setUp(self): - self.source_mesh = UnitSquareMesh(1, 1, diagonal="left") - self.target_mesh = UnitSquareMesh(1, 1, diagonal="right") - - def sinusoid(self, source=True): - x, y = SpatialCoordinate(self.source_mesh if source else self.target_mesh) - return sin(pi * x) * sin(pi * y) - - def test_notimplemented_error(self): - Vs = FunctionSpace(self.source_mesh, "CG", 1) - Vt = FunctionSpace(self.target_mesh, "CG", 1) - with self.assertRaises(NotImplementedError) as cm: - project(2 * Function(Vs), Vt) - msg = "Can only currently project Functions and Cofunctions." - self.assertEqual(str(cm.exception), msg) - - def test_no_sub_source_space(self): - Vs = FunctionSpace(self.source_mesh, "CG", 1) - Vt = FunctionSpace(self.target_mesh, "CG", 1) - Vt = Vt * Vt - with self.assertRaises(ValueError) as cm: - project(Function(Vs), Function(Vt)) - msg = "Target space has multiple components but source space does not." - self.assertEqual(str(cm.exception), msg) - - def test_no_sub_source_space_adjoint(self): - Vs = FunctionSpace(self.source_mesh, "CG", 1) - Vt = FunctionSpace(self.target_mesh, "CG", 1) - Vt = Vt * Vt - with self.assertRaises(ValueError) as cm: - project(Cofunction(Vs.dual()), Cofunction(Vt.dual())) - msg = "Source space has multiple components but target space does not." - self.assertEqual(str(cm.exception), msg) - - def test_no_sub_target_space(self): - Vs = FunctionSpace(self.source_mesh, "CG", 1) - Vt = FunctionSpace(self.target_mesh, "CG", 1) - Vs = Vs * Vs - with self.assertRaises(ValueError) as cm: - project(Function(Vs), Function(Vt)) - msg = "Source space has multiple components but target space does not." - self.assertEqual(str(cm.exception), msg) - - def test_no_sub_target_space_adjoint(self): - Vs = FunctionSpace(self.source_mesh, "CG", 1) - Vt = FunctionSpace(self.target_mesh, "CG", 1) - Vs = Vs * Vs - with self.assertRaises(ValueError) as cm: - project(Cofunction(Vs.dual()), Cofunction(Vt.dual())) - msg = "Target space has multiple components but source space does not." - self.assertEqual(str(cm.exception), msg) - - def test_wrong_number_sub_spaces(self): - P1 = FunctionSpace(self.source_mesh, "CG", 1) - P0 = FunctionSpace(self.target_mesh, "DG", 0) - Vs = P1 * P1 * P1 - Vt = P0 * P0 - with self.assertRaises(ValueError) as cm: - project(Function(Vs), Function(Vt)) - msg = "Inconsistent numbers of components in source and target spaces: 3 vs. 2." - self.assertEqual(str(cm.exception), msg) - - def test_project_same_space(self): - Vs = FunctionSpace(self.source_mesh, "CG", 1) - source = assemble(interpolate(self.sinusoid(), Vs)) - target = Function(Vs) - project(source, target) - expected = source - self.assertAlmostEqual(errornorm(expected, target), 0) - - def test_project_same_space_adjoint(self): - Vs = FunctionSpace(self.source_mesh, "CG", 1) - source = assemble(interpolate(self.sinusoid(), Vs)) - source = function2cofunction(source) - target = Cofunction(Vs.dual()) - project(source, target) - expected = source - self.assertAlmostEqual(errornorm(expected, target), 0) - - def test_project_same_space_mixed(self): - P1 = FunctionSpace(self.source_mesh, "CG", 1) - Vs = P1 * P1 - source = Function(Vs) - s1, s2 = source.subfunctions - s1.interpolate(self.sinusoid()) - s2.interpolate(-self.sinusoid()) - target = Function(Vs) - project(source, target) - expected = source - self.assertAlmostEqual(errornorm(expected, target), 0) - - def test_project_same_space_mixed_adjoint(self): - P1 = FunctionSpace(self.source_mesh, "CG", 1) - Vs = P1 * P1 - source = Function(Vs) - s1, s2 = source.subfunctions - s1.interpolate(self.sinusoid()) - s2.interpolate(-self.sinusoid()) - source = function2cofunction(source) - target = Cofunction(Vs.dual()) - project(source, target) - expected = source - self.assertAlmostEqual(errornorm(expected, target), 0) - - def test_project_same_mesh(self): - Vs = FunctionSpace(self.source_mesh, "CG", 1) - Vt = FunctionSpace(self.source_mesh, "DG", 0) - source = assemble(interpolate(self.sinusoid(), Vs)) - target = Function(Vt) - project(source, target) - expected = Function(Vt).project(source) - self.assertAlmostEqual(errornorm(expected, target), 0) - - def test_project_same_mesh_adjoint(self): - Vs = FunctionSpace(self.source_mesh, "CG", 1) - Vt = FunctionSpace(self.source_mesh, "DG", 0) - source = assemble(interpolate(self.sinusoid(), Vs)) - target = Cofunction(Vt.dual()) - project(function2cofunction(source), target) - expected = function2cofunction(Function(Vt).project(source)) - self.assertAlmostEqual(errornorm(expected, target), 0) - - def test_project_same_mesh_mixed(self): - P1 = FunctionSpace(self.source_mesh, "CG", 1) - P0 = FunctionSpace(self.source_mesh, "DG", 0) - Vs = P1 * P1 - Vt = P0 * P0 - source = Function(Vs) - s1, s2 = source.subfunctions - s1.interpolate(self.sinusoid()) - s2.interpolate(-self.sinusoid()) - target = Function(Vt) - project(source, target) - expected = Function(Vt) - e1, e2 = expected.subfunctions - e1.project(s1) - e2.project(s2) - self.assertAlmostEqual(errornorm(expected, target), 0) - - def test_project_same_mesh_mixed_adjoint(self): - P1 = FunctionSpace(self.source_mesh, "CG", 1) - P0 = FunctionSpace(self.source_mesh, "DG", 0) - Vs = P1 * P1 - Vt = P0 * P0 - source = Function(Vs) - s1, s2 = source.subfunctions - s1.interpolate(self.sinusoid()) - s2.interpolate(-self.sinusoid()) - target = Cofunction(Vt.dual()) - project(function2cofunction(source), target) - expected = Function(Vt) - e1, e2 = expected.subfunctions - e1.project(s1) - e2.project(s2) - self.assertAlmostEqual(errornorm(expected, target), 0) diff --git a/test/test_math.py b/test/test_math.py index d61aa790..50f6ade9 100644 --- a/test/test_math.py +++ b/test/test_math.py @@ -1,8 +1,10 @@ -from firedrake import UnitTriangleMesh -from goalie.math import * +import unittest + import numpy as np import scipy as sp -import unittest +from firedrake import UnitTriangleMesh + +from goalie.math import * class TestBessel(unittest.TestCase): diff --git a/test/test_mesh_seq.py b/test/test_mesh_seq.py index 8ac31590..4cbadf5c 100644 --- a/test/test_mesh_seq.py +++ b/test/test_mesh_seq.py @@ -1,14 +1,23 @@ """ Testing for the mesh sequence objects. """ -from goalie.mesh_seq import MeshSeq -from goalie.time_partition import TimePartition, TimeInterval -from firedrake import * -from pyadjoint.block_variable import BlockVariable + import re import unittest from unittest.mock import patch +from firedrake import ( + Function, + FunctionSpace, + UnitCubeMesh, + UnitSquareMesh, + UnitTriangleMesh, +) +from pyadjoint.block_variable import BlockVariable + +from goalie.mesh_seq import MeshSeq +from goalie.time_partition import TimeInterval, TimePartition + class TestGeneric(unittest.TestCase): """ diff --git a/test/test_metric_normalise.py b/test/test_metric.py similarity index 63% rename from test/test_metric_normalise.py rename to test/test_metric.py index 081c6741..e877596c 100644 --- a/test/test_metric_normalise.py +++ b/test/test_metric.py @@ -1,31 +1,67 @@ """ Test metric normalisation functionality. """ -from firedrake import * -from goalie import * + +import unittest + +import numpy as np +import pytest from animate.metric import RiemannianMetric +from animate.utility import errornorm +from firedrake import * +from parameterized import parameterized from sensors import * from utility import uniform_mesh -import numpy as np -from parameterized import parameterized -import pytest -import unittest + +from goalie import * -class TestMetricNormalisation(unittest.TestCase): +class BaseClasses: + """ + Base classes for unit testing. + """ + + class MetricTestCase(unittest.TestCase): + """ + Base class for metric unit tests. + """ + + def setUp(self): + self.P1_ten = TensorFunctionSpace(uniform_mesh(2, 1), "CG", 1) + + @property + def simple_metric(self): + return RiemannianMetric(self.P1_ten) + + +class TestEnforceVariableConstraints(BaseClasses.MetricTestCase): + """ + Unit tests for enforcing variable constraints. + """ + + def test_defaults(self): + enforce_variable_constraints(self.simple_metric) + + def test_noniterable(self): + hmin, hmax, amax = 1.0e-05, 1.0, 1.0e05 + iterable = enforce_variable_constraints( + [self.simple_metric], h_min=[hmin], h_max=[hmax], a_max=[amax] + ) + noniterable = enforce_variable_constraints( + self.simple_metric, h_min=hmin, h_max=hmax, a_max=amax + ) + self.assertAlmostEqual(errornorm(iterable[0], noniterable[0]), 0) + + +class TestMetricNormalisation(BaseClasses.MetricTestCase): """ Unit tests for metric normalisation. """ def setUp(self): + super().setUp() self.time_partition = TimeInterval(1.0, 1.0, "u") - @property - def simple_metric(self): - mesh = uniform_mesh(2, 1) - P1_ten = TensorFunctionSpace(mesh, "CG", 1) - return RiemannianMetric(P1_ten) - def test_time_partition_length_error(self): time_partition = TimePartition(1.0, 2, [0.5, 0.5], "u") mp = {"dm_plex_metric_target_complexity": 1.0} @@ -132,3 +168,47 @@ def test_consistency(self, sensor, degree): # Check that the metrics coincide self.assertAlmostEqual(errornorm(M, M_st), 0) + + +class TestRampComplexity(unittest.TestCase): + """ + Unit tests for :func:`ramp_complexity`. + """ + + def test_base_non_positive_error(self): + with self.assertRaises(ValueError) as cm: + ramp_complexity(0, 100, 1) + msg = "Base complexity must be positive, not 0." + self.assertEqual(str(cm.exception), msg) + + def test_target_non_positive_error(self): + with self.assertRaises(ValueError) as cm: + ramp_complexity(100, 0, 1) + msg = "Target complexity must be positive, not 0." + self.assertEqual(str(cm.exception), msg) + + def test_iteration_non_positive_error(self): + with self.assertRaises(ValueError) as cm: + ramp_complexity(100, 1000, -1) + msg = "Current iteration must be non-negative, not -1." + self.assertEqual(str(cm.exception), msg) + + def test_num_iterations_non_positive_error(self): + with self.assertRaises(ValueError) as cm: + ramp_complexity(100, 1000, 0, num_iterations=-1) + msg = "Number of iterations must be non-negative, not -1." + self.assertEqual(str(cm.exception), msg) + + def test_instant_ramp(self): + C = ramp_complexity(100, 1000, 0, num_iterations=0) + self.assertEqual(C, 1000) + + def test_ramp(self): + base = 100 + target = 1000 + niter = 3 + for i in range(niter): + C = ramp_complexity(base, target, i, num_iterations=niter) + self.assertAlmostEqual(C, base + i * (target - base) / niter) + C = ramp_complexity(base, target, niter, num_iterations=niter) + self.assertEqual(C, target) diff --git a/test/test_metric_utils.py b/test/test_metric_utils.py deleted file mode 100644 index 362b2613..00000000 --- a/test/test_metric_utils.py +++ /dev/null @@ -1,46 +0,0 @@ -from goalie import * -import unittest - - -class TestRampComplexity(unittest.TestCase): - """ - Unit tests for :func:`ramp_complexity`. - """ - - def test_base_non_positive_error(self): - with self.assertRaises(ValueError) as cm: - ramp_complexity(0, 100, 1) - msg = "Base complexity must be positive, not 0." - self.assertEqual(str(cm.exception), msg) - - def test_target_non_positive_error(self): - with self.assertRaises(ValueError) as cm: - ramp_complexity(100, 0, 1) - msg = "Target complexity must be positive, not 0." - self.assertEqual(str(cm.exception), msg) - - def test_iteration_non_positive_error(self): - with self.assertRaises(ValueError) as cm: - ramp_complexity(100, 1000, -1) - msg = "Current iteration must be non-negative, not -1." - self.assertEqual(str(cm.exception), msg) - - def test_num_iterations_non_positive_error(self): - with self.assertRaises(ValueError) as cm: - ramp_complexity(100, 1000, 0, num_iterations=-1) - msg = "Number of iterations must be non-negative, not -1." - self.assertEqual(str(cm.exception), msg) - - def test_instant_ramp(self): - C = ramp_complexity(100, 1000, 0, num_iterations=0) - self.assertEqual(C, 1000) - - def test_ramp(self): - base = 100 - target = 1000 - niter = 3 - for i in range(niter): - C = ramp_complexity(base, target, i, num_iterations=niter) - self.assertAlmostEqual(C, base + i * (target - base) / niter) - C = ramp_complexity(base, target, niter, num_iterations=niter) - self.assertEqual(C, target) diff --git a/test/test_options.py b/test/test_options.py index 4543537d..1a9a18d5 100644 --- a/test/test_options.py +++ b/test/test_options.py @@ -1,8 +1,10 @@ -from firedrake import TensorFunctionSpace +import unittest + from animate.adapt import RiemannianMetric -from goalie.options import * +from firedrake import TensorFunctionSpace from utility import uniform_mesh -import unittest + +from goalie.options import * class TestAdaptParameters(unittest.TestCase): diff --git a/test/test_time_partition.py b/test/test_time_partition.py index 0fc24327..6a0abf88 100644 --- a/test/test_time_partition.py +++ b/test/test_time_partition.py @@ -1,9 +1,11 @@ """ Testing for the time partition objects. """ -from goalie.time_partition import TimePartition, TimeInterval, TimeInstant + import unittest +from goalie.time_partition import TimeInstant, TimeInterval, TimePartition + class TestSetup(unittest.TestCase): """ diff --git a/test/test_utility.py b/test/test_utility.py index 97268194..097eb9cb 100644 --- a/test/test_utility.py +++ b/test/test_utility.py @@ -2,17 +2,15 @@ Test utility functions. """ -from firedrake import * -from firedrake.norms import norm as fnorm -from firedrake.norms import errornorm as ferrnorm -from goalie import * -from utility import uniform_mesh import os import pathlib -from parameterized import parameterized import shutil import unittest +from firedrake import * +from utility import uniform_mesh + +from goalie import * pointwise_norm_types = [["l1"], ["l2"], ["linf"]] integral_scalar_norm_types = [["L1"], ["L2"], ["L4"], ["H1"], ["HCurl"]] @@ -23,226 +21,6 @@ # --------------------------- -class TestVTK(unittest.TestCase): - """ - Test the subclass of Firedrake's :class:`VTKFile`. - """ - - def setUp(self): - mesh = UnitSquareMesh(1, 1) - self.fs = FunctionSpace(mesh, "CG", 1) - pwd = os.path.dirname(__file__) - self.fname = os.path.join(pwd, "tmp.pvd") - - def tearDown(self): - fname = os.path.splitext(self.fname)[0] - for ext in (".pvd", "_0.vtu", "_1.vtu"): - if os.path.exists(fname + ext): - os.remove(fname + ext) - - def test_adaptive(self): - file = VTKFile(self.fname) - self.assertTrue(os.path.exists(self.fname)) - self.assertTrue(file._adaptive) - - def test_different_fnames(self): - f = Function(self.fs, name="f") - g = Function(self.fs, name="g") - file = VTKFile(self.fname) - file.write(f) - file.write(g) - self.assertEqual("f", g.name()) - - def test_different_lengths(self): - f = Function(self.fs, name="f") - g = Function(self.fs, name="g") - file = VTKFile(self.fname) - file.write(f) - with self.assertRaises(ValueError) as cm: - file.write(f, g) - msg = "Writing different number of functions: expected 1, got 2." - self.assertEqual(str(cm.exception), msg) - - -class TestMassMatrix(unittest.TestCase): - """ - Unit tests for :func:`assemble_mass_matrix`. - """ - - @parameterized.expand([["L2"], ["H1"]]) - def test_tiny(self, norm_type): - mesh = uniform_mesh(2, 1) - V = FunctionSpace(mesh, "DG", 0) - M = assemble_mass_matrix(V, norm_type=norm_type) - expected = np.array([[0.5, 0], [0, 0.5]]) - got = M.convert("dense").getDenseArray() - self.assertTrue(np.allclose(expected, got)) - - def test_norm_type_error(self): - mesh = uniform_mesh(2, 1) - V = FunctionSpace(mesh, "DG", 0) - with self.assertRaises(ValueError) as cm: - assemble_mass_matrix(V, norm_type="HDiv") - msg = "Norm type 'HDiv' not recognised." - self.assertEqual(str(cm.exception), msg) - - -class TestNorm(unittest.TestCase): - """ - Unit tests for :func:`norm`. - """ - - def setUp(self): - self.mesh = uniform_mesh(2, 4) - self.x, self.y = SpatialCoordinate(self.mesh) - V = FunctionSpace(self.mesh, "CG", 1) - self.f = assemble(interpolate(self.x**2 + self.y, V)) - - def test_boundary_error(self): - with self.assertRaises(NotImplementedError) as cm: - norm(self.f, norm_type="l1", boundary=True) - msg = "lp errors on the boundary not yet implemented." - self.assertEqual(str(cm.exception), msg) - - def test_l1(self): - expected = np.sum(np.abs(self.f.dat.data)) - got = norm(self.f, norm_type="l1") - self.assertAlmostEqual(expected, got) - - def test_l2(self): - expected = np.sqrt(np.sum(self.f.dat.data**2)) - got = norm(self.f, norm_type="l2") - self.assertAlmostEqual(expected, got) - - def test_linf(self): - expected = np.max(self.f.dat.data) - got = norm(self.f, norm_type="linf") - self.assertAlmostEqual(expected, got) - - def test_Linf(self): - expected = np.max(self.f.dat.data) - got = norm(self.f, norm_type="Linf") - self.assertAlmostEqual(expected, got) - - def test_notimplemented_lp_error(self): - with self.assertRaises(NotImplementedError) as cm: - norm(self.f, norm_type="lp") - msg = "lp norm of order p not supported." - self.assertEqual(str(cm.exception), msg) - - def test_L0_error(self): - with self.assertRaises(ValueError) as cm: - norm(self.f, norm_type="L0") - msg = "'L0' norm does not make sense." - self.assertEqual(str(cm.exception), msg) - - def test_notimplemented_Lp_error(self): - with self.assertRaises(ValueError) as cm: - norm(self.f, norm_type="Lp") - msg = "Don't know how to interpret 'Lp' norm." - self.assertEqual(str(cm.exception), msg) - - @parameterized.expand(integral_scalar_norm_types) - def test_consistency_firedrake(self, norm_type): - expected = fnorm(self.f, norm_type=norm_type) - got = norm(self.f, norm_type=norm_type) - self.assertAlmostEqual(expected, got) - - def test_consistency_hdiv(self): - V = VectorFunctionSpace(self.mesh, "CG", 1) - x, y = SpatialCoordinate(self.mesh) - f = assemble(interpolate(as_vector([y * y, -x * x]), V)) - expected = fnorm(f, norm_type="HDiv") - got = norm(f, norm_type="HDiv") - self.assertAlmostEqual(expected, got) - - def test_invalid_norm_type_error(self): - with self.assertRaises(ValueError) as cm: - norm(self.f, norm_type="X") - msg = "Unknown norm type 'X'." - self.assertEqual(str(cm.exception), msg) - - @parameterized.expand([("L1", 0.25), ("L2", 0.5)]) - def test_condition_integral(self, norm_type, expected): - self.f.assign(1.0) - condition = conditional(And(self.x < 0.5, self.y < 0.5), 1.0, 0.0) - val = norm(self.f, norm_type=norm_type, condition=condition) - self.assertAlmostEqual(val, expected) - - -class TestErrorNorm(unittest.TestCase): - """ - Unit tests for :func:`errornorm`. - """ - - def setUp(self): - self.mesh = uniform_mesh(2, 4) - self.x, self.y = SpatialCoordinate(self.mesh) - V = FunctionSpace(self.mesh, "CG", 1) - self.f = assemble(interpolate(self.x**2 + self.y, V)) - self.g = assemble(interpolate(self.x + self.y**2, V)) - - def test_shape_error(self): - with self.assertRaises(RuntimeError) as cm: - errornorm(self.f, self.mesh.coordinates) - msg = "Mismatching rank between u and uh." - self.assertEqual(str(cm.exception), msg) - - def test_not_function_error(self): - with self.assertRaises(TypeError) as cm: - errornorm(self.f, 1.0) - msg = "uh should be a Function, is a ''." - self.assertEqual(str(cm.exception), msg) - - def test_not_function_error_lp(self): - with self.assertRaises(TypeError) as cm: - errornorm(1.0, self.f, norm_type="l1") - msg = "u should be a Function, is a ''." - self.assertEqual(str(cm.exception), msg) - - def test_mixed_space_invalid_norm_error(self): - V = self.f.function_space() * self.f.function_space() - with self.assertRaises(NotImplementedError) as cm: - errornorm(Function(V), Function(V), norm_type="L1") - msg = "Norm type 'L1' not supported for mixed spaces." - self.assertEqual(str(cm.exception), msg) - - @parameterized.expand(scalar_norm_types) - def test_zero_scalar(self, norm_type): - err = errornorm(self.f, self.f, norm_type=norm_type) - self.assertAlmostEqual(err, 0.0) - - def test_zero_hdiv(self): - V = VectorFunctionSpace(self.mesh, "CG", 1) - x, y = SpatialCoordinate(self.mesh) - f = assemble(interpolate(as_vector([y * y, -x * x]), V)) - err = errornorm(f, f, norm_type="HDiv") - self.assertAlmostEqual(err, 0.0) - - @parameterized.expand(integral_scalar_norm_types) - def test_consistency_firedrake(self, norm_type): - expected = ferrnorm(self.f, self.g, norm_type=norm_type) - got = errornorm(self.f, self.g, norm_type=norm_type) - self.assertAlmostEqual(expected, got) - - def test_consistency_hdiv(self): - V = VectorFunctionSpace(self.mesh, "CG", 1) - x, y = SpatialCoordinate(self.mesh) - f = assemble(interpolate(as_vector([y * y, -x * x]), V)) - g = assemble(interpolate(as_vector([x * y, 1.0]), V)) - expected = ferrnorm(f, g, norm_type="HDiv") - got = errornorm(f, g, norm_type="HDiv") - self.assertAlmostEqual(expected, got) - - @parameterized.expand([("L1", 0.25), ("L2", 0.5)]) - def test_condition_integral(self, norm_type, expected): - self.f.assign(1.0) - self.g.assign(0.0) - condition = conditional(And(self.x < 0.5, self.y < 0.5), 1.0, 0.0) - val = errornorm(self.f, self.g, norm_type=norm_type, condition=condition) - self.assertAlmostEqual(val, expected) - - class TestEffectivityIndex(unittest.TestCase): """ Unit tests for :func:`effectivity_index`. diff --git a/test/utility.py b/test/utility.py index b139d7ba..5afc7eab 100644 --- a/test/utility.py +++ b/test/utility.py @@ -1,10 +1,12 @@ """ Functions used frequently for testing. """ + import firedrake -from goalie.metric import RiemannianMetric import ufl +from goalie.metric import RiemannianMetric + def uniform_mesh(dim, n, l=1, **kwargs): args = [n] * dim + [l] diff --git a/test_adjoint/conftest.py b/test_adjoint/conftest.py index a6b299b4..8a582c51 100644 --- a/test_adjoint/conftest.py +++ b/test_adjoint/conftest.py @@ -3,9 +3,11 @@ **Disclaimer: some functions copied from firedrake/src/tests/conftest.py """ + +from subprocess import check_call + import pyadjoint import pytest -from subprocess import check_call def parallel(item): @@ -83,9 +85,10 @@ def pytest_runtest_setup(item): if MPI.COMM_WORLD.size > 1: # Turn on source hash checking - from firedrake import parameters from functools import partial + from firedrake import parameters + def _reset(check): parameters["pyop2_options"]["check_src_hashes"] = check diff --git a/test_adjoint/examples/burgers.py b/test_adjoint/examples/burgers.py index ddb17414..e866ca57 100644 --- a/test_adjoint/examples/burgers.py +++ b/test_adjoint/examples/burgers.py @@ -10,7 +10,6 @@ from firedrake import * from firedrake.__future__ import interpolate - # Problem setup n = 32 mesh = UnitSquareMesh(n, n, diagonal="left") diff --git a/test_adjoint/examples/point_discharge2d.py b/test_adjoint/examples/point_discharge2d.py index 57d91098..034fe433 100644 --- a/test_adjoint/examples/point_discharge2d.py +++ b/test_adjoint/examples/point_discharge2d.py @@ -15,10 +15,11 @@ Paris: R&D, Electricite de France, p. 134 (2014). """ -from firedrake import * -from goalie.math import bessk0 + import numpy as np +from firedrake import * +from goalie.math import bessk0 # Problem setup n = 0 diff --git a/test_adjoint/examples/point_discharge3d.py b/test_adjoint/examples/point_discharge3d.py index e9124f31..82e8397d 100644 --- a/test_adjoint/examples/point_discharge3d.py +++ b/test_adjoint/examples/point_discharge3d.py @@ -24,10 +24,11 @@ Proceedings of the 28th International Meshing Roundtable (2020). """ -from firedrake import * -from goalie.math import bessk0 + import numpy as np +from firedrake import * +from goalie.math import bessk0 # Problem setup n = 0 diff --git a/test_adjoint/examples/steady_flow_past_cyl.py b/test_adjoint/examples/steady_flow_past_cyl.py index 7864a961..71f78ac6 100644 --- a/test_adjoint/examples/steady_flow_past_cyl.py +++ b/test_adjoint/examples/steady_flow_past_cyl.py @@ -11,10 +11,10 @@ https://nbviewer.jupyter.org/github/firedrakeproject/firedrake/blob/master/docs/notebooks/06-pde-constrained-optimisation.ipynb """ -from firedrake import * -from firedrake.__future__ import interpolate import os +from firedrake import * +from firedrake.__future__ import interpolate mesh = Mesh(os.path.join(os.path.dirname(__file__), "mesh-with-hole.msh")) fields = ["up"] diff --git a/test_adjoint/test_adjoint.py b/test_adjoint/test_adjoint.py index a7985cb7..1b2649bd 100644 --- a/test_adjoint/test_adjoint.py +++ b/test_adjoint/test_adjoint.py @@ -2,15 +2,17 @@ Test adjoint drivers. """ -from firedrake import * -from goalie_adjoint import * -import pyadjoint -import pytest import importlib import os import sys import unittest +import pyadjoint +import pytest +from animate.utility import errornorm, norm +from firedrake import Cofunction, UnitTriangleMesh + +from goalie_adjoint import * sys.path.append(os.path.join(os.path.dirname(__file__), "examples")) diff --git a/test_adjoint/test_demos.py b/test_adjoint/test_demos.py index c18499ac..fe8a0114 100644 --- a/test_adjoint/test_demos.py +++ b/test_adjoint/test_demos.py @@ -1,14 +1,15 @@ """ Checks that all demo scripts run. """ + import glob import os -from os.path import splitext -import pytest import shutil import subprocess import sys +from os.path import splitext +import pytest # Set environment variable to specify that tests are being run so that we can # cut down the length of the demos diff --git a/test_adjoint/test_fp_iteration.py b/test_adjoint/test_fp_iteration.py index 46cf8a4b..56bb5f87 100644 --- a/test_adjoint/test_fp_iteration.py +++ b/test_adjoint/test_fp_iteration.py @@ -1,10 +1,12 @@ -from firedrake import * -from goalie_adjoint import * import abc -from parameterized import parameterized import unittest from unittest.mock import MagicMock +from firedrake import * +from parameterized import parameterized + +from goalie_adjoint import * + def constant_qoi(mesh_seq, solutions, index): R = FunctionSpace(mesh_seq[index], "R", 0) diff --git a/test_adjoint/test_mesh_seq.py b/test_adjoint/test_mesh_seq.py index ee38b1cc..d3871484 100644 --- a/test_adjoint/test_mesh_seq.py +++ b/test_adjoint/test_mesh_seq.py @@ -1,17 +1,21 @@ """ Testing for the mesh sequence objects. """ + +import logging +import unittest + +import pyadjoint +import pytest +from animate.utility import norm from firedrake import * -from goalie_adjoint import * +from parameterized import parameterized + +from goalie.go_mesh_seq import GoalOrientedMeshSeq from goalie.log import * from goalie.mesh_seq import MeshSeq -from goalie.go_mesh_seq import GoalOrientedMeshSeq from goalie.time_partition import TimeInterval -from parameterized import parameterized -import logging -import pyadjoint -import pytest -import unittest +from goalie_adjoint import * class TestGetSolveBlocks(unittest.TestCase): diff --git a/test_adjoint/test_utils.py b/test_adjoint/test_utils.py index b572eef7..7be2c310 100644 --- a/test_adjoint/test_utils.py +++ b/test_adjoint/test_utils.py @@ -1,8 +1,10 @@ +import unittest + +import numpy as np from firedrake import * -from goalie_adjoint import * + from goalie.adjoint import annotate_qoi -import numpy as np -import unittest +from goalie_adjoint import * class TestAdjointUtils(unittest.TestCase): From a20c97535d381ffbab7aaf6a99e4715d7ea23df5 Mon Sep 17 00:00:00 2001 From: ddundo Date: Tue, 23 Apr 2024 09:32:45 +0000 Subject: [PATCH 09/10] #28: Remove get_bcs in demos --- demos/bubble_shear-goal_oriented.py | 13 +++---------- demos/bubble_shear.py | 17 +++++------------ 2 files changed, 8 insertions(+), 22 deletions(-) diff --git a/demos/bubble_shear-goal_oriented.py b/demos/bubble_shear-goal_oriented.py index adeddde0..d21d8c65 100644 --- a/demos/bubble_shear-goal_oriented.py +++ b/demos/bubble_shear-goal_oriented.py @@ -37,13 +37,6 @@ def get_function_spaces(mesh): return {"c": FunctionSpace(mesh, "CG", 1)} -def get_bcs(mesh_seq): - def bcs(index): - return [DirichletBC(mesh_seq.function_spaces["c"][index], 0.0, "on_boundary")] - - return bcs - - def ball_initial_condition(x, y): ball_r0, ball_x0, ball_y0 = 0.15, 0.5, 0.65 r = sqrt(pow(x - ball_x0, 2) + pow(y - ball_y0, 2)) @@ -61,7 +54,7 @@ def get_initial_condition(mesh_seq): def get_solver(mesh_seq): def solver(index, ic): tp = mesh_seq.time_partition - t_start, t_end = tp.subintervals[index] + t_end = tp.subintervals[index][1] dt = tp.timesteps[index] time = mesh_seq.get_time(index) @@ -82,7 +75,8 @@ def solver(index, ic): # We pass both the concentration and velocity Functions to get_form form_fields = {"c": (c, c_), "u": (u, u_)} F = mesh_seq.form(index, form_fields)["c"] - nlvp = NonlinearVariationalProblem(F, c, bcs=mesh_seq.bcs(index)) + bcs = DirichletBC(mesh_seq.function_spaces["c"][index], 0.0, "on_boundary") + nlvp = NonlinearVariationalProblem(F, c, bcs=bcs) nlvs = NonlinearVariationalSolver(nlvp, ad_block_tag="c") # Time integrate from t_start to t_end @@ -259,7 +253,6 @@ def adaptor(mesh_seq, solutions, indicators): meshes, get_function_spaces=get_function_spaces, get_initial_condition=get_initial_condition, - get_bcs=get_bcs, get_form=get_form, get_solver=get_solver, get_qoi=get_qoi, diff --git a/demos/bubble_shear.py b/demos/bubble_shear.py index 0d0c0bac..3fa5607d 100644 --- a/demos/bubble_shear.py +++ b/demos/bubble_shear.py @@ -60,20 +60,13 @@ def get_function_spaces(mesh): return {"c": FunctionSpace(mesh, "CG", 1)} -# We proceed similarly with prescribing initial and boundary conditions. At :math:`t=0`, -# we initialise the tracer concentration :math:`c_0 = c(x, y, 0)` to be :math:`1` inside +# 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 this is a discontinuous function # which will not be represented well on a coarse uniform mesh. :: -def get_bcs(mesh_seq): - def bcs(index): - return [DirichletBC(mesh_seq.function_spaces["c"][index], 0.0, "on_boundary")] - - return bcs - - def ball_initial_condition(x, y): ball_r0, ball_x0, ball_y0 = 0.15, 0.5, 0.65 r = sqrt(pow(x - ball_x0, 2) + pow(y - ball_y0, 2)) @@ -136,7 +129,7 @@ def form(index, form_fields): def get_solver(mesh_seq): def solver(index, ic): tp = mesh_seq.time_partition - t_start, t_end = tp.subintervals[index] + t_end = tp.subintervals[index][1] dt = tp.timesteps[index] time = mesh_seq.get_time(index) @@ -157,7 +150,8 @@ def solver(index, ic): # We pass both the concentration and velocity Functions to get_form form_fields = {"c": (c, c_), "u": (u, u_)} F = mesh_seq.form(index, form_fields)["c"] - nlvp = NonlinearVariationalProblem(F, c, bcs=mesh_seq.bcs(index)) + bcs = DirichletBC(mesh_seq.function_spaces["c"][index], 0.0, "on_boundary") + nlvp = NonlinearVariationalProblem(F, c, bcs=bcs) nlvs = NonlinearVariationalSolver(nlvp, ad_block_tag="c") # Time integrate from t_start to t_end @@ -203,7 +197,6 @@ def solver(index, ic): meshes, get_function_spaces=get_function_spaces, get_initial_condition=get_initial_condition, - get_bcs=get_bcs, get_form=get_form, get_solver=get_solver, ) From 29e959f5b82f28c496433347c3492008a868354c Mon Sep 17 00:00:00 2001 From: ddundo Date: Tue, 23 Apr 2024 09:37:29 +0000 Subject: [PATCH 10/10] #28: Recompile forms at every timestep in indicate_errors --- goalie/go_mesh_seq.py | 26 ++++++-------------------- 1 file changed, 6 insertions(+), 20 deletions(-) diff --git a/goalie/go_mesh_seq.py b/goalie/go_mesh_seq.py index 821bf5ea..0666bcbb 100644 --- a/goalie/go_mesh_seq.py +++ b/goalie/go_mesh_seq.py @@ -3,7 +3,6 @@ """ from collections.abc import Iterable -from inspect import signature import numpy as np import ufl @@ -159,9 +158,6 @@ def indicate_errors( enriched_mesh_seq.solve_adjoint(**solver_kwargs) tp = self.time_partition - # check whether get_form contains a kwarg indicating time-dependent constants - timedep_const = "err_ind_time" in signature(enriched_mesh_seq.form).parameters - FWD, ADJ = "forward", "adjoint" FWD_OLD = "forward" if self.steady else "forward_old" ADJ_NEXT = "adjoint" if self.steady else "adjoint_next" @@ -183,29 +179,19 @@ def indicate_errors( u_star_next[f] = Function(fs_e) u_star_e[f] = Function(fs_e) - # Get forms for each equation in enriched space - if not timedep_const: - forms = enriched_mesh_seq.form(i, mapping) - if not isinstance(forms, dict): - raise TypeError( - "The function defined by get_form should return a dictionary" - f", not type '{type(forms)}'." - ) - # Loop over each strongly coupled field for f in self.fields: # Loop over each timestep solutions = self.solutions.extract(layout="field") indicators = self.indicators.extract(layout="field") for j in range(len(solutions[f]["forward"][i])): + # Recompile the form with updated time-dependent constants + time.assign( + tp.subintervals[i][0] + + (j + 1) * tp.timesteps[i] * tp.num_timesteps_per_export[i] + ) + forms = enriched_mesh_seq.form(i, mapping) # Update fields - if timedep_const: - # recompile the form with updated time-dependent constants - time.assign( - tp.subintervals[i][0] - + (j + 1) * tp.timesteps[i] * tp.num_timesteps_per_export[i] - ) - forms = enriched_mesh_seq.form(i, mapping) transfer(self.solutions[f][FWD][i][j], u[f]) transfer(self.solutions[f][FWD_OLD][i][j], u_[f]) transfer(self.solutions[f][ADJ][i][j], u_star[f])