Code Basics

Code Basics

Basic Notation

The code solves partial differential equations of the general form:

\[\partial_t u = \mathcal{L}u + \mathcal{N}(u)\ .\]

(Note: ODEs are special cases of the above. Thus the code also solves ODEs.)

We decompose the right hand side of the above in a linear part ($\mathcal{L}u$) and a nonlinear part ($\mathcal{N}(u)$). The time steppers treat the linear and nonlinear parts differently.

Boundary conditions in all spatial dimensions are periodic. That allows us to expand all variables using a Fourier decomposition. For example, a variable $\phi(x, t)$ that depends in one spatial dimension is expanded as:

\[\phi(x, t) = \sum_{k} \widehat{\phi}(k, t)\,e^{\mathrm{i} k x}\ ,\]

where wavenumbers $k$ take the values $\tfrac{2\pi}{L_x}[0,\pm 1,\pm 2,\dots]$. The equation is time-stepped forward in Fourier space. That way $u$ becomes the array with all Fourier coefficients of the solution.

The coefficients for the linear operator $\mathcal{L}$ are stored in an array called LC. The term $\mathcal{N}(u)$ is computed for by calling the function calcN!.

Abstract SuperTypes

The code is divided along conceptual lines into problem-agnostic and problem-specific components. Files that contain problem-agnostic parts of the code are stored in /src. Files in /src define the domain, 'AbstractTypes' that supertype problem-specific types, and time-stepper types and routines. Problem-specific modules are stored in /src/physics.

Below is a list of all Abstract Supertypes used by the code:

Grids and time-steppers are generic and work for any problem of that particular dimension. State and Problem just gathers things together. Thus, to write a solver for a new physical problem you only need to prescribe params, vars, the coefficients of the linear part, LC, and function calcN!.

Source code organization

The code is divided along conceptual lines into problem-agnostic and problem-specific components. Files that contain problem-agnostic parts of the code are stored in /src. Files in /src define the domain, 'AbstractTypes' that supertype problem-specific types, and time-stepper types and routines. Problem-specific modules are stores in /src/physics.

Here's an overview of the code structure:

Basic steps for solving a problem: step through an example script

To illustrate the basic steps for solving a problem consider the 1D Kuramoto-Sivashinsky equation for $u(x, t)$:

\[\partial_t u + \partial_x^4 u + \partial_x^2 u + u\partial_x u = 0\ ,\]

which in Fourier base reads:

\[\partial_t \widehat{u} = \underbrace{(- k_x^4 + k_x^2) \widehat{u}}_{\mathcal{L}\widehat{u}} + \underbrace{\widehat{ -u\partial_x u }}_{\mathcal{N}(\widehat{u})}\ .\]

The steps to construct an AbstractProblem for the above are:

The example script found in examples/kuramotosivashinsky/trefethenexample.jl demonstrates the above steps needed to construct an AbstractProblem. The prob is constructed by calling prob = InitialValueProblem(nx=nx, Lx=Lx, dt=dt, stepper="ETDRK4"). Looking into the InitialValueProblem function we can see the above steps:

function InitialValueProblem(;
     nx = 256,
     Lx = 2π,
     dt = 0.01,
stepper = "RK4"
)

g  = OneDGrid(nx, Lx)
pr = Params()
vs = Vars(g)
eq = Equation(pr, g)
ts = FourierFlows.autoconstructtimestepper(stepper, dt, eq.LC, g)

FourierFlows.Problem(g, vs, pr, eq, ts)
end

The OneDGrid function is called for the grid. Within grid the wavenumber array is constructed:

i1 = 0:Int(nx/2)
i2 = Int(-nx/2+1):-1
k = Array{T}(2π/Lx*cat(1, i1, i2))
kr = Array{T}(2π/Lx*cat(1, i1))

For real-valued fields we use rfft and thus only positive wavenumbers are involved: array kr. E.g., for nx=8 and Lx=2π the wavenumber grids are: k = [0, 1, 2, 3, 4, -3, -2, -1] and kr = [0, 1, 2, 3, 4].

The construction of the grids only works for even number of grid points. Moreover, since the code relies on the $\mathrm{FFT}$ algorithm, we suggest you use a power of 2 as the number of grid points, since then $\mathrm{FFT}$ is most efficient.

Function Vars(g) initialize variables u, ux, and uux as real valued arrays of length nx and variables uh, uxh, and uuxh as complex valued arrays of length nkr = Int(nx/2+1) (the same length as kr). As a general convention variable names with h denote the Fourier transforms of the corresponding variable (h stands for 'hat').

The array LC is constructed by Equation function

function Equation(p, g)
  LC = @. g.kr^2 - g.kr^4
  FourierFlows.Equation(LC, calcN!)
end

Also eq includes function calcN! which computes the nonlinear term $\mathcal{N}(\widehat{u})$:

function calcN!(N, sol, t, s, v, p, g)
  @. v.uh = sol
  @. v.uxh = im*g.kr*sol
  A_mul_B!(v.u, g.irfftplan, v.uh)
  A_mul_B!(v.ux, g.irfftplan, v.uxh)
  @. v.uux = v.u*v.ux
  A_mul_B!(v.uuxh, g.rfftplan, v.uux)
  @. N = -v.uuxh
  dealias!(N, g)
  nothing
end

The time-stepper is constructed and stored as ts. Finally, all supertypes are gathered together as an AbstractProblem.

Tutorials