Vectorial Point Source Term (Stokes equations)

General
, ,
1

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.

2

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.

3

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()

dipole
dipole1920×1449 204 KB

4

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.

5

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.

6

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.

7

Hi, I’m still working on the same problem and I think the example code and plot I posted before still qualifies as a minimal example. My problem now is that when I place my dipole, that right now it’s surrounded by 8 additional point forces, close to the edges of the domain the solution does not converge. The program times out after 10000 iterations.

One difference of my current code and the one posted is that now I’m defining the mesh like shown below and I’m refining my mesh in the center of the dipole since I need more resolution close to all my point sources:

p0 = Point([-0.5, -2.5, -0.5])  
p1 = Point([0.5, 2.5, 0.5])  
nx = 10
ny = 50
nz = 10
mesh = BoxMesh(p0, p1, nx, ny, nz)
for _ in range(2):
    self.mesh = self.refine_mesh(Point( (0.45+0.55)/2, 0.5, 0.5), radius = 0.2)

with:

def refine_mesh(self, point_force, radius):
        cell_markers = MeshFunction("bool", self.mesh, self.mesh.topology().dim())
        cell_markers.set_all(False)
        for cell in cells(self.mesh):
            if cell.midpoint().distance(point_force) < radius:
                cell_markers[cell] = True
        return refine(self.mesh, cell_markers)

I guess the program is crashing when close to the boundary because non of my point sources are nodes of the mesh and maybe because the refining is not enough. If I don’t apply the refinement it converges but my solution is garbage.
I’ve been trying for weeks to either switch to DOLFINx, where I’ve heard unstructured meshes are better implemented, or to create a mesh using gmsh where my point sources form part of the mesh, but I’m having trouble finding how to modify my code to run it in DOLFINx or importing the mesh I created with gmsh to legacy-dolfin. I’ve read a lot of threads here and I just replicate other people’s errors again and again.
Can someone help me with this?
Thanks.

8

Have you considered
https://scientificcomputing.github.io/scifem/examples/point_source.html

9

I recently revisited the tutorial you mentioned, either in a post or while answering someone else’s question. I liked it when I first saw it, but to be honest, it left me with more questions than answers.

I only discovered the FEniCS project last year, and even after going through all the tutorials (legacy-version), I’m still unclear about how everything works—especially when it comes to solving coupled systems like the Stokes equations. It seems to involve more than just constructing a and L, and I haven’t been able to find clear explanations of the process.

I’ve also struggled with the transition to DOLFINx. While the DOLFINx tutorial is helpful, it doesn’t seem to cover all the documentation, and beyond the tutorial, I haven’t found many additional resources. I’ve thought about diving directly into the implementation .py files, but I feel like there must be a better way to learn.

Could you point me to resources or documentation that cover transitioning to DOLFINx, resolving the Stokes equations, and finding all available functions for mesh generation?

I’ve read so many posts that you’ve answered by now, that I truly feel you’ve already helped me a lot. Thanks again for your help and for your work in the FEniCS project.

10

We have a specific demo for the Stokes equation at:
https://docs.fenicsproject.org/dolfinx/v0.9.0/python/demos/demo_stokes.html
and recently I’ve used a similar technique at:

There is even more «basic» documentation at
http://jsdokken.com/FEniCS-workshop/README.html
and in condensed version at: FEniCS tutorial for Sorbonne University — FEniCS Tutorial @ Sorbonne