Quadrature Family - Nonlinear Poisson Example

Dear All,

I was studying the use of quadrature elements in Fenics and tried to solve the nonlinear poisson example - \nabla \cdot (q(u) \nabla u) = f where q(u) = (1+u)^2 using the implementation suggested here and an approach using quadrature elements (see link to ufl doc)

In the ufl doc it is stated that

Clearly, using the standard FiniteElement (for user-defined functions C and \sigma_0) leads to a poor convergence whereas the QuadratureElement still leads to quadratic convergence

It seems I can’t reproduce these results. I get equally good convergence results whether I use “Quadrature”, “DG” or even “CG” element families.

Sorry in advance, but I don’t really know how to shorten my code. I’ll put the script here so you might be able to try it out yourself.

Is the statement in the ufl doc wrong or did I miss something here?

Any help is much appreciated.

#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""Solve the nonlinear poisson eq.
     - nabla · (q(u) nabla u) = f

    nonlinear_poisson.py [options] NELEMENTS M ORDER
    nonlinear_poisson.py -h | --help

    NELEMENTS           Number of elements in each spatial direction.
    M                   q(u) = (1 + u)^M
    ORDER               Order of Lagrange polynomial.

    --family=FAMILY     Element family to be used for history variables. [default: Quadrature]
    --summary           Print summary of results to the screen.
    --save              Save (append) summary of results to file.

import sys
import dolfin as df
import numpy as np
import warnings
from docopt import docopt
from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning

INFO = 20
warnings.simplefilter("ignore", QuadratureRepresentationDeprecationWarning)

def local_project(v, V, u=None):
    """project v onto V and store the values in u"""
    QUAD_DEG = 2
    df.parameters["form_compiler"]["representation"] = "quadrature"
    metadata = {"quadrature_degree": QUAD_DEG, "quadrature_scheme": "default"}
    dx = df.dx(metadata=metadata)
    v_trial = df.TrialFunction(V)
    v_test = df.TestFunction(V)
    a_proj = df.inner(v_trial, v_test) * dx
    b_proj = df.inner(v, v_test) * dx
    solver = df.LocalSolver(a_proj, b_proj)
    if u is None:
        u = df.Function(V)
        df.parameters["form_compiler"]["representation"] = "uflacs"
        return u
        df.parameters["form_compiler"]["representation"] = "uflacs"

def parse_arguments(args):
    args = docopt(__doc__, args)
    args['NELEMENTS'] = int(args['NELEMENTS'])
    args['M'] = int(args['M'])
    args['ORDER'] = int(args['ORDER'])
    assert args['--family'] in {"Quadrature", "DG", "CG"}
    if args['--save']:
        assert args['--summary'], "Can't save without --summary option."
    return args

