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 

    a = dot(div(grad(u)), div(grad(v)))*dx \
       + inner(jump(v), avg(dot(grad(div(grad(u))), n)))*dS \
       + lamb1*inner(jump(u), avg(dot(grad(div(grad(v))), n)))*dS \
       - inner(avg(div(grad(u))), jump(grad(v),n))*dS \
       - lamb2*inner(jump(grad(u),n), avg(div(grad(v))))*dS \
       + inner(alpha_1*jump(u), jump(v))*dS \
       + inner(beta_1*jump(grad(u), n),jump(grad(v), n))*dS \
       + inner(v, dot(grad(div(grad(u))),n))*ds \
       + lamb1*inner(u, dot(grad(div(grad(v))),n))*ds \
       + inner(alpha_1*u, v)*ds \
       - inner(div(grad(u)), dot(grad(v),n))*ds \
       - lamb2*inner(dot(grad(u),n), div(grad(v)))*ds \
       + inner(beta_1*dot(grad(u), n),dot(grad(v), n))*ds

       

# Define linear form(1) 
    L = f*v*dx \
       + lamb1*inner(g_0, dot(grad(div(grad(v))),n))*ds \
       - lamb2*inner(g_1, div(grad(v)))*ds \
       + alpha_1*g_0*v*ds \
       + inner(beta_1*g_1, dot(grad(v),n))*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
Please login to add an answer/comment or follow this question.

Similar posts:
Search »