Hello everyone! Now the program can run smoothly, but there are obvious errors. For example, the L2 error order of the 8th and 16th order has a great deviation from the H1 error order. Does anyone know where this problem lies and how to solve it?
from __future__ import print_function
import matplotlib.pyplot as plt
from fenics import *
from ufl import nabla_div
from dolfin import *
from sympy import cos,sin
def NS(N):
T = 1.0
num_steps = pow(N, 2)
dt = T / num_steps
# Creat Mesh and function spaces 创造函数空间
mesh = UnitSquareMesh(N, N)
f_degree = 3
my_degree = 1
#普通网格
#V = VectorFunctionSpace(mesh, 'P', 2)
#Q = FunctionSpace(mesh, 'P', 1)
#W = FunctionSpace(mesh, V * Q)
# p2-p1
#V = VectorElement("Lagrange", mesh.ufl_cell(), my_degree + 1)
#Q = FiniteElement("CG", mesh.ufl_cell(), 1)
#W = FunctionSpace(mesh, V * Q )
# p1b-p1
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
B = FiniteElement("Bubble", mesh.ufl_cell(), 3)
V= VectorElement(P1 + B)
Q = P1
W = FunctionSpace(mesh, V * Q)
Vh1 = FunctionSpace(mesh, V)
Vh2 = FunctionSpace(mesh, Q)
# Define boundary conditions 定义边界条件
t=0
u1_ex = """10*pow(x[0],2)*pow((x[0]-1),2)*x[1]*(x[1]-1)*(2*x[1]-1)"""
u2_ex = """10*x[0]*(x[0]-1)*(2*x[0]-1)*pow(x[1],2)*pow((x[1]-1),2)"""
#u3_ex = """0.0"""
u_D = Expression(('10*pow(x[0],2)*pow((x[0]-1),2)*x[1]*(x[1]-1)*(2*x[1]-1)*cos(t)','10*x[0]*(x[0]-1)*(2*x[0]-1)*pow(x[1],2)*pow((x[1]-1),2)*cos(t)'), degree=f_degree, t=t)
# P条件
p1_ex = """10*(2*x[0]-1)*(2*x[1]-1)"""
p_D = Expression('10*(2*x[0]-1)*(2*x[1]-1)*cos(t)', degree=f_degree, t=t)
# Define expressions used in variational forms 定义用于变分形式的表达式
f1_ex = """20*cos(t)*(2*x[1] - 1) - 40*pow(x[0] ,2)*cos(t)*pow((x[0] - 1) ,2)*(x[1] - 1) - 20*pow(x[0] ,2)*cos(t)*(2*x[1] - 1)*pow((x[0] - 1),2) - 40*pow( x[0],2)*x[1]*cos(t)*pow((x[0] - 1) ,2) - 20*pow(x[0] ,2)*x[1]*cos(t)*(2*x[1] - 1)*(x[1] - 1) - 20*x[1]*cos(t)*(2*x[1] - 1)*pow((x[0] - 1) ,2)*(x[1] - 1) - 40*x[0]*x[1]*cos(t)*(2*x[0] - 2)*(2*x[1] - 1)*(x[1] - 1) - 10*pow(x[0] ,2)*x[1]*sin(t)*(2*x[1] - 1)*pow( (x[0] - 1),2)*(x[1] - 1) - 10*x[0]*pow( x[1],2)*cos(t)*(2*x[0] - 1)*(x[0] - 1)*pow((x[1] - 1) ,2)* (10*pow( x[0],2)*cos(t)*(2*x[1] - 1)*pow((x[0] - 1) ,2)*(x[1] - 1) + 20*pow(x[0] ,2)*x[1]*cos(t)*pow((x[0] - 1) ,2)*(x[1] - 1) + 10*pow(x[0] ,2)*x[1]*cos(t)*(2*x[1] - 1)*pow((x[0] - 1) ,2)) + 10*pow(x[0] ,2)*x[1]*cos(t)*(2*x[1] - 1)*(20*x[0]*x[1]*cos(t)*(2*x[1] - 1)*pow( (x[0] - 1),2)*(x[1] - 1) + 10*pow(x[0] ,2)*x[1]*cos(t)*(2*x[0] - 2)*(2*x[1] - 1)*(x[1] - 1))*pow((x[0] - 1) ,2)*(x[1] - 1)"""
# 20*cos(t)*(2*y - 1)- 40*x^2*cos(t)*(x - 1)^2*(y - 1)- 20*x^2*cos(t)*(2*y - 1)*(x - 1)^2- 40*x^2*y*cos(t)*(x - 1)^2 - 20*x^2*y*cos(t)*(2*y - 1)*(y - 1) - 20*y*cos(t)*(2*y - 1)*(x - 1)^2*(y - 1) - 40*x*y*cos(t)*(2*x - 2)*(2*y - 1)*(y - 1) - 10*x^2*y*sin(t)*(2*y - 1)*(x - 1)^2*(y - 1) - 10*x*y^2*cos(t)*(2*x - 1)*(x - 1)*(y - 1)^2*(10*x^2*cos(t)*(2*y - 1)*(x - 1)^2*(y - 1) + 20*x^2*y*cos(t)*(x - 1)^2*(y - 1) + 10*x^2*y*cos(t)*(2*y - 1)*(x - 1)^2) + 10*x^2*y*cos(t)*(2*y - 1)*(20*x*y*cos(t)*(2*y - 1)*(x - 1)^2*(y - 1) + 10*x^2*y*cos(t)*(2*x - 2)*(2*y - 1)*(y - 1))*(x - 1)^2*(y - 1)
f2_ex = """2*cos(t)*(20*x[0] - 10) + 40*pow(x[1] ,2)*cos(t)*(x[0] - 1)*pow((x[1] - 1) ,2) + 20*pow(x[1] ,2)*cos(t)*(2*x[0] - 1)*pow((x[1] - 1) ,2) + 40*x[0]*pow(x[1] ,2)*cos(t)*pow( (x[1] - 1),2) + 20*x[0]*pow(x[1] ,2)*cos(t)*(2*x[0] - 1)*(x[0] - 1) + 20*x[0]*cos(t)*(2*x[0] - 1)*(x[0] - 1)*pow((x[1] - 1) ,2) + 40*x[0]*x[1]*cos(t)*(2*x[0] - 1)*(2*x[1] - 2)*(x[0] - 1) + 10*x[0]*pow(x[1] ,2)*sin(t)*(2*x[0] - 1)*(x[0] - 1)*pow((x[1] - 1) ,2) - 10*pow(x[0] ,2)*x[1]*cos(t)*(2*x[1] - 1)*pow((x[0] - 1) ,2)*(x[1] - 1)*(10*pow(x[1] ,2)*cos(t)*(2*x[0] - 1)*(x[0] - 1)*pow((x[1] - 1) ,2) + 20*x[0]*pow(x[1] ,2)*cos(t)*(x[0] - 1)*pow( (x[1] - 1),2) + 10*x[0]*pow( x[1],2)*cos(t)*(2*x[0] - 1)*pow((x[1] - 1) ,2)) + 10*x[0]*pow(x[1] ,2)*cos(t)*(2*x[0] - 1)*(20*x[0]*x[1]*cos(t)*(2*x[0] - 1)*(x[0] - 1)*pow((x[1] - 1) ,2) + 10*x[0]*pow(x[1] ,2)*cos(t)*(2*x[0] - 1)*(2*x[1] - 2)*(x[0] - 1))*(x[0] - 1)*(x[1] - 1)*pow( (x[1] - 1),2)"""
# 2*cos(t)*(20*x - 10) + 40*y^2*cos(t)*(x - 1)*(y - 1)^2 + 20*y^2*cos(t)*(2*x - 1)*(y - 1)^2 + 40*x*y^2*cos(t)*(y - 1)^2 + 20*x*y^2*cos(t)*(2*x - 1)*(x - 1) + 20*x*cos(t)*(2*x - 1)*(x - 1)*(y - 1)^2 + 40*x*y*cos(t)*(2*x - 1)*(2*y - 2)*(x - 1) + 10*x*y^2*sin(t)*(2*x - 1)*(x - 1)*(y - 1)^2 - 10*x^2*y*cos(t)*(2*y - 1)*(x - 1)^2*(y - 1)*(10*y^2*cos(t)*(2*x - 1)*(x - 1)*(y - 1)^2 + 20*x*y^2*cos(t)*(x - 1)*(y - 1)^2 + 10*x*y^2*cos(t)*(2*x - 1)*(y - 1)^2) + 10*x*y^2*cos(t)*(2*x - 1)*(20*x*y*cos(t)*(2*x - 1)*(x - 1)*(y - 1)^2 + 10*x*y^2*cos(t)*(2*x - 1)*(2*y - 2)*(x - 1))*(x - 1)*(y - 1)^2
# f1_ex = """0.0"""
f = Expression((f1_ex, f2_ex), degree=f_degree,t=t)
#Dirichlet条件
#def boundary(x):
#tol = 1E-12
#return abs(x[0]) < tol or abs(x[1]) < tol or abs(x[0] - 1) < tol or abs(x[1] - 1) < tol
#bcu = DirichletBC(W.sub(0), u_D, boundary)
#bcu = [bc]
bc = DirichletBC(W.sub(0), u_D, "on_boundary")
bcu = [bc]
# Define functions for solutions at previous and current time stepsfrom _future_import print function
#从_future_import打印函数中定义解决方案在过去和当前时间步骤的函数
# Define trial and test functions 试探和测试函数
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
u_n = interpolate(u_D, Vh1)
u_ = Function(Vh1)
p_ = Function(Vh2)
#U = Function(W)
W_ex = Expression((u1_ex, u2_ex, p1_ex), degree=f_degree)
u_s = interpolate(W_ex, W)
Uk = Function(W)
# Time-stepping 时域
t=0
for N in range(num_steps):
# Update current time 更新当前时间
#初始值变分块
#u1_ex = """10*pow(x[0],2)*pow((x[0]-1),2)*x[1]*(x[1]-1)*(2*x[1]-1)*cos(t)"""
#u2_ex = """10*x[0]*(x[0]-1)*(2*x[0]-1)*pow(x[1],2)*pow((x[1]-1),2)*cos(t)"""
u_D = Expression(('10*pow(x[0],2)*pow((x[0]-1),2)*x[1]*(x[1]-1)*(2*x[1]-1)*cos(t)',\
'10*x[0]*(x[0]-1)*(2*x[0]-1)*pow(x[1],2)*pow((x[1]-1),2)*cos(t)'), degree=f_degree, t=t)
# P条件
#p1_ex = """10*(2*x[0]-1)*(2*x[1]-1)*cos(t)"""
p_D = Expression('10*(2*x[0]-1)*(2*x[1]-1)*cos(t)', degree=f_degree, t=t)
k = Constant(dt)
F=(dot((u - u_n), v) + k*inner(grad(u),grad(v)) + k*inner(dot(u_n,grad(u)),v) - div(u)*q \
+ k*div(v)*p + 10e-8*inner(p,q))*dx - k*inner(f,v)*dx
#F = inner((u - u_n) / k, v) * dx - inner(f, v) * dx + inner(grad(u), grad(v)) * dx - inner(grad(p), v) * dx \
#+ inner(grad(u_n) * u, v) * dx \
#+ inner(grad(q), u) * dx
a = lhs(F)
L = rhs(F)
u_D.t = t
f.t = t
t += dt
A = assemble(a)
[bc.apply(A) for bc in bcu]
b = assemble(L)
[bc.apply(b) for bc in bcu]
#解法一
#solve(A == b, u_s, bcu)
#(u_ , p_)=u_s.split()
#解法二
solve(A, u_s.vector(), b)
#u_n.assign(u_)
(u_, p_) = u_s.split()
#误差计算
u_error_L2 = errornorm(u_D, u_, 'L2')
u_error_H1 = errornorm(u_D, u_, 'H1')
p_error_L2 = errornorm(p_D, p_, 'L2')
p_error_H1 = errornorm(p_D, p_, 'H1')
return [u_error_L2 , u_error_H1, p_error_L2, p_error_H1]
#for N in [2, 4, 8, 16]:
out = []
for N in [2,4,8]:
t = NS(N)
out = out + t
print(out[0:4]) # 2 [u_error_L2 , u_error_H1, p_error_L2, p_error_H1]
print(out[4:8]) # 4 [u_error_L2 , u_error_H1, p_error_L2, p_error_H1]
print(out[8:12]) # 8 [u_error_L2 , u_error_H1, p_error_L2, p_error_H1]
Here is results
N=2,4,8
[0.04048438418096054, 0.3722916114590978, 4.097405149150013, 21.0971489140043]
[0.026808636805903896, 0.21622793229724555, 3.8975503602723465, 19.575437769175945]
[0.023356153171252336, 0.17572340833172742, 3.685405072757736, 18.462012322798884]
thanks everyone