Patch 1

.travis.yml

@@ -4,7 +4,7 @@ os:
   - linux
   - osx
 julia:
-  - 1.0
+  - 1.1
   - nightly
 matrix:
   allow_failures:
@@ -22,7 +22,7 @@ after_success:
 jobs:
   include:
     - stage: "Documentation"
-      julia: 1.0
+      julia: 1.1
       os: linux
       install:
         # mkdocs-material isn't available at the default py2 on travis-ci for now

Project.toml

@@ -1,17 +1,19 @@
 name = "JuEQ"
 uuid = "70fec03b-432a-5687-8a7d-e45297f2e4e6"
 authors = ["Pengcheng Shi"]
-version = "0.1.6"
+version = "0.2.0"
 
 [deps]
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SharedArrays = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
 SimpleTraits = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
+Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 
 [extras]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"

docs/src/examples/bp1.jl

@@ -1,13 +1,12 @@
 # !!! note
 #     This example is from [Benchmark Problem 1](https://www.scec.org/publication/8214) (hence referred as BP1).
 #
-# ### Define Parameters
 
 # First, we load the package
 using JuEQ
 using Plots
 
-# Instead of using SI unit, we refactor ours into the follow:
+# The prerequisite parameters in this benchmark are list below:
 
 ms2mmyr = 365 * 86400 * 1e3 # convert velocity from m/s to mm/yr
 ρ = 2670.0 # density [kg/m³]
@@ -25,84 +24,71 @@ H = 15.0 # vw zone [km]
 h = 3.0 # vw-vs changing zone [km]
 Wf = 40.0 # fault depth [km]
 Δz = 100.0e-3 # grid size interval [km]
-tf = 400.0; # simulation time [yr]
+tf = 200.0; # simulation time [yr]
+nothing
 
 # !!! warning
-#     Make sure your units are consistent across the whole variable space. Pontenial imporvement may
-#     incoporate [Unitful.jl](https://github.com/ajkeller34/Unitful.jl) package.
+#     Make sure your units are consistent across the whole variable space.
 
 # Then we arrive at some parameters that are implicit by above:
 
 μ = vs^2 * ρ / 1e5 / 1e6 # shear modulus [bar·km/mm]
 λ = μ # poisson material
 η = μ / 2(vs * 1e-3 * 365 * 86400)
 ngrid = round(Int, Wf / Δz); # number of grids
+nothing
 
-# Now, we start to construct our model using parameters above. First, we create a 'fault' by specifying fault type and depth:
+# First, create a **fault space** by specifying fault type, depth
+# and with the desired discretization interval.
 # !!! tip
 #     Here, we do not need to provide `dip` for strike-slip fault as it automatically choose `90`. See [`fault`](@ref).
 
-# ### Construct Model
+fa = fault(Val(:CSFS), STRIKING(), 40.0, Δz)
+nothing
 
-fa = fault(StrikeSlipFault, Wf);
+# Then, provide the material properties w.r.t. our 'fault space'. The properties contain two part. First one is
+# **Fault Properties** which includes *elastic modulus*, *radiation damping term*, *plate rate*, *reference friction* and *reference velocity*.
 
-# Next, we generate the grid regarding the `fault` we just created by giving number of grids:
-# !!! note
-#     This package use `ξ` for denoting downdip coordinate and `x` for along-strike one. See [`discretize`](@ref).
-
-gd = discretize(fa; nξ=ngrid);
+faprop = init_fault_prop(λ, μ, η, vpl, f0, v0)
+nothing
 
-# Next, we construct the required frictional parameter profile:
+# Second one is **Rate-and-State Friction Properties** which includes *a*, *b*, *L*, *σ*. The default friction formula is **Regularized Form**
+# which is stored in the field `flf` and is denoted as `RForm()`. The default state evolution law is `DieterichStateLaw()`.
 
-z = -gd.ξ
-az = fill(a0, size(z))
-az[z .≥ (H + h)] .= amax
-az[H .< z .< H + h] = a0 .+ (amax - a0) / (h / Δz) * collect(1: Int(h / Δz));
-
-# Then, we provide the required initial condition satisfying uniform slip distribution over the depth:
+frprop = init_friction_prop(fa)
+fill!(frprop.a, a0)
+frprop.a[-fa.mesh.z .≥ (H + h)] .= amax
+frprop.a[H .< -fa.mesh.z .< H + h] .= a0 .+ (amax - a0) / (h / Δz) * collect(1: Int(h / Δz))
+fill!(frprop.b, b0)
+fill!(frprop.L, L)
+fill!(frprop.σ, σ)
+nothing
 
