Simple Runge-Kutta Scheme for nonlinear PDE

Question

I am new to this community as well as to scientific programming. I programmed a simple 4th order Runge-Kutta for the 1-D Cahn-Hilliard Equation for some first simple calculations on pattern forming systems. It turned out to be extremely sensitive to the change of the size of the spatial grid. With 100 grid elements the time step needs to be like 100 times smaller in order to converge ... This is the code i used:

import numpy as np
import matplotlib.pyplot as plt

# Parameters

k_c = 1.
eps = 0.5
L=12. # Domainsize

size = 50 # Gridsize
dx = L/size

T = 100 # total time
dt = .001 # time step
n = int(T/dt)

# Random Initial Conditions

U = np.random.normal(0,0.01,size) # Normalverteilte Zufallswerte um 0

# Laplace-Operator

def laplacian(Z):
    Zleft = np.roll(Z,1)
    Zright = np.roll(Z,-1)
    Zcenter = Z
    return (Zright + Zleft -2*Zcenter)/dx**2

def DoubleLaplacian(Z):
    Z2left = np.roll(Z,2)
    Z2right = np.roll(Z,-2)
    Zleft = np.roll(Z,1)
    Zright = np.roll(Z,-1)
    Zcenter = Z
    return (Z2left - 4* Zleft + 6*Zcenter - 4*Zright + Z2right)/dx**4

# Equation Body

def equation(U):
     return laplacian(U**3)-2*eps/(k_c**2)*laplacian(U)-eps*(k_c**4)*DoubleLaplacian(U)

# Runge-Kutta-Step

def RK4Step(U, dt, equation):
    k1 = dt*equation(U)
    U_temp = U + k1/2; k2 = dt*equation(U_temp);
    U_temp = U + k2/2; k3 = dt*equation(U_temp);
    U_temp = U + k3; k4 = dt*equation(U_temp);
    U += (k1 + 2*k2 + 2*k3 + k4)/6

# Integration
for i in range(n):
    RK4Step(U, dt, equation)

Since I am new to this kind of stuff I wanted to ask whether there is a conceptual fault in my thinking or coding or it is simply due to the nonlinearity. I would be as well pretty happy to get some hints how to improve the performance of the code.

Answer 1

Profiling your code

import cProfile, pstats, StringIO
pr = cProfile.Profile()
pr.enable()
# Integration
for i in range(n):
    RK4Step(U, dt, equation)
pr.disable()
s = StringIO.StringIO()
sortby = 'cumulative'
ps = pstats.Stats(pr, stream=s).sort_stats(sortby)
ps.print_stats()
print s.getvalue()

gives the following :

         27300018 function calls in 103.787 seconds

   Ordered by: cumulative time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        2    0.000    0.000  103.787   51.894 /Users/Leo/anaconda/lib/python2.7/site-packages/IPython/core/interactiveshell.py:3005(run_code)
        1    0.466    0.466  103.787  103.787 <ipython-input-2-3cf9531a6a7d>:5(<module>)
   100000    4.605    0.000  103.317    0.001 <ipython-input-1-3a5916441c4c>:44(RK4Step)
   400000    9.648    0.000   98.712    0.000 <ipython-input-1-3a5916441c4c>:39(equation)
  3200000   16.343    0.000   63.233    0.000 /Users/Leo/anaconda/lib/python2.7/site-packages/numpy/core/numeric.py:1279(roll)
   800000   13.794    0.000   45.859    0.000 <ipython-input-1-3a5916441c4c>:23(laplacian)
   400000   12.037    0.000   43.205    0.000 <ipython-input-1-3a5916441c4c>:29(DoubleLaplacian)
  3200000   12.589    0.000   12.589    0.000 {numpy.core.multiarray.concatenate}
  6400000   12.371    0.000   12.371    0.000 {numpy.core.multiarray.arange}

You can see how Laplacian and Double laplacian take about 88% of the computing time.

You can either rewrite your code in Fortran or C and wrap it with f2py or rewrite those functions in cython.

Another way is compiling your code using Nuitka, this decreases runtime to 64%.

You also see that roll is what takes most of the time. So you can start working on an optimised version of that.

         25600000 function calls in 64.852 seconds

   Ordered by: cumulative time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
  3200000   17.494    0.000   64.852    0.000 /usr/local/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/numpy/core/numeric.py:1279(roll)
  6400000   12.644    0.000   12.644    0.000 {numpy.core.multiarray.arange}
  3200000   12.635    0.000   12.635    0.000 {numpy.core.multiarray.concatenate}
  3200000    4.499    0.000   10.167    0.000 /usr/local/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/numpy/core/numeric.py:464(asanyarray)
  3200000    6.472    0.000    6.472    0.000 {method 'take' of 'numpy.ndarray' objects}
  3200000    5.668    0.000    5.668    0.000 {numpy.core.multiarray.array}
  3200000    5.440    0.000    5.440    0.000 {method 'reshape' of 'numpy.ndarray' objects}

Finally, find a way to get your CFL number to define your $\Delta t$