1d_fd_simple_shear_transient

commit 5e5b568af441c1457133afd2deb7e09867fbb31e
parent 3bd8e7c3f474481aa9a925151da29945bac8e8e5
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date:   Thu,  4 Apr 2019 21:53:21 +0200

Fix array copying during jacobi iterations

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

1 file changed, 36 insertions(+), 26 deletions(-)
diff --git a/1d_fd_simple_shear.jl b/1d_fd_simple_shear.jl
@@ -9,26 +9,41 @@ import PyPlot
 const G = 9.81
 
 ## Wall parameters
-const P_wall = 1.0e3  # normal stress [Pa]
-# const v_x_wall = 1e-0 # wall velocity along x [m/s]
+
+# A larger normal stress deepens the deformational depth
+const P_wall = 100.0e3  # normal stress [Pa]
+
 const v_x_bot = 0.0   # bottom velocity along x [m/s]
-#const μ_wall = 0.392  # stress ratio at top wall
-const μ_wall = 0.52  # stress ratio at top wall
+const μ_wall = 0.40  # stress ratio at top wall
 
 ## Nodes along z
 const nz = 100
 
 ## Material properties
+
+# lower values of A mean that the velocity curve can have sharper curves, e.g.
+# at the transition from μ ≈ μ_s
 const A = 0.48        # nonlocal amplitude [-]
+
+
+# lower values of b mean larger shear velocity for a given stress ratio above
+# yield
 const b = 0.9377      # rate dependence beyond yield [-]
+
 const μ_s = 0.3819    # static yield friction coefficient [-]
 const ϕ = 0.38        # porosity [-]
+
+# lower values of d mean that the shear velocity curve can have sharper curves,
+# e.g. at the transition from μ ≈ μ_s
 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
+#const L_z = 20.0*d
+#const L_z = 200.0*d
+const L_z = 0.5
 
 ## Spatial coordiantes
 const z = collect(range(origo_z, L_z, length=nz))
@@ -103,13 +118,15 @@ end
 
 ## 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, μ, p, i, Δz, verbose=false)
+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_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)
+    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 -----------------")
@@ -117,20 +134,18 @@ function poisson_solver_1d_iteration(g_in, g_out, μ, p, i, Δz, verbose=false)
         println("   g_local: $g_local_")
         println("      g_in: $(g_in[i+1])")
         println("     g_out: $(g_out[i+1])")
-        println("  res_norm: $res_norm")
+        println("    r_norm: $(r_norm[i])")
     end
-    return res_norm
 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,
-                                             max_iter=10,
+                                             rel_tol=1e-6,
+                                             max_iter=10_000,
                                              verbose=false)
 
     # allocate second array of g for Jacobian solution procedure
-    g_out = g
+    g_out = zero(g)
 
     if verbose
         println("g: ")
@@ -141,7 +156,7 @@ function implicit_1d_jacobian_poisson_solver(g, p, μ, Δz,
     end
 
     # array of normalized residuals
-    local r_norm = zero(g)
+    r_norm = zero(p)
 
     for iter=1:max_iter
 
@@ -152,16 +167,10 @@ function implicit_1d_jacobian_poisson_solver(g, p, μ, Δz,
         end
 
         # perform a single jacobi iteration in each cell
-        # if iter%2 == 0
-            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
-        #     for iz=length(g)-1:2
-        #         r_norm[iz] = poisson_solver_1d_iteration(g, g_out,
-        #                                                  f, iz, Δz)
-        #     end
-        # end
+        for iz=1:length(p)
+            poisson_solver_1d_iteration(g, g_out, r_norm,
+                                        μ, p, iz, Δz)
+        end
         r_norm_max = maximum(r_norm)
 
         if verbose
@@ -184,12 +193,13 @@ function implicit_1d_jacobian_poisson_solver(g, p, μ, Δz,
 
         # stop iterating if the relative tolerance is satisfied
         if r_norm_max <= rel_tol
-            if verbose
+            #if verbose
                 @info ".. Solution converged after $iter iterations"
-            end
+            #end
             return
         end
     end
+    @error "Solution did not converge after $iter (r_norm_max = $r_norm_max)"
 end
 
 function plot_profile(z, v, label, filename)