# Tutorial 3: Problems and Solvers

**Overview:** This tutorials covers the basics of setting up and solving problems using Dedalus. Dedalus solves symbolically specified initial value, boundary value, and eigenvalue problems.

```
[1]:
```

```
import numpy as np
import matplotlib.pyplot as plt
from dedalus import public as de
```

```
[2]:
```

```
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
```

## 3.1: Problems

### Problem formulations

Dedalus standardizes the formulation of all initial value problems by taking systems of symbolically specified equations boundary conditions and manipulating them into the following generic form:

where \(\mathcal{M}\) and \(\mathcal{L}\) are matrices of linear differential operators, \(\mathcal{X}\) is a state vector of the unknown fields, \(\mathcal{F}\) is a vector of general nonlinear expressions. This form encapsulates prognostic/evolution equations (equations containing time derivatives) and diagnostic/algebraic constraints (equations without time derivatives), both in the interior of the domain and on the boundary.

Importantly, the left-hand side (LHS) of the equations must be time-independent, first-order in temporal derivatives, and linear in the problem variables. The right-hand side (RHS) of the equations may contain nonlinear and time-dependent terms, but no temporal derivatives. In addition, the LHS must be first-order in non-separable (e.g. Chebyshev) derivatives, and may only contain non-constant coefficients that vary in the non-periodic dimensions (others must be placed on the RHS).

In addition to initial value problems, there are similar standard forms for generalized eigenvalue problems (\(\sigma \mathcal{M} \cdot \mathcal{X} + \mathcal{L} \cdot \mathcal{X} = 0\)), linear boundary value problems (\(\mathcal{L} \cdot \mathcal{X} = \mathcal{G}\)), and nonlinear boundary value problems (\(\mathcal{L} \cdot \mathcal{X} = \mathcal{F}(\mathcal{X})\)). These four problem types are represented by the `IVP`

, `EVP`

, `LBVP`

, and `NLBVP`

problem classes,
respectively.

To create a problem object, you must provide a domain object and the names of the variables that appear in the equations. Let’s start setting up the complex Ginzburg-Landau equation (CGLE) for a variable \(u(x,t)\) on a finite interval \(x \in [0, 300]\) with homogeneous Dirichlet boundary conditions:

The CGLE only has one primative variable but is second-order in its spatial derivatives, so we’ll have to introduce an extra variable to reduce the equation to first-order. We also pick a dealiasing factor of 2 to correctly dealias the cubic nonlinearity.

```
[3]:
```

```
# Build bases and domain
x_basis = de.Chebyshev('x', 1024, interval=(0, 300), dealias=2)
domain = de.Domain([x_basis], grid_dtype=np.complex128)
# Build problem
problem = de.IVP(domain, variables=['u', 'ux'])
```

### Variable metadata

Metadata for the problem variables can be specified through the `meta`

attribute of the problem object, and indexing by variable name, axis, and property, respectively.

The most common metadata to set here is the `parity`

option when using the `SinCos`

basis. The parity for each field and for each `SinCos`

basis needs to be set to either `+1`

or `-1`

to indicate that the field has even (cosine) or odd (sine) parity around the interval endpoints for the corresponding dimension.

Here we aren’t using a `SinCos`

basis, but if we were, we could set the variable parities as:

```
[4]:
```

```
# problem.meta['u']['x']['parity'] = +1
# problem.meta['ux']['x']['parity'] = -1
```

### Parameters and non-constant coefficients

Before adding equations to the problem, we first add any parameters, defined as fields or scalars used in the equations but not part of the state vector of problem variables, to the `problem.parameters`

dictionary.

For constant/scalar parameters, like we have here, we simply add the desired numerical values to the parameters dictionary.

```
[5]:
```

```
problem.parameters['b'] = 0.5
problem.parameters['c'] = -1.76
```

For non-constant coefficients (NCCs), we pass a field object with the desired data. Dedalus only accepts NCCs that couple the non-periodic dimensions, i.e. are constant along all `Fourier`

and `SinCos`

dimensions, so that those dimensions remain linearly uncoupled. To inform the parser that a NCC will not couple these directions, you must explicitly add some metadata to the NCC fields indiciating that they are constant along any periodic directions.

We don’t have NCCs or periodic dimensions here, but we’ll sketch the process here anyways. Consider a 3D problem on a Fourier (x), SinCos (y), and Chebyshev (z) domain. Here’s how we would add a simple non-constant coefficient in z to a problem:

```
[6]:
```

```
# ncc = Field(domain, name='c')
# ncc['g'] = z**2
# ncc.meta['x', 'y']['constant'] = True
# problem.parameters['c'] = ncc
```

### Substitutions

To simplify equation entry, you can define substitutions which effectively act as string-replacement rules during the parsing process. Substitutions can be used to provide short aliases to quantities computed from the problem variables. They can also define shortcut functions similar to Python lambda functions, but with normal mathematical syntax.

Here we’ll create a simple substitution for computing the squared magnitude of a variable:

```
[7]:
```

```
# Function-like substitution using dummy variables
problem.substitutions["mag_sq(A)"] = "A * conj(A)"
```

