Convergence rate w.r.t. to a finer-grid solution

dolfinx
,
1

I use code from dolfinx tutorial and interpolation test to find convergence rates with respect to a finer-grid solution. It works perfectly for the Poisson problem given in the tutorial, but gives suspicious rates for a plate being pulled by a fixed displacement. Am I making some mistake here, or is it expected? Here is the code (change problem to run for the two problems ):

import numpy as np
import ufl
from dolfinx import default_scalar_type
from dolfinx.fem import (
    Function,
    assemble_scalar,
    create_interpolation_data,
    dirichletbc,
    form,
    functionspace,
    locate_dofs_geometrical,
    locate_dofs_topological,
)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import (
    CellType,
    create_unit_square,
    locate_entities_boundary,
)
from mpi4py import MPI
from ufl import (
    SpatialCoordinate,
    TestFunction,
    TrialFunction,
    div,
    dot,
    dx,
    grad,
    inner,
)

problem = "plate"  # 'poisson' | 'plate'


def u_ex(mod):
    return lambda x: mod.cos(2 * mod.pi * x[0]) * mod.cos(2 * mod.pi * x[1])


u_numpy = u_ex(np)
u_ufl = u_ex(ufl)


def map_uh_coarse_to_fine(uh_coarse, V_fine):
    V_coarse = uh_coarse.function_space
    domain_fine = V_fine.mesh
    domain_fine_cell_map = domain_fine.topology.index_map(domain_fine.topology.dim)
    num_cells_on_proc = (
        domain_fine_cell_map.size_local + domain_fine_cell_map.num_ghosts
    )
    cells = np.arange(num_cells_on_proc, dtype=np.int32)
    interpolation_data = create_interpolation_data(V_fine, V_coarse, cells)
    uh_coarse_on_fine = Function(V_fine)
    uh_coarse_on_fine.interpolate_nonmatching(
        uh_coarse, cells, interpolation_data=interpolation_data
    )
    uh_coarse_on_fine.x.scatter_forward()
    return uh_coarse_on_fine


def L2_err(uh_coarse, uh_fine):
    uh_coarse_on_fine = map_uh_coarse_to_fine(uh_coarse, uh_fine.function_space)
    comm = uh_coarse.function_space.mesh.comm
    error = form((uh_coarse_on_fine - uh_fine) ** 2 * dx)
    E = np.sqrt(comm.allreduce(assemble_scalar(error), MPI.SUM))
    if comm.rank == 0:
        return E


def solve_poisson(N=10, degree=1):
    mesh = create_unit_square(MPI.COMM_WORLD, N, N)
    x = SpatialCoordinate(mesh)
    f = -div(grad(u_ufl(x)))
    V = functionspace(mesh, ("Lagrange", degree))
    u = TrialFunction(V)
    v = TestFunction(V)
    a = inner(grad(u), grad(v)) * dx
    L = f * v * dx
    u_bc = Function(V)
    u_bc.interpolate(u_numpy)
    facets = locate_entities_boundary(
        mesh, mesh.topology.dim - 1, lambda x: np.full(x.shape[1], True)
    )
    dofs = locate_dofs_topological(V, mesh.topology.dim - 1, facets)
    bcs = [dirichletbc(u_bc, dofs)]
    default_problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "cg"})
    return default_problem.solve(), u_ufl(x)


def solve_plate(N=10, degree=1):
    uD = 0.1
    E = 10
    nu = 0.3
    lam = E * nu / (1 - nu**2)  # For plane-stress
    mu = E / 2 / (1 + nu)

    # Mesh
    mesh = create_unit_square(MPI.COMM_WORLD, N, N, CellType.quadrilateral)
    tdim = mesh.topology.dim
    fdim = tdim - 1

    # Function space
    V = functionspace(mesh, ("P", degree, (2,)))

    # Boundary conditions
    bottom_dofs = locate_dofs_geometrical(V, lambda x: np.isclose(x[1], 0))
    bottom_bc = dirichletbc(np.zeros((2,)), bottom_dofs, V)
    top_facets = locate_entities_boundary(mesh, fdim, lambda x: np.isclose(x[1], 1))
    top_ydofs = locate_dofs_topological(V.sub(1), fdim, top_facets)
    top_bc = dirichletbc(uD, top_ydofs, V.sub(1))
    bcs = [bottom_bc, top_bc]

    # Variational form
    def strain(u):
        return ufl.sym(grad(u))

    def stress(u):
        return lam * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * strain(u)

    u = TrialFunction(V)
    v = TestFunction(V)
    a = inner(stress(u), strain(v)) * dx
    f = Function(V)
    L = dot(f, v) * dx

    # Solve
    problem = LinearProblem(a, L, bcs)
    uh = problem.solve()

    return uh


def solve_problem(name, N, degree):
    if name == "poisson":
        return solve_poisson(N, degree)[0]
    elif name == "plate":
        return solve_plate(N, degree)
    else:
        assert False, 'Unknown problem name "{name}"!'


# Solve on a fine grid
uh_fine = solve_problem(problem, N=100, degree=1)

# Compute convergence rates
Ns = [4, 8, 16, 32, 64]
Es = np.zeros(len(Ns), dtype=default_scalar_type)
hs = np.zeros(len(Ns), dtype=np.float64)
for i, N in enumerate(Ns):
    uh_coarse = solve_problem(problem, N, degree=1)
    comm = uh_coarse.function_space.mesh.comm
    Es[i] = L2_err(uh_coarse, uh_fine)
    hs[i] = 1.0 / Ns[i]
    if comm.rank == 0:
        print(f"h: {hs[i]:.2e} Error: {Es[i]:.2e}")

rates = np.log(Es[1:] / Es[:-1]) / np.log(hs[1:] / hs[:-1])
if comm.rank == 0:
    print(f"Rates: {rates}")

for problem=’poisson’ I get Rates: [1.61683649 1.9088231 2.02968856 2.15457504], but for problem=’plate’ I get Rates: [1.3495696 1.13147778 0.90571456 0.16758422].

Thanks.

2

When considering convergence rates, it is very important to ensure that the solver converge to a tolerance that is sufficient to measure convergence with.
Usually, this is done by enforcing a direct solver, rather than the iterative (which is the PETSc default).
i.e.

    petsc_options = {"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps",
                     "ksp_monitor": None, "ksp_error_if_not_converged": True}

    problem = LinearProblem(a, L, bcs=bcs, petsc_options=petsc_options)

results in

Rates: [1.41327707 1.45397808 1.60882883 2.09653075]
2 Likes
3

Thanks Dokken. I wouldn’t have got it by myself, especially because Poisson problem was working with the same code. Could you please direct me to some source where I can read about petsc_options in detail?

4

PETSc is a toolkit with its own documentation. You can read about how it works and how it’s designed to use options for initialising its own features. For example, the documentation for KSPSetTypelists its options database key, with lots of types to choose from.

1 Like