0

I am trying to implement two numerical solutions. A forward Euler and a second order Runge-Kutta for the unsteady 2D heat equation with periodic boundary conditions. I am using a 3 point central difference in both cases for the spatial discretization.

This is my code:

def analytical(npx,npy,t,alpha=0.1): #Function to create analytical solution data
    X = numpy.linspace(0,1,npx,endpoint=False) #X spatial range
    Y = numpy.linspace(0,1,npy,endpoint=False) #Y Spatial Range
    uShape = (1,numpy.shape(X)[0],numpy.shape(Y)[0]) #Shape of data array
    u = numpy.zeros(uShape) #Allocate data array
    m = 2 #m and n = 2 makes the function periodic in 0->1 
    n = 2
    for i,x in enumerate(X): #Looping through x and y to produce analytical solution
        for j,y in enumerate(Y):
            u[0,i,j] = numpy.sin(m*pi*x)*numpy.sin(n*pi*y)*numpy.exp(-(m*m+n*n)*pi*pi*alpha*t)
    return u,X,Y

def numericalHeatComparisonFE(): #Numerical function for forward euler
    arraysize = 10 #Size of simulation array in x and y
    d = 0.1 #value of pde coefficient alpha*dt/dx/dx
    alpha = 1 #thermal diffusivity
    dx = 1/arraysize #getting spatial step
    dt = float(d*dx**2/alpha) #getting time step
    T,x,y = analytical(arraysize,arraysize,0,alpha=alpha) #get analytical solution
    numerical = numpy.zeros((2,)+T.shape) #create numerical data array
    ns = numerical.shape #shape of numerical array aliased
    numerical[0,:,:,:] = T[:,:,:] # assign initial conditions to first element in numerical
    error = [] #create empty error list for absolute error
    courant = alpha*dt/dx/dx #technically not the courant number but the coefficient for PDE
    for i in range(1,20):#looping through twenty times for testing - solving FE each step
        T,x,y = analytical(arraysize,arraysize,i*dt,alpha=alpha)
        for idx,idy in numpy.ndindex(ns[2:]):
            dxx = numerical[0,0,idx-1,idy]+numerical[0,0,(idx+1)%ns[-2],idy]-2*numerical[0,0,idx,idy] #X direction diffusion
            dyy = numerical[0,0,idx,idy-1]+numerical[0,0,idx,(idy+1)%ns[-1]]-2*numerical[0,0,idx,idy] #Y direction diffusion
            numerical[1,0,idx,idy] = courant*(dxx+dyy)+numerical[0,0,idx,idy] #Update formula
        error.append(numpy.amax(numpy.absolute(numerical[1,:,:,:]-T[:,:,:])))#Add max error to error list
        numerical[0,:,:,:] = numerical[1,:,:,:] #Update initial condition
    print(numpy.amax(error))

def numericalHeatComparisonRK2():
    arraysize = 10 #Size of simulation array in x and y
    d = 0.1 #value of pde coefficient alpha*dt/dx/dx
    alpha = 1 #thermal diffusivity
    dx = 1/arraysize #getting spatial step
    dt = float(d*dx**2/alpha) #getting time step
    T,x,y = analytical(arraysize,arraysize,0,alpha=alpha) #get analytical solution
    numerical = numpy.zeros((3,)+T.shape) #create numerical data array
    ns = numerical.shape #shape of numerical array aliased
    numerical[0,:,:,:] = T[:,:,:] # assign initial conditions to first element in numerical
    error = [] #create empty error list for absolute error
    courant = alpha*dt/dx/dx #technically not the courant number but the coefficient for PDE
    for i in range(1,20): #Test twenty time steps -RK2
        T,x,y = analytical(arraysize,arraysize,i*dt,alpha=alpha)
        for idx,idy in numpy.ndindex(ns[2:]): #Intermediate step looping through indices
            #Intermediate
            dxx = numerical[0,0,idx-1,idy]+numerical[0,0,(idx+1)%ns[-2],idy]-2*numerical[0,0,idx,idy]
            dyy = numerical[0,0,idx,idy-1]+numerical[0,0,idx,(idy+1)%ns[-1]]-2*numerical[0,0,idx,idy]
            numerical[1,0,idx,idy] = 0.5*courant*(dxx+dyy)+numerical[0,0,idx,idy]
            
        for idx,idy in numpy.ndindex(ns[2:]): #Update step looping through indices
            #RK Step
            dxx = numerical[1,0,idx-1,idy]+numerical[1,0,(idx+1)%ns[-2],idy]-2*numerical[1,0,idx,idy]
            dyy = numerical[1,0,idx,idy-1]+numerical[1,0,idx,(idy+1)%ns[-1]]-2*numerical[1,0,idx,idy]
            numerical[2,0,idx,idy] = courant*(dxx+dyy)+numerical[0,0,idx,idy]
        error.append(numpy.amax(numpy.absolute(numerical[2,:,:,:]-T[:,:,:]))) #Add maximum error to list
        numerical[0,:,:,:] = numerical[2,:,:,:] #Update initial conditions
    print(numpy.amax(error))

