Steady State Stokes Equations with Iterative Solver

mathematics
1

I’m following the demo for the Stokes equations with an iterative solver and I can get that to run. However, if I try to reformulate the problem as below, it appears that my velocity boundary condition gets applied correctly on the elements that share a face with the boundary, but the pressure and the remainder of the velocity solution remain close to 0. My question is, is there a way to modify the existing code below to solve the problem correctly?

from fenics import *
mesh = UnitCubeMesh(25, 25, 25)

CG2 = VectorElement("CG", mesh.ufl_cell(), 3)
CG1 = FiniteElement("CG", mesh.ufl_cell(), 2)
TH = MixedElement([CG2, CG1])
W = FunctionSpace(mesh, TH)
noslip = DirichletBC(W.sub(0), (0.0, 0.0, 0.0), "on_boundary && ( x[1] < DOLFIN_EPS | x[1] > 1.0 - DOLFIN_EPS | x[2] < DOLFIN_EPS | x[2] > 1 - DOLFIN_EPS)")
inflow = DirichletBC(W.sub(0), (1.0, 0.0, 0.0), "x[0] < DOLFIN_EPS")
ambient_pressure = DirichletBC(W.sub(1), 0.0, "x[0] > 1 - DOLFIN_EPS")
bcs = [noslip,
       inflow,
       ambient_pressure]

mu = Constant(0.89e-3)
rho = Constant(997.0)
f = Constant((0.0, 0.0, 0.0))
nu = Constant(mu.values()[0]/rho.values()[0])

w = Function(W)
(u_vel, u_pre) = split(w)
(v_vel, v_pre) = TestFunctions(W)
F = nu*inner(grad(u_vel), grad(v_vel))*dx + div(v_vel)*u_pre*dx + v_pre*div(u_vel)*dx - inner(f, v_vel)*dx
solve(F == 0, w, bcs, solver_parameters={"newton_solver":
                                        {"relative_tolerance": 1e-4,
                                        'linear_solver': 'gmres',
                                        'preconditioner': 'amg'}})

u_vel, u_pre = w.split(deepcopy=True)

# Store velocity field data
ufile = File(f"test/velocity_stokes.pvd")
ufile << u_vel

# Store pressure field data
pfile = File(f"test/pressure_stokes.pvd")
pfile << u_pre
2

Your problem has a saddle point. Where’s the preconditioning matrix (the matrix P in your linked example code)?

1 Like
3

Hi Nate–thanks for your reply. Your question is kinda the information I’m looking for. The code above is in a format that seems pretty consistent between a lot of the FEniCS examples (i.e. solve(F == 0, w, bcs, ...) ), whereas they KrylovSolver creation in the example is called in a different format. I don’t mind using the format used in the linked example, but I am interested if I can achieve the same result while retaining the format of the code listed above. I may be way off base with trying to do things this way, and if I am, let me know!

I should also mention that eventually I want to solve the stationary Navier-Stokes equations with nonlinear term/s and I’m not sure what the best way to do this is, starting from the linked Stokes example.