commit f65899ee85de5ed0234d7b2164c22ff80b25a978
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Thu, 4 Apr 2019 15:55:32 +0200
Implement correct jacobi iteration call
Diffstat:
| A | 1d_fd_simple_shear.jl | | | 205 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
1 file changed, 205 insertions(+), 0 deletions(-)
diff --git a/1d_fd_simple_shear.jl b/1d_fd_simple_shear.jl
@@ -0,0 +1,205 @@
+#!/usr/bin/env julia
+
+ENV["MPLBACKEND"] = "Agg"
+import PyPlot
+
+# Simulation settings
+
+## Gravitational acceleration magnitude
+const G = 9.81
+
+## Wall parameters
+const P_wall = 10.0e3 # normal stress [Pa]
+# const v_x_wall = 1e-0 # wall velocity along x [m/s]
+const v_x_bot = 0.0 # bottom velocity along x [m/s]
+const μ_wall = 0.382 # stress ratio at top wall
+
+## Nodes along z
+const nz = 10
+
+## Material properties
+const A = 0.48 # nonlocal amplitude [-]
+const b = 0.9377 # rate dependence beyond yield [-]
+const μ_s = 0.3819 # static yield friction coefficient [-]
+const ϕ = 0.38 # porosity [-]
+const d = 1e-3 # representative grain size [m]
+const ρ_s = 2.485e3 # grain material density [kg/m^3]
+
+## Domain size
+const origo_z = 0.0
+const L_z = 20.0*d
+
+## Spatial coordiantes
+const z = collect(range(origo_z, L_z, length=nz))
+const Δz = z[2] - z[1]
+
+## Allocate other arrays
+μ = zero(z) # local stress ratio
+p = zero(z) # local pressure
+v_x = zero(z) # local shear velocity
+g_ghost = zeros(size(z)[1]+2) # local fluidity with ghost nodes
+
+
+# Function definitions
+
+## Shear plastic strain rate (eq 2)
+γ_dot_p(g, μ) = g.*μ
+
+## Stress ratio at the upper wall
+#μ_wall(τ_wall, P_wall) = τ_wall.*P_wall
+
+## Normal stress
+p_lithostatic(P_wall, z) = P_wall .+ (1 - ϕ).*ρ_s.*G.*(L_z .- z)
+
+## local cooperativity length
+ξ(μ) = A*d./sqrt.(abs.(μ .- μ_s))
+
+## Local fluidity
+function g_local(p, μ)
+ if μ <= μ_s
+ return 0.0
+ else
+ return sqrt(p./ρ_s.*d.^2.0) .* (μ .- μ_s)./(b.*μ)
+ end
+end
+
+## Update ghost nodes for g from current values
+## BC: neumann (dg/dx = 0)
+function g_set_bc_neumann(g_ghost)
+ g_ghost[1] = g_ghost[2]
+ g_ghost[end] = g_ghost[end-1]
+end
+
+## Compute shear stress from velocity gradient using finite differences
+function shear_stress(v, Δz)
+ τ = zero(v)
+
+ # compute inner values with central differences
+ for i=2:length(v)-1
+ τ[i] = (v[i+1] - v[i-1])/(2.0*Δz)
+ end
+
+ # use forward/backward finite differences at edges
+ τ[1] = (v[2] - v[1])/Δz
+ τ[end] = (v[end] - v[end-1])/Δz
+
+ return τ
+end
+
+friction(τ, p) = τ./(p .+ 1e-16)
+
+## Fluidity
+#fluidity(p, μ, g) = 1.0./(ξ.(μ).^2.0) .* (g .- g_local.(p, μ))
+
+# ## set up the forcing function for the Poisson equation. The forcing is the
+# ## right-hand side of the fluidity equation
+# function forcing(p, μ, g_ghost)
+# return fluidity.(μ, g_ghost[2:end-1], p)
+# end
+
+## 1d Laplacian operator (2nd order central finite difference)
+#laplacian(phi, i, dx) = (phi[i-1] - 2.0*phi[i] + phi[i+1])/dx
+
+## A single iteration for solving a Poisson equation Laplace(phi) = f on a
+## Cartesian grid. The function returns the normalized residual value
+# TODO f is not used?
+function poisson_solver_1d_iteration(g_in, g_out, μ, p, i, Δz)
+
+ coorp_term = Δz^2.0/(2.0*ξ(μ[i])^2.0)
+ g_out[i] = 1.0/(1.0 + coorp_term) * (coorp_term*g_local(p[i], μ[i])
+ + g_in[i+1]/2.0 + g_in[i-1]/2.0)
+
+ return (g_out[i] - g_in[i])^2.0 / (g_out[i]^2.0 + 1e-16)
+end
+
+## Iteratively solve the system laplace(phi) = f
+function implicit_1d_jacobian_poisson_solver(g, p, μ, Δz,
+ rel_tol=1e-4,
+ #max_iter=10_000,
+ max_iter=10,
+ verbose=true)
+
+ # allocate second array of g for Jacobian solution procedure
+ g_out = g
+
+ # transfer ghost node values
+ g_out[1] = g[1]
+ g_out[end] = g[end]
+
+ # array of normalized residuals
+ local r_norm = zero(g)
+
+ for iter=1:max_iter
+
+ # perform a single jacobi iteration in each cell
+ # if iter%2 == 0
+ for iz=2:length(g)-2
+ r_norm[iz] = poisson_solver_1d_iteration(g, g_out, μ, p, iz, Δz)
+ end
+ # else # reverse direction every second time
+ # for iz=length(g)-1:2
+ # r_norm[iz] = poisson_solver_1d_iteration(g, g_out,
+ # f, iz, Δz)
+ # end
+ # end
+
+ r_norm_max = maximum(r_norm)
+
+ # Flip-flop arrays
+ tmp = g
+ g = g_out
+ g_out = tmp
+
+ if verbose
+ @info ".. Relative normalized error: $r_norm_max"
+ end
+
+ # stop iterating if the relative tolerance is satisfied
+ if r_norm_max <= rel_tol
+ if verbose
+ @info ".. Solution converged after $iter iterations"
+ end
+ return
+ end
+ end
+end
+
+function plot_profile(z, v, label, filename)
+ PyPlot.clf()
+ PyPlot.plot(v, z, "+k")
+ #PyPlot.gca().invert_yaxis()
+ PyPlot.xlabel(label)
+ PyPlot.ylabel("\$z\$ [m]")
+ PyPlot.tight_layout()
+ PyPlot.savefig(filename)
+ PyPlot.close()
+end
+
+init_μ(μ_wall, ϕ, ρ_s, G, z, P_wall) =
+ μ_wall./(1.0 .+ (1.0-ϕ)*ρ_s*G.*(L_z .- z)./P_wall)
+
+
+# Main
+
+## Set velocity BCs
+v_x[1] = v_x_bot
+
+## Calculate stressses
+p = p_lithostatic(P_wall, z)
+μ = init_μ(μ_wall, ϕ, ρ_s, G, z, P_wall)
+τ = shear_stress(v_x, Δz)
+
+#f = forcing(p, μ, g_ghost)
+#g_ghost[2:end-1] = f
+
+g_set_bc_neumann(g_ghost)
+
+implicit_1d_jacobian_poisson_solver(g_ghost, p, μ, Δz)
+
+v_x = γ_dot_p(g_ghost[2:end-1], μ)
+
+plot_profile(z, v_x, "Shear velocity, \$v_x\$ [m/s]", "1d_fd_simple_shear_v_x.png")
+plot_profile(z, μ, "Stress ratio, μ [-]", "1d_fd_simple_shear_mu.png")
+plot_profile(z, p, "Normal stress, \$p\$ [Pa]", "1d_fd_simple_shear_p.png")
+#plot_profile(z, f, "Forcing, \$f\$", "1d_fd_simple_shear_f.png")
+plot_profile(z, g_ghost[2:end-1], "Fluidity, \$g\$", "1d_fd_simple_shear_g.png")
