.. SPDX-FileCopyrightText: 2022 SeisSol Group SPDX-License-Identifier: BSD-3-Clause SPDX-LicenseComments: Full text under /LICENSE and /LICENSES/ SPDX-FileContributor: Author lists in /AUTHORS and /CITATION.cff .. _energy_output: Energy output ============= Introduction ------------ The energy output computes the energy of the simulation. It is divided into multiple parts: Energy in the water layer ~~~~~~~~~~~~~~~~~~~~~~~~~ The water layer is modelled as an acoustic medium over the domain :math:`\Omega_a`. The constitutive behaviour is governed by the relation :math:`p = -K \nabla \cdot \mathbf{u}`, with bulk modulus :math:`K`, pressure :math:`p`, and displacement :math:`\mathbf{u}`. **Gravitational potential energy** is associated with the deformed sea surface :math:`\Gamma_\mathrm{free}` (the top boundary of :math:`\Omega_a`): .. math:: W_\mathrm{grav} = \int_{\Gamma_\mathrm{free}} \frac{1}{2} \rho\, g\, \eta^2 \,\mathrm{d}S with :math:`\rho` the density, :math:`g` the gravitational acceleration and :math:`\eta` the sea-surface elevation. It is the work needed to deform the sea surface against gravity, measured from the undisturbed state, and is quadratic in :math:`\eta` under the same linearisation that gives the free-surface boundary condition with gravitational restoring force [LottoDunham2015]_. Reported as ``gravitational_potential_energy``. **Acoustic potential energy** is the energy stored in the compression of the fluid: .. math:: W_\mathrm{ac,pot} = \int_{\Omega_a} \frac{1}{2K}\, p^2 \,\mathrm{d}\mathbf{x} with :math:`p` the acoustic pressure. Reported as ``acoustic_potential_energy``. **Acoustic kinetic energy:** .. math:: W_\mathrm{ac,kin} = \int_{\Omega_a} \frac{1}{2} \rho\, v_i\, v_i \,\mathrm{d}\mathbf{x} with :math:`\rho` the density and :math:`v_i` the velocity. Reported as ``acoustic_kinetic_energy``. **Acoustic energy.** In the acoustics literature the term denotes the sum of the two: .. math:: W_\mathrm{ac} = W_\mathrm{ac,pot} + W_\mathrm{ac,kin} so neither part carries that name on its own. SeisSol writes the two parts to separate columns and prints their sum, with the share of each part next to it. **Mechanical energy of the water layer:** .. math:: W_\mathrm{water} = W_\mathrm{grav} + W_\mathrm{ac} **Dissipation.** The acoustic medium is modelled as inviscid and non-dissipative: .. math:: \dot{D}_\mathrm{ac} = 0 **Energy balance:** .. math:: \frac{\mathrm{d}}{\mathrm{d}t} W_\mathrm{water} = P_\mathrm{ext,ac} - \Phi_{\Gamma} where :math:`P_\mathrm{ext,ac}` is the power of external sources in the acoustic domain and :math:`\Phi_{\Gamma}` is the energy flux through the acoustic-elastic coupling interface :math:`\Gamma`. At this interface, continuity of normal velocity and normal traction ensures that energy leaving the water layer enters the solid domain and vice versa: .. math:: \Phi_{\Gamma} = \int_{\Gamma} p\, v_i\, n_i \,\mathrm{d}S with :math:`n_i` the outward normal of :math:`\Omega_a`. In a closed coupled system without external sources, the total energy :math:`W_\mathrm{water} + W_\mathrm{kin} + W_\mathrm{strain}` (summed over both domains) is conserved; any drift in the discrete balance indicates numerical dissipation or energy injection by the scheme. Energy in the Earth ~~~~~~~~~~~~~~~~~~~ For comprehensive derivations of the energy expressions for elastic, anelastic, anisotropic and poroelastic media used throughout this section, see [Carcione2014]_. The Earth domain is denoted :math:`\Omega_e`. Kinetic energy ^^^^^^^^^^^^^^ For single-phase solid material models (elastic, anisotropic, viscoelastic), the kinetic energy of the solid is: .. math:: W_\mathrm{kin} = \int_{\Omega_e} \frac{1}{2} \rho\, v_i\, v_i \,\mathrm{d}\mathbf{x} with :math:`\rho` the density and :math:`v_i` the velocity. The poroelastic case is two-phase (solid skeleton + pore fluid) and requires additional terms for solid-fluid coupling and relative fluid motion; the corresponding expression is given in the poroelastic subsection below. Strain energy ^^^^^^^^^^^^^ **General form.** The strain energy stored in the deformed solid is: .. math:: W_\mathrm{strain} = \int_{\Omega_e} \frac{1}{2} \epsilon_{ij}\, \sigma_{ij} \,\mathrm{d}\mathbf{x} with :math:`\epsilon_{ij}` the strain tensor and :math:`\sigma_{ij}` the stress tensor. Together with :math:`W_\mathrm{kin}` it makes up the mechanical energy of the solid; the two are reported as ``elastic_strain_energy`` and ``elastic_kinetic_energy``. **Isotropic elastic materials.** For isotropic materials with constitutive law :math:`\sigma_{ij} = \lambda \delta_{ij} \epsilon_{kk} + 2\mu \epsilon_{ij}`, the strain can be expressed purely in terms of stress by inverting the constitutive relation: .. math:: \epsilon_{ij} = \frac{1}{2\mu}\left(\sigma_{ij} - \frac{\lambda}{3\lambda+2\mu}\,\sigma_{kk}\,\delta_{ij}\right) Substituting into the energy expression yields: .. math:: W_\mathrm{strain} = \int_{\Omega_e} \frac{1}{4\mu} \left(\sigma_{ij}\, \sigma_{ij} - \frac{\lambda}{3\lambda+2\mu}\, \sigma_{kk}^2\right) \mathrm{d}\mathbf{x} with :math:`\lambda` and :math:`\mu` the Lamé parameters. Written out explicitly (:math:`\sigma_{kk}^2` is shorthand for :math:`\sigma_{kk}\,\sigma_{ll} = (\mathrm{tr}\,\boldsymbol{\sigma})^2`): .. math:: W_\mathrm{strain} = \int_{\Omega_e} \frac{1}{4\mu}\Big(\sigma_{xx}^2 + \sigma_{yy}^2 + \sigma_{zz}^2 + 2\sigma_{xy}^2 + 2\sigma_{xz}^2 + 2\sigma_{yz}^2 - \frac{\lambda}{3\lambda+2\mu}\left(\sigma_{xx}+\sigma_{yy}+\sigma_{zz}\right)^2\Big)\,\mathrm{d}\mathbf{x} **Anisotropic elastic materials.** For anisotropic materials with constitutive law :math:`\sigma_{ij} = c_{ijkl}\,\epsilon_{kl}`, the strain is obtained by introducing the compliance tensor :math:`s_{ijkl}`, defined as the inverse of the stiffness tensor on the space of symmetric second-order tensors: .. math:: c_{ijmn}\, s_{mnkl} = \tfrac{1}{2}\left(\delta_{ik}\,\delta_{jl} + \delta_{il}\,\delta_{jk}\right) so that :math:`\epsilon_{ij} = s_{ijkl}\,\sigma_{kl}`. The elastic energy becomes: .. math:: W_\mathrm{strain} = \int_{\Omega_e} \frac{1}{2}\, s_{ijkl}\, \sigma_{ij}\, \sigma_{kl} \,\mathrm{d}\mathbf{x} In Voigt notation, with the stress vector :math:`\boldsymbol{\sigma} = (\sigma_{xx}, \sigma_{yy}, \sigma_{zz}, \sigma_{yz}, \sigma_{xz}, \sigma_{xy})^T`, the :math:`6 \times 6` Voigt stiffness matrix :math:`\mathbf{C}` (with components :math:`C_{IJ}` corresponding to :math:`c_{ijkl}` under the standard index pairing) and its inverse :math:`\mathbf{S} = \mathbf{C}^{-1}`, the energy reads: .. math:: W_\mathrm{strain} = \int_{\Omega_e} \frac{1}{2}\, \boldsymbol{\sigma}^T \mathbf{S}\, \boldsymbol{\sigma} \,\mathrm{d}\mathbf{x} This relation uses the engineering convention in which the off-diagonal entries of the Voigt strain vector are :math:`2\epsilon_{ij}` (engineering shear strains); only under this convention does the matrix identity :math:`\mathbf{S} = \mathbf{C}^{-1}` reproduce the tensor relation :math:`\epsilon_{ij} = s_{ijkl}\,\sigma_{kl}`. .. warning:: The Voigt vector above is ordered :math:`(xx, yy, zz, yz, xz, xy)`, which is also the ordering the material parameters :math:`c_{IJ}` follow. SeisSol's *quantity* vector, and hence the stress moments the energy output works with, is ordered :math:`(xx, yy, zz, xy, yz, xz)` -- the last three are permuted. The compliance therefore has to be assembled in the quantity ordering before it is contracted with the stress moments. Getting this wrong is silent for any material with an isotropic shear block, and only shows up for genuinely anisotropic ones. The isotropic result is recovered when :math:`c_{ijkl} = \lambda\,\delta_{ij}\,\delta_{kl} + \mu\,(\delta_{ik}\,\delta_{jl} + \delta_{il}\,\delta_{jk})`, the standard isotropic elasticity tensor in the Lamé parameters :math:`\lambda` and :math:`\mu`. **Dissipation.** Both isotropic and anisotropic elastic materials are non-dissipative: .. math:: \dot{D} = 0 **Energy balance.** The total mechanical energy is exactly conserved: .. math:: \frac{\mathrm{d}}{\mathrm{d}t}\left(W_\mathrm{kin} + W_\mathrm{strain}\right) = P_\mathrm{ext} where :math:`P_\mathrm{ext}` is the power of external sources (boundary tractions, body forces, fault slip). Without external sources: .. math:: W_\mathrm{kin}(t) + W_\mathrm{strain}(t) = W_\mathrm{kin}(0) + W_\mathrm{strain}(0) = \mathrm{const} This conservation property makes the energy balance a useful diagnostic for numerical accuracy: any drift in :math:`W_\mathrm{kin} + W_\mathrm{strain}` indicates numerical dissipation or energy injection by the scheme. Viscoelastic energy ^^^^^^^^^^^^^^^^^^^ SeisSol implements the generalised Maxwell body with *strain-rate-like* anelastic variables, which is worth stating explicitly because much of the literature uses stress-like memory variables instead. The system solved is .. math:: \dot{\sigma}_{ij} = c_{ijkl}^U\,\dot{\epsilon}_{kl} - \sum_{l=1}^{L} d_{ijkl}^{(l)}\,\vartheta_{kl}^{(l)} .. math:: \dot{\vartheta}_{ij}^{(l)} = \omega_l\left(\dot{\epsilon}_{ij} - \vartheta_{ij}^{(l)}\right) with the unrelaxed stiffness :math:`c^U` (this is what ``material.lambda`` and ``material.mu`` hold once ``fitAttenuation`` has run), the relaxation frequencies :math:`\omega_l` (``material.omega``), and the modulus defect :math:`d^{(l)}` of mechanism :math:`l`. The latter is stored in ``material.theta`` as .. math:: \theta^{(l)} = \bigl(-(\delta\lambda_l + 2\,\delta\mu_l),\; -\delta\lambda_l,\; -2\,\delta\mu_l\bigr) so that :math:`d^{(l)}` is the isotropic stiffness built from :math:`(\delta\lambda_l, \delta\mu_l)`, and the source tensor assembled by ``getTransposedSourceCoefficientTensor`` is exactly :math:`-d^{(l)}`. Note that :math:`\delta\lambda_l` and :math:`\delta\mu_l` are fitted **independently** (from two separate least-squares solves in ``fitAttenuation``); in particular :math:`d^{(l)} \neq \omega_l\, c^U`, so formulations that assume a single scalar weight per mechanism do not apply here. The relaxed stiffness is .. math:: c^R = c^U - \sum_{l=1}^{L} d^{(l)}, \qquad \lambda^R = \lambda^U - \sum_l \delta\lambda_l, \qquad \mu^R = \mu^U - \sum_l \delta\mu_l. **The anelastic variables are strain rates.** The relaxation equation forces :math:`[\vartheta] = [\dot{\epsilon}]`. Writing :math:`\xi^{(l)}` for their time integral and :math:`e^{(l)} := \epsilon - \xi^{(l)}` for the spring strain of Maxwell branch :math:`l`, one gets :math:`\dot{\xi}^{(l)} = \omega_l\, e^{(l)}`, hence .. math:: e_{ij}^{(l)} = \frac{\vartheta_{ij}^{(l)}}{\omega_l}. This is what makes the whole energy budget computable without any history: every quantity below is a function of the instantaneous state :math:`(\sigma, \vartheta^{(1..L)})` alone. **Recovering the strain.** The total strain is not a state variable, but follows from the relaxed constitutive relation: .. math:: \sigma^R_{ij} := \sigma_{ij} - \sum_{l=1}^{L} \frac{1}{\omega_l}\, d_{ijkl}^{(l)}\,\vartheta_{kl}^{(l)} \;=\; c^R_{ijkl}\,\epsilon_{kl}. **Stored energy.** The free energy splits into the strain energy of the relaxed spring and the strain energy of the :math:`L` Maxwell springs, :math:`W_\mathrm{stored} = W_\mathrm{strain} + W_\mathrm{ane}`, with .. math:: W_\mathrm{strain} = \frac{1}{2}\int_{\Omega_e} c^R_{ijkl}\,\epsilon_{ij}\,\epsilon_{kl}\,\mathrm{d}\mathbf{x}, \qquad W_\mathrm{ane} = \frac{1}{2}\int_{\Omega_e} \sum_{l=1}^{L} d^{(l)}_{ijkl}\,e^{(l)}_{ij}\,e^{(l)}_{kl}\,\mathrm{d}\mathbf{x} These are reported as ``elastic_strain_energy`` and ``anelastic_strain_energy``. For :math:`L = 0` the first reduces to the purely elastic expression, as it must, and the second vanishes. Both belong to the same reporting group as the kinetic energy, so the terminal line reads .. code-block:: none Elastic mechanical energy: , kinetic x % , strain y % , anelastic z % with the total being :math:`W_\mathrm{kin} + W_\mathrm{stored}`. Reporting only the equilibrium spring would understate the stored energy by the branch contribution, which grows with the attenuation. .. warning:: :math:`\tfrac{1}{2}\,\sigma_{ij}\,s^U_{ijkl}\,\sigma_{kl}` -- the elastic formula evaluated with the unrelaxed compliance -- is **not** the stored energy, and neither is the same expression with the relaxed compliance. The :math:`\sigma^R` correction cannot be absorbed into a choice of moduli; both variants are wrong by several percent, with a sign that depends on the loading history. **Dissipation rate.** The power absorbed by the dashpots: .. math:: \dot{D} = \int_{\Omega_e} \sum_{l=1}^{L} \omega_l\; d^{(l)}_{ijkl}\,e^{(l)}_{ij}\,e^{(l)}_{kl}\,\mathrm{d}\mathbf{x} \;=\; \int_{\Omega_e} \sum_{l=1}^{L} \frac{1}{\omega_l}\; d^{(l)}_{ijkl}\,\vartheta^{(l)}_{ij}\,\vartheta^{(l)}_{kl}\,\mathrm{d}\mathbf{x} \;\geq\; 0 This is reported as ``viscous_dissipation_rate``, in watts. It is an *instantaneous rate*, not a cumulative energy: the time integral .. math:: W_\mathrm{diss}(t) = \int_0^t \dot{D}(t')\,\mathrm{d}t' is left to postprocessing, since accumulating it in the solver would require carrying an extra state through the time-stepping kernels. **Volumetric/deviatoric evaluation.** All of the above are isotropic quadratic forms, so the implementation evaluates them in a volumetric/deviatoric split and never inverts a :math:`6\times 6` matrix. With :math:`\Delta K_l = \delta\lambda_l + \tfrac{2}{3}\delta\mu_l`: .. math:: \mathrm{tr}\,\sigma^R = \mathrm{tr}\,\sigma - \sum_l \frac{3\,\Delta K_l}{\omega_l}\,\mathrm{tr}\,\vartheta^{(l)}, \qquad \mathrm{dev}\,\sigma^R = \mathrm{dev}\,\sigma - \sum_l \frac{2\,\delta\mu_l}{\omega_l}\,\mathrm{dev}\,\vartheta^{(l)} .. math:: d^{(l)}_{ijkl}\,a_{ij}\,a_{kl} = \Delta K_l\,(\mathrm{tr}\,a)^2 + 2\,\delta\mu_l\,\|\mathrm{dev}\,a\|^2 **Admissibility.** The least-squares fit in ``fitAttenuation`` does not constrain the signs of its coefficients, so a poor combination of :math:`Q` and frequency band can yield a fit that reproduces the target :math:`Q` while being thermodynamically inconsistent (a branch with :math:`\delta\mu_l \leq 0`, or non-positive relaxed moduli). The reported dissipation would then be negative. SeisSol checks this once during setup and warns; the remedy is a different central frequency or frequency ratio. **Energy balance:** .. math:: \frac{\mathrm{d}}{\mathrm{d}t}\left(W_\mathrm{kin} + W_\mathrm{stored}\right) = P_\mathrm{ext} - \dot{D} Without external sources, the total energy decreases monotonically due to internal relaxation. **Viscoacoustic materials.** Without a shear modulus, only the bulk modulus relaxes. The viscoacoustic material therefore carries one scalar memory variable per mechanism, driven by the divergence of the velocity -- i.e. the trace of its tensor counterpart above -- and everything above reduces to its volumetric part. The same holds for a cell with :math:`\mu = 0` in a viscoelastic build. The Maxwell branch springs then hold potential energy of the sound field, :math:`W_\mathrm{ac,ane}`, which is reported as ``anelastic_potential_energy`` in the acoustic group: .. code-block:: none Acoustic mechanical energy: , kinetic x % , potential y % , anelastic z % Poroelastic energy ^^^^^^^^^^^^^^^^^^ For poroelastic materials following Biot's theory [Biot1962]_, the constitutive equations are: .. math:: \sigma_{ij} = c_{ijkl}^d\,\epsilon_{kl} - \alpha\, p\, \delta_{ij} .. math:: p = M\left(\zeta - \alpha\, \epsilon_{kk}\right) with drained stiffness :math:`c_{ijkl}^d`, Biot coefficient :math:`\alpha`, Biot modulus :math:`M`, pore pressure :math:`p`, and fluid content increment :math:`\zeta`. **Strain energy.** The total stored energy has a skeleton and a fluid contribution: .. math:: W_p = \frac{1}{2}\int_{\Omega_e}\left(\sigma_{ij}\,\epsilon_{ij} + p\,\zeta\right)\mathrm{d}\mathbf{x} In quadratic form over the state variables :math:`(\epsilon_{ij},\,\zeta)`, substituting the constitutive relations gives the compact form: .. math:: W_p = \frac{1}{2}\int_{\Omega_e}\left(c_{ijkl}^d\,\epsilon_{ij}\,\epsilon_{kl} + M\,\bigl(\zeta - \alpha\,\epsilon_{kk}\bigr)^2\right)\mathrm{d}\mathbf{x} Expressed purely in stresses and pore pressure, using the Biot effective stress :math:`\sigma'_{ij} = \sigma_{ij} + \alpha\, p\, \delta_{ij}` and the drained compliance :math:`s_{ijkl}^d` (defined analogously to :math:`s_{ijkl}` above, with :math:`c_{ijkl}^d` in place of :math:`c_{ijkl}`): .. math:: W_p = \frac{1}{2}\int_{\Omega_e}\left(s_{ijkl}^d\,\sigma'_{ij}\,\sigma'_{kl} + \frac{p^2}{M}\right)\mathrm{d}\mathbf{x} For isotropic drained material (Lamé parameters :math:`\lambda_d, \mu`): .. math:: s_{ijkl}^d\,\sigma'_{ij}\,\sigma'_{kl} = \frac{1}{2\mu}\left(\sigma'_{ij}\,\sigma'_{ij} - \frac{\lambda_d}{3\lambda_d+2\mu}\,{\sigma'_{kk}}^2\right) **Kinetic energy.** With solid velocity :math:`v_i^s`, Darcy velocity :math:`w_i = \phi(\dot{u}_i^f - \dot{u}_i^s)`, total density :math:`\rho = (1-\phi)\rho_s + \phi\rho_f`, and tortuosity tensor :math:`T_{ij}`: .. math:: W_\mathrm{kin} = \frac{1}{2}\int_{\Omega_e}\left(\rho\, v_i^s\, v_i^s + 2\rho_f\, v_i^s\, w_i + \frac{\rho_f}{\phi}\, T_{ij}\, w_i\, w_j\right)\mathrm{d}\mathbf{x} The three terms represent the kinetic energy of the skeleton, the solid-fluid coupling, and the fluid motion relative to the skeleton, respectively. Both differ from their single-phase counterparts and are reported under their own names, ``poroelastic_strain_energy`` and ``poroelastic_kinetic_energy``. **Dissipation rate.** Viscous dissipation due to Darcy friction: .. math:: \dot{D} = \int_{\Omega_e} \frac{\eta_f}{\kappa}\, w_i\, w_i \,\mathrm{d}\mathbf{x} with dynamic fluid viscosity :math:`\eta_f` and permeability :math:`\kappa`. For anisotropic permeability, replace :math:`\frac{\eta_f}{\kappa}\, w_i\, w_i` by :math:`\eta_f\, \kappa_{ij}^{-1}\, w_i\, w_j`. **Energy balance:** .. math:: \frac{\mathrm{d}}{\mathrm{d}t}\left(W_\mathrm{kin} + W_p\right) = P_\mathrm{ext} - \dot{D} Without external sources, the total energy decreases monotonically due to viscous fluid friction. .. _darcy_dissipation_stiff: .. warning:: For typical permeabilities, the Darcy friction is stiff: the relative fluid velocity :math:`w_i` relaxes within :math:`\rho_f T \kappa / (\phi \eta_f)`, usually microseconds and far below the time step. The solver resolves that relaxation in its element-local predictor, but the state at an output time still carries the part of :math:`w_i` that the last flux update added, which the next predictor removes again within the relaxation time. ``darcy_dissipation_rate`` evaluated from that state is therefore not the dissipation the scheme actually performs: it grows with the length of the last, shortened time step before the output, and can exceed the actual dissipation by orders of magnitude. In this regime -- seismic frequencies far below the Biot characteristic frequency, where the physical Darcy dissipation of the waves is small -- do not integrate the column into an energy budget. It is meaningful when the relaxation time is long compared to the time step. Earthquake source energy ~~~~~~~~~~~~~~~~~~~~~~~~~ **Total frictional work** done by the stress change (cf. eq. 3 in [MaArchuleta2006]_): .. math:: W_\mathrm{total} = -\int_{0}^{t_f} \int_{\Sigma} \Delta\boldsymbol{\sigma}(t) \cdot \Delta\dot{\mathbf{u}}(t) \,\mathrm{d}S\,\mathrm{d}t with :math:`\Sigma` the fault surface, :math:`\Delta\boldsymbol{\sigma}(t) = \boldsymbol{\sigma}(t) - \boldsymbol{\sigma}(0)` the shear traction change, :math:`\Delta\dot{\mathbf{u}}(t)` the fault slip rate, and :math:`t_f` the end time of the simulation. **Static frictional work** done by the stress change (cf. eq. 4 in [MaArchuleta2006]_): .. math:: W_\mathrm{static} = -\int_{\Sigma} \frac{1}{2} \Delta\boldsymbol{\sigma}(t_f) \cdot \Delta\mathbf{u}(t_f) \,\mathrm{d}S **Radiated energy** can then be computed as (cf. eq. 5 in [MaArchuleta2006]_): .. math:: E_\mathrm{r} = W_\mathrm{total} - W_\mathrm{static} Other quantities ~~~~~~~~~~~~~~~~ **Potency:** .. math:: \int_{\Sigma} \Delta u_\mathrm{acc}(t_f) \,\mathrm{d}S with :math:`\Delta u_\mathrm{acc}` the accumulated fault slip (scalar). **Seismic moment:** .. math:: \int_{\Sigma} \mu\, \Delta u_\mathrm{acc}(t_f) \,\mathrm{d}S with :math:`\mu` the shear modulus (second Lamé parameter). **Plastic moment:** .. math:: \int_{\Omega_e} \mu\, \eta_p \,\mathrm{d}\mathbf{x} with :math:`\mu` the shear modulus and :math:`\eta_p` a scalar measure of accumulated plastic strain (off-fault material damage). **Total momentum:** .. math:: \int_{\Omega} \rho\, v_i \,\mathrm{d}\mathbf{x} for every material, and :math:`\int_{\Omega} (\rho\, v_i^s + \rho_f\, w_i) \,\mathrm{d}\mathbf{x}` (the momentum of the mixture) for a poroelastic one. It is reported as ``momentumX``, ``momentumY`` and ``momentumZ``. Without sources, and as long as no wave has left the domain, it is conserved. Names in the output ~~~~~~~~~~~~~~~~~~~ Each of the quantities above is written to its own column of the CSV file. Which columns exist depends on the equation system that SeisSol was built with, and, for the acoustic/elastic split, on the materials present in the mesh: a cell with :math:`\mu = 0` contributes to the acoustic columns, every other cell to the solid ones. +------------------------------------+---------------------------+------------------------------------------------+ | CSV column | Symbol | Quantity | +====================================+===========================+================================================+ | ``gravitational_potential_energy`` | :math:`W_\mathrm{grav}` | Potential energy of the deformed sea surface | +------------------------------------+---------------------------+------------------------------------------------+ | ``acoustic_potential_energy`` | :math:`W_\mathrm{ac,pot}` | Energy stored in the compression of the fluid | | | | cells | +------------------------------------+---------------------------+------------------------------------------------+ | ``acoustic_kinetic_energy`` | :math:`W_\mathrm{ac,kin}` | Kinetic energy of the fluid cells | +------------------------------------+---------------------------+------------------------------------------------+ | ``elastic_strain_energy`` | :math:`W_\mathrm{strain}` | Strain energy of the solid cells; for a | | | | viscoelastic material, that of the equilibrium | | | | spring | +------------------------------------+---------------------------+------------------------------------------------+ | ``elastic_kinetic_energy`` | :math:`W_\mathrm{kin}` | Kinetic energy of the solid cells | +------------------------------------+---------------------------+------------------------------------------------+ | ``anelastic_strain_energy`` | :math:`W_\mathrm{ane}` | Strain energy held in the Maxwell branches | | | | (viscoelastic only) | +------------------------------------+---------------------------+------------------------------------------------+ | ``anelastic_potential_energy`` | :math:`W_\mathrm{ac,ane}` | Potential energy held in the Maxwell branches | | | | of the fluid cells (viscoelastic and | | | | viscoacoustic only) | +------------------------------------+---------------------------+------------------------------------------------+ | ``poroelastic_strain_energy`` | :math:`W_p` | Strain energy of skeleton and pore fluid | | | | (poroelastic only) | +------------------------------------+---------------------------+------------------------------------------------+ | ``poroelastic_kinetic_energy`` | :math:`W_\mathrm{kin}` | Kinetic energy of skeleton and pore fluid, | | | | including their coupling (poroelastic only) | +------------------------------------+---------------------------+------------------------------------------------+ | ``viscous_dissipation_rate`` | :math:`\dot{D}` | Power absorbed by the Maxwell dashpots | +------------------------------------+---------------------------+------------------------------------------------+ | ``darcy_dissipation_rate`` | :math:`\dot{D}` | Power absorbed by Darcy friction in the pore | | | | fluid | +------------------------------------+---------------------------+------------------------------------------------+ Every column holds one contribution only. The terminal output, in contrast, prints one line per group -- ``Acoustic mechanical energy`` for :math:`W_\mathrm{ac,pot} + W_\mathrm{ac,kin} + W_\mathrm{ac,ane}`, ``Elastic mechanical energy`` for :math:`W_\mathrm{kin} + W_\mathrm{strain} + W_\mathrm{ane}`, and ``Poroelastic mechanical energy`` for :math:`W_\mathrm{kin} + W_p` -- followed by the share each contribution has in that sum. So a name that stands for a sum never also stands for one of its parts. Summary of energy balances ~~~~~~~~~~~~~~~~~~~~~~~~~~ +---------------------+-------------------------------------+---------------------------------------------+ | Material model | Total energy | Dissipation mechanism | +=====================+=====================================+=============================================+ | Acoustic | :math:`W_\mathrm{grav} + | None internally; energy exchanged via | | | W_\mathrm{ac,pot} + | interface flux :math:`\Phi_\Gamma` | | | W_\mathrm{ac,kin}` | | +---------------------+-------------------------------------+---------------------------------------------+ | Elastic (iso/aniso) | :math:`W_\mathrm{kin} + | None (energy conserved) | | | W_\mathrm{strain}` | | +---------------------+-------------------------------------+---------------------------------------------+ | Viscoelastic | :math:`W_\mathrm{kin} + | Internal relaxation of the Maxwell | | | W_\mathrm{strain} + W_\mathrm{ane}` | branches, ``viscous_dissipation_rate`` | +---------------------+-------------------------------------+---------------------------------------------+ | Viscoacoustic | :math:`W_\mathrm{grav} + | Internal relaxation of the Maxwell | | | W_\mathrm{ac,pot} + | branches, ``viscous_dissipation_rate`` | | | W_\mathrm{ac,kin} + | | | | W_\mathrm{ac,ane}` | | +---------------------+-------------------------------------+---------------------------------------------+ | Poroelastic | :math:`W_\mathrm{kin} + W_p` | Viscous Darcy friction in the pore fluid, | | | | ``darcy_dissipation_rate`` | +---------------------+-------------------------------------+---------------------------------------------+ In all dissipative cases the energy balance reads: .. math:: \frac{\mathrm{d}}{\mathrm{d}t}\left(W_\mathrm{kin} + W_\mathrm{pot}\right) = P_\mathrm{ext} - \dot{D} where :math:`W_\mathrm{pot}` is the respective potential/stored energy and :math:`\dot{D} \geq 0` the dissipation rate. .. note:: The dissipation columns are reported as instantaneous **rates** (unit W), not as cumulative energies. Closing the balance therefore means integrating the reported rate over time in postprocessing and comparing against the change in :math:`W_\mathrm{kin} + W_\mathrm{pot}`. The accuracy of that integration is limited by the output interval, so a coarse ``EnergyOutputInterval`` will not close the budget to machine precision even for a perfectly conservative scheme. For a stiff dissipation, the sampled rate is not representative at all; see :ref:`the warning on the Darcy dissipation `. For coupled acoustic-elastic simulations, the global energy balance over both domains reads: .. math:: \frac{\mathrm{d}}{\mathrm{d}t}\left(W_\mathrm{water} + W_\mathrm{kin} + W_\mathrm{pot}\right) = P_\mathrm{ext} - \dot{D} The interface flux :math:`\Phi_\Gamma` cancels in the sum since the energy leaving one domain enters the other. Configuration ------------- .. code-block:: Fortran &Output OutputFile = 'output/conv' EnergyOutput = 1 EnergyTerminalOutput = 1 EnergyTerminalPrecision = 6 EnergyOutputInterval = 0.05 ComputeVolumeEnergiesEveryOutput = 4 ! Compute volume energies only once every ComputeVolumeEnergiesEveryOutput * EnergyOutputInterval / Energy output ~~~~~~~~~~~~~ Controlled via ``EnergyOutput``. | 0 : no output | 1 : csv output For the example configuration, the output is written in the file "output/conv-energy.csv". Terminal output ~~~~~~~~~~~~~~~ Controlled via ``EnergyTerminalOutput``. | 0 : no output | 1 : additional output to stdout Additionally, the energy can be written to stdout. This can be useful for debugging. To increase the precision of the terminal output, adjust the ``EnergyTerminalPrecision`` field to the desired accuracy (a value of about 15 is sufficient to check for machine accuracy). Note that ``EnergyOutput`` needs to be enabled for the terminal output to work. Output interval ~~~~~~~~~~~~~~~ The output interval is controlled by ``EnergyOutputInterval``. If the output interval is not specified, the energy will be computed at the start of the simulation and at the end of the simulation. Postprocessing and plotting ~~~~~~~~~~~~~~~~~~~~~~~~~~~ The code below suggests a way to process and plot variables of the energy output file: .. code-block:: Python import pandas as pd import numpy as np import matplotlib.pylab as plt df = pd.read_csv("prefix-energy.csv") df = df.pivot_table(index="time", columns="variable", values="measurement") df["seismic_moment_rate"] = np.gradient(df["seismic_moment"], df.index[1]) df.plot(y="seismic_moment_rate", use_index=True) # if ComputeVolumeEnergiesEveryOutput > 1 volume_output = df.dropna() volume_output.plot(y="elastic_strain_energy", use_index=True) plt.show() References ---------- .. [Biot1962] Biot, M. A. (1962). Mechanics of deformation and acoustic propagation in porous media. *Journal of Applied Physics*, 33(4), 1482–1498. doi: `10.1063/1.1728759 `_ .. [Carcione2014] Carcione, J. M. (2014). *Wave Fields in Real Media: Theory and Numerical Simulation of Wave Propagation in Anisotropic, Anelastic, Porous and Electromagnetic Media* (3rd ed.). Elsevier. ISBN 978-0-08-099999-9. .. [LottoDunham2015] Lotto, G. C. and Dunham, E. M. (2015). High-order finite difference modeling of tsunami generation in a compressible ocean from offshore earthquakes. *Computational Geosciences*, 19(2), 327–340. doi: `10.1007/s10596-015-9472-0 `_ .. [MaArchuleta2006] Ma, S. and Archuleta, R. J. (2006). Radiated seismic energy based on dynamic rupture models of faulting. *Journal of Geophysical Research: Solid Earth*, 111, B05315. doi: `10.1029/2005JB004055 `_