# Step-1: A first spatial problem¶

In this step, we create a spatial problem and play with it a little bit. No SDC, no multi-level here so far.

## Part A: Spatial problem setup¶

We start as simple as possible by creating a first, very simple problem using one of the `problem_classes`, namely `HeatEquation_1D_FD`. This basically consists only of the matrix `A`, which represents a finite difference discretization of the 1D Laplacian. This is tested for one particular example.

Important things to note:

• Many (most) parameters for pySDC are passed using dictionaries

• Data types are encapsulated: the real values are stored in `values`, meta-information can be stored separately in the data structure

• A quick peak into `HeatEquation_1D_FD` reveals that the `init` and the `params.nvars` attribute contain the same values (namely `nvars`). Yet, sometimes one or the other is used here (and throughout the code). The reason is this: the data structure (`mesh` in this case) needs some form of standard initialization. For this, pySDC uses the `init` attribute each problem class has. In our simple example, this is he same as `nvars`, but e.g. for Finite Elements, the function space is stored in the `init` attribute. Thus, whenever a data type is created, use `init`, otherwise do whatever looks nice.

• We happily pass classes around so that they can be used to instantiate things within subroutines

```import numpy as np
from pathlib import Path

from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_unforced

def main():
"""
A simple test program to set up a spatial problem and play with it
"""
# instantiate problem
prob = heatNd_unforced(
nvars=1023,  # number of degrees of freedom
nu=0.1,  # diffusion coefficient
freq=4,  # frequency for the test value
bc='dirichlet-zero',  # boundary conditions
)

# run accuracy test, get error back
err = run_accuracy_check(prob)

Path("data").mkdir(parents=True, exist_ok=True)
f = open('data/step_1_A_out.txt', 'w')
out = 'Error of the spatial accuracy test: %8.6e' % err
f.write(out)
print(out)
f.close()

assert err <= 2e-04, "ERROR: the spatial accuracy is higher than expected, got %s" % err

def run_accuracy_check(prob):
"""
Routine to check the error of the Laplacian vs. its FD discretization

Args:
prob: a problem instance

Returns:
the error between the analytic Laplacian and the computed one of a given function
"""

# create x values, use only inner points
xvalues = np.array([(i + 1) * prob.dx for i in range(prob.nvars)])

# create a mesh instance and fill it with a sine wave
u = prob.dtype_u(init=prob.init)
u[:] = np.sin(np.pi * prob.freq * xvalues)

# create a mesh instance and fill it with the Laplacian of the sine wave
u_lap = prob.dtype_u(init=prob.init)
u_lap[:] = -((np.pi * prob.freq) ** 2) * prob.nu * np.sin(np.pi * prob.freq * xvalues)

# compare analytic and computed solution using the eval_f routine of the problem class
err = abs(prob.eval_f(u, 0) - u_lap)

return err

if __name__ == "__main__":
main()
```

Results:

```Error of the spatial accuracy test: 1.981784e-04
```

## Part B: Spatial accuracy check¶

We now do a more thorough test of the accuracy of the spatial discretization. We loop over a set of `nvars`, compute the Laplacian of our test vector and look at the error. Then, the order of accuracy in space is checked by looking at the numbers and by looking at “points on a line”.

Important things to note:

• You better test your operators.. use nvars > 2**16 and things will break!

• Add your results into a dictionary for later usage. Use IDs to find the data! Also, use headers to store meta-information.

