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 464550e2aba5f3dcec7454d023843dae46aafc24
parent 7631c7da68bfda4e8924457ca4d7615cd2fe1862
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date:   Mon,  1 Jul 2019 10:24:47 +0200

Save pressure solution as dp_f/dt terms, reduce memory operations

Diffstat:
Mfluid.c | 149++++++++++++++++++++++++++++++++++++++++---------------------------------------
1 file changed, 75 insertions(+), 74 deletions(-)

diff --git a/fluid.c b/fluid.c @@ -67,7 +67,6 @@ darcy_pressure_change_1d(const int i, const double* phi, const double* k, const double dz, - const double dt, const double beta_f, const double mu_f) { @@ -96,7 +95,7 @@ darcy_pressure_change_1d(const int i, i, phi[i], i, div_k_grad_p); #endif - return dt/(beta_f*phi[i]*mu_f)*div_k_grad_p; + return 1.0/(beta_f*phi[i]*mu_f)*div_k_grad_p; } int @@ -104,15 +103,16 @@ darcy_solver_1d(struct simulation* sim, const int max_iter, const double rel_tol) { - int i, iter; + int i, iter, solved; double epsilon, theta, p_f, p_f_top, r_norm_max; - double *dp_f_expl, *dp_f_impl, *p_f_ghost_out, *r_norm; + double *dp_f_dt_expl, *dp_f_dt_impl, *p_f_ghost_out, *r_norm; /* choose integration method, parameter in [0.0; 1.0] * epsilon = 0.0: explicit * epsilon = 0.5: Crank-Nicolson * epsilon = 1.0: implicit */ - epsilon = 0.5; + /* epsilon = 0.5; */ + epsilon = 0.0; /* choose relaxation factor, parameter in ]0.0; 1.0] * theta in ]0.0; 1.0]: underrelaxation @@ -133,94 +133,95 @@ darcy_solver_1d(struct simulation* sim, set_bc_neumann(sim->p_f_ghost, sim->nz, -1); /* compute explicit solution to pressure change */ - dp_f_expl = zeros(sim->nz); + solved = 1; + dp_f_dt_expl = zeros(sim->nz); for (i=0; i<sim->nz; ++i) - dp_f_expl[i] = darcy_pressure_change_1d(i, - sim->nz, - sim->p_f_ghost, - sim->phi, - sim->k, - sim->dz, - sim->dt, - sim->beta_f, - sim->mu_f); - - /* compute implicit solution to pressure change */ - dp_f_impl = zeros(sim->nz); - p_f_ghost_out = zeros(sim->nz+2); - r_norm = zeros(sim->nz); - for (iter=0; iter<max_iter; ++iter) { - - /* set fluid BCs (2 of 2) */ - set_bc_dirichlet(sim->p_f_ghost, sim->nz, +1, p_f_top); - sim->p_f_ghost[idx1g(sim->nz-1)] = p_f_top; /* Include top node in BC */ - set_bc_neumann(sim->p_f_ghost, sim->nz, -1); + dp_f_dt_expl[i] = darcy_pressure_change_1d(i, + sim->nz, + sim->p_f_ghost, + sim->phi, + sim->k, + sim->dz, + sim->beta_f, + sim->mu_f); + + if (epsilon > 1e-3) { + solved = 0; + /* compute implicit solution to pressure change */ + dp_f_dt_impl = zeros(sim->nz); + p_f_ghost_out = zeros(sim->nz+2); + r_norm = zeros(sim->nz); + for (iter=0; iter<max_iter; ++iter) { + + /* set fluid BCs (2 of 2) */ + set_bc_dirichlet(sim->p_f_ghost, sim->nz, +1, p_f_top); + sim->p_f_ghost[idx1g(sim->nz-1)] = p_f_top; + set_bc_neumann(sim->p_f_ghost, sim->nz, -1); #ifdef DEBUG - puts(".. p_f_ghost after BC:"); print_array(sim->p_f_ghost, sim->nz+2); + puts(".. p_f_ghost after BC:"); + print_array(sim->p_f_ghost, sim->nz+2); #endif - /* for (int i=0; i<sim->nz; ++i) */ - for (i=0; i<sim->nz-1; ++i) - dp_f_impl[i] = darcy_pressure_change_1d(i, - sim->nz, - sim->p_f_ghost, - sim->phi, - sim->k, - sim->dz, - sim->dt, - sim->beta_f, - sim->mu_f); - for (i=0; i<sim->nz-1; ++i) { + /* for (int i=0; i<sim->nz; ++i) */ + for (i=0; i<sim->nz-1; ++i) + dp_f_dt_impl[i] = darcy_pressure_change_1d(i, + sim->nz, + sim->p_f_ghost, + sim->phi, + sim->k, + sim->dz, + sim->beta_f, + sim->mu_f); + for (i=0; i<sim->nz-1; ++i) { #ifdef DEBUG - printf("dp_f_expl[%d] = %g\ndp_f_impl[%d] = %g\n", - i, dp_f_expl[i], i, dp_f_impl[i]); + printf("dp_f_expl[%d] = %g\ndp_f_impl[%d] = %g\n", + i, dp_f_expl[i], i, dp_f_impl[i]); #endif - p_f = sim->p_f_ghost[idx1g(i)]; + p_f = sim->p_f_ghost[idx1g(i)]; - p_f_ghost_out[idx1g(i)] = p_f - + (1.0 - epsilon)*dp_f_expl[i] - + epsilon*dp_f_impl[i]; + p_f_ghost_out[idx1g(i)] = p_f + + ((1.0 - epsilon)*dp_f_dt_expl[i] + + epsilon*dp_f_dt_impl[i])*sim->dt; - /* apply relaxation */ - p_f_ghost_out[idx1g(i)] = p_f*(1.0 - theta) - + p_f_ghost_out[idx1g(i)]*theta; + /* apply relaxation */ + p_f_ghost_out[idx1g(i)] = p_f*(1.0 - theta) + + p_f_ghost_out[idx1g(i)]*theta; - r_norm[i] = (p_f_ghost_out[idx1g(i)] - p_f)/(p_f + 1e-16); - } + r_norm[i] = (p_f_ghost_out[idx1g(i)] - p_f)/(p_f + 1e-16); + } - r_norm_max = max(r_norm, sim->nz); + r_norm_max = max(r_norm, sim->nz); #ifdef DEBUG - puts(".. p_f_ghost_out:"); print_array(p_f_ghost_out, sim->nz+2); + puts(".. p_f_ghost_out:"); + print_array(p_f_ghost_out, sim->nz+2); #endif - copy_values(p_f_ghost_out, sim->p_f_ghost, sim->nz+2); + copy_values(p_f_ghost_out, sim->p_f_ghost, sim->nz+2); #ifdef DEBUG - puts(".. p_f_ghost after update:"); - print_array(sim->p_f_ghost, sim->nz+2); + puts(".. p_f_ghost after update:"); + print_array(sim->p_f_ghost, sim->nz+2); #endif - if (r_norm_max <= rel_tol) { - set_bc_dirichlet(sim->p_f_ghost, sim->nz, +1, p_f_top); - sim->p_f_ghost[idx1g(sim->nz-1)] = p_f_top; /* top node in BC */ - set_bc_neumann(sim->p_f_ghost, sim->nz, -1); - free(dp_f_expl); - free(dp_f_impl); - free(p_f_ghost_out); - free(r_norm); + if (r_norm_max <= rel_tol) { #ifdef DEBUG - printf(".. Solution converged after %d iterations\n", iter); + printf(".. Solution converged after %d iterations\n", iter); #endif - return 0; - } + solved = 1; + break; + } + } + free(dp_f_dt_impl); + free(p_f_ghost_out); + free(r_norm); } - free(dp_f_expl); - free(dp_f_impl); - free(p_f_ghost_out); - free(r_norm); - fprintf(stderr, "darcy_solver_1d: "); - fprintf(stderr, "Solution did not converge after %d iterations\n", iter); - fprintf(stderr, ".. Residual normalized error: %f\n", r_norm_max); - return 1; + if (!solved) { + fprintf(stderr, "darcy_solver_1d: "); + fprintf(stderr, "Solution did not converge after %d iterations\n", + iter); + fprintf(stderr, ".. Residual normalized error: %f\n", r_norm_max); + } + free(dp_f_dt_expl); + return solved - 1; }