refactor stress components (introducing OptimalStress_Szz)

doc/maps.rst

@@ -206,24 +206,37 @@ Provides values by evaluating another easi tree.
   evaluate this intermediate variable before executing the
   "!STRESS\_STR\_DIP\_SLIP\_AM" component.
 
-OptimalStress
--------------
+OptimalStress_Szz
+------------------
+
+This component generates a stress tensor which maximizes shear traction
+on the optimally oriented plane defined by the strike and dip angles, along the rake angle orientation.
+Such optimally oriented plane can be a virtual plane if it does not correspond to any segment of the fault system. 
+
+The principal stress magnitudes are prescribed by:
+
+- the relative prestress ratio R (where :math:`R=1/(1+S)`, with S the relative fault strength), 
+- the stress shape ratio s2ratio :math:`= (s_2-s_3)/(s_1-s_3)`, where :math:`s_1`, :math:`s_2` and :math:`s_3` are the maximum, intermediate, and minimum compressive stress, respetively, and 
+- the :math:`s_{zz}` component of the stress tensor
 
-This function allows computing the stress which would result in faulting
-in the rake direction on the optimally oriented plane defined by strike
-and dip angles (this can be only a virtual plane if such optimal
-orientation does not correspond to any segment of the fault system). The
-principal stress magnitudes are prescribed by the relative prestress
-ratio R (where :math:`R=1/(1+S)`), the effective confining stress
-(effectiveConfiningStress :math:`= Tr(sii)/3`) and the stress shape ratio
-:math:`s2ratio = (s_2-s_3)/(s_1-s_3)`, where :math:`s_1>s_2>s_3` are the principal stress
-magnitudes, following the procedure described in Ulrich et al.
-(2019), methods section 'Initial Stress'. To prescribe R, static and dynamic friction
-(mu\_s and mu\_d) as well as cohesion are required. 
+
+To prescribe R, static and dynamic friction (mu\_s and mu\_d) as well as cohesion are required. 
+The procedure is described in Ulrich et al. (2019), methods section 'Initial Stress' 
+(note that in Ulrich et al. (2019), as in the deprecated `OptimalStress` component, the effectiveConfiningStress :math:`= Tr(s_{ii})/3` is prescribed instead of :math:`s_{zz}`).
+
+The principal stresses are oriented relatively to the strike, dip and rake angles as follow:
+
+.. image:: figs/OptimalStress_orientation.png
+    :width: 50%
+    :align: center
+
+The red plane, that contains :math:`s_1` and :math:`s_3` is normal to the optimally oriented fault plane and contains the rake vector.
+:math:`s_2` (not represented) is normal to the red plane.
+Note that the OptimalStress_Szz component is only valid for an ENU system.
 
 .. code-block:: YAML
 
-        components: !OptimalStress
+        components: !OptimalStress_Szz
           constants:
             mu_d:      <double>
             mu_s:      <double>
@@ -232,23 +245,23 @@ magnitudes, following the procedure described in Ulrich et al.
             rake:      <double>
             cohesion:  <double>
             s2ratio:   <double>
-            R:         <double>
-            effectiveConfiningStress: <double>
+            R0:         <double>
+            s_zz: <double>
 
 :Domain:
   *inherited*
 :Codomain:
-  stress components (s\_xx, s\_yy, s\_zz, s\_xy, s\_yz, and s\_xz)
+  stress components (b\_xx, b\_yy, b\_zz, b\_xy, b\_yz, and b\_xz)
 
 AndersonianStress
 -----------------
 
-This function allows computing Andersonian stresses (for which one principal axis of the stress tensor is vertical).
+This component allows computing Andersonian stress tensors (for which one principal axis of the stress tensor is vertical).
 
 The principal stress orientations are defined by SH_max (measured from North, positive eastwards), the direction of maximum horizontal compressive stress.
 
 S_v defines which of the principal stresses :math:`s_i` is vertical where :math:`s_1>s_2>s_3`.
-S_v = 1, 2 or 3 should be used if the vertical principal stress is the maximum, intermediate or minimum compressive stress.
+S_v = 1, 2, or 3 should be used when the vertical principal stress is the maximum, intermediate, or minimum compressive stress, respectively.
 Assuming mu_d=0.6, S_v = 1 favours normal faulting on a 60° dipping fault plane striking SH_max,
 S_v = 2 favours strike-slip faulting on a vertical fault plane making an angle of 30° with SH_max and
 S_v = 3 favours reverse faulting on a 30° dipping fault plane striking SH_max.