```from pathlib import Path
import matplotlib

matplotlib.use('Agg')

from collections import namedtuple
import matplotlib.pylab as plt
import numpy as np
import os.path

from pySDC.implementations.datatype_classes.mesh import mesh
from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_unforced

# setup id for gathering the results (will sort by nvars)
ID = namedtuple('ID', 'nvars')

def main():
"""
A simple test program to check order of accuracy in space for a simple test problem
"""

# initialize problem parameters
problem_params = dict()
problem_params['nu'] = 0.1  # diffusion coefficient
problem_params['freq'] = 4  # frequency for the test value
problem_params['bc'] = 'dirichlet-zero'  # boundary conditions

# create list of nvars to do the accuracy test with
nvars_list = [2**p - 1 for p in range(4, 15)]

# run accuracy test for all nvars
results = run_accuracy_check(nvars_list=nvars_list, problem_params=problem_params)

# compute order of accuracy
order = get_accuracy_order(results)

Path("data").mkdir(parents=True, exist_ok=True)
f = open('data/step_1_B_out.txt', 'w')
for l in range(len(order)):
out = 'Expected order: %2i -- Computed order %4.3f' % (2, order[l])
f.write(out + '\n')
print(out)
f.close()

# visualize results
plot_accuracy(results)

assert os.path.isfile('data/step_1_accuracy_test_space.png'), 'ERROR: plotting did not create file'

assert all(np.isclose(order, 2, rtol=0.06)), "ERROR: spatial order of accuracy is not as expected, got %s" % order

def run_accuracy_check(nvars_list, problem_params):
"""
Routine to check the error of the Laplacian vs. its FD discretization

Args:
nvars_list: list of nvars to do the testing with
problem_params: dictionary containing the problem-dependent parameters

Returns:
a dictionary containing the errors and a header (with nvars_list)
"""

results = {}
# loop over all nvars
for nvars in nvars_list:
# setup problem
problem_params['nvars'] = nvars
prob = heatNd_unforced(**problem_params)

# create x values, use only inner points
xvalues = np.array([(i + 1) * prob.dx for i in range(prob.nvars)])

# create a mesh instance and fill it with a sine wave
u = prob.u_exact(t=0)

# create a mesh instance and fill it with the Laplacian of the sine wave
u_lap = prob.dtype_u(init=prob.init)
u_lap[:] = -((np.pi * prob.freq) ** 2) * prob.nu * np.sin(np.pi * prob.freq * xvalues)

# compare analytic and computed solution using the eval_f routine of the problem class
err = abs(prob.eval_f(u, 0) - u_lap)

# get id for this nvars and put error into dictionary
id = ID(nvars=nvars)
results[id] = err

# add nvars_list to dictionary for easier access later on
results['nvars_list'] = nvars_list

return results

def get_accuracy_order(results):
"""
Routine to compute the order of accuracy in space

Args:
results: the dictionary containing the errors

Returns:
the list of orders
"""

# retrieve the list of nvars from results
assert 'nvars_list' in results, 'ERROR: expecting the list of nvars in the results dictionary'
nvars_list = sorted(results['nvars_list'])

order = []
# loop over two consecutive errors/nvars pairs
for i in range(1, len(nvars_list)):
# get ids
id = ID(nvars=nvars_list[i])
id_prev = ID(nvars=nvars_list[i - 1])

# compute order as log(prev_error/this_error)/log(this_nvars/old_nvars) <-- depends on the sorting of the list!
tmp = np.log(results[id_prev] / results[id]) / np.log(nvars_list[i] / nvars_list[i - 1])
order.append(tmp)

return order

def plot_accuracy(results):
"""
Routine to visualize the errors as well as the expected errors

Args:
results: the dictionary containing the errors
"""

# retrieve the list of nvars from results
assert 'nvars_list' in results, 'ERROR: expecting the list of nvars in the results dictionary'
nvars_list = sorted(results['nvars_list'])

# Set up plotting parameters
params = {
'legend.fontsize': 20,
'figure.figsize': (12, 8),
'axes.labelsize': 20,
'axes.titlesize': 20,
'xtick.labelsize': 16,
'ytick.labelsize': 16,
'lines.linewidth': 3,
}
plt.rcParams.update(params)

# create new figure
plt.figure()
# take x-axis limits from nvars_list + some spacing left and right
plt.xlim([min(nvars_list) / 2, max(nvars_list) * 2])
plt.xlabel('nvars')
plt.ylabel('abs. error')
plt.grid()

# get guide for the order of accuracy, i.e. the errors to expect
# get error for first entry in nvars_list
id = ID(nvars=nvars_list)
base_error = results[id]
# assemble optimal errors for 2nd order method and plot
order_guide_space = [base_error / (2 ** (2 * i)) for i in range(0, len(nvars_list))]
plt.loglog(nvars_list, order_guide_space, color='k', ls='--', label='2nd order')

min_err = 1e99
max_err = 0e00
err_list = []
# loop over nvars, get errors and find min/max error for y-axis limits
for nvars in nvars_list:
id = ID(nvars=nvars)
err = results[id]
min_err = min(err, min_err)
max_err = max(err, max_err)
err_list.append(err)
plt.loglog(nvars_list, err_list, ls=' ', marker='o', markersize=10, label='experiment')

plt.ylim([min_err / 10, max_err * 10])
plt.legend(loc=1, ncol=1, numpoints=1)

# save plot as PDF, beautify
fname = 'data/step_1_accuracy_test_space.png'
plt.savefig(fname, bbox_inches='tight')

return None

if __name__ == "__main__":
main()
```

Results:

```Expected order:  2 -- Computed order 1.888
Expected order:  2 -- Computed order 1.949
Expected order:  2 -- Computed order 1.976
Expected order:  2 -- Computed order 1.988
Expected order:  2 -- Computed order 1.994
Expected order:  2 -- Computed order 1.997
Expected order:  2 -- Computed order 1.999
Expected order:  2 -- Computed order 1.999
Expected order:  2 -- Computed order 1.999
Expected order:  2 -- Computed order 1.976
``` ## Part C: Collocation problem setup¶

