cngf-pf

commit 872d50f8ad089e76a93a5415628f580a3cc5b51b
parent 95cdb6dd109d1a82e9ae6d0f20e46d2929b5dd80
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date:   Fri,  5 Apr 2019 10:01:39 +0200

Include fixed bottom particles in plot, add entr command to makefile

Diffstat:
M
1d_fd_simple_shear.jl
 | 
27
++++++++++++++++++---------
M
Makefile
 | 
10
++++++++--

2 files changed, 26 insertions(+), 11 deletions(-)
diff --git a/1d_fd_simple_shear.jl b/1d_fd_simple_shear.jl
@@ -30,7 +30,8 @@ nz = 100
 ### 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.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
@@ -43,7 +44,7 @@ b = 0.9377      # Henann and Kamrin 2016
 
 ### porosity [-]
 #ϕ = 0.38        # Henann and Kamrin 2016
-ϕ = 0.28        # Damsgaard et al 2013
+ϕ = 0.25        # Damsgaard et al 2013
 
 # representative grain size [m]
 # lower values of d mean that the shear velocity curve can have sharper curves,
@@ -103,11 +104,11 @@ end
 
 ## Update ghost nodes for g from current values
 ## BC: Dirichlet (g = 0)
-function set_bc_dirichlet(g_ghost, boundary, value=0.0)
+function set_bc_dirichlet(g_ghost, boundary; value=0.0, idx_offset=0)
     if boundary == "-z"
-        g_ghost[1] = value
+        g_ghost[1+idx_offset] = value
     elseif boundary == "+z"
-        g_ghost[end] = value
+        g_ghost[end-idx_offset] = value
     else
         @error "boundary '$boundary' not understood"
     end
@@ -155,7 +156,7 @@ end
 
 ## Iteratively solve the system laplace(phi) = f
 function implicit_1d_jacobian_poisson_solver(g, p, μ, Δz,
-                                             rel_tol=1e-6,
+                                             rel_tol=1e-5,
                                              max_iter=10_000,
                                              verbose=false)
 
@@ -172,18 +173,26 @@ function implicit_1d_jacobian_poisson_solver(g, p, μ, Δz,
 
     # array of normalized residuals
     r_norm = zero(p)
+    r_norm_max = 0.0
 
     for iter=1:max_iter
 
         set_bc_neumann(g, "+z")
         set_bc_dirichlet(g, "-z")
+
+        # Damsgaard et al 2013 include fixed grains in plot
+        for i=1:9
+            set_bc_dirichlet(g, "-z", idx_offset=i)
+        end
+
         if verbose
             println("g after BC: ")
             println(g)
         end
 
         # perform a single jacobi iteration in each cell
-        for iz=1:length(p)
+        #for iz=1:length(p)
+        for iz=10:length(p)
             poisson_solver_1d_iteration(g, g_out, r_norm,
                                         μ, p, iz, Δz)
         end
@@ -215,7 +224,7 @@ function implicit_1d_jacobian_poisson_solver(g, p, μ, Δz,
             return
         end
     end
-    @error "Solution did not converge after $iter (r_norm_max = $r_norm_max)"
+    @error "Solution didn't converge after $max_iter iterations ($r_norm_max)"
 end
 
 function plot_profile(z, v, label, filename)
@@ -265,7 +274,7 @@ 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.savefig("1d_fd_simple_shear.png")
 PyPlot.close()
 
 end # end let
diff --git a/Makefile b/Makefile
@@ -1,2 +1,8 @@
-1d_fd_simple_shear_v_x.png: 1d_fd_simple_shear.jl
-	julia --banner=no --color=yes $<
+JULIA=julia --banner=no --color=yes
+
+.PHONY: run
+run: 1d_fd_simple_shear.jl
+	echo "$<" | entr -s '$(JULIA) "$<"'
+
+1d_fd_simple_shear.png: 1d_fd_simple_shear.jl
+	$(JULIA) $<