@@ -258,8 +271,7 @@ the vertical stress sig_zz and the stress shape ratio
 :math:`s2ratio = (s_2-s_3)/(s_1-s_3)`, where :math:`s_1>s_2>s_3` are the principal stress
 magnitudes, following the procedure described in Ulrich et al.
 (2019), methods section 'Initial Stress'. To prescribe S, static and dynamic friction
-(mu\_s and mu\_d) as well as cohesion are required. 
-
+(mu\_s and mu\_d), as well as cohesion, are required. 
 
 
 .. code-block:: YAML
@@ -278,41 +290,9 @@ magnitudes, following the procedure described in Ulrich et al.
 :Domain:
   *inherited*
 :Codomain:
-  stress components (s\_xx, s\_yy, s\_zz, s\_xy, s\_yz, and s\_xz)
-
+  stress components (b\_xx, b\_yy, b\_zz, b\_xy, b\_yz, and b\_xz)
 
 
-STRESS\_STR\_DIP\_SLIP\_AM (deprecated)
----------------------------------------
-
-This routine is now replaced by the more complete and exact
-'OptimalStress' routine. It is nevertheless preserved in the code for
-being able to run the exact setup we use for the Sumatra SC paper (Uphoff
-et al., 2017). It is mostly similar with the 'OptimalStress' routine,
-but instead of a rake parameter, the direction of slip can only be pure
-strike-slip and pure dip-slip faulting (depending on the parameter
-DipSlipFaulting). In this routine the s\_zz component of the stress
-tensor is prescribed (and not the confining stress tr(sii)/3) as in
-'OptimalStress'.
-
-.. code-block:: YAML
-
-        components: !STRESS_STR_DIP_SLIP_AM
-          constants:
-            mu_d:      <double>
-            mu_s:      <double>
-            strike:    <double>
-            dip:       <double>
-            DipSlipFaulting: <double> (0 or 1)
-            cohesion:  <double>
-            s2ratio:   <double>
-
-:Domain:
-  *inherited*
-:Codomain:
-  stress components (s\_xx, s\_yy, s\_zz, s\_xy, s\_yz, and s\_xz)
-:Example:
-  `120_initial_stress <https://github.com/SeisSol/easi/blob/9e93f35fbacc950d00534643c59a64dff306a381/examples/120_initial_stress.yaml#L44>`__
 
 SpecialMap
 ----------
@@ -376,3 +356,158 @@ Evaluates application-defined functions.
 
   The domain of !Special is now i1, i3 and the codomain is o1, o2.
   i2 is constant and has the value 3.