Here, we set up our first collocation problem using one of the `collocation_classes`, namely `GaussRadau_Right`. Using the spatial operator, we can build the collocation problem which in this case is linear. This fully coupled system is then solved directly and the error is compared to the exact solution.

Important things to note:

• The collocation matrix `Q` is and will be always relative to the temporal domain [0,1]. Use `dt` to scale it appropriately.

• Although convenient to analyze, the matrix formulation is not suited for larger (in space or time) computations. All remaining computations in pySDC are thus based on decoupling space and time operators (i.e. no Kronecker product).

• We can use the `u_exact` routine here to return the values at any given point in time. It is recommended to include either an initialization routine or the exact solution (if applicable) into the problem class.

• This is where the fun with parameters comes in: How many DOFs do we need in space? How large/small does dt have to be? What frequency in space is fine? …

```import numpy as np
import scipy.sparse as sp
from pathlib import Path

from pySDC.core.Collocation import CollBase
from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_unforced

def main():
"""
A simple test program to create and solve a collocation problem directly
"""

# instantiate problem
prob = heatNd_unforced(
nvars=1023,  # number of degrees of freedom
nu=0.1,  # diffusion coefficient
freq=4,  # frequency for the test value
bc='dirichlet-zero',  # boundary conditions
)

# instantiate collocation class, relative to the time interval [0,1]

# set time-step size (warning: the collocation matrices are relative to [0,1], see above)
dt = 0.1

# solve collocation problem
err = solve_collocation_problem(prob=prob, coll=coll, dt=dt)

Path("data").mkdir(parents=True, exist_ok=True)
f = open('data/step_1_C_out.txt', 'w')
out = 'Error of the collocation problem: %8.6e' % err
f.write(out + '\n')
print(out)
f.close()

assert err <= 4e-04, "ERROR: did not get collocation error as expected, got %s" % err

def solve_collocation_problem(prob, coll, dt):
"""
Routine to build and solve the linear collocation problem

Args:
prob: a problem instance
coll: a collocation instance
dt: time-step size

Return:
the analytic error of the solved collocation problem
"""

# shrink collocation matrix: first line and column deals with initial value, not needed here
Q = coll.Qmat[1:, 1:]

# build system matrix M of collocation problem
M = sp.eye(prob.nvars * coll.num_nodes) - dt * sp.kron(Q, prob.A)

# get initial value at t0 = 0
u0 = prob.u_exact(t=0)
# fill in u0-vector as right-hand side for the collocation problem
u0_coll = np.kron(np.ones(coll.num_nodes), u0)
# get exact solution at Tend = dt
uend = prob.u_exact(t=dt)

# solve collocation problem directly
u_coll = sp.linalg.spsolve(M, u0_coll)

# compute error
err = np.linalg.norm(u_coll[-prob.nvars :] - uend, np.inf)

return err

if __name__ == "__main__":
main()
```

Results:

```Error of the collocation problem: 3.803471e-04
```

## Part D: Collocation accuracy test¶

As for the spatial operator, we now test the accuracy in time. This time, we loop over a set of `dt`, compute the solution to the collocation problem and look at the error.

Important things to note:

• We take a large number of DOFs in space, since we need to beat 5th order in time with a 2nd order stencil in space.

• Orders of convergence are not as stable as for the space-only test. One of the problems of this example is that we are actually trying to compute 0 very, very thoroughly…

