1d_fd_simple_shear

Continuum model for granular flows with pore-pressure dynamics
git clone git://src.adamsgaard.dk/1d_fd_simple_shear
Log | Files | Refs | README | LICENSE

commit 6cb6fc35f0aa16cadba133b4bd63f7376c7cdc58
parent 1620e64a74db3220031e60593953ae346764d138
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date:   Fri,  5 Apr 2019 20:08:43 +0200

Implement final fixes, minimize stdout output

Diffstat:
M1d_fd_simple_shear.jl | 13++-----------
Marrays.c | 17++++++++++++++---
Marrays.h | 3+++
Mmain.c | 6+++---
Msimulation.c | 16+++++++++++-----
5 files changed, 33 insertions(+), 22 deletions(-)

diff --git a/1d_fd_simple_shear.jl b/1d_fd_simple_shear.jl @@ -140,7 +140,6 @@ function poisson_solver_1d_iteration(g_in, g_out, r_norm, verbose=true) 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) r_norm[i] = (g_out[i+1] - g_in[i+1])^2.0 / (g_out[i+1]^2.0 + 1e-16) @@ -148,7 +147,7 @@ function poisson_solver_1d_iteration(g_in, g_out, r_norm, if verbose println("-- $i -----------------") println("coorp_term: $coorp_term") - println(" g_local: $g_local_") + println(" g_local: $(g_local(p[i], μ[i]))") println(" g_in: $(g_in[i+1])") println(" g_out: $(g_out[i+1])") println(" r_norm: $(r_norm[i])") @@ -178,6 +177,7 @@ function implicit_1d_jacobian_poisson_solver(g, p, μ, Δz, r_norm_max = 0.0 for iter=1:max_iter + println("\n@@@ ITERATION $iter @@@") set_bc_dirichlet(g, "-z") set_bc_neumann(g, "+z") @@ -200,15 +200,6 @@ function implicit_1d_jacobian_poisson_solver(g, p, μ, Δz, 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 diff --git a/arrays.c b/arrays.c @@ -48,9 +48,8 @@ double* linspace(const double lower, const double upper, const int n) { double *x = malloc(n*sizeof(double)); double dx = (upper - lower)/(double)(n-1); - for (int i=0; i<n; ++i) { + for (int i=0; i<n; ++i) x[i] = lower + dx*i; - printf("x[%d] = %f\n", i, x[i]);} return x; } @@ -101,5 +100,17 @@ double min(const double* a, const int n) void print_array(const double* a, const int n) { for (int i=0; i<n; ++i) - printf("%g\n", a[i]); + printf("%.17g\n", a[i]); +} + +void print_arrays(const double* a, const double* b, const int n) +{ + for (int i=0; i<n; ++i) + printf("%.17g\t%.17g\n", a[i], b[i]); +} + +void copy_values(const double* in, double* out, const int n) +{ + for (int i=0; i<n; ++i) + out[i] = in[i]; } diff --git a/arrays.h b/arrays.h @@ -26,5 +26,8 @@ double max(const double* a, const int n); double min(const double* a, const int n); void print_array(const double* a, const int n); +void print_arrays(const double* a, const double* b, const int n); + +void copy_values(const double* in, double* out, const int n); #endif diff --git a/main.c b/main.c @@ -18,9 +18,11 @@ int main(int argc, char** argv) { init_friction(&sim); compute_cooperativity_length(&sim); +#ifdef DEBUG puts("\n## Before solver"); puts(".. p:"); print_array(sim.p, sim.nz); puts(".. mu:"); print_array(sim.mu, sim.nz); +#endif if (implicit_1d_jacobian_poisson_solver(&sim, 10000, 1e-5)) exit(1); @@ -28,9 +30,7 @@ int main(int argc, char** argv) { compute_shear_strain_rate_plastic(&sim); compute_shear_velocity(&sim); - puts("\n## After solver"); - puts(".. g:"); print_array(sim.g_ghost, sim.nz+2); - puts(".. v_x:"); print_array(sim.v_x, sim.nz); + print_arrays(sim.z, sim.v_x, sim.nz); free_arrays(&sim); return 0; diff --git a/simulation.c b/simulation.c @@ -42,6 +42,7 @@ void compute_shear_strain_rate_plastic(struct simulation* sim) void compute_shear_velocity(struct simulation* sim) { + // TODO: implement iterative solver // Dirichlet BC at bottom sim->v_x[0] = sim->v_x_bot; @@ -91,7 +92,7 @@ void set_bc_neumann(double* g_ghost, const int nz, const int boundary) if (boundary == -1) g_ghost[0] = g_ghost[1]; else if (boundary == +1) - g_ghost[idx1g(nz)+1] = g_ghost[idx1g(nz)]; + g_ghost[nz+1] = g_ghost[nz]; else { fprintf(stderr, "set_bc_neumann: Unknown boundary %d\n", boundary); exit(1); @@ -137,12 +138,14 @@ void poisson_solver_1d_cell_update( r_norm[i] = pow(g_out[gi] - g_in[gi], 2.0) / (pow(g_out[gi], 2.0) + 1e-16); +#ifdef DEBUG 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[i]); +#endif } int implicit_1d_jacobian_poisson_solver( @@ -157,6 +160,9 @@ int implicit_1d_jacobian_poisson_solver( double r_norm_max = NAN; for (iter=0; iter<max_iter; ++iter) { +#ifdef DEBUG + printf("\n@@@ ITERATION %d @@@\n", iter); +#endif set_bc_dirichlet(sim->g_ghost, sim->nz, -1, 0.0); set_bc_neumann(sim->g_ghost, sim->nz, +1); @@ -177,14 +183,14 @@ int implicit_1d_jacobian_poisson_solver( sim->d); r_norm_max = max(r_norm, sim->nz); - double* tmp = sim->g_ghost; - sim->g_ghost = g_ghost_out; - g_ghost_out = tmp; + copy_values(g_ghost_out, sim->g_ghost, sim->nz+2); if (r_norm_max <= rel_tol) { + set_bc_dirichlet(sim->g_ghost, sim->nz, -1, 0.0); + set_bc_neumann(sim->g_ghost, sim->nz, +1); free(g_ghost_out); free(r_norm); - printf(".. Solution converged after %d iterations\n", iter); + /* printf(".. Solution converged after %d iterations\n", iter); */ return 0; } }