if __name__ == "__main__":
    numericalHeatComparisonFE()
    numericalHeatComparisonRK2()

when running the code, I expect that the maximum error for the RK2 should be less than that of the FE but I get 0.0021498590913591187 for the FE and 0.011325197051528346 for the RK2. I have searched the code pretty thoroughly and haven't found any glaring typos or errors. I feel it has to be something minor that I am missing but I can't seem to find it. If you happen to spot an error or know something I don't help or a comment would be appreciated.

Thanks!

sokato
  • 354
  • 4
  • 14
  • Before trying to see anything in this undocumented code: 1. Why do you expect the RK2 error to be less than the FE error? 2. Did you verify your methods by checking their orders? 3. What norms do you use to measure your errors? – Peter Meisrimel Jul 13 '20 at 08:12
  • @PeterMeisrimel well Runge-Kutta second order is second order in time and Forward Euler is first order in time so running both with the same step sizes I would expect RK2 to have less error or at least produce less error as the time step is decreased. The central difference is second order in space for both. I am not looking at error norms but absolute error from the analytical solution developed in this code. Yes, the code is undocumented but I felt it straightforward enough that it could be understood mainly from variable and function names. I can add comments though if that will help you. – sokato Jul 13 '20 at 20:03
  • 1
    Sure, these methods have first and second order, but you might want to numerically verify that your implementation also shows the expected orders. Order is a statement for stepsize going to zero, i.e., it only tells you something about how the error decreases for dt resp. dx to zero, but does not make a statement for any given dt or dx. For norms, it would be more common to use the discrete L2 norm, since you have a discrete function in space. With regards to documentation, one can follow the structure, but the devil is mostly in details and little indices. – Peter Meisrimel Jul 13 '20 at 20:13
  • @PeterMeisrimel I added some comments to hopefully clarify things, if there is more that I can add to aide understanding please let me know. I think verifying the order is a good suggestion as well as the L2 norm. I have seen that used pretty frequently in papers. I will give those a shot and see what they show me. Thank you for your suggestions. – sokato Jul 13 '20 at 20:23
  • So, if I understand this correctly now, you take a fixed dx and dt and compare with an analytical solution in time in space, thus measuring both the error in space and time. To get the error in time only, you need a sufficiently fine resolution in space. I.e., given dt, you get the time-integration error by taking increasing spatial resolutions, dx -> 0, and check where the error stagnates. This will be your time-integration error. – Peter Meisrimel Jul 14 '20 at 10:36
  • @PeterMeisrimel, I plotted the order of accuracy for both yesterday and got the RK2 to be 1st order on a loglog as dt->0 which isn't correct. I did this by fixing the coefficient d=0.1 and reducing the array size then solving for the appropriate dt to maintain d=0.1. As for you previous comment, I was taking a fixed dt and dx and comparing them directly though yes. – sokato Jul 14 '20 at 17:10
  • If you fix d and reduce arraysize (dx increasing) you do not get dt->0? The typical way to verify the order of your time-integration method is to first fix alpha and then dx, then you solve for decreasing dt. This ensures you are solving **the same problem** for different dt. Sidenote: I am not 100% sure how stability concerns work out for convergence test, but I guess this is another plus point for implicit methods, which I would recommend here. – Peter Meisrimel Jul 14 '20 at 21:24
  • @PeterMeisrimel, I mixed up my wording. I meant increase arraysize i.e. reducing dx. Fixing alpha and dx and reducing dt is probably better anyways so it shows how the error changes only with dt. Implicit methods are great in that regard but the solver I am working on specifically requires an explicit method for the other methods to work. This is a unit test within that solver. – sokato Jul 15 '20 at 17:43
  • Do you satisfy some CFL condition on your time and space step size ? – Thibault Cimic Aug 17 '20 at 10:13
  • @ThibaultCimic I was satisfying a CFL condition with nu=0.1 if I remember correctly but I decided to move away from RK2 anyway because it wasn't entirely necessary for my work. I just seems like it should work. Thanks for the comment though. – sokato Aug 31 '20 at 17:57

0 Answers0