cngf-pf

commit bbd24afac544c4afc9772f126eec9c3ba2c602e4
parent 2e2fa53fd466f6e8a224249b5f353cf0029f257a
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date:   Sat,  6 Apr 2019 12:40:28 +0200

Use Dirichlet BC for top

Diffstat:
M
1d_fd_simple_shear.png
 | 
0

M
julia/1d_fd_simple_shear.jl
 | 
16
+++-------------
M
julia/1d_fd_simple_shear.png
 | 
0

M
simulation.c
 | 
7
+++++--

4 files changed, 8 insertions(+), 15 deletions(-)
diff --git a/1d_fd_simple_shear.png b/1d_fd_simple_shear.png
Binary files differ.
diff --git a/julia/1d_fd_simple_shear.jl b/julia/1d_fd_simple_shear.jl
@@ -175,10 +175,11 @@ function implicit_1d_jacobian_poisson_solver(g, p, μ, Δz,
     r_norm_max = 0.0
 
     for iter=1:max_iter
-        println("\n@@@ ITERATION $iter @@@")
+        #println("\n@@@ ITERATION $iter @@@")
 
         set_bc_dirichlet(g, "-z")
-        set_bc_neumann(g, "+z")
+        set_bc_dirichlet(g, "+z")
+        #set_bc_neumann(g, "+z")
 
         if verbose
             println("g after BC: ")
@@ -245,12 +246,6 @@ for P_wall in P_wall_
     ## calculate stresses
     p = p_lithostatic(P_wall, z)
     μ = init_μ(μ_wall, ϕ, ρ_s, G, z, P_wall)
-    println(".. z:")
-    println(z)
-    println(".. p:")
-    println(p)
-    println(".. mu:")
-    println(μ)
 
     ## solve for fluidity
     implicit_1d_jacobian_poisson_solver(g_ghost, p, μ, Δz)
@@ -259,11 +254,6 @@ for P_wall in P_wall_
     γ_dot = γ_dot_p(g_ghost[2:end-1], μ)
     v_x = shear_velocity(γ_dot, Δz, v_x_bot)
 
-    println(".. g:")
-    println(g_ghost)
-    println(".. v_x:")
-    println(v_x)
-
     ## plot results
     P = Int(round(P_wall/1e3))
     PyPlot.plot(v_x/maximum(v_x), z, "+-", label="\$P_{wall}\$ = $P kPa")
diff --git a/julia/1d_fd_simple_shear.png b/julia/1d_fd_simple_shear.png
Binary files differ.
diff --git a/simulation.c b/simulation.c
@@ -163,9 +163,12 @@ int implicit_1d_jacobian_poisson_solver(
 #ifdef DEBUG
         printf("\n@@@ ITERATION %d @@@\n", iter);
 #endif
-
+        /* Dirichlet BCs resemble fixed particle velocities */
         set_bc_dirichlet(sim->g_ghost, sim->nz, -1, 0.0);
-        set_bc_neumann(sim->g_ghost, sim->nz, +1);
+        set_bc_dirichlet(sim->g_ghost, sim->nz, +1, 0.0);
+
+        /* Neumann BCs resemble free surfaces */
+        /* set_bc_neumann(sim->g_ghost, sim->nz, +1); */
 
         for (int i=0; i<sim->nz; ++i)
             poisson_solver_1d_cell_update(