commit e94f455468778b12e109147a3c3eb2732b554a9a
parent 464550e2aba5f3dcec7454d023843dae46aafc24
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Mon, 1 Jul 2019 11:10:36 +0200
Fix conditional operations depending on epsilon value
Diffstat:
| M | fluid.c | | | 69 | +++++++++++++++++++++++++++++++++++++++++++-------------------------- |
1 file changed, 43 insertions(+), 26 deletions(-)
diff --git a/fluid.c b/fluid.c
@@ -112,7 +112,8 @@ darcy_solver_1d(struct simulation* sim,
* epsilon = 0.5: Crank-Nicolson
* epsilon = 1.0: implicit */
/* epsilon = 0.5; */
- epsilon = 0.0;
+ /* epsilon = 0.0; */
+ epsilon = 1.0;
/* choose relaxation factor, parameter in ]0.0; 1.0]
* theta in ]0.0; 1.0]: underrelaxation
@@ -132,21 +133,23 @@ darcy_solver_1d(struct simulation* sim,
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);
- /* compute explicit solution to pressure change */
- solved = 1;
- dp_f_dt_expl = zeros(sim->nz);
- for (i=0; i<sim->nz; ++i)
- 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;
+ solved = 0;
+
+ if (epsilon < 1.0) {
+ /* compute explicit solution to pressure change */
+ dp_f_dt_expl = zeros(sim->nz);
+ for (i=0; i<sim->nz; ++i)
+ 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 > 0.0) {
/* compute implicit solution to pressure change */
dp_f_dt_impl = zeros(sim->nz);
p_f_ghost_out = zeros(sim->nz+2);
@@ -175,14 +178,17 @@ darcy_solver_1d(struct simulation* sim,
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]);
+ i, dp_f_dt_expl[i], i, dp_f_dt_impl[i]);
#endif
p_f = sim->p_f_ghost[idx1g(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;
+ + epsilon*dp_f_dt_impl[i]*sim->dt;
+
+ if (epsilon < 1.0)
+ p_f_ghost_out[idx1g(i)] += (1.0 - epsilon)
+ *dp_f_dt_expl[i]*sim->dt;
/* apply relaxation */
p_f_ghost_out[idx1g(i)] = p_f*(1.0 - theta)
@@ -214,14 +220,25 @@ darcy_solver_1d(struct simulation* sim,
free(dp_f_dt_impl);
free(p_f_ghost_out);
free(r_norm);
+ 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);
+ }
+ } else {
+ for (i=0; i<sim->nz; ++i)
+ sim->p_f_ghost[idx1g(i)] += dp_f_dt_expl[i]*sim->dt;
+ solved = 1;
+#ifdef DEBUG
+ puts(".. dp_f_dt_expl:");
+ print_array(dp_f_dt_expl, sim->nz);
+ puts(".. p_f_ghost after explicit solution:");
+ print_array(sim->p_f_ghost, sim->nz+2);
+#endif
}
- 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);
+ if (epsilon < 1.0)
+ free(dp_f_dt_expl);
return solved - 1;
}
