Dear all,
i tried to implement Chorin's method to solve NS problem for square driven cavity, but i'm not sure is correct.
from dolfin import *
mesh = UnitSquareMesh(129,129)
V = VectorFunctionSpace(mesh, "Lagrange", 2)
Q = FunctionSpace(mesh, "Lagrange", 1)
Define trial and test functions
u = TrialFunction(V)
p = TrialFunction(Q)
v = TestFunction(V)
q = TestFunction(Q)
Functions
u0 = Function(V)
u1 = Function(V)
p1 = Function(Q)
boundary conditions
def cavityns(x, on_boundary): return x[0] < DOLFIN_EPS or x[0] > 1- DOLFIN_EPS or x[1] < DOLFIN_EPS
def cavitytop(x, on_boundary): return x[1] > 1.0 - DOLFIN_EPS
noslip
noslip = Constant((0.0, 0.0))
bc0 = DirichletBC(V, noslip, cavityns)
lid
lid = Constant((1.0, 0.0))
bcl = DirichletBC(V, lid, cavitytop)
pressure
pref = DirichletBC(Q, 0, "x[0] < DOLFIN_EPS && x[1] < DOLFIN_EPS", "pointwise")
collect BC
bcu = [bc0, bcl]
bcp = [pref]
bc = [bcu, bcp]
parameters
n = FacetNormal(mesh)
nu = Constant(0.01)
f = Constant((0.0, 0.0))
Set parameter values
dt = 0.1
T = 1
Define coefficients
k = Constant(dt)
Tentative velocity step
a1 = (1/k)inner(u, v)dx + nuinner(grad(u), grad(v))dx
a1 = a1 + inner(grad(u)u0, v)dx
L1 = -inner(grad(u0)u0, v)dx + (1/k)inner(u0, v)dx
L1 = L1 + inner(f, v)*dx
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 - kinner(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
Create files for storing solution
ufile = File("results/velocity.pvd")
pfile = File("results/pressure.pvd")
#
Iterazioni
#
eps = 1.0
tol = 1.0e-5
niter = 0
maxiter = 100
#
Iterations
#
u_old = Function(V)
u_new = Function(V)
Time-stepping
t = dt
while t <= T:
begin("Computing tentative velocity")
while eps > tol and niter < maxiter:
a1 = (1/k)inner(u, v)dx + nuinner(grad(u), grad(v))dx
tri = (u_old[j]Dx(u[i], j)v[i])dx #trilinear form
L1 = inner(f, v)dx + (1/k)inner(u0, v)dx
A = assemble(a1)
b1 = assemble(L1)
C = assemble(tri)
A1 = A + C
[bc.apply(A1, b1) for bc in bcu]
solve(A1, u_new.vector(), b1, "bicgstab", "default")
u_ex= interpolate(u_ex, V)
eps = errornorm(u_old,u_new)
print "iter=%d: norm=%g" % (niter, eps)
u_old.assign(u_new)
niter += 1
end()
u1.assign(u_new)
niter = 0
eps = 1
# Pressure correction
begin("Computing pressure correction")
a2 = inner(grad(p), grad(q))dx
L2 = -(1/k)div(u1)qdx
A2 = assemble(a2)
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)
end()
# Velocity correction
begin("Computing velocity correction")
a3 = inner(u, v)*dx
L3 = inner(u1, v)*dx - k*inner(grad(p1), v)*dx
A3 = assemble(a3)
b3 = assemble(L3)
[bc.apply(A3, b3) for bc in bcu]
solve(A3, u1.vector(), b3, "bicgstab", "default")
end()
# Plot solution
plot(p1, title="Pressure", rescale=True)
plot(u1, title="Velocity", rescale=True)
# Save to file
ufile << u1
pfile << p1
# Move to next time step
u0.assign(u1)
t += dt
print("t =", t)
Hold plot
interactive()
Best reguards,
Francesco
