further simplified BP1-qd and test code

Thrase · Apr 5, 2024 · 86269c2 · 86269c2
1 parent dfa9ae0
commit 86269c2
Show file tree

Hide file tree

Showing 8 changed files with 213 additions and 66 deletions.
diff --git a/Manifest.toml b/Manifest.toml
@@ -1,8 +1,8 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.10.0"
+julia_version = "1.10.2"
 manifest_format = "2.0"
-project_hash = "e903746fa488e3703286eac6764293ec01fd8c83"
+project_hash = "0b8ccef874ca12289c9bdafef2290b49668af449"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "41c37aa88889c171f1300ceac1313c06e891d245"
@@ -64,6 +64,12 @@ weakdeps = ["SparseArrays"]
 [[deps.Artifacts]]
 uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
 
+[[deps.AxisAlgorithms]]
+deps = ["LinearAlgebra", "Random", "SparseArrays", "WoodburyMatrices"]
+git-tree-sha1 = "01b8ccb13d68535d73d2b0c23e39bd23155fb712"
+uuid = "13072b0f-2c55-5437-9ae7-d433b7a33950"
+version = "1.1.0"
+
 [[deps.BandedMatrices]]
 deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra", "PrecompileTools"]
 git-tree-sha1 = "27baf04c642465b4289179f29bb7127f0673d4f1"
@@ -131,6 +137,16 @@ git-tree-sha1 = "f641eb0a4f00c343bbc32346e1217b86f3ce9dad"
 uuid = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
 version = "0.5.1"
 
+[[deps.ChainRulesCore]]
+deps = ["Compat", "LinearAlgebra"]
+git-tree-sha1 = "575cd02e080939a33b6df6c5853d14924c08e35b"
+uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+version = "1.23.0"
+weakdeps = ["SparseArrays"]
+
+    [deps.ChainRulesCore.extensions]
+    ChainRulesCoreSparseArraysExt = "SparseArrays"
+
 [[deps.CloseOpenIntervals]]
 deps = ["Static", "StaticArrayInterface"]
 git-tree-sha1 = "70232f82ffaab9dc52585e0dd043b5e0c6b714f1"
@@ -195,7 +211,7 @@ weakdeps = ["Dates", "LinearAlgebra"]
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
-version = "1.0.5+1"
+version = "1.1.0+0"
 
 [[deps.ConcreteStructs]]
 git-tree-sha1 = "f749037478283d372048690eb3b5f92a79432b34"
@@ -337,15 +353,12 @@ deps = ["LinearAlgebra", "Statistics", "StatsAPI"]
 git-tree-sha1 = "66c4c81f259586e8f002eacebc177e1fb06363b0"
 uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
 version = "0.10.11"
+weakdeps = ["ChainRulesCore", "SparseArrays"]
 
     [deps.Distances.extensions]
     DistancesChainRulesCoreExt = "ChainRulesCore"
     DistancesSparseArraysExt = "SparseArrays"
 
-    [deps.Distances.weakdeps]
-    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-    SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-
 [[deps.Distributed]]
 deps = ["Random", "Serialization", "Sockets"]
 uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
@@ -656,6 +669,16 @@ version = "2024.0.2+0"
 deps = ["Markdown"]
 uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 
+[[deps.Interpolations]]
+deps = ["Adapt", "AxisAlgorithms", "ChainRulesCore", "LinearAlgebra", "OffsetArrays", "Random", "Ratios", "Requires", "SharedArrays", "SparseArrays", "StaticArrays", "WoodburyMatrices"]
+git-tree-sha1 = "88a101217d7cb38a7b481ccd50d21876e1d1b0e0"
+uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+version = "0.15.1"
+weakdeps = ["Unitful"]
+
+    [deps.Interpolations.extensions]
+    InterpolationsUnitfulExt = "Unitful"
+
 [[deps.IrrationalConstants]]
 git-tree-sha1 = "630b497eafcc20001bba38a4651b327dcfc491d2"
 uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