+# Next, construct the initial condition and ODE problem using Okada's Green's function.
 τ0 = σ * amax * asinh(vinit / 2v0 * exp((f0 + b0 * log(v0 / vinit)) / amax)) + η * vinit
-τz = fill(τ0, size(z))
-θz = @. L / v0 * exp(az / b0 * log(2v0 / vinit * sinh((τz - η * vinit) / az / σ)) - f0 / b0)
-vz = fill(vinit, size(z))
-u0 = hcat(vz, θz);
-
-# Let's simulate only the first 200 years:
+τz = fill(τ0, size(fa.mesh.z))
+θz = @. L / v0 * exp(frprop.a / b0 * log(2v0 / vinit * sinh((τz - η * vinit) / frprop.a / σ)) - f0 / b0)
+vz = fill(vinit, size(fa.mesh.ξ))
+δz = fill(0.0, size(fa.mesh.ξ))
+u0 = hcat(vz, θz, δz);
+prob = assemble(Val(:okada), fa, faprop, frprop, u0, (0.0, tf))
+nothing
 
-tspan = (0., 200.);
+# Check our depth profile now.
 
-# Finally, we provide the material properties w.r.t. our 'fault', 'grid' as well as other necessary parameters predefined
-# using the same grid size & dimension:
-
-mp = properties(;fault=fa, grid=gd, parameters=[:a=>az, :b=>b0, :L=>L, :σ=>σ, :η=>η, :k=>[:λ=>λ, :μ=>μ], :vpl=>vpl, :f0=>f0, :v0=>v0]);
-
-# !!! tip
-#     Check [`properties`](@ref) for extended options.
+plot([frprop.a, frprop.b], fa.mesh.z, label=["a", "b"], yflip=true, ylabel="Depth (km)")
 
-# Check our profile now:
+# Afterwards, solve ODE thanks to [DifferentialEquations.jl](https://github.com/JuliaDiffEq/DifferentialEquations.jl)
 
-plot([mp.a, mp.b], z, label=["a", "b"], yflip=true, ylabel="Depth (km)")
-
-# We then contruct the `ODEProblem` as following by stating which state evolution law to use and frcitonal law form,
-# plus initial condition and simulation time:
-
-prob = EarthquakeCycleProblem(gd, mp, u0, tspan; se=DieterichStateLaw(), fform=RForm());
-
-# ### Solve Model
-# We then solve the ODEs:
-
-sol = solve(prob, Tsit5(), reltol=1e-6, abstol=1e-6);
-
-# !!! tip
-#     For details of solving options, see [here](http://docs.juliadiffeq.org/latest/basics/common_solver_opts.html).
+sol = solve(prob, Tsit5(), reltol=1e-6, abstol=1e-6)
+nothing
 
 # !!! tip
-#     Raise the accuracy option if you get instability when solving these ODEs.
+#     Raise the accuracy option or switch to other algorithms if you get instability when solving these ODEs.
 
-# ### Results
-# The first event happens at around 196 year:
+# Finally, check the results. The first event happens at around 196 year:
 
-maxv = max_velocity(sol)
+maxv = JuEQ.max_velocity(sol);
 plot(sol.t, log10.(maxv / ms2mmyr), xlabel="Time (year)", ylabel="Max Velocity (log10 (m/s))", xlims=(190, 200), label="")
 
 # !!! note

docs/src/examples/otfsync.jl

@@ -9,8 +9,7 @@
 #     using JuEQ
 #     ```
 
-# ### Define Parameters
-# First, we load the package and define some basic parameters:
+# First, list all the essential parameters:
 
 using JuEQ
 using Plots
@@ -21,94 +20,88 @@ cs = 3044.0 # m/s
 vpl = 100.0 # mm/yr
 v0 = 3.2e4 # mm/yr
 f0 = 0.6;
-
-# Then we come to parameters implicit by above:
-
 μ = 0.3 # Bar·km/mm
 λ = μ # poisson material
 α = (λ + μ) / (λ + 2μ)
-η = μ / 2(cs * 1e-3 * 365 * 86400); # Bar·yr/mm
+η = μ / 2(cs * 1e-3 * 365 * 86400) # Bar·yr/mm
+nothing
 
-# Create a fault:
+# First, create a fault space.
 
-fa = fault(StrikeSlipFault, (80., 10.));
+fa = fault(Val(:CSFS), STRIKING(), (80., 10.), (0.5, 0.5))
+nothing
 
