after time loop, the solution can't hit exact solution


182
views
1
6 months ago by
Dear all,
I want to implement this :
To make an exact solution at first:

for problem with pure neumann boundary conditions:

where


Time discretisation is :

rewrite it to solve

and parameters are:

and initional solution is 

and other parameters are constants, you can see them in the code


Here is my code:
from __future__ import print_function

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


# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["optimize"] = True

# Make mesh ghosted for evaluation of DG terms
parameters["ghost_mode"] = "shared_facet"


def solver(Nx, Ny,
           degree,
           dt,
           s,
           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 Lagrange
    elements of specified degree and u = u_D on the boundary.
    """
    #dt     = 0.01  # time step
    t  = 0*dt
    tt = t + 0.5*dt
    l1 = 2.0
    l2 = 2.0  
    TT = t + s*dt
    # Create mesh and define function space
    mesh = RectangleMesh(Point(-1*pi, -1*pi), Point(1*pi, 1*pi), Nx, Ny)
    V = FunctionSpace(mesh, 'P', degree)


    # Create mesh and define function space
    #W = FunctionSpace(mesh,"Lagrange",degree)
    R  = FunctionSpace(mesh,"R",0)
    W =  V * R

    # Define trial and test functions
    (u, c)  = TrialFunction(W)
    (v, d)  = TestFunction(W)

    # Define functions
    #u   = Function(W)  # current solution
    w = Function(W)
    w0  = Function(W)  # solution from previous converged step
    w1  = Function(W)

    # Split mixed functions
    #u, c  = split(w)
    (u0, c0)  = split(w0)
    (u1, c1)  = split(w1)

    # Initional conditions
    class InitialConditions(Expression):
        #def __init__(self, **kwargs):
        #   random.seed(2 + MPI.rank(mpi_comm_world()))
        def eval(self, values, x):
            values[0] = cos(l1*x[0])*cos(l2*x[1])*(t*t+0.01)

    u_init = InitialConditions(degree=2)

    u0 = interpolate(u_init, V)
    u1 = u0 

# Some parameters
# which is then used in the definition of the variational forms::
    class Source(Expression):
        def eval(self, values, x):
            values[0] = cos(l1*x[0])*cos(l2*x[1])
    M = Constant(1.0)
    alpha = Constant(1.0)        
    K = Constant(2.0)      
    g = Constant(0.0)        
    z = Source(degree=2)
    fplus1 = 2*tt*z + M*(K*(l1*l1+l2*l2)*z*(tt*tt+0.01)+alpha*z*(tt*tt+0.01)+g)


    g_1 = Constant(0.0)  #Neumann boundary condition
    #g_2 = Constant(0.0)
    A = -K*(M*dt)/2
    B = (alpha*M*dt)/2+1
    

# varitional forms

    a = -A*dot(grad(v),grad(u))*dx + B*u*v*dx + c*v*dx +u*d*dx

    L = -A*v*g_1*ds - A*dot(v,div(grad(u0)))*dx - (B-2)*v*u0*dx - M*dt*g*v*dx + dt*fplus1*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
    w   = Function(W)
    solve(a == L, w, solver_parameters={'linear_solver':'mumps'})
    (u, c) = split(w)
#Output file
    file = File("output.pvd", "compressed")
    
# Step in time
    T = TT 
    for i in range(s-1):
        t += dt
        tt += dt
        #w1.vector()[:] = w0.vector()
        #w0.vector()[:] = w.vector()
        u_last = u
        (u, c)  = TrialFunction(W)
        (v, d)  = TestFunction(W)

    # Define functions
    #u   = Function(W)  # current solution
        w = Function(W)
        w0  = Function(W)  # solution from previous converged step
        #w1  = Function(W)

    # Split mixed functions
    #u, c  = split(w)
        (u0, c0)  = split(w0)
        #(u1, c1)  = split(w1)
    #c0, mu0 = split(u0)
        
        
        #u1 = u0
        u0 = u_last
        fplus1 = 2*tt*z + M*(K*(l1*l1+l2*l2)*z*(tt*tt+0.1)+alpha*z*(tt*tt+0.1)+g)

        
        a = -A*dot(grad(v),grad(u))*dx + B*u*v*dx + c*v*dx +u*d*dx
        L = -A*v*g_1*ds - A*dot(v,div(grad(u0)))*dx - (B-2)*v*u0*dx - M*dt*g*v*dx + dt*fplus1*v*dx
        w   = Function(W)
        solve(a == L, w, solver_parameters={'linear_solver':'mumps'})
        (u, c) = split(w)
        u = w.split()[0]
        file << u
        #plot (u)

    return u

def compute_errors(u_e, u, n, degree):
    """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
    mesh = RectangleMesh(Point(-1*pi, -1*pi), Point(1*pi, 1*pi), n, n)    
    V = FunctionSpace(mesh,"Lagrange",degree)
    R  = FunctionSpace(mesh,"R",0)
    W =  V * R
    #W = u.function_space()
    u_0 = u#.split()[0]
    u_e0 = u_e#.split()[0]
    V = FunctionSpace(mesh,"Lagrange",degree) 
    #V = u_0.function_space()
    # Explicit computation of L2 norm
    error = (u_0 - u_e0)**2*dx(mesh)
    E1 = sqrt(abs(assemble(error)))

    # Explicit interpolation of u_e onto the same space as u
    u_e0_ = interpolate(u_e0, V)
    error = (u_0 - u_e0_)**2*dx(mesh)
    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', degree+6)
    u_e0_ = interpolate(u_e0, Ve)
    error = (u_0 - u_e0)**2*dx(mesh)
    E3 = sqrt(abs(assemble(error)))

    # Infinity norm based on nodal values
    u_e0_ = interpolate(u_e0, V)
    #E4 = abs(u_e0_.vector().array() - u_0.vector().array()).max()
    E4 = E3
    # L2 norm
    E5 = errornorm(u_e0_, u_0, norm_type='L2', degree_rise=3)

    # H1 seminorm
    E6 = errornorm(u_e0, u_0, 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(dt, s,
                              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(2, max_degree + 1)
    for degree in degrees:
        n = 8  # coarsest mesh division
        mesh = RectangleMesh(Point(-1*pi, -1*pi), Point(1*pi, 1*pi), n, n)    
        V = FunctionSpace(mesh,"Lagrange",degree)
        R  = FunctionSpace(mesh,"R",0)
        W =  V * R
        u_ex = Function(V)
        TT = s*dt
        l1 = Constant(2.0)
        l2 = Constant(2.0)
        class exactsolution(Expression):
    #def __init__(self, **kwargs):
     #   random.seed(2 + MPI.rank(mpi_comm_world()))
            def eval(self, values, x):
                values[0] = cos(l1*x[0])*cos(l2*x[1])*(TT*TT+0.01)
              
        
            #def value_shape(self):
                #return (2,)
        uex = exactsolution(degree=degree+6)
#u.interpolate(u_init)
        u_ex.interpolate(uex)
        file2 = File("periodic_exact.pvd", "compressed")
#u_exact = interpolate(uex, VV)
        file2 << u_ex
        #plot (u_ex)
        h[degree] = []
        E[degree] = []
        for i in range(num_levels):
            h[degree].append(1.0 / n)
            u = solver(n, n, degree, dt, s, linear_solver='Krylov', abs_tol=1E-16, rel_tol=1E-13, max_iter=10000000)
            plot (u)
            #u_ex0 = u_ex.split()[0]            
            #u_0 = u.split()[0] 
            errors = compute_errors(u_ex, u, n, degree)
            E[degree].append(errors)
            print('2 x (%d x %d) P%d mesh, %d unknowns, E3 = %f' %
              (n, n, degree, u.function_space().dim(), errors['interpolate(u, Ve) - interpolate(u_e, Ve)']))
            n *= 2
            mesh = RectangleMesh(Point(-1*pi, -1*pi), Point(1*pi, 1*pi), n, n)    
            V = FunctionSpace(mesh,"Lagrange",degree)
            R  = FunctionSpace(mesh,"R",0)
            W =  V * R
            u_ex = Function(V)
            TT = s*dt
            class exactsolution(Expression):
    #def __init__(self, **kwargs):
     #   random.seed(2 + MPI.rank(mpi_comm_world()))
                def eval(self, values, x):
                    values[0] = cos(l1*x[0])*cos(l2*x[1])*(TT*TT+0.01)
                    #values[1] = -((2*2/pi)*(2*2/pi)+(2*3/pi)*(2*3/pi))*cos(2*2*x[0]/pi)*cos(2*3*x[1]/pi)*TT*TT #n1=2,n2=3
        
                #def value_shape(self):
                    #return (2,)
            uex = exactsolution(degree=8)
#u.interpolate(u_init)
            u_ex.interpolate(uex)
            file2 = File("periodic_exact.pvd", "compressed")
#u_exact = interpolate(uex, VV)
            file2 << u_ex
           # print (TT)
            plot (u_ex)
    # 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"
    dt = 0.01
    s = 10
   
    # Compute and print convergence rates
    etypes, degrees, rates = compute_convergence_rates(dt, s)
    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()​

And I make a plot for numerical solution and exact solution, but they are different, and I don't find the mistake. Please give me some help. Thank you!
Community: FEniCS Project
Please login to add an answer/comment or follow this question.

Similar posts:
Search »