Reorganize examples

README.md

@@ -55,23 +55,31 @@ include("test/runtests.jl");
 ```
 
 # Getting Started:
-To run the code for the <a href="https://strike.scec.org/cvws/seas/download/SEAS_BP1_QD.pdf">first SEAS benchmark problem 
 
-```
-localARGS = ["examples/bp1-qd.dat"]
-include("examples/stripped_qd_driver.jl");
+### Benchmark Problem 1: 
+To run the example code for the <a href="https://strike.scec.org/cvws/seas/download/SEAS_BP1_QD.pdf">first SEAS benchmark problem 
 
 ```
-To run the code for <a href="https://strike.scec.org/cvws/seas/download/SEAS_BP6.pdf">BP6
-
+localARGS = ["examples/BP1/bp1-qd.dat"]
+include("examples/BP1/stripped_qd_driver.jl");
 ```
-localARGS = ["examples/bp6.dat"]
-include("examples/BP6_driver.jl");
 
+Note that the parameters set in the input file "examples/bp1-qd.dat" are slightly modified from the benchmark description to allow for fast simulation on most personal computers. The driver file calls a stripped down version of the code (for training purposes).
+
+### Benchmark Problem 6: 
+To run the example code for the <a href="https://strike.scec.org/cvws/seas/download/SEAS_BP6.pdf">sixth SEAS benchmark problem
 
 ```
-Note that the parameters set in the input file "examples/bp1-qd.dat" are slightly modified from the benchmark description to allow for fast simulation on most personal computers. The driver file calls a stripped down version of the code (for training purposes).
-
+localARGS = ["examples/BP6/bp6.dat"]
+include("examples/BP6/BP6_driver.jl");
+```
+
+### Benchmark Problem 5:
+To run the example code for the <a href="https://strike.scec.org/cvws/seas/download/SEAS_BP5.pdf">fifth SEAS benchmark problem</a> you will need GPU resources to solve in 3 dimensions, without a GPU it will be able to complete in a timely manner
+```
+include("examples/BP5/BP5-QD.jl");
+```
+
 # Documentation: 
 <a href="https://thrase.github.io/Thrase.jl/dev/"> More to come soon...
 

docs/make.jl

@@ -16,8 +16,9 @@ makedocs(;
     ),
     pages=[
         "Home" => "index.md",
-        "Numerical Method" => "method.md",
         "Install and Quick Start" => "quickstart.md",
+        "Examples" => "examples.md",
+        "Numerical Methods" => "method.md",
         "Input Parameters" => "parameters.md",
     ],
 )

docs/src/examples.md

@@ -0,0 +1,28 @@
+
+Once you have installed the Thrase project you should be able to run any of the below examples. If you want to change any input parameters you can do so by updating the .dat files for each benchmark problem.
+
+# Benchmark Problem 1
+To run the code for the first <a href="https://strike.scec.org/cvws/seas/download/SEAS_BP1_QD.pdf">SEAS benchmark problem
+
+```
+localARGS = ["examples/bp1-qd.dat"]
+include("examples/stripped_qd_driver.jl");
+```
+
+Note that the parameters set in the input file "examples/bp1-qd.dat" are slightly modified from the benchmark description to allow for fast simulation on most personal computers. The driver file calls a stripped down version of the code (for training purposes).
+
+
+# Benchmark Problem 6
+To run the example code for the <a href="https://strike.scec.org/cvws/seas/download/SEAS_BP6.pdf">sixth SEAS benchmark problem
+
+```
+localARGS = ["examples/BP6/bp6.dat"]
+include("examples/BP6/BP6_driver.jl");
+```
+
+# Benchmark Problem 5 (GPU):
+To run the example code for the <a href="https://strike.scec.org/cvws/seas/download/SEAS_BP5.pdf">fifth SEAS benchmark problem</a> you will need GPU resources to solve in 3 dimensions, without a GPU it will not be able to complete in a timely manner
+```
+include("examples/BP5/BP5-QD.jl");
+```
+The setup of the domains and coefficients for the BP5 example solutions is done in "examples/BP5/domain.jl" and "examples/BP5/domain_256.jl" for the different problem sizes respectively. To change values from the default you only need to update the constructor values in these files (or any domain size variables in the file).

