### Why the error don't converge when refine as nx=512?

448
views
0
8 months ago by
Hi,

I am testing the code ft10_poisson_extended.py of this page And I only change maximum refinement level as nx=ny=512, and run the code as below

from __future__ import print_function

from fenics import *
import numpy as np

def solver(kappa, f, u_D, Nx, Ny,
degree=1,
linear_solver='Krylov',
abs_tol=1E-5,
rel_tol=1E-3,
max_iter=1000):
"""
Solve -div(kappa*grad(u) = f on (0, 1) x (0, 1) with 2*Nx*Ny Lagrange
elements of specified degree and u = u_D on the boundary.
"""

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

# Define boundary condition
def boundary(x, on_boundary):
return on_boundary

bc = DirichletBC(V, u_D, boundary)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
L = f*v*dx

# 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, bc, 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, u_D, kappa,
max_degree=3, 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(1, 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, u_D, n, n, degree, linear_solver='direct')
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[1][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
u_e = Expression('sin(omega*pi*x[0])*sin(omega*pi*x[1])',
degree=6, omega=omega)
f = 2*omega**2*pi**2*u_e
u_D = Constant(0)
kappa = Constant(1)

# Compute and print convergence rates
etypes, degrees, rates = compute_convergence_rates(u_e, f, u_D, 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()​

The results are:

Solving linear variational problem.
2 x (256 x 256) P3 mesh, 591361 unknowns, E5 = 1.87663e-11
Solving linear variational problem.
2 x (512 x 512) P3 mesh, 2362369 unknowns, E5 = 2.01626e-11

L2 norm
P1: 1.97, 1.99, 2.0, 2.0, 2.0, 2.0
P2: 3.0, 3.0, 3.0, 3.0, 3.0, 2.99
P3: 4.04, 4.02, 4.01, 4.0, 3.95, -0.1

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.
Community: FEniCS Project

5
8 months ago by
You've hit machine precision.
Dear Nate,

Thank you very much for replying so quick!
But the double precision is approximately 1.11e-16, why have I hit machine precision when the result is just approximately 1.87663e-11?
Thank you!

Chucui1016
written 8 months ago by chucui1016
You need to bear in mind you're computing something like $\sqrt{\epsilon}$ϵ
written 8 months ago by Nate
OK, thank you very much!
written 8 months ago by chucui1016