(Using legacy FEniCS): I am trying to numerically solve a two-species advection-diffusion system involving four nonlinear PDEs (two for each species), on a rectangular domain with periodic conditions on the left+right wall, and no flux top and bottom. I setup the mesh, elements and functions spaces as:
Nx, Ny = 40, 20
mesh = RectangleMesh(Point(0.0,0.0), Point(1.0, 0.5), Nx, Ny)
class PeriodicLeftRight(SubDomain):
def inside(self, x, on_bdry):
return on_bdry and near(x[0], 0.0)
def map(self, x, y):
y[0] = x[0] - 1.0
y[1] = x[1]
per_map = PeriodicLeftRight()
P1 = FiniteElement("CG", mesh.ufl_cell(), 1)
P1_vec = VectorElement("CG", mesh.ufl_cell(), 1)
Z = FunctionSpace(mesh, MixedElement([P1, P1, P1_vec, P1_vec]),
constrained_domain=per_map)
z = Function(Z)
z_prev = Function(Z)
rho1, rho2, q1, q2 = split(z)
rho1_n, rho2_n, q1_n, q2_n = split(z_prev)
Here rho1, rho2 are the densities of the two species, and q1, q2 can be thought of as “momentum”.
I set-up the following variational form (with previous simple regularisations due to a singular term in the PDEs that I don’t include here) :
F = 0
F += ((rho1 - rho1_n) / dt * v_rho1
- dot(q1, grad(v_rho1))
+ epsilon * dot(rho1 * grad_phi, grad(v_rho1))) * dx
F += ((rho2 - rho2_n) / dt * v_rho2
- dot(q2, grad(v_rho2))
+ epsilon * dot(rho2 * grad_phi, grad(v_rho2))) * dx
F += (dot((q1 - q1_n) / dt, v_q1)
- inner(outer(q1, q1) / (rho1_pos + DEN), grad(v_q1))
+ epsilon * inner(outer(grad_phi, q1), grad(v_q1)) )*dx
F += dot(q1, n) * dot((q1 - epsilon * rho1_pos * grad_phi) / (rho1_pos + DEN), v_q1) * ds_wall(1)
F += (dot((q2 - q2_n) / dt, v_q2)
- inner(outer(q2, q2) / (rho2_pos + DEN), grad(v_q2))
+ epsilon * inner(outer(grad_phi, q2), grad(v_q2)))*dx
F += dot(q2, n) * dot((q2 - epsilon * rho2_pos * grad_phi) / (rho2_pos + DEN), v_q2) * ds_wall(1)
J = derivative(F, z)
I solve this with snes solver, with linesearch “cp”, linear solver as “gmres” and preconditioner “hypre_amg”.
I then setup the initial configuration with:
A_block = (0.0, 0.45, 0.1, 0.5)
B_block = (0.55, 1.0, 0.0, 0.4)
class BlockScalarA(UserExpression):
def __init__(self, blocks, value_inside=0.5, **kw):
super().__init__(**kw); self.blocks = blocks; self.v = value_inside
def eval(self, value, x):
value[0] = self.v if any(b[0] < x[0] < b[1] and b[2] < x[1] < b[3]
for b in self.blocks) else 0.0
def value_shape(self): return ()
class BlockScalarB(UserExpression):
def __init__(self, blocks, value_inside=0.5, **kw):
super().__init__(**kw); self.blocks = blocks; self.v = value_inside
def eval(self, value, x):
value[0] = self.v if any(b[0] < x[0] < b[1] and b[2] < x[1] < b[3]
for b in self.blocks) else 0.0
def value_shape(self): return ()
class BlockVector(UserExpression):
def __init__(self, blocks, vx, vy=0.0, **kw):
super().__init__(**kw); self.blocks = blocks; self.vx = vx; self.vy = vy
def eval(self, value, x):
if any(b[0] < x[0] < b[1] and b[2] < x[1] < b[3] for b in self.blocks):
value[0], value[1] = self.vx, self.vy
else:
value[0] = value[1] = 0.0
def value_shape(self): return (2,)
Then do the time-stepping and animate.
Block A corresponds to species 1, Block B corresponds to species 2.
Now to the problem I am facing: I want Block A to move to the right and Block B to move to the left. Now Block A moves as desired, but Block B just remains stationary and diffuses. When A collides with Block B though, the collision happens as expected and nothing weird seems to happen.
To try and see what was going wrong, I then set both Block A and B to move to the right (i.e both positive vx) - this worked as desired, and the periodic walls also functioned properly.
But when I set both to move to the left (negative vx), both remained stationary and only diffused. It seems the problem is that a negative vx leads to no motion, and I would appreciate any help or directions in how to fix this.