From 24546e99300e02f1d119093eff9193b48f97da3b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Sierro?= Date: Thu, 11 Nov 2021 12:14:15 +0100 Subject: [PATCH] solver not working --- src/scgenerator/solver.py | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/src/scgenerator/solver.py b/src/scgenerator/solver.py index 076a3db..9caef11 100644 --- a/src/scgenerator/solver.py +++ b/src/scgenerator/solver.py @@ -140,14 +140,15 @@ class ConservedQuantityIntegrator(Integrator): ) new_qty = self.conserved_quantity(new_state) - self.current_error = np.abs(new_qty - self.last_qty) / self.last_qty + self.current_error = np.abs(new_qty - self.last_qty) + error_ok = self.last_qty * self.tolerated_error - if self.current_error > 2 * self.tolerated_error: + if self.current_error > 2 * error_ok: h_next_step = h * 0.5 - elif self.tolerated_error < self.current_error <= 2.0 * self.tolerated_error: + elif error_ok < self.current_error <= 2.0 * error_ok: h_next_step = h / size_fac break - elif self.current_error < 0.1 * self.tolerated_error: + elif self.current_error < 0.1 * error_ok: h_next_step = h * size_fac break else: @@ -381,9 +382,9 @@ def RK4IP_step( expD = np.exp(h * 0.5 * init_linear) A_I = expD * spectrum - k1 = expD * init_nonlinear - k2 = nonlinear_operator(init_state.replace(A_I + k1 * 0.5 * h)) - k3 = nonlinear_operator(init_state.replace(A_I + k2 * 0.5 * h)) - k4 = nonlinear_operator(init_state.replace(expD * (A_I + h * k3))) + k1 = h * expD * init_nonlinear + k2 = h * nonlinear_operator(init_state.replace(A_I + k1 * 0.5)) + k3 = h * nonlinear_operator(init_state.replace(A_I + k2 * 0.5)) + k4 = h * nonlinear_operator(init_state.replace(expD * (A_I + k3))) - return expD * (A_I + h / 6 * (k1 + 2 * k2 + 2 * k3)) + h / 6 * k4 + return expD * (A_I + k1 / 6 + k2 / 3 + k3 / 3) + k4 / 6