Energy output

Introduction

The energy output computes the energy of the simulation. It is divided into multiple parts:

Energy in the water layer (= acoustic medium, \(\Omega_a\))
- Gravitational energy
  \(\int_{\Omega_a} \frac{1}{2} \rho g \eta^2 \,\mathbf{dx}\), with \(\rho\) the density, \(g\) the gravitational acceleration and \(\eta\) the sea-surface elevation.
- Acoustic energy
  \(\int_{\Omega_a} \frac{1}{2K} p^2 \,\mathbf{dx}\), with \(p\) the acoustic pressure and \(K\) the compressibility.
- Acoustic kinetic energy
  \(\int_{\Omega_a} \frac{1}{2} \rho v^2 \,\mathbf{dx}\), with \(\rho\) the density and \(v\) the velocity.
Energy in Earth \(\Omega_e\)
- Elastic kinetic energy
  \(\int_{\Omega_e} \frac{1}{2} \rho v^2 \,\mathbf{dx}\), with \(\rho\) the density and \(v\) the velocity.
- Elastic energy
  \(\int_{\Omega_e} \frac{1}{2} \epsilon_{ij} \sigma_{kl} \,\mathbf{dx}\), with \(\epsilon_{ij}\) the strain tensor and \(\sigma_{ij}\) the stress tensor. It reduces for isotropic materials to \(\int_{\Omega_e} \frac{1}{2\mu} (\sigma_{ij} \sigma_{ij} -\frac{\lambda}{3\lambda+2\mu} \sigma_{kk}^2)\,\mathbf{dx}\), with \(\lambda\) and \(\mu\) the Lame coefficients.
Earthquake source energy
- Total frictional work done by the stress change
  \(W_\mathrm{total} = -\int_{0}^{t_f} \int_{\Sigma} \Delta\mathbf{\sigma}(t) \cdot \Delta\mathbf{\dot{u}}(t) \,\mathbf{dx}dt\), with \(\Sigma\) the fault surface, \(\Delta\mathbf{\sigma}(t) = \mathbf{\sigma}(t) - \mathbf{\sigma}(0)\) the shear traction change, and \(\Delta\mathbf{\dot{u}}(t)\) the fault slip rate, and \(t_f\) the end time of the simulation (see eq. 3 in Ma and Archuleta, 2006).
- Static frictional work done by the stress change
  \(W_\mathrm{static} = -\int_{\Sigma} \frac{1}{2} \mathbf{\Delta\sigma}(t_f) \cdot \mathbf{\Delta u}(t_f) \,\mathbf{dx}\) (see eq. 4 in Ma and Archuleta, 2006).
- Radiated energy can then be computed with:
  \(E_\mathrm{r} = W_\mathrm{total} - W_\mathrm{static}\) (see eq. 5 in Ma and Archuleta, 2006).
Potency
\(\int_{\Sigma} \Delta u_\mathrm{acc}(t_f) \,\mathbf{dx}\), with \(\Delta u_\mathrm{acc}\) the accumulated fault slip (scalar).
Seismic moment
\(\int_{\Sigma} \mu \Delta u_\mathrm{acc}(t_f) \,\mathbf{dx}\), with \(\mu\) the second Lame coefficient.
Plastic moment
\(\int_{\Omega_e} \mu \eta \,\mathbf{dx}\), with \(\mu\) the second Lame coefficient and eta a scalar quantity measuring the accumulated material damage.

Currently, the output is only supported for the elastic wave equation.

Configuration

&Output
OutputFile = 'output/conv'
EnergyOutput = 1
EnergyTerminalOutput = 1
EnergyOutputInterval = 0.05
ComputeVolumeEnergiesEveryOutput = 4 ! Compute volume energies only once every ComputeVolumeEnergiesEveryOutput * EnergyOutputInterval
/

Energy output

0 : no output

1 : csv output

For the example configuration, the output is written in the file “output/conv-energy.csv”.

Terminal output

0 : no output

1 : additional output to stdout

Additionally, the energy can be written to stdout. This can be useful for debugging.

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:

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_energy", use_index=True)

plt.show()