strain softening

docker/mandyoc/Dockerfile

@@ -32,6 +32,7 @@ COPY --chown=${USER}:${USER} . ${MANDYOC_DIR}
 RUN mkdir -p ${HOME}/.local/bin && \
     cd ${MANDYOC_DIR} && \
     make clear && \
+    make all && \
     make install
 
 # =============================================================================

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:
 

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
    **********************

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                 

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}
+}

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)

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]];
+				}
 			}
 		}
 

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){

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;