Merge pull request #849 from SeisSol/thomas/fix_fault_energies

refactor code calculating fault energies and fix associated documentation

Documentation/energy-output.rst

@@ -21,12 +21,15 @@ The energy output computes the energy of the simulation. It is divided into mult
     - Elastic energy 
         :math:`\int_{\Omega_e} \frac{1}{2} \epsilon_{ij} \sigma_{kl} \,\mathbf{dx}`, with  :math:`\epsilon_{ij}` the strain tensor and :math:`\sigma_{ij}` the stress tensor. It reduces for isotropic materials to :math:`\int_{\Omega_e} \frac{1}{2\mu} (\sigma_{ij} \sigma_{ij} -\frac{\lambda}{3\lambda+2\mu} \sigma_{kk}^2)\,\mathbf{dx}`, with :math:`\lambda` and :math:`\mu` the Lame coefficients.
 - Earthquake source energy
-    - Total frictional work 
-        :math:`\int_{0}^{t_f} \int_{\Sigma} \frac{1}{2} \Delta\mathbf{\sigma}(t) \cdot \Delta\mathbf{\dot{u}}(t) \,\mathbf{dx}dt` with :math:`\Sigma` the fault surface, :math:`\Delta\mathbf{\sigma}(t)` the shear traction change, and :math:`\Delta\mathbf{\dot{u}}(t)` the fault slip rate, and :math:`t_f` the end time of the simulation.
-    - Static frictional work 
-        :math:`\int_{\Sigma} \frac{1}{2} \mathbf{\Delta\sigma}(t_f) \cdot \mathbf{\Delta u}(t_f) \,\mathbf{dx}`.
+    - Total frictional work done by the stress change
+        :math:`W_\mathrm{total} = -\int_{0}^{t_f} \int_{\Sigma} \Delta\mathbf{\sigma}(t) \cdot \Delta\mathbf{\dot{u}}(t) \,\mathbf{dx}dt`, with :math:`\Sigma` the fault surface, :math:`\Delta\mathbf{\sigma}(t) = \mathbf{\sigma}(t) - \mathbf{\sigma}(0)` the shear traction change, and :math:`\Delta\mathbf{\dot{u}}(t)` the fault slip rate, and :math:`t_f` the end time of the simulation (see eq. 3 in Ma and Archuleta, 2006). 
+    - Static frictional work done by the stress change
+        :math:`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:
+        :math:`E_\mathrm{r} = W_\mathrm{total} - W_\mathrm{static}` (see eq. 5 in Ma and Archuleta, 2006).
+
 - Seismic moment
-        :math:`\int_{\Sigma} \frac{1}{2} \mu \Delta u_\mathrm{acc}(t_f) \,\mathbf{dx}`, with :math:`\mu` the second Lame coefficient and :math:`\Delta u_\mathrm{acc}` the accumulated fault slip (scalar).
+        :math:`\int_{\Sigma} \mu \Delta u_\mathrm{acc}(t_f) \,\mathbf{dx}`, with :math:`\mu` the second Lame coefficient and :math:`\Delta u_\mathrm{acc}` the accumulated fault slip (scalar).
 - Plastic moment 
     :math:`\int_{\Omega_e} \mu \eta  \,\mathbf{dx}`, with :math:`\mu` the second Lame coefficient and \eta a scalar quantity measuring the accumulated material damage.
 

generated_code/DynamicRupture.py

@@ -120,17 +120,17 @@ def interpolateQGenerator(i,h):
   # where the normal points from the plus side to the minus side
   QInterpolatedPlus = OptionalDimTensor('QInterpolatedPlus', aderdg.Q.optName(), aderdg.Q.optSize(), aderdg.Q.optPos(), gShape, alignStride=True)
   QInterpolatedMinus = OptionalDimTensor('QInterpolatedMinus', aderdg.Q.optName(), aderdg.Q.optSize(), aderdg.Q.optPos(), gShape, alignStride=True)
-  slipRateInterpolated = Tensor('slipRateInterpolated', (numberOfPoints,3), alignStride=True)
+  slipInterpolated = Tensor('slipInterpolated', (numberOfPoints,3), alignStride=True)
   tractionInterpolated = Tensor('tractionInterpolated', (numberOfPoints,3), alignStride=True)
-  frictionalEnergy = Tensor('frictionalEnergy', ())
-  timeWeight = Scalar('timeWeight')
+  staticFrictionalWork = Tensor('staticFrictionalWork', ())
+  minusSurfaceArea = Scalar('minusSurfaceArea')
   spaceWeights = Tensor('spaceWeights', (numberOfPoints,), alignStride=True)
 
   computeTractionInterpolated = tractionInterpolated['kp'] <= QInterpolatedMinus['kq'] * aderdg.tractionMinusMatrix['qp'] + QInterpolatedPlus['kq'] * aderdg.tractionPlusMatrix['qp']
   generator.add('computeTractionInterpolated', computeTractionInterpolated)
 
-  accumulateFrictionalEnergy = frictionalEnergy[''] <= frictionalEnergy[''] + timeWeight * tractionInterpolated['kp'] * slipRateInterpolated['kp'] * spaceWeights['k']
-  generator.add('accumulateFrictionalEnergy', accumulateFrictionalEnergy)
+  accumulateStaticFrictionalWork = staticFrictionalWork[''] <= staticFrictionalWork[''] + minusSurfaceArea * tractionInterpolated['kp'] * slipInterpolated['kp'] * spaceWeights['k']
+  generator.add('accumulateStaticFrictionalWork', accumulateStaticFrictionalWork)
 
   ## Dynamic Rupture Precompute
   qPlus = OptionalDimTensor('Qplus', aderdg.Q.optName(), aderdg.Q.optSize(), aderdg.Q.optPos(), gShape, alignStride=True)

src/DynamicRupture/FrictionLaws/FrictionSolverCommon.h

@@ -470,10 +470,7 @@ inline void computeFrictionEnergy(
   auto* qIPlus = reinterpret_cast<QInterpolatedShapeT>(qInterpolatedPlus);
   auto* qIMinus = reinterpret_cast<QInterpolatedShapeT>(qInterpolatedMinus);
 
-  const auto aPlus = impAndEta.etaP * impAndEta.invZp;
   const auto bPlus = impAndEta.etaS * impAndEta.invZs;
-
-  const auto aMinus = impAndEta.etaP * impAndEta.invZpNeig;
   const auto bMinus = impAndEta.etaS * impAndEta.invZsNeig;
 
   using Range = typename NumPoints<Type>::Range;
@@ -501,14 +498,13 @@ inline void computeFrictionEnergy(
       slip[1][i] += timeWeight * interpolatedSlipRate2;
       slip[2][i] += timeWeight * interpolatedSlipRate3;
 
-      const real interpolatedTraction11 = aPlus * qIMinus[o][XX][i] + aMinus * qIPlus[o][XX][i];
-      const real interpolatedTraction12 = bPlus * qIMinus[o][XY][i] + bMinus * qIPlus[o][XY][i];
-      const real interpolatedTraction13 = bPlus * qIMinus[o][XZ][i] + bMinus * qIPlus[o][XZ][i];
+      const real interpolatedTraction12 = bPlus * qIMinus[o][T1][i] + bMinus * qIPlus[o][T1][i];
+      const real interpolatedTraction13 = bPlus * qIMinus[o][T2][i] + bMinus * qIPlus[o][T2][i];
 
       const auto spaceWeight = spaceWeights[i];
-      const auto weight = -1.0 * timeWeight * spaceWeight * doubledSurfaceArea;
-      frictionalEnergy[i] += weight * (interpolatedTraction11 * interpolatedSlipRate1 +
-                                       interpolatedTraction12 * interpolatedSlipRate2 +
+
+      const auto weight = -timeWeight * spaceWeight * doubledSurfaceArea;
+      frictionalEnergy[i] += weight * (interpolatedTraction12 * interpolatedSlipRate2 +
                                        interpolatedTraction13 * interpolatedSlipRate3);
     }
   }

src/Initializer/typedefs.hpp

@@ -407,11 +407,20 @@ struct DRGodunovData {
   real TinvT[seissol::tensor::TinvT::size()];
   real tractionPlusMatrix[seissol::tensor::tractionPlusMatrix::size()];
   real tractionMinusMatrix[seissol::tensor::tractionMinusMatrix::size()];
+  // When integrating quantities over the fault (e.g. mu*slip for the seismic moment)
+  // we need to integrate over each physical element.
+  // The integration is effectively done in the reference element, and the scaling factor of
+  // the transformation, the surface Jacobian (e.g. |n^e(\chi)| in eq. (35) of Uphoff et al. (2023))
+  // is incorporated. This explains the factor 2 (doubledSurfaceArea)
+  //
+  // Uphoff, C., May, D. A., & Gabriel, A. A. (2023). A discontinuous Galerkin method for
+  // sequences of earthquakes and aseismic slip on multiple faults using unstructured curvilinear
+  // grids. Geophysical Journal International, 233(1), 586-626.
   double doubledSurfaceArea;
 };
 
 struct DREnergyOutput {
-  real slip[seissol::tensor::slipRateInterpolated::size()];
+  real slip[seissol::tensor::slipInterpolated::size()];
   real accumulatedSlip[seissol::dr::misc::numPaddedPoints];
   real frictionalEnergy[seissol::dr::misc::numPaddedPoints];
 };

src/ResultWriter/EnergyOutput.cpp

@@ -89,12 +89,11 @@ void EnergyOutput::simulationStart() {
   syncPoint(0.0);
 }
 
-real EnergyOutput::computeStaticWork(
-    const real* degreesOfFreedomPlus,
-    const real* degreesOfFreedomMinus,
-    const DRFaceInformation& faceInfo,
-    const DRGodunovData& godunovData,
-    const real slip[seissol::tensor::slipRateInterpolated::size()]) {
+real EnergyOutput::computeStaticWork(const real* degreesOfFreedomPlus,
+                                     const real* degreesOfFreedomMinus,
+                                     const DRFaceInformation& faceInfo,
+                                     const DRGodunovData& godunovData,
+                                     const real slip[seissol::tensor::slipInterpolated::size()]) {
   real points[NUMBER_OF_SPACE_QUADRATURE_POINTS][2];
   real spaceWeights[NUMBER_OF_SPACE_QUADRATURE_POINTS];
   seissol::quadrature::TriangleQuadrature(points, spaceWeights, CONVERGENCE_ORDER + 1);
@@ -127,12 +126,12 @@ real EnergyOutput::computeStaticWork(
   trKrnl.execute();
 
   real staticFrictionalWork = 0.0;
-  dynamicRupture::kernel::accumulateFrictionalEnergy feKrnl;
-  feKrnl.slipRateInterpolated = slip;
+  dynamicRupture::kernel::accumulateStaticFrictionalWork feKrnl;
+  feKrnl.slipInterpolated = slip;
   feKrnl.tractionInterpolated = tractionInterpolated;
   feKrnl.spaceWeights = spaceWeights;
-  feKrnl.frictionalEnergy = &staticFrictionalWork;
-  feKrnl.timeWeight = -0.5 * godunovData.doubledSurfaceArea;
+  feKrnl.staticFrictionalWork = &staticFrictionalWork;
+  feKrnl.minusSurfaceArea = -0.5 * godunovData.doubledSurfaceArea;
   feKrnl.execute();
 
   return staticFrictionalWork;
@@ -154,7 +153,16 @@ void EnergyOutput::computeDynamicRuptureEnergies() {
     seissol::model::IsotropicWaveSpeeds* waveSpeedsMinus = it->var(dynRup->waveSpeedsMinus);
 
 #if defined(_OPENMP) && !defined(__NVCOMPILER)
-#pragma omp parallel for reduction(+ : totalFrictionalWork, staticFrictionalWork, seismicMoment) default(none) shared(it, drEnergyOutput, faceInformation, timeDerivativeMinus, timeDerivativePlus, godunovData, waveSpeedsPlus, waveSpeedsMinus)
+#pragma omp parallel for reduction(                                                                \
+        + : totalFrictionalWork, staticFrictionalWork, seismicMoment) default(none)                \
+    shared(it,                                                                                     \
+               drEnergyOutput,                                                                     \
+               faceInformation,                                                                    \
+               timeDerivativeMinus,                                                                \
+               timeDerivativePlus,                                                                 \
+               godunovData,                                                                        \
+               waveSpeedsPlus,                                                                     \
+               waveSpeedsMinus)
 #endif
     for (unsigned i = 0; i < it->getNumberOfCells(); ++i) {
       if (faceInformation[i].plusSideOnThisRank) {
@@ -171,13 +179,13 @@ void EnergyOutput::computeDynamicRuptureEnergies() {
                       waveSpeedsPlus[i].sWaveVelocity;
         real muMinus = waveSpeedsMinus[i].density * waveSpeedsMinus[i].sWaveVelocity *
                        waveSpeedsMinus[i].sWaveVelocity;
-        real mu = 2.0 * muPlus * muMinus / (muPlus + muMinus);
+        real mu = muPlus * muMinus / (muPlus + muMinus);
         real seismicMomentIncrease = 0.0;
         for (unsigned k = 0; k < seissol::dr::misc::numberOfBoundaryGaussPoints; ++k) {
           seismicMomentIncrease += drEnergyOutput[i].accumulatedSlip[k];
         }
-        seismicMomentIncrease *= 0.5 * godunovData[i].doubledSurfaceArea * mu /
-                                 seissol::dr::misc::numberOfBoundaryGaussPoints;
+        seismicMomentIncrease *=
+            godunovData[i].doubledSurfaceArea * mu / seissol::dr::misc::numberOfBoundaryGaussPoints;
         seismicMoment += seismicMomentIncrease;
       }
     }
@@ -200,7 +208,13 @@ void EnergyOutput::computeVolumeEnergies() {
   // Note: Default(none) is not possible, clang requires data sharing attribute for g, gcc forbids
   // it
 #if defined(_OPENMP) && !defined(__NVCOMPILER)
-#pragma omp parallel for schedule(static) reduction(+ : totalGravitationalEnergyLocal, totalAcousticEnergyLocal, totalAcousticKineticEnergyLocal, totalElasticEnergyLocal, totalElasticKineticEnergyLocal, totalPlasticMoment) shared(elements, vertices, lts, ltsLut, global)
+#pragma omp parallel for schedule(static) reduction(+ : totalGravitationalEnergyLocal,             \
+                                                        totalAcousticEnergyLocal,                  \
+                                                        totalAcousticKineticEnergyLocal,           \
+                                                        totalElasticEnergyLocal,                   \
+                                                        totalElasticKineticEnergyLocal,            \
+                                                        totalPlasticMoment)                        \
+    shared(elements, vertices, lts, ltsLut, global)
 #endif
   for (std::size_t elementId = 0; elementId < elements.size(); ++elementId) {
     real volume = MeshTools::volume(elements[elementId], vertices);

src/ResultWriter/EnergyOutput.h

@@ -63,7 +63,7 @@ class EnergyOutput : public Module {
                          const real* degreesOfFreedomMinus,
                          DRFaceInformation const& faceInfo,
                          DRGodunovData const& godunovData,
-                         const real slip[seissol::tensor::slipRateInterpolated::size()]);
+                         const real slip[seissol::tensor::slipInterpolated::size()]);
 
   void computeDynamicRuptureEnergies();