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 https://fenicsproject.org/pub/tutorial/html/._ftut1020.html
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)
    a = kappa*dot(grad(u), grad(v))*dx
    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

1 Answer


5
8 months ago by
Nate  
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  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »