### Why the error of solving biharmonic equation by nonsymmetric interior penalty DGFEM (NIPG) cannot converge at n=512 or finer?

186
views
0
6 months ago by
I have a similar question as https://www.allanswered.com/post/pgqbp/why-the-error-dont-converge-when-refine-as-nx512/  but the reason may not be hitting the machine precision.
The problem is to solve:

Now I need to solve it by discontinuous galerkin finite element method (DGFEM), and the paper is here:
File attached: hp-version interior penalty DGFEMs for the biharmonic equation.pdf (370.78 KB)

The bilinear form is:

where
and the linear form is:

The exact solution is:

so,

and the parameters in bilinear form and linear form are chosen from the paper:

Here are my codes for p=2,
λ1=-1,λ2=-1, which is called nonsymmetric interior penalty DGFEM (NIPG) in the paper

from __future__ import print_function

from fenics import *
#from boxfield import *
import numpy as np

def solver(kappa, f, g_0, g_1, Nx, Ny,
degree=1,
linear_solver='Krylov',
abs_tol=1E-16,
rel_tol=1E-13,
max_iter=1000000):
"""
Solve div(kappa*grad(u) = f on (0, 1) x (0, 1) with 2*Nx*Ny DG
elements of specified degree and u = u_D on the boundary.
"""
# Define class marking Dirichlet boundary (x = 0 or x = 1)
class DirichletBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary #and near(x[0]*(1 - x[0]), 0)

# Define class marking Neumann boundary (y = 0 or y = 1)
class NeumanBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary #and near(x[1]*(1 - x[1]), 0)

# Create mesh and define function space
mesh = UnitSquareMesh(Nx, Ny)
V = FunctionSpace(mesh, 'DG', degree)

# Define normal vector and mesh size
n = FacetNormal(mesh)
h = CellSize(mesh)
h_avg = (h('+') + h('-'))/2

# Mark facets of the mesh
boundaries = FacetFunction('size_t', mesh, 0)
#NeumanBoundary().mark(boundaries)
#DirichletBoundary().mark(boundaries)
NeumanBoundary().mark(boundaries, 2)
DirichletBoundary().mark(boundaries, 1)

u0 = Constant(0.0)
bc = DirichletBC(V, u0, DirichletBoundary())
#bc = DirichletBC(V, Constant(0.0), DirichletBoundary)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)

# Define outer surface measure aware of Dirichlet and Neumann boundaries
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

# Define parameters
alpha = 4.0
gamma = 8.0
lamb1 = -1.0
lamb2 = -1.0

alpha_1 = 10*degree*degree*degree*degree*Nx*Nx*Nx
beta_1 = 10*degree*degree*Nx

+ inner(alpha_1*jump(u), jump(v))*dS \
+ inner(alpha_1*u, v)*ds \

# Define linear form(1)
L = f*v*dx \
+ alpha_1*g_0*v*ds \

# Set linear solver parameters
prm = LinearVariationalSolver.default_parameters()
if linear_solver == 'Krylov':
prm.linear_solver = 'gmres'
prm.preconditioner = 'ilu'
prm.krylov_solver.absolute_tolerance = abs_tol
prm.krylov_solver.relative_tolerance = rel_tol
prm.krylov_solver.maximum_iterations = max_iter
else:
prm.linear_solver = 'lu'

# Compute solution
u = Function(V)
#solve(a == L, u, bc, solver_parameters=prm)
solve(a == L, u, solver_parameters={'linear_solver':'mumps'})
return u

def compute_errors(u_e, u):
"""Compute various measures of the error u - u_e, where
u is a finite element Function and u_e is an Expression."""

# Get function space
V = u.function_space()

# Explicit computation of L2 norm
error = (u - u_e)**2*dx
E1 = sqrt(abs(assemble(error)))

# Explicit interpolation of u_e onto the same space as u
u_e_ = interpolate(u_e, V)
error = (u - u_e_)**2*dx
E2 = sqrt(abs(assemble(error)))

# Explicit interpolation of u_e to higher-order elements.
# u will also be interpolated to the space Ve before integration
Ve = FunctionSpace(V.mesh(), 'P', 5)
u_e_ = interpolate(u_e, Ve)
error = (u - u_e)**2*dx
E3 = sqrt(abs(assemble(error)))

# Infinity norm based on nodal values
u_e_ = interpolate(u_e, V)
E4 = abs(u_e_.vector().array() - u.vector().array()).max()

