From b7e358cee08f3630b1851fc6ceccd4abe88187a9 Mon Sep 17 00:00:00 2001
From: Etienne <etienne.gavazzi@uit.no>
Date: Fri, 15 Dec 2023 16:30:35 +0000
Subject: [PATCH] fix #22

---
 src/main.jl        |  2 +-
 src/utilitaries.jl | 50 +++++++++++++++++++++++++++++++++-------------
 2 files changed, 37 insertions(+), 15 deletions(-)

diff --git a/src/main.jl b/src/main.jl
index 8ae3b0e..3be4dbd 100644
--- a/src/main.jl
+++ b/src/main.jl
@@ -79,7 +79,7 @@ function calculate_e_transport(altitude_max, θ_lims, E_max, B_angle_to_zenith,
 
     ## And save the simulation parameters in it
     save_parameters(altitude_max, θ_lims, E_max, B_angle_to_zenith, t_sampling, t, n_loop,
-        Nthreads, CFL_number, INPUT_OPTIONS, savedir)
+        CFL_number, INPUT_OPTIONS, savedir)
     save_neutrals(h_atm, n_neutrals, ne, Te, savedir)
 
     # Initialize arrays for the ionization collisions part of the energy degradation
diff --git a/src/utilitaries.jl b/src/utilitaries.jl
index 15765c6..3b492a4 100644
--- a/src/utilitaries.jl
+++ b/src/utilitaries.jl
@@ -20,6 +20,32 @@ function v_of_E(E)
 	return v
 end
 
+# function CFL_criteria(t, h_atm, v, CFL_number=64)
+#     dt = t[2] - t[1]
+#     dz = h_atm[2] - h_atm[1]
+
+#     # The Courant-Freidrichs-Lewy (CFL) number normally hase to be small (<4) to ensure numerical
+#     # stability. However, as a Crank-Nicolson scheme is always stable, we can take a bigger CFL. We
+#     # should be careful about numerical accuracy though.
+#     # For Gaussian inputs (or similar), it seems that the CFL can be set to 64 without major effects
+#     # on the results, while reducing computational time tremendously
+#     CFL = v * dt / dz
+#     n_factors = 2 .^ collect(0:22)
+#     iFactor = 1
+#     # This while loop effectively reduces dt by a factor of 2 at each iteration and check if the new
+#     # CFL is < 64. If not, it continues reducing dt.
+#     t_finer = t
+#     while (CFL > CFL_number) && (iFactor < length(n_factors))
+#         t_finer = range(t[1], t[end], length(t) * n_factors[iFactor] + 1 - n_factors[iFactor])
+#         dt = t_finer[2] - t_finer[1]
+#         CFL = v * dt / dz
+#         iFactor += 1
+#     end
+#     CFL_factor = n_factors[max(1, iFactor - 1)]
+
+#     return t_finer, CFL_factor
+# end
+
 function CFL_criteria(t, h_atm, v, CFL_number=64)
     dt = t[2] - t[1]
     dz = h_atm[2] - h_atm[1]
@@ -29,19 +55,15 @@ function CFL_criteria(t, h_atm, v, CFL_number=64)
     # should be careful about numerical accuracy though.
     # For Gaussian inputs (or similar), it seems that the CFL can be set to 64 without major effects
     # on the results, while reducing computational time tremendously
-    CFL = v * dt / dz
-    n_factors = 2 .^ collect(0:22)
-    iFactor = 1
-    # This while loop effectively reduces dt by a factor of 2 at each iteration and check if the new
-    # CFL is < 64. If not, it continues reducing dt.
-    t_finer = t
-    while (CFL > CFL_number) && (iFactor < length(n_factors))
-        t_finer = range(t[1], t[end], length(t) * n_factors[iFactor] + 1 - n_factors[iFactor])
-        dt = t_finer[2] - t_finer[1]
-        CFL = v * dt / dz
-        iFactor += 1
-    end
-    CFL_factor = n_factors[max(1, iFactor - 1)]
+
+    # Calculate the maximum dt that still satisfies the CFL criteria.
+    dt_max = CFL_number * dz / v
+
+    # Then find the smallest integer dt can be divided by and still be smaller than dt_max.
+    CFL_factor = ceil(Int, dt / dt_max) # that integer is the CFL_factor
+
+    # Now make a t_finer enumbe
+    t_finer = range(t[1], t[end], CFL_factor * (length(t) - 1) + 1)
 
     return t_finer, CFL_factor
 end
@@ -116,7 +138,7 @@ function save_parameters(altitude_max, θ_lims, E_max, B_angle_to_zenith, t_samp
         write(f, "t = $t \n")
         write(f, "n_loop = $n_loop \n")
         write(f, "\n")
-        write(f, "Nthreads = $Nthreads \n")
+        # write(f, "Nthreads = $Nthreads \n")
         write(f, "CFL_number = $CFL_number")
         write(f, "\n")
         write(f, "input_options = $INPUT_OPTIONS \n")