strain softening (#106)

* strain softening from interfaces.txt working for cohesion and friction angle but not for weakening seed * read weakening_seed from interfaces.txt and if negative, use initial random strain * calc_visc checks if strain softenins is from interfaces.txt or command line; corrected where 'p' was used insted of 'layer_array[p]' * small if correction * fixed problem where interfaces.txt required strain softening parameters regardless of the strain softening input * uses weakening_min as weakening_seed if no parameters are given in the interfaces.txt file; also uses Beaumont values for cohesion and friction angle * small changes * add missing 'make all' command to build mandyoc image (#110) * fix typo in line break character (#111) * change order of parse options function (#112) * strain softening from interfaces.txt working for cohesion and friction angle but not for weakening seed * read weakening_seed from interfaces.txt and if negative, use initial random strain * calc_visc checks if strain softenins is from interfaces.txt or command line; corrected where 'p' was used insted of 'layer_array[p]' * small if correction * fixed problem where interfaces.txt required strain softening parameters regardless of the strain softening input * uses weakening_min as weakening_seed if no parameters are given in the interfaces.txt file; also uses Beaumont values for cohesion and friction angle * small changes * fixed previous issues where the interfaces files was read incorrectly and changed code to fit the parse_options() commign before the reader() * calc_visc conflict solving * solved conflict with SMSwarm_3d.cpp * fix check of number of layers read from seed command and interfaces file (#114) * strain softening from interfaces.txt working for cohesion and friction angle but not for weakening seed * read weakening_seed from interfaces.txt and if negative, use initial random strain * calc_visc checks if strain softenins is from interfaces.txt or command line; corrected where 'p' was used insted of 'layer_array[p]' * small if correction * fixed problem where interfaces.txt required strain softening parameters regardless of the strain softening input * uses weakening_min as weakening_seed if no parameters are given in the interfaces.txt file; also uses Beaumont values for cohesion and friction angle * small changes * fixed previous issues where the interfaces files was read incorrectly and changed code to fit the parse_options() commign before the reader() * strain softening from interfaces.txt working for cohesion and friction angle but not for weakening seed * read weakening_seed from interfaces.txt and if negative, use initial random strain * calc_visc checks if strain softenins is from interfaces.txt or command line; corrected where 'p' was used insted of 'layer_array[p]' * small if correction * fixed problem where interfaces.txt required strain softening parameters regardless of the strain softening input * solved small conflict * fixed weakening arrays construction * default weakening start value is 0.0 * Update src/reader.cpp to truncate friction angle values Co-authored-by: Rafael Silva <rafaelmds@users.noreply.github.com> * set default weakening_seed to be set by rand_r() * documentation + reference changing --------- Co-authored-by: Rafael Silva <rafaelmds@users.noreply.github.com>
ggciag · Apr 13, 2023 · 66afb7c · 66afb7c
1 parent 8ce1f00
commit 66afb7c
Show file tree

Hide file tree

Showing 11 changed files with 330 additions and 214 deletions.
diff --git a/docs/files/implementation.rst b/docs/files/implementation.rst
@@ -283,7 +283,7 @@ The plastic deformation can be calculated using the Byerlee Law :cite:`byerlee19
 
     \tau_{yield} = c_{0}+\mu \rho g z
 
-where :math:`c_{0}` in the internal cohesion, :math:`\mu` is the friction coefficient, :math:`\rho` is the density and :math:`z` is the depth.
+where :math:`c_{0}` in the internal cohesion of the rock, :math:`\mu` is the friction coefficient, :math:`\rho` is the density and :math:`z` is the depth.
 
 Alternatively, the user can choose to use the the Druker-Prager criterion :cite:`drucker-prager1952`, which is presented in :eq:`drucker-prager`.
 
@@ -306,7 +306,7 @@ The ductile rheology can be simulated using the Frank-Kamenetskii approximation,
 
     \eta_{visc}(T) = C \eta_r b^* \exp{(-\gamma T)}
 
-where :math:`\eta_r` is the reference viscosity, :math:`C`is a compositional factor to scale the effective viscosity, and :math:`b^*` and :math:`\gamma = E_a / RT^2_b` are constants, which in turn, :math:`E_a` is the activation energy, :math:`R` is the gas constant, and :math:`T_b` is the basal temperature.
+where :math:`\eta_r` is the reference viscosity, :math:`C` is a compositional factor to scale the effective viscosity, and :math:`b^*` and :math:`\gamma = E_a / RT^2_b` are constants, which in turn, :math:`E_a` is the activation energy, :math:`R` is the gas constant, and :math:`T_b` is the basal temperature.
 
 Additionally, the rheology can also be considered to follow a power law, as a function of the temperature :math:`T`, compositional factor :math:`C`, pressure :math:`P` and strain rate :math:`\varepsilon` as follows:
 

diff --git a/docs/files/input.rst b/docs/files/input.rst
@@ -79,7 +79,7 @@ The ``interfaces.txt`` file, for a 2-D grid, starts with seven lines of variable
 2-D Initial interfaces
 **********************
 
-Below, the example corresponds to a 2-D grid with two interfaces and, therefore, three lithological units. The first column contains the vertical positions **in meters** of every grid node :math:`y_m` that corresponds to the  **deepest** interface boundary, starting at :math:`x=x_0` on line 8, and ending at :math:`x=x_{nx-1}` on line :math:`nx+7`. The second column contains the vertical position of every :math:`y_m` that corresponds to the second interface boundary. When defining the interfaces, it is rather common for them to "touch". Because of that, all the interfaces must be provided in a "tetris" manner, where interfaces that are collinear in parts fit the interface below.
+Below, the example corresponds to a 2-D grid with two interfaces and, therefore, three lithological units. The first column contains the vertical positions **in meters** of every grid node :math:`y_m` that corresponds to the **deepest** interface boundary, starting at :math:`x=x_0` on line 8, and ending at :math:`x=x_{nx-1}` on line :math:`nx+7`. The second column contains the vertical position of every :math:`y_m` that corresponds to the second interface boundary. When defining the interfaces, it is rather common for them to "touch". Because of that, all the interfaces must be provided in a "tetris" manner, where interfaces that are collinear in parts fit the interface below.
 
 .. note::
    Because the interfaces are defined linearly between the nodes, it is important to define them properly, so every point inside the grid can be attributed to a lithological unit.
@@ -90,6 +90,16 @@ Below, the example corresponds to a 2-D grid with two interfaces and, therefore,
 
 The values of the variables in the first seven lines are the values of the lithological units bound by the interfaces. For the 2-D grid with two interfaces, the first values of each variable refers to the lithological unit below the first interface, the second value of each variable refers to the lithological unit between the first and second interfaces, and the third value of each variable refers to the lithological unit above the second interface. If more interfaces are added, there will be more units bounded by an upper and a lower interface.
 
+To consider different values of internal cohesion :math:`c_0` and internal angle of friction :math:`\varphi` for each layer, the interface file can be written as the file below, where *weakening_seed*, *cohesion_min*, *cohesion_max*, *friction_angle_min*, and *friction_angle_max* are given.
+
+.. literalinclude:: src/new-initial-interfaces-2d.txt
+   :language: text
+   :linenos:
+
+The *weakening_seed* value for each layer represents its initial accumulated strain value (seed), and when it is negative, a random seed will be used instead. When *cohesion_min* and *cohesion_max* are different, that corresponding layer will be linearly softened between those values. Notice that a constant cohesion can be set if the values are simply equal. The same logic aplies to *friction_angle_min* and *friction_angle_max*.
+
+The default interval of accumualted strain where strain softening occurs is 0.05 and 1.05, cf. :cite:`salazarmora2018`. To use different values, the user must provide them in the ``param.txt`` withe the keywords *weakeing_min* and *weakening_max*.
+
 ..
    3-D Initial interfaces
    **********************

diff --git a/docs/files/parameter-file.rst b/docs/files/parameter-file.rst
@@ -332,6 +332,19 @@ The parameter file ``param.txt`` contains the information that is necessary for
         unit: :raw-html:`<br />`
         definition: set if radiogenic heating is active 
 
+#. Strain softening parameters
+
+    * weakeing_min                 
+        default: 0.05 :raw-html:`<br />`
+        type: real number :raw-html:`<br />`
+        unit: :raw-html:`<br />`
+        definition: value of the accumulated strain where strain softening starts :cite:`salazarmora2018`
+    * weakening_max             
+        default: 1.05 :raw-html:`<br />`
+        type: real number :raw-html:`<br />`
+        unit: :raw-html:`<br />`
+        definition: value of the accumulated strain where strain softening stops :cite:`salazarmora2018`
+
 #. Velocity boundary conditions
 
     * top_normal_velocity                 

diff --git a/docs/files/refs.bib b/docs/files/refs.bib
@@ -259,3 +259,12 @@ @article{gerya2003
   publisher={Elsevier}
 }
 
+@article{salazarmora2018,
+  title = {The Wilson Cycle and Effects of Tectonic Structural Inheritance on Rifted Passive Margin Formation},
+  author = {Salazar-Mora, Claudio A. and Huismans, Ritske S. and Fossen, Haakon and Egydio-Silva, Marcos},
+  journal = {Tectonics},  
+  volume = {37},
+  number = {9},
+  pages = {3085-3101},
+  year = {2018}
+}
diff --git a/docs/files/src/new-initial-interfaces-2d.txt b/docs/files/src/new-initial-interfaces-2d.txt
@@ -0,0 +1,19 @@
+C    				1.0    1.0    1.0
+rho  				3000.0 3000.0 3000.0
+H    				0.0    0.0    0.0
+A    				0.0    0.0    0.0
+n    				0.0    0.0    0.0
+Q    				0.0    0.0    0.0
+V    				0.0    0.0    0.0
+weakeing_seed		1.0    -1.0   0.5
+cohesion_min		5.0e6  5.0e6  5.0e6
+cohesion_max		30.0e7 5.0e6  1.0e7
+friction_angle_min	0.0    0.2    10.0
+friction_angle_max	0.0    20.0   10.0
+y_1(x_0)    y_2(x_0)
+y_1(x_1)    y_2(x_1)
+y_1(x_2)    y_2(x_2)
+y_1(x_3)    y_2(x_3)
+y_1(x_4)    y_2(x_4)
+y_1(x_5)    y_2(x_5)
+y_1(x_6)    y_2(x_6)
diff --git a/src/DMSwarm_2d.cpp b/src/DMSwarm_2d.cpp
@@ -64,6 +64,14 @@ extern PetscInt seed_layer_size;
 extern PetscBool seed_layer_set;
 extern PetscBool strain_seed_constant;
 
+extern PetscBool weakening_from_interfaces_file;
+
+// From param.txt and interfaces.txt
+extern PetscInt WITH_NON_LINEAR;
+extern PetscInt PLASTICITY;
+extern PetscScalar *weakening_seed;
+extern PetscScalar weakening_min;
+
 extern PetscReal *strain_seed_layer;
 extern PetscInt strain_seed_layer_size;
 extern PetscBool strain_seed_layer_set;
@@ -442,31 +450,13 @@ PetscErrorCode createSwarm_2d()
 					layer_array[p] = n_interfaces;
 				}
 
+				// Start initial strain_array with random value and overwrite it (if set)
 				rand_r(&seed_strain);
-				strain_array[p]=random_initial_strain*(float)rand_r(&seed_strain)/RAND_MAX;
-
-				if (seed_layer_set == PETSC_TRUE) {
-					if (seed_layer_size == 1 && layer_array[p] == seed_layer[0]) {
-						if (strain_seed_constant == PETSC_TRUE) {
-							strain_array[p] = strain_seed_layer[0];
-						} else {
-							strain_array[p] += strain_seed_layer[0];
-						}
-					}
-					else {
-						for (int k = 0; k < seed_layer_size; k++) {
-							if (layer_array[p] == seed_layer[k]) {
-								if (strain_seed_constant == PETSC_TRUE) {
-									strain_array[p] = strain_seed_layer[k];
-								} else {
-									strain_array[p] += strain_seed_layer[k];
-								}
-							}
-						}
-					}
-				}
-
+				strain_array[p] = random_initial_strain*(float)rand_r(&seed_strain)/RAND_MAX;
 
+				if (weakening_seed[layer_array[p]] >= 0) {
+					strain_array[p] = weakening_seed[layer_array[p]];
+				}
 			}
 		}
 

diff --git a/src/DMSwarm_3d.cpp b/src/DMSwarm_3d.cpp
@@ -65,6 +65,10 @@ extern PetscInt print_step_files;
 
 extern PetscReal random_initial_strain;
 
+// From param.txt and interfaces.txt
+extern PetscInt WITH_NON_LINEAR;
+extern PetscInt PLASTICITY;
+extern PetscScalar *weakening_seed;
 
 extern PetscInt *seed_layer;
 extern PetscInt seed_layer_size;
@@ -475,28 +479,12 @@ PetscErrorCode createSwarm_3d()
 					//printf("entrei!\n");
 				}
 
+				// Start initial strain_array with random value and overwrite it (if set)
 				rand_r(&seed_strain);
-				strain_array[p]=random_initial_strain*(float)rand_r(&seed_strain)/RAND_MAX;
-
-				if (seed_layer_set == PETSC_TRUE) {
-					if (seed_layer_size == 1 && layer_array[p] == seed_layer[0]) {
-						if (strain_seed_constant == PETSC_TRUE) {
-							strain_array[p] = strain_seed_layer[0];
-						} else {
-							strain_array[p] += strain_seed_layer[0];
-						}
-					}
-					else {
-						for (int k = 0; k < seed_layer_size; k++) {
-							if (layer_array[p] == seed_layer[k]) {
-								if (strain_seed_constant == PETSC_TRUE) {
-									strain_array[p] = strain_seed_layer[k];
-								} else {
-									strain_array[p] += strain_seed_layer[k];
-								}
-							}
-						}
-					}
+				strain_array[p] = random_initial_strain*(float)rand_r(&seed_strain)/RAND_MAX;
+
+				if (weakening_seed[layer_array[p]] >= 0) {
+					strain_array[p] = weakening_seed[layer_array[p]];
 				}
 				/////!!!!
 				/*if (rank==0){

diff --git a/src/DMSwarm_move.cpp b/src/DMSwarm_move.cpp
@@ -18,7 +18,7 @@ typedef struct {
 } Stokes3d;
 
 double calc_visco_ponto(double T, double P, double x, double z,double geoq_ponto,double e2_inva,double strain_cumulate,
-						double A, double n_exp, double QE, double VE);
+						double A, double n_exp, double QE, double VE, PetscInt layer_number);
 
 extern DM dms;
 
@@ -383,7 +383,8 @@ PetscErrorCode moveSwarm(int dimensions, PetscReal dt)
 
 
 			rarray[p] = calc_visco_ponto(tp,Pp,cx,cz,inter_geoq[layer_array[p]],E2_invariant,strain_fac[p],
-										inter_A[layer_array[p]], inter_n[layer_array[p]], inter_Q[layer_array[p]], inter_V[layer_array[p]]);
+										inter_A[layer_array[p]], inter_n[layer_array[p]], inter_Q[layer_array[p]], inter_V[layer_array[p]], layer_array[p]);
+			// PetscPrintf(PETSC_COMM_WORLD, "p: %d, strain: %E\n", p, strain_fac[p]);
 
 
 
@@ -575,8 +576,10 @@ PetscErrorCode moveSwarm(int dimensions, PetscReal dt)
 
 
 			rarray[p] = calc_visco_ponto(tp,Pp,cx,cz,inter_geoq[layer_array[p]],E2_invariant,strain_fac[p],
-										inter_A[layer_array[p]], inter_n[layer_array[p]], inter_Q[layer_array[p]], inter_V[layer_array[p]]);
+										inter_A[layer_array[p]], inter_n[layer_array[p]], inter_Q[layer_array[p]], inter_V[layer_array[p]], layer_array[p]);
 
+			// PetscPrintf(PETSC_COMM_WORLD, "p: %d, strain: %E\n", p, strain_fac[p]);
+
 			array[3*p  ] += dt * vx;
 			array[3*p+1] += dt * vy;
 			array[3*p+2] += dt * vz;