Summary
I’m trying to learn about local a posteriori error estimation, and automatic mesh refinement, using the Poisson problem as a test case.
I have some questions about how to best define the local problem using fenicsx. I’ve found previous discussions regarding local problems in the form of
integrating over an interior surface, as well as previous implementations of error control, but these regard the previous version of dolfin, and make use of some deprecated mesh-manipulation functions (ex: MeshFunction(), facet()).
Right now I’m trying to integrate some form on each cell, and its corresponding edges. My specific questions are:
- What is the most efficient way to define the local problem? Currently I am redefining/compiling it for each cell, within a loop. This seems too costly to repeat for each adapted mesh.
- Does my setup of integrals over the edges of each cell look correct here? I’ve tried to assign the index of each edge as its meshtag, and then define a measure on each edge. Each of these is then mapped to a cell using the domain topology cell-to-edges indexing.
Code
Below is a MWE. This version does not include the mesh adaptation code, just the residual/error calculations.
from mpi4py import MPI # Parallel communication
from dolfinx import mesh # Mesh tools
import dolfinx as dolx
from dolfinx import fem
from dolfinx.fem import FunctionSpace
import numpy as np
import ufl
from ufl import (grad, div, dot, inner, jump, FacetNormal, dx,
TrialFunction, TestFunction)
from petsc4py.PETSc import ScalarType
import matplotlib.pyplot as plt
# Domain and Function Space
nx = 10
ny = 10
domain = mesh.create_unit_square(MPI.COMM_WORLD, nx, ny, mesh.CellType.triangle)
func_space = ("CG", 2)
V = FunctionSpace(domain, func_space)
#% Define Solution/Basis functions
u = TrialFunction(V)
v = TestFunction(V)
#% BC
uD = fem.Function(V)
bound_func = lambda x: 1 + x[0]**2 + 2 * x[1]**2
uD.interpolate(bound_func)
# Define Boundary
# Create facet to cell connectivity required to determine boundary facets
tdim = domain.topology.dim # Topology dimension
fdim = tdim - 1 # Facet dimensio
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
# Specify BC
boundary_dofs = fem.locate_dofs_topological(V, fdim, boundary_facets)
bc = fem.dirichletbc(uD, boundary_dofs)
#% Compute exact solution
V2 = fem.FunctionSpace(domain, ("CG", 2)) # Error space is second order?
uex = fem.Function(V2)
uex.interpolate(bound_func)
#% Define Variational Form
f = fem.Constant(domain, ScalarType(-6))
a = -dot(grad(u), grad(v)) * dx
L = f * v * dx
# Use solver
problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
#% Compute Error
V2 = fem.FunctionSpace(domain, ("CG", 2)) # Error space is second order?
uex = fem.Function(V2)
uex.interpolate(bound_func)
# L2
L2_expr = fem.form(ufl.inner(uh - uex, uh - uex) * ufl.dx) # L2 norm in ufl
error_local = fem.assemble_scalar(L2_expr) # Compute value from ufl
error_gather = domain.comm.allreduce(error_local, op=MPI.SUM) # Gather results from other processors
error_L2 = np.sqrt(error_gather)
# Max point-wise
error_max = np.max(np.abs(uD.x.array-uh.x.array))
# Only print the error on one process
if domain.comm.rank == 0:
print(f"Error_L2 : {error_L2:.2e}")
print(f"Error_max : {error_max:.2e}")
#%% Mesh manipulation (local problem definition)
# Init
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(tdim, 2) # Compute cells
c_to_e = domain.topology.connectivity(tdim, 1) # Map cells to edges
e_to_c = domain.topology.connectivity(1, 2) # Map edges to cells
e_to_v = domain.topology.connectivity(1, 0) # Map edges to vertices
# Tag each cell
num_cells = domain.topology.index_map(tdim).size_local + domain.topology.index_map(tdim).num_ghosts
tag_list = np.arange(num_cells, dtype=np.int32)
cell_geo = dolx.cpp.mesh.entities_to_geometry(domain, tdim, tag_list, False)
cell_tags = mesh.meshtags(domain, tdim, tag_list, tag_list)
# Tag each edge
num_edges = domain.topology.index_map(1).size_local + domain.topology.index_map(1).num_ghosts
tag_list = np.arange(num_edges, dtype=np.int32)
edge_geo = dolx.cpp.mesh.entities_to_geometry(domain, 1, tag_list, False)
edge_tags = mesh.meshtags(domain, 1, tag_list, tag_list)
# Define local measures based on previous tags, for local problem definitions
dA = ufl.Measure("dx", domain, subdomain_data=cell_tags)
dE = ufl.Measure("dS", domain, subdomain_data=edge_tags)
#%% Calculating Residual
# Some operators for the form
n = FacetNormal(domain)
lap_uh = div(grad(uh))
# Entity Sizes
h_arr = dolx.cpp.mesh.h(domain, tdim, range(num_cells)) # Cell sizes
L_arr = dolx.cpp.mesh.h(domain, tdim-1, range(num_edges))
# Compute Local Residuals
Ri_arr = np.zeros(num_cells) # Interior Residuals (cells)
Re_arr = np.zeros(num_cells) # Exterior Residuals (edges)
eta_arr = np.zeros(num_cells) # Error Predictor
E_arr = np.zeros(num_cells) # Actual Error
for itr in range(num_cells):
# Compute Interior Residual
Ri_form = ( inner(lap_uh, lap_uh) - f )*dA(itr)
Ri_val = fem.assemble_scalar( fem.form(Ri_form) )
Ri_arr[itr] = Ri_val
# Compute Exterior Residual
Re_val = 0
for itrE in c_to_e.links(itr):
Re_form = -jump( dot(grad(uh), n) )*dE(itrE)
Re_val += np.sqrt(L_arr[itrE])*fem.assemble_scalar( fem.form(Re_form) )
Re_arr[itr] = Re_val
# Error Estimator
eta_arr[itr] = np.sqrt( h_arr[itr]**2 * Ri_val**2 + Re_val**2 )
# True Error
L2_expr = fem.form(ufl.inner(uh - uex, uh - uex) * dA(itr)) # L2 norm in ufl
error_local = fem.assemble_scalar(L2_expr) # Compute value from ufl
error_gather = domain.comm.allreduce(error_local, op=MPI.SUM) # Gather results from other processors
error_L2 = np.sqrt(error_gather)
E_arr[itr] = error_L2
#%% Plot Residual
plt.plot( abs(eta_arr) / np.max(eta_arr))
plt.plot( abs(E_arr) / np.max(E_arr))
plt.xlabel('Cell Number')
plt.ylabel('Normalized Errors')
plt.grid()
plt.legend( ['Residual', 'L2 Error'] )