ABC-flow integration and visualization

Question

I am working my way through the book "A First Course in Computational Fluid Dynamics". One of the exercises is to integrate the ABC-flow equations using Runge-Kutta and visualize it.

Here is my naive implementation in Python.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

This implements the ABC-flow computation:

def abc_flow( t, p ):
    x = p[ 0 ]
    y = p[ 1 ]
    z = p[ 2 ]

    u = A * np.sin( z ) + C * np.cos( y )
    v = B * np.sin( x ) + A * np.cos( z )
    w = C * np.sin( y ) + B * np.cos( x )

    return np.array( [ u, v, w ] )

This is the 4th order Runge-Kutta implementation (I know I could use a ready-made solver, this is just for learning purposes):

def runge_kutta( f, p, t ):
    k0 = f( t, p )
    k1 = f( t + dt * 0.5, p + ( dt * 0.5 ) * k0 )
    k2 = f( t + dt * 0.5, p + ( dt * 0.5 ) * k1 )
    k3 = f( t + dt, p + dt * k2 )

    p_next = p + ( dt / 6.0 ) * ( k0 + 2 * k1 + 2 * k2 + k3 )

    return p_next

These are the parameters as suggested in the book:

L = 2 * np.pi

A = 1
B = 1
C = 0.65

Finally, this is the integration loop:

N = 100000
dt = 0.01

if __name__ == '__main__':
    p = np.array( [ 3, np.pi, np.pi ] )

    x = [ p[ 0 ] ]
    y = [ p[ 1 ] ]
    z = [ p[ 2 ] ]

    for _ in range( N ):
        p = runge_kutta( abc_flow, p, 0 )

        x.append( p[ 0 ] )
        y.append( p[ 1 ] )
        z.append( p[ 2 ] )

The book asks to plot the particle position from each side of the box, which I have implemented as follows:

    fig = plt.figure( tight_layout=True )
    gs = gridspec.GridSpec( 3, 1 )

    ax = fig.add_subplot( gs[0, 0] )
    ax.plot( x, y )

    ax = fig.add_subplot( gs[1, 0] )
    ax.plot( x, z )

    ax = fig.add_subplot( gs[2, 0] )
    ax.plot( y, z )

    plt.show()

This is the image I get:

I think I am either missing a simple bug or my understanding of the problem is incorrect.

EDIT: adding the expected results: and problem statement:

Answer 1

Your solution looks correct. Note that the reference pictures display the solutions modulo $2\pi$, as also the equations are invariant under shifts by $2\pi$. In plotting these reduced solutions one has to make sure to cut out the jumps at the wrap-arounds.

See https://scicomp.stackexchange.com/questions/29823/poincare-map-for-arnold-beltrami-childress-magnetic-field-in-python, https://stackoverflow.com/questions/53792164/how-to-implement-a-method-to-generate-poincar%c3%a9-sections-for-a-non-linear-system for other discussions of this topic.

---

With the paramters and initial point used, you get one tunnel among the many stable looking components. Using 20 random initial points gives the projections

One can get "true 3D" with a stereoscopic image. There will be slight irritations as matplotlib does not implement sufficient 3D Z-ordering.

and finally the section with the plane $\{y=0\}$