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 < DOLFIN_EPS and x > -DOLFIN_EPS and on_boundary) # Map right boundary (H) to left boundary (G) def map(self, x, y): y = x - 1.0 y = x # 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)', 'sin(x)'), 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, 2) + 9*pow(x, 2)'), degree=2) t=0 F = u*v*dx - u_n*v*dx - dt/2*(dot(grad(u), grad(v)))*dx - (dt/2)*f*u*v*dx - dt/2*(dot(grad(u_n), grad(v)))*dx -(dt/2)*f*u_n*v*dx + u*v*dx + (dt/2)*(dot(grad(u), grad(v)))*dx + (dt/2)*(u*f*v)*dx - u_n*v*dx + (dt/2)*(dot(grad(u_n), grad(v)))*dx + (dt/2)*f*u_n*v*dx - 5*dt*(u**2 + u_n**2 + u**2 + u_n**2)*(u+u_n)*v*dx + 5*dt*(u**2 + u_n**2 + u**2 + u_n**2)*(u+u_n)*v*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.