Hello,
Let’s summarize my issue, I’ve been trying to solve the Toner-Tu equation, basically the Navier-Stokes equation with a nonlinear term, and it worked pretty well. Previously I was solving with Dirichlet boundary conditions.
Now I want to implement Periodic Boundary Conditions. I’ve faced two major issues. First, the solution returned ain’t periodic. Then, after multiple time steps, I a peak in density appears on the border and few steps later the simulations diverge. I join an Image of what’s happening then…
I work with mixed space and it wasn’t a problem before. I join some parts of my code in case you have any ideas of where it could come from…
The implementation of the periodic boundary conditions:
class PeriodicBoundary(SubDomain):
def inside(self, x, on_boundary):
return bool((near(x[0], SIZE/2) or near(x[1], SIZE/2))
and not(near(x[0], -SIZE/2) or near(x[1], SIZE/2)) and on_boundary)
def map(self, x, y):
if near(x[0], SIZE/2) and near(x[0], SIZE/2):
y[0] = x[0]-SIZE
y[1] = x[1]-SIZE
elif near(x[0], SIZE/2):
y[0] = x[0]-SIZE
y[1] = x[1]
elif near(x[1], SIZE/2) and not(near(x[0], -SIZE/2) or near(x[0], SIZE/2)):
y[0] = x[0]
y[1] = x[1]-SIZE
else:
y[0] = x[0]-100*SIZE
y[1] = x[1]-100*SIZE
Then I implement the domain:
V1 = VectorElement("Lagrange", mesh.ufl_cell(), DEG)
V2 = FiniteElement("Lagrange", mesh.ufl_cell(), DEG)
V = VectorFunctionSpace(mesh, 'P', DEG)
Vs = FunctionSpace(mesh, 'P', DEG)
# Make a mixed space
TH = V1*V2
if CL_P:
# Create periodic boundary condition
pbc = PeriodicBoundary()
W = FunctionSpace(mesh, TH, constrained_domain=pbc)
V = VectorFunctionSpace(mesh, 'P', DEG, constrained_domain=pbc)
Vs = FunctionSpace(mesh, 'P', DEG, constrained_domain=pbc)
And even if it’s ugly the form that worked pretty well before:
lamb = Constant(LAMBDA) # SD
alph = Constant(1e2) # s-1
pho_c = Constant(3e-3) # SD
pho = Constant(0.1) # SD
a2 = Constant(alph*(pho-pho_c)) # s-1
a4 = Constant(10) # mm-2 s
D = Constant(1e-2) # mm2 s-1
D2 = Constant(4e-6*FAC_D_RHO) # mm2 s-1
sigma = Constant(5) # mm2 s-2
Dt = Constant(DT)
# Declare the form
F = dot((u-u_1)/Dt, v)*dx \
+lamb*dot(dot(u, nabla_grad(u)), v)*dx \
+D*(dot(grad(u[0]), grad(v[0]))+dot(grad(u[1]), grad(v[1])))*dx \
-D*(dot(grad(u[0]),n)*v[0]+dot(grad(u[1]),n)*v[1])*ds\
-(alph*(p-pho_c)-a4*dot(u, u))*dot(u, v)*dx \
+sigma*dot(grad(p), v)*dx \
+(((p-p_1)/Dt+div(p*u))*q+D2*dot(grad(p), grad(q)))*dx \
-D2*dot(grad(p),n)*q*ds
np.savetxt("./Solutions/coordo.txt", np.reshape(np.array(x.vector()), (-1, 2)))
for i in range(0, NUM_STEP):
plotfigure = False
try:
if CL_P:
solve(F == 0, w)
else:
solve(F == 0, w, bcu)
except RuntimeError:
break
assign(w_1, w)
Any help would be appreciated.
Thanks!
