[cig-commits] r19277 - seismo/2D/SPECFEM2D/trunk/src/specfem2D

[dkomati1 at geodynamics.org]

Author: dkomati1
Date: 2011-12-06 17:42:53 -0800 (Tue, 06 Dec 2011)
New Revision: 19277

Modified:
   seismo/2D/SPECFEM2D/trunk/src/specfem2D/attenuation_compute_param.c
   seismo/2D/SPECFEM2D/trunk/src/specfem2D/gmat01.f90
Log:
fixed a division by zero problem in the attenuation routine when called from elements with no attenuation (e.g. acoustic elements)


Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/attenuation_compute_param.c
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/attenuation_compute_param.c	2011-12-07 00:29:58 UTC (rev 19276)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/attenuation_compute_param.c	2011-12-07 01:42:53 UTC (rev 19277)
@@ -72,12 +72,12 @@
     printf("T2 > T1\n");
     exit; }
 
-  if (target_Q1 <= 0.0001) {
-    printf("Q1 cannot be negative or null\n");
+  if (target_Q1 <= -0.0001) {
+    printf("Q1 cannot be negative\n");
     exit; }
 
-  if (target_Q2 <= 0.0001) {
-    printf("Q2 cannot be negative or null\n");
+  if (target_Q2 <= -0.0001) {
+    printf("Q2 cannot be negative\n");
     exit; }
 
   if (n < 1) {
@@ -112,6 +112,9 @@
     if (nu == 1) { Q_value = target_Q1 ; }
     if (nu == 2) { Q_value = target_Q2 ; }
 
+/* no need to compute these parameters if there is no attenuation; it could lead to a division by zero in the code */
+    if (Q_value > 0.00001) {
+
     tau_s = dvector(1, n);
     tau_e = dvector(1, n);
 
@@ -153,6 +156,8 @@
 
   }
 
+  }
+
 }
 
 void   plot_modulus(f1, f2, n, m, mR, Q, tau_e, tau_s ,xmgr)

Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/gmat01.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/gmat01.f90	2011-12-07 00:29:58 UTC (rev 19276)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/gmat01.f90	2011-12-07 01:42:53 UTC (rev 19277)
@@ -253,9 +253,6 @@
         aniso_array(4,n) = c33
         aniso_array(5,n) = c35
         aniso_array(6,n) = c55
-! dummy Q values, trick to avoid a bug in attenuation_model
-        QKappa_array(n) = 15
-        Qmu_array(n) = 15
         porosity_array(n) = 0.d0
      elseif (indic == 3) then
         density_array(1,n) = density(1)
@@ -274,8 +271,6 @@
         poroelastcoef(2,3,n) = mu_fr
         poroelastcoef(3,3,n) = lambdaplus2mu_fr
         poroelastcoef(4,3,n) = zero
-        QKappa_array(n) = 10.d0 ! dummy for attenuation_model
-        Qmu_array(n) = Qmu
      endif
 
      !