Unable to plot squared norm of solution over domain

mesh
1

While trying to solve the time-dependent Schrodinger equation with a initial condition and dirichlet boundary conditions over a rectangular shaped domain, we’re trying to plot the values of |u|^2 on the domain. This is done with a VectorFunctionSpace were u[0] corresponds to the real part of the equation and u[1] corresponds to the imaginary part, as seen in the residual F below. We want to plot the u[0]^2+u[1]^2 in a contour plot.

from fenics import *
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt

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), 50, 50)
V = VectorFunctionSpace(mesh, 'Lagrange', 1, dim=2)
# Define initial value
u_0 = Expression(  ( 'exp(-2*(pow(x[0], 2)))+ exp(-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)

QQ = FunctionSpace(mesh, "R", 0, 2)

#Residual
F = u[0]*v[0]*dx - dt/2*dot(grad(u[1]),grad(v[0]))*dx - dt/2*f[0]*u[1]*v[0]*dx - u_n[0]*v[0]*dx - dt/2*dot(grad(u_n[1]),grad(v[0]))*dx - dt/2*f[0]*u_n[1]*v[0]*dx + u[1]*v[1]*dx + dt/2*dot(grad(u[0]),grad(v[1]))*dx + dt/2*f[0]*u[0]*v[1]*dx - u[1]*v[1]*dx + dt/2*dot(grad(u[0]), grad(v[1]))*dx + dt/2*f[0]*u[0]*v[1]*dx  
J = derivative(F, u)

t=0
for n in range(num_steps):  
    # Update current time
    t += dt
    # Compute solution
    solve(F == 0, u, bc, J=J)
    #print(norm_u)
    unormsq = project(inner(u,u),QQ) 
    # Update previous solution
    u_n.assign(u)

plot(unormsq)
plt.savefig('testunorm2.png')

the right contour plot is what we expect to get (given from our supervisor) and the left plot is what we get.

plots
plots1178×750 137 KB

Before we’re even able to conclude if the solution is correct we’re not able to plot the solution with a contour plot on the mesh. How would we go about plotting this data?

2

Maybe this is what you need:

u0, u1 = u.split()
plot(u0**2+u1**2)
plt.show()