New question #155893 on DOLFIN:
https://answers.launchpad.net/dolfin/+question/155893
I have created code for the Steady-state Lid Driven (enclosed flow) for Navier
Stokes. I have a convergence error in the Newton Solver. Is this because I have
no pressure BC? Should I set a pressure to the entire domain or something like
this?
The code is below:
#!/usr/bin/python
from dolfin import *
import numpy, sys
# MESHING
mesh = Rectangle(0.0, -0.5, 1.0, 0.5, 10, 10)
# FUNCTION SPACE
scalar = FunctionSpace(mesh, "CG", 1) #Pressure
vector = VectorFunctionSpace(mesh, "CG", 2) #Velocity
system = vector * scalar #Mixed Function Space
D = mesh.topology().dim()
# BOUNDARIES
right = compile_subdomains('x[0] == 1.0')
left = compile_subdomains('x[0] == 0.0')
top = compile_subdomains('x[1] == 0.5')
bottom = compile_subdomains('x[1] == -0.5')
subdomains = MeshFunction("uint", mesh, D-1)
subdomains.set_all(0)
# INFLOW VELOCITY BC FOR TOP
bc1 = DirichletBC(system.sub(0), Constant((1.0, 0.0)), top)
# NO-SLIP BC FOR TOP-BOTTOM
noslip = Constant((0.0, 0.0))
bc0 = DirichletBC(system.sub(0), noslip, left)
bc3 = DirichletBC(system.sub(0), noslip, bottom)
bc2 = DirichletBC(system.sub(0), noslip, right)
# STORAGE OF BCs
bcs = [bc1, bc2, bc0, bc3]
# VARIATIONAL PROBLEM
class NavierStokes(NonlinearProblem):
def __init__(self, a, L, bc, exter_B):
NonlinearProblem.__init__(self)
self.L = L
self.a = a
self.bc = bc
self.exter = exter_B
self.reset_sparsity = True
def F(self, b, x):
assemble(self.L, tensor=b, exterior_facet_domains=self.exter)
for condition in self.bc: condition.apply(b, x)
def J(self, A, x):
assemble(self.a, tensor=A, exterior_facet_domains=self.exter)
for condition in self.bc: condition.apply(A)
self.reset_sparsity = False
parameters["form_compiler"]["cpp_optimize"] = True
ts = TestFunction(system)
ds = TrialFunction(system)
s = Function(system)
v, q = split(ts)
dv, dq = split(ds)
u , p = split(s)
f = Constant((0, 0))
nu = 0.1
F = (nu*inner(grad(v), grad(u)) - div(v)*p + q*div(u) + v[i]*u[j]*Dx(u[i], j) -
inner(v,f))*dx # Navier-Stokes
J = derivative(F,s,ds)
# SET UP PDE
problem = NavierStokes(J, F, bcs, subdomains)
solver = NewtonSolver("lu")
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6
# SOLVE PDE
solver.solve(problem,s.vector())
(U, P) = s.split()
# FILES FOR PARAVIEW
ufile_pvd = File("velocity.pvd")
ufile_pvd << U
pfile_pvd = File("pressure.pvd")
pfile_pvd << P
# PLOTS
plot(U)
plot(P)
interactive()
--
You received this question notification because you are a member of
DOLFIN Team, which is an answer contact for DOLFIN.
_______________________________________________
Mailing list: https://launchpad.net/~dolfin
Post to : dolfin@lists.launchpad.net
Unsubscribe : https://launchpad.net/~dolfin
More help : https://help.launchpad.net/ListHelp
