Setting Dirichlet Boundary Conditions from Previous (Adjacent) Simulations?

Hi all,

TL;DR:

  1. Why isn’t the pressure boundary condition, with values from a previous solution, being applied on the border?
  2. For near(x[0], 2.0), why is x[0] = 1.95 being called?
  3. How can I fix my code, i.e. apply the pressure boundary condition on the border?

I am trying to simulate fluid flow (Navier-Stokes without time dependency) in a long pipe by splitting it into several smaller simulations, as seen in the image below. I am using Dirichlet Boundary Conditions for velocity on the left and pressure on the right of each pipe.

In the first simulation, I assume a velocity and pressure boundary. The second simulation is connected to the first spatially, so I want to use the velocity profile on the right side of the first simulation as the profile on the left side of the second simulation, again setting pressure to a constant. This is working more-or-less as intended.

Figure 1: Splitting a long pipe into two sections, using adjacent edges from previous simulations as Dirichlet boundary conditions on future ones

In the third simulation, I want to use the pressure from the left side of the second simulation as the pressure profile on the right side of the third simulation. However, this part is not running correctly. It seems the pressure is still being set to 0, and the values being called to apply the boundary condition are not on the boundary.

My question has 3 parts.

  1. Why isn’t the pressure boundary condition being applied on the border?
  2. In simulation 3, why is x[0] = 1.95 being called?
  3. How can I fix this to run as intended, i.e. apply the pressure boundary condition on the border?

Below is my best attempt at a MWE. I am running FEniCS 2018.2.0.dev0 on Ubuntu 18.04.1.

from fenics import *

class borderVelBC(UserExpression):
    def getPrevSoln(self, prevSol):
        self.prevSol = prevSol
        
    def eval(self, values, x):
        try:
            val = self.prevSol(x[0], x[1])
        except:
            val = [0, 0]
        values[0] = val[0]
        values[1] = val[1]
        
    def value_shape(self):
        return (2,)
    
class borderPreBC(UserExpression):
    def getPrevSoln(self, prevSol):
        self.prevSol = prevSol
        
    def eval(self, values, x):
        try:
            val = self.prevSol(x[0], x[1])
            print(str(val) + " at " + str(x))
        except:
            print("attempt invalid at " + str(x))
            val = -99
        values = val
        
########## Simulation 1 ##########
lefBound = 0.0
rigBound = 2.0
botBound = 0.0
topBound = 1.0

mesh = RectangleMesh(Point(lefBound, botBound), Point(rigBound, topBound), 40, 20)

V = VectorElement("CG", mesh.ufl_cell(), 2)
Q = FiniteElement("CG", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, MixedElement((V, Q)))

noSlip = DirichletBC(W.sub(0), Constant((0.0, 0.0)), "near(x[1], 0.0) || near(x[1], 1.0)")
leftBC = DirichletBC(W.sub(0), Constant((0.01, 0.0)), "near(x[0], 0.0) && x[1] < 0.9 && x[1] > 0.1")
rightBC = DirichletBC(W.sub(1), Constant(0.0), "near(x[0], 2.0)")

bcs = [noSlip, leftBC, rightBC]

w = Function(W)

(u, p) = split(w)
(v, q) = TestFunctions(W)

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

solve(F == 0, w, bcs = bcs)

(u_soln, p_soln) = w.split(deepcopy=True)

File("u1.pvd") << u_soln
File("p1.pvd") << p_soln

prevUSoln = u_soln
prevPSoln = p_soln

########## Simulation 2 ##########
lefBound = 2.0
rigBound = 4.0
botBound = 0.0
topBound = 1.0

mesh = RectangleMesh(Point(lefBound, botBound), Point(rigBound, topBound), 40, 20)

V = VectorElement("CG", mesh.ufl_cell(), 2)
Q = FiniteElement("CG", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, MixedElement((V, Q)))

noSlip = DirichletBC(W.sub(0), Constant((0.0, 0.0)), "near(x[1], 0.0) || near(x[1], 1.0)")
rightBC = DirichletBC(W.sub(1), Constant(0.0), "near(x[0], 4.0)")

bcFromSim = borderVelBC()
bcFromSim.getPrevSoln(prevUSoln)
leftBC = DirichletBC(W.sub(0), bcFromSim, "near(x[0], 2.0)")

bcs = [noSlip, leftBC, rightBC]

w = Function(W)

(u, p) = split(w)
(v, q) = TestFunctions(W)

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

solve(F == 0, w, bcs = bcs)

(u_soln, p_soln) = w.split(deepcopy=True)

File("u2.pvd") << u_soln
File("p2.pvd") << p_soln

prevUSoln = u_soln
prevPSoln = p_soln

########## Sim. 1 Correction ##########
lefBound = 0.0
rigBound = 2.0
botBound = 0.0
topBound = 1.0

mesh = RectangleMesh(Point(lefBound, botBound), Point(rigBound, topBound), 40, 20)

V = VectorElement("CG", mesh.ufl_cell(), 2)
Q = FiniteElement("CG", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, MixedElement((V, Q)))

noSlip = DirichletBC(W.sub(0), Constant((0.0, 0.0)), "near(x[1], 0.0) || near(x[1], 1.0)")
leftBC = DirichletBC(W.sub(0), Constant((0.01, 0.0)), "near(x[0], 0.0) && x[1] < 0.9 && x[1] > 0.1")

bcFromSim = borderPreBC()
bcFromSim.getPrevSoln(prevPSoln)
rightBC = DirichletBC(W.sub(1), bcFromSim, "near(x[0], 2.0)")

bcs = [noSlip, leftBC, rightBC]

w = Function(W)

(u, p) = split(w)
(v, q) = TestFunctions(W)

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

solve(F == 0, w, bcs = bcs)

(u_soln, p_soln) = w.split(deepcopy=True)

File("u1correction.pvd") << u_soln
File("p1correction.pvd") << p_soln

Additional images (new users can only post one image)

Figure 2: Pressure values of split simulations. All color mappings (red = high, blue = low) are in the same range. We expect the third simulation’s right edge and second simulation’s left edge to have the same values, which they do not.

Figure 3: Velocity magnitude values of split simulations. These results are expected.

Figure 4: Some output from running this code. Note that in Sim. 1 Correction, the BC is applied where near(x[0], 2.0) but is still being queried at x[0] = 1.95.
image

I found the issue. In the borderPreBC class, I was using values = val.
It should be values[0] = val.