Everybody, this is a program that I’ve been working on for a long time, and since I’m new to Fenics, Mr.Dokken has helped me run it. Now I want to ask for your help. The program out of L2, H1 errors are very large, error order is also wrong. I guess the problem may be a circular problem, but it still outputs errors after correcting according to the textbook. Here is the code:
import math
from dolfin import *
from fenics import dx
#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
def boundary(x, on_boundary): # define the Dirichlet boundary
return on_boundary
def E_L2_H1(n):
#定义网格
T = 1.0
num_steps = pow(n, 2)
dt = T / num_steps
mesh=UnitSquareMesh(n,n)
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)
#初始时刻条件
#U条件 Define manufactured exact solutions : OK 定义已有的精确解
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)"""
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=3, t=0)
# 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=3, t=0)
# force term of momentum equation 动量方程的力项
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)"""
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)"""
f = Expression((f1_ex, f2_ex), degree=3, t=0)
#Dirichlet条件
bc = DirichletBC(W.sub(0), u_D , boundary)
W_ex=Expression((u1_ex,u2_ex,p1_ex),degree=3)
Uw=interpolate(W_ex,W)
#Uw= Function(W)
#U=Function(W)
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
(uk, pk) = Uw.split()
#变分
k = Constant(dt)
f = (dot((u - uk), v) + k* inner(grad(u), grad(v)) + k * inner(dot(uk, grad(u) ), v) - div(u) * q + k * div(v) * p + 10e-8 * inner(p, q)) * dx - k * inner(f, v) * dx
a = lhs(f)
L = rhs(f)
#U.assign(Uw)
A = assemble(a)
#时间循环
U = Function(W)
t = 0
for n in range(num_steps):
#while (t <= T):
t += dt
u_D.t = t
b = assemble( L )
bc.apply(A, b)
#U = Function(W)
solve(A, U.vector(), b)
#solve(A == b , U , bc)
(u, p) = U.split()
Uw.assign(U)
#误差计算
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,num_steps]
def EOrder(e_L2_f,h_f,e_L2_a,h_a):
error_order=(math.log(e_L2_f,2)-math.log(e_L2_a,2))/(math.log(h_f,2)-math.log(h_a,2))
return error_order
n1=8
n2=16
[a1,b1,c1,d1,e1]=E_L2_H1(n1)
[a2,b2,c2,d2,e2]=E_L2_H1(n2)
print('---')
print('%d 等分u_error_L2=%.15g'%(n1,a1))
print('%d 等分u_error_H1=%.15g'%(n1,c1))
print('---')
print('%d 等分p_error_L2=%.15g'%(n1,b1))
print('%d 等分p_error_H1=%.15g'%(n1,d1))
print('---')
print('%d 等分u_error_L2=%.15g'%(n2,a2))
print('%d 等分u_error_H1=%.15g'%(n2,c2))
print('---')
print('%d 等分p_error_L2=%.15g'%(n2,b2))
print('%d 等分p_error_H1=%.15g'%(n2,d2))
print('---')
print('u L2error_order：',EOrder(a1,e1,a2,e2))
print('u H1error_order：',EOrder(c1,e1,c2,e2))
print('---')
print('p L2error_order：',EOrder(b1,e1,b2,e2))
print('p H1error_order：',EOrder(d1,e1,d2,e2))
If you need additional information, please leave a message, it is greatly appreciated!