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.
Best,
Philipp
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""Solve the nonlinear poisson eq.
- nabla · (q(u) nabla u) = f
Usage:
nonlinear_poisson.py [options] NELEMENTS M ORDER
nonlinear_poisson.py -h | --help
Arguments:
NELEMENTS Number of elements in each spatial direction.
M q(u) = (1 + u)^M
ORDER Order of Lagrange polynomial.
Options:
--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
WARNING = 30
INFO = 20
PROGRESS = 16
df.set_log_level(WARNING)
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)
solver.factorize()
if u is None:
u = df.Function(V)
solver.solve_local_rhs(u)
df.parameters["form_compiler"]["representation"] = "uflacs"
return u
else:
solver.solve_local_rhs(u)
df.parameters["form_compiler"]["representation"] = "uflacs"
return
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,
quad_scheme="default")
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)
else:
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)
else:
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()
u0.assign(u)
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'],
element=V.ufl_element())
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"""
print(summary)
if args['--save']:
print("Writing summary results to file.")
try:
with open("summary.out", "a") as outstream:
outstream.write(summary)
except:
with open("summary.out", "w") as outstream:
outstream.write(summary)
if __name__ == "__main__":
main(sys.argv[1:])