diff --git a/demos/bbm/demo_bbm_galerkin.py.rst b/demos/bbm/demo_bbm_galerkin.py.rst new file mode 100644 index 00000000..eb9fef67 --- /dev/null +++ b/demos/bbm/demo_bbm_galerkin.py.rst @@ -0,0 +1,112 @@ +Solving the Benjamin-Bona-Mahony equation with Galerkin-in-Time +=============================================================== + +This demo solves the Benjamin-Bona-Mahony equation: + +.. math:: + + u_t + u_x + u u_x - u_{txx} = 0 + +typically posed on :math:`\mathbb{R}` or a bounded interval with periodic +boundaries. + +It is interesting because it is nonlinear (:math:`u u_x`) and also a Sobolev-type equation, with spatial derivatives acting on a time derivative. We can obtain a weak form in the standard way: + +.. math:: + + (u_t, v) + (u_x, v) + (u u_x, v) + (u_{tx}, v_x) = 0 + +BBM is known to have a Hamiltonian structure, and there are three canonical polynomial invariants: + +.. math:: + + I_1 & = \int u \, dx + + I_2 & = \int u^2 + (u_x)^2 \, dx + + I_3 & = \int (u_x)^2 + \tfrac{1}{3} u^3 \, dx + +We are mainly interested in accuracy and in conserving these quantities reasonably well. + + +Firedrake imports:: + + from firedrake import * + from irksome import Dt, GalerkinTimeStepper, MeshConstant + +This function seems to be left out of UFL, but solitary wave solutions for BBM need it:: + + def sech(x): + return 2 / (exp(x) + exp(-x)) + +Set up problem parameters, etc:: + + N = 1000 + L = 100 + h = L / N + msh = PeriodicIntervalMesh(N, L) + + MC = MeshConstant(msh) + t = MC.Constant(0) + dt = MC.Constant(10*h) + + x, = SpatialCoordinate(msh) + +Here is the true solution, which is right-moving solitary wave (but not a soliton):: + + c = Constant(0.5) + + center = 30.0 + delta = -c * center + + uexact = 3 * c**2 / (1-c**2) * sech(0.5 * (c * x - c * t / (1 - c ** 2) + delta))**2 + +Create the approximating space and project true solution:: + + V = FunctionSpace(msh, "CG", 1) + u = project(uexact, V) + v = TestFunction(V) + + F = (inner(Dt(u), v) * dx + + inner(u.dx(0), v) * dx + + inner(u * u.dx(0), v) * dx + + inner((Dt(u)).dx(0), v.dx(0)) * dx) + +For a 1d problem, we don't worry much about efficient solvers.:: + + params = {"mat_type": "aij", + "ksp_type": "preonly", + "pc_type": "lu"} + +The Galerkin-in-Time approach should have symplectic properties:: + + stepper = GalerkinTimeStepper(F, 2, t, dt, u, + solver_parameters=params) + +UFL for the mathematical invariants and containers to track them over time:: + + I1 = u * dx + I2 = (u**2 + (u.dx(0))**2) * dx + I3 = ((u.dx(0))**2 - u**3 / 3) * dx + + I1s = [] + I2s = [] + I3s = [] + + +Time-stepping loop, keeping track of :math:`I` values:: + + tfinal = 18.0 + while (float(t) < tfinal): + if float(t) + float(dt) > tfinal: + dt.assign(tfinal - float(t)) + stepper.advance() + + I1s.append(assemble(I1)) + I2s.append(assemble(I2)) + I3s.append(assemble(I3)) + + print('%.15f %.15f %.15f %.15f' % (float(t), I1s[-1], I2s[-1], I3s[-1])) + t.assign(float(t) + float(dt)) + + print(errornorm(uexact, u) / norm(uexact)) diff --git a/demos/mixed_wave/demo_RTwave_PEP.py.rst b/demos/mixed_wave/demo_RTwave_PEP.py.rst index 60f9258c..6207e495 100644 --- a/demos/mixed_wave/demo_RTwave_PEP.py.rst +++ b/demos/mixed_wave/demo_RTwave_PEP.py.rst @@ -1,5 +1,5 @@ -Solving the Mixed Wave Equation: Energy conservation -==================================================== +Solving the Mixed Wave Equation: Energy conservation and explicit RK +==================================================================== Let :math:`\Omega` be the unit square with boundary :math:`\Gamma`. We write the wave equation as a first-order system of PDE: diff --git a/demos/mixed_wave/demo_RTwave_galerkin.py.rst b/demos/mixed_wave/demo_RTwave_galerkin.py.rst new file mode 100644 index 00000000..591143cb --- /dev/null +++ b/demos/mixed_wave/demo_RTwave_galerkin.py.rst @@ -0,0 +1,130 @@ +Solving the Mixed Wave Equation: Energy conservation, Multigrid, Galerkin-in-Time +================================================================================= + +Let :math:`\Omega` be the unit square with boundary :math:`\Gamma`. We write +the wave equation as a first-order system of PDE: + +.. math:: + + u_t + \nabla p & = 0 + + p_t + \nabla \cdot u & = 0 + +together with homogeneous Dirichlet boundary conditions + +.. math:: + + p = 0 \quad \textrm{on}\ \Gamma + +In this form, at each time, :math:`u` is a vector-valued function in the Soboleve space :math:`H(\mathrm{div})` and `p` is a scalar-valued function. If we select appropriate test functions :math:`v` and :math:`w`, then we can arrive at the weak form + +.. math:: + + (u_t, v) - (p, \nabla \cdot v) & = 0 + + (p_t, w) + (\nabla \cdot u, w) & = 0 + +Note that in mixed formulations, the Dirichlet boundary condition is weakly +enforced via integration by parts rather than strongly in the definition of +the approximating space. + +In this example, we will use the next-to-lowest order Raviart-Thomas elements +for the velocity variable :math:`u` and discontinuous piecewise linear +polynomials for the scalar variable :math:`p`. + +Here is some typical Firedrake boilerplate and the construction of a simple +mesh and the approximating spaces. We are going to use a multigrid preconditioner for each timestep, so we create a MeshHierarchy as well:: + + from firedrake import * + from irksome import Dt, MeshConstant, GalerkinTimeStepper + + N = 10 + + base = UnitSquareMesh(N, N) + mh = MeshHierarchy(base, 2) + msh = mh[-1] + V = FunctionSpace(msh, "RT", 2) + W = FunctionSpace(msh, "DG", 1) + Z = V*W + +Now we can build the initial condition, which has zero velocity and a sinusoidal displacement:: + + x, y = SpatialCoordinate(msh) + up0 = project(as_vector([0, 0, sin(pi*x)*sin(pi*y)]), Z) + u0, p0 = split(up0) + + +We build the variational form in UFL:: + + v, w = TestFunctions(Z) + F = inner(Dt(u0), v)*dx + inner(div(u0), w) * dx + inner(Dt(p0), w)*dx - inner(p0, div(v)) * dx + +Energy conservation is an important principle of the wave equation, and we can +test how well the spatial discretization conserves energy by creating a +UFL expression and evaluating it at each time step:: + + E = 0.5 * (inner(u0, u0)*dx + inner(p0, p0)*dx) + +The time and time step variables:: + + MC = MeshConstant(msh) + t = MC.Constant(0.0) + dt = MC.Constant(1.0/N) + +Here, we experiment with a multigrid preconditioner for the CG(2)-in-time discretization:: + + mgparams = {"mat_type": "aij", + "snes_type": "ksponly", + "ksp_type": "fgmres", + "ksp_rtol": 1e-8, + "pc_type": "mg", + "mg_levels": { + "ksp_type": "chebyshev", + "ksp_max_it": 2, + "ksp_convergence_test": "skip", + "pc_type": "python", + "pc_python_type": "firedrake.PatchPC", + "patch": { + "pc_patch": { + "save_operators": True, + "partition_of_unity": False, + "construct_type": "star", + "construct_dim": 0, + "sub_mat_type": "seqdense", + "dense_inverse": True, + "precompute_element_tensors": None}, + "sub": { + "ksp_type": "preonly", + "pc_type": "lu"}}}, + "mg_coarse": { + "pc_type": "lu", + "pc_factor_mat_solver_type": "mumps"} + } + + + stepper = GalerkinTimeStepper(F, 2, t, dt, up0, + solver_parameters=mgparams) + + +And, as with the heat equation, our time-stepping logic is quite simple. At easch time step, we print out the energy in the system:: + + print("Time Energy") + print("==============") + + while (float(t) < 1.0): + if float(t) + float(dt) > 1.0: + dt.assign(1.0 - float(t)) + + stepper.advance() + + t.assign(float(t) + float(dt)) + print("{0:1.1e} {1:5e}".format(float(t), assemble(E))) + +If all is well with the world, the energy will be nearly identical (up +to roundoff error) at each time step because the Galerkin-in-time methods +are symplectic and applied to a linear Hamiltonian system. + +We can also confirm that the multigrid preconditioner is effective, by computing the average number of linear iterations per time-step:: + + (steps, nl_its, linear_its) = stepper.solver_stats() + print(f"The average number of multigrid iterations per time-step is {linear_its/steps}.") diff --git a/demos/preconditioning/demo_heat_mg_dg.py.rst b/demos/preconditioning/demo_heat_mg_dg.py.rst new file mode 100644 index 00000000..efaf456b --- /dev/null +++ b/demos/preconditioning/demo_heat_mg_dg.py.rst @@ -0,0 +1,123 @@ +Solving the heat equation with monolithic multigrid and Discontinuous Galerkin-in-Time +====================================================================================== + +This reprise of the heat equation demo uses a monolithic multigrid +algorithm suggested by Patrick Farrell to perform time advancement, +but now applies it to the Discontinuous Galerkin-in-Time discretization. + +We consider the heat equation on :math:`\Omega = [0,10] +\times [0,10]`, with boundary :math:`\Gamma`: giving rise to the weak form + +.. math:: + + (u_t, v) + (\nabla u, \nabla v) = (f, v) + +This demo implements an example used by Solin with a particular choice +of :math:`f` given below + +We perform similar imports and setup as before:: + + from firedrake import * + from irksome import DiscontinuousGalerkinTimeStepper, Dt, MeshConstant + from ufl.algorithms.ad import expand_derivatives + + +However, we need to set up a mesh hierarchy to enable geometric multigrid +within Firedrake:: + + N = 128 + x0 = 0.0 + x1 = 10.0 + y0 = 0.0 + y1 = 10.0 + + from math import log + coarseN = 8 # size of coarse grid + nrefs = log(N/coarseN, 2) + assert nrefs == int(nrefs) + nrefs = int(nrefs) + base = RectangleMesh(coarseN, coarseN, x1, y1) + mh = MeshHierarchy(base, nrefs) + msh = mh[-1] + +From here, setting up the function space, manufactured solution, etc, +are just as for the regular heat equation demo:: + + V = FunctionSpace(msh, "CG", 1) + + MC = MeshConstant(msh) + dt = MC.Constant(10.0 / N) + t = MC.Constant(0.0) + + x, y = SpatialCoordinate(msh) + S = Constant(2.0) + C = Constant(1000.0) + B = (x-Constant(x0))*(x-Constant(x1))*(y-Constant(y0))*(y-Constant(y1))/C + R = (x * x + y * y) ** 0.5 + uexact = B * atan(t)*(pi / 2.0 - atan(S * (R - t))) + rhs = expand_derivatives(diff(uexact, t)) - div(grad(uexact)) + + u = Function(V) + u.interpolate(uexact) + v = TestFunction(V) + F = inner(Dt(u), v)*dx + inner(grad(u), grad(v))*dx - inner(rhs, v)*dx + bc = DirichletBC(V, 0, "on_boundary") + +And now for the solver parameters. Note that we are solving a +block-wise system with all stages coupled together. This performs a +monolithic multigrid with pointwise block Jacobi preconditioning:: + + mgparams = {"mat_type": "aij", + "snes_type": "ksponly", + "ksp_type": "gmres", + "ksp_monitor_true_residual": None, + "pc_type": "mg", + "mg_levels": { + "ksp_type": "chebyshev", + "ksp_max_it": 1, + "ksp_convergence_test": "skip", + "pc_type": "python", + "pc_python_type": "firedrake.PatchPC", + "patch": { + "pc_patch": { + "save_operators": True, + "partition_of_unity": False, + "construct_type": "star", + "construct_dim": 0, + "sub_mat_type": "seqdense", + "dense_inverse": True, + "precompute_element_tensors": None}, + "sub": { + "ksp_type": "preonly", + "pc_type": "lu"}}}, + "mg_coarse": { + "pc_type": "lu", + "pc_factor_mat_solver_type": "mumps"} + } + +These solver parameters work just fine in the :class:`.DiscontinuousGalerkinTimeStepper`:: + + stepper = DiscontinuousGalerkinTimeStepper(F, 2, t, dt, u, bcs=bc, + solver_parameters=mgparams) + +And we can advance the solution in time in typical fashion:: + + while (float(t) < 1.0): + if (float(t) + float(dt) > 1.0): + dt.assign(1.0 - float(t)) + stepper.advance() + print(float(t), flush=True) + t.assign(float(t) + float(dt)) + +After the solve, we can retrieve some statistics about the solver:: + + steps, nonlinear_its, linear_its = stepper.solver_stats() + + print("Total number of timesteps was %d" % (steps)) + print("Average number of nonlinear iterations per timestep was %.2f" % (nonlinear_its/steps)) + print("Average number of linear iterations per timestep was %.2f" % (linear_its/steps)) + +Finally, we print out the relative :math:`L^2` error:: + + print() + print(norm(u-uexact)/norm(uexact)) diff --git a/docs/source/index.rst b/docs/source/index.rst index 66f37728..c46a02ff 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -66,6 +66,7 @@ to deploy more efficient methods: demos/demo_heat_pc.py demos/demo_heat_mg.py + demos/demo_heat_mg_dg.py demos/demo_heat_rana.py demos/demo_dirk_parameters.py @@ -76,6 +77,13 @@ We now have support for DIRKs: demos/demo_heat_dirk.py +and for Galerkin-in-Time: + +.. toctree:: + :maxdepth: 1 + + demos/demo_RTwave_galerkin.py + and for explicit schemes: .. toctree:: @@ -109,13 +117,20 @@ Or check out two IMEX-type methods for the monodomain equations: Advanced demos -------------- -There's also an example solving a Sobolev-type equation +There's also an example solving a Sobolev-type equation with symplectic RK methods: .. toctree:: :maxdepth: 1 demos/demo_bbm.py +and with a Galerkin-in-Time approach: + +.. toctree:: + :maxdepth: 1 + + demos/demo_bbm_galerkin.py + Finally, if you feel you must bypass the :py:class:`.TimeStepper` abstraction, we have some examples how to interact with Irksome at a slightly lower level: diff --git a/irksome/__init__.py b/irksome/__init__.py index f69722cd..3bbc92da 100644 --- a/irksome/__init__.py +++ b/irksome/__init__.py @@ -22,3 +22,5 @@ from .stepper import TimeStepper # noqa: F401 from .tools import MeshConstant # noqa: F401 from .wso_dirk_tableaux import WSODIRK # noqa: F401 +from .galerkin_stepper import GalerkinTimeStepper # noqa: F401 +from .discontinuous_galerkin_stepper import DiscontinuousGalerkinTimeStepper # noqa: F401 diff --git a/irksome/discontinuous_galerkin_stepper.py b/irksome/discontinuous_galerkin_stepper.py new file mode 100644 index 00000000..8443cca7 --- /dev/null +++ b/irksome/discontinuous_galerkin_stepper.py @@ -0,0 +1,296 @@ +from functools import reduce +from FIAT import (Bernstein, DiscontinuousElement, + DiscontinuousLagrange, + IntegratedLegendre, Lagrange, + make_quadrature, ufc_simplex) +from operator import mul +from ufl.classes import Zero +from ufl.constantvalue import as_ufl +from .bcs import stage2spaces4bc +from .manipulation import extract_terms, strip_dt_form +from .stage import getBits +from .tools import getNullspace, replace +import numpy as np +from firedrake import as_vector, Constant, dot, TestFunction, Function, NonlinearVariationalProblem as NLVP, NonlinearVariationalSolver as NLVS +from firedrake.dmhooks import pop_parent, push_parent + + +def getFormDiscGalerkin(F, L, Q, t, dt, u0, bcs=None, nullspace=None): + + """Given a time-dependent variational form, trial and test spaces, and + a quadrature rule, produce UFL for the Discontinuous Galerkin-in-Time method. + + :arg F: UFL form for the semidiscrete ODE/DAE + :arg L: A :class:`FIAT.FiniteElement` for the test and trial functions in time + :arg Q: A :class:`FIAT.QuadratureRule` for the time integration + :arg t: a :class:`Function` on the Real space over the same mesh as + `u0`. This serves as a variable referring to the current time. + :arg dt: a :class:`Function` on the Real space over the same mesh as + `u0`. This serves as a variable referring to the current time step. + The user may adjust this value between time steps. + :arg u0: a :class:`Function` referring to the state of + the PDE system at time `t` + :arg bcs: optionally, a :class:`DirichletBC` object (or iterable thereof) + containing (possibly time-dependent) boundary conditions imposed + on the system. + :arg nullspace: A list of tuples of the form (index, VSB) where + index is an index into the function space associated with `u` + and VSB is a :class: `firedrake.VectorSpaceBasis` instance to + be passed to a `firedrake.MixedVectorSpaceBasis` over the + larger space associated with the Runge-Kutta method + + On output, we return a tuple consisting of four parts: + + - Fnew, the :class:`Form` corresponding to the DG-in-Time discretized problem + - UU, the :class:`Function` representing the stages to be solved for + - `bcnew`, a list of :class:`firedrake.DirichletBC` objects to be posed + on the Galerkin-in-time solution, + - 'nspnew', the :class:`firedrake.MixedVectorSpaceBasis` object + that represents the nullspace of the coupled system + """ + assert Q.ref_el.get_spatial_dimension() == 1 + assert L.get_reference_element() == Q.ref_el + + v = F.arguments()[0] + V = v.function_space() + assert V == u0.function_space() + + vecconst = Constant + + num_fields = len(V) + num_stages = L.space_dimension() + + Vbig = reduce(mul, (V for _ in range(num_stages))) + + VV = TestFunction(Vbig) + UU = Function(Vbig) + + qpts = Q.get_points() + qwts = Q.get_weights() + + tabulate_basis = L.tabulate(1, qpts) + basis_vals = tabulate_basis[(0,)] + basis_dvals = tabulate_basis[(1,)] + + element = L + if isinstance(element, DiscontinuousElement): + element = element._element + # sort dofs geometrically by entity location + edofs = element.entity_dofs() + perm = [*edofs[0][0], *edofs[1][0], *edofs[0][1]] + basis_vals = basis_vals[perm] + basis_dvals = basis_dvals[perm] + + # mass matrix later for BC + mmat = np.multiply(basis_vals, qwts) @ basis_vals.T + + # L2 projector + proj = Constant(np.linalg.solve(mmat, np.multiply(basis_vals, qwts))) + + u0bits, vbits, VVbits, UUbits = getBits(num_stages, num_fields, + u0, UU, v, VV) + + split_form = extract_terms(F) + dtless = strip_dt_form(split_form.time) + + Fnew = Zero() + + basis_vals = vecconst(basis_vals) + basis_dvals = vecconst(basis_dvals) + + qpts = vecconst(qpts.reshape((-1,))) + qwts = vecconst(qwts) + + # Terms with time derivatives + for i in range(num_stages): + repl = {} + for j in range(num_fields): + repl[vbits[j]] = VVbits[i][j] + for ii in np.ndindex(vbits[j].ufl_shape): + repl[vbits[j][ii]] = VVbits[i][j][ii] + F_i = replace(dtless, repl) + + # now loop over quadrature points + for q in range(len(qpts)): + repl = {t: t + dt * qpts[q]} + for k in range(num_fields): + d_tosub = sum(basis_dvals[ell, q] * UUbits[ell][k] for ell in range(num_stages)) / dt + repl[u0bits[k]] = d_tosub + + for ii in np.ndindex(u0bits[k].ufl_shape): + d_tosub = sum(basis_dvals[ell, q] * UUbits[ell][k][ii] + for ell in range(num_stages)) / dt + repl[u0bits[k][ii]] = d_tosub + Fnew += dt * qwts[q] * basis_vals[i, q] * replace(F_i, repl) + + # jump terms + repl = {} + for k in range(num_fields): + repl[u0bits[k]] = UUbits[0][k] - u0bits[k] + repl[vbits[k]] = VVbits[0][k] + for ii in np.ndindex(u0bits[k].ufl_shape): + repl[u0bits[k][ii]] = UUbits[0][k][ii] - u0bits[k][ii] + repl[vbits[k][ii]] = VVbits[0][k][ii] + Fnew += replace(dtless, repl) + + # handle the rest of the terms + for i in range(num_stages): + repl = {} + for j in range(num_fields): + repl[vbits[j]] = VVbits[i][j] + for ii in np.ndindex(vbits[j].ufl_shape): + repl[vbits[j][ii]] = VVbits[i][j][ii] + F_i = replace(split_form.remainder, repl) + + # now loop over quadrature points + for q in range(len(qpts)): + repl = {t: t + dt * qpts[q]} + for k in range(num_fields): + tosub = sum(basis_vals[ell, q] * UUbits[ell][k] for ell in range(num_stages)) + repl[u0bits[k]] = tosub + + for ii in np.ndindex(u0bits[k].ufl_shape): + tosub = sum(basis_vals[ell, q] * UUbits[ell][k][ii] + for ell in range(num_stages)) + repl[u0bits[k][ii]] = tosub + Fnew += dt * qwts[q] * basis_vals[i, q] * replace(F_i, repl) + + # Oh, honey, is it the boundary conditions? + if bcs is None: + bcs = [] + bcsnew = [] + for bc in bcs: + bcarg = as_ufl(bc._original_arg) + bcblah_at_qp = np.zeros((len(qpts),), dtype="O") + for q in range(len(qpts)): + tcur = t + qpts[q] * dt + bcblah_at_qp[q] = replace(bcarg, {t: tcur}) + bc_func_for_stages = dot(proj, as_vector(bcblah_at_qp)) + for i in range(num_stages): + Vbigi = stage2spaces4bc(bc, V, Vbig, i) + bcsnew.append(bc.reconstruct(V=Vbigi, g=bc_func_for_stages[i])) + + return Fnew, UU, bcsnew, getNullspace(V, Vbig, num_stages, nullspace) + + +class DiscontinuousGalerkinTimeStepper: + """Front-end class for advancing a time-dependent PDE via a Discontinuous Galerkin + in time method + + :arg F: A :class:`ufl.Form` instance describing the semi-discrete problem + F(t, u; v) == 0, where `u` is the unknown + :class:`firedrake.Function and `v` is the + :class:firedrake.TestFunction`. + :arg order: an integer indicating the order of the DG space to use + (with order == 0 corresponding to DG(0)-in-time) + :arg t: a :class:`Function` on the Real space over the same mesh as + `u0`. This serves as a variable referring to the current time. + :arg dt: a :class:`Function` on the Real space over the same mesh as + `u0`. This serves as a variable referring to the current time step. + The user may adjust this value between time steps. + :arg u0: A :class:`firedrake.Function` containing the current + state of the problem to be solved. + :arg bcs: An iterable of :class:`firedrake.DirichletBC` containing + the strongly-enforced boundary conditions. Irksome will + manipulate these to obtain boundary conditions for each + stage of the method. + :arg basis_type: A string indicating the finite element family (either + `'Lagrange'` or `'Bernstein'`) or the Lagrange variant for the + test/trial spaces. Defaults to equispaced Lagrange elements. + :arg quadrature: A :class:`FIAT.QuadratureRule` indicating the quadrature + to be used in time, defaulting to GL with order+1 points + :arg solver_parameters: A :class:`dict` of solver parameters that + will be used in solving the algebraic problem associated + with each time step. + :arg appctx: An optional :class:`dict` containing application context. + This gets included with particular things that Irksome will + pass into the nonlinear solver so that, say, user-defined preconditioners + have access to it. + :arg nullspace: A list of tuples of the form (index, VSB) where + index is an index into the function space associated with + `u` and VSB is a :class: `firedrake.VectorSpaceBasis` + instance to be passed to a + `firedrake.MixedVectorSpaceBasis` over the larger space + associated with the Runge-Kutta method + """ + def __init__(self, F, order, t, dt, u0, bcs=None, basis_type=None, + quadrature=None, + solver_parameters=None, appctx=None, nullspace=None): + assert order >= 0 + self.u0 = u0 + self.orig_bcs = bcs + self.t = t + self.dt = dt + self.order = order + self.basis_type = basis_type + + V = u0.function_space() + self.num_fields = len(V) + + ufc_line = ufc_simplex(1) + + if order == 0: + self.el = DiscontinuousLagrange(ufc_line, 0) + elif basis_type == "Bernstein": + self.el = DiscontinuousElement(Bernstein(ufc_line, order)) + elif basis_type == "integral": + self.el = DiscontinuousElement(IntegratedLegendre(ufc_line, order)) + else: + # Let recursivenodes handle the general case + variant = None if basis_type == "Lagrange" else basis_type + self.el = DiscontinuousElement(Lagrange(ufc_line, order, variant=variant)) + + if quadrature is None: + quadrature = make_quadrature(ufc_line, order+1) + self.quadrature = quadrature + assert np.size(quadrature.get_points()) >= order+1 + + self.num_steps = 0 + self.num_nonlinear_iterations = 0 + self.num_linear_iterations = 0 + + bigF, UU, bigBCs, bigNSP = \ + getFormDiscGalerkin(F, self.el, + quadrature, t, dt, u0, bcs, nullspace) + + self.UU = UU + self.bigBCs = bigBCs + problem = NLVP(bigF, UU, bigBCs) + appctx_irksome = {"F": F, + "t": t, + "dt": dt, + "u0": u0, + "bcs": bcs, + "nullspace": nullspace} + if appctx is None: + appctx = appctx_irksome + else: + appctx = {**appctx, **appctx_irksome} + + push_parent(u0.function_space().dm, UU.function_space().dm) + self.solver = NLVS(problem, + appctx=appctx, + solver_parameters=solver_parameters, + nullspace=bigNSP) + pop_parent(u0.function_space().dm, UU.function_space().dm) + + def advance(self): + """Advances the system from time `t` to time `t + dt`. + Note: overwrites the value `u0`.""" + push_parent(self.u0.function_space().dm, self.UU.function_space().dm) + self.solver.solve() + pop_parent(self.u0.function_space().dm, self.UU.function_space().dm) + + self.num_steps += 1 + self.num_nonlinear_iterations += self.solver.snes.getIterationNumber() + self.num_linear_iterations += self.solver.snes.getLinearSolveIterations() + + u0 = self.u0 + u0bits = u0.subfunctions + UUs = self.UU.subfunctions + + for i, u0bit in enumerate(u0bits): + u0bit.assign(UUs[self.num_fields*(self.order)+i]) + + def solver_stats(self): + return (self.num_steps, self.num_nonlinear_iterations, self.num_linear_iterations) diff --git a/irksome/galerkin_stepper.py b/irksome/galerkin_stepper.py new file mode 100644 index 00000000..48709f3f --- /dev/null +++ b/irksome/galerkin_stepper.py @@ -0,0 +1,272 @@ +from functools import reduce +from FIAT import (Bernstein, DiscontinuousElement, DiscontinuousLagrange, + IntegratedLegendre, Lagrange, Legendre, + make_quadrature, ufc_simplex) +from operator import mul +from ufl.classes import Zero +from ufl.constantvalue import as_ufl +from .bcs import bc2space, stage2spaces4bc +from .deriv import TimeDerivative +from .stage import getBits +from .tools import getNullspace, replace +import numpy as np +from firedrake import as_vector, dot, Constant, TestFunction, Function, NonlinearVariationalProblem as NLVP, NonlinearVariationalSolver as NLVS +from firedrake.dmhooks import pop_parent, push_parent + + +def getFormGalerkin(F, L_trial, L_test, Q, t, dt, u0, bcs=None, nullspace=None): + + """Given a time-dependent variational form, trial and test spaces, and + a quadrature rule, produce UFL for the Galerkin-in-Time method. + + :arg F: UFL form for the semidiscrete ODE/DAE + :arg L_trial: A :class:`FIAT.FiniteElement` for the trial functions in time + :arg L_test: A :class:`FIAT.FinteElement` for the test functions in time + :arg Q: A :class:`FIAT.QuadratureRule` for the time integration + :arg t: a :class:`Function` on the Real space over the same mesh as + `u0`. This serves as a variable referring to the current time. + :arg dt: a :class:`Function` on the Real space over the same mesh as + `u0`. This serves as a variable referring to the current time step. + The user may adjust this value between time steps. + :arg u0: a :class:`Function` referring to the state of + the PDE system at time `t` + :arg bcs: optionally, a :class:`DirichletBC` object (or iterable thereof) + containing (possibly time-dependent) boundary conditions imposed + on the system. + :arg nullspace: A list of tuples of the form (index, VSB) where + index is an index into the function space associated with `u` + and VSB is a :class: `firedrake.VectorSpaceBasis` instance to + be passed to a `firedrake.MixedVectorSpaceBasis` over the + larger space associated with the Runge-Kutta method + + On output, we return a tuple consisting of four parts: + + - Fnew, the :class:`Form` corresponding to the Galerkin-in-Time discretized problem + - UU, the :class:`Function` representing the stages to be solved for + - `bcnew`, a list of :class:`firedrake.DirichletBC` objects to be posed + on the Galerkin-in-time solution, + - 'nspnew', the :class:`firedrake.MixedVectorSpaceBasis` object + that represents the nullspace of the coupled system + """ + assert L_test.get_reference_element() == Q.ref_el + assert L_trial.get_reference_element() == Q.ref_el + assert Q.ref_el.get_spatial_dimension() == 1 + assert L_trial.get_order() == L_test.get_order() + 1 + + v = F.arguments()[0] + V = v.function_space() + assert V == u0.function_space() + + num_fields = len(V) + num_stages = L_test.space_dimension() + + Vbig = reduce(mul, (V for _ in range(num_stages))) + VV = TestFunction(Vbig) + UU = Function(Vbig) # u0 + this are coefficients of the Galerkin polynomial + + qpts = Q.get_points() + qwts = Q.get_weights() + + tabulate_trials = L_trial.tabulate(1, qpts) + trial_vals = tabulate_trials[(0,)] + trial_dvals = tabulate_trials[(1,)] + test_vals = L_test.tabulate(0, qpts)[(0,)] + + # sort dofs geometrically by entity location + edofs = L_trial.entity_dofs() + trial_perm = [*edofs[0][0], *edofs[1][0], *edofs[0][1]] + trial_vals = trial_vals[trial_perm] + trial_dvals = trial_dvals[trial_perm] + + # mass-ish matrix later for BC + mmat = np.multiply(test_vals, qwts) @ trial_vals[1:].T + + # L2 projector + proj = Constant(np.linalg.solve(mmat, np.multiply(test_vals, qwts))) + + u0bits, vbits, VVbits, UUbits = getBits(num_stages, num_fields, + u0, UU, v, VV) + + Fnew = Zero() + + trial_vals = Constant(trial_vals) + trial_dvals = Constant(trial_dvals) + test_vals = Constant(test_vals) + qpts = Constant(qpts.reshape((-1,))) + qwts = Constant(qwts) + + for i in range(num_stages): + repl = {} + for j in range(num_fields): + repl[vbits[j]] = VVbits[i][j] + for ii in np.ndindex(vbits[j].ufl_shape): + repl[vbits[j][ii]] = VVbits[i][j][ii] + F_i = replace(F, repl) + + # now loop over quadrature points + for q in range(len(qpts)): + repl = {t: t + dt * qpts[q]} + for k in range(num_fields): + tosub = u0bits[k] * trial_vals[0, q] + d_tosub = u0bits[k] * trial_dvals[0, q] + for ell in range(num_stages): + tosub += UUbits[ell][k] * trial_vals[1 + ell, q] + d_tosub += UUbits[ell][k] * trial_dvals[1 + ell, q] + + repl[u0bits[k]] = tosub + repl[TimeDerivative(u0bits[k])] = d_tosub / dt + + for ii in np.ndindex(u0bits[k].ufl_shape): + tosub = u0bits[k][ii] * trial_vals[0, q] + d_tosub = u0bits[k][ii] * trial_dvals[0, q] + for ell in range(num_stages): + tosub += UUbits[ell][k][ii] * trial_vals[1 + ell, q] + d_tosub += UUbits[ell][k][ii] * trial_dvals[1 + ell, q] + repl[u0bits[k][ii]] = tosub + repl[TimeDerivative(u0bits[k][ii])] = d_tosub / dt + Fnew += dt * qwts[q] * test_vals[i, q] * replace(F_i, repl) + + # Oh, honey, is it the boundary conditions? + if bcs is None: + bcs = [] + bcsnew = [] + for bc in bcs: + u0_sub = bc2space(bc, u0) + bcarg = as_ufl(bc._original_arg) + bcblah_at_qp = np.zeros((len(qpts),), dtype="O") + for q in range(len(qpts)): + tcur = t + qpts[q] * dt + bcblah_at_qp[q] = replace(bcarg, {t: tcur}) - u0_sub * trial_vals[0, q] + bc_func_for_stages = dot(proj, as_vector(bcblah_at_qp)) + for i in range(num_stages): + Vbigi = stage2spaces4bc(bc, V, Vbig, i) + bcsnew.append(bc.reconstruct(V=Vbigi, g=bc_func_for_stages[i])) + + return Fnew, UU, bcsnew, getNullspace(V, Vbig, num_stages, nullspace) + + +class GalerkinTimeStepper: + """Front-end class for advancing a time-dependent PDE via a Galerkin + in time method + + :arg F: A :class:`ufl.Form` instance describing the semi-discrete problem + F(t, u; v) == 0, where `u` is the unknown + :class:`firedrake.Function and `v` is the + :class:firedrake.TestFunction`. + :arg order: an integer indicating the order of the DG space to use + (with order == 1 corresponding to CG(1)-in-time for the trial space) + :arg t: a :class:`Function` on the Real space over the same mesh as + `u0`. This serves as a variable referring to the current time. + :arg dt: a :class:`Function` on the Real space over the same mesh as + `u0`. This serves as a variable referring to the current time step. + The user may adjust this value between time steps. + :arg u0: A :class:`firedrake.Function` containing the current + state of the problem to be solved. + :arg bcs: An iterable of :class:`firedrake.DirichletBC` containing + the strongly-enforced boundary conditions. Irksome will + manipulate these to obtain boundary conditions for each + stage of the method. + :arg basis_type: A string indicating the finite element family (either + `'Lagrange'` or `'Bernstein'`) or the Lagrange variant for the + test/trial spaces. Defaults to equispaced Lagrange elements. + :arg quadrature: A :class:`FIAT.QuadratureRule` indicating the quadrature + to be used in time, defaulting to GL with order points + :arg solver_parameters: A :class:`dict` of solver parameters that + will be used in solving the algebraic problem associated + with each time step. + :arg appctx: An optional :class:`dict` containing application context. + This gets included with particular things that Irksome will + pass into the nonlinear solver so that, say, user-defined preconditioners + have access to it. + :arg nullspace: A list of tuples of the form (index, VSB) where + index is an index into the function space associated with + `u` and VSB is a :class: `firedrake.VectorSpaceBasis` + instance to be passed to a + `firedrake.MixedVectorSpaceBasis` over the larger space + associated with the Runge-Kutta method + """ + def __init__(self, F, order, t, dt, u0, bcs=None, basis_type=None, + quadrature=None, + solver_parameters=None, appctx=None, nullspace=None): + assert order >= 1 + self.u0 = u0 + self.orig_bcs = bcs + self.t = t + self.dt = dt + self.order = order + self.basis_type = basis_type + + V = u0.function_space() + self.num_fields = len(V) + + ufc_line = ufc_simplex(1) + if basis_type == "Bernstein": + self.trial_el = Bernstein(ufc_line, order) + if order == 1: + self.test_el = DiscontinuousLagrange(ufc_line, 0) + else: + self.test_el = DiscontinuousElement( + Bernstein(ufc_line, order-1)) + elif basis_type == "integral": + self.trial_el = IntegratedLegendre(ufc_line, order) + self.test_el = Legendre(ufc_line, order-1) + else: + # Let recursivenodes handle the general case + variant = None if basis_type == "Lagrange" else basis_type + self.trial_el = Lagrange(ufc_line, order, variant=variant) + self.test_el = DiscontinuousLagrange(ufc_line, order-1, variant=variant) + + if quadrature is None: + quadrature = make_quadrature(ufc_line, order) + self.quadrature = quadrature + assert np.size(quadrature.get_points()) >= order + + self.num_steps = 0 + self.num_nonlinear_iterations = 0 + self.num_linear_iterations = 0 + + bigF, UU, bigBCs, bigNSP = \ + getFormGalerkin(F, self.trial_el, self.test_el, + quadrature, t, dt, u0, bcs, nullspace) + + self.UU = UU + self.bigBCs = bigBCs + problem = NLVP(bigF, UU, bigBCs) + appctx_irksome = {"F": F, + "t": t, + "dt": dt, + "u0": u0, + "bcs": bcs, + "nullspace": nullspace} + if appctx is None: + appctx = appctx_irksome + else: + appctx = {**appctx, **appctx_irksome} + + push_parent(u0.function_space().dm, UU.function_space().dm) + self.solver = NLVS(problem, + appctx=appctx, + solver_parameters=solver_parameters, + nullspace=bigNSP) + pop_parent(u0.function_space().dm, UU.function_space().dm) + + def advance(self): + """Advances the system from time `t` to time `t + dt`. + Note: overwrites the value `u0`.""" + push_parent(self.u0.function_space().dm, self.UU.function_space().dm) + self.solver.solve() + pop_parent(self.u0.function_space().dm, self.UU.function_space().dm) + + self.num_steps += 1 + self.num_nonlinear_iterations += self.solver.snes.getIterationNumber() + self.num_linear_iterations += self.solver.snes.getLinearSolveIterations() + + u0 = self.u0 + u0bits = u0.subfunctions + UUs = self.UU.subfunctions + + for i, u0bit in enumerate(u0bits): + u0bit.assign(UUs[self.num_fields*(self.order-1)+i]) + + def solver_stats(self): + return (self.num_steps, self.num_nonlinear_iterations, self.num_linear_iterations) diff --git a/irksome/getForm.py b/irksome/getForm.py index 4a8e2b77..1c3c91d3 100644 --- a/irksome/getForm.py +++ b/irksome/getForm.py @@ -33,7 +33,7 @@ def getForm(F, butch, t, dt, u0, bcs=None, bc_type=None, splitting=AI, :arg u0: a :class:`Function` referring to the state of the PDE system at time `t` :arg bcs: optionally, a :class:`DirichletBC` object (or iterable thereof) - containing (possible time-dependent) boundary conditions imposed + containing (possibly time-dependent) boundary conditions imposed on the system. :arg bc_type: How to manipulate the strongly-enforced boundary conditions to derive the stage boundary conditions. Should @@ -181,6 +181,6 @@ def bc2gcur(bc, i): gdat = BCStageData(Vsp, gcur, u0, u0_mult, i, t, dt) bcnew.append(bc.reconstruct(V=Vbigi, g=gdat)) - nspnew = getNullspace(V, Vbig, butch, nullspace) + nspnew = getNullspace(V, Vbig, num_stages, nullspace) return Fnew, w, bcnew, nspnew diff --git a/irksome/imex.py b/irksome/imex.py index f2aed394..04c324b3 100644 --- a/irksome/imex.py +++ b/irksome/imex.py @@ -1,6 +1,6 @@ import FIAT import numpy as np -from firedrake import (Function, NonlinearVariationalProblem, +from firedrake import (Constant, Function, NonlinearVariationalProblem, NonlinearVariationalSolver, TestFunction, as_ufl, dx, inner, split) from firedrake.dmhooks import pop_parent, push_parent @@ -18,19 +18,18 @@ def riia_explicit_coeffs(k): of a RadauIIA-IMEX method.""" U = FIAT.ufc_simplex(1) L = FIAT.GaussRadau(U, k - 1) - Q = FIAT.make_quadrature(L.ref_el, 2*k) - c = np.asarray([list(ell.pt_dict.keys())[0][0] - for ell in L.dual.nodes]) Q = FIAT.make_quadrature(L.ref_el, 2*k) qpts = Q.get_points() qwts = Q.get_weights() A = np.zeros((k, k)) - for i in range(k): - qpts_i = 1 + qpts * c[i] - qwts_i = qwts * c[i] - Lvals_i = L.tabulate(0, qpts_i)[0, ] + for i, ell in enumerate(L.dual.nodes): + pt, = ell.pt_dict + ci, = pt + qpts_i = 1 + qpts * ci + qwts_i = qwts * ci + Lvals_i = L.tabulate(0, qpts_i)[(0,)] A[i, :] = Lvals_i @ qwts_i return A @@ -42,18 +41,15 @@ def getFormExplicit(Fexp, butch, u0, UU, t, dt, splitting=None): which really just differ by which constants are in them.""" v = Fexp.arguments()[0] V = v.function_space() - msh = V.mesh() Vbig = UU.function_space() VV = TestFunction(Vbig) num_stages = butch.num_stages num_fields = len(V) - MC = MeshConstant(msh) - vc = np.vectorize(lambda c: MC.Constant(c)) Aexp = riia_explicit_coeffs(num_stages) - Aprop = vc(Aexp) - Ait = vc(butch.A) - C = vc(butch.c) + Aprop = Constant(Aexp) + Ait = Constant(butch.A) + C = Constant(butch.c) u0bits, vbits, VVbits, UUbits = getBits(num_stages, num_fields, u0, UU, v, VV) @@ -103,7 +99,7 @@ def getFormExplicit(Fexp, butch, u0, UU, t, dt, splitting=None): Fit += dt * replace(Fexp, repl) # dense contribution to propagator - Ablah = vc(np.linalg.solve(butch.A, Aexp)) + Ablah = Constant(np.linalg.solve(butch.A, Aexp)) for i in range(num_stages): # replace test function diff --git a/irksome/stage.py b/irksome/stage.py index 68ebb515..3549627e 100644 --- a/irksome/stage.py +++ b/irksome/stage.py @@ -5,11 +5,11 @@ import numpy as np from FIAT import Bernstein, ufc_simplex -from firedrake import (Function, NonlinearVariationalProblem, +from firedrake import (Constant, + Function, NonlinearVariationalProblem, NonlinearVariationalSolver, TestFunction, dx, inner, split) from firedrake.petsc import PETSc -from numpy import vectorize from ufl.classes import Zero from ufl.constantvalue import as_ufl @@ -25,20 +25,13 @@ def isiterable(x): def split_field(num_fields, u): - return np.array((u,) if num_fields == 1 else split(u), dtype="O") + ubits = np.array(split(u), dtype="O") + return ubits def split_stage_field(num_stages, num_fields, UU): - if num_fields == 1: - if num_stages == 1: # single-stage method - UUbits = np.reshape(np.array((UU,), dtype='O'), (num_stages, num_fields)) - else: # multi-stage methods - UUbits = np.zeros((len(split(UU)),), dtype='O') - for (i, x) in enumerate(split(UU)): - UUbits[i] = np.zeros((1,), dtype='O') - UUbits[i][0] = x - else: - UUbits = np.reshape(np.asarray(split(UU), dtype="O"), (num_stages, num_fields)) + UUbits = np.reshape(np.asarray(split(UU), dtype="O"), + (num_stages, num_fields)) return UUbits @@ -78,7 +71,7 @@ def getFormStage(F, butch, u0, t, dt, bcs=None, splitting=None, vandermonde=None If none is provided, we assume it is the identity, working in the Lagrange basis. :arg bcs: optionally, a :class:`DirichletBC` object (or iterable thereof) - containing (possible time-dependent) boundary conditions imposed + containing (possibly time-dependent) boundary conditions imposed on the system. :arg bc_constraints: optionally, a dictionary mapping (some of) the boundary conditions in `bcs` to triples of the form (params, lower, upper) indicating @@ -129,7 +122,7 @@ def getFormStage(F, butch, u0, t, dt, bcs=None, splitting=None, vandermonde=None u0, ZZ, v, VV) MC = MeshConstant(V.mesh()) - vecconst = np.vectorize(lambda c: MC.Constant(c)) + vecconst = Constant C = vecconst(butch.c) A = vecconst(butch.A) @@ -146,13 +139,7 @@ def getFormStage(F, butch, u0, t, dt, bcs=None, splitting=None, vandermonde=None Vander_col = Vander[1:, 0] Vander0 = Vander[1:, 1:] - v0u0 = np.zeros((num_stages, num_fields), dtype="O") - for i in range(num_stages): - for j in range(num_fields): - v0u0[i, j] = Vander_col[i] * u0bits[j] - - if num_fields == 1: - v0u0 = np.reshape(v0u0, (-1,)) + v0u0 = np.outer(Vander_col, u0bits) UUbits = v0u0 + Vander0 @ ZZbits @@ -174,9 +161,7 @@ def getFormStage(F, butch, u0, t, dt, bcs=None, splitting=None, vandermonde=None for j in range(num_fields): repl[u0bits[j]] = UUbits[i][j] - u0bits[j] repl[vbits[j]] = VVbits[i][j] - - # Also get replacements right for indexing. - for j in range(num_fields): + # Also get replacements right for indexing. for ii in np.ndindex(u0bits[j].ufl_shape): repl[u0bits[j][ii]] = UUbits[i][j][ii] - u0bits[j][ii] repl[vbits[j][ii]] = VVbits[i][j][ii] @@ -208,7 +193,7 @@ def getFormStage(F, butch, u0, t, dt, bcs=None, splitting=None, vandermonde=None Fnew += A[i, j] * dt * replace(Ftmp, repl) elif splitting == IA: - Ainv = np.vectorize(lambda c: MC.Constant(c))(np.linalg.inv(butch.A)) + Ainv = vecconst(np.linalg.inv(butch.A)) # time derivative part gets inverse of Butcher matrix. for i in range(num_stages): @@ -281,7 +266,7 @@ def getFormStage(F, butch, u0, t, dt, bcs=None, splitting=None, vandermonde=None bcsnew.extend(bcnew_cur) - nspacenew = getNullspace(V, Vbig, butch, nullspace) + nspacenew = getNullspace(V, Vbig, num_stages, nullspace) # only form update stuff if we need it # which means neither stiffly accurate nor Vandermonde @@ -289,8 +274,9 @@ def getFormStage(F, butch, u0, t, dt, bcs=None, splitting=None, vandermonde=None unew = Function(V) Fupdate = inner(unew - u0, v) * dx - B = vectorize(lambda c: MC.Constant(c))(butch.b) - C = vectorize(lambda c: MC.Constant(c))(butch.c) + vecconst = np.vectorize(MC.Constant) + B = vecconst(butch.b) + C = vecconst(butch.c) for i in range(num_stages): repl = {t: t + C[i] * dt} diff --git a/irksome/tools.py b/irksome/tools.py index e9bf6c1d..16b1d0a6 100644 --- a/irksome/tools.py +++ b/irksome/tools.py @@ -9,20 +9,19 @@ from irksome.deriv import TimeDerivative -def getNullspace(V, Vbig, butch, nullspace): +def getNullspace(V, Vbig, num_stages, nullspace): """ Computes the nullspace for a multi-stage method. :arg V: The :class:`FunctionSpace` on which the original time-dependent PDE is posed. :arg Vbig: The multi-stage :class:`FunctionSpace` for the stage problem - :arg butch: The :class:`ButcherTableau` defining the RK method + :arg num_stages: The number of stages in the RK method :arg nullspace: The nullspace for the original problem. On output, we produce a :class:`MixedVectorSpaceBasis` defining the nullspace for the multistage problem. """ - num_stages = butch.num_stages num_fields = len(V) if nullspace is None: nspnew = None diff --git a/tests/test_disc_galerkin.py b/tests/test_disc_galerkin.py new file mode 100644 index 00000000..80c45b33 --- /dev/null +++ b/tests/test_disc_galerkin.py @@ -0,0 +1,168 @@ +from math import isclose + +import pytest +from firedrake import * +from irksome import Dt, MeshConstant, DiscontinuousGalerkinTimeStepper +from irksome import TimeStepper, RadauIIA +from ufl.algorithms.ad import expand_derivatives +import FIAT + + +@pytest.mark.parametrize("order", [0, 1, 2]) +@pytest.mark.parametrize("basis_type", ["Lagrange", "Bernstein", "integral"]) +def test_1d_heat_dirichletbc(order, basis_type): + # Boundary values + u_0 = Constant(2.0) + u_1 = Constant(3.0) + + N = 10 + x0 = 0.0 + x1 = 10.0 + msh = IntervalMesh(N, x1) + V = FunctionSpace(msh, "CG", 1) + MC = MeshConstant(msh) + dt = MC.Constant(1.0 / N) + t = MC.Constant(0.0) + (x,) = SpatialCoordinate(msh) + + # Method of manufactured solutions copied from Heat equation demo. + S = Constant(2.0) + C = Constant(1000.0) + B = (x - Constant(x0)) * (x - Constant(x1)) / C + R = (x * x) ** 0.5 + # Note end linear contribution + uexact = ( + B * atan(t) * (pi / 2.0 - atan(S * (R - t))) + + u_0 + + ((x - x0) / x1) * (u_1 - u_0) + ) + rhs = expand_derivatives(diff(uexact, t)) - div(grad(uexact)) + u = Function(V) + u.interpolate(uexact) + v = TestFunction(V) + F = ( + inner(Dt(u), v) * dx + + inner(grad(u), grad(v)) * dx + - inner(rhs, v) * dx + ) + bc = [ + DirichletBC(V, u_1, 2), + DirichletBC(V, u_0, 1), + ] + + luparams = {"mat_type": "aij", "ksp_type": "preonly", "pc_type": "lu"} + + stepper = DiscontinuousGalerkinTimeStepper( + F, order, t, dt, u, bcs=bc, basis_type=basis_type, + solver_parameters=luparams + ) + + t_end = 2.0 + while float(t) < t_end: + if float(t) + float(dt) > t_end: + dt.assign(t_end - float(t)) + stepper.advance() + t.assign(float(t) + float(dt)) + # Check solution and boundary values + assert errornorm(uexact, u) / norm(uexact) < 10.0 ** -3 + assert isclose(u.at(x0), u_0) + assert isclose(u.at(x1), u_1) + + +@pytest.mark.parametrize("order", [0, 1, 2]) +def test_1d_heat_neumannbc(order): + N = 20 + msh = UnitIntervalMesh(N) + V = FunctionSpace(msh, "CG", 1) + MC = MeshConstant(msh) + dt = MC.Constant(1.0 / N) + t = MC.Constant(0.0) + (x,) = SpatialCoordinate(msh) + butcher_tableau = RadauIIA(order+1) + + uexact = cos(pi*x)*exp(-(pi**2)*t) + rhs = expand_derivatives(diff(uexact, t)) - div(grad(uexact)) + u_Radau = Function(V) + u = Function(V) + u_Radau.interpolate(uexact) + u.interpolate(uexact) + + v = TestFunction(V) + F = ( + inner(Dt(u), v) * dx + + inner(grad(u), grad(v)) * dx + - inner(rhs, v) * dx + ) + F_Radau = replace(F, {u: u_Radau}) + + luparams = {"mat_type": "aij", "ksp_type": "preonly", "pc_type": "lu"} + + ufc_line = FIAT.ufc_simplex(1) + quadrature = FIAT.quadrature.RadauQuadratureLineRule(ufc_line, order+1) + + stepper = DiscontinuousGalerkinTimeStepper( + F, order, t, dt, u, quadrature=quadrature, + solver_parameters=luparams + ) + stepper_Radau = TimeStepper( + F_Radau, butcher_tableau, t, dt, u_Radau, solver_parameters=luparams + ) + + t_end = 1.0 + while float(t) < t_end: + if float(t) + float(dt) > t_end: + dt.assign(t_end - float(t)) + stepper.advance() + stepper_Radau.advance() + t.assign(float(t) + float(dt)) + assert (errornorm(u_Radau, u) / norm(u)) < 1.e-10 + + +@pytest.mark.parametrize("order", [1, 2, 3]) +def test_1d_heat_homogeneous_dirichletbc(order): + N = 20 + msh = UnitIntervalMesh(N) + V = FunctionSpace(msh, "CG", 1) + MC = MeshConstant(msh) + dt = MC.Constant(1.0 / N) + t = MC.Constant(0.0) + (x,) = SpatialCoordinate(msh) + butcher_tableau = RadauIIA(order+1) + + uexact = sin(pi*x)*exp(-(pi**2)*t) + rhs = expand_derivatives(diff(uexact, t)) - div(grad(uexact)) + bcs = DirichletBC(V, uexact, "on_boundary") + u_Radau = Function(V) + u = Function(V) + u_Radau.interpolate(uexact) + u.interpolate(uexact) + + v = TestFunction(V) + F = ( + inner(Dt(u), v) * dx + + inner(grad(u), grad(v)) * dx + - inner(rhs, v) * dx + ) + F_Radau = replace(F, {u: u_Radau}) + + luparams = {"mat_type": "aij", "ksp_type": "preonly", "pc_type": "lu"} + + ufc_line = FIAT.ufc_simplex(1) + quadrature = FIAT.quadrature.RadauQuadratureLineRule(ufc_line, order+1) + + stepper = DiscontinuousGalerkinTimeStepper( + F, order, t, dt, u, bcs=bcs, quadrature=quadrature, + solver_parameters=luparams + ) + stepper_Radau = TimeStepper( + F_Radau, butcher_tableau, t, dt, u_Radau, bcs=bcs, solver_parameters=luparams + ) + + t_end = 1.0 + while float(t) < t_end: + if float(t) + float(dt) > t_end: + dt.assign(t_end - float(t)) + stepper.advance() + stepper_Radau.advance() + t.assign(float(t) + float(dt)) + assert (errornorm(u_Radau, u) / norm(u)) < 1.e-10 diff --git a/tests/test_galerkin.py b/tests/test_galerkin.py new file mode 100644 index 00000000..0c4de9ab --- /dev/null +++ b/tests/test_galerkin.py @@ -0,0 +1,218 @@ +from math import isclose + +import pytest +from firedrake import * +from irksome import Dt, MeshConstant, GalerkinTimeStepper +from irksome import TimeStepper, GaussLegendre +from ufl.algorithms.ad import expand_derivatives +from FIAT import make_quadrature, ufc_simplex + + +@pytest.mark.parametrize("order", [1, 2, 3]) +@pytest.mark.parametrize("basis_type", ["Lagrange", "Bernstein", "integral"]) +def test_1d_heat_dirichletbc(order, basis_type): + # Boundary values + u_0 = Constant(2.0) + u_1 = Constant(3.0) + + N = 10 + x0 = 0.0 + x1 = 10.0 + msh = IntervalMesh(N, x1) + V = FunctionSpace(msh, "CG", 1) + MC = MeshConstant(msh) + dt = MC.Constant(1.0 / N) + t = MC.Constant(0.0) + (x,) = SpatialCoordinate(msh) + + # Method of manufactured solutions copied from Heat equation demo. + S = Constant(2.0) + C = Constant(1000.0) + B = (x - Constant(x0)) * (x - Constant(x1)) / C + R = (x * x) ** 0.5 + # Note end linear contribution + uexact = ( + B * atan(t) * (pi / 2.0 - atan(S * (R - t))) + + u_0 + + ((x - x0) / x1) * (u_1 - u_0) + ) + rhs = expand_derivatives(diff(uexact, t)) - div(grad(uexact)) + u = Function(V) + u.interpolate(uexact) + v = TestFunction(V) + F = ( + inner(Dt(u), v) * dx + + inner(grad(u), grad(v)) * dx + - inner(rhs, v) * dx + ) + bc = [ + DirichletBC(V, u_1, 2), + DirichletBC(V, u_0, 1), + ] + + luparams = {"mat_type": "aij", "ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"} + + stepper = GalerkinTimeStepper( + F, order, t, dt, u, bcs=bc, basis_type=basis_type, + solver_parameters=luparams + ) + + t_end = 2.0 + while float(t) < t_end: + if float(t) + float(dt) > t_end: + dt.assign(t_end - float(t)) + stepper.advance() + t.assign(float(t) + float(dt)) + # Check solution and boundary values + assert errornorm(uexact, u) / norm(uexact) < 10.0 ** -3 + assert isclose(u.at(x0), u_0) + assert isclose(u.at(x1), u_1) + + +@pytest.mark.parametrize("order", [1, 2, 3]) +@pytest.mark.parametrize("num_quad_points", [3, 4]) +def test_1d_heat_neumannbc(order, num_quad_points): + N = 20 + msh = UnitIntervalMesh(N) + V = FunctionSpace(msh, "CG", 1) + MC = MeshConstant(msh) + dt = MC.Constant(1.0 / N) + t = MC.Constant(0.0) + (x,) = SpatialCoordinate(msh) + butcher_tableau = GaussLegendre(order) + + uexact = cos(pi*x)*exp(-(pi**2)*t) + rhs = expand_derivatives(diff(uexact, t)) - div(grad(uexact)) + u_GL = Function(V) + u = Function(V) + u_GL.interpolate(uexact) + u.interpolate(uexact) + + v = TestFunction(V) + F = ( + inner(Dt(u), v) * dx + + inner(grad(u), grad(v)) * dx + - inner(rhs, v) * dx + ) + F_GL = replace(F, {u: u_GL}) + + luparams = {"mat_type": "aij", "ksp_type": "preonly", "pc_type": "lu"} + + ufc_line = ufc_simplex(1) + quadrature = make_quadrature(ufc_line, num_quad_points) + + stepper = GalerkinTimeStepper( + F, order, t, dt, u, quadrature=quadrature, + solver_parameters=luparams + ) + stepper_GL = TimeStepper( + F_GL, butcher_tableau, t, dt, u_GL, solver_parameters=luparams + ) + + t_end = 1.0 + while float(t) < t_end: + if float(t) + float(dt) > t_end: + dt.assign(t_end - float(t)) + stepper.advance() + stepper_GL.advance() + t.assign(float(t) + float(dt)) + assert (errornorm(u_GL, u) / norm(u)) < 1.e-10 + + +@pytest.mark.parametrize("order", [1, 2, 3]) +def test_1d_heat_homogeneous_dirichletbc(order): + N = 20 + msh = UnitIntervalMesh(N) + V = FunctionSpace(msh, "CG", 1) + MC = MeshConstant(msh) + dt = MC.Constant(1.0 / N) + t = MC.Constant(0.0) + (x,) = SpatialCoordinate(msh) + butcher_tableau = GaussLegendre(order) + + uexact = sin(pi*x)*exp(-(pi**2)*t) + rhs = expand_derivatives(diff(uexact, t)) - div(grad(uexact)) + bcs = DirichletBC(V, uexact, "on_boundary") + u_GL = Function(V) + u = Function(V) + u_GL.interpolate(uexact) + u.interpolate(uexact) + + v = TestFunction(V) + F = ( + inner(Dt(u), v) * dx + + inner(grad(u), grad(v)) * dx + - inner(rhs, v) * dx + ) + F_GL = replace(F, {u: u_GL}) + + luparams = {"mat_type": "aij", "ksp_type": "preonly", "pc_type": "lu"} + + stepper = GalerkinTimeStepper( + F, order, t, dt, u, bcs=bcs, + solver_parameters=luparams + ) + stepper_GL = TimeStepper( + F_GL, butcher_tableau, t, dt, u_GL, bcs=bcs, solver_parameters=luparams + ) + + t_end = 1.0 + while float(t) < t_end: + if float(t) + float(dt) > t_end: + dt.assign(t_end - float(t)) + stepper.advance() + stepper_GL.advance() + t.assign(float(t) + float(dt)) + assert (errornorm(u_GL, u) / norm(u)) < 1.e-10 + + +def galerkin_wave(n, deg, alpha, order): + N = 2**n + msh = UnitIntervalMesh(N) + + params = {"snes_type": "ksponly", + "ksp_type": "preonly", + "mat_type": "aij", + "pc_type": "lu"} + + V = FunctionSpace(msh, "CG", deg) + W = FunctionSpace(msh, "DG", deg - 1) + Z = V * W + + x, = SpatialCoordinate(msh) + + MC = MeshConstant(msh) + t = MC.Constant(0.0) + dt = MC.Constant(alpha / N) + + up = project(as_vector([0, sin(pi*x)]), Z) + u, p = split(up) + + v, w = TestFunctions(Z) + + F = (inner(Dt(u), v)*dx + inner(u.dx(0), w) * dx + + inner(Dt(p), w)*dx - inner(p, v.dx(0)) * dx) + + E = 0.5 * (inner(u, u)*dx + inner(p, p)*dx) + + stepper = GalerkinTimeStepper(F, order, t, dt, up, + solver_parameters=params) + + energies = [] + + while (float(t) < 1.0): + if (float(t) + float(dt) > 1.0): + dt.assign(1.0 - float(t)) + stepper.advance() + t.assign(float(t) + float(dt)) + energies.append(assemble(E)) + + return np.array(energies) + + +@pytest.mark.parametrize('N', [2**j for j in range(2, 3)]) +@pytest.mark.parametrize('deg', (2,)) +@pytest.mark.parametrize('order', (1, 2)) +def test_wave_eq_galerkin(deg, N, order): + energy = galerkin_wave(N, deg, 0.3, order) + assert np.allclose(energy[1:], energy[:-1])