r"""Module for solving a 1D conservation law via DG.
Adapted from a Discontinuous Galerkin (DG) solver written
by Per Olof-Persson.
Check out an example `notebook`_ using these utilities to
solve the problem.
.. _notebook: http://nbviewer.jupyter.org/github/dhermes/\
berkeley-m273-s2016/blob/master/assignment1/dg1.ipynb
"""
import numpy as np
from numpy.polynomial import legendre
import six
_RK_STEPS = (4, 3, 2, 1)
# pylint: disable=too-few-public-methods
[docs]class MathProvider(object):
"""Mutable settings for doing math.
For callers that wish to swap out the default behavior -- for example,
to use high-precision values instead of floats -- this class can
just be monkey patched on the module.
The module assumes through-out that solution data is in NumPy arrays,
but the data inside those arrays may be any type.
.. note::
The callers assume :data:`exp_func` is a vectorized exponential
that can act on a NumPy array containing elements of the relevant
type.
.. note::
The :data:`zeros` constructor should also be able to take the
``order`` argument (and should produce a NumPy array).
"""
exp_func = staticmethod(np.exp)
linspace = staticmethod(np.linspace)
num_type = staticmethod(float)
mat_inv = staticmethod(np.linalg.inv)
solve = staticmethod(np.linalg.solve)
zeros = staticmethod(np.zeros)
# pylint: enable=too-few-public-methods
[docs]def gauss_lobatto_points(start, stop, num_points):
r"""Get the node points for Gauss-Lobatto quadrature.
Using :math:`n` points, this quadrature is accurate to degree
:math:`2n - 3`. The node points are :math:`x_1 = -1`,
:math:`x_n = 1` and the interior are :math:`n - 2` roots of
:math:`P'_{n - 1}(x)`.
Though we don't compute them here, the weights are
:math:`w_1 = w_n = \frac{2}{n(n - 1)}` and for the interior points
.. math::
w_j = \frac{2}{n(n - 1) \left[P_{n - 1}\left(x_j\right)\right]^2}
This is in contrast to the scheme used in Gaussian quadrature, which
use roots of :math:`P_n(x)` as nodes and use the weights
.. math::
w_j = \frac{2}{\left(1 - x_j\right)^2
\left[P'_n\left(x_j\right)\right]^2}
.. note::
This method is **not** generic enough to accommodate non-NumPy
types as it relies on the :mod:`numpy.polynomial.legendre`.
:type start: float
:param start: The beginning of the interval.
:type stop: float
:param stop: The end of the interval.
:type num_points: int
:param num_points: The number of points to use.
:rtype: :class:`numpy.ndarray`
:returns: 1D array, the interior quadrature nodes.
"""
p_n_minus1 = [0] * (num_points - 1) + [1]
inner_nodes = legendre.legroots(legendre.legder(p_n_minus1))
# Utilize symmetry about 0.
inner_nodes = 0.5 * (inner_nodes - inner_nodes[::-1])
if start != -1.0 or stop != 1.0:
# [-1, 1] --> [0, 2] --> [0, stop - start] --> [start, stop]
inner_nodes = start + (inner_nodes + 1.0) * 0.5 * (stop - start)
return np.hstack([[start], inner_nodes, [stop]])
[docs]def get_legendre_matrix(points, max_degree=None):
r"""Evaluate Legendre polynomials at a set of points.
If our points are :math:`x_0, \ldots, x_p`, this computes
.. math::
\left[ \begin{array}{c c c c}
L_0(x_0) & L_1(x_0) & \cdots & L_d(x_0) \\
L_0(x_1) & L_1(x_1) & \cdots & L_d(x_p) \\
\vdots & \vdots & \ddots & \vdots \\
L_0(x_p) & L_1(x_p) & \cdots & L_d(x_p)
\end{array}\right]
by utilizing the recurrence
.. math::
n L_n(x) = (2n - 1) x L_{n - 1}(x) - (n - 1) L_{n - 2}(x)
:type points: :class:`numpy.ndarray`
:param points: 1D array. The points at which to evaluate Legendre
polynomials.
:type max_degree: int
:param max_degree: (Optional) The maximum degree of Legendre
polynomial to use. Defaults to one less than the
number of points (which will produce a square
output).
:rtype: :class:`numpy.ndarray`
:returns: The 2D array containing the Legendre polynomials evaluated
at our input points.
"""
num_points, = np.shape(points)
if max_degree is None:
max_degree = num_points - 1
# Use Fortran order since we operate on columns.
result = MathProvider.zeros((num_points, max_degree + 1), order='F')
result[:, 0] = MathProvider.num_type(1.0)
result[:, 1] = points
for degree in six.moves.xrange(2, max_degree + 1):
result[:, degree] = (
(2 * degree - 1) * points * result[:, degree - 1] -
(degree - 1) * result[:, degree - 2]) / degree
return result
def _find_matrices_helper(vals1, vals2):
r"""Helper for :func:`find_matrices`.
Computes a shoelace-like product of two vectors :math:`u, v`
via
.. math::
u_0 (v_1 + v_3 + \cdots) + u_1 (v_2 + v_4 + \cdots) +
u_{p - 1} v_p
:type vals1: :class:`numpy.ndarray`.
:param vals1: The vector :math`u`.
:type vals2: :class:`numpy.ndarray`.
:param vals2: The vector :math`v`.
:rtype: float
:returns: The shoelace-like product of the vectors.
"""
result = 0
for i, val in enumerate(vals1):
result += val * np.sum(vals2[i + 1::2])
return result
[docs]def get_evenly_spaced_points(start, stop, num_points):
"""Get points on an interval that are evenly spaced.
This is intended to be used to give points on a reference
interval when using DG on the 1D problem.
:type start: float
:param start: The beginning of the interval.
:type stop: float
:param stop: The end of the interval.
:type num_points: int
:param num_points: The number of points to use on the interval.
:rtype: :class:`numpy.ndarray`
:returns: The evenly spaced points on the interval.
"""
return MathProvider.linspace(start, stop, num_points)
[docs]def find_matrices(p_order, points_on_ref_int=None):
r"""Find mass and stiffness matrices.
We do this on the reference interval :math:`\left[-1, 1\right]`.
By default we use the evenly spaced points
.. math::
x_0 = -1, x_1 = -(p - 2)/p, \ldots, x_p = 1
but the set of nodes to use on the reference interval can be specified
via the ``points_on_ref_int`` argument. With our points, we
compute the polynomials :math:`\varphi_j(x)` such that
:math:`\varphi_j\left(x_i\right) = \delta_{ij}`. We do this by
writing
.. math::
\varphi_j(x) = \sum_{n = 0}^p c_n^{(j)} L_n(x)
where :math:`L_n(x)` is the Legendre polynomial of degree :math:`n`.
With this representation, we need to solve
.. math::
\left[ \begin{array}{c c c c}
L_0(x_0) & L_1(x_0) & \cdots & L_p(x_0) \\
L_0(x_1) & L_1(x_1) & \cdots & L_p(x_p) \\
\vdots & \vdots & \ddots & \vdots \\
L_0(x_p) & L_1(x_p) & \cdots & L_p(x_p)
\end{array}\right]
\left[ \begin{array}{c c c c}
c_0^{(0)} & c_0^{(1)} & \cdots & c_0^{(p)} \\
c_1^{(0)} & c_1^{(1)} & \cdots & c_1^{(p)} \\
\vdots & \vdots & \ddots & \vdots \\
c_p^{(0)} & c_p^{(1)} & \cdots & c_p^{(p)}
\end{array}\right]
= \left(\delta_{ij}\right) = I_{p + 1}
Then use these to compute the mass matrix
.. math::
M_{ij} = \int_{-1}^1 \varphi_i(x) \varphi_j(x) \, dx
and the stiffness matrix
.. math::
K_{ij} = \int_{-1}^1 \varphi_i'(x) \varphi_j(x) \, dx
Utilizing the fact that
.. math::
\left\langle L_n, L_m \right\rangle =
\int_{-1}^1 L_n(x) L_m(x) \, dx =
\frac{2}{2n + 1} \delta_{nm}
we can compute
.. math::
M_{ij} = \left\langle \varphi_i, \varphi_j \right\rangle =
\sum_{n, m} \left\langle c_n^{(i)} L_n,
c_m^{(j)} L_m \right\rangle =
\sum_{n = 0}^p \frac{2}{2n + 1} c_n^{(i)} c_n^{(j)}.
Similarly
.. math::
\left\langle L_n'(x), L_m(x) \right\rangle =
\begin{cases}
2 & \text{ if } n > m \text{ and }
n - m \equiv 1 \bmod{2} \\
0 & \text{ otherwise}.
\end{cases}
gives
.. math::
\begin{align*}
K_{ij} &= \left\langle \varphi_i', \varphi_j \right\rangle
= \sum_{n, m} \left\langle c_n^{(i)} L_n',
c_m^{(j)} L_m \right\rangle \\
&= 2 \left(c_0^{(j)} \left(c_1^{(i)} + c_3^{(i)} +
\cdots\right)
+ c_1^{(j)} \left(c_2^{(i)} + c_4^{(i)} +
\cdots\right)
+ \cdots
+ c_{p - 1}^{(j)} c_p^{(i)}\right) \\
\end{align*}
(For more general integrals, one might use Gaussian quadrature.
The largest degree integrand :math:`\varphi_i \varphi_j` has
degree :math:`2 p` so this would require :math:`n = p + 1` points
to be exact up to degree :math:`2(p + 1) - 1 = 2p + 1`.)
:type p_order: int
:param p_order: The degree of precision for the method.
:type points_on_ref_int: :data:`function <types.FunctionType>`
:param points_on_ref_int: (Optional) The method used to partition the
reference interval :math:`\left[0, h\right]`
into node points. Defaults to
:func:`get_evenly_spaced_points`.
:rtype: tuple
:returns: Pair of mass and stiffness matrices, square
:class:`numpy.ndarray` of dimension ``p_order + 1``.
"""
if points_on_ref_int is None:
points_on_ref_int = get_evenly_spaced_points
# Find the coefficients of the L_n(x) for each basis function.
x_vals = points_on_ref_int(-1, 1, p_order + 1)
coeff_mat = MathProvider.mat_inv(get_legendre_matrix(x_vals))
# Populate the mass and stiffness matrices.
legendre_norms = (MathProvider.num_type(2.0) /
np.arange(1, 2 * p_order + 2, 2))
mass_mat = MathProvider.zeros((p_order + 1, p_order + 1))
stiffness_mat = MathProvider.zeros((p_order + 1, p_order + 1))
for i in six.moves.xrange(p_order + 1):
mass_mat_base = legendre_norms * coeff_mat[:, i]
for j in six.moves.xrange(i, p_order + 1):
mass_mat[i, j] = np.sum(mass_mat_base * coeff_mat[:, j])
stiffness_mat[i, j] = 2.0 * _find_matrices_helper(
coeff_mat[:, j], coeff_mat[:, i])
if j > i:
mass_mat[j, i] = mass_mat[i, j]
stiffness_mat[j, i] = -stiffness_mat[i, j]
return mass_mat, stiffness_mat
[docs]def low_storage_rk(ode_func, u_val, dt):
r"""Update an ODE solutuon with an order 2/4 Runge-Kutta function.
The method is given by the following Butcher array:
.. math::
\begin{array}{c | c c c c}
0 & 0 & & & \\
1/4 & 1/4 & 0 & & \\
1/3 & 0 & 1/3 & 0 & \\
1/2 & 0 & 0 & 1/2 & 0 \\
\hline
& 0 & 0 & 0 & 1
\end{array}
It is advantageous because the updates :math:`k_j` can be over-written at
each step, since they are never re-used.
One can see that this method is **order 2** for general
:math:`\dot{u} = f(u)` by verifying that not all order 3 node
conditions are satisfied. For example:
.. math::
\frac{1}{3} \neq \sum_i b_i c_i^2 = 0 + 0 + 0 +
1 \cdot \left(\frac{1}{2}\right)^2
However, for linear ODEs, the method is **order 4**. To see this, note
that the test problem :math:`\dot{u} = \lambda u` gives the stability
function
.. math::
R\left(\lambda \Delta t\right) = R(z) =
1 + z + \frac{z^2}{2} + \frac{z^3}{6} + \frac{z^4}{24}
which matches the Taylor series for :math:`e^z` to order 4.
See `Problem Set 3`_ from Persson's Math 228A for more details.
.. _Problem Set 3: http://persson.berkeley.edu/228A/ps3.pdf
:type ode_func: callable
:param ode_func: The RHS in the ODE :math:`\dot{u} = f(u)`.
:type u_val: :class:`numpy.ndarray`
:param u_val: The input to :math:`f(u)`.
:type dt: float
:param dt: The timestep to use.
:rtype: :class:`numpy.ndarray`
:returns: The updated solution value.
"""
u_prev = u_val
# NOTE: We may not actually need to copy this if the ``ode_func`` does
# not mutate ``u_val`` in place, but we copy anyhow to avoid any
# accidental mutation.
u_curr = u_val.copy()
for irk in _RK_STEPS:
u_curr = u_prev + dt / irk * ode_func(u_curr)
return u_curr
[docs]def get_node_points(num_points, p_order, step_size=None,
points_on_ref_int=None):
r"""Return node points to split unit interval for DG.
:type num_points: int
:param num_points: The number :math:`n` of intervals to divide
:math:`\left[0, 1\right]` into.
:type p_order: int
:param p_order: The degree of precision for the method.
:type step_size: float
:param step_size: (Optional) The step size :math:`1 / n`.
:type points_on_ref_int: :data:`function <types.FunctionType>`
:param points_on_ref_int: (Optional) The method used to partition the
reference interval :math:`\left[0, h\right]`
into node points. Defaults to
:func:`get_evenly_spaced_points`.
:rtype: :class:`numpy.ndarray`
:returns: The :math:`x`-values for the node points, with
``p_order + 1`` rows and :math:`n` columns. The columns
correspond to each sub-interval and the rows correspond
to the node points within each sub-interval.
"""
if step_size is None:
step_size = MathProvider.num_type(1.0) / num_points
if points_on_ref_int is None:
points_on_ref_int = get_evenly_spaced_points
interval_starts = MathProvider.linspace(0, 1 - step_size, num_points)
# Split the first interval [0, h] in ``p_order + 1`` points.
first_interval = points_on_ref_int(0, step_size, p_order + 1)
# Broadcast the values with ``first_interval`` as rows and
# columns as ``interval_starts``.
return (first_interval[:, np.newaxis] +
interval_starts[np.newaxis, :])
[docs]def get_gaussian_like_initial_data(node_points):
r"""Get the default initial solution data.
In this case it is
.. math::
u(x, 0) = \exp\left(-\left(\frac{x - \frac{1}{2}}{0.1}
\right)^2\right)
:type node_points: :class:`numpy.ndarray`
:param node_points: Points at which evaluate the initial data function.
:rtype: :class:`numpy.ndarray`
:returns: The :math:`u`-values at each node point.
"""
return MathProvider.exp_func(- 25 * (2 * node_points - 1)**2)
[docs]class DG1Solver(object):
r"""Discontinuous Galerkin (DG) solver for the 1D conservation law.
.. math::
\frac{\partial u}{\partial t} + \frac{\partial u}{\partial x} = 0
on the interval :math:`\left[0, 1\right]`. By default, uses
Gaussian-like initial data
.. math::
u(x, 0) = \exp\left(-\left(\frac{x - \frac{1}{2}}{0.1}
\right)^2\right)
but :math:`u(x, 0)` can be specified via ``get_initial_data``.
We represent our solution via the :math:`(p + 1) \times n` rectangular
matrix:
.. math::
\mathbf{u} = \left[ \begin{array}{c c c c}
u_0^1 & u_0^2 & \cdots & u_0^n \\
u_1^1 & u_1^2 & \cdots & u_1^n \\
\vdots & \vdots & \ddots & \vdots \\
u_p^1 & u_p^2 & \cdots & u_p^n
\end{array}\right]
where each column represents one of :math:`n` sub-intervals and each row
represents one of the :math:`p + 1` node points within each sub-interval.
:type num_intervals: int
:param num_intervals: The number :math:`n` of intervals to divide
:math:`\left[0, 1\right]` into.
:type p_order: int
:param p_order: The degree of precision for the method.
:type total_time: float
:param total_time: The amount of time to run the solver for (starts at
:math:`t = 0`).
:type dt: float
:param dt: The timestep to use in the solver.
:type get_initial_data: callable
:param get_initial_data: (Optional) The function to use to evaluate
:math:`u(x, 0)` at the points in our solution.
Defaults to
:func:`get_gaussian_like_initial_data`.
:type points_on_ref_int: :data:`function <types.FunctionType>`
:param points_on_ref_int: (Optional) The method used to partition the
reference interval :math:`\left[0, h\right]`
into node points. Defaults to
:func:`get_evenly_spaced_points`.
"""
def __init__(self, num_intervals, p_order, total_time, dt,
get_initial_data=None, points_on_ref_int=None):
self.num_intervals = num_intervals
self.p_order = p_order
self.total_time = total_time
self.dt = dt
self.current_step = 0
self.points_on_ref_int = points_on_ref_int
# Computed values.
self.num_steps = int(round(self.total_time / self.dt))
self.step_size = MathProvider.num_type(1.0) / self.num_intervals
self.node_points = get_node_points(
self.num_intervals, self.p_order,
step_size=self.step_size,
points_on_ref_int=self.points_on_ref_int)
# The solution: u(x, t).
if get_initial_data is None:
get_initial_data = get_gaussian_like_initial_data
self.solution = get_initial_data(self.node_points)
(self.mass_mat,
self.stiffness_mat) = self.get_mass_and_stiffness_matrices()
[docs] def get_mass_and_stiffness_matrices(self):
"""Get the mass and stiffness matrices for the current solver.
Uses pre-computed mass matrix and stiffness matrix for :math:`p = 1`,
:math:`p = 2` and :math:`p = 3` degree polynomials and computes
the matrices on the fly for larger :math:`p`.
Depends on the sub-interval width ``h`` and the order of accuracy
``p_order``.
:rtype: tuple
:returns: Pair of mass and stiffness matrices, both with
``p_order + 1`` rows and columns.
"""
mass_mat, stiffness_mat = find_matrices(
self.p_order, points_on_ref_int=self.points_on_ref_int)
# We are solving on [0, 1] but ``find_matrices`` is
# on [-1, 1], and the mass matrix is translation invariant
# but scales with interval length.
scaled_mass_mat = (self.step_size *
MathProvider.num_type(0.5) * mass_mat)
return scaled_mass_mat, stiffness_mat
[docs] def ode_func(self, u_val):
r"""Compute the right-hand side for the ODE.
When we write
.. math::
M \dot{\mathbf{u}} = K \mathbf{u} +
\left[ \begin{array}{c c c c c}
u_p^2 & u_p^3 & \cdots & u_p^n & u_p^1 \\
0 & 0 & \cdots & 0 & 0 \\
\vdots & \vdots & \ddots & \vdots & \vdots \\
0 & 0 & \cdots & 0 & 0 \\
-u_p^1 & -u_p^2 & \cdots & -u_p^{n-1} & -u_p^n \\
\end{array}\right]
we specify a RHS :math:`f(u)` via solving the system.
:type u_val: :class:`numpy.ndarray`
:param u_val: The input to :math:`f(u)`.
:rtype: :class:`numpy.ndarray`
:returns: The value of the slope function evaluated at ``u_val``.
"""
# NOTE: Assumes ``stiffness_mat`` and ``u_val`` are NumPy arrays.
# In NumPy, MATMUL tries to use CBLAS if the data type is
# NPY_FLOAT/NPY_CFLOAT/NPY_DOUBLE/NPY_CDOUBLE and if not
# it falls back to the naive way of computing the product
# element-wise.
rhs = np.dot(self.stiffness_mat, u_val)
# First we modify
# K u^k --> K u^k - up^k ep
# so we just take the final row of ``u``
# [up^0, up^1, ..., up^{n-1}]
# and subtract it from the last component of ``r``.
rhs[-1, :] -= u_val[-1, :]
# Then we modify
# K u^k - up^k ep --> K u^k - up^k ep + up^{k-1} e0
# with the assumption that up^{-1} = up^{n-1}, i.e. we
# assume the solution is periodic, so we just roll
# the final row of ``u`` around to
# [up^1, ..., up^{n-1}, up^0]
# and add it to the first component of ``r``.
rhs[0, :] += np.roll(u_val[-1, :], shift=1)
return MathProvider.solve(self.mass_mat, rhs)
[docs] def update(self):
r"""Update the solution for a single time step.
We use :meth:`ode_func` to compute :math:`\dot{u} = f(u)` and
pair it with an RK method (:func:`low_storage_rk`) to compute
the updated value.
"""
self.solution = low_storage_rk(self.ode_func, self.solution, self.dt)
self.current_step += 1