def main(args):
    args = parse_arguments(args)
    mesh = df.UnitSquareMesh(args['NELEMENTS'], args['NELEMENTS'])
    V = df.FunctionSpace(mesh, "Lagrange", args['ORDER'])
    if args['--family'] == "Quadrature":
        QUAD_DEG = 2
        metadata = {"quadrature_degree": QUAD_DEG, "quadrature_scheme": "default"}
        dx = df.dx(metadata=metadata)
        QElement = df.FiniteElement(args['--family'], mesh.ufl_cell(), degree=QUAD_DEG, quad_scheme="default")
        Q = df.FunctionSpace(mesh, QElement)
        QVectorElement = df.VectorElement(args['--family'], mesh.ufl_cell(), degree=QUAD_DEG,
        QV = df.FunctionSpace(mesh, QVectorElement)
    elif args['--family'] == "DG":
        dx = df.dx
        Q = df.FunctionSpace(mesh, args['--family'], args['ORDER']-1)
        QV = df.VectorFunctionSpace(mesh, args['--family'], args['ORDER']-1)
    elif args['--family'] == "CG":
        dx = df.dx
        Q = df.FunctionSpace(mesh, args['--family'], args['ORDER'])
        QV = df.VectorFunctionSpace(mesh, args['--family'], args['ORDER'])
        assert Q.ufl_element() == V.ufl_element()

    tol = 1e-14
    def left_boundary(x, on_boundary):
        return on_boundary and abs(x[0]) < tol

    def right_boundary(x, on_boundary):
        return on_boundary and abs(x[0]-1) < tol

    Gamma_0 = df.DirichletBC(V, df.Constant(0.0), left_boundary)
    Gamma_1 = df.DirichletBC(V, df.Constant(1.0), right_boundary)
    bcs = [Gamma_0, Gamma_1]

    # Compute Initial guess for q(u)=1
    u = df.TrialFunction(V)
    v = df.TestFunction(V)
    a = df.inner(df.grad(u), df.grad(v)) * dx
    f = df.Constant(0.0)
    L = f * v * dx
    A, b = df.assemble_system(a, L, bcs)
    u0 = df.Function(V)
    df.solve(A, u0.vector(), b)

    Gamma_0_du = df.DirichletBC(V, df.Constant(0.0), left_boundary)
    Gamma_1_du = df.DirichletBC(V, df.Constant(0.0), right_boundary)
    bcs_du = [Gamma_0_du, Gamma_1_du]

    def q(u, m=args['M']):
        return (1+u)**m

    def Dq(u, m=args['M']):
        return m*(1+u)**(m-1)

    du = df.TrialFunction(V)# u = u0 + omega*du
    C = df.Function(Q)
    sigma0 = df.Function(QV)

    def compute_C(C, space, u, family):
        if family == "DG" or family == "CG":
            C.assign(df.project((1.0 + u)**2, space))
        elif family == "Quadrature":
            local_project(q(u), space, C)
            assert False

    def compute_sigma(sigma, space, u, family):
        if family == "DG" or family == "CG":
            sigma.assign(df.project((1.0 + u)**2*df.grad(u), space))
        elif family == "Quadrature":
            local_project(q(u)*df.grad(u), space, sigma)
            assert False

    # initialize C and sigma0
    compute_C(C, Q, u0, args['--family'])
    compute_sigma(sigma0, QV, u0, args['--family'])

    a = df.inner(C*df.grad(du), df.grad(v))*dx + df.inner(Dq(u0)*du*df.grad(u0), df.grad(v))*dx
    L = - df.inner(sigma0, df.grad(v))*dx

    du = df.Function(V)
    u = df.Function(V)
    omega = 1.0
    eps = 1.0
    tol = 1e-5
    niter = 0
    maxiter = 25
    while eps > tol and niter < maxiter:
        niter += 1
        A, b = df.assemble_system(a, L, bcs_du)
        df.solve(A, du.vector(), b)
        eps = np.linalg.norm(du.vector()[:], ord=np.inf)
        print("Residual: ", eps)
        # solution and history update
        u.vector()[:] = u0.vector() + omega*du.vector()
        compute_C(C, Q, u0, args['--family'])
        compute_sigma(sigma0, QV, u0, args['--family'])

    if args['--summary']:
        convergence = f"convergence after {niter} newton iterations."
        if niter >= maxiter:
            convergence = "no " + convergence
        # Find max error
        u_exact = df.Expression('pow((pow(2, m+1)-1)*x[0] + 1, 1.0/(m+1)) - 1', m=args['M'],
        u_e = df.interpolate(u_exact, V)
        diff = np.abs(u_e.vector()[:] - u.vector()[:]).max()
        summary = f"""Nonlinear Poisson Results:
            Number of DOFs:     {V.dim()},
            Degree of V:        {args['ORDER']},
            family:             {args['--family']},
            convergence:        {convergence},
            nonlinearity:       (1 + u)**{args['M']},
            max error:          {diff}\n"""
    if args['--save']:
        print("Writing summary results to file.")
            with open("summary.out", "a") as outstream:
            with open("summary.out", "w") as outstream:

if __name__ == "__main__":