-# Generate grids:
+# Next, establish frictional and fault space parameters:
 
-gd = discretize(fa; nx=160, nξ=20, buffer_ratio=1);
+frprop = init_friction_prop(fa)
+faprop = init_fault_prop(λ, μ, η, vpl, f0, v0)
 
-# !!! tip
-#     It is recommended (from Yajing Liu's personal communication) to add buffer zones adjacent the horizontal edges
-#     to immitate *zero* dislocation at the ridge region.
-#     Basically, it affects how the stiffness tensor are periodically summed. To what extent it alters the results remains further testing.
+fill!(frprop.a, 0.015)
+fill!(frprop.b, 0.0115)
+fill!(frprop.L, 12.0)
 
-#     Under the hood, it shall impose buffer areas on both sides of along-strike, each of which has a length of `bufferratio/2*fa[:x]`.
-#     Thus, the stiffness contributions falling into those buffer zone shall be neglected, which is equivalent to impose zero-slip correspondingly.
+left_patch = @. -25. ≤ fa.mesh.x ≤ -5.
+right_patch = @. 5. ≤ fa.mesh.x ≤ 25.
+vert_patch = @. -6. ≤ fa.mesh.z ≤ -1
 
-# Time for us to establish frictional parameters profile:
+frprop.b[xor.(left_patch, right_patch), vert_patch] .= 0.0185
 
-a = 0.015 .* ones(gd.nx, gd.nξ)
-b = 0.0115 .* ones(gd.nx, gd.nξ)
-left_patch = @. -25. ≤ gd.x ≤ -5.
-right_patch = @. 5. ≤ gd.x ≤ 25.
-vert_patch = @. -6. ≤ gd.z ≤ -1.
-b[xor.(left_patch, right_patch), vert_patch] .= 0.0185
-amb = a - b
 σmax = 500.
