commit bfb86ef9008e86da090c6d6c84e165517d0a85b0
parent 7fe3756e69502a715e4554615f16304d6ef76565
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Wed, 18 Nov 2020 13:55:10 +0100
preallocate arrays for iterative solvers
Diffstat:
| M | fluid.c | | | 135 | ++++++++++++++++++++++++++++++++++++++++--------------------------------------- |
| M | simulation.c | | | 97 | ++++++++++++++++++++++++++++++++++++++----------------------------------------- |
| M | simulation.h | | | 6 | ++++++ |
3 files changed, 121 insertions(+), 117 deletions(-)
diff --git a/fluid.c b/fluid.c
@@ -111,21 +111,24 @@ set_fluid_bcs(double *p_f_ghost, struct simulation *sim, const double p_f_top)
static double
darcy_pressure_change_1d(const int i,
- const int nz,
- const double *p_f_ghost_in,
- const double *phi,
- const double *phi_dot,
- const double *k,
- const double dz,
- const double beta_f,
- const double alpha,
- const double mu_f,
- const double D)
+ const int nz,
+ const double *p_f_ghost_in,
+ const double *phi,
+ const double *phi_dot,
+ const double *k,
+ const double dz,
+ const double beta_f,
+ const double alpha,
+ const double mu_f,
+ const double D)
{
double k_, div_k_grad_p, k_zn, k_zp;
if (D > 0.0)
- return D*(p_f_ghost_in[i + 2] - 2.0 * p_f_ghost_in[i + 1] + p_f_ghost_in[i])/(dz*dz);
+ return D * (p_f_ghost_in[i + 2]
+ - 2.0 * p_f_ghost_in[i + 1]
+ + p_f_ghost_in[i])
+ / (dz * dz);
else {
@@ -171,17 +174,13 @@ darcy_solver_1d(struct simulation *sim,
{
int i, iter, solved = 0;
double epsilon, p_f_top, r_norm_max = NAN;
- double *p_f_dot_expl, *p_f_dot_impl, *p_f_ghost_new;
- double *tmp, *r_norm, *old_val;
-
- p_f_dot_impl = zeros(sim->nz);
- p_f_dot_expl = zeros(sim->nz);
- p_f_ghost_new = zeros(sim->nz + 2);
+ double *tmp;
/* choose integration method, parameter in [0.0; 1.0]
* epsilon = 0.0: explicit
* epsilon = 0.5: Crank-Nicolson
- * epsilon = 1.0: * implicit */
+ * epsilon = 1.0: implicit */
+ /* epsilon = 0.5; */
epsilon = 0.5;
if (isnan(sim->p_f_mod_pulse_time))
@@ -205,84 +204,80 @@ darcy_solver_1d(struct simulation *sim,
/* set fluid BCs (1 of 2) */
- copy_values(sim->p_f_ghost, p_f_ghost_new, sim->nz + 2);
- set_fluid_bcs(p_f_ghost_new, sim, p_f_top);
+ set_fluid_bcs(sim->p_f_ghost, sim, p_f_top);
+ copy_values(sim->p_f_ghost, sim->new_ghost, sim->nz + 2);
+ /* explicit solution to pressure change */
if (epsilon < 1.0) {
- /* compute explicit solution to pressure change */
for (i = 0; i < sim->nz; ++i)
- p_f_dot_expl[i] = darcy_pressure_change_1d(i,
- sim->nz,
- sim->p_f_ghost,
- sim->phi,
- sim->phi_dot,
- sim->k,
- sim->dz,
- sim->beta_f,
- sim->alpha,
- sim->mu_f,
- sim->D);
+ sim->p_f_dot_expl[i] = darcy_pressure_change_1d(i,
+ sim->nz,
+ sim->p_f_ghost,
+ sim->phi,
+ sim->phi_dot,
+ sim->k,
+ sim->dz,
+ sim->beta_f,
+ sim->alpha,
+ sim->mu_f,
+ sim->D);
}
+
+ /* implicit solution with Jacobian iterations */
if (epsilon > 0.0) {
- /* compute implicit solution with Jacobian iterations */
- r_norm = zeros(sim->nz);
- old_val = empty(sim->nz);
for (iter = 0; iter < max_iter; ++iter) {
- copy_values(sim->p_f_dot, old_val, sim->nz);
+ copy_values(sim->p_f_dot_impl, sim->old_val, sim->nz);
/* set fluid BCs (2 of 2) */
- set_fluid_bcs(p_f_ghost_new, sim, p_f_top);
+ set_fluid_bcs(sim->new_ghost, sim, p_f_top);
#ifdef DEBUG
puts(".. p_f_ghost after BC:");
- print_array(p_f_ghost_new, sim->nz + 2);
+ print_array(sim->new_ghost, sim->nz + 2);
#endif
for (i = 0; i < sim->nz - 1; ++i)
- p_f_dot_impl[i] = darcy_pressure_change_1d(i,
- sim->nz,
- p_f_ghost_new,
- sim->phi,
- sim->phi_dot,
- sim->k,
- sim->dz,
- sim->beta_f,
- sim->alpha,
- sim->mu_f,
- sim->D);
+ sim->p_f_dot_impl[i] = darcy_pressure_change_1d(i,
+ sim->nz,
+ sim->new_ghost,
+ sim->phi,
+ sim->phi_dot,
+ sim->k,
+ sim->dz,
+ sim->beta_f,
+ sim->alpha,
+ sim->mu_f,
+ sim->D);
for (i = 0; i < sim->nz - 1; ++i) {
- sim->p_f_dot[i] = epsilon * p_f_dot_impl[i] + (1.0 - epsilon) * p_f_dot_expl[i];
- p_f_ghost_new[i + 1] += sim->p_f_dot[i] * sim->dt;
- r_norm[i] = fabs(residual(sim->p_f_dot[i], old_val[i]));
+ sim->new_ghost[i + 1] = sim->p_f_ghost[i + 1]
+ + sim->p_f_dot_impl[i] * sim->dt;
+ sim->p_f_dot_impl_r_norm[i] = fabs(residual(
+ sim->p_f_dot_impl[i], sim->old_val[i]));
}
- r_norm_max = max(r_norm, sim->nz);
+ r_norm_max = max(sim->p_f_dot_impl_r_norm, sim->nz);
#ifdef DEBUG
puts(".. p_f_ghost_new:");
- print_array(p_f_ghost_new, sim->nz + 2);
+ print_array(sim->new_ghost, sim->nz + 2);
#endif
- tmp = p_f_ghost_new;
- p_f_ghost_new = sim->p_f_ghost;
+ tmp = sim->new_ghost;
+ sim->new_ghost = sim->p_f_ghost;
sim->p_f_ghost = tmp;
#ifdef DEBUG
puts(".. p_f_ghost after update:");
- print_array(p_f_ghost_new, sim->nz + 2);
+ print_array(sim->new_ghost, sim->nz + 2);
#endif
if (r_norm_max <= rel_tol) {
#ifdef DEBUG
- printf(".. Solution converged after %d iterations\n", iter);
+ printf(".. Iterative solution converged after %d iterations\n", iter);
#endif
solved = 1;
break;
}
}
- free(p_f_dot_impl);
- free(p_f_dot_expl);
- free(p_f_ghost_new);
- free(r_norm);
if (!solved) {
fprintf(stderr, "darcy_solver_1d: ");
fprintf(stderr, "Solution did not converge after %d iterations\n",
@@ -290,15 +285,21 @@ darcy_solver_1d(struct simulation *sim,
fprintf(stderr, ".. Residual normalized error: %f\n", r_norm_max);
}
} else {
- for (i = 0; i < sim->nz; ++i)
- sim->p_f_dot[i] = p_f_dot_expl[i];
solved = 1;
+ }
+
+ for (i = 0; i < sim->nz; ++i)
+ sim->p_f_dot[i] = epsilon * sim->p_f_dot_impl[i]
+ + (1.0 - epsilon) * sim->p_f_dot_expl[i];
+
#ifdef DEBUG
+ printf(".. epsilon = %s\n", epsilon);
+ puts(".. p_f_dot_impl:");
+ print_array(sim->p_f_dot_impl, sim->nz);
puts(".. p_f_dot_expl:");
- print_array(p_f_dot_expl, sim->nz);
+ print_array(sim->p_f_dot_expl, sim->nz);
#endif
- free(p_f_dot_expl);
- }
+
set_fluid_bcs(sim->p_f_ghost, sim, p_f_top);
return solved - 1;
diff --git a/simulation.c b/simulation.c
@@ -141,27 +141,33 @@ prepare_arrays(struct simulation *sim)
free(sim->phi);
free(sim->k);
- sim->z = linspace(sim->origo_z, /* spatial coordinates */
+ sim->z = linspace(sim->origo_z,
sim->origo_z + sim->L_z,
sim->nz);
- sim->dz = sim->z[1] - sim->z[0]; /* cell spacing */
- sim->mu = zeros(sim->nz); /* stress ratio */
- sim->mu_c = zeros(sim->nz); /* critical-state stress ratio */
- sim->sigma_n_eff = zeros(sim->nz); /* effective normal stress */
- sim->sigma_n = zeros(sim->nz); /* normal stess */
- sim->p_f_ghost = zeros(sim->nz + 2); /* fluid pressure with ghost nodes */
- sim->p_f_dot = zeros(sim->nz); /* fluid pressure change */
- sim->phi = zeros(sim->nz); /* porosity */
- sim->phi_c = zeros(sim->nz); /* critical-state porosity */
- sim->phi_dot = zeros(sim->nz); /* rate of porosity change */
- sim->k = zeros(sim->nz); /* permeability */
- sim->xi = zeros(sim->nz); /* cooperativity length */
- sim->gamma_dot_p = zeros(sim->nz); /* shear velocity */
- sim->v_x = zeros(sim->nz); /* shear velocity */
- sim->g_local = zeros(sim->nz); /* local fluidity */
- sim->g_ghost = zeros(sim->nz + 2); /* fluidity with ghost nodes */
- sim->I = zeros(sim->nz); /* inertia number */
- sim->tan_psi = zeros(sim->nz); /* tan(dilatancy_angle) */
+ sim->dz = sim->z[1] - sim->z[0];
+ sim->mu = zeros(sim->nz);
+ sim->mu_c = zeros(sim->nz);
+ sim->sigma_n_eff = zeros(sim->nz);
+ sim->sigma_n = zeros(sim->nz);
+ sim->p_f_ghost = zeros(sim->nz + 2);
+ sim->p_f_dot = zeros(sim->nz);
+ sim->p_f_dot_expl = zeros(sim->nz);
+ sim->p_f_dot_impl = zeros(sim->nz);
+ sim->p_f_dot_impl_r_norm = zeros(sim->nz);
+ sim->phi = zeros(sim->nz);
+ sim->phi_c = zeros(sim->nz);
+ sim->phi_dot = zeros(sim->nz);
+ sim->k = zeros(sim->nz);
+ sim->xi = zeros(sim->nz);
+ sim->gamma_dot_p = zeros(sim->nz);
+ sim->v_x = zeros(sim->nz);
+ sim->g_local = zeros(sim->nz);
+ sim->g_ghost = zeros(sim->nz + 2);
+ sim->g_r_norm = zeros(sim->nz);
+ sim->I = zeros(sim->nz);
+ sim->tan_psi = zeros(sim->nz);
+ sim->old_val = empty(sim->nz);
+ sim->new_ghost = empty(sim->nz + 2);
}
void
@@ -174,6 +180,9 @@ free_arrays(struct simulation *sim)
free(sim->sigma_n);
free(sim->p_f_ghost);
free(sim->p_f_dot);
+ free(sim->p_f_dot_expl);
+ free(sim->p_f_dot_impl);
+ free(sim->p_f_dot_impl_r_norm);
free(sim->k);
free(sim->phi);
free(sim->phi_c);
@@ -183,8 +192,11 @@ free_arrays(struct simulation *sim)
free(sim->v_x);
free(sim->g_local);
free(sim->g_ghost);
+ free(sim->g_r_norm);
free(sim->I);
free(sim->tan_psi);
+ free(sim->old_val);
+ free(sim->new_ghost);
}
static void
@@ -642,11 +654,11 @@ implicit_1d_jacobian_poisson_solver(struct simulation *sim,
{
int iter, i;
double r_norm_max = NAN;
- double *tmp, *g_ghost_out = empty(sim->nz + 2), *r_norm = empty(sim->nz);
+ double *tmp;
for (iter = 0; iter < max_iter; ++iter) {
#ifdef DEBUG
- printf("\n@@@ ITERATION %d @@@\n", iter);
+ printf("\n@@@ %s: ITERATION %d @@@\n", __func__, iter);
#endif
/* Dirichlet BCs resemble fixed particle velocities */
set_bc_dirichlet(sim->g_ghost, sim->nz, -1, 0.0);
@@ -659,31 +671,28 @@ implicit_1d_jacobian_poisson_solver(struct simulation *sim,
poisson_solver_1d_cell_update(i,
sim->g_ghost,
sim->g_local,
- g_ghost_out,
- r_norm,
+ sim->new_ghost,
+ sim->g_r_norm,
sim->dz,
sim->xi);
- r_norm_max = max(r_norm, sim->nz);
+ r_norm_max = max(sim->g_r_norm, sim->nz);
- tmp = g_ghost_out;
- g_ghost_out = sim->g_ghost;
+ tmp = sim->new_ghost;
+ sim->new_ghost = sim->g_ghost;
sim->g_ghost = tmp;
if (r_norm_max <= rel_tol) {
set_bc_dirichlet(sim->g_ghost, sim->nz, -1, 0.0);
set_bc_dirichlet(sim->g_ghost, sim->nz, +1, 0.0);
- free(g_ghost_out);
- free(r_norm);
- /* printf(".. Solution converged after %d
- * iterations\n", iter); */
+#ifdef DEBUG
+ printf(".. Solution converged after %d iterations\n", iter + 1);
+#endif
return 0;
}
}
- free(g_ghost_out);
- free(r_norm);
fprintf(stderr, "implicit_1d_jacobian_poisson_solver: ");
- fprintf(stderr, "Solution did not converge after %d iterations\n", iter);
+ fprintf(stderr, "Solution did not converge after %d iterations\n", iter + 1);
fprintf(stderr, ".. Residual normalized error: %f\n", r_norm_max);
return 1;
}
@@ -753,24 +762,19 @@ temporal_increment(struct simulation *sim)
int
coupled_shear_solver(struct simulation *sim,
- const int max_iter,
- const double rel_tol)
+ const int max_iter,
+ const double rel_tol)
{
int i, coupled_iter, stress_iter = 0;
- double *r_norm, *old_val;
double r_norm_max, vel_res_norm = NAN, mu_wall_orig = sim->mu_wall;
- if (sim->transient) {
- r_norm = empty(sim->nz);
- old_val = empty(sim->nz);
- }
do { /* stress iterations */
coupled_iter = 0;
do { /* coupled iterations */
if (sim->transient) {
- copy_values(sim->phi_dot, old_val, sim->nz);
+ copy_values(sim->phi_dot, sim->old_val, sim->nz);
/* step 1 */
compute_inertia_number(sim); /* Eq. 1 */
@@ -828,14 +832,12 @@ coupled_shear_solver(struct simulation *sim,
/* stable porosity change field == coupled solution converged */
if (sim->transient) {
for (i = 0; i < sim->nz; ++i)
- r_norm[i] = fabs(residual(sim->phi_dot[i], old_val[i]));
- r_norm_max = max(r_norm, sim->nz);
+ sim->g_r_norm[i] = fabs(residual(sim->phi_dot[i], sim->old_val[i]));
+ r_norm_max = max(sim->g_r_norm, sim->nz);
if (r_norm_max <= rel_tol)
break;
if (coupled_iter++ >= max_iter) {
- free(r_norm);
- free(old_val);
fprintf(stderr, "coupled_shear_solver: ");
fprintf(stderr, "Transient solution did not converge "
"after %d iterations\n", coupled_iter);
@@ -872,11 +874,6 @@ coupled_shear_solver(struct simulation *sim,
if (!isnan(sim->v_x_limit))
sim->mu_wall = mu_wall_orig;
- if (sim->transient) {
- free(r_norm);
- free(old_val);
- }
-
temporal_increment(sim);
return 0;
diff --git a/simulation.h b/simulation.h
@@ -107,6 +107,9 @@ struct simulation {
double *sigma_n; /* normal stress [Pa] */
double *p_f_ghost; /* fluid pressure [Pa] */
double *p_f_dot; /* fluid pressure change [Pa/s] */
+ double *p_f_dot_expl; /* fluid pressure change (explicit solution) [Pa/s] */
+ double *p_f_dot_impl; /* fluid pressure change (implicit solution) [Pa/s] */
+ double *p_f_dot_impl_r_norm; /* normalized residual fluid pressure change [-] */
double *k; /* hydraulic permeability [m^2] */
double *phi; /* porosity [-] */
double *phi_c; /* critical-state porosity [-] */
@@ -116,8 +119,11 @@ struct simulation {
double *v_x; /* shear velocity [m/s] */
double *g_local; /* local fluidity */
double *g_ghost; /* fluidity with ghost nodes */
+ double *g_r_norm; /* normalized residual of fluidity field */
double *I; /* inertia number [-] */
double *tan_psi; /* tan(dilatancy_angle) [-] */
+ double *old_val; /* temporary storage for iterative solvers */
+ double *new_ghost; /* temporary storage for iterative solvers */
};
void init_sim(struct simulation *sim);
