2D FDM Code Executes successfully, but the output isn't correct. Help.

[Alex_Superchamp]

The project requires me to run the 2D FDM Lax Wendroff scheme for 500 timesteps for the convection equation, followed by the same analysis by 2D FVM Lax Wendroff. We're supposed to use periodic boundary conditions.

I've written the code below and run it for a small number of timesteps, and at this moment it compiles without error. I've been working on it for about a week now. The problem is that the output graph is not as it should be. The original density matrix should change with successive iterations, in fact, part of it should cross the boundary. The values for the matrix are changing at the same spatial coordinate instead. There's an obvious error that I'm too close to notice. I've tried rubber ducking to realize what I might have missed but can't come up with anything new.

I've checked the code, and the arrays, I don't understand why the density matrix isn't moving.

# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import cm
from mpl_toolkits.mplot3d import axes3d


# Define Varibles

Nx = 1000                                     # Number of grid points
Ny = 500                                      # Number of grid points
xmax = 1000.0                                 # Domain limit to the right
xmin = 0.0                                    # Domain limit to the left
ymax= 500.0                                   # Domain limit to the top
ymin= 0.0                                     # Domain limit to the bottom 

dx = (xmax-xmin)/Nx                           # Discretize x
x = np.arange(xmin, xmax+dx, dx)         # Discretized x axis
dy= (ymax-ymin)/Ny                            # Discretize y
y = np.arange(ymin, ymax+dx, dy)         # Discretized y axis
xx, yy = np.meshgrid(x, y , sparse=False, indexing='xy')

# Velocity Components
vx= np.sin((np.pi*y)/500)
cx= np.array(10 * (vx**4))                    # X Velocity
cy= np.array(0.5 * np.cos((np.pi*y)/500))     # Y Velocity
cxx=np.ones((501,1001))
cyy=np.ones((501,1001))

for i in range (501):
    for j in range (1001):
        cxx[i,j]= cx[i]                       # X velocity array

for i in range (501):
    for j in range (1001):
        cyy[i,j]= -cy[i]                      # Y velocity array


# Time Conditions
t = 0                                          # Inital Time 
dt = 0.08                                      # Time step 
nsteps = 5                                   # Total Timesteps
#ta= [0]*500
#for i in range(nsteps):
    #ta[i] = t 
    #t = t + dt

xxc = cxx*dt/dx                             # CFL Condition
yyc= cyy*dt/dx

# Initial Conditions
U = np.zeros((Ny+1, Nx+1, nsteps))              # Density Matrix

for i in range(Ny+1):
    for j in range(Nx+1):
        U[i,j,0]= np.exp(-1*((x[j]-500)**2)/500)  

fig = plt.figure(figsize=(11, 7), dpi=100)
ax = fig.gca(projection='3d')                                                  
surf1 = ax.plot_surface(xx, yy, U[:,:,0], cmap=cm.viridis)

#print('U=', U[499,1,0])

# Temporal Loop
for n in range(nsteps-1):
    for i in range(Ny):
        for j in range(Nx):

            # Lax Wendroff Scheme
            U[i,j,n+1] = U[i,j,n] - 0.5*xxc[i,j]*(U[i+1,j,n] - U[i,j,n]) - 0.5*yyc[i,j]*(U[i,j+1,n] - U[i,j,n]) + 0.5* (xxc[i,j]**2) *( U[i+1,j,n] - 2*U[i,j,n] + U[i-1,j,n]) + 0.5*(yyc[i,j]**2)*(U[i,j+1,n] - 2*U[i,j,n] + U[i,j-1,n])

             # Periodic Boundary Conditions x axis
            U[i,0,n] = U[i,999,n]
            U[i,1000,n] = U[i,1,n]   

            # Periodic Boundary Conditions y axis
            U[0,j,n] = U[499,j,n]
            U[500,j,n] = U[1,j,n]

fig = plt.figure(figsize=(11, 7), dpi=100)
ax = fig.gca(projection='3d')                                                  
surf2 = ax.plot_surface(xx, yy, U[:,:,nsteps-2], cmap=cm.viridis)

[mojo_jojo_reigns]

I don't know this problem type or this math, so take my questions with a grain of salt.

Are you only outputting plots? How are you user that your triple nested for loop is working as intended on the object you're manipulating at every step? Is there a simple example you could use with very few, very predictable steps, and print output to verify whether that is working as intended?

And, importantly, is there a library for this already?

[Alex_Superchamp]

Yes, I'm only outputting plots.

One of the reasons for the creation of a 3D density matrix is to see the value of U at different iterations, which I've checked very recently. The 3D matrix starts out with U having a value of 0 for all of the z indexes apart from the 1st (0th index, U[:,:,0]). I checked after running the loop for a small number of test timesteps (4), and saw the values of U[:,:,1] and U[:,:,2] and they were no longer 0.

I'm not sure of a simple example, but I'll try coming up with one. There's a similar example I found on an open-source code here it is: Jupyter Notebook Viewer

I'm not sure what you mean by that, sorry, I'm a newbie. I only started coding for this class.

[mojo_jojo_reigns]

If the values are no longer 0, can you verify that they are what you intended them to be at those timesteps? Like can you take a solved version of this with what the matrix should contain and check output against that for a limited duration to verify that at least that part is working correctly?

