Vectorial Point Source Term (Stokes equations)

Hi. I want to solve the Stokes equations in a unit cube with no slip boundary conditions in all faces of the cube and one point force in the center of the cube (pointing, for example, along the z axis). This is to solve the Green’s function of the Stokes flow. Bellow I’m sharing the code (using legacy dolfin) that I wrote for that purpose but I’m missing the part were I assign a vectorial point force to the center of the cube.

I’ve seen how you can create a PointSource like: delta = PointSource(V, point, 1.0) and then apply it to the rhs like: delta.apply(bb) but that’s not quite what I want to do. I don’t even know what is happening when I apply that point source that is a single scalar. What I want to have is f = (1.0, 0.0, 0.0) but applied to a single point in my domain. A rhs equal to:

f * n_vector * 3D_DiracDelta(x_vector -x0_vector)

where f is a scalar and n_vector is a unit vector representing the orientation of the point force in space.

Here is my code missing the point source term:

from dolfin import *

mesh = UnitCubeMesh(10, 10, 10)
# Define function spaces
V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)
W = FunctionSpace(mesh, MixedElement([V.ufl_element(), Q.ufl_element()]))

# Boundaries
def all_boundaries(x, on_boundary):
return x[0] > (1.0 - DOLFIN_EPS) or x[0] < (DOLFIN_EPS) or x[1] > (1.0 - DOLFIN_EPS) or x[1] < (DOLFIN_EPS) or x[2] > (1.0 - DOLFIN_EPS) or x[2] < (DOLFIN_EPS)
def top_y(x, on_boundary): return x[1] > (1.0 - DOLFIN_EPS)

noslip = Constant((0.0, 0.0, 0.0))
bc0 = DirichletBC(W.sub(0), noslip, all_boundaries)
bcs = [bc0]

# lid = Constant((1.0, 0.0, 0.0))
# bc1 = DirichletBC(W.sub(0), lid, top_y)
# bcs = [bc0, bc1]

# Define variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
f = Constant((0.0, 0.0, 0.0))

a = inner(grad(u), grad(v))*dx + div(v)*p*dx + q*div(u)*dx
L = inner(f, v)*dx

b = inner(grad(u), grad(v))*dx + p*q*dx

A, bb = assemble_system(a, L, bcs)
P, btmp = assemble_system(b, L, bcs)

solver = KrylovSolver("minres", "amg")
solver.set_operators(A, P)

U = Function(W)
solver.solve(U.vector(), bb)

u, p = U.split()

I’m new to Fenics. I started trying to use Fenicsx for this but as a newby everything looked more complicated there and even experienced users were having trouble with “scalar” point sources. Thanks in advance.

I haven’t tested this, but you may want to have a try with delta = PointSource(W.sub(0).sub(0), point, 1.0), where the first sub(0) is to get the velocity space out of the mixed space, and the second sub(0) is to get the x-component.

Thank you @francesco-ballarin. It worked exactly like you described. Here is the updated code and a plot of a cut-plane at z=0.5 of two of these point forces forming a dipole.
Now I have the problem of using real units. The intensity of the point source should be ~ 1e-12 and the dimensions or the BoxMesh different size lengths should be of the order of 1e-6: force in pN and spacial dimensions in um. If I just change this dimensions in the code I get an error saying:

line 121, in solve_stokes_equations_compact solver.solve(U.vector(), bb)
Error: Unable to solve linear system using PETSc Krylov solver.
*** Reason: Solution failed to converge in 0 iterations (PETSc reason DIVERGED_PC_FAILED, residual norm ||r|| = 0.000000e+00).
*** Where: This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0


*** DOLFIN version: 2019.2.0.13.dev0
*** Git changeset: ubuntu

Any advice here?
Here is the corrected code including @francesco-ballarin answer for the point source and a plot (described above):

from dolfin import *

mesh = UnitCubeMesh(10, 10, 10)
# Define function spaces
V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)
W = FunctionSpace(mesh, MixedElement([V.ufl_element(), Q.ufl_element()]))

# Boundaries
def all_boundaries(x, on_boundary):
return x[0] > (1.0 - DOLFIN_EPS) or x[0] < (DOLFIN_EPS) or x[1] > (1.0 - DOLFIN_EPS) or x[1] < (DOLFIN_EPS) or x[2] > (1.0 - DOLFIN_EPS) or x[2] < (DOLFIN_EPS)
def top_y(x, on_boundary): return x[1] > (1.0 - DOLFIN_EPS)

noslip = Constant((0.0, 0.0, 0.0))
bc0 = DirichletBC(W.sub(0), noslip, all_boundaries)
bcs = [bc0]

# lid = Constant((1.0, 0.0, 0.0))
# bc1 = DirichletBC(W.sub(0), lid, top_y)
# bcs = [bc0, bc1]

# Define the directional Point Sources
point0 = Point(0.45, 0.5, 0.5)
delta0 = PointSource(W.sub(0).sub(0), point0, -1.0)

point1 = Point(0.55, 0.5, 0.5)
delta1 = PointSource(W.sub(0).sub(0), point1, 1.0)

# Define variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
f = Constant((0.0, 0.0, 0.0))

a = inner(grad(u), grad(v))*dx + div(v)*p*dx + q*div(u)*dx
L = inner(f, v)*dx

b = inner(grad(u), grad(v))*dx + p*q*dx

A, bb = assemble_system(a, L, bcs)

# Apply point sources
delta0.apply(bb)
delta1.apply(bb)

P, btmp = assemble_system(b, L, bcs)

solver = KrylovSolver("minres", "amg")
solver.set_operators(A, P)

U = Function(W)
solver.solve(U.vector(), bb)

u, p = U.split()

Just make your problem dimensionless, that’s how the Reynolds number is introduced in the Navier-Stokes problem (i.e., the nonlinear version of the problem you are solving).

In principle you could define your mesh so that the box length is 1e-6 instead of 1 (as in UnitCubeMesh), but you would need to be careful (for instance: double check the value of DOLFIN_EPS, since it may very well be comparable to 1e-6). I wouldn’t bother, since it’s so simple to make the problem dimensionless and not having to fight with numbers like 1e-12, which are very close to the machine epsilon.

Thank you @francesco-ballarin. I made my problem dimensionless as suggested. Working with real units was making the code crash. Now everything is good.
Sorry for the new questions in the same post but, do you know how can I use the solution for “u” and “p” from my code and incorporate it as initial guess for my next iteration of the same problem so the solution converges quicker.
Let me explain myself better. In my simulations the center of mass of the dipole shown in the uploaded image moves tiny increments following certain equation of motion. What I want to do is to compute the velocity field generated by the dipole using the code uploaded before, and use this solution as initial guess at the next position calculation. In theory that should speed up things.

Any ideas? Thanks.

The nonlinear solver uses as initial guess whatever is stored in the solution dolfinx.fem.Functionobject . If you use the same object for the first solver and the second solver, the second one will automatically use the converged solution of the first one as initial guess. If you have further doubts about this please make a new topic.