diff --git a/Manifest.toml b/Manifest.toml index bc332c1..78d6aa3 100644 --- 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 index bca6a43..66500f2 100644 --- 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"] \ No newline at end of file +test = ["Test"] diff --git a/examples/bp1-qd.dat b/examples/bp1-qd.dat index 5cbf190..fd95583 100644 --- 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 index 2d507f0..e708e4e 100644 --- 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) diff --git a/src/odefun.jl b/src/odefun.jl index 2ed06af..84248d8 100644 --- a/src/odefun.jl +++ b/src/odefun.jl @@ -7,6 +7,86 @@ using Printf using DelimitedFiles +function odefun_stripped(dψV, ψδ, p, t) + + Vp = p.Vp + A = p.A + u = p.u + Δτ = p.Δτ + τf = p.τf + b = p.b + μshear = p.μshear + RSa = p.RSa + RSb = p.RSb + σn = p.σn + η = p.η + RSV0 = p.RSV0 + τ0 = p.τ0 + RSDc = p.RSDc + RSf0 = p.RSf0 + δNp = p.δNp + N = p.N + F = p.F + τ = p.τ + x = p.x + z = p.z + HfI_FT = p.HfI_FT + Lx = p.Lx + Lz = p.Lz + + current_time = t ./ 31556926 + print("TIME [YRS] = $(current_time).\n") + + ψ = @view ψδ[ (1:δNp) ] + δ = ψδ[ δNp .+ (1:N+1) ] + + + bdry_vec_strip!(b, F, τ, x, z, δ ./ 2, (t .* Vp./2)*ones(size(z)), zeros(size(x)), Lx, Lz) + + # solve for displacements everywhere in domain + u[:] = A \ b + + # set up rates of change for state and slip + dψ = @view dψV[ (1:δNp) ] + V = @view dψV[ δNp .+ (1:N+1)] + + dψ .= 0 # initialize values to 0 + V .= 0 # initialize values to 0 + + # Update the fault data + Δτ .= 0 + Δτ .= -computetraction_stripped(HfI_FT, τ, u, δ) + τf .= τ0 .+ Δτ + + # Do safe-guarded Newton at every node in rate-and-state friction zone in order to solve for slip rate V. + for n = 1:δNp + ψn = ψ[n] + an = RSa[n] + + τn = Δτ[n] + τ0 + + VR = abs(τn / η) + VL = -VR + Vn = V[n] + obj_rs(V) = rateandstate(V, ψn, σn, τn, η, an, RSV0) + (Vn, _, iter) = newtbndv(obj_rs, VL, VR, Vn; ftol = 1e-9, + atolx = 1e-9, rtolx = 1e-9) + + + + V[n] = Vn # update slip rate + + dψ[n] = (RSb * RSV0 / RSDc) * (exp((RSf0 - ψn) / RSb) - abs(Vn) / RSV0) # update aging law + + end + + V[δNp+1:N+1] .= Vp # for all points below rate-and-state region, set slip rate to plate rate Vp. + + + nothing +end + + function odefun(dψV, ψδ, p, t) @@ -60,7 +140,7 @@ function odefun(dψV, ψδ, p, t) Δτ .= 0 lf1 = 1 # fault is at face 1 - Δτ .= -computetraction_stripped(HfI_FT, τ, lf1, u, δ) + Δτ .= -computetraction_stripped(HfI_FT, τ, u, δ) τf .= τz0 .+ Δτ @@ -93,4 +173,4 @@ end -export odefun \ No newline at end of file +export odefun, odefun_stripped \ No newline at end of file diff --git a/src/ops_stripped.jl b/src/ops_stripped.jl index b481ba7..b370ed0 100644 --- a/src/ops_stripped.jl +++ b/src/ops_stripped.jl @@ -160,6 +160,32 @@ function get_operators(p, Nx, Nz, μ; xc = (-1, 1), zc = (-1, 1)) return (M̃ , F, τ, H̃, HfI_FT) end +function bdry_vec_strip!(g, F, τ, x, z, slip_data, remote_data, free_surface_data, Lx, Lz) + + + g[:] .= 0 + + # fault (Dirichlet): + vf = slip_data + g[:] -= F[1] * vf + + # FACE 2 (Dirichlet): + vf = remote_data + g[:] -= F[2] * vf + + # FACE 3 (Neumann): + gN = free_surface_data + vf = gN ./ diag(τ[3]) + g[:] -= F[3] * vf + + # FACE 4 (Neumann): + gN = free_surface_data + vf = gN ./ diag(τ[4]) + g[:] -= F[4] * vf + + +end + function bdry_vec_mod!(g, F, τ, x, z, bc_Dirichlet, bc_Neumann, Lx, Lz) @@ -186,9 +212,9 @@ function bdry_vec_mod!(g, F, τ, x, z, bc_Dirichlet, bc_Neumann, Lx, Lz) end -function computetraction_stripped(HfI_FT, τ, lf, u, δ) - HfI_FT = HfI_FT[lf] - τf = τ[lf] +function computetraction_stripped(HfI_FT, τ, u, δ) + HfI_FT = HfI_FT[1] + τf = τ[1] return (HfI_FT * u + τf * (δ .- δ / 2)) end @@ -240,5 +266,5 @@ function computetraction_stripped(HfI_FT, τ, lf, u, δ) return (x, f, -maxiter) end - export get_operators, bdry_vec_mod!, computetraction_stripped + export get_operators, bdry_vec_mod!, bdry_vec_strip!, computetraction_stripped export rateandstate, newtbndv \ No newline at end of file diff --git a/src/utils.jl b/src/utils.jl index 1babda6..4d1b0cc 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -3,6 +3,7 @@ using SparseArrays using LinearAlgebra using DelimitedFiles using DifferentialEquations +using Interpolations function create_text_files(pth, flt_loc, flt_loc_indices, stations, station_strings, station_indices, t, RSVinit, δ, τz0, θ) diff --git a/test/runtests.jl b/test/runtests.jl index acb4f15..c67fb98 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,8 +4,8 @@ using Printf @testset "Thrase.jl" begin try - localARGS = ["../examples/test.dat"] - include("../examples/stripped_qd_driver.jl"); + localARGS = ["examples/test.dat"] + include("examples/stripped_qd_driver.jl"); catch print("cannot run bp1-qd") end