tMove main loop into simulation.c - cngf-pf - continuum model for granular flows with pore-pressure dynamics (renamed from 1d_fd_simple_shear)
HTML git clone git://src.adamsgaard.dk/cngf-pf
DIR Log
DIR Files
DIR Refs
DIR README
DIR LICENSE
---
DIR commit f74ed98b4c7ea4bde059ced6673690a20f85e016
DIR parent 699306346b7a372e537e8b3edd7598c7ec419411
HTML Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Mon, 6 Apr 2020 13:50:25 +0200
Move main loop into simulation.c
Diffstat:
M 1d_fd_simple_shear.c | 102 +++----------------------------
M fluid.c | 2 +-
M simulation.c | 109 +++++++++++++++++++++++++++++++
M simulation.h | 2 ++
4 files changed, 119 insertions(+), 96 deletions(-)
---
DIR diff --git a/1d_fd_simple_shear.c b/1d_fd_simple_shear.c
t@@ -12,16 +12,8 @@
#include "parameter_defaults.h"
/* relative tolerance criteria for granular fluidity solver */
-#define RTOL_GRANULAR 1e-5
-#define MAX_ITER_GRANULAR 10000
-
-/* relative tolerance criteria for fluid-pressure solver */
-#define RTOL_DARCY 1e-5
-#define MAX_ITER_DARCY 10000
-
-/* relative tolerance criteria when shear velocity is restricted */
-#define RTOL_STRESS 1e-3
-#define MAX_ITER_STRESS 20000
+#define RTOL 1e-5
+#define MAX_ITER_1D_FD_SIMPLE_SHEAR 10000
/* uncomment to print time spent per time step to stdout */
/*#define BENCHMARK_PERFORMANCE*/
t@@ -77,8 +69,8 @@ int
main(int argc, char* argv[])
{
int i, normalize;
- unsigned long iter, stressiter;
- double new_phi, new_k, filetimeclock, res_norm, mu_wall_orig;
+ unsigned long iter;
+ double new_phi, new_k, filetimeclock, mu_wall_orig;
#ifdef BENCHMARK_PERFORMANCE
clock_t t_begin, t_end;
double t_elapsed;
t@@ -87,7 +79,7 @@ main(int argc, char* argv[])
#ifdef __OpenBSD__
if (pledge("stdio wpath cpath", NULL) == -1) {
fprintf(stderr, "error: pledge failed");
- exit(1);
+ exit(2);
}
#endif
t@@ -255,94 +247,14 @@ main(int argc, char* argv[])
filetimeclock = 0.0;
iter = 0;
mu_wall_orig = sim.mu_wall;
- res_norm = NAN;
while (sim.t <= sim.t_end) {
#ifdef BENCHMARK_PERFORMANCE
t_begin = clock();
#endif
- stressiter = 0;
- do {
- printf("\niter %lu, sim.t = %g s\n", iter, sim.t);
-
- if (sim.transient) {
- /* step 1 */
- compute_inertia_number(); /* Eq. 1 */
- for (i=0; i<sim.nz; i++)
- printf("sim.I[%d] = %g\n", i, sim.I[i]);
- /* step 2 */
- compute_critical_state_porosity(); /* Eq. 2 */
- for (i=0; i<sim.nz; i++)
- printf("sim.phi_c[%d] = %g\n", i, sim.phi_c[i]);
- /* step 3 */
- compute_tan_dilatancy_angle(); /* Eq. 5 */
- for (i=0; i<sim.nz; i++)
- printf("sim.tan_psi[%d] = %g\n", i, sim.tan_psi[i]);
- }
- compute_critical_state_friction(); /* Eq. 7 */
- for (i=0; i<sim.nz; i++)
- printf("sim.mu_c[%d] = %g\n", i, sim.mu_c[i]);
-
- /* step 4 */
- if (sim.transient) {
- compute_porosity_change(); /* Eq. 3 */
- for (i=0; i<sim.nz; i++)
- printf("sim.phi_dot[%d] = %g\n", i, sim.phi_dot[i]);
- compute_permeability(); /* Eq. 6 */
- for (i=0; i<sim.nz; i++)
- printf("sim.k[%d] = %g\n", i, sim.k[i]);
- }
- compute_friction(); /* Eq. 4 */
- for (i=0; i<sim.nz; i++)
- printf("sim.mu[%d] = %g\n", i, sim.mu[i]);
-
- /* step 5, Eq. 13 */
- if (sim.fluid)
- if (darcy_solver_1d(MAX_ITER_DARCY, RTOL_DARCY))
- exit(1);
-
- /* step 6 */
- compute_effective_stress(); /* Eq. 9 */
-
- /* step 7 */
- compute_local_fluidity(); /* Eq. 10 */
- compute_cooperativity_length(); /* Eq. 12 */
-
- /* step 8, Eq. 11 */
- if (implicit_1d_jacobian_poisson_solver(MAX_ITER_GRANULAR,
- RTOL_GRANULAR))
- exit(1);
-
- /* step 9 */
- compute_shear_strain_rate_plastic(); /* Eq. 8 */
- compute_inertia_number();
- compute_shear_velocity();
-
- if (!isnan(sim.v_x_limit) || !isnan(sim.v_x_fix)) {
- if (!isnan(sim.v_x_limit)) {
- res_norm = (sim.v_x_limit - sim.v_x[sim.nz-1])
- /(sim.v_x[sim.nz-1] + 1e-12);
- if (res_norm > 0.0)
- res_norm = 0.0;
- } else {
- res_norm = (sim.v_x_fix - sim.v_x[sim.nz-1])
- /(sim.v_x[sim.nz-1] + 1e-12);
- }
- sim.mu_wall *= 1.0 + (res_norm*1e-2);
- }
-
- if (++stressiter > MAX_ITER_STRESS) {
- fprintf(stderr, "error: stress solution did not converge:\n");
- fprintf(stderr, "v_x=%g, v_x_fix=%g, v_x_limit=%g, "
- "res_norm=%g, mu_wall=%g\n",
- sim.v_x[sim.nz-1], sim.v_x_fix, sim.v_x_limit,
- res_norm, sim.mu_wall);
- return 10;
- }
-
- } while ((!isnan(sim.v_x_fix) || !isnan(sim.v_x_limit))
- && fabs(res_norm) > RTOL_STRESS);
+ if (coupled_shear_solver(MAX_ITER_1D_FD_SIMPLE_SHEAR, RTOL))
+ exit(10);
#ifdef BENCHMARK_PERFORMANCE
t_end = clock();
DIR diff --git a/fluid.c b/fluid.c
t@@ -32,7 +32,7 @@ set_largest_fluid_timestep(const double safety)
fprintf(stderr, "error: cell spacing negative (%g) in cell %d\n",
dx[i], i);
free(dx);
- exit(10);
+ exit(1);
}
if (dx[i] < dx_min) dx_min = dx[i];
}
DIR diff --git a/simulation.c b/simulation.c
t@@ -3,9 +3,22 @@
#include <math.h>
#include "arrays.h"
#include "simulation.h"
+#include "fluid.h"
/* #define SHOW_PARAMETERS */
+/* relative tolerance criteria for granular fluidity solver */
+#define RTOL_GRANULAR 1e-5
+#define MAX_ITER_GRANULAR 10000
+
+/* relative tolerance criteria for fluid-pressure solver */
+#define RTOL_DARCY 1e-5
+#define MAX_ITER_DARCY 10000
+
+/* relative tolerance criteria when shear velocity is restricted */
+#define RTOL_STRESS 1e-3
+#define MAX_ITER_STRESS 20000
+
void
prepare_arrays()
{
t@@ -621,3 +634,99 @@ print_output(FILE* fp, const int norm)
free(v_x_out);
}
+
+int
+coupled_shear_solver(const int max_iter,
+ const double rel_tol)
+{
+ int coupled_iter, stress_iter;
+
+ double res_norm;
+ res_norm = NAN;
+
+ stress_iter = 0;
+ do { /* stress iterations */
+
+ coupled_iter = 0;
+ do { /* coupled iterations */
+
+ if (sim.transient) {
+ /* step 1 */
+ compute_inertia_number(); /* Eq. 1 */
+ /* step 2 */
+ compute_critical_state_porosity(); /* Eq. 2 */
+ /* step 3 */
+ compute_tan_dilatancy_angle(); /* Eq. 5 */
+ }
+ compute_critical_state_friction(); /* Eq. 7 */
+
+ /* step 4 */
+ if (sim.transient) {
+ compute_porosity_change(); /* Eq. 3 */
+ compute_permeability(); /* Eq. 6 */
+ }
+ compute_friction(); /* Eq. 4 */
+
+ /* step 5, Eq. 13 */
+ if (sim.fluid)
+ if (darcy_solver_1d(MAX_ITER_DARCY, RTOL_DARCY))
+ exit(11);
+
+ /* step 6 */
+ compute_effective_stress(); /* Eq. 9 */
+
+ /* step 7 */
+ compute_local_fluidity(); /* Eq. 10 */
+ compute_cooperativity_length(); /* Eq. 12 */
+
+ /* step 8, Eq. 11 */
+ if (implicit_1d_jacobian_poisson_solver(MAX_ITER_GRANULAR,
+ RTOL_GRANULAR))
+ exit(12);
+
+ /* step 9 */
+ compute_shear_strain_rate_plastic(); /* Eq. 8 */
+ compute_inertia_number();
+ compute_shear_velocity();
+
+#ifdef DEBUG
+ for (i=0; i<sim.nz) {
+ printf("\nsim.t = %g s\n", sim.t);
+ printf("sim.I[%d] = %g\n", i, sim.I[i]);
+ printf("sim.phi_c[%d] = %g\n", i, sim.phi_c[i]);
+ printf("sim.tan_psi[%d] = %g\n", i, sim.tan_psi[i]);
+ printf("sim.mu_c[%d] = %g\n", i, sim.mu_c[i]);
+ printf("sim.phi_dot[%d] = %g\n", i, sim.phi_dot[i]);
+ printf("sim.k[%d] = %g\n", i, sim.k[i]);
+ printf("sim.mu[%d] = %g\n", i, sim.mu[i]);
+ }
+#endif
+
+ if (!isnan(sim.v_x_limit) || !isnan(sim.v_x_fix)) {
+ if (!isnan(sim.v_x_limit)) {
+ res_norm = (sim.v_x_limit - sim.v_x[sim.nz-1])
+ /(sim.v_x[sim.nz-1] + 1e-12);
+ if (res_norm > 0.0)
+ res_norm = 0.0;
+ } else {
+ res_norm = (sim.v_x_fix - sim.v_x[sim.nz-1])
+ /(sim.v_x[sim.nz-1] + 1e-12);
+ }
+ sim.mu_wall *= 1.0 + (res_norm*1e-2);
+ }
+
+ if (++stress_iter > MAX_ITER_STRESS) {
+ fprintf(stderr, "error: stress solution did not converge:\n");
+ fprintf(stderr, "v_x=%g, v_x_fix=%g, v_x_limit=%g, "
+ "res_norm=%g, mu_wall=%g\n",
+ sim.v_x[sim.nz-1], sim.v_x_fix, sim.v_x_limit,
+ res_norm, sim.mu_wall);
+ return 10;
+ }
+ } while (0);
+
+ } while ((!isnan(sim.v_x_fix) || !isnan(sim.v_x_limit))
+ && fabs(res_norm) > RTOL_STRESS);
+ return 0;
+}
+
DIR diff --git a/simulation.h b/simulation.h
t@@ -153,4 +153,6 @@ int implicit_1d_jacobian_poisson_solver(const int max_iter,
void write_output_file(const int normalize);
void print_output(FILE* fp, const int normalize);
+int coupled_shear_solver(const int max_iter, const double rel_tol);
+
#endif