```from pathlib import Path
import matplotlib

matplotlib.use('Agg')

from collections import namedtuple
import matplotlib.pylab as plt
import numpy as np
import os.path
import scipy.sparse as sp

from pySDC.core.Collocation import CollBase
from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_unforced

# setup id for gathering the results (will sort by dt)
ID = namedtuple('ID', 'dt')

def main():
"""
A simple test program to compute the order of accuracy in time
"""

# initialize problem parameters
problem_params = {}
problem_params['nu'] = 0.1  # diffusion coefficient
problem_params['freq'] = 4  # frequency for the test value
problem_params['nvars'] = 16383  # number of DOFs in space
problem_params['bc'] = 'dirichlet-zero'  # boundary conditions

# instantiate problem
prob = heatNd_unforced(**problem_params)

# instantiate collocation class, relative to the time interval [0,1]

# assemble list of dt
dt_list = [0.1 / 2**p for p in range(0, 4)]

# run accuracy test for all dt
results = run_accuracy_check(prob=prob, coll=coll, dt_list=dt_list)

# get order of accuracy
order = get_accuracy_order(results)

Path("data").mkdir(parents=True, exist_ok=True)
f = open('data/step_1_D_out.txt', 'w')
for l in range(len(order)):
out = 'Expected order: %2i -- Computed order %4.3f' % (5, order[l])
f.write(out + '\n')
print(out)
f.close()

# visualize results
plot_accuracy(results)

assert os.path.isfile('data/step_1_accuracy_test_coll.png')

assert all(np.isclose(order, 2 * coll.num_nodes - 1, rtol=0.4)), (
"ERROR: did not get order of accuracy as expected, got %s" % order
)

def run_accuracy_check(prob, coll, dt_list):
"""
Routine to build and solve the linear collocation problem

Args:
prob: a problem instance
coll: a collocation instance
dt_list: list of time-step sizes

Return:
the analytic error of the solved collocation problem
"""

results = {}
# loop over all nvars
for dt in dt_list:
# shrink collocation matrix: first line and column deals with initial value, not needed here
Q = coll.Qmat[1:, 1:]

# build system matrix M of collocation problem
M = sp.eye(prob.nvars * coll.num_nodes) - dt * sp.kron(Q, prob.A)

# get initial value at t0 = 0
u0 = prob.u_exact(t=0)
# fill in u0-vector as right-hand side for the collocation problem
u0_coll = np.kron(np.ones(coll.num_nodes), u0)
# get exact solution at Tend = dt
uend = prob.u_exact(t=dt)

# solve collocation problem directly
u_coll = sp.linalg.spsolve(M, u0_coll)

# compute error
err = np.linalg.norm(u_coll[-prob.nvars :] - uend, np.inf)
# get id for this dt and store error in results
id = ID(dt=dt)
results[id] = err

# add list of dt to results for easier access
results['dt_list'] = dt_list
return results

def get_accuracy_order(results):
"""
Routine to compute the order of accuracy in time

Args:
results: the dictionary containing the errors

Returns:
the list of orders
"""

# retrieve the list of dt from results
assert 'dt_list' in results, 'ERROR: expecting the list of dt in the results dictionary'
dt_list = sorted(results['dt_list'], reverse=True)

order = []
# loop over two consecutive errors/dt pairs
for i in range(1, len(dt_list)):
# get ids
id = ID(dt=dt_list[i])
id_prev = ID(dt=dt_list[i - 1])

# compute order as log(prev_error/this_error)/log(this_dt/old_dt) <-- depends on the sorting of the list!
tmp = np.log(results[id] / results[id_prev]) / np.log(dt_list[i] / dt_list[i - 1])
order.append(tmp)

return order

def plot_accuracy(results):
"""
Routine to visualize the errors as well as the expected errors

Args:
results: the dictionary containing the errors
"""

# retrieve the list of nvars from results
assert 'dt_list' in results, 'ERROR: expecting the list of dts in the results dictionary'
dt_list = sorted(results['dt_list'])

# Set up plotting parameters
params = {
'legend.fontsize': 20,
'figure.figsize': (12, 8),
'axes.labelsize': 20,
'axes.titlesize': 20,
'xtick.labelsize': 16,
'ytick.labelsize': 16,
'lines.linewidth': 3,
}
plt.rcParams.update(params)

# create new figure
plt.figure()
# take x-axis limits from nvars_list + some spacning left and right
plt.xlim([min(dt_list) / 2, max(dt_list) * 2])
plt.xlabel('dt')
plt.ylabel('abs. error')
plt.grid()

# get guide for the order of accuracy, i.e. the errors to expect
# get error for first entry in nvars_list
id = ID(dt=dt_list)
base_error = results[id]
# assemble optimal errors for 5th order method and plot
order_guide_space = [base_error * (2 ** (5 * i)) for i in range(0, len(dt_list))]
plt.loglog(dt_list, order_guide_space, color='k', ls='--', label='5th order')

min_err = 1e99
max_err = 0e00
err_list = []
# loop over nvars, get errors and find min/max error for y-axis limits
for dt in dt_list:
id = ID(dt=dt)
err = results[id]
min_err = min(err, min_err)
max_err = max(err, max_err)
err_list.append(err)
plt.loglog(dt_list, err_list, ls=' ', marker='o', markersize=10, label='experiment')

plt.ylim([min_err / 10, max_err * 10])
plt.legend(loc=2, ncol=1, numpoints=1)

# save plot as PDF, beautify
fname = 'data/step_1_accuracy_test_coll.png'
plt.savefig(fname, bbox_inches='tight')

return None

if __name__ == "__main__":
main()
```

Results:

```Expected order:  5 -- Computed order 4.791
Expected order:  5 -- Computed order 5.364
Expected order:  5 -- Computed order 5.662
``` 