For a Non-linear Schrodinger flow in a rectangular channel I want to implement periodic boundaries at the left and right side so that what flows out at the left side is the same as what flows in at the right side. At the top and bottom sides a homogenous Neumann boundary should be imposed. This should be a zero-flux BC. It is my understanding that if we have a homogenous neumann BC, or in other words a “do-nothing” boundary condition, nothing should be implemented in the FEniCS code. With this given information, it is correct to assume that only the periodic BC should be implemented in the FEniCS code? Will the Neumann BC be imposed correctly simply by doing nothing? The zero flux should only be imposed at the top and bottom sides, not left and right. The results we get looks very strange in the GIF simulation shown below.
My question is, is the boundary conditions correctly implemented in the code below and if so, what could be the problem causing the strange behaviour of the flow?
It can be noted that convergence tests has been done and the mass and energy is conserved with the Crank-Nicolsson method.
class PeriodicBoundary(SubDomain):
# Left boundary is "target domain" G
def inside(self, x, on_boundary):
return bool(x[0] < DOLFIN_EPS and x[0] > -DOLFIN_EPS and on_boundary)
# Map right boundary (H) to left boundary (G)
def map(self, x, y):
y[0] = x[0] - 1.0
y[1] = x[1]
# Create periodic boundary condition
pbc = PeriodicBoundary()
T = 1.0 # final time
num_steps = 100 # number of time steps
dt = T / num_steps # time step size
k=20
channel = Rectangle(Point(-2, 0), Point(2.5, 2))
cylinder = Ellipse(Point(0.8, 0.75), 0.5, 0.05)
domain = channel - cylinder
mesh = generate_mesh(domain, 64)
V = VectorFunctionSpace(mesh, 'P', 2, dim=2, constrained_domain=pbc)
# Define initial value
u_0 = Expression( ( 'cos(x[0])', 'sin(x[0])'), degree=1)
u_n = interpolate(u_0, V)
un0, un1 = u_n.split()
def boundary(x, on_boundary):
return on_boundary
n = FacetNormal(mesh)
u = Function(V)
v = TestFunction(V)
f = Expression(('(1/4)*pow(x[0], 2) + 9*pow(x[1], 2)'), degree=2)
t=0
F = u[0]*v[1]*dx - u_n[0]*v[1]*dx - dt/2*(dot(grad(u[1]), grad(v[1])))*dx - (dt/2)*f*u[1]*v[1]*dx - dt/2*(dot(grad(u_n[1]), grad(v[1])))*dx -(dt/2)*f*u_n[1]*v[1]*dx + u[1]*v[0]*dx + (dt/2)*(dot(grad(u[0]), grad(v[0])))*dx + (dt/2)*(u[0]*f*v[0])*dx - u_n[1]*v[0]*dx + (dt/2)*(dot(grad(u_n[0]), grad(v[0])))*dx + (dt/2)*f*u_n[0]*v[0]*dx - 5*dt*(u[0]**2 + u_n[0]**2 + u[1]**2 + u_n[1]**2)*(u[1]+u_n[1])*v[1]*dx + 5*dt*(u[0]**2 + u_n[0]**2 + u[1]**2 + u_n[1]**2)*(u[0]+u_n[0])*v[0]*dx
vtkfile_8 = File('GFlow/solution.pvd')
J = derivative(F, u)
for n in range(num_steps):
# Update current time
t += dt
# Compute solution
solve(F == 0, u, J=J)
u0, u1 = u.split()
vtkfile_8 << (u, t)
# Update previous solution
u_n.assign(u)
u0, u1 = u.split()
plot(u0**2+u1**2)
pl = plot(u0**2+u1**2)
plt.savefig('GFlowfigs/unormsq_'+ str(n)+'.png')
In the GIF below the density is plotted.