-σ = [min(σmax, 15. + 180. * z) for z in -gd.z]
-σ = Matrix(repeat(σ, 1, gd.nx)')
-L = 12.;
+σ = [min(σmax, 15. + 180. * z) for z in -fa.mesh.ξ]
+σ = Matrix(repeat(σ, 1, fa.mesh.nx)')
+L = 12.
 
-# Check our profile:
+frprop.σ .= σ
+nothing
 
-p1 = heatmap(amb',
-    xticks=(0: 10/gd.Δx: gd.nx, -fa.span[1]/2: 10: fa.span[1]/2),
-    yticks=(0: 5/gd.Δξ: gd.nξ, 0: -5: -fa.span[2]),
-    yflip=true, color=:isolum, aspect_ratio=2, title="a-b"
+# Make sure our profile match our expectation:
+
+p1 = plot((frprop.a .- frprop.b)', seriestype=:heatmap,
+    xticks=(collect(1: 40: fa.mesh.nx+1), [-40, -20, 0, 20, 40]),
+    yticks=(collect(1: 5: fa.mesh.nξ+1), [0, 5, 10, 15, 20]),
+    yflip=true, color=:isolum, aspect_ratio=2, title="a-b",
     );
 
 p2 = heatmap(σ',
-    xticks=(0: 10/gd.Δx: gd.nx, -fa.span[1]/2: 10: fa.span[1]/2),
-    yticks=(0: 5/gd.Δξ: gd.nξ, 0: -5: -fa.span[2]),
+    xticks=(collect(1: 40: fa.mesh.nx+1), [-40, -20, 0, 20, 40]),
+    yticks=(collect(1: 5: fa.mesh.nξ+1), [0, 5, 10, 15, 20]),
     yflip=true, color=:isolum, aspect_ratio=2, title="\\sigma"
     );
 
 plot(p1, p2, layout=(2, 1))
 
-# ### Construct Model
-# Construct our material property profile:
-
-mp = properties(fa, gd, [:a=>a, :b=>b, :L=>L, :σ=>σ, :η=>η, :k=>[:λ=>λ, :μ=>μ], :vpl=>vpl, :f0=>f0, :v0=>v0]);
+# Then, provide the initial condition and assemble the ODEs:
 
-# Provide the initial condition:
-
-vinit = vpl .* ones(gd.nx, gd.nξ)
+vinit = vpl .* ones(fa.mesh.nx, fa.mesh.nξ)
 θ0 = L ./ vinit ./ 1.1
-u0 = cat(vinit, θ0, dims=3);
+u0 = cat(vinit, θ0, zeros(Float64, fa.mesh.nx, fa.mesh.nξ), dims=3)
+prob = assemble(Val(:okada), fa, faprop, frprop, u0, (0., 18.), buffer_ratio=1)
+nothing
 
-# Get our ODEs problem:
+# !!! tip
+#     It is recommended (from Yajing Liu's personal communication) to add buffer zones adjacent the horizontal edges
+#     to immitate *zero* dislocation at the ridge region.
+#     Basically, it affects how the stiffness tensor are periodically summed. To what extent it alters the results remains further testing.
 
-prob = EarthquakeCycleProblem(gd, mp, u0, (0., 18.); se=DieterichStateLaw(), fform=CForm());
+#     Under the hood, it shall impose buffer areas on both sides of along-strike, each of which has a length of `bufferratio/2*fa[:x]`.
+#     Thus, the stiffness contributions falling into those buffer zone shall be neglected, which is equivalent to impose zero-slip correspondingly.
 
-# ### Solve Model
-# Solve the model:
+# Afterwards, solve ODEs problem:
 
-sol = solve(prob, Tsit5(), reltol=1e-6, abstol=1e-6);
+sol = solve(prob, Tsit5(), reltol=1e-6, abstol=1e-6)
+nothing
 
-# ### Sanity Check of Results
-# Take a look at the max velocity:
+# Last, take a look at the max velocity time series:
 
-maxv = max_velocity(sol)
+maxv = JuEQ.max_velocity(sol)
 plot(sol.t, log10.(maxv / ms2mmyr), xlabel="Time (year)", ylabel="Max Velocity (log10 (m/s))", label="")
 
-# View some snapshots to see the rupture (quasi-dynamic) patterns:
+# And view some snapshots of ruptures (quasi-dynamic) patterns:
 
 ind = argmax(maxv)
 myplot = (ind) -> heatmap(log10.(sol.u[ind][:,:,1]./ms2mmyr)',
-    xticks=(0: 10/gd.Δx: gd.nx, -fa.span[1]/2: 10: fa.span[1]/2),
-    yticks=(0: 5/gd.Δξ: gd.nξ, 0: -5: -fa.span[2]),
+    xticks=(collect(1: 40: fa.mesh.nx+1), [-40, -20, 0, 20, 40]),
+    yticks=(collect(1: 5: fa.mesh.nξ+1), [0, 5, 10, 15, 20]),
     yflip=true, color=:isolum, aspect_ratio=2, title="t = $(sol.t[ind])")
-
-snaps = [myplot(i) for i in ind-700: 200: ind+500]
-
-plot(snaps..., layout=(length(snaps), 1), size=(600, 1800))
+snaps = myplot(ind)
+plot(snaps)

docs/src/index.md

@@ -18,7 +18,7 @@ Features to be implemented:
 Get the latest version with Julia's package manager:
 
 ```julia
-] add https://github.com/shipengcheng1230/JuEQ.jl
+(v1.1) pkg> add https://github.com/shipengcheng1230/JuEQ.jl
 ```
 
 To load the package:

src/JuEQ.jl

@@ -3,6 +3,8 @@ module JuEQ
 using Reexport
 
 @reexport using DifferentialEquations
+@reexport using Sundials
+@reexport using HDF5
 
 using SimpleTraits
 using Parameters
@@ -14,27 +16,18 @@ using Base.Threads
 using LinearAlgebra
 using SharedArrays
 
-include("rsf.jl")
-include("fault.jl")
-include("bem.jl")
-include("constructor.jl")
-include("utils.jl")
-include("dc3d.jl")
 include("config.jl")
+include("rsf.jl")
+include("faulttype.jl")
+include("geometry.jl")
 
-export NormalFault, ThrustFault, StrikeSlipFault
-export PlaneFaultDomain
-export fault
-
-export DieterichStateLaw, RuinaStateLaw, PrzStateLaw
-export CForm, RForm
-export friction
-export MaterialProperties
-
-export discretize, properties, stiffness_tensor
-export EarthquakeCycleProblem
-export dc3d_okada
+include("faultspace.jl")
+include("greensfun.jl")
+include("gfoperator.jl")
+include("properties.jl")
+include("assemble.jl")
 
-export max_velocity, moment_magnitude, DECallbackSaveToFile
+const UTILSDIR = abspath(joinpath(@__DIR__, "utils"))
+map(x -> include(joinpath(UTILSDIR, x)), filter!(x -> endswith(x, ".jl"), readdir(UTILSDIR)))
 
-end
+end # module