tFirst untested implementation of fluid solver - cngf-pf - continuum model for … | |
git clone git://src.adamsgaard.dk/cngf-pf | |
Log | |
Files | |
Refs | |
README | |
LICENSE | |
--- | |
commit d6d5abd9c90f013b39ad99462f4bdf419562942a | |
parent 0a594e0ee9565c0f780e85e1af04e76d6c4b92cf | |
Author: Anders Damsgaard <[email protected]> | |
Date: Tue, 16 Apr 2019 10:36:21 +0200 | |
First untested implementation of fluid solver | |
Diffstat: | |
A fluid.c | 131 +++++++++++++++++++++++++++++… | |
A fluid.h | 11 +++++++++++ | |
M main.c | 9 +++++++-- | |
M simulation.c | 2 +- | |
4 files changed, 150 insertions(+), 3 deletions(-) | |
--- | |
diff --git a/fluid.c b/fluid.c | |
t@@ -0,0 +1,131 @@ | |
+#include <stdlib.h> | |
+#include "simulation.h" | |
+#include "arrays.h" | |
+ | |
+static double darcy_pressure_change_1d( | |
+ const int i, | |
+ const int nz, | |
+ const double* p_f_ghost_in, | |
+ const double* phi, | |
+ const double* k, | |
+ const double dz, | |
+ const double dt, | |
+ const double beta_f, | |
+ const double mu_f) | |
+{ | |
+ const double p = p_f_ghost_in[idx1g(i)]; | |
+ double p_zn = p_f_ghost_in[idx1g(i-1)]; | |
+ double p_zp = p_f_ghost_in[idx1g(i+1)]; | |
+ | |
+ const double k_ = k[i]; | |
+ const double k_zn = k[i-1]; | |
+ const double k_zp = k[i+1]; | |
+ | |
+ double div_k_grad_p = | |
+ (2.*k_zp*k_/(k_zp + k_) * | |
+ (p_zp - p)/dz | |
+ - | |
+ 2.*k_zn*k_/(k_zn + k_) * | |
+ (p - p_zn)/dz)/dz; | |
+ | |
+ /* return delta p */ | |
+ return dt/(beta_f*phi[i]*mu_f)*div_k_grad_p; | |
+} | |
+ | |
+int darcy_solver_1d( | |
+ struct simulation* sim, | |
+ const int max_iter, | |
+ const double rel_tol) | |
+{ | |
+ | |
+ /* compute explicit solution to pressure change */ | |
+ double* dp_f_expl = zeros(sim->nz); | |
+ for (int 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); | |
+ | |
+ /* choose integration method, parameter in [0.0; 1.0] | |
+ * epsilon = 0.0: explicit | |
+ * epsilon = 0.5: Crank-Nicolson | |
+ * epsilon = 1.0: implicit */ | |
+ const double 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: overrrelaxation */ | |
+ const double theta = 0.05; | |
+ | |
+ double p_f; | |
+ | |
+ /* compute implicit solution to pressure change */ | |
+ int iter; | |
+ double* dp_f_impl = zeros(sim->nz); | |
+ double* p_f_ghost_out = zeros(sim->nz+2); | |
+ double* r_norm = zeros(sim->nz); | |
+ double r_norm_max = NAN; | |
+ for (iter=0; iter<max_iter; ++iter) { | |
+ | |
+ /* Dirichlet BCs resemble fixed particle velocities */ | |
+ set_bc_dirichlet(sim->p_f_ghost, sim->nz, +1, 0.0); | |
+ | |
+ /* Neumann BCs resemble free surfaces */ | |
+ set_bc_neumann(sim->p_f_ghost, sim->nz, -1); | |
+ | |
+ for (int i=0; i<sim->nz; ++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 (int i=0; i<sim->nz; ++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]; | |
+ | |
+ /* 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_max = max(r_norm, sim->nz); | |
+ | |
+ copy_values(p_f_ghost_out, sim->p_f_ghost, sim->nz+2); | |
+ | |
+ if (r_norm_max <= rel_tol) { | |
+ set_bc_dirichlet(sim->p_f_ghost, sim->nz, +1, 0.0); | |
+ 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); | |
+ printf(".. Solution converged after %d iterations\n", iter); | |
+ return 0; | |
+ } | |
+ } | |
+ | |
+ 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; | |
+} | |
diff --git a/fluid.h b/fluid.h | |
t@@ -0,0 +1,11 @@ | |
+#ifndef FLUID_H_ | |
+#define FLUID_H_ | |
+ | |
+#include "simulation.h" | |
+ | |
+int darcy_solver_1d( | |
+ struct simulation* sim, | |
+ const int max_iter, | |
+ const double rel_tol); | |
+ | |
+#endif | |
diff --git a/main.c b/main.c | |
t@@ -4,8 +4,9 @@ | |
#include <getopt.h> | |
#include "simulation.h" | |
+#include "fluid.h" | |
-#define VERSION "0.1" | |
+#define VERSION "0.2" | |
#define PROGNAME "1d_fd_simple_shear" | |
#include "parameter_defaults.h" | |
t@@ -200,8 +201,12 @@ int main(int argc, char* argv[]) | |
while (sim.t <= sim.t_end) { | |
init_normal_stress(&sim); | |
- if (sim.fluid) | |
+ if (sim.fluid) { | |
init_water_pressure(&sim); | |
+ if (darcy_solver_1d(&sim, 10000, 1e-5)) | |
+ exit(1); | |
+ } | |
+ | |
compute_effective_stress(&sim); | |
init_friction(&sim); | |
compute_cooperativity_length(&sim); | |
diff --git a/simulation.c b/simulation.c | |
t@@ -371,7 +371,7 @@ int implicit_1d_jacobian_poisson_solver( | |
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); | |
+ 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); */ |