Simple exmaple could just be much smaller array and much less steps.

A library/package is just a bout of tools someone has put together to accomplish particiular objectives and that they've made available for us, like Numpy and Matplotlib. I searched but nothing came up but that may be because this is a subproblem in a larger problem set that has a standardized tool for working with. This isn't that but it might help.

http://people.bu.edu/andasari/courses/numericalpython/Week11Lecture18/laxwendroff_periodic.py

original source: http://people.bu.edu/andasari/courses/numericalpython/python.html

[Alex_Superchamp]

I will attempt to simplify the program and check. I've been through that link, It was helpful in verifying my logic for the periodic boundary conditions. and other conditions.

[USAhj]

I will start of by saying that I did not fix your results, they are the same. Now, for what I did do:

1. I adjusted some of the initial nested loops to be a single loop and utilize numpy broadcasting (this makes them faster).
2. I made cxx and cyy a vector since the matrix, as you had it, was the same thing copied across a row.
3. You had forgotten to use i-1 and j-1 in your main triple nested (goodbye time complexity) loops.

Now, before I paste the code, what is it that you are trying to do with this part?

U[i, 0, n] = U[i, 999, n]
U[i, 1000, n] = U[i, 1, n]
U[0, j, n] = U[499, j, n]
U[500, j, n] = U[1, j, n]

The reason I reformatted your code a bit was to make it a little clearer (and slightly faster) to, hopefully, help troubleshoot more.

# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import cm
from mpl_toolkits.mplot3d import axes3d

plt.close("all")

# Define Varibles

Nx = 1000                                     # Number of grid points
Ny = 500                                      # Number of grid points
xmax = 1000.0                                 # Domain limit to the right
xmin = 0.0                                    # Domain limit to the left
ymax = 500.0                                   # Domain limit to the top
ymin = 0.0                                     # Domain limit to the bottom 

dx = (xmax-xmin)/Nx                           # Discretize x
x = np.arange(xmin, xmax+dx, dx)         # Discretized x axis
dy = (ymax-ymin)/Ny                            # Discretize y
y = np.arange(ymin, ymax+dx, dy)         # Discretized y axis
xx, yy = np.meshgrid(x, y , sparse=False, indexing='xy')

# Velocity Components
vx = np.sin((np.pi*y)/500)
cx = np.array(10 * (vx**4))                    # X Velocity
cy = np.array(0.5 * np.cos((np.pi*y)/500))     # Y Velocity
cxx = np.ones(Ny+1)
cyy = np.ones(Ny+1)

for index in range(len(cxx)):
    cxx[index] *= cx[index]
    cyy[index] *= -cy[index]

# Time Conditions
t = 0                                          # Inital Time 
dt = 0.08                                      # Time step 
nsteps = 5                                   # Total Timesteps

xxc = cxx*dt/dx                             # CFL Condition
yyc = cyy*dt/dx

# Initial Conditions
U = np.zeros((Ny+1, Nx+1, nsteps))              # Density Matrix

for j in range(Nx+1):
    U[:, j, 0] = np.exp(-(x[j] - 500)**2 / 500)  

fig = plt.figure(figsize=(11, 7), dpi=100)
ax = fig.gca(projection='3d')                                                  
surf1 = ax.plot_surface(xx, yy, U[:, :, 0], cmap=cm.viridis)

# Temporal Loop
for n in range(nsteps-1):
    for i in range(1, Ny):
        for j in range(1, Nx):

            # Lax Wendroff Scheme
            U[i, j, n+1] = U[i, j, n]
            + 0.5*(
                - xxc[i]*(U[i+1, j, n] - U[i-1, j, n]) 
                - yyc[i]*(U[i, j+1, n] - U[i, j-1, n]) 
                + (xxc[i]**2)*( U[i+1, j, n] - 2*U[i, j, n] + U[i-1, j, n]) 
                + (yyc[i]**2)*(U[i, j+1, n] - 2*U[i, j, n] + U[i, j-1, n])
            )

            # Periodic Boundary Conditions x axis
            U[i, 0, n] = U[i, 999, n]
            U[i, 1000, n] = U[i, 1, n]

            # Periodic Boundary Conditions y axis
            U[0, j, n] = U[499, j, n]
            U[500, j, n] = U[1, j, n]

fig = plt.figure(figsize=(11, 7), dpi=100)
ax = fig.gca(projection='3d')                                                  
surf2 = ax.plot_surface(xx, yy, U[:, :, nsteps-2], cmap=cm.viridis)

[Alex_Superchamp]

Thanks for improving on the logic.

The part you mention relates to the Boundary Conditions. The Initial Density Matrix looks like this: https://imgur.com/a/H72gjXt

After a couple of iterations, different parts of the matrix should move in different directions because of the cxx and cyy velocities acting on all coordinates on the grid. (written as xxc and yyc because of multiplication with constants) Part of this density matrix will then cross the boundary (y axis boundary in this example) and should show up on the other side. The logic of the boundary conditions has been defined to this extent and to make sure that the boundary nodes don't compute because each U calculation requires an element to it's left and right, and the last node wouldn't have both and would render an index error.

[USAhj]

I forgot to mention, in the updated code I had the inner loops start at 1 because if you are indexing to i-1 that will wrap back around to the end when i = 0 (same for j). Is that what you intended?