docs/src/quickstart.md

@@ -31,12 +31,3 @@ julia --project=.
 ```
 include("test/runtests.jl");
 ```
-
-# Getting Started
-To run the code for the first [SEAS benchmark problem](https://strike.scec.org/cvws/seas/download/SEAS_BP1_QD.pdf)
-
-```
-localARGS = ["examples/bp1-qd.dat"]
-include("examples/stripped_qd_driver.jl");
-```
-Note that the parameters set in the input file "examples/bp1-qd.dat" are slightly modified from the benchmark description to allow for fast simulation on most personal computers. The driver file calls a stripped down version of the code (for training purposes).

examples/stripped_qd_driver.jl → examples/BP1/stripped_qd_driver.jl

@@ -99,7 +99,7 @@ function main()
           station_ind[i]
         end
         return Integer.(station_ind)
-      end
+    end
 
     station_indices = find_station_index(stations, z)
     station_strings = ["000", "025", "005", "075", "100", "125", "150", "175", "200", "250", "300", "350"] # "125" corresponds to 12.5 km down dip; these are necessary for writing to files
@@ -146,7 +146,7 @@ function main()
     prob = ODEProblem(odefun_stripped, ψδ, tspan, odeparam)
 
     # Set call-back function so that files are written to after successful time steps only.
-    cb_fun = SavingCallback((ψδ, t, i) -> write_to_file(pth, ψδ, t, i, z, flt_loc, flt_loc_indices,station_strings, station_indices, odeparam, "BP1_", 0.1 * year_seconds), SavedValues(Float64, Float64))
+    cb_fun = SavingCallback((ψδ, t, i) -> write_to_file(pth, ψδ, t, i, z, flt_loc, flt_loc_indices,station_strings, station_indices, odeparam, sim_years, "BP1_", 0.1 * year_seconds), SavedValues(Float64, Float64))
 
     # Make text files to store on-fault time series and slip data,
     # Also initialize with initial data:

examples/BP5/BP5-QD.jl

@@ -0,0 +1,149 @@
+# loading helper functions for calculations
+include("../../src/helper.jl")
+
+domain_file = "domain.jl"
+include(domain_file)
+
+# loading odefun defined for ODEProblem
+include("../../src/utils_MG.jl")
+include("../../src/odefun_BP5.jl")
+
+
+CUDA.allowscalar(false)
+
+function main()
+    time_string = Dates.format(now(), "yyyymmddHHMM")
+    @show BP5_coeff.Δz
+
+    # setting b = b0
+    b = BP5_coeff.b0
+    Vzero = 1e-20
+
+    @unpack_namedtuple odeparam
+
+    # calling create_ψVδ to create variables used for odefun, create_ψVδ defined in domain.jl
+    dψV, ψδ = create_ψVδ()
+    # calling create_view to create "views" for dψ, V, ψ, δ, create_view defined in domain.jl
+    dψ, V, ψ, δ = create_view(dψV, ψδ)
+    θ = zeros(size(ψ))
+
+    for i in 1:fN2
+        for j in 1:fN3
+            index = i + (j - 1) * fN2
+            x2 = (i - 1) * BP5_coeff.Δz / 1000 - BP5_coeff.lf / 2
+            x3 = (j - 1) * BP5_coeff.Δz / 1000
+            RSas[index] = a_func(x2, x3, BP5_coeff)
+        end
+    end
+
+    # initializing \boldsymbol{V} over the entire region
+    for i in 1:Ny
+        for j in 1:Nz
+            index = i + (j - 1) * Ny
+            V[2*index-1] = BP5_coeff.Vinit
+            V[2*index] = Vzero
+        end
+    end
+
+
+
+    # initializing \boldsymbol{τ}^0 for the entire domain
+    V_norm = norm([BP5_coeff.Vinit, Vzero])
+    τ = @view τb[1:2:length(τb)]
+    τz = @view τb[2:2:length(τb)]
+
+    # only using \tau values for the RS zone
+    for i in 1:fN2
+        for j in 1:fN3
+            index = i + (j - 1) * fN2
+            tau_index = RS_filter_2D_nzind[index]
+            τ0 = BP5_coeff.σn * RSas[index] * asinh((BP5_coeff.Vinit / (2 * BP5_coeff.V0) *
+                                                     exp((BP5_coeff.f0 + BP5_coeff.b0 * log(BP5_coeff.V0 / BP5_coeff.Vinit)) /
+                                                         RSas[index])) + η * BP5_coeff.Vinit)
+            τ[tau_index] = τ0 * BP5_coeff.Vinit / V_norm
+            τz[tau_index] = τ0 * Vzero / V_norm
+
+            θ0 = BP5_coeff.L / BP5_coeff.V0 * exp(RSas[index] / BP5_coeff.b0 *
+                                                  log(2 * BP5_coeff.V0 / BP5_coeff.Vinit * sinh((τ0 - η * BP5_coeff.Vinit) / (RSas[index] * BP5_coeff.σn)))
+                                                  -
+                                                  BP5_coeff.f0 / BP5_coeff.b0)
+            ψ0 = BP5_coeff.f0 + BP5_coeff.b0 * log(BP5_coeff.V0 * θ0 / BP5_coeff.L)
+            ψ[index] = ψ0
+            θ[index] = θ0
+        end
+    end
+
+
+    for i in 1:fN2
+        for j in 1:fN3
+            index = i + (j - 1) * fN2
+            tau_index = RS_filter_2D_nzind[index]
+            if index in VW_favorable_filter_RS_nzind
+                τ0 = BP5_coeff.σn * RSas[index] * asinh((0.03 / (2 * BP5_coeff.V0) *
+                                                         exp((BP5_coeff.f0 + BP5_coeff.b0 * log(BP5_coeff.V0 / BP5_coeff.Vinit)) /
+                                                             RSas[index])) + η * 0.03)
+                τ[tau_index] = τ0
+                V2_v[index] = 0.03
+                RSLs[index] = 0.13
+            end
+        end
+    end
+
+
+
+    path_time = "$path/$time_string/"
+    try
+        mkdir(path_time)
+    catch
+        # folder already exists
+    end
+
+    tspan = (0, BP5_coeff.sim_years * year_seconds)
+    prob = ODEProblem(odefun_BP5, ψδ, tspan, odeparam) 
+
+    if !isfile(path_time * "restart_value.json")
+        global ctr[] = 1
+        create_text_files_BP5(path_time, station_strings, station_indices, δ, τb, θ, 0)
+        tspan = (0, BP5_coeff.sim_years * year_seconds)
+        prob = ODEProblem(odefun_BP5, ψδ, tspan, odeparam)
+    else
+        contents = read(path_time * "restart_values.json", String)
+        json_data = JSON.parse(contents)
+        ψ_restart = json_data["ψ"]
+        δ_restart = json_data["δ"]
+        t = json_data["t"]
+        V_restart = json_data["V"]
+        tspan = (t, BP5_coeff.sim_years * year_seconds)
+        ψ .= ψ_restart
+        δ .= δ_restart
+        V .= V_restart
+        prob = ODEProblem(odefun_BP5, ψδ, tspan, odeparam)
+    end
+
+
+    # Creating output
+    callback_func = SavingCallback(
+        (ψδ, t, i) -> write_to_file_BP5(path_time, ψδ, t, i, odeparam, station_strings, station_indices), SavedValues(Float64, Float64))
+
+
+
+    function stepcheck(_, odeparam, _)
+        if odeparam.reject_step[1]
+            odeparam.reject_step[1] = false
+            println("reject")
+            return true
+        end
+        return false
+    end
+
+    # abstol = 1e-6, reltol = 1e-8, gamma = 0.8 works
+    # try larger abstol and reltol for better performance[]
+    # sol = solve(prob, Tsit5(); isoutofdomain=stepcheck, dt=0.001, dtmin=1e-8, abstol=1e-6, reltol=1e-8, gamma = 0.8, save_everystep=false,
+    #     callback=callback_func) 
+
+    sol = solve(prob, Tsit5(); isoutofdomain=stepcheck, dt=0.001, dtmin=1e-8, abstol=1e-6, reltol=1e-6, save_everystep=false,
+        callback=callback_func) 
+end
+
+
+main()