+
+
+Deprecated components
+---------------------
+
+OptimalStress
+^^^^^^^^^^^^^^
+
+This component does the same as the component `OptimalStress_Szz`, but prescribes the effectiveConfiningStress :math:`= Tr(s_{ii})/3`,
+rather than the :math:`s_{zz}` component of the stress tensor.
+Note that the OptimalStress component is only valid for an ENU system.
+
+.. code-block:: YAML
+
+        components: !OptimalStress
+          constants:
+            mu_d:      <double>
+            mu_s:      <double>
+            strike:    <double>
+            dip:       <double>
+            rake:      <double>
+            cohesion:  <double>
+            s2ratio:   <double>
+            R:         <double>
+            effectiveConfiningStress: <double>
+
+:Domain:
+  *inherited*
+:Codomain:
+  stress components (b\_xx, b\_yy, b\_zz, b\_xy, b\_yz, and b\_xz)
+
+The following code:
+
+.. code-block:: YAML
+
+    !OptimalStress
+        constants:
+            mu_d:     a
+            mu_s:     b
+            strike:   c
+            dip:      d
+            rake:     e
+            cohesion: f
+            s2ratio:  g
+            R:        h
+            effectiveConfiningStress:     i
+
+can be (in most case accurately) approximated by:
+
+.. code-block:: YAML
+
+    !EvalModel
+    parameters: [b_xx,b_yy,b_zz,b_xy, b_xz, b_yz]
+    model: !ConstantMap
+        map:
+                mu_d: a
+        components: !OptimalStress_Szz
+              constants:
+                mu_s:      b
+                strike:    c
+                dip:       d
+                rake:      e
+                cohesion:  f
+                s2ratio:   g
+                R0:        h
+                s_zz:      i
+    components: !FunctionMap
+        map:
+           b_xx: return b_xx*(3.*b_zz)/(b_xx+b_yy+b_zz);
+           b_yy: return b_yy*(3.*b_zz)/(b_xx+b_yy+b_zz);
+           b_zz: return b_zz*(3.*b_zz)/(b_xx+b_yy+b_zz);
+           b_xy: return b_xy*(3.*b_zz)/(b_xx+b_yy+b_zz);
+           b_xz: return b_xz*(3.*b_zz)/(b_xx+b_yy+b_zz);
+           b_yz: return b_yz*(3.*b_zz)/(b_xx+b_yy+b_zz);
+
+
+Note that this is not fully equivalent if `cohesion` is not small compared to `effectiveConfiningStress`.
+
+
+STRESS\_STR\_DIP\_SLIP\_AM
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. code-block:: YAML
+
+        components: !STRESS_STR_DIP_SLIP_AM
+          constants:
+            mu_d:      <double>
+            mu_s:      <double>
+            strike:    <double>
+            dip:       <double>
+            DipSlipFaulting: <double> (0 or 1)
+            cohesion:  <double>
+            s2ratio:   <double>
+            s_zz:      <double>
+            R:         <double>
+            s2ratio:   <double>
+
+:Domain:
+  *inherited*
+:Codomain:
+  normalized stress components (b\_xx, b\_yy, b\_zz, b\_xy, b\_yz, and b\_xz, with b\_zz=1)
+:Example:
+  `120_initial_stress <https://github.com/SeisSol/easi/blob/9e93f35fbacc950d00534643c59a64dff306a381/examples/120_initial_stress.yaml#L44>`__
+
+This component, used for example in the setup of the Sumatra SC paper 
+(Uphoff et al., 2017) is deprecated and can be substituted by the more complete 'OptimalStress_Szz' routine. 
+While in the 'OptimalStress_Szz' routine, a rake parameter defines the direction of maximized shear traction 
+on the optimally oriented fault, such direction is here defined by the parameter
+DipSlipFaulting (1 for pure dip-slip, 0 for pure strike-slip).
+Another difference with the OptimalStress component is that STRESS\_STR\_DIP\_SLIP\_AM returns a normalized stress tensor.
+Note that the STRESS\_STR\_DIP\_SLIP\_AM component is only valid for an ENU system.
+
+The following code:
+
+.. code-block:: YAML
+
+    !STRESS_STR_DIP_SLIP_AM
+        constants:
+            mu_d:     a
+            mu_s:     b
+            strike:   c
+            dip:      d
+            DipSlipFaulting:  1.0
+            cohesion: e
+            s2ratio:  f
+            R:        g
+            s_zz:     h
+
+is equivalent to:
+
+.. code-block:: YAML
+
+    !EvalModel
+    parameters: [b_xx,b_yy,b_zz,b_xy, b_xz, b_yz]
+    model: !ConstantMap
+        map:
+                mu_d: a
+        components: !OptimalStress_Szz
+              constants:
+                mu_s:      b
+                strike:    c
+                dip:       d
+                rake:     90.0
+                cohesion:  e
+                s2ratio:   f
+                R0:        g
+                s_zz:      h
+    components: !FunctionMap
+        map:
+           b_xx: return b_xx/b_zz;
+           b_yy: return b_yy/b_zz;
+           b_zz: return 1;
+           b_xy: return b_xy/b_zz;
+           b_xz: return b_xz/b_zz;
+           b_yz: return b_yz/b_zz;

include/easi/YAMLParser.h

@@ -115,6 +115,7 @@ YAMLParser::YAMLParser(unsigned dimDomain, AsagiReader* externalAsagiReader, cha
   // Specials
   registerSpecial<STRESS_STR_DIP_SLIP_AM>("!STRESS_STR_DIP_SLIP_AM");
   registerSpecial<OptimalStress>("!OptimalStress");
+  registerSpecial<OptimalStress_Szz>("!OptimalStress_Szz");
   registerSpecial<AndersonianStress>("!AndersonianStress");
 }