@@ -937,16 +960,12 @@ deps = ["ArrayInterface", "CPUSummary", "CloseOpenIntervals", "DocStringExtensio
 git-tree-sha1 = "0f5648fbae0d015e3abe5867bca2b362f67a5894"
 uuid = "bdcacae8-1622-11e9-2a5c-532679323890"
 version = "0.12.166"
+weakdeps = ["ChainRulesCore", "ForwardDiff", "SpecialFunctions"]
 
     [deps.LoopVectorization.extensions]
     ForwardDiffExt = ["ChainRulesCore", "ForwardDiff"]
     SpecialFunctionsExt = "SpecialFunctions"
 
-    [deps.LoopVectorization.weakdeps]
-    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-    ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
-    SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
-
 [[deps.MKL_jll]]
 deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl"]
 git-tree-sha1 = "72dc3cf284559eb8f53aa593fe62cb33f83ed0c0"
@@ -1084,7 +1103,7 @@ version = "1.3.5+1"
 [[deps.OpenBLAS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
 uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
-version = "0.3.23+2"
+version = "0.3.23+4"
 
 [[deps.OpenLibm_jll]]
 deps = ["Artifacts", "Libdl"]
@@ -1293,6 +1312,16 @@ git-tree-sha1 = "043da614cc7e95c703498a491e2c21f58a2b8111"
 uuid = "e6cf234a-135c-5ec9-84dd-332b85af5143"
 version = "1.5.3"
 
+[[deps.Ratios]]
+deps = ["Requires"]
+git-tree-sha1 = "1342a47bf3260ee108163042310d26f2be5ec90b"
+uuid = "c84ed2f1-dad5-54f0-aa8e-dbefe2724439"
+version = "0.4.5"
+weakdeps = ["FixedPointNumbers"]
+
+    [deps.Ratios.extensions]
+    RatiosFixedPointNumbersExt = "FixedPointNumbers"
+
 [[deps.RecipesBase]]
 deps = ["PrecompileTools"]
 git-tree-sha1 = "5c3d09cc4f31f5fc6af001c250bf1278733100ff"
@@ -1514,13 +1543,11 @@ deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_j
 git-tree-sha1 = "e2cfc4012a19088254b3950b85c3c1d8882d864d"
 uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
 version = "2.3.1"
+weakdeps = ["ChainRulesCore"]
 
     [deps.SpecialFunctions.extensions]
     SpecialFunctionsChainRulesCoreExt = "ChainRulesCore"
 
-    [deps.SpecialFunctions.weakdeps]
-    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-
 [[deps.Static]]
 deps = ["IfElse"]
 git-tree-sha1 = "f295e0a1da4ca425659c57441bcb59abb035a4bc"
@@ -1543,15 +1570,12 @@ deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"]
 git-tree-sha1 = "f68dd04d131d9a8a8eb836173ee8f105c360b0c5"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
 version = "1.9.1"
+weakdeps = ["ChainRulesCore", "Statistics"]
 
     [deps.StaticArrays.extensions]
     StaticArraysChainRulesCoreExt = "ChainRulesCore"
     StaticArraysStatisticsExt = "Statistics"
 
-    [deps.StaticArrays.weakdeps]
-    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-    Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
-
 [[deps.StaticArraysCore]]
 git-tree-sha1 = "36b3d696ce6366023a0ea192b4cd442268995a0d"
 uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
@@ -1780,6 +1804,12 @@ git-tree-sha1 = "93f43ab61b16ddfb2fd3bb13b3ce241cafb0e6c9"
 uuid = "2381bf8a-dfd0-557d-9999-79630e7b1b91"
 version = "1.31.0+0"
 
+[[deps.WoodburyMatrices]]
+deps = ["LinearAlgebra", "SparseArrays"]
+git-tree-sha1 = "c1a7aa6219628fcd757dede0ca95e245c5cd9511"
+uuid = "efce3f68-66dc-5838-9240-27a6d6f5f9b6"
+version = "1.0.0"
+
 [[deps.XML2_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Zlib_jll"]
 git-tree-sha1 = "801cbe47eae69adc50f36c3caec4758d2650741b"

diff --git a/Project.toml b/Project.toml
@@ -6,6 +6,7 @@ version = "0.1.0"
 [deps]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
@@ -19,4 +20,4 @@ julia = "1"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test"]
+test = ["Test"]
diff --git a/examples/bp1-qd.dat b/examples/bp1-qd.dat
@@ -11,20 +11,20 @@ stride_space = 1
 stride_time = 5
 #
 # computational domain size
-# x-domain is (x1, x2) (km)
-x1 = 0
-x2 = 80
-# z-domain is (z1, z2) (km)
-z1 = 0
-z2 = 80
+# x-domain is (Lx1, Lx2) (km)
+Lx1 = 0
+Lx2 = 80
+# z-domain is (Lz1, Lz2) (km)
+Lz1 = 0
+Lz2 = 80
 #
 # number of grid points in each dimension 
-# make sure to resolve h* and l_b - should be written as integer.
+# make sure to resolve h* and l_b; should be written as integer.
 Nx = 128
 Nz = 128
 #
 # how many years to simulate
-sim_years = 15
+sim_years = 500
 #
 # loading rate
 Vp = 1e-9

diff --git a/examples/stripped_qd_driver.jl b/examples/stripped_qd_driver.jl
@@ -35,29 +35,30 @@ function main()
 
     Lx = xc[2]
     Lz = zc[2]
+
     # create operators
     (M̃, F, τ, H̃, HfI_FT) = get_operators(SBPp, Nx, Nz, μ; xc = xc, zc = zc)
-    #(M̃, F, τ, H̃, HfI_FT) = get_operators(SBPp, Nr, Ns, μ)
 
-    # factor with Cholesky
-    M = cholesky(Symmetric(M̃))
+    # factor matrix with Cholesky
+    A = cholesky(Symmetric(M̃))
 
-    # initialize vector g that stores boundary data
+    # initialize time and vector b that stores boundary data (linear system will be Au = b)
     t = 0
-    g = zeros((Nx+1) * (Nz+1))
+    b = zeros((Nx+1) * (Nz+1))
 
+    # initial slip vector
     δ = zeros(Nz+1)
 
-    bc_Dirichlet = (lf, x, z) -> (2-lf) .* (0 * x .+ 0 .* z) + (lf-1) .* (0 .* x .+ 0 .* z)
-    bc_Neumann   = (lf, x, z, nx, ny) -> zeros(size(x))
-
-    bdry_vec_mod!(g, F, τ, x, z, bc_Dirichlet, bc_Neumann, Lx, Lz)
+    # fill in initial boundary data into b
+    bdry_vec_strip!(b, F, τ, x, z, δ ./ 2, (t .* Vp./2)*ones(size(z)), zeros(size(x)), Lx, Lz)
 
-    u = M \ g
+    u = A \ b # solve linear system with a backsolve to obtain initial displacement u(x, z, 0)
 
+    # Find z-index δNp corresponding to base of rate-and-state friction zone RSWf
     (mm, δNp) = findmin(abs.(RSWf .- z))
     @assert z[δNp] ≈ RSWf
 
+    # initialize change in shear stress due to quasi-static deformation
     Δτ = zeros(Nz+1)
 
     # Assemble fault variables/data
@@ -67,26 +68,29 @@ function main()
             min(1, max(0, (RSH1 - z[n])/(RSH1 - RSH2)))
     end
 
-
-    τz0 = σn * RSamax * asinh(RSVinit / (2 * RSV0) *
+    # Set pre-stress according to benchmark
+    τ0 = σn * RSamax * asinh(RSVinit / (2 * RSV0) *
                                     exp.((RSf0 + RSb * log.(RSV0 / RSVinit)) /
                                         RSamax)) + η * RSVinit
 
 
+    # Set initial state variable according to benchmark
     θ = (RSDc ./ RSV0) .* exp.((RSa ./ RSb) .* log.((2 .* RSV0 ./ RSVinit) .*
-        sinh.((τz0 .- η .* RSVinit) ./ (RSa .* σn))) .- RSf0 ./ RSb)
-
+        sinh.((τ0 .- η .* RSVinit) ./ (RSa .* σn))) .- RSf0 ./ RSb)
 
-    ψ0 = RSf0 .+ RSb .* log.(RSV0 .* θ ./ RSDc)
+    # Initialize psi version of state variable
+    ψ = RSf0 .+ RSb .* log.(RSV0 .* θ ./ RSDc)
 
 
-    ψδ = zeros(δNp + Nz + 1)  #because length(ψ) = δNp,  length(δ) = N+1
-    ψδ[1:δNp] .= ψ0
-
+    # Set initial condition for index 1 DAE - this is a stacked vector of psi, followed by slip
+    ψδ = zeros(δNp + Nz + 1)  #because length(ψ) = δNp,  length(δ) = Nz+1
+    ψδ[1:δNp] .= ψ
+    ψδ[δNp+1:δNp + Nz + 1] .= δ
 
+    # Set fault station locations (depths) specified in benchmark
     stations = [0, 2.5, 5, 7.5, 10, 12.5, 15, 17.5, 20, 25, 30, 35] # km
 
-    # Next part of code related to code output
+    # Function that finds the depth-index corresponding to a station location
     function find_station_index(stations, grid_points)
         numstations = length(stations)
         station_ind = zeros(numstations)
@@ -98,28 +102,28 @@ function main()
       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. 
+    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
 
 
-    # Fault locations to record slip evolution output:
+    # Benchmark description also requests slip evolution output recorded at fault locations every ~stride_space spatial nodes:
     flt_loc = z[1:stride_space:δNp]  # physical stations (units of km)
     flt_loc_indices = find_station_index(flt_loc, z)
 
-    # set up parameters sent to right hand side
+    # set up parameters sent to the right hand side of the DAE:
     odeparam = (reject_step = [false],
                 Vp=Vp,
-                M = M,
+                A = A,
                 u=u,
                 Δτ = Δτ,
-                τf = τz0*ones(Nz+1),
-                g = g,
+                τf = τ0*ones(Nz+1),
+                b = b,
                 μshear=μshear,
                 RSa=RSa,
                 RSb=RSb,
                 σn=σn,
                 η=η,
                 RSV0=RSV0,
-                τz0=τz0,
+                τ0=τ0,
                 RSDc=RSDc,
                 RSf0=RSf0,
                 δNp = δNp,
@@ -134,16 +138,21 @@ function main()
                 save_stride_fields = stride_time # save every save_stride_fields time steps
                 )
 
-    #dψV = zeros(δNp + Nz + 1)
+    # Set time span over which to solve:
     tspan = (0, sim_years * year_seconds)
-    prob = ODEProblem(odefun, ψδ, tspan, odeparam)
+
+
+    # Set up ODE problem corresponding to DAE
+    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))
 
     # Make text files to store on-fault time series and slip data,
     # Also initialize with initial data:
-    create_text_files(pth, flt_loc, flt_loc_indices, stations, station_strings, station_indices, 0, RSVinit, δ, τz0, θ)
+    create_text_files(pth, flt_loc, flt_loc_indices, stations, station_strings, station_indices, 0, RSVinit, δ, τ0, θ)
 
+    # Solve DAE using Tsit5()
     sol = solve(prob, Tsit5(); dt=0.2,
              abstol = 1e-5, reltol = 1e-5, save_everystep=true, gamma = 0.2,
              internalnorm=(x, _)->norm(x, Inf), callback=cb_fun)