# L2 norm
E5 = errornorm(u_e, u, norm_type='L2', degree_rise=3)

# H1 seminorm
E6 = errornorm(u_e, u, norm_type='H10', degree_rise=3)

# Collect error measures in a dictionary with self-explanatory keys
errors = {'u - u_e': E1,
'u - interpolate(u_e, V)': E2,
'interpolate(u, Ve) - interpolate(u_e, Ve)': E3,
'infinity norm (of dofs)': E4,
'L2 norm': E5,
'H10 seminorm': E6}

return errors

def compute_convergence_rates(u_e, f, g_0, g_1, kappa,
max_degree=2, num_levels=7):
"Compute convergences rates for various error norms"

h = {}  # discretization parameter: h[degree][level]
E = {}  # error measure(s): E[degree][level][error_type]

# Iterate over degrees and mesh refinement levels
degrees = range(2, max_degree + 1)
for degree in degrees:
n = 8  # coarsest mesh division
h[degree] = []
E[degree] = []
for i in range(num_levels):
h[degree].append(1.0 / n)
u = solver(kappa, f, g_0, g_1, n, n, degree, linear_solver='Krylov', abs_tol=1E-16, rel_tol=1E-13, max_iter=100000000)
errors = compute_errors(u_e, u)
E[degree].append(errors)
print('2 x (%d x %d) P%d mesh, %d unknowns, E5 = %g' %
(n, n, degree, u.function_space().dim(), errors['L2 norm']))
n *= 2

# Compute convergence rates
from math import log as ln  # log is a fenics name too
etypes = list(E[2][0].keys())
rates = {}
for degree in degrees:
rates[degree] = {}
for error_type in sorted(etypes):
rates[degree][error_type] = []
for i in range(1, num_levels):
Ei = E[degree][i][error_type]
Eim1 = E[degree][i - 1][error_type]
r = ln(Ei / Eim1) / ln(h[degree][i] / h[degree][i - 1])
rates[degree][error_type].append(round(r, 2))

return etypes, degrees, rates

def demo_convergence_rates():
"Compute convergence rates in various norms for P1, P2, P3"

# Define exact solution and coefficients
omega = 1.0

class Source(Expression):
def eval(self, values, x):
values[0] = sin(omega*pi*x[0])*sin(omega*pi*x[1])*sin(omega*pi*x[0])*sin(omega*pi*x[1])
class Source1(Expression):
def eval(self, values, x):
values[0] = 24*pi**4*sin(pi*x[1])*sin(pi*x[0])*sin(pi*x[1])*sin(pi*x[0])-16*pi**4*cos(pi*x[0])*sin(pi*x[1])*cos(pi*x[0])*sin(pi*x[1])+8*pi**4*cos(pi*x[0])*cos(pi*x[1])*cos(pi*x[0])*cos(pi*x[1])-16*pi**4*cos(pi*x[1])*sin(pi*x[0])*cos(pi*x[1])*sin(pi*x[0])
u_e = Source(degree=3+3)
f = Source1(degree=3+3)

g_0 = Constant(0)
g_1 = Constant(0)
kappa = Constant(1)

# Compute and print convergence rates
etypes, degrees, rates = compute_convergence_rates(u_e, f, g_0, g_1, kappa)
for error_type in etypes:
print('\n' + error_type)
for degree in degrees:
print('P%d: %s' % (degree, str(rates[degree][error_type])[1:-1]))

if __name__ == '__main__':

# Run demo
demo_convergence_rates()

# Hold plot
interactive()

But the results are like this
NIPG:
2 x (8 x 8) P2 mesh, 768 unknowns, E5 = 0.0627282
2 x (16 x 16) P2 mesh, 3072 unknowns, E5 = 0.0230552
2 x (32 x 32) P2 mesh, 12288 unknowns, E5 = 0.00670823
2 x (64 x 64) P2 mesh, 49152 unknowns, E5 = 0.00176518
2 x (128 x 128) P2 mesh, 196608 unknowns, E5 = 0.000449102
2 x (256 x 256) P2 mesh, 786432 unknowns, E5 = 0.000114316
2 x (512 x 512) P2 mesh, 3145728 unknowns, E5 = 5.04492e-05

and the convergence rate of L2 norm is:
L2 norm
P2: 1.44, 1.78, 1.93, 1.97, 1.97, 1.18

Why this kind of error happen? And I have interpolate u_e to a higher-order space, but the error still don't converge at nx=512 or finer. And this error is only  5.04492e-5, it can't hit the machine precision?

Community: FEniCS Project