Advanced usage¶
This page contains short examples that demonstrate advanced usage of PyMGRIT. The source code for these and more examples is available in the examples folder.
Advanced multigrid hierarchy¶
example_time_integrators.py and example_heat_1d_bdf2.py
PyMGRIT allows using different application classes and/or time integration schemes in the multigrid hierarchy.
Example 1 shows how to implement different time integration methods in an application class.
Example 2 shows how to use multiple application classes in the multigrid hierarchy.
Example 1 One application class with different time integration methods
The member function step() of an application class can carry out a time integration step based on different time integration methods. The Dahlquist application class implements the following time integration schemes:
Backward Euler
Forward Euler
Trapezoidal rule
Implicit mid-point rule
that can be controlled by the member variable method:
def step(self, u_start: VectorDahlquist, t_start: float, t_stop: float) -> VectorDahlquist:
"""
Time integration routine for Dahlquist's test problem:
BE: Backward Euler
FE: Forward Euler
TR: Trapezoidal rule
MR: implicit Mid-point rule
:param u_start: approximate solution for the input time t_start
:param t_start: time associated with the input approximate solution u_start
:param t_stop: time to evolve the input approximate solution to
:return: approximate solution for the input time t_stop
"""
z = (t_stop - t_start) * -1 # Note: lambda = -1
if self.method == 'BE':
tmp = 1 / (1 - z) * u_start.get_values()
elif self.method == 'FE':
tmp = (1 + z) * u_start.get_values()
elif self.method == 'TR':
tmp = (1 + z / 2) / (1 - z / 2) * u_start.get_values()
elif self.method == 'MR':
k1 = -1 / (1 - z / 2) * u_start.get_values()
tmp = u_start.get_values() + (t_stop - t_start) * k1
return VectorDahlquist(tmp)
The corresponding example (example_time_integrators.py) creates a two-level hierarchy for Dahlquist’s test problem, using the implicit mid-point rule on the fine grid (level 0) and backward Euler on the coarse grid (level 1):
# Create Dahlquist's test problem using implicit mid-point rule time integration
dahlquist_lvl0 = Dahlquist(t_start=0, t_stop=5, nt=101, method='MR')
# Create Dahlquist's test problem using backward Euler time integration
dahlquist_lvl1 = Dahlquist(t_start=0, t_stop=5, nt=51, method='BE')
# Setup an MGRIT solver and solve the problem
mgrit = Mgrit(problem=[dahlquist_lvl0, dahlquist_lvl1])
info = mgrit.solve()
Example 2 Two application classes
In the second example, we use two application classes for implementing two different time integration methods for the 1D heat equation example:
Note: The vector class used in both application classes contains the solution at two consecutive time points.
The corresponding example (example_heat_1d_bdf2.py) constructs a three-level multigrid hierarchy for the 1D heat equation example using the BDF2 application class on the fine grid (level 0) and the BDF1 application class on the first and second coarse grids (levels 1 and 2):
import numpy as np
from pymgrit.core.mgrit import Mgrit
from pymgrit.heat.heat_1d_2pts_bdf1 import Heat1DBDF1
from pymgrit.heat.heat_1d_2pts_bdf2 import Heat1DBDF2
def rhs(x, t):
"""
Right-hand side of 1D heat equation example problem at a given space-time point (x,t),
-sin(pi*x)(sin(t) - a*pi^2*cos(t)), a = 1
:param x: spatial grid point
:param t: time point
:return: right-hand side of 1D heat equation example problem at point (x,t)
"""
return - np.sin(np.pi * x) * (np.sin(t) - 1 * np.pi ** 2 * np.cos(t))
def init_cond(x):
"""
Initial condition of 1D heat equation example,
u(x,0) = sin(pi*x)
:param x: spatial grid point
:return: initial condition of 1D heat equation example problem
"""
return np.sin(np.pi * x)
# Time interval
t_start = 0
t_stop = 2
nt = 512 # number of time points excluding t_start
dtau = t_stop / nt # time-step size
# Time points are grouped in pairs of two consecutive time points
# => (nt/2) + 1 pairs
# Note: * Each pair is associated with the time value of its first point.
# * The second value of the last pair (associated with t_stop) is not used.
# * The spacing within each pair is the same (= dt) on all grid levels.
t_interval = np.linspace(t_start, t_stop, int(nt / 2 + 1))
heat0 = Heat1DBDF2(x_start=0, x_end=1, nx=1001, a=1, dtau=dtau, rhs=rhs, init_cond=init_cond,
t_interval=t_interval)
heat1 = Heat1DBDF1(x_start=0, x_end=1, nx=1001, a=1, dtau=dtau, rhs=rhs, init_cond=init_cond,
t_interval=heat0.t[::2])
heat2 = Heat1DBDF1(x_start=0, x_end=1, nx=1001, a=1, dtau=dtau, rhs=rhs, init_cond=init_cond,
t_interval=heat1.t[::2])
# Setup three-level MGRIT solver and solve the problem
problem = [heat0, heat1, heat2]
mgrit = Mgrit(problem=problem)
info = mgrit.solve()
Spatial coarsening¶
This example demonstrates how to use the transfer parameter transfer of the MGRIT solver to apply spatial coarsening on different levels of the time-grid hierarchy for solving a 1D heat equation problem (see Heat Equation).
The first step is to import all necessary PyMGRIT classes (and numpy for later use):
import numpy as np
from pymgrit.heat.heat_1d import Heat1D # 1D Heat equation problem
from pymgrit.heat.heat_1d import VectorHeat1D # 1D Heat equation vector class
from pymgrit.core.mgrit import Mgrit # MGRIT solver
from pymgrit.core.grid_transfer import GridTransfer # Parent grid transfer class
from pymgrit.core.grid_transfer_copy import GridTransferCopy # Copy transfer class
Then, we define the class GridTransferHeat for the 1D heat equation:
class GridTransferHeat(GridTransfer):
"""
Grid Transfer for the Heat Equation.
Interpolation: Linear interpolation
Restriction: Full weighting
"""
def __init__(self):
"""
Constructor.
:rtype: GridTransferHeat object
"""
super().__init__()
The grid transfer class must contain the two member functions restriction() and interpolation().
The function restriction() receives a VectorHeat1D object and returns another VectorHeat1D object that contains the restricted solution vector:
def restriction(self, u: VectorHeat1D) -> VectorHeat1D:
"""
Restrict input vector u using standard full weighting restriction.
Note: In the 1d heat equation example, we consider homogeneous Dirichlet BCs in space.
The Heat1D vector class only stores interior points.
:param u: approximate solution vector
:return: input solution vector u restricted to a coarse grid
"""
# Get values at interior points
sol = u.get_values()
# Create array for restricted values
ret_array = np.zeros(int((len(sol) - 1) / 2))
# Full weighting restriction
for i in range(len(ret_array)):
ret_array[i] = sol[2 * i] * 1 / 4 + sol[2 * i + 1] * 1 / 2 + sol[2 * i + 2] * 1 / 4
# Create and return a VectorHeat1D object with the restricted values
ret = VectorHeat1D(len(ret_array))
ret.set_values(ret_array)
return ret
Similarly, we define the function interpolation() as follows:
def interpolation(self, u: VectorHeat1D) -> VectorHeat1D:
"""
Interpolate input vector u using linear interpolation.
Note: In the 1d heat equation example, we consider homogeneous Dirichlet BCs in space.
The Heat1D vector class only stores interior points.
:param u: approximate solution vector
:return: input solution vector u interpolated to a fine grid
"""
# Get values at interior points
sol = u.get_values()
# Create array for interpolated values
ret_array = np.zeros(int(len(sol) * 2 + 1))
# Linear interpolation
for i in range(len(sol)):
ret_array[i * 2] += 1 / 2 * sol[i]
ret_array[i * 2 + 1] += sol[i]
ret_array[i * 2 + 2] += 1 / 2 * sol[i]
# Create and return a VectorHeat1D object with interpolated values
ret = VectorHeat1D(len(ret_array))
ret.set_values(ret_array)
return ret
Now, we construct a multigrid hierarchy for the 1d heat example. Here, we set up the following hierarchy:
level 0: 129 time points, 17 points in space
level 1: 65 time points, 9 points in space
level 2: 33 time points, 5 points in space
level 3: 17 time points, 5 points in space
Note: In this example, it is not possible to use PyMGRIT’s core function simple_setup_problem(), since the number of spatial grid points changes in the multigrid hiearchy:
def rhs(x, t):
return - np.sin(np.pi * x) * (np.sin(t) - 1 * np.pi ** 2 * np.cos(t))
def init_cond(x):
return np.sin(np.pi * x)
heat0 = Heat1D(x_start=0, x_end=2, nx=2 ** 4 + 1, a=1, rhs=rhs, init_cond=init_cond, t_start=0, t_stop=2,
nt=2 ** 7 + 1)
heat1 = Heat1D(x_start=0, x_end=2, nx=2 ** 3 + 1, a=1, rhs=rhs, init_cond=init_cond, t_interval=heat0.t[::2])
heat2 = Heat1D(x_start=0, x_end=2, nx=2 ** 2 + 1, a=1, rhs=rhs, init_cond=init_cond, t_interval=heat1.t[::2])
heat3 = Heat1D(x_start=0, x_end=2, nx=2 ** 2 + 1, a=1, rhs=rhs, init_cond=init_cond, t_interval=heat2.t[::2])
problem = [heat0, heat1, heat2, heat3]
Before we can set up the MGRIT solver, we have to define the grid transfer between all two consecutive levels in the multigrid hierarchy. These grid transfers are specified by a list of grid transfer objects of length (#levels -1). For our four-level example, this list is of length three with two objects of the new class GridTransferHeat for the transfer between levels 0 and 1 as well as between levels 1 and 2 and an object of PyMGRIT’s core class GridTransferCopy for the transfer between levels 2 and 3:
transfer = [GridTransferHeat(), GridTransferHeat(), GridTransferCopy()]
Finally, we set up the MGRIT solver and solve the problem:
mgrit = Mgrit(problem=problem, transfer=transfer)
info = mgrit.solve()
Complete code:
import numpy as np
from pymgrit.heat.heat_1d import Heat1D # 1D Heat equation problem
from pymgrit.heat.heat_1d import VectorHeat1D # 1D Heat equation vector class
from pymgrit.core.mgrit import Mgrit # MGRIT solver
from pymgrit.core.grid_transfer import GridTransfer # Parent grid transfer class
from pymgrit.core.grid_transfer_copy import GridTransferCopy # Copy transfer class
# Create class for the grid transfer between spatial grids.
# Note: The class must inherit from PyMGRIT's core GridTransfer class.
class GridTransferHeat(GridTransfer):
"""
Grid Transfer class for the Heat Equation.
Interpolation: Linear interpolation
Restriction: Full weighting
"""
def __init__(self):
"""
Constructor.
:rtype: GridTransferHeat object
"""
super().__init__()
# Define restriction operator
def restriction(self, u: VectorHeat1D) -> VectorHeat1D:
"""
Restrict input vector u using standard full weighting restriction.
Note: In the 1d heat equation example, we consider homogeneous Dirichlet BCs in space.
The Heat1D vector class only stores interior points.
:param u: approximate solution vector
:return: input solution vector u restricted to a coarse grid
"""
# Get values at interior points
sol = u.get_values()
# Create array for restricted values
ret_array = np.zeros(int((len(sol) - 1) / 2))
# Full weighting restriction
for i in range(len(ret_array)):
ret_array[i] = sol[2 * i] * 1 / 4 + sol[2 * i + 1] * 1 / 2 + sol[2 * i + 2] * 1 / 4
# Create and return a VectorHeat1D object with the restricted values
ret = VectorHeat1D(len(ret_array))
ret.set_values(ret_array)
return ret
# Define interpolation operator
def interpolation(self, u: VectorHeat1D) -> VectorHeat1D:
"""
Interpolate input vector u using linear interpolation.
Note: In the 1d heat equation example, we consider homogeneous Dirichlet BCs in space.
The Heat1D vector class only stores interior points.
:param u: approximate solution vector
:return: input solution vector u interpolated to a fine grid
"""
# Get values at interior points
sol = u.get_values()
# Create array for interpolated values
ret_array = np.zeros(int(len(sol) * 2 + 1))
# Linear interpolation
for i in range(len(sol)):
ret_array[i * 2] += 1 / 2 * sol[i]
ret_array[i * 2 + 1] += sol[i]
ret_array[i * 2 + 2] += 1 / 2 * sol[i]
# Create and return a VectorHeat1D object with interpolated values
ret = VectorHeat1D(len(ret_array))
ret.set_values(ret_array)
return ret
# Construct a four-level multigrid hierarchy for the 1d heat example
# * use a coarsening factor of 2 in time on all levels
# * apply spatial coarsening by a factor of 2 on the first two levels
heat0 = Heat1D(x_start=0, x_end=2, nx=2 ** 4 + 1, a=1, t_start=0, t_stop=2, nt=2 ** 7 + 1)
heat1 = Heat1D(x_start=0, x_end=2, nx=2 ** 3 + 1, a=1, t_interval=heat0.t[::2])
heat2 = Heat1D(x_start=0, x_end=2, nx=2 ** 2 + 1, a=1, t_interval=heat1.t[::2])
heat3 = Heat1D(x_start=0, x_end=2, nx=2 ** 2 + 1, a=1, t_interval=heat2.t[::2])
problem = [heat0, heat1, heat2, heat3]
# Specify a list of grid transfer operators of length (#levels - 1) for the transfer between two consecutive levels
# * Use the new class GridTransferHeat to apply spatial coarsening for transfers between the first three levels
# * Use PyMGRIT's core class GridTransferCopy for the transfer between the last two levels (no spatial coarsening)
transfer = [GridTransferHeat(), GridTransferHeat(), GridTransferCopy()]
# Setup four-level MGRIT solver and solve the problem
mgrit = Mgrit(problem=problem, transfer=transfer)
info = mgrit.solve()
Convergence criteria¶
example_convergence_criterion.py
In this example, we define a customized version of PyMGRIT’s MGRIT solver that uses a different convergence criterion. A two-level variant of this customized solver is then applied to compute an Arenstorf orbit.
All it takes to define a customized MGRIT solver is to create a new class that inherits from PyMGRIT’s core Mgrit class and to define in this child class the method convergence_criterion(), overwriting Mgrit’s member routine convergence_criterion() that is called in the algorithm after each iteration. In this example, we implement a convergence criterion based on the maximum norm of the relative difference at C-points of successive iterates:
class MgritCustomized(Mgrit):
"""
Customized MGRIT class.
Use maximum norm of the relative difference at C-points of two successive
iterates as convergence criterion.
"""
def __init__(self, *args, **kwargs) -> None:
"""
Cumstomized MGRIT constructor.
"""
# Call parent constructor
super().__init__(*args, **kwargs)
# New member variable for saving the C-point values of the last iteration
self.last_it = []
# Initialize the new member variable
self.convergence_criterion(iteration=0)
def convergence_criterion(self, iteration: int) -> None:
"""
Stopping criterion based on achieving a maximum relative difference at C-points
of two successive iterates below the specified stopping tolerance.
Note: The stopping tolerance is specified when setting up the solver.
:param iteration: Iteration number
"""
# Create list in the first function call
if len(self.last_it) != len(self.index_local_c[0]):
self.last_it = np.zeros((len(self.index_local_c[0]), len(self.u[0][0].get_values())))
new = np.zeros_like(self.last_it)
j = 0
tmp = 0
# If process has a C-point
if self.index_local_c[0].size > 0:
# Loop over all C-points of the process
for i in np.nditer(self.index_local_c[0]):
new[j] = self.u[0][i].get_values()
j = j + 1
# Compute relative difference between two iterates
tmp = 100 * np.max(
np.abs(np.abs(np.divide((new - self.last_it), new, out=np.zeros_like(self.last_it), where=new != 0))))
# Communicate the local value
tmp = self.comm_time.allgather(tmp)
# Take maximum norm
self.conv[iteration] = np.max(np.abs(tmp))
self.last_it = np.copy(new)
We can use the new class MgritCustomized to set up an MGRIT solver with the new convergence criterion and solve our problem in the usual way:
# Create two-level time-grid hierarchy for the ODE system describing Arenstorf orbits
ahrenstorf_lvl_0 = ArenstorfOrbit(t_start=0, t_stop=17.06521656015796, nt=10001)
ahrenstorf_lvl_1 = ArenstorfOrbit(t_interval=ahrenstorf_lvl_0.t[::100])
# Set up customized MGRIT solver and solve the problem.
# Note: Setting the solver tolerance to 1 means that iterations stop
# if the maximum relative change at C-points of all four variables of the ODE system
# is smaller than 1%.
info = MgritCustomized(problem=[ahrenstorf_lvl_0, ahrenstorf_lvl_1], tol=1).solve()