cngf-pf

commit f3eb6b6c4944f8836c67355a082a8dc92d2bca04
parent 76e6b4ac97282ee09d512c61cda430007884b345
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date:   Wed, 14 Jan 2026 21:23:30 +0100

perf(fluid): precompute reciprocal terms in Darcy solver to avoid repeated divisions

Diffstat:
M
fluid.c
 | 
16
++++++++++------

1 file changed, 10 insertions(+), 6 deletions(-)
diff --git a/fluid.c b/fluid.c
@@ -181,6 +181,7 @@ darcy_pressure_change_1d_impl(const int i,
                               const double D)
 {
 	double k_, div_k_grad_p, k_zn, k_zp, rhs_term, omega = 1.0;
+	double inv_fluid_term, inv_porosity_term;
 
 	if (D > 0.0)
 		return D * (p_f_ghost_in[i + 2]
@@ -197,19 +198,22 @@ darcy_pressure_change_1d_impl(const int i,
 		else
 			k_zp = k[i + 1];
 
-		rhs_term = dt / ((alpha + beta_f * phi[i]) * mu_f)
+		/* Precompute reciprocal terms to avoid repeated divisions */
+		inv_fluid_term = 1.0 / ((alpha + beta_f * phi[i]) * mu_f);
+		inv_porosity_term = 1.0 / ((alpha + beta_f * phi[i]) * (1.0 - phi[i]));
+
+		rhs_term = dt * inv_fluid_term
 		           * ((2.0 * k_zp * k_ / (k_zp + k_) / (dz * dz))
 		              + (2.0 * k_zn * k_ / (k_zn + k_) / (dz * dz)));
 
 		p_f_ghost_out[i + 1] = 1.0 / (1.0 + rhs_term)
 		                       * (p_f_old_val[i + 1] + dt
-		                       * (1.0 / ((alpha + beta_f * phi[i]) * mu_f)
+		                       * (inv_fluid_term
 		                       * (2.0 * k_zp * k_ / (k_zp + k_)
 		                       * (p_f_ghost_in[i + 2]) / dz
 		                       - 2.0 * k_zn * k_ / (k_zn + k_)
 		                       * -p_f_ghost_in[i] / dz) / dz
-		                       - 1.0 / ((alpha + beta_f * phi[i])
-		                       * (1.0 - phi[i])) * phi_dot[i]));
+		                       - inv_porosity_term * phi_dot[i]));
 		p_f_ghost_out[i + 1] = omega * p_f_ghost_out[i + 1]
 		                       + (1.0 - omega) * p_f_ghost_in[i + 1];
 
@@ -228,8 +232,8 @@ darcy_pressure_change_1d_impl(const int i,
 		       k_zn, k_, k_zp);
 #endif
 		/* use the values from the next time step as the time derivative for this iteration */
-		return 1.0 / ((alpha + beta_f * phi[i]) * mu_f) * div_k_grad_p
-		       - 1.0 / ((alpha + beta_f * phi[i]) * (1.0 - phi[i])) * phi_dot[i];
+		return inv_fluid_term * div_k_grad_p
+		       - inv_porosity_term * phi_dot[i];
 	}
 }