From 8006a351cca3d4de48d8dfa40c26a472eb39157d Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Fri, 24 Apr 2020 12:11:33 -0500 Subject: [PATCH 1/2] Add initial advection demo --- experiments/advection.py | 474 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 474 insertions(+) create mode 100644 experiments/advection.py diff --git a/experiments/advection.py b/experiments/advection.py new file mode 100644 index 0000000..f0dc376 --- /dev/null +++ b/experiments/advection.py @@ -0,0 +1,474 @@ +import contextlib +import numpy as np +import numpy.linalg as la +import numpy.polynomial.legendre as leg +import pytest + + +__doc__ = """ +Notation convention for operator shapes +======================================= + +* m - number of elements in the discretization +* n - number of volume degrees of freedom per element +""" + +import functools + + +memoized = functools.lru_cache(maxsize=None) + + +def ortholegvander(x, deg): + """See numpy.polynomial.legendre.legvander(). Uses an orthonormal basis.""" + result = leg.legvander(x, deg) + factors = np.array([np.sqrt((2*n+1)/2) for n in range(0, 1 + deg)]) + return result * factors + + +def ortholegder(c): + """See numpy.polynomial.legendre.legder(). Uses an orthonormal basis.""" + fw_factors = np.array([np.sqrt((2*n+1)/2) for n in range(len(c))]) + derivs = leg.legder(c * fw_factors) + return derivs / fw_factors[:len(derivs)] + + +def ortholegval(x, c): + """See numpy.polynomial.legendre.legval(). Uses an orthonormal basis.""" + factors = np.array([np.sqrt((2*n+1)/2) for n in range(len(c))]) + return leg.legval(x, c * factors) + + +class DGDiscr1D(object): + """A one-dimensional Discontinuous Galerkin discretization.""" + + def __init__(self, left, right, nelements, nnodes): + """ + Inputs: + left - left endpoint + right - right endpoint + nelements - number of discretization panels + nnodes - number of degrees of freedom per panel + """ + self.left = left + self.right = right + self.nelements = nelements + self.nnodes = nnodes + + @property + @memoized + def ref_nodes(self): + """Return reference nodes for a single element. + + Signature: ->(n,) + """ + nodes, _ = leg.leggauss(self.nnodes) + return nodes + + @property + @memoized + def ref_weights(self): + """Return reference quadrature weights for a single element. + + Signature: ->(n,) + """ + _, weights = leg.leggauss(self.nnodes) + return weights + + def zeros(self): + """Return a zero solution. + + Signature: ->(n*m,) + """ + return np.zeros(self.nnodes * self.nelements) + + @property + def h(self): + """Return the element size. + + Signature: ->() + """ + return self.elements[0,1] - self.elements[0,0] + + def nodes(self): + """Return the vector of node coordinates. + + Signature: ->(n*m,) + """ + centers = (self.elements[:,0] + self.elements[:,1]) / 2 + radii = (self.elements[:,1] - self.elements[:,0]) / 2 + return ((self.ref_nodes[:,np.newaxis] * radii) + centers).T.ravel() + + @property + @memoized + def vdm(self): + """Return the elementwise Vandermonde (modal-to-nodal) matrix. + + Signature: ->(n, n) + """ + return ortholegvander(self.ref_nodes, self.nnodes - 1) + + @property + @memoized + def _ref_mass(self): + """Return the (volume) mass matrix for the reference element. + + Signature: ->(n, n) + """ + return la.inv(self.vdm @ self.vdm.T) + + @property + @memoized + def mass(self): + """Return the elementwise volume mass matrix. + + Signature: ->(n, n) + """ + h = (self.right - self.left) / self.nelements + return (h/2) * self._ref_mass + + @property + @memoized + def inv_mass(self): + """Return the inverse of the elementwise volume mass matrix. + + Signature: ->(n, n) + """ + return la.inv(self.mass) + + @property + @memoized + def face_mass(self): + """Return the face mass matrix. + + The face mass matrix combines the effects of applying the face mass + operator on each face and interpolating the output to the volume nodes. + + Signature: ->(n, 2) + """ + return self.interp.T + + @property + @memoized + def diff(self): + """Return the elementwise differentiation matrix. + + Signature: ->(n, n) + """ + VrT = [] + for row in np.eye(self.nnodes): + deriv = ortholegder(row) + VrT.append(ortholegval(self.ref_nodes, deriv)) + Vr = np.vstack(VrT).T + return Vr @ la.inv(self.vdm) + + @property + @memoized + def stiffness(self): + """Return the stiffness matrix. + + Signature: ->(n, n) + """ + return (self._ref_mass @ self.diff) + + @property + @memoized + def interp(self): + """Return the volume-to-face interpolation matrix. + + Signature: ->(2, n) + """ + return ortholegvander([-1, 1], self.nnodes - 1) @ la.inv(self.vdm) + + @property + @memoized + def elements(self): + """Return the list of elements, each given by their left/right boundaries. + + Signature: ->(m, 2) + """ + h = (self.right - self.left) / self.nelements + return np.array(list(zip( + np.linspace(self.left, self.right, self.nelements, endpoint=False), + np.linspace(h + self.left, self.right, self.nelements)))) + + @property + def dg_ops(self): + """Return a context manager yielding a DGOps1D instance. + """ + return contextlib.contextmanager(lambda: (yield DGOps1DRef(self))) + + @property + def normals(self): + """Return the face unit normals. + + Signature: ->(m, 2) + """ + result = np.zeros((self.nelements, 2)) + result[:,0] = -1 + result[:,1] = 1 + return result + + +def interpolate(discr, vec, nodes): + """Return an interpolated solution at *nodes*. + + Input: + discr - a DGDiscr1D instance + vec - vector of nodal values at degrees of freedom + nodes - vector of nodes to interpolate to + + Signature: (m*n,) -> (len(nodes),) + """ + elements = discr.elements + nelements = discr.nelements + nnodes = discr.nnodes + inv_vdm = la.inv(discr.vdm) + + sorter = np.argsort(nodes) + sorted_nodes = nodes[sorter] + result = [] + + indices = np.searchsorted(sorted_nodes, elements) + for i, (start, end) in enumerate(indices): + if i == 0: + start = 0 + elif i == nelements - 1: + end = len(nodes) + + center = (elements[i][0] + elements[i][1]) / 2 + radius = (elements[i][1] - elements[i][0]) / 2 + element_nodes = sorted_nodes[start:end] + mapped_nodes = (element_nodes - center) / radius + + modal_vals = inv_vdm @ vec[i * nnodes:(i + 1) * nnodes] + nodal_vals = ortholegvander(mapped_nodes, nnodes - 1) @ modal_vals + result.append(nodal_vals) + + result = np.hstack(result) + unsorter = np.arange(len(nodes))[sorter] + return result[unsorter] + + +def integrate(discr, soln): + """Return the integral of the solution. + + Signature: (n*m,) -> () + """ + soln = soln.reshape((discr.nelements, discr.nnodes)) + h = discr.elements[0][1] - discr.elements[0][0] + weights = discr.ref_weights * h / 2 + return np.sum(soln * weights) + + +def elementwise(mat, vec): + """Apply a matrix to rows of the input representing per-element + degrees of freedom. + + Inputs: + mat: Shape (a, b) + vec: Shape (c, b) + + Signature: (a, b), (c, b) -> (c, a) + """ + return np.einsum("ij,kj->ki", mat, vec) + + +class AbstractDGOps1D(object): + + def __init__(self, discr): + self.discr = discr + + def interp(self, vec): + """Apply elementwise volume-to-face interpolation. + + Signature: (m, n) -> (m, 2) + """ + raise NotImplementedError + + def inv_mass(self, vec): + """Apply the elementwise inverse mass matrix. + + Signature: (m, n) -> (m, n) + """ + raise NotImplementedError + + def stiffness(self, vec): + """Apply the elementwise stiffness matrix. + + Signature: (m, n) -> (m, n) + """ + raise NotImplementedError + + def face_mass(self, vec): + """Apply the elementwise face mass matrix. + + Signature: (m, 2) -> (m, n) + """ + raise NotImplementedError + + def face_swap(self, vec): + """Swap values at opposite faces. + + Signature: (m, 2) -> (m, 2) + """ + raise NotImplementedError + + +def elementwise(mat, vec): + """Apply a matrix to rows of the input representing per-element + degrees of freedom. + + Inputs: + mat: Shape (a, b) + vec: Shape (c, b) + + Signature: (a, b), (c, b) -> (c, a) + """ + return np.einsum("ij,kj->ki", mat, vec) + + +class DGOps1DRef(AbstractDGOps1D): + """A reference NumPy implementation of the AbstractDGOps1D interface.""" + + def interp(self, vec): + return elementwise(self.discr.interp, vec) + + def inv_mass(self, vec): + return elementwise(self.discr.inv_mass, vec) + + def stiffness(self, vec): + return elementwise(self.discr.stiffness, vec) + + def face_mass(self, vec): + return elementwise(self.discr.face_mass, vec) + + def face_swap(self, vec): + result = np.zeros_like(vec) + result[:,0] = np.roll(vec[:,1], +1) + result[:,1] = np.roll(vec[:,0], -1) + return result + + +class DGOps1D(AbstractDGOps1D): + pass + + +class AdvectionOperator(object): + """A class representing a DG advection operator.""" + + def __init__(self, discr, c, flux_type): + """ + Inputs: + discr: an instance of DGDiscr1D + c: advection speed parameter + flux_type: "upwind" or "central" + """ + self.discr = discr + self.c = c + assert flux_type in ("upwind", "central") + self.flux_type = flux_type + + def weak_flux(self, dg, vec): + """Apply the flux, weak form. + + Inputs: + dg: a DGOps1D instance + vec: vector of nodal values at the faces + + Signature: (m, 2) -> (m, 2) + """ + if self.flux_type == "central": + flux = (vec + dg.face_swap(vec)) / 2 + + elif self.flux_type == "upwind": + swp = dg.face_swap(vec) + if self.c >= 0: + flux = np.stack((vec[:,0], swp[:,1]), axis=1) + else: + flux = np.stack((swp[:,0], vec[:,1]), axis=1) + + flux *= self.c + flux *= self.discr.normals + + return flux + + def strong_flux(self, dg, vec): + """Apply the flux, strong form. + + Inputs: + dg: a DGOps1D instance + vec: vector of nodal values at the faces + + Signature: (m, 2) -> (m, 2) + """ + return self.c * self.discr.normals * vec - self.weak_flux(dg, vec) + + def apply(self, vec): + """Main operator implementation. + + Signature: (m, n) -> (m, n) + """ + with self.discr.dg_ops() as dg: + pt1 = dg.face_mass(self.strong_flux(dg, dg.interp(vec))) + pt2 = self.c * dg.stiffness(vec) + return -dg.inv_mass(pt1 - pt2) + + def __call__(self, vec): + """Apply the DG advection operator to the vector of degrees of freedom. + + Signature: (m*n,) -> (m*n,) + """ + vec = vec.reshape((self.discr.nelements, self.discr.nnodes)) + return self.apply(vec).reshape((-1,)) + + +def rk4(rhs, initial, t_initial, t_final, dt): + """RK4 integrator. + + Inputs: + - rhs: a callable that takes arguments (t, y) + - initial: initial value + - t_initial: initial time + - t_final: final time + - dt: step size + + Returns: + The solution computed at the final time. + """ + t = t_initial + sol = initial + + while t < t_final: + dt = min(dt, t_final - t) + s0 = rhs(t, sol) + s1 = rhs(t + dt/2, sol + dt/2 * s0) + s2 = rhs(t + dt/2, sol + dt/2 * s1) + s3 = rhs(t + dt, sol + dt * s2) + sol = sol + dt / 6 * (s0 + 2 * s1 + 2 * s2 + s3) + t += dt + + return sol + + +def test_rk4(): + assert np.isclose(rk4(lambda t, y: -y, 1, 0, 5, 0.01), np.exp(-5)) + + +@pytest.mark.parametrize("order", (3, 4, 5)) +@pytest.mark.parametrize("flux_type", ("central", "upwind")) +def test_advection_convergence(order, flux_type): + errors = [] + hs = [] + + for nelements in (8, 12, 16, 20): + discr = DGDiscr1D(0, 2*np.pi, nelements=nelements, nnodes=order) + u_initial = np.sin(discr.nodes()) + op = AdvectionOperator(discr, c=1, flux_type=flux_type) + u = rk4(lambda t, y: op(y), u_initial, t_initial=0, t_final=np.pi, dt=0.01) + u_ref = -u_initial + hs.append(discr.h) + errors.append(integrate(discr, (u - u_ref)**2)**0.5) + + eoc, _ = np.polyfit(np.log(hs), np.log(errors), 1) + assert eoc >= order - 0.1, eoc -- GitLab From 13ae382bbf3c1649de2011a71f237c396e0036df Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Sat, 25 Apr 2020 21:56:52 -0500 Subject: [PATCH 2/2] Run pytest --- experiments/advection.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/experiments/advection.py b/experiments/advection.py index f0dc376..52d6eff 100644 --- a/experiments/advection.py +++ b/experiments/advection.py @@ -472,3 +472,8 @@ def test_advection_convergence(order, flux_type): eoc, _ = np.polyfit(np.log(hs), np.log(errors), 1) assert eoc >= order - 0.1, eoc + + +if __name__ == "__main__": + from pytest import main + main([__file__]) -- GitLab