### Equation entry

Equations and boundary conditions are then entered in plain text, optionally with condition strings specifying which separable modes (indexed by `nx`

and `ny`

for separable axes named `x`

and `y`

, etc.) that equation applies to.

The parsing namespace basically consists of:

The variables, parameters, and substitutions defined in the problem

The axis names representing the individual basis grids, e.g.

`'x'`

The derivative operators for each basis, named as e.g.

`'dx'`

The

`differentiate`

,`integrate`

, and`interpolate`

factories as`'d'`

,`'integ'`

, and`'interp'`

`'left'`

and`'right'`

as aliases to interpolation at the Chebyshev endpoints, if presentTime and temporal derivatives as

`'t'`

and`'dt'`

Simple mathematical functions (logarithmic and trigonometric), e.g.

`'sin'`

,`'exp'`

, …

Let’s enter the CGLE (as two first-order equations) and the boundary conditions:

```
[8]:
```

```
# Add main equation, with linear terms on the LHS and nonlinear terms on the RHS
problem.add_equation("dt(u) - u - (1 + 1j*b)*dx(ux) = - (1 + 1j*c) * mag_sq(u) * u")
# Add auxiliary equation defining the first-order reduction
problem.add_equation("ux - dx(u) = 0")
# Add boundary conditions
problem.add_equation("left(u) = 0")
problem.add_equation("right(u) = 0")
```

## 3.2: Solvers

### Building a solver

Each problem type (IVP, EVP, LBVP, and NLBVP) has a corresponding solver class that actually performs the solution steps for a corresponding problem. Solvers are simply built using the `problem.build_solver`

method.

For IVPs, we select a timestepping method when building the solver. Several multistep and Runge-Kutta IMEX schemes are available in the timesteppers.py module, and can be selected by name.

```
[9]:
```

```
# Build solver
solver = problem.build_solver('RK222')
```

```
2022-07-24 12:55:03,925 pencil 0/1 INFO :: Building pencil matrix 1/1 (~100%) Elapsed: 0s, Remaining: 0s, Rate: 5.5e+01/s
```

### Setting stop criteria

For IVPs, stopping criteria for halting time evolution are specified by setting the `solver.stop_iteration`

, `solver.stop_wall_time`

(seconds since solver instantiation), and/or `solver.stop_sim_time`

attributes.

Let’s stop at \(t = 500\) in simulation units:

```
[10]:
```

```
# Stopping criteria
solver.stop_sim_time = 500
solver.stop_wall_time = np.inf
solver.stop_iteration = np.inf
```

### Setting initial conditions

The fields representing the problem variables can be accessed with a dictionary-like interface through the `solver.state`

system. For IVPs and nonlinear BVPs, initial conditions are set by directly modifying the state variable data before running a simulation.

```
[11]:
```

```
# Reference local grid and state fields
x = domain.grid(0)
u = solver.state['u']
ux = solver.state['ux']
# Setup a sine wave
u.set_scales(1)
u['g'] = 1e-3 * np.sin(5 * np.pi * x / 300)
u.differentiate('x', out=ux);
```

### Solving/iterating a problem

IVPs are iterated using the `solver.step`

method with a provided timestep. EVPs are solved using the `solver.solve_dense`

or `solver.solve_sparse`

methods. LBVPs are solved using the `solver.solve`

method. NLBVPs are iterated using the `solver.newton_iteration`

method.

The logic controlling the main-loop of a Dedalus IVP simulation occurs explicitly in the simulation script. The `solver.proceed`

property can change from `True`

to `False`

once any of the specified stopping criteria have been met. Let’s timestep our problem until the halting condition is reached, copying the grid values of `u`

every few iterations. This should take just a few seconds to run.

```
[12]:
```

```
# Setup storage
u.set_scales(1)
u_list = [np.copy(u['g'])]
t_list = [solver.sim_time]
# Main loop
dt = 0.05
while solver.ok:
solver.step(dt)
if solver.iteration % 10 == 0:
u.set_scales(1)
u_list.append(np.copy(u['g']))
t_list.append(solver.sim_time)
if solver.iteration % 1000 == 0:
print('Completed iteration {}'.format(solver.iteration))
# Convert storage lists to arrays
u_array = np.array(u_list)
t_array = np.array(t_list)
```

```
Completed iteration 1000
Completed iteration 2000
Completed iteration 3000
Completed iteration 4000
Completed iteration 5000
Completed iteration 6000
Completed iteration 7000
Completed iteration 8000
Completed iteration 9000
Completed iteration 10000
2022-07-24 12:55:10,822 solvers 0/1 INFO :: Simulation stop time reached.
```

Now let’s make a space-time plot of the magnitude of the solution:

```
[13]:
```

```
# Plot solution
plt.figure(figsize=(6, 7), dpi=100)
plt.pcolormesh(x, t_array, np.abs(u_array), shading='nearest')
plt.colorbar()
plt.xlabel('x')
plt.ylabel('t')
plt.title('Hole-defect chaos in the CGLE: |u|')
plt.tight_layout()
```