JOSS paper draft and documentation updates

README.md

@@ -11,20 +11,26 @@ seismic imaging project.
 ## Quickstart
 
 Seigen requires the installation of Firedrake and must be run from
-within the Firedrake virtual environment. To first install Firedrake
+within the Firedrake virtual environment (venv). To first install Firedrake
 please follow these [instructions](http://www.firedrakeproject.org/download.html#).
 
-Once Firedrake is installed and the virtual environment is activated, you can install
+Once Firedrake is installed and the venv is activated, you can install
 Seigen using the following commands:
 
 ```
 git clone https://github.com/opesci/seigen.git
 pip install -e seigen
 ```
 
+## Documentation
+
+The documentation can be built by navigating to the `docs` directory and running `make html` at the command line. A new `build` folder will then be created containing the documentation in HTML format.
+
+Details regarding the model equations and the discretisation and solution procedure are provided, along with a description of the pulse propagation application (see `tests/pulse`). This acts as an example of how to set up a simulation in Seigen.
+
 ## Licence
 
-Seigen is open-source software and is released under the [MIT License](https://opensource.org/licenses/MIT). See the file ``LICENSE.md`` for more information.
+Seigen is open-source software and is released under the [MIT License](https://opensource.org/licenses/MIT). See the file `LICENSE.md` for more information.
 
 ## Contact
 

codemeta.json

@@ -0,0 +1,17 @@
+{
+  "@context": "https://raw.githubusercontent.com/mbjones/codemeta/master/codemeta.jsonld",
+  "@type": "Code",
+  "author": [
+
+  ],
+  "identifier": "",
+  "codeRepository": "https://github.com/opesci/seigen",
+  "datePublished": "2018-06-06",
+  "dateModified": "2018-06-06",
+  "dateCreated": "2018-06-06",
+  "description": "Seigen is a finite element-based elastic wave equation solver for seimological problems.",
+  "keywords": "seismology, seismic modelling, numerical modelling, code generation, finite element method, elastic wave equation",
+  "license": "MIT",
+  "title": "Seigen",
+  "version": "v0.9-beta"
+}

docs/source/discretisation.rst

@@ -0,0 +1,87 @@
+Discretisation and solution procedure
+=====================================
+
+The model equations are discretised using the Galerkin finite element method, and are advanced forward in time using a fourth-order leapfrog scheme. These are described below.
+
+Spatial discretisation
+----------------------
+
+The discretisation of the governing equations in space is performed using the Galerkin finite element method. As a first step, the weak form of the elastic wave equation is derived by multiplying through by a test function :math:`\mathbf{w} \in H^1(\Omega)^d`, where :math:`H^1(\Omega)^d` is the first Hilbertian Sobolev space of dimension :math:`d`, and integrating over the domain :math:`\Omega` by parts on the RHS. For the velocity equation in the velocity-stress formulation, this yields
+
+.. math:: \int_\Omega{\mathbf{w}\cdot\rho\frac{\partial\mathbf{u}}{\partial t}}\ \mathrm{dV} = -\int_\Omega{\mathbb{T}\cdot\nabla\mathbf{w}}\ \mathrm{dV} + \int_{\partial\Omega_e}{(\mathbf{w}\cdot\mathbb{T}|_{\partial\Omega_e})\cdot\mathbf{n}_e}\ \mathrm{dS} + \int_{\partial\Omega}{(\mathbf{w}\cdot\mathbb{T})\cdot\mathbf{n}}\ \mathrm{ds}
+    :label: weak_form_velocity
+
+where :math:`\mathbf{n}` is the vector normal to the domain boundary, :math:`\partial\Omega`, and :math:`\mathbf{n}_e` is the vector normal to the facet of element :math:`e`, :math:`\partial\Omega_e`. Given this weak form, a solution :math:`\mathbf{u} \in H^1(\Omega)^d` is sought for all test functions :math:`\mathbf{w} \in H^1(\Omega)^d`.
+
+Discrete representations for both :math:`\mathbf{w}` and :math:`\mathbf{u}`, given by a linear combination of basis functions :math:`\{\phi_i\}_{i=1}^{N_e}`, which are discontinuous across element facets, are then used to replace them in :eq:`weak_form_velocity`:
+
+.. math:: \mathbf{w} = \sum_{i=1}^{N_e} \phi_i\mathbf{w}_i,
+    :label: discrete_test_function
+
+.. math:: \mathbf{u} = \sum_{j=1}^{N_e} \phi_j\mathbf{u}_j,
+    :label: discrete_velocity_trial_function
+
+where :math:`N` is the total number of solution nodes. The stress field :math:`\mathbb{T}` is also represened by its own set of basis functions :math:`\{\psi_i\}_{i=1}^{N_e}`
+
+.. math:: \mathbb{T} = \sum_{k=1}^{N_e} \psi_k\mathbb{T}_k.
+    :label: discrete_stress_trial_function
+
+Note that the underlying Firedrake framework permits the use of a range of basis function types and high-order function spaces, and it is assumed in Seigen that both the velocity and stress fields share the same function space.
+
+By substituting :eq:`discrete_test_function`, :eq:`discrete_velocity_trial_function` and :eq:`discrete_stress_trial_function` into :eq:`weak_form_velocity`, and applying the fact that :math:`\mathbf{w}_i` are arbitrary, this yields
+
+.. math:: \sum_{j=1}^{N_e}\int_{\Omega_e}{\phi_i\cdot\rho\phi_j}\ \mathrm{dV}\ \frac{\partial\mathbf{u}_j}{\partial t} = -\sum_{k=1}^{N_e}\int_{\Omega_e}{\psi_k\cdot\nabla\phi_i}\ \mathrm{dV}\ \mathbb{T}_k + \sum_{k=1}^{N_e}\int_{\partial\Omega_e}{(\phi_i\cdot\mathbb{T}|_{\partial\Omega_e})\cdot\mathbf{n}_e}\ \mathrm{dS} + \sum_{k=1}^{N_e}\int_{\partial\Omega}{(\phi_i\cdot\mathbb{T})\cdot\mathbf{n}}\ \mathrm{ds}
+    :label: weak_form_velocity
+
+for all :math:`\phi_i` in element :math:`e`. The volume integrals and interior facet integrals have been restricted to an individual element :math:`e` here because the discontinuous formulation results in each element becoming its own independent problem.
+
+First-order upwinding
+~~~~~~~~~~~~~~~~~~~~~
+
+As a result of the discontinuous nature of the solution fields and test functions at element boundaries, the values of :math:`\mathbb{T}|_{\partial\Omega_e}` and :math:`\mathbf{u}|_{\partial\Omega_e}` in the facet integral terms must be treated carefully. For this work we use a first-order upwinding scheme such that the average values of :math:`\mathbb{T}` and :math:`\mathbf{u}` across the facets are used.
+
+The discrete system
+~~~~~~~~~~~~~~~~~~~
+
+The end result is a discrete system of size :math:`N \times N` for both the velocity and stress equations:
+
+.. math:: \mathbf{M}\frac{\partial\mathbf{u}}{\partial t} = \mathbf{D}\mathbb{T},
+
+where
+
+.. math:: \mathbf{M}_{ij} = \int_{\Omega_e}{\phi_i\rho\phi_j\ \mathrm{dV}},
+
+.. math:: \mathbf{D}_{ij} = -\int_{\Omega_e}{\psi_j\cdot\nabla\phi_i}\ \mathrm{dV} + \int_{\partial\Omega_e}{(\phi_i\cdot\mathbb{T}|_{\partial\Omega_e})\cdot\mathbf{n}_e}\ \mathrm{dS} + \int_{\partial\Omega}{(\phi_i\cdot\mathbb{T})\cdot\mathbf{n}}\ \mathrm{ds}
+
+The coefficients of the discrete representation of :math:`\mathbf{u}` must be solved for using a numerical solution method. A similar procedure can be performed on the stress equation. 
+
+Temporal discretisation
+-----------------------
+
+A fourth-order explicit leap-frog scheme is used to treat the time derivatives in the elastic wave equation. This is based on a truncated Taylor series expansion of the velocity and stress fields whilst staggering their solutions by half a time unit. Hence, the velocity is first solved to obtain a solution at time level :math:`n+1` using information about the stress at time :math:`n+\frac{1}{2}`. This new velocity solution is then in turn used to solve for the stress field at time :math:`n+\frac{3}{2}`. Mathematically, this is written as
+
+.. math:: \mathbf{u}^{n+1} = \mathbf{u}^{n} + \Delta t \frac{\partial\mathbf{u}^{n+\frac{1}{2}}}{\partial t} + \frac{\Delta t^3}{24}\frac{\partial^3\mathbf{u}^{n+\frac{1}{2}}}{\partial t^3}
+
+.. math:: \mathbb{T}^{n+\frac{3}{2}} = \mathbb{T}^{n + \frac{1}{2}} + \Delta t \frac{\partial\mathbb{T}^{n+1}}{\partial t} + \frac{\Delta t^3}{24}\frac{\partial^3\mathbb{T}^{n+1}}{\partial t^3}
+
+where the higher-order terms (:math:`O(\Delta t^5)`) have been neglected. The remaining time derivatives can be evaluated using auxiliary fields which need to be solved for at each time-step:
+
+.. math:: \mathbf{u}^{n+1} = \mathbf{u}^{n} + \Delta t \mathbf{u}^{n+\frac{1}{2}}_{\bigstar} + \frac{\Delta t^3}{24}\mathbf{u}^{n+\frac{1}{2}}_{\bigstar\bigstar}
+
+.. math:: \mathbb{T}^{n+\frac{3}{2}} = \mathbb{T}^{n + \frac{1}{2}} + \Delta t \mathbb{T}^{n+1}_{\bigstar} + \frac{\Delta t^3}{24} \mathbb{T}^{n+1}_{\bigstar\bigstar}
+
+These auxiliary fields are defined as
+
+.. math:: \mathbf{u}^{n+\frac{1}{2}}_{\bigstar} = f(\mathbb{T}^{n + \frac{1}{2}})
+.. math:: \mathbb{T}^{n + \frac{1}{2}}_{\bullet} = g(\mathbf{u}^{n+\frac{1}{2}}_{\bigstar})
+.. math:: \mathbf{u}^{n+\frac{1}{2}}_{\bigstar\bigstar} = f(\mathbb{T}^{n + \frac{1}{2}}_{\bullet})
+
+and
+
+.. math:: \mathbb{T}^{n+1}_{\bigstar} = g(\mathbf{u}^{n+1})
+.. math:: \mathbf{u}^{n+1}_{\bullet} = f(\mathbb{T}^{n+1}_{\bigstar})
+.. math:: \mathbb{T}^{n+1}_{\bigstar\bigstar} = g(\mathbf{u}^{n+1}_{\bullet})
+
+where :math:`f` and :math:`g` are the right-hand sides of the velocity and stress equations (from the velocity-stress formulation of the elastic wave equation), respectively.
+
+Both implicit and explicit solvers are available, encapsulated within the ``ElasticLF4`` class. The former solves the individual finite element variational problems (defined above) implicitly using a linear solver, while the latter explicitly solves the variational problems by assembling RHS vectors and multiplying them with the according global inverse mass matrix.

docs/source/equations.rst

@@ -26,7 +26,7 @@ Absorption term
 An optional sponge/absorption region can be specified through an additional term, :math:`-\sigma\mathbf{u}`, added to the velocity equation.
 
 Source term
-~~~~~~~~~~~~~~~
+~~~~~~~~~~~
 
 An optional source can be specified through an additional term, :math:`\mathbf{S}`, added to the velocity equation.
 

docs/source/example.rst

@@ -1,7 +1,28 @@
-Example
-=======
+Example: Pulse propagation
+==========================
 
-The problem of a one-dimensional travelling pulse wave will be used here as an example of how to set up a simulation in Seigen. This involves a single Python file (``pulse_1d_lf4.py``) which describes the problem to be solved, and can be found in the ``tests/pulse`` directory.
+A one-dimensional travelling pulse wave will be considered here as an example of how to set up a simulation in Seigen. This involves a single Python file (``pulse_1d_lf4.py``) which describes the problem to be solved, and can be found in the ``tests/pulse`` directory.
+
+Background
+----------
+
+This test case considered the propagation of a smooth wave across 0 :math:`\leq x \leq` 4 m. Not only does this allow one to check that the wave is propagating at the expected velocity, it also demonstrates that the wave does not dissipate over time as a result of the non-dissipative nature of the numerical scheme.
+
+The initial condition for the :math:`x`-component of velocity and :math:`xx`-component of stress are given by
+
+.. math::  \mathbf{u}_x = e^{-50(x-1)^2}
+
+.. math:: \mathbb{T}_{xx} = -e^{-50(x-1)^2}
+
+An absorption region was applied at both ends of the domain (0 :math:`\leq x \leq` 0.5 m and 3.5 :math:`\leq x \leq` 4.0 m), with an absorption coefficient of :math:`\sigma = 100`, to model an infinite domain and to prevent any spurious reflections off the end boundaries.
+
+The physical parameters were set to :math:`\mu` = 0.5 Pa, :math:`\lambda` = 0.25 Pa and :math:`\rho` = 1.0 :math:`\mathrm{kgm}^{-3}`. For discretisation purposes, a structured mesh of characteristic element length :math:`\Delta x` = 0.01 m was used. The time-step :math:`\Delta t` of 0.0025 s was set, corresponding to a maximum Courant number of 0.25, and the solutions were advanced forward in time until :math:`T` = 2.0 s. Both the velocity and stress fields were represented by piecewise-discontinuous linear (also known as P1-DG) polynomials, since this simulation did not require such a high degree of accuracy to demonstrate the robustness of the numerical scheme.
+
+The results in the figures below show how the pulse propagates from left to right, whilst displaying excellent agreement with the analytical solution derived by `Delcourte et al. (2009) <https://doi.org/10.1051/proc/2009020>`_. The :math:`x`-component of the velocity field is shown at :math:`t` = 0.25, 1.0 and 2.0 (from left to right). The pulse wave travels at the correct velocity of 1 :math:`\mathrm{ms}^{-1}`. The difference between the analytical and numerical solution, :math:`\epsilon`, is scaled by a factor of 10 to emphasise the insignificant dissipation over time.
+
+.. image::  ux_0.25.png
+.. image::  ux_1.00.png
+.. image::  ux_2.00.png
 
 Modules
 -------
@@ -67,7 +88,7 @@ Additional terms can be added to the solver by setting the ``absorption`` and/or
    elastic.absorption = Expression("x[0] >= 3.5 || x[0] <= 0.5 ? 100.0 : 0")
 
 Initial conditions
-~~~~~~~~~~~~~~~~~~
+------------------
 
 Expressions for velocity and stress at time t = 0 should be provided via the ``u0`` and ``s0`` solver attributes. Again, this uses ``Expresssion`` UFL objects. However, note that these need to be interpolated onto the DG function spaces used by the velocity and stress fields.
 

docs/source/index.rst

@@ -13,7 +13,7 @@ Contents:
 
    introduction
    equations
-   solvers
+   discretisation
    example
 
 Indices and tables