commit d18eb7e0f5a9d7247f3a2641501f4ca795800d70
parent 3b38b840aa892559966075831c33ae15f8d9d74e
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Sat, 7 Feb 2026 09:17:02 +0100
refactor(simulation): reuse persistent solver workspace buffers
Diffstat:
3 files changed, 50 insertions(+), 38 deletions(-)
diff --git a/fluid.c b/fluid.c
@@ -93,7 +93,14 @@ static void set_fluid_bcs(double *p_f_ghost, struct simulation *sim,
int darcy_solver_1d(struct simulation *sim) {
int i, solved = 0;
- double p_f_top, *k_n, *phi_n;
+ double p_f_top;
+ double *a = sim->tdma_a;
+ double *b = sim->tdma_b;
+ double *c = sim->tdma_c;
+ double *d = sim->tdma_d;
+ double *x = sim->tdma_x;
+ double *k_n = sim->darcy_k_n;
+ double *phi_n = sim->darcy_phi_n;
for (i = 0; i < sim->nz; ++i)
sim->p_f_dot_impl[i] = sim->p_f_dot_expl[i] = 0.0;
@@ -119,15 +126,6 @@ int darcy_solver_1d(struct simulation *sim) {
#ifdef DEBUG
printf("\nIMPLICIT SOLVER IN %s\n", __func__);
#endif
- double *a = empty(sim->nz);
- double *b = empty(sim->nz);
- double *c = empty(sim->nz);
- double *d = empty(sim->nz);
- double *x = empty(sim->nz);
-
- k_n = zeros(sim->nz);
- phi_n = zeros(sim->nz);
-
/* Predictor step for nonlinear coefficients */
if (sim->transient)
for (i = 0; i < sim->nz; ++i) {
@@ -214,7 +212,8 @@ int darcy_solver_1d(struct simulation *sim) {
c[sim->nz - 1] = 0.0;
/* Solve */
- tridiagonal_solver(x, a, b, c, d, sim->nz);
+ tridiagonal_solver(x, a, b, c, d, sim->tdma_c_prime, sim->tdma_d_prime,
+ sim->nz);
/* Store result in p_f_dot (rate) */
for (i = 0; i < sim->nz; ++i) {
@@ -223,14 +222,6 @@ int darcy_solver_1d(struct simulation *sim) {
sim->p_f_dot_impl[i] = sim->p_f_dot[i];
}
- free(a);
- free(b);
- free(c);
- free(d);
- free(x);
- free(k_n);
- free(phi_n);
-
add_darcy_iters(1);
solved = 1;
}
diff --git a/simulation.c b/simulation.c
@@ -150,6 +150,15 @@ void prepare_arrays(struct simulation *sim) {
sim->I = zeros(sim->nz);
sim->tan_psi = zeros(sim->nz);
sim->old_val = empty(sim->nz);
+ sim->tdma_a = empty(sim->nz);
+ sim->tdma_b = empty(sim->nz);
+ sim->tdma_c = empty(sim->nz);
+ sim->tdma_d = empty(sim->nz);
+ sim->tdma_x = empty(sim->nz);
+ sim->tdma_c_prime = empty(sim->nz);
+ sim->tdma_d_prime = empty(sim->nz);
+ sim->darcy_k_n = empty(sim->nz);
+ sim->darcy_phi_n = empty(sim->nz);
}
void free_arrays(struct simulation *sim) {
@@ -177,6 +186,15 @@ void free_arrays(struct simulation *sim) {
free(sim->I);
free(sim->tan_psi);
free(sim->old_val);
+ free(sim->tdma_a);
+ free(sim->tdma_b);
+ free(sim->tdma_c);
+ free(sim->tdma_d);
+ free(sim->tdma_x);
+ free(sim->tdma_c_prime);
+ free(sim->tdma_d_prime);
+ free(sim->darcy_k_n);
+ free(sim->darcy_phi_n);
}
static void warn_parameter_value(const char message[], const double value,
@@ -562,7 +580,8 @@ double residual(double new_val, double old_val) {
}
void tridiagonal_solver(double *x, const double *a, const double *b,
- const double *c, double *d, int n) {
+ const double *c, const double *d, double *c_prime,
+ double *d_prime, int n) {
/*
* TDMA (Thomas Algorithm) solver for tridiagonal system
* a: sub-diagonal (means a[i] is coeff for x[i-1])
@@ -578,8 +597,6 @@ void tridiagonal_solver(double *x, const double *a, const double *b,
* caller or zeroed)
*/
int i;
- double *c_prime = empty(n);
- double *d_prime = empty(n);
/* Forward sweep */
c_prime[0] = c[0] / b[0];
@@ -596,9 +613,6 @@ void tridiagonal_solver(double *x, const double *a, const double *b,
for (i = n - 2; i >= 0; i--) {
x[i] = d_prime[i] - c_prime[i] * x[i + 1];
}
-
- free(c_prime);
- free(d_prime);
}
static int implicit_1d_sor_poisson_solver(struct simulation *sim) {
@@ -609,11 +623,11 @@ static int implicit_1d_sor_poisson_solver(struct simulation *sim) {
*/
int i;
int n = sim->nz;
- double *a = empty(n);
- double *b = empty(n);
- double *c = empty(n);
- double *d = empty(n);
- double *x = empty(n);
+ double *a = sim->tdma_a;
+ double *b = sim->tdma_b;
+ double *c = sim->tdma_c;
+ double *d = sim->tdma_d;
+ double *x = sim->tdma_x;
/* Set up TDMA arrays */
for (i = 0; i < n; i++) {
@@ -637,7 +651,7 @@ static int implicit_1d_sor_poisson_solver(struct simulation *sim) {
a[0] = 0.0;
c[n - 1] = 0.0;
- tridiagonal_solver(x, a, b, c, d, n);
+ tridiagonal_solver(x, a, b, c, d, sim->tdma_c_prime, sim->tdma_d_prime, n);
/* Copy result back to ghost array (indices 1 to n) */
set_bc_dirichlet(sim->g_ghost, sim->nz, -1, 0.0);
@@ -646,12 +660,6 @@ static int implicit_1d_sor_poisson_solver(struct simulation *sim) {
sim->g_ghost[i + 1] = x[i];
}
- free(a);
- free(b);
- free(c);
- free(d);
- free(x);
-
g_stats.poisson_iters += 1; /* Count as 1 iteration for stats */
return 0;
}
diff --git a/simulation.h b/simulation.h
@@ -125,6 +125,18 @@ struct simulation {
double *I; /* inertia number [-] */
double *tan_psi; /* tan(dilatancy_angle) [-] */
double *old_val; /* temporary storage for iterative solvers */
+
+ /* Persistent solver workspace (size nz unless noted otherwise), allocated
+ * once in prepare_arrays() and reused by Darcy/Poisson TDMA paths. */
+ double *tdma_a; /* shared TDMA sub-diagonal coefficients */
+ double *tdma_b; /* shared TDMA diagonal coefficients */
+ double *tdma_c; /* shared TDMA super-diagonal coefficients */
+ double *tdma_d; /* shared TDMA RHS coefficients */
+ double *tdma_x; /* shared TDMA solution vector */
+ double *tdma_c_prime; /* shared TDMA forward-sweep scratch */
+ double *tdma_d_prime; /* shared TDMA forward-sweep scratch */
+ double *darcy_k_n; /* Darcy predictor permeability workspace */
+ double *darcy_phi_n; /* Darcy predictor porosity workspace */
};
void init_sim(struct simulation *sim);
@@ -158,6 +170,7 @@ void print_solver_stats(FILE *fp);
void reset_solver_stats(void);
void tridiagonal_solver(double *x, const double *a, const double *b,
- const double *c, double *d, int n);
+ const double *c, const double *d, double *c_prime,
+ double *d_prime, int n);
#endif
