cngf-pf

commit 95cdb6dd109d1a82e9ae6d0f20e46d2929b5dd80
parent c1aa4dce15eba37cc67913aa232337e9205feb70
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date:   Fri,  5 Apr 2019 09:15:14 +0200

Add Dirichlet BC and improve BC calls

Diffstat:
M
1d_fd_simple_shear.jl
 | 
58
++++++++++++++++++++++++++++++++++++++++++++++------------

1 file changed, 46 insertions(+), 12 deletions(-)
diff --git a/1d_fd_simple_shear.jl b/1d_fd_simple_shear.jl
@@ -90,10 +90,27 @@ function g_local(p, μ)
 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]
+## 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)
+    if boundary == "-z"
+        g_ghost[1] = value
+    elseif boundary == "+z"
+        g_ghost[end] = value
+    else
+        @error "boundary '$boundary' not understood"
+    end
 end
 
 ## Compute shear stress from velocity gradient using finite differences
@@ -158,7 +175,8 @@ function implicit_1d_jacobian_poisson_solver(g, p, μ, Δz,
 
     for iter=1:max_iter
 
-        g_set_bc_neumann(g)
+        set_bc_neumann(g, "+z")
+        set_bc_dirichlet(g, "-z")
         if verbose
             println("g after BC: ")
             println(g)
@@ -215,23 +233,39 @@ init_μ(μ_wall, ϕ, ρ_s, G, z, P_wall) =
 
 
 # Main
+PyPlot.clf()
 
 for P_wall in P_wall_
-    ## Set velocity BCs
-    v_x[1] = v_x_bot
 
-    ## Calculate stressses
+    ## 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
     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_P$(P_wall).png")
-    plot_profile(z, μ, "Stress ratio, μ [-]", "1d_fd_simple_shear_mu_P$(P_wall).png")
-    plot_profile(z, p, "Normal stress, \$p\$ [Pa]", "1d_fd_simple_shear_p_P$(P_wall).png")
-    plot_profile(z, g_ghost[2:end-1], "Fluidity, \$g\$", "1d_fd_simple_shear_g_P$(P_wall).png")
+
+    ## 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_v_x.png")
+PyPlot.close()
+
 end # end let