tRemove Julia implementation - cngf-pf - continuum model for granular flows wit… | |
git clone git://src.adamsgaard.dk/cngf-pf | |
Log | |
Files | |
Refs | |
README | |
LICENSE | |
--- | |
commit edc6affeb0cb2691c0db6ea8ef5043f3a43b1e3e | |
parent bfa41cc2e556c1fb56b794cac4505371890dc86b | |
Author: Anders Damsgaard <[email protected]> | |
Date: Mon, 15 Apr 2019 13:05:05 +0200 | |
Remove Julia implementation | |
Diffstat: | |
D julia/1d_fd_simple_shear.jl | 277 -----------------------------… | |
D julia/1d_fd_simple_shear.png | 0 | |
D julia/Makefile | 7 ------- | |
3 files changed, 0 insertions(+), 284 deletions(-) | |
--- | |
diff --git a/julia/1d_fd_simple_shear.jl b/julia/1d_fd_simple_shear.jl | |
t@@ -1,277 +0,0 @@ | |
-#!/usr/bin/env julia | |
- | |
-ENV["MPLBACKEND"] = "Agg" | |
-import PyPlot | |
- | |
-let | |
- | |
-# Simulation settings | |
- | |
-## Gravitational acceleration magnitude | |
-G = 9.81 | |
- | |
-## Wall parameters | |
- | |
-### Effective normal stress exerted by top wall | |
-# A larger normal stress deepens the deformational depth | |
-P_wall_ = [10, 20, 40, 60, 80, 120] .* 1e3 # normal stress [Pa] | |
- | |
-### bottom velocity along x [m/s] | |
-v_x_bot = 0.0 | |
- | |
-# stress ratio at top wall | |
-μ_wall = 0.40 | |
- | |
-### Nodes along z | |
-nz = 100 | |
- | |
-## Material properties | |
- | |
-### nonlocal amplitude [-] | |
-# lower values of A mean that the velocity curve can have sharper curves, e.g. | |
-# at the transition from μ ≈ μ_s | |
-#A = 0.48 # Henann and Kamrin 2016 | |
-A = 0.40 # Loose fit to Damsgaard et al 2013 | |
- | |
-### rate dependence beyond yield [-] | |
-# lower values of b mean larger shear velocity for a given stress ratio above | |
-# yield | |
-b = 0.9377 # Henann and Kamrin 2016 | |
- | |
-### bulk and critical state static yield friction coefficient [-] | |
-#μ_s = 0.3819 # Henann and Kamrin 2016 | |
-μ_s = atan(deg2rad(22.0)) # Damsgaard et al 2013 | |
- | |
-### porosity [-] | |
-#ϕ = 0.38 # Henann and Kamrin 2016 | |
-ϕ = 0.25 # Damsgaard et al 2013 | |
- | |
-# representative grain size [m] | |
-# lower values of d mean that the shear velocity curve can have sharper curves, | |
-# e.g. at the transition from μ ≈ μ_s | |
-#d = 1e-3 # Henann and Kamrin 2016 | |
-d = 0.04 # Damsgaard et al 2013 | |
- | |
-### grain material density [kg/m^3] | |
-#ρ_s = 2.485e3 # Henann and Kamrin 2016 | |
-ρ_s = 2.6e3 # Damsgaard et al 2013 | |
- | |
-## Spatial settings | |
-origo_z = 0.0 | |
-#L_z = 20.0*d # Henann and Kamrin 2016 | |
-L_z = 0.7 # Damsgaard et al 2013 | |
-z = collect(range(origo_z, L_z, length=nz)) | |
-Δ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.*μ | |
- | |
-## 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 set_bc_neumann(g_ghost, boundary) | |
- if boundary == "-z" | |
- g_ghost[1] = g_ghost[2] | |
- elseif boundary == "+z" | |
- g_ghost[end] = g_ghost[end-1] | |
- else | |
- @error "boundary '$boundary' not understood" | |
- end | |
-end | |
- | |
-## Update ghost nodes for g from current values | |
-## BC: Dirichlet (g = 0) | |
-function set_bc_dirichlet(g_ghost, boundary; value=0.0, idx_offset=0) | |
- if boundary == "-z" | |
- g_ghost[1+idx_offset] = value | |
- elseif boundary == "+z" | |
- g_ghost[end-idx_offset] = value | |
- else | |
- @error "boundary '$boundary' not understood" | |
- end | |
-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) | |
- | |
-## A single iteration for solving a Poisson equation Laplace(phi) = f on a | |
-## Cartesian grid. The function returns the normalized residual value | |
-function poisson_solver_1d_iteration(g_in, g_out, r_norm, | |
- μ, p, i, Δz, | |
- verbose=false) | |
- | |
- coorp_term = Δz^2.0/(2.0*ξ(μ[i])^2.0) | |
- g_out[i+1] = 1.0/(1.0 + coorp_term) * (coorp_term*g_local(p[i], μ[i]) | |
- + g_in[i+1+1]/2.0 + g_in[i-1+1]/2.0) | |
- r_norm[i] = (g_out[i+1] - g_in[i+1])^2.0 / (g_out[i+1]^2.0 + 1e-16) | |
- | |
- if verbose | |
- println("-- $i -----------------") | |
- println("coorp_term: $coorp_term") | |
- println(" g_local: $(g_local(p[i], μ[i]))") | |
- println(" g_in: $(g_in[i+1])") | |
- println(" g_out: $(g_out[i+1])") | |
- println(" r_norm: $(r_norm[i])") | |
- end | |
-end | |
- | |
-## Iteratively solve the system laplace(phi) = f | |
-function implicit_1d_jacobian_poisson_solver(g, p, μ, Δz, | |
- rel_tol=1e-5, | |
- max_iter=10_000, | |
- verbose=false) | |
- | |
- # allocate second array of g for Jacobian solution procedure | |
- g_out = zero(g) | |
- | |
- if verbose | |
- println("g: ") | |
- println(g) | |
- println() | |
- println("g_out: ") | |
- println(g_out) | |
- end | |
- | |
- # array of normalized residuals | |
- r_norm = zero(p) | |
- r_norm_max = 0.0 | |
- | |
- for iter=1:max_iter | |
- #println("\n@@@ ITERATION $iter @@@") | |
- | |
- set_bc_dirichlet(g, "-z") | |
- set_bc_dirichlet(g, "+z") | |
- #set_bc_neumann(g, "+z") | |
- | |
- if verbose | |
- println("g after BC: ") | |
- println(g) | |
- end | |
- | |
- # perform a single jacobi iteration in each cell | |
- for iz=1:length(p) | |
- poisson_solver_1d_iteration(g, g_out, r_norm, | |
- μ, p, iz, Δz) | |
- 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 | |
- @info ".. Solution converged after $iter iterations" | |
- return | |
- end | |
- end | |
- @error "Solution didn't converge after $max_iter iterations ($r_norm_max)" | |
-end | |
- | |
-function shear_velocity(γ_dot, Δz, v_x_bot) | |
- | |
- v_x = zero(γ_dot) | |
- | |
- # BC | |
- v_x[1] = v_x_bot | |
- | |
- for i=2:length(v_x) | |
- v_x[i] = v_x[i-1] + γ_dot[i]*Δz | |
- end | |
- | |
- return v_x | |
-end | |
- | |
-function plot_profile(z, v, label, filename) | |
- PyPlot.figure(figsize=[4,4]) | |
- PyPlot.plot(v, z, "+k") | |
- 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 | |
- | |
-for P_wall in P_wall_ | |
- | |
- ## calculate stresses | |
- p = p_lithostatic(P_wall, z) | |
- μ = init_μ(μ_wall, ϕ, ρ_s, G, z, P_wall) | |
- | |
- ## solve for fluidity | |
- implicit_1d_jacobian_poisson_solver(g_ghost, p, μ, Δz) | |
- | |
- ## calculate shear velocitiesj | |
- γ_dot = γ_dot_p(g_ghost[2:end-1], μ) | |
- v_x = shear_velocity(γ_dot, Δz, v_x_bot) | |
- | |
- ## plot results | |
- P = Int(round(P_wall/1e3)) | |
- PyPlot.plot(v_x/maximum(v_x), z, "+-", label="\$P_{wall}\$ = $P kPa") | |
- # plot_profile(z, v_x, "Shear velocity, \$v_x\$ [m/s]", | |
- # "1d_fd_simple_shear_v_x_P$(P)kPa.png") | |
- # plot_profile(z, μ, "Stress ratio, μ [-]", | |
- # "1d_fd_simple_shear_mu_P$(P)kPa.png") | |
- # plot_profile(z, p, "Normal stress, \$p\$ [Pa]", | |
- # "1d_fd_simple_shear_p_P$(P)kPa.png") | |
- # plot_profile(z, g_ghost[2:end-1], "Fluidity, \$g\$", | |
- # "1d_fd_simple_shear_g_P$(P)kPa.png") | |
-end | |
- | |
-PyPlot.xlabel("Normalized shear displacement, [m]") | |
-PyPlot.ylabel("Vertical position, \$z\$ [m]") | |
-PyPlot.legend() | |
-PyPlot.tight_layout() | |
-PyPlot.savefig("1d_fd_simple_shear.png") | |
-PyPlot.close() | |
- | |
-end # end let | |
diff --git a/julia/1d_fd_simple_shear.png b/julia/1d_fd_simple_shear.png | |
Binary files differ. | |
diff --git a/julia/Makefile b/julia/Makefile | |
t@@ -1,7 +0,0 @@ | |
-JULIA=julia --banner=no --color=yes | |
-.PHONY: run-julia | |
-run-julia: 1d_fd_simple_shear.jl | |
- echo "$<" | entr -s '$(JULIA) "$<"' | |
- | |
-1d_fd_simple_shear.png: 1d_fd_simple_shear.jl | |
- $(JULIA) $< |