Hello,
I am currently trying to solve a steady Stokes problem. And I am testing different methods and preconditioners. I posted below the code I am using.
Solving this problem using the default solve(a == L, w, bcs) gives expected result but when I want to use CG method, the velocity field is incorrect. What is wrong here? Thank you!
from mshr import *
from dolfin import *
# Adjust log level
set_log_level(PROGRESS)
# Set linear algebra backend
parameters["linear_algebra_backend"] = "PETSc"
# Set global DOLFIN parameters
flags = ["-O3", "-ffast-math", "-march=native"]
parameters["form_compiler"]["quadrature_degree"] = 4
parameters["form_compiler"]["representation"] = "uflacs"
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["cpp_optimize_flags"] = " ".join(flags)
parameters["form_compiler"]["precision"] = 8
TEST=1
if TEST:
cylinder = Cylinder(Point(0, 0, 0), Point(0, 0, 3.0), 0.8, 0.7)
geometry = cylinder
mesh = generate_mesh (geometry , 10)
d=0.1
# coordinates
V = FunctionSpace(mesh, "CG", 2)
ndim = V.dim()
nd = mesh.geometry().dim()
dof_coordinates = V.dofmap().tabulate_all_coordinates(mesh)
dof_coordinates.resize((ndim, nd))
dof_z = dof_coordinates[:, 2]
facetZ_max=max(dof_z)
facetZ_min=min(dof_z)
# Construct facet markers
boundaries = FacetFunction("size_t", mesh)
class Inflow(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x[2] < facetZ_min+d
class Outflow(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x[2] > facetZ_max-d
class WallOut(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x[2] >= facetZ_min+0. and x[2] <= facetZ_max-0.
WallOut=WallOut()
Inflow=Inflow()
Outflow=Outflow()
boundaries.set_all(0)
WallOut.mark(boundaries, 3)
Inflow.mark(boundaries, 1)
Outflow.mark(boundaries, 2)
# Measurements
n = FacetNormal(mesh)
ds = Measure('ds')[boundaries]
# Define function spaces
P2 = VectorFunctionSpace(mesh, "CG", 2)
P1 = FunctionSpace(mesh, "CG", 1)
W = MixedFunctionSpace([P2,P1])
# No-slip boundary conditions
noslip = Constant(( 0.0 , 0.0 , 0.0 ))
bc0 = DirichletBC(W.sub(0), noslip, boundaries, 3)
bcs = [bc0]
# Define variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
p_out=0; p_in=Constant(10)
a = inner(grad(u), grad(v))*dx + div(v)*p*dx + q*div(u)*dx
L = - inner(p_out, dot(v, n))*ds(2, domain=mesh) - inner(p_in, dot(v, n))*ds(1, domain=mesh)
(A, b) = assemble_system(a, L, bcs)
krylov_method="cg"
precond="amg"
solver = KrylovSolver(krylov_method, precond)
solver.parameters["relative_tolerance"] = 1.0e-8
solver.parameters["absolute_tolerance"] = 1.0e-6
solver.parameters["monitor_convergence"] = True
solver.parameters["maximum_iterations"] = 1000
print("Solver: " + krylov_method + "... Preconditioner: " + precond)
# Compute solution
w = Function(W)
solver.solve(A, w.vector(), b)
# solve(a == L, w, bcs)
# Check mass conservation
(u, p) = w.split(True)
#inlet flow
form = dot(n,u)*ds(1)
ans1 = assemble(form)
print("\nInlet flow: %.10g" % ans1)
#outlet flow
form = dot(n,u)*ds(2)
ans2 = assemble(form)
print("\nOutlet flow: %.10g" % ans2)
plot(u)
interactive()
