commit 2562762c17ca9c5f80b30c7016dca66bba884f0a
parent 78f850afa0d533bca82ae05a01988ccb4b08df8b
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Wed, 14 Jan 2026 08:39:38 +0100
refactor(simulation): replace Jacobi with SOR solver and add performance tracking
Replace Jacobi iteration with SOR (Successive Over-Relaxation, omega=1.2) for
faster convergence in the Poisson solver. Add solver_stats tracking for
benchmarking (iteration counts, timesteps, elapsed time).
Diffstat:
4 files changed, 87 insertions(+), 33 deletions(-)
diff --git a/fluid.c b/fluid.c
@@ -3,6 +3,7 @@
#include <err.h>
#include "simulation.h"
#include "arrays.h"
+#include "fluid.h"
void
hydrostatic_fluid_pressure_distribution(struct simulation *sim)
@@ -372,6 +373,7 @@ darcy_solver_1d(struct simulation *sim,
#ifdef DEBUG
printf(".. Iterative solution converged after %d iterations\n", iter);
#endif
+ add_darcy_iters(iter + 1);
solved = 1;
break;
}
@@ -383,6 +385,7 @@ darcy_solver_1d(struct simulation *sim,
fprintf(stderr, "Solution did not converge after %d iterations\n",
iter);
fprintf(stderr, ".. Residual normalized error: %f\n", r_norm_max);
+ add_darcy_iters(iter + 1);
}
} else
solved = 1;
diff --git a/fluid.h b/fluid.h
@@ -13,4 +13,6 @@ int darcy_solver_1d(struct simulation *sim,
const int max_iter,
const double rel_tol);
+void add_darcy_iters(int iters);
+
#endif
diff --git a/simulation.c b/simulation.c
@@ -14,6 +14,17 @@
/* tolerance criteria when in velocity driven or velocity limited mode */
#define RTOL_VELOCITY 1e-3
+/* solver statistics for benchmarking */
+struct solver_stats {
+ long poisson_iters;
+ long darcy_iters;
+ long coupled_iters;
+ long stress_iters;
+ long timesteps;
+};
+
+static struct solver_stats g_stats = {0, 0, 0, 0, 0};
+
/* lower limit for effective normal stress sigma_n_eff for granular solver */
#define SIGMA_N_EFF_MIN 1e-1
@@ -628,38 +639,45 @@ residual(double new_val, double old_val)
}
static void
-poisson_solver_1d_cell_update(int i,
- const double *g_in,
- const double *g_local,
- double *g_out,
- double *r_norm,
- const double dz,
- const double *xi)
-{
- double coorp_term = dz * dz / (2.0 * pow(xi[i], 2.0));
-
- g_out[i + 1] = 1.0 / (1.0 + coorp_term)
- * (coorp_term * g_local[i] + g_in[i + 2] / 2.0
- + g_in[i] / 2.0);
- r_norm[i] = fabs(residual(g_out[i + 1], g_in[i + 1]));
+poisson_solver_1d_sor_cell_update(int i,
+ double *g,
+ const double *g_local,
+ double *r_norm,
+ const double dz,
+ const double *xi,
+ const double omega)
+{
+ double coorp_term = dz * dz / (2.0 * xi[i] * xi[i]);
+ double g_old = g[i + 1];
+
+ /* Gauss-Seidel update (uses already-updated g[i] from this iteration) */
+ double g_gs = 1.0 / (1.0 + coorp_term)
+ * (coorp_term * g_local[i] + g[i + 2] / 2.0
+ + g[i] / 2.0);
+
+ /* SOR relaxation */
+ g[i + 1] = (1.0 - omega) * g_old + omega * g_gs;
+ r_norm[i] = fabs(residual(g[i + 1], g_old));
#ifdef DEBUG
printf("-- %d --------------\n", i);
printf("coorp_term: %g\n", coorp_term);
printf(" g_local: %g\n", g_local[i]);
- printf(" g_in: %g\n", g_in[i + 1]);
- printf(" g_out: %g\n", g_out[i + 1]);
+ printf(" g_old: %g\n", g_old);
+ printf(" g_gs: %g\n", g_gs);
+ printf(" g_new: %g\n", g[i + 1]);
printf(" r_norm: %g\n", r_norm[i]);
#endif
}
static int
-implicit_1d_jacobian_poisson_solver(struct simulation *sim,
- const int max_iter,
- const double rel_tol)
+implicit_1d_sor_poisson_solver(struct simulation *sim,
+ const int max_iter,
+ const double rel_tol)
{
int iter, i;
- double r_norm_max = NAN, *tmp;
+ double r_norm_max = NAN;
+ const double omega = 1.2; /* SOR relaxation parameter */
for (iter = 0; iter < max_iter; ++iter) {
#ifdef DEBUG
@@ -672,33 +690,32 @@ implicit_1d_jacobian_poisson_solver(struct simulation *sim,
/* Neumann BCs resemble free surfaces */
/* set_bc_neumann(sim->g_ghost, sim->nz, +1, 0.0, sim->dz); */
+ /* SOR updates g_ghost in-place, using already-updated values */
for (i = 0; i < sim->nz; ++i)
- poisson_solver_1d_cell_update(i,
- sim->g_ghost,
- sim->g_local,
- sim->tmp_ghost,
- sim->g_r_norm,
- sim->dz,
- sim->xi);
+ poisson_solver_1d_sor_cell_update(i,
+ sim->g_ghost,
+ sim->g_local,
+ sim->g_r_norm,
+ sim->dz,
+ sim->xi,
+ omega);
r_norm_max = max(sim->g_r_norm, sim->nz);
- tmp = sim->tmp_ghost;
- sim->tmp_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);
#ifdef DEBUG
printf(".. Solution converged after %d iterations\n", iter + 1);
#endif
+ g_stats.poisson_iters += iter + 1;
return 0;
}
}
- fprintf(stderr, "implicit_1d_jacobian_poisson_solver: ");
+ fprintf(stderr, "implicit_1d_sor_poisson_solver: ");
fprintf(stderr, "Solution did not converge after %d iterations\n", iter + 1);
fprintf(stderr, ".. Residual normalized error: %f\n", r_norm_max);
+ g_stats.poisson_iters += iter + 1;
return 1;
}
@@ -825,7 +842,7 @@ coupled_shear_solver(struct simulation *sim,
compute_cooperativity_length(sim); /* Eq. 12 */
/* step 8, Eq. 11 */
- if (implicit_1d_jacobian_poisson_solver(sim,
+ if (implicit_1d_sor_poisson_solver(sim,
MAX_ITER_GRANULAR,
rel_tol))
exit(12);
@@ -896,6 +913,10 @@ coupled_shear_solver(struct simulation *sim,
temporal_increment(sim);
+ g_stats.coupled_iters += coupled_iter;
+ g_stats.stress_iters += stress_iter;
+ g_stats.timesteps++;
+
return 0;
}
@@ -949,3 +970,28 @@ set_coupled_fluid_transient_timestep(struct simulation *sim, const double safety
if (sim->dt > safety * dt)
sim->dt = safety * dt;
}
+
+void
+print_solver_stats(FILE *fp)
+{
+ fprintf(fp, "solver_stats: timesteps=%ld poisson_iters=%ld "
+ "darcy_iters=%ld coupled_iters=%ld stress_iters=%ld\n",
+ g_stats.timesteps, g_stats.poisson_iters,
+ g_stats.darcy_iters, g_stats.coupled_iters, g_stats.stress_iters);
+}
+
+void
+reset_solver_stats(void)
+{
+ g_stats.poisson_iters = 0;
+ g_stats.darcy_iters = 0;
+ g_stats.coupled_iters = 0;
+ g_stats.stress_iters = 0;
+ g_stats.timesteps = 0;
+}
+
+void
+add_darcy_iters(int iters)
+{
+ g_stats.darcy_iters += iters;
+}
diff --git a/simulation.h b/simulation.h
@@ -162,4 +162,7 @@ void set_coupled_fluid_transient_timestep(struct simulation *sim, const double s
double find_flux(const struct simulation *sim);
+void print_solver_stats(FILE *fp);
+void reset_solver_stats(void);
+
#endif
