cngf-pf

commit 4def6527e85a80d0dc2b52f51aa376c782d296d4
parent f65899ee85de5ed0234d7b2164c22ff80b25a978
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date:   Thu,  4 Apr 2019 16:29:15 +0200

Add working script

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

1 file changed, 60 insertions(+), 34 deletions(-)
diff --git a/1d_fd_simple_shear.jl b/1d_fd_simple_shear.jl
@@ -9,13 +9,14 @@ import PyPlot
 const G = 9.81
 
 ## Wall parameters
-const P_wall = 10.0e3 # normal stress [Pa]
+const P_wall = 1.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
+#const μ_wall = 0.392  # stress ratio at top wall
+const μ_wall = 0.52  # stress ratio at top wall
 
 ## Nodes along z
-const nz = 10
+const nz = 100
 
 ## Material properties
 const A = 0.48        # nonlocal amplitude [-]
@@ -71,22 +72,22 @@ function g_set_bc_neumann(g_ghost)
 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
+# 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)
+#friction(τ, p) = τ./(p .+ 1e-16)
 
 ## Fluidity
 #fluidity(p, μ, g) = 1.0./(ξ.(μ).^2.0) .* (g .- g_local.(p, μ))
@@ -102,38 +103,57 @@ 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
-# TODO f is not used?
-function poisson_solver_1d_iteration(g_in, g_out, μ, p, i, Δz)
+function poisson_solver_1d_iteration(g_in, g_out, μ, p, i, Δz, verbose=false)
 
     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)
+    g_local_ = g_local(p[i], μ[i])
+    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)
+    res_norm = (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_")
+        println("      g_in: $(g_in[i+1])")
+        println("     g_out: $(g_out[i+1])")
+        println("  res_norm: $res_norm")
+    end
+    return res_norm
 end
 
 ## Iteratively solve the system laplace(phi) = f
 function implicit_1d_jacobian_poisson_solver(g, p, μ, Δz,
-                                             rel_tol=1e-4,
+                                             rel_tol=1e-5,
                                              #max_iter=10_000,
                                              max_iter=10,
-                                             verbose=true)
+                                             verbose=false)
 
     # 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]
+    if verbose
+        println("g: ")
+        println(g)
+        println()
+        println("g_out: ")
+        println(g_out)
+    end
 
     # array of normalized residuals
     local r_norm = zero(g)
 
     for iter=1:max_iter
 
+        g_set_bc_neumann(g)
+        if verbose
+            println("g after BC: ")
+            println(g)
+        end
+
         # perform a single jacobi iteration in each cell
         # if iter%2 == 0
-            for iz=2:length(g)-2
+            for iz=1:length(p)
                 r_norm[iz] = poisson_solver_1d_iteration(g, g_out, μ, p, iz, Δz)
             end
         # else  # reverse direction every second time
@@ -142,9 +162,17 @@ function implicit_1d_jacobian_poisson_solver(g, p, μ, Δz,
         #                                                  f, iz, Δz)
         #     end
         # end
-
         r_norm_max = maximum(r_norm)
 
+        if verbose
+            println("g after jacobi iteration $iter: ")
+            println(g)
+            println("r_norm:")
+            println(r_norm)
+            println("r_norm_max:")
+            println(r_norm_max)
+        end
+
         # Flip-flop arrays
         tmp = g
         g = g_out
@@ -187,13 +215,11 @@ 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)
+#τ = 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], μ)