### 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