For our Bachelors thesis we're trying to solve the Schrodinger equation $i\partial_tu = -\nabla^2u+Vu$ in FEniCS. Given the domain $[-5, 5]^2$ with an initial value of $u_0(x, y)=e^{(-2(x^2+y^2))}$ and with $V(x,y)=\frac{x^2}{4}+9y^2$ at time $T=2$ the solution $|u|^2$ should look like the following image given by our supervisor.
This is however the result we get:
The equation is dicretized with Crank-Nicolson as follows: (where $u=p+qi$)
For the real part:
$\frac{\Delta t}{2}(\nabla^2p^{n+1}-Vp^{n+1})-q^{n+1}=-q^n + \frac{\Delta t}{2}(Vp^n-\nabla^2p^n)$
for the imaginary part:
$\frac{\Delta t}{2}(\nabla^2q^{n+1}-Vq^{n+1})+p^{n+1}=-p^n + \frac{\Delta t}{2}(\nabla^2q^n-Vq^n)$
In the code below, $u[0]$ is $p$ and $u[1]$ is $q$
We're having trouble creating the correct residual when solving for $u$.
In FEniCS the equation is written in weak form and when multiplying with the test function we change sign of $dot(grad(), grad())$ to accommodate for multiplication with $i$. This is were we probably are at fault because we do not know how to construct the residual. When changing the sign of the gradient for the imaginary case the solution does not converge.
T = 2.0 # final time
num_steps = 200 # number of time steps
dt = T / num_steps # time step size
mesh = RectangleMesh(Point(-5, -5), Point(5, 5), 55, 55)
V = VectorFunctionSpace(mesh, 'Lagrange', 1, dim=2)
# Define initial value
u_0 = Expression( ( 'exp(-2*(pow(x[0], 2)+pow(x[1], 2)))', '0'),
degree=1)
u_n = interpolate(u_0, V)
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V,(0, 0), boundary)
u = Function(V)
v = TestFunction(V)
f = Expression(('(1/4)*pow(x[0], 2) + 9*pow(x[1], 2)', '0'), degree=2)
t=0
#Residual
ReVL = -u[1]*v[0]*dx + (dt/2)*(-dot(grad(u[0]), grad(v[0])) -
f[0]*u[0]*v[0])*dx
ReHL = -u_n[1]*v[0]*dx + (dt/2)*(f[0]*u_n[0]*v[0] + dot(grad(u_n[0]),
grad(v[0])))*dx
ImVL = u[0]*v[1]*dx + (dt/2)*(-dot(grad(u[1]), grad(v[1])) -
f[0]*u[1]*v[1])*dx
ImHL = u_n[0]*v[1]*dx + dt/2*(f[0]*u_n[1]*v[1] + dot(grad(u_n[1]),
grad(v[1])))*dx
FReal = ReVL - ReHL
FIm = ImVL - ImHL
Fny = FReal + FIm
J = derivative(Fny, u)
for n in range(num_steps):
# Update current time
t += dt
# Compute solution
solve(Fny == 0, u, bc, J=J)
# Update previous solution
u_n.assign(u)
We've also tried our luck doing as the example found in this link without success: Gross-Pitaevskii in FEniCS
What are we doing wrong here?
