commit 92b06adca3562c68372cf40983590ddea9c3f17e
parent b61f9f0ce52b135939c72bb5294f607638bd9dd1
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Wed, 14 Jan 2026 23:01:20 +0100
refactor(fluid): remove unused epsilon and untrack merge_request.md
Diffstat:
4 files changed, 10 insertions(+), 90 deletions(-)
diff --git a/fluid.c b/fluid.c
@@ -196,23 +196,12 @@ static double darcy_pressure_change_1d_impl(
}
}
-int darcy_solver_1d(struct simulation *sim, const int max_iter,
- const double rel_tol) {
- int i, iter, solved = 0;
- double epsilon, p_f_top, omega, r_norm_max = NAN, *k_n, *phi_n;
+int darcy_solver_1d(struct simulation *sim) {
+ int i, solved = 0;
+ double p_f_top, *k_n, *phi_n;
copy_values(sim->p_f_dot, sim->p_f_dot_old, sim->nz);
- /* choose integration method, parameter in [0.0; 1.0]
- * epsilon = 0.0: explicit
- * epsilon = 0.5: Crank-Nicolson
- * epsilon = 1.0: implicit */
- epsilon = 0.5;
-
- /* underrelaxation parameter in ]0.0; 1.0],
- * where omega = 1.0: no underrelaxation */
- omega = 1.0;
-
for (i = 0; i < sim->nz; ++i)
sim->p_f_dot_impl[i] = sim->p_f_dot_expl[i] = 0.0;
@@ -232,19 +221,8 @@ int darcy_solver_1d(struct simulation *sim, const int max_iter,
set_fluid_bcs(sim->p_f_ghost, sim, p_f_top);
set_fluid_bcs(sim->p_f_next_ghost, sim, p_f_top);
- /* explicit solution to pressure change */
- if (epsilon < 1.0) {
-#ifdef DEBUG
- printf("\nEXPLICIT SOLVER IN %s\n", __func__);
-#endif
- for (i = 0; i < sim->nz - 1; ++i)
- sim->p_f_dot_expl[i] = darcy_pressure_change_1d(
- i, sim->nz, sim->p_f_ghost, sim->phi, sim->phi_dot, sim->k, sim->dz,
- sim->beta_f, sim->alpha, sim->mu_f, sim->D);
- }
-
/* implicit solution with TDMA */
- if (epsilon > 0.0) {
+ {
#ifdef DEBUG
printf("\nIMPLICIT SOLVER IN %s\n", __func__);
#endif
@@ -327,22 +305,8 @@ int darcy_solver_1d(struct simulation *sim, const int max_iter,
/* Apply Boundary Conditions to Linear System */
/* Bottom (i=0): Neumann. p_{-1} = p_0 + C. */
- double bc_neumann_val;
- if (sim->D > 0.0) {
- /* For D > 0, we still use the ghost nodes set by set_fluid_bcs.
- p_ghost[0] = p_ghost[1] + ...
- The 'a[0]*p_{-1}' logic requires a[0] to be valid.
- For D>0, a[0] = -coeff.
- Equation: a[0]*p_{-1} + b[0]*p_0 + c[0]*p_1 = d[0].
- Substitute p_{-1} = p_0 + C.
- a[0]*p_0 + a[0]*C + b[0]*p_0 + c[0]*p_1 = d[0].
- (a[0]+b[0])*p_0 + c[0]*p_1 = d[0] - a[0]*C.
- Same algebra applies.
- */
- bc_neumann_val = sim->phi[0] * sim->rho_f * sim->G * sim->dz;
- } else {
- bc_neumann_val = sim->phi[0] * sim->rho_f * sim->G * sim->dz;
- }
+ /* Bottom (i=0): Neumann. p_{-1} = p_0 + C. */
+ double bc_neumann_val = sim->phi[0] * sim->rho_f * sim->G * sim->dz;
b[0] += a[0];
d[0] -= a[0] * bc_neumann_val;
@@ -376,8 +340,6 @@ int darcy_solver_1d(struct simulation *sim, const int max_iter,
add_darcy_iters(1);
solved = 1;
- } else {
- /* Explicit mode skipped */
}
for (i = 0; i < sim->nz - 1; ++i)
@@ -387,7 +349,6 @@ int darcy_solver_1d(struct simulation *sim, const int max_iter,
set_fluid_bcs(sim->p_f_ghost, sim, p_f_top);
set_fluid_bcs(sim->p_f_next_ghost, sim, p_f_top);
#ifdef DEBUG
- printf(".. epsilon = %g\n", epsilon);
puts(".. p_f_dot_expl:");
print_array(sim->p_f_dot_expl, sim->nz);
puts(".. p_f_dot_impl:");
diff --git a/fluid.h b/fluid.h
@@ -9,9 +9,7 @@ void hydrostatic_fluid_pressure_distribution(struct simulation *sim);
int set_largest_fluid_timestep(struct simulation *sim, const double safety);
-int darcy_solver_1d(struct simulation *sim,
- const int max_iter,
- const double rel_tol);
+int darcy_solver_1d(struct simulation *sim);
void add_darcy_iters(int iters);
diff --git a/merge_request.md b/merge_request.md
@@ -1,37 +0,0 @@
-# Merge Request: Numerical Performance Optimizations
-
-## Overview
-This MR implements significant numerical performance optimizations for the `cngf-pf` solver. The primary change involves replacing iterative solvers for the Poisson and Darcy equations with direct Tridiagonal Matrix Algorithm (TDMA) solvers.
-
-**Branch**: `feat/numerical-optimization`
-
-## Changes
-
-### 1. Solvers
-- **Poisson Equation (`simulation.c`)**: Replaced the iterative Successive Over-Relaxation (SOR) solver with a direct TDMA solver.
-- **Darcy Equation (`fluid.c`)**: Replaced the iterative implicit solver (Crank-Nicolson) with a Fully Implicit direct TDMA solver.
- - *Note*: This change improves stability and performance but modifies the time integration scheme from 2nd order (CN) to 1st order (Backward Euler), resulting in unconditional stability.
- - **Correction**: Updated Top BC handling to correctly respect hydrostatic corrections (ghost nodes).
- - **Feature**: Added support for constant diffusivity mode (`-D`) in the TDMA path.
-
-### 2. Math & Build
-- **Micro-optimizations**: Replaced `pow(x, 2.0)` and `pow(x, 3.0)` with direct multiplication.
-- **Loop Fusion**: Fused adjacent loops in `coupled_shear_solver` to improve cache locality.
-- **Compiler Flags**: Added `-O3 -march=native -flto` to `Makefile`.
-
-## Performance Results
-Benchmarks run on the `wet_large` test case (most intensive):
-
-| Metric | Baseline | Optimized | Improvement |
-| :--- | :--- | :--- | :--- |
-| **Elapsed Time** | 0.115s | 0.015s | **7.7x Speedup** |
-| **Poisson Iterations** | 24,435 | 3,600 | **6.8x Reduction** |
-| **Darcy Iterations** | 10,405 | 3,599 | **2.9x Reduction** |
-
-## Verification
-- **Correctness**: Validated using `test/accuracy.sh`.
-- **Regression Tests**: All standard tests passed (`make test`). Standard output files (`test/*.std`) were updated to reflect the improved precision and the solver change.
-- **Stability**: The transition to a Fully Implicit Darcy solver provides unconditional linear stability.
-
-## Artifacts
-- [Walkthrough](walkthrough.md): Detailed verification report and benchmark logs.
diff --git a/simulation.c b/simulation.c
@@ -609,9 +609,7 @@ void tridiagonal_solver(double *x, const double *a, const double *b,
free(d_prime);
}
-static int implicit_1d_sor_poisson_solver(struct simulation *sim,
- const int max_iter,
- const double rel_tol) {
+static int implicit_1d_sor_poisson_solver(struct simulation *sim) {
/*
* Replaced SOR solver with direct TDMA solver.
* System: -0.5*g[i-1] + (1 + C)*g[i] - 0.5*g[i+1] = C*g_local[i]
@@ -786,7 +784,7 @@ int coupled_shear_solver(struct simulation *sim, const int max_iter,
/* step 5, Eq. 13 */
if (sim->fluid && (sim->t > 0))
- if (darcy_solver_1d(sim, MAX_ITER_DARCY, 1e-3 * rel_tol))
+ if (darcy_solver_1d(sim))
exit(11);
/* step 6 */
@@ -797,7 +795,7 @@ int coupled_shear_solver(struct simulation *sim, const int max_iter,
compute_cooperativity_length(sim); /* Eq. 12 */
/* step 8, Eq. 11 */
- if (implicit_1d_sor_poisson_solver(sim, MAX_ITER_GRANULAR, rel_tol))
+ if (implicit_1d_sor_poisson_solver(sim))
exit(12);
/* step 9 */
