From 0da95db32af509e1bebca6378ecf6e83fecb9d76 Mon Sep 17 00:00:00 2001
From: Adam Beall <abea0020@mu00186821x.gateway>
Date: Sat, 19 Feb 2022 19:27:29 +1100
Subject: [PATCH] added single asperity mode

---
 README.md |  2 ++
 plot.py   |  2 --
 visc.py   | 47 +++++++++++++++++++++++++++++++----------------
 3 files changed, 33 insertions(+), 18 deletions(-)

diff --git a/README.md b/README.md
index 5f2f19f..543e092 100644
--- a/README.md
+++ b/README.md
@@ -6,6 +6,8 @@ by Adam Beall, Martijn van den Ende, Jean-Paul Ampuero and Ake Fagereng. These s
 The script 'visc.py' runs the earthquake cycle model, where the inputs can be varied to reproduce any of the heterogenous fault models in the study.
 Use 'plot.py' to colalte the synthetic earthquake catalogue and visualise the model output.
 
+The script 'visc.py' can also be used to run the single asperity test model, by setting 'visc_dist' to 'single'.
+
 
 Requirements:
 - QDYN (and associated requirements), with the viscosity branch:
diff --git a/plot.py b/plot.py
index 3957699..2aca4ce 100644
--- a/plot.py
+++ b/plot.py
@@ -101,7 +101,6 @@
 arrL = []
 etathresh = 60.0e6 * W / 1e-9
 stress_drop = 0e6 * W / 1e-9
-print("%.2e" %etathresh)
 L = 0.0
 
 arrActivepatchesX = [] 
@@ -202,7 +201,6 @@
 
 for i,t in enumerate(arrTime):
     arrI[i] = np.argmin( np.abs(time*t_yr-t))
-    print(arrI[i])
 
 
 for i in arrI:
diff --git a/visc.py b/visc.py
index ab7f386..9de5975 100755
--- a/visc.py
+++ b/visc.py
@@ -26,7 +26,11 @@ def load():
     # ---- Key paramaters
     # Shear-zone thickness
     params['W'] = 100.0
-    params['visc_contrast'] = 100.0
+    # Shear-zone viscosity min/max in Pa s
+    params['visc_min'] = 1e18
+    params['visc_max'] = 1e20
+
+
     # seed for patch size and viscosity random generation 
     np.random.seed(3)  # model-set A
     #np.random.seed(5)  # model-set B
@@ -36,10 +40,13 @@ def load():
     # Model duration
     params['maxtime'] = 1000 * 3600 * 24 * 365.0    # model time in seconds
 
-    # options are log, powerlaw or bimodal
+    # options are log, powerlaw, bimodal or single
+    # single = single asperity
     params['visc_dist'] = 'log' 
 
 
+    # Patch width - only used for single asperity model
+    params['asp_width'] = 5e3 
 
 
     # Instantiate the QDYN class object
@@ -113,8 +120,8 @@ def load():
 
     # Generate random viscosities
 
-    eta_min = 18   # min. log(eta)
-    eta_max = eta_min + np.log10(params['visc_contrast'])   # max. log(eta)
+    eta_min = np.log10(params['visc_min'])
+    eta_max = np.log10(params['visc_max'])
 
 
     xrand = np.random.random(200)
@@ -128,25 +135,33 @@ def load():
         arrStressRand += 10.**eta_min
     elif params['visc_dist'] == 'powerlaw':
         arrStressRand = (eta_min**D + (eta_max**D - eta_min**D)*xrand)**(1./D)
-    else:
+    elif params['visc_dist'] != 'single':
         raise NameError('No viscosity distribution chosen')
 
 
     # map patches to fault x co-ordinates
-    arrWCumul = np.cumsum(arrW)
-
     arrX = np.linspace(0,params['L'],params['res'])
-    arrPert = np.zeros(len(arrX))
     arrMaterial = np.zeros(arrX.size)
 
-    pCount = 0
-    for i,xi in enumerate(arrX):
-        if xi > arrWCumul[pCount+1]:
-            pCount += 1
-        material = pCount % 2
-        arrMaterial[i] =  arrStressRand[pCount+1]
-            
-    arrW = arrW[:pCount]
+    if params['visc_dist'] == 'single':
+        x_l = 0.5 * (params['L'] - params['asp_width'])
+        x_r = 0.5 * (params['L'] + params['asp_width'])
+        for i,xi in enumerate(arrX):
+            if xi < x_l or xi>x_r:
+                arrMaterial[i] = 10. ** eta_min
+            else:
+                arrMaterial[i] = 10. ** eta_max
+
+    else:
+        arrWCumul = np.cumsum(arrW)
+
+        pCount = 0
+        for i,xi in enumerate(arrX):
+            if xi > arrWCumul[pCount+1]:
+                pCount += 1
+            material = pCount % 2
+            arrMaterial[i] =  arrStressRand[pCount+1]
+                
 
 
     # Set viscosity