Time-dependent Navier–Stokes equations

[Kan2046]

Hello everyone, here is a student from China who is trying to use Fenics to solve the timeNS problem, but the program still cannot run after debugging for a long time. Can anyone find out what the problem is?

from future import print_function
import matplotlib.pyplot as plt
from fenics import *
import sympy as sym
import numpy as np
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

t=0

# 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, MixedFunctionSpace([V , 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 trial and test functions 试探和测试函数
#u = TrialFunction(V)
#v = TestFunction(V)
#pt = TrialFunction(Q)
#qt = TestFunction(Q)

# Define functions for solutions at previous and current time stepsfrom _future_import print function
#从_future_import打印函数中定义解决方案在过去和当前时间步骤的函数

(u, p) = TrialFunction(V)
(v, q) = TestFunction(Q)
u_n = interpolate(u_D, V)
#u_s=Function(V)
#u_ = Function(V)
#p_n = Function(Q)
#p_ = Function(Q)
#u_s = Function(V)
u_ = Function(V)
p_ = Function(Q)

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

for N in range(num_steps):
    # Update current time 更新当前时间
    t=0
    #初始值变分块
    #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)

    #Uk.assign(u_s)
    #(u, p) = TrialFunction(W)
    #(v, q) = TestFunction(W)
    #(u_n, pk) = Uk.split()

    k = Constant(dt)

    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)


    #t += dt
    u_D.t = t
    f.t = t

    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_.vector(), b)
    #u_n.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')

t += dt

return u_error_L2,u_error_H1,p_error_L2,p_error_H1

for N in [2, 4, 8, 16]:
print("{0:d} u_error_L2 ={1:.15g}".format(N, NS(N)))
print("{0:d} u_error_H1 ={1:.15g}".format(N, NS(N)))

print("{0:d} p_error_L2 ={1:.15g}".format(N, NS(N)))
print("{0:d} p_error_H1 ={1:.15g}".format(N, NS(N)))

[Kan2046]

AND this is the code error

Traceback (most recent call last):
File “/mnt/c/Users/ALSACE/PycharmProjects/pythonProject3/main.py”, line 156, in
print("{0:d} u_error_L2 ={1:.15g}".format(N, NS(N)))
File “/mnt/c/Users/ALSACE/PycharmProjects/pythonProject3/main.py”, line 80, in NS
(u, p) = TrialFunction(V)
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/function/argument.py”, line 100, in TrialFunction
return Argument(V, 1, part)
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/function/argument.py”, line 53, in init
raise TypeError(_ufl_dolfin_difference_message)
TypeError: \ When constructing an Argument, TestFunction or TrialFunction, you
must to provide a FunctionSpace and not a FiniteElement. The
FiniteElement class provided by ufl only represents an abstract finite
element space and is only used in standalone .ufl files, while the
FunctionSpace provides a full discrete function space over a given
mesh and should be used in dolfin programs in Python.

[dokken]

Please format your code with 3x` encapsulation such that your code is rendered as

```
def test(x):
    return 2*x[0]

