commit 1620e64a74db3220031e60593953ae346764d138
parent bf84e746592960c537206a8913ffbef07110f228
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Fri, 5 Apr 2019 16:11:36 +0200
Perform detailed comparison between implementations
Diffstat:
8 files changed, 61 insertions(+), 18 deletions(-)
diff --git a/1d_fd_simple_shear.jl b/1d_fd_simple_shear.jl
@@ -14,7 +14,8 @@ G = 9.81
### Effective normal stress exerted by top wall
# A larger normal stress deepens the deformational depth
-P_wall_ = [10, 20, 40, 60, 80, 120] .* 1e3 # normal stress [Pa]
+#P_wall_ = [10, 20, 40, 60, 80, 120] .* 1e3 # normal stress [Pa]
+P_wall_ = [10] .* 1e3 # normal stress [Pa]
### bottom velocity along x [m/s]
v_x_bot = 0.0
@@ -23,7 +24,7 @@ v_x_bot = 0.0
μ_wall = 0.40
### Nodes along z
-nz = 100
+nz = 10
## Material properties
@@ -136,7 +137,7 @@ end
## Cartesian grid. The function returns the normalized residual value
function poisson_solver_1d_iteration(g_in, g_out, r_norm,
μ, p, i, Δz,
- verbose=false)
+ verbose=true)
coorp_term = Δz^2.0/(2.0*ξ(μ[i])^2.0)
g_local_ = g_local(p[i], μ[i])
@@ -158,6 +159,7 @@ end
function implicit_1d_jacobian_poisson_solver(g, p, μ, Δz,
rel_tol=1e-5,
max_iter=10_000,
+ #max_iter=1,
verbose=false)
# allocate second array of g for Jacobian solution procedure
@@ -177,13 +179,13 @@ function implicit_1d_jacobian_poisson_solver(g, p, μ, Δz,
for iter=1:max_iter
- set_bc_neumann(g, "+z")
set_bc_dirichlet(g, "-z")
+ set_bc_neumann(g, "+z")
# Damsgaard et al 2013 include fixed grains in plot
- for i=1:9
- set_bc_dirichlet(g, "-z", idx_offset=i)
- end
+ # for i=1:9
+ # set_bc_dirichlet(g, "-z", idx_offset=i)
+ # end
if verbose
println("g after BC: ")
@@ -191,8 +193,8 @@ function implicit_1d_jacobian_poisson_solver(g, p, μ, Δz,
end
# perform a single jacobi iteration in each cell
- #for iz=1:length(p)
- for iz=10: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
@@ -227,6 +229,20 @@ function implicit_1d_jacobian_poisson_solver(g, p, μ, Δz,
@error "Solution didn't converge after $max_iter iterations ($r_norm_max)"
end
+function shear_velocity(γ_dot, Δz, v_x_bot)
+
+ v_x = zero(γ_dot)
+
+ # BC
+ v_x[1] = v_x_bot
+
+ for i=2:length(v_x)
+ v_x[i] = v_x[i-1] + γ_dot[i]*Δz
+ end
+
+ return v_x
+end
+
function plot_profile(z, v, label, filename)
PyPlot.figure(figsize=[4,4])
PyPlot.plot(v, z, "+k")
@@ -248,13 +264,24 @@ 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)
## calculate shear velocitiesj
- v_x = γ_dot_p(g_ghost[2:end-1], μ)
+ γ_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))
diff --git a/1d_fd_simple_shear.png b/1d_fd_simple_shear.png
Binary files differ.
diff --git a/1d_fd_simple_shear_damsgaard2013.h b/1d_fd_simple_shear_damsgaard2013.h
@@ -2,6 +2,7 @@
#define ONED_FD_SIMPLE_SHEAR_
#include <math.h>
+#include <stdio.h>
#include "arrays.h"
#include "simulation.h"
@@ -17,6 +18,7 @@ struct simulation init_sim(void)
sim.P_wall = 10e3; /* larger normal stress deepens the shear depth */
sim.mu_wall = 0.40;
+ sim.v_x_bot = 0.0;
sim.nz = 10;
@@ -57,7 +59,8 @@ void init_friction(struct simulation* sim)
{
for (int i=0; i<sim->nz; ++i)
sim->mu[i] = sim->mu_wall /
- (1.0 + (1.0 - sim->phi)*sim->rho_s*sim->G*(sim->L_z - sim->z[i]));
+ (1.0 + (1.0 - sim->phi)*sim->rho_s*sim->G*(sim->L_z - sim->z[i])/
+ sim->P_wall);
}
#endif
diff --git a/1d_fd_simple_shear_henann_kamrin2016.h b/1d_fd_simple_shear_henann_kamrin2016.h
@@ -16,6 +16,7 @@ struct simulation init_sim(void)
sim.P_wall = 10e3; /* larger normal stress deepens the deformational depth */
sim.mu_wall = 0.40;
+ sim.v_x_bot = 0.0;
sim.nz = 100;
diff --git a/arrays.c b/arrays.c
@@ -47,9 +47,10 @@ unsigned int idx1g(const unsigned int i)
double* linspace(const double lower, const double upper, const int n)
{
double *x = malloc(n*sizeof(double));
- double dx = (upper - lower)/(double)n;
- for (int i=0; i<n; ++i)
+ double dx = (upper - lower)/(double)(n-1);
+ for (int i=0; i<n; ++i) {
x[i] = lower + dx*i;
+ printf("x[%d] = %f\n", i, x[i]);}
return x;
}
diff --git a/main.c b/main.c
@@ -26,6 +26,7 @@ int main(int argc, char** argv) {
exit(1);
compute_shear_strain_rate_plastic(&sim);
+ compute_shear_velocity(&sim);
puts("\n## After solver");
puts(".. g:"); print_array(sim.g_ghost, sim.nz+2);
diff --git a/simulation.c b/simulation.c
@@ -40,6 +40,15 @@ void compute_shear_strain_rate_plastic(struct simulation* sim)
shear_strain_rate_plastic(sim->g_ghost[idx1g(i)], sim->mu[i]);
}
+void compute_shear_velocity(struct simulation* sim)
+{
+ // Dirichlet BC at bottom
+ sim->v_x[0] = sim->v_x_bot;
+
+ for (int i=1; i<sim->nz; ++i)
+ sim->v_x[i] = sim->v_x[i-1] + sim->gamma_dot_p[i]*sim->dz;
+}
+
double cooperativity_length(
const double A,
const double d,
@@ -126,14 +135,14 @@ void poisson_solver_1d_cell_update(
+ g_in[gi+1]/2.0
+ g_in[gi-1]/2.0);
- r_norm[i] = pow(g_out[gi] - g_in[gi], 1.0) / (pow(g_out[gi], 2.0) + 1e-16);
+ r_norm[i] = pow(g_out[gi] - g_in[gi], 2.0) / (pow(g_out[gi], 2.0) + 1e-16);
printf("-- %d --------------\n", i);
printf("coorp_term: %g\n", coorp_term);
printf(" g_local: %g\n", local_fluidity(p[i], mu[i], mu_s, b, rho_s, d));
printf(" g_in: %g\n", g_in[gi]);
printf(" g_out: %g\n", g_out[gi]);
- printf(" r_norm: %g\n", r_norm[gi]);
+ printf(" r_norm: %g\n", r_norm[i]);
}
int implicit_1d_jacobian_poisson_solver(
@@ -144,14 +153,14 @@ int implicit_1d_jacobian_poisson_solver(
double* g_ghost_out = empty(sim->nz);
double* r_norm = empty(sim->nz);
- set_bc_neumann(sim->g_ghost, sim->nz, +1);
- set_bc_dirichlet(sim->g_ghost, sim->nz, -1, 0.0);
-
int iter;
double r_norm_max = NAN;
for (iter=0; iter<max_iter; ++iter) {
+ set_bc_dirichlet(sim->g_ghost, sim->nz, -1, 0.0);
+ set_bc_neumann(sim->g_ghost, sim->nz, +1);
+
for (int i=0; i<sim->nz; ++i)
poisson_solver_1d_cell_update(
i,
diff --git a/simulation.h b/simulation.h
@@ -78,6 +78,7 @@ void set_bc_dirichlet(
void compute_cooperativity_length(struct simulation* sim);
void compute_shear_strain_rate_plastic(struct simulation* sim);
+void compute_shear_velocity(struct simulation* sim);
int implicit_1d_jacobian_poisson_solver(
struct simulation* sim,
