tperform temporal integration of fluid pressure outside of fluid solver - cngf-… | |
git clone git://src.adamsgaard.dk/cngf-pf | |
Log | |
Files | |
Refs | |
README | |
LICENSE | |
--- | |
commit 97e88f765d0b872280c99f2cd5753ed1afefe420 | |
parent 3dfa0a14ecc48bdd9021dbe83072c161af72cbf8 | |
Author: Anders Damsgaard <[email protected]> | |
Date: Wed, 18 Nov 2020 11:31:56 +0100 | |
perform temporal integration of fluid pressure outside of fluid solver | |
Diffstat: | |
M fluid.c | 82 +++++++++++++----------------… | |
M simulation.c | 16 ++++++++-------- | |
2 files changed, 41 insertions(+), 57 deletions(-) | |
--- | |
diff --git a/fluid.c b/fluid.c | |
t@@ -98,11 +98,11 @@ square_pulse(const double time, | |
} | |
static void | |
-set_fluid_bcs(struct simulation *sim, const double p_f_top) | |
+set_fluid_bcs(double *p_f_ghost, struct simulation *sim, const double p_f_top) | |
{ | |
- set_bc_dirichlet(sim->p_f_ghost, sim->nz, +1, p_f_top); | |
- sim->p_f_ghost[sim->nz] = p_f_top; /* Include top node in BC */ | |
- set_bc_neumann(sim->p_f_ghost, | |
+ set_bc_dirichlet(p_f_ghost, sim->nz, +1, p_f_top); | |
+ p_f_ghost[sim->nz] = p_f_top; /* Include top node in BC */ | |
+ set_bc_neumann(p_f_ghost, | |
sim->nz, | |
-1, | |
sim->phi[0] * sim->rho_f * sim->G, | |
t@@ -170,9 +170,13 @@ darcy_solver_1d(struct simulation *sim, | |
const double rel_tol) | |
{ | |
int i, iter, solved = 0; | |
- double epsilon, theta, p_f_top, r_norm_max = NAN; | |
- double *dp_f_dt_expl, *p_f_ghost_old, *dp_f_dt_impl, *p_f_ghost_new; | |
- double *tmp, *r_norm; | |
+ 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); | |
/* choose integration method, parameter in [0.0; 1.0] | |
* epsilon = 0.0: explicit | |
t@@ -180,12 +184,6 @@ darcy_solver_1d(struct simulation *sim, | |
* epsilon = 1.0: * implicit */ | |
epsilon = 0.5; | |
- /* choose relaxation factor, parameter in ]0.0; 1.0] theta in | |
- * ]0.0; 1.0]: underrelaxation theta = 1.0: Gauss-Seidel theta > | |
- * 1.0: overrelaxation */ | |
- /* TODO: values other than 1.0 do not work! */ | |
- theta = 1.0; | |
- | |
if (isnan(sim->p_f_mod_pulse_time)) | |
p_f_top = sim->p_f_top | |
+ sine_wave(sim->t, | |
t@@ -205,14 +203,15 @@ darcy_solver_1d(struct simulation *sim, | |
sim->p_f_mod_freq, | |
sim->p_f_mod_pulse_time); | |
+ | |
/* set fluid BCs (1 of 2) */ | |
- set_fluid_bcs(sim, p_f_top); | |
+ copy_values(sim->p_f_ghost, p_f_ghost_new, sim->nz + 2); | |
+ set_fluid_bcs(p_f_ghost_new, sim, p_f_top); | |
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, | |
+ p_f_dot_expl[i] = darcy_pressure_change_1d(i, | |
sim->nz, | |
sim->p_f_gh… | |
sim->phi, | |
t@@ -226,25 +225,23 @@ darcy_solver_1d(struct simulation *sim, | |
} | |
if (epsilon > 0.0) { | |
/* compute implicit solution with Jacobian iterations */ | |
- dp_f_dt_impl = zeros(sim->nz); | |
- p_f_ghost_new = zeros(sim->nz + 2); | |
r_norm = zeros(sim->nz); | |
- p_f_ghost_old = empty(sim->nz + 2); | |
- copy_values(sim->p_f_ghost, p_f_ghost_old, sim->nz + 2); | |
+ old_val = empty(sim->nz); | |
for (iter = 0; iter < max_iter; ++iter) { | |
+ copy_values(sim->p_f_dot, old_val, sim->nz); | |
/* set fluid BCs (2 of 2) */ | |
- set_fluid_bcs(sim, p_f_top); | |
+ set_fluid_bcs(p_f_ghost_new, sim, p_f_top); | |
#ifdef DEBUG | |
puts(".. p_f_ghost after BC:"); | |
- print_array(sim->p_f_ghost, sim->nz + 2); | |
+ print_array(p_f_ghost_new, sim->nz + 2); | |
#endif | |
for (i = 0; i < sim->nz - 1; ++i) | |
- dp_f_dt_impl[i] = darcy_pressure_change_1d(i, | |
+ p_f_dot_impl[i] = darcy_pressure_change_1d(i, | |
sim… | |
- sim… | |
+ p_f… | |
sim… | |
sim… | |
sim… | |
t@@ -255,22 +252,13 @@ 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… | |
- i, dp_f_dt_expl[i], i, dp_f_dt_impl[i]); | |
-#endif | |
- p_f_ghost_new[i + 1] = p_f_ghost_old[i + 1] | |
- + epsilon * dp_f_dt_impl[i] * sim->dt; | |
+ sim->p_f_dot[i] = epsilon * p_f_dot_impl[i]; | |
if (epsilon < 1.0) | |
- p_f_ghost_new[i + 1] += (1.0 - epsilon) | |
- * dp_f_dt_expl[i] * sim->dt; | |
- | |
- p_f_ghost_new[i + 1] = p_f_ghost_old[i + 1] * … | |
- + p_f_ghost_new[i + 1] * theta; | |
+ sim->p_f_dot[i + 1] += (1.0 - epsilon)… | |
- r_norm[i] = fabs(residual(p_f_ghost_new[i + 1], | |
- sim->p_f_ghost[i + 1… | |
+ p_f_ghost_new[i + 1] += sim->p_f_dot[i] * sim-… | |
+ r_norm[i] = fabs(residual(sim->p_f_dot[i], old… | |
} | |
r_norm_max = max(r_norm, sim->nz); | |
t@@ -284,7 +272,7 @@ darcy_solver_1d(struct simulation *sim, | |
sim->p_f_ghost = tmp; | |
#ifdef DEBUG | |
puts(".. p_f_ghost after update:"); | |
- print_array(sim->p_f_ghost, sim->nz + 2); | |
+ print_array(p_f_ghost_new, sim->nz + 2); | |
#endif | |
if (r_norm_max <= rel_tol) { | |
t@@ -295,10 +283,10 @@ darcy_solver_1d(struct simulation *sim, | |
break; | |
} | |
} | |
- free(dp_f_dt_impl); | |
+ free(p_f_dot_impl); | |
+ free(p_f_dot_expl); | |
free(p_f_ghost_new); | |
free(r_norm); | |
- free(p_f_ghost_old); | |
if (!solved) { | |
fprintf(stderr, "darcy_solver_1d: "); | |
fprintf(stderr, "Solution did not converge after %d it… | |
t@@ -307,19 +295,15 @@ darcy_solver_1d(struct simulation *sim, | |
} | |
} else { | |
for (i = 0; i < sim->nz; ++i) | |
- sim->p_f_ghost[i + 1] += dp_f_dt_expl[i] * sim->dt; | |
+ sim->p_f_dot[i] = p_f_dot_expl[i]; | |
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); | |
+ puts(".. p_f_dot_expl:"); | |
+ print_array(p_f_dot_expl, sim->nz); | |
#endif | |
+ free(p_f_dot_expl); | |
} | |
- set_fluid_bcs(sim, p_f_top); | |
- | |
- if (epsilon < 1.0) | |
- free(dp_f_dt_expl); | |
+ set_fluid_bcs(sim->p_f_ghost, sim, p_f_top); | |
return solved - 1; | |
} | |
diff --git a/simulation.c b/simulation.c | |
t@@ -744,9 +744,9 @@ temporal_increment(struct simulation *sim) | |
for (i = 0; i < sim->nz; ++i) | |
sim->phi[i] += sim->phi_dot[i] * sim->dt; | |
- /*if (sim->fluid) | |
+ if (sim->fluid) | |
for (i = 1; i <= sim->nz; ++i) | |
- sim->p_f_ghost[i] += sim->p_f_dot[i] * sim->dt;*/ | |
+ sim->p_f_ghost[i] += sim->p_f_dot[i] * sim->dt; | |
sim->t += sim->dt; | |
} | |
t@@ -757,12 +757,12 @@ coupled_shear_solver(struct simulation *sim, | |
const double rel_tol) | |
{ | |
int i, coupled_iter, stress_iter = 0; | |
- double *r_norm, *oldval; | |
+ 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); | |
- oldval = empty(sim->nz); | |
+ old_val = empty(sim->nz); | |
} | |
do { /* stress iterations */ | |
t@@ -770,7 +770,7 @@ coupled_shear_solver(struct simulation *sim, | |
do { /* coupled iterations */ | |
if (sim->transient) { | |
- copy_values(sim->phi_dot, oldval, sim->nz); | |
+ copy_values(sim->phi_dot, old_val, sim->nz); | |
/* step 1 */ | |
compute_inertia_number(sim); /* Eq. 1 */ | |
t@@ -828,14 +828,14 @@ coupled_shear_solver(struct simulation *sim, | |
/* stable porosity change field == coupled solution co… | |
if (sim->transient) { | |
for (i = 0; i < sim->nz; ++i) | |
- r_norm[i] = fabs(residual(sim->phi_dot… | |
+ r_norm[i] = fabs(residual(sim->phi_dot… | |
r_norm_max = max(r_norm, sim->nz); | |
if (r_norm_max <= rel_tol) | |
break; | |
if (coupled_iter++ >= max_iter) { | |
free(r_norm); | |
- free(oldval); | |
+ free(old_val); | |
fprintf(stderr, "coupled_shear_solver:… | |
fprintf(stderr, "Transient solution di… | |
"after %d iterations\n", coupl… | |
t@@ -874,7 +874,7 @@ coupled_shear_solver(struct simulation *sim, | |
if (sim->transient) { | |
free(r_norm); | |
- free(oldval); | |
+ free(old_val); | |
} | |
temporal_increment(sim); |