diff --git a/doc/figs/OptimalStress_orientation.png b/doc/figs/OptimalStress_orientation.png new file mode 100644 index 0000000..c41b4a1 Binary files /dev/null and b/doc/figs/OptimalStress_orientation.png differ diff --git a/doc/maps.rst b/doc/maps.rst index 381d898..7a62097 100644 --- a/doc/maps.rst +++ b/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: mu_s: @@ -232,23 +245,23 @@ magnitudes, following the procedure described in Ulrich et al. rake: cohesion: s2ratio: - R: - effectiveConfiningStress: + R0: + s_zz: :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: - mu_s: - strike: - dip: - DipSlipFaulting: (0 or 1) - cohesion: - s2ratio: - -:Domain: - *inherited* -:Codomain: - stress components (s\_xx, s\_yy, s\_zz, s\_xy, s\_yz, and s\_xz) -:Example: - `120_initial_stress `__ 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: + mu_s: + strike: + dip: + rake: + cohesion: + s2ratio: + R: + effectiveConfiningStress: + +: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: + mu_s: + strike: + dip: + DipSlipFaulting: (0 or 1) + cohesion: + s2ratio: + s_zz: + R: + s2ratio: + +: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 `__ + +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; diff --git a/include/easi/YAMLParser.h b/include/easi/YAMLParser.h index 964f424..e80be63 100644 --- a/include/easi/YAMLParser.h +++ b/include/easi/YAMLParser.h @@ -115,6 +115,7 @@ YAMLParser::YAMLParser(unsigned dimDomain, AsagiReader* externalAsagiReader, cha // Specials registerSpecial("!STRESS_STR_DIP_SLIP_AM"); registerSpecial("!OptimalStress"); + registerSpecial("!OptimalStress_Szz"); registerSpecial("!AndersonianStress"); } diff --git a/include/easi/component/OptimalStress_Szz.h b/include/easi/component/OptimalStress_Szz.h new file mode 100644 index 0000000..569ca59 --- /dev/null +++ b/include/easi/component/OptimalStress_Szz.h @@ -0,0 +1,113 @@ +/** + * @file + * This file is part of SeisSol. + * + * @author Thomas Ulrich (ulrich AT geophysik.uni-muenchen.de) + * + * @section LICENSE + * Copyright (c) 2018, SeisSol Group + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * 3. Neither the name of the copyright holder nor the names of its + * contributors may be used to endorse or promote products derived from this + * software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * + * @section DESCRIPTION + **/ +#ifndef EASI_COMPONENT_OPTIMALSTRESS_SZZ_H_ +#define EASI_COMPONENT_OPTIMALSTRESS_SZZ_H_ + +#include +#include "easi/util/MagicStruct.h" + +namespace easi { + struct OptimalStress_Szz { + struct in { + double mu_d, mu_s, strike, dip, rake, s_zz, cohesion, R0, s2ratio; + }; + in i; + + struct out { + double b_xx, b_yy, b_zz, b_xy, b_yz, b_xz; + }; + out o; + + inline void evaluate(); + }; +} + +SELF_AWARE_STRUCT(easi::OptimalStress_Szz::in, mu_d, mu_s, strike, dip, rake, s_zz, cohesion, R0, s2ratio) +SELF_AWARE_STRUCT(easi::OptimalStress_Szz::out, b_xx, b_yy, b_zz, b_xy, b_yz, b_xz) + +// COMPUTE NORMALIZED STRESS FOLLOWING THE METHOD OF Ulrich et al. (2018) +void easi::OptimalStress_Szz::evaluate() { + double const pi = 3.141592653589793; + // most favorable direction (A4, AM2003) + double Phi = pi/4.0-0.5*atan(i.mu_s); + double strike_rad = i.strike*pi/180.0; + double dip_rad = i.dip*pi/180.0; + double rake_rad = i.rake*pi/180.0; + double s2 = sin(2.0*Phi); + double c2 = cos(2.0*Phi); + double alpha = (2.0*i.s2ratio-1.0)/3.0; + + double cd = cos(-Phi); + double sd = sin(-Phi); + + double cr = cos(rake_rad); + double sr = sin(rake_rad); + + double ci = cos(dip_rad); + double si = sin(dip_rad); + + double cs = cos(strike_rad); + double ss = sin(strike_rad); + + // a1, a2 and a3 are the coefficient of o.b_zz = -a1 s11 -a2 s22 -a3 s33; + double a1 = pow(si * sr * cd + ci * sd, 2.0); + double a2 = si * si * cr * cr; + double a3 = pow(-si * sr * sd + ci * cd, 2.0); + double a = a1 + a2 + a3; + double b = a1 + (2.0*i.s2ratio-1.0)*a2 - a3; + double mufac = (i.mu_d + i.R0 * (i.mu_s-i.mu_d)); + double ds = (std::fabs(i.s_zz) + i.R0 * a * std::fabs(i.cohesion) / mufac) / (b + a * c2 + a *s2 / mufac); + double sm = (std::fabs(i.s_zz) - b * ds) / a; + + //sii are all positive + double s11 = sm + ds; + double s22= sm - ds + 2.0*ds*i.s2ratio; + double s33 = sm - ds; + + + o.b_xx = -((ci * ci * s11 * sr * sr + s33 * si * si) * cd * cd + 2 * ci * sd * si * sr * (s33 - s11) * cd + (s33 * sd * sd * sr * sr + cr * cr * s22) * ci * ci + sd * sd * si * si * s11) * cs * cs + 2 * cr * cs * (cd * cd * ci * sr * s11 + sd * si * (s33 - s11) * cd - ci * sr * (-s33 * sd * sd + s22)) * ss - ss * ss * (cd * cd * cr * cr * s11 + cr * cr * s33 * sd * sd + s22 * sr * sr); + o.b_xy = cr * (cd * cd * ci * sr * s11 + sd * si * (s33 - s11) * cd - ci * sr * (-s33 * sd * sd + s22)) * cs * cs + ss * ((-cd * cd * s11 + ci * ci * s22 - s33 * sd * sd) * cr * cr + (ci * ci * s11 * sr * sr + s33 * si * si) * cd * cd + 2 * ci * sd * si * sr * (s33 - s11) * cd + s33 * ci * ci * sd * sd * sr * sr + sd * sd * si * si * s11 - s22 * sr * sr) * cs - ss * ss * cr * (cd * cd * ci * sr * s11 + sd * si * (s33 - s11) * cd - ci * sr * (-s33 * sd * sd + s22)); + o.b_xz = cd * cs * sd * sr * (s33 - s11) * si * si - ((cs * (-s11 * sr * sr + s33) * ci + cr * sr * ss * s11) * cd * cd - cs * ((s33 * sr * sr - s11) * sd * sd + cr * cr * s22) * ci - cr * sr * ss * (-s33 * sd * sd + s22)) * si - cd * ci * sd * (ci * cs * sr - cr * ss) * (s33 - s11); + o.b_yy = -((ci * ci * s11 * sr * sr + s33 * si * si) * cd * cd + 2 * ci * sd * si * sr * (s33 - s11) * cd + (s33 * sd * sd * sr * sr + cr * cr * s22) * ci * ci + sd * sd * si * si * s11) * ss * ss - 2 * cr * cs * (cd * cd * ci * sr * s11 + sd * si * (s33 - s11) * cd - ci * sr * (-s33 * sd * sd + s22)) * ss - cs * cs * (cd * cd * cr * cr * s11 + cr * cr * s33 * sd * sd + s22 * sr * sr); + o.b_yz = -cd * sd * sr * ss * (s33 - s11) * si * si - ((-ss * (-s11 * sr * sr + s33) * ci + cr * cs * sr * s11) * cd * cd + ss * ((s33 * sr * sr - s11) * sd * sd + cr * cr * s22) * ci - cr * cs * sr * (-s33 * sd * sd + s22)) * si + cd * ci * sd * (ci * sr * ss + cr * cs) * (s33 - s11); + o.b_zz = -(cd * cd * s11 * sr * sr + s33 * sd * sd * sr * sr + cr * cr * s22) * si * si - 2 * cd * ci * sd * sr * (s11 - s33) * si - ci * ci * (cd * cd * s33 + s11 * sd * sd); + +} + +#endif diff --git a/include/easi/component/Special.h b/include/easi/component/STRESS_STR_DIP_SLIP_AM.h similarity index 100% rename from include/easi/component/Special.h rename to include/easi/component/STRESS_STR_DIP_SLIP_AM.h diff --git a/include/easi/parser/YAMLComponentParsers.h b/include/easi/parser/YAMLComponentParsers.h index 2b79fee..0e51b86 100644 --- a/include/easi/parser/YAMLComponentParsers.h +++ b/include/easi/parser/YAMLComponentParsers.h @@ -52,8 +52,9 @@ #include "easi/component/EvalModel.h" #include "easi/component/LayeredModelBuilder.h" #include "easi/component/SpecialMap.h" -#include "easi/component/Special.h" +#include "easi/component/STRESS_STR_DIP_SLIP_AM.h" #include "easi/component/OptimalStress.h" +#include "easi/component/OptimalStress_Szz.h" #include "easi/component/AndersonianStress.h" #ifdef USE_ASAGI