# This is a comment
from dolfin import *
```

[Kan2046]

I re-uploaded the code so that it can be rendered in its entirety

from __future__ import print_function
import matplotlib.pyplot as plt
from fenics import *
import sympy as sym
import numpy as np
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

    t=0

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

    # 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 trial and test functions 试探和测试函数
    #u = TrialFunction(V)
    #v = TestFunction(V)
    #pt = TrialFunction(Q)
    #qt = TestFunction(Q)

    # Define functions for solutions at previous and current time stepsfrom _future_import print function
    #从_future_import打印函数中定义解决方案在过去和当前时间步骤的函数

    (u, p) = TrialFunction(V)
    (v, q) = TestFunction(Q)
    u_n = interpolate(u_D, V)
    #u_s=Function(V)
    #u_ = Function(V)
    #p_n = Function(Q)
    #p_ = Function(Q)
    #u_s = Function(V)
    u_ = Function(V)
    p_ = Function(Q)

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

    for N in range(num_steps):
        # Update current time 更新当前时间
        t=0
        #初始值变分块
        #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)

        #Uk.assign(u_s)
        #(u, p) = TrialFunction(W)
        #(v, q) = TestFunction(W)
        #(u_n, pk) = Uk.split()

        k = Constant(dt)

        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)


        #t += dt
        u_D.t = t
        f.t = t

        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_.vector(), b)
        #u_n.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')

    t += dt

    return u_error_L2,u_error_H1,p_error_L2,p_error_H1


for N in [2, 4, 8, 16]:
    print("{0:d} u_error_L2 ={1:.15g}".format(N, NS(N)))
    print("{0:d} u_error_H1 ={1:.15g}".format(N, NS(N)))

    print("{0:d} p_error_L2 ={1:.15g}".format(N, NS(N)))
    print("{0:d} p_error_H1 ={1:.15g}".format(N, NS(N)))

[Kan2046]

Thanks Mr.dokken suggest !!!

[dokken]

The first issue is:

@Kan2046

(u, p) = TrialFunction(V)
    (v, q) = TestFunction(Q)

which should be

    (u, p) = TrialFunctions(W)
    (v, q) = TestFunctions(W)

Secondly, you get into a similar issue when trying to interpolate onto a finite element, not a function space, in:

@Kan2046

u_n = interpolate(u_D, V)

[Kan2046]

Thanks Mr.dokken. But in first issue, I try change V,Q to W.

File "/mnt/c/Users/ALSACE/PycharmProjects/pythonProject3/main.py", line 80, in NS
    (u, p) = TrialFunction(W)
ValueError: too many values to unpack (expected 2)

It make confuses to me.

The second issue,The way I understand it, the u sub n computed in the variational is of this form. The space V is combined with the function u_D to form a computable function. Do you mean it is not acceptable to use it in calculation? If there are other mistakes, please kindly advise.

[dokken]

@Kan2046

TrialFunction

As I said above, it should be Trialfunctions, not TrialFunction

@dokken

(u, p) = TrialFunctions(W)

@Kan2046

The second issue,The way I understand it, the u sub n computed in the variational is of this form. The space V is combined with the function u_D to form a computable function. Do you mean it is not acceptable to use it in calculation? If there are other mistakes, please kindly advise.

The second issue is that V is not a FunctionSpace, but a FiniteElement. You would have to create a FunctionSpace based on V if you want to do interpolation of any data.

[Kan2046]

Dear Mr.dokken
Thanks for yousuggest about FunctionSpace, and I already solve other similar issue.
But after the new code runs, there is a new error. The new code

import matplotlib.pyplot as plt
from fenics import *
import sympy as sym
import numpy as np
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

    t=0

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

    for N in range(num_steps):
        # Update current time 更新当前时间
        t=0
        #初始值变分块
        #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 = 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)

        #t += dt
        u_D.t = t
        f.t = t

        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_.vector(), b)
        #u_n.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')

    t += dt

    return u_error_L2,u_error_H1,p_error_L2,p_error_H1


for N in [2, 4, 8, 16]:
    print("{0:d} u_error_L2 ={1:.15g}".format(N, NS(N)))
    print("{0:d} u_error_H1 ={1:.15g}".format(N, NS(N)))

    print("{0:d} p_error_L2 ={1:.15g}".format(N, NS(N)))
    print("{0:d} p_error_H1 ={1:.15g}".format(N, NS(N)))

New error

  File "/mnt/c/Users/ALSACE/PycharmProjects/pythonProject3/main.py", line 125, in NS
    solve(A == b, u_s, bcu)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/la/__init__.py", line 76, in __eq__
    return self.get_local() == value.get_local()
AttributeError: 'dolfin.cpp.la.Matrix' object has no attribute 'get_local'

Looking forward and thanks for your reply

[dokken]

@Kan2046

solve(A == b, u_s, bcu)

I am not sure what kind of syntax you are trying to use here. As A is a matrix, and b a vector, you should use the solve command that you have commented:

@Kan2046

     #解法二
        #solve(A, u_.vector(), b)
        #u_n.assign(u_)

[Kan2046]

@dokken

am not sure what kind of syntax you are trying to use here. As A is a matrix, and b a vector, you should use the solve command that you have commented

非常非常感谢您的指导！希望我的谢意可以跨过欧亚大陆。
Thank you very much for your guidance! I hope my gratitude goes across Eurasia.

[Kan2046]

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

[dokken]

There are many things to consider.
1: you are looking at very coarse meshes. Consider increasing the number of elements further to make sure the code converges to your analytical solution
2. is your time step sufficiently small?
3. have you considered any of the benchmarks covered in: The Navier-Stokes equations — FEniCSx tutorial

As your code has quite complicated expressions that can be prone to typos.

[Kan2046]

Dear Mr.dokken
Thank you again for your reply.
I tried dense grid again [8, 16, 32] and the results is bad.
But the result is still disappointing.
Time = 1, which should be small enough.And h=1/N^2, this should enough?
I’ll study the problem carefully again.

What I want to ask is, will there be a problem of program structure that will lead to this result? Is it also related to the iterative form?

Thank you for your reply.

[dokken]

I do not understand why you have:

@Kan2046

10e-8*inner(p,q))*dx

In your variational form. You should pin the pressure at a single degree of freedom, see: Pointsource on specific points python - #2 by kamensky

There are various sources that have already considered solving the Navier stokes equation with dolfin, including:

978-3-319-52462-7.pdf

Chapter 3.4

You can also find a similar problem in:
https://www.sciencedirect.com/science/article/abs/pii/S0045782520303145
Solving the Taylor-Green vortex problem.
The source code is available at:

Source code for: A multimesh finite element method for the Navier-Stokes...

[Kan2046]

Dear Mr.Dokken
When correcting my code, I met a question, why some NS programs contain the exact solution of u and p, and can calculate the exact value of f, but f is still written in the program as
F = constant ((0.0))
This is the code format for community questions, right? Or some solution to the equation?

Thank you for your reply!

[dokken]

You would need to reference where this is the case. If you are using a manufactured solution, f should be a function of u_ex and p_ex.

[sun1130]

I think you should check whether your exact solution u is divergence-free, i.e, \mathrm{div}\,u=0.

[Kan2046]

Thanks Mr.Sun1130 reply. Let me learn new knowledge.

[Kan2046]

This is my new code, which is still used to solve the timeNS problem, but this time I changed the previous mixed space into an independent space to calculate the code. I carefully checked the previous mistakes made by Mr.Dokken, but the result is still unsatisfactory, and the L2 error order is still not 2. Here is my new code, I hope someone can help me! Thank you all very much!

ut - delta(u) + (u.nable)*u + nable(p) = f
div(u) = 0
time-dependent Navier–Stokes equations

from __future__ import print_function
from fenics import *
from mshr import *
import numpy as np
from dolfin import *
import math

def NS(N):
    # 定义网格
    T = 1.0
    num_steps = pow(N, 2)
    dt = T / num_steps

    mesh = UnitSquareMesh(N, N)

    #定义函数空间 p2-p1
    #V = VectorFunctionSpace(mesh, 'P', 2)
    #Q = FunctionSpace(mesh, 'P', 1)

    #p1b-p1
    P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
    B = FiniteElement("Bubble", mesh.ufl_cell(), 3)
    Vh1 = VectorElement(P1 + B)
    Vh2 = P1
    V = FunctionSpace(mesh, Vh1)
    Q = FunctionSpace(mesh, Vh2)

    t=0
    #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)*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((u1_ex,u2_ex), degree=3, t=t)

    # P条件
    p1_ex = """10*(2*x[0]-1)*(2*x[1]-1)*cos(t)"""
    p_D = Expression(p1_ex, degree=3, t=t)

    # 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=t)
    #f=Constant((0,0))

    # Define trial and test functions
    u = TrialFunction(V)
    p = TrialFunction(Q)
    v = TestFunction(V)
    q = TestFunction(Q)

    #定义边界条件
    noslip = DirichletBC(V, u_D,
                         "on_boundary && \
                           (x[0] < DOLFIN_EPS | x[1] < DOLFIN_EPS | \
                           (x[0] > 1.0- DOLFIN_EPS && x[1] > 1.0 - DOLFIN_EPS))")              #黏附条件
    inflow = DirichletBC(Q, p_D, "x[1] > 1.0 - DOLFIN_EPS")
    outflow = DirichletBC(Q, 0, "x[0] > 1.0 - DOLFIN_EPS")
    bcu = [noslip]
    bcp = [inflow, outflow]

    # Create functions
    u0 = Function(V)
    u1 = Function(V)
    p1 = Function(Q)

    # Define coefficients
    k = Constant(dt)

    # Tentative velocity step
    F1 = (1 / k) * inner(u - u0, v) * dx + inner(grad(u0) * u0, v) * dx + \
         inner(grad(u), grad(v)) * dx - inner(f, v) * dx
    a1 = lhs(F1)
    L1 = rhs(F1)

    # Pressure update
    a2 = inner(grad(p), grad(q)) * dx
    L2 = -(1 / k) * div(u1) * q * dx

    # Velocity update
    a3 = inner(u, v) * dx
    L3 = inner(u1, v) * dx - k * inner(grad(p1), v) * dx

    # Assemble matrices
    A1 = assemble(a1)
    A2 = assemble(a2)
    A3 = assemble(a3)

    # Use amg preconditioner if available
    prec = "amg" if has_krylov_solver_preconditioner("amg") else "default"

    # Use nonzero guesses - essential for CG with non-symmetric BC
    parameters['krylov_solver']['nonzero_initial_guess'] = True

    # Time-stepping
    t = dt
    for N in range(num_steps):
        # Update pressure boundary condition
        p_D.t = t

        # Compute tentative velocity step
        b1 = assemble(L1)
        [bc.apply(A1, b1) for bc in bcu]
        solve(A1, u1.vector(), b1, "bicgstab", "default")

        # Pressure correction
        b2 = assemble(L2)
        [bc.apply(A2, b2) for bc in bcp]
        [bc.apply(p1.vector()) for bc in bcp]
        solve(A2, p1.vector(), b2, "bicgstab", prec)

        # Velocity correction
        b3 = assemble(L3)
        [bc.apply(A3, b3) for bc in bcu]
        solve(A3, u1.vector(), b3, "bicgstab", "default")

        # Move to next time step
        u0.assign(u1)

        u_error_L2 = errornorm(u_D, u1, 'L2')
        u_error_H1 = errornorm(u_D, u1, 'H1')

        p_error_L2 = errornorm(p_D, p1, 'L2')
        p_error_H1 = errornorm(p_D, p1, 'H1')

        t += dt

    return  [u_error_L2 , u_error_H1, p_error_L2, p_error_H1,N]

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]=NS(n1)
[a2,b2,c2,d2,e2]=NS(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))

Thank you again.