Give Δt to step forward TPM

Astroshaper · Oct 14, 2023 · d52b82c · d52b82c
1 parent ac3e319
commit d52b82c
Show file tree

Hide file tree

Showing 3 changed files with 34 additions and 21 deletions.
diff --git a/src/TPM.jl b/src/TPM.jl
@@ -359,7 +359,8 @@ function run_TPM!(stpm::SingleTPM, ephem, savepath)
         ProgressMeter.next!(p; showvalues)
 
         nₜ == length(ephem.time) && break  # Stop to update the temperature at the final step
-        update_temperature!(stpm)
+        Δt = ephem.time[nₜ+1] - ephem.time[nₜ]
+        update_temperature!(stpm, Δt)
     end
 
     jldsave(savepath; stpm, ephem, surf_temps, forces, torques)
@@ -427,7 +428,8 @@ function run_TPM!(btpm::BinaryTPM, ephem, savepath)
 
         ## Update temperature distribution
         nₜ == length(ephem.time) && break  # Stop to update the temperature at the final step
-        update_temperature!(btpm)
+        Δt = ephem.time[nₜ+1] - ephem.time[nₜ]
+        update_temperature!(btpm, Δt)
     end
 
     jldsave(savepath; btpm, ephem, surf_temps, forces, torques)

diff --git a/src/heat_conduction.jl b/src/heat_conduction.jl
@@ -5,38 +5,38 @@
 # ****************************************************************
 
 """
-    update_temperature!(stpm::SingleTPM)
+    update_temperature!(stpm::SingleTPM, Δt)
 
 Calculate the temperature for the next time step based on 1D heat conduction equation.
 
 # Arguments
 - `stpm` : Thermophysical model for a single asteroid
 """
-function update_temperature!(stpm::SingleTPM)
+function update_temperature!(stpm::SingleTPM, Δt)
     if stpm.SOLVER isa ForwardEulerSolver
-        forward_euler!(stpm)
+        forward_euler!(stpm, Δt)
     elseif stpm.SOLVER isa BackwardEulerSolver
-        backward_euler!(stpm)
+        backward_euler!(stpm, Δt)
     elseif stpm.SOLVER isa CrankNicolsonSolver
-        crank_nicolson!(stpm)
+        crank_nicolson!(stpm, Δt)
     else
         error("The solver is not implemented.")
     end
 end
 
 
 """
-    update_temperature!(btpm::BinaryTPM)
+    update_temperature!(btpm::BinaryTPM, Δt)
 
 Calculate the temperature for the next time step based on 1D heat conductivity equation.
 
 # Arguments
 - `btpm` : Thermophysical model for a binary asteroid
 - `nₜ`   : Index of the current time step
 """
-function update_temperature!(btpm::BinaryTPM)
-    update_temperature!(btpm.pri)
-    update_temperature!(btpm.sec)
+function update_temperature!(btpm::BinaryTPM, Δt)
+    update_temperature!(btpm.pri, Δt)
+    update_temperature!(btpm.sec, Δt)
 end
 
 # ****************************************************************
@@ -45,20 +45,30 @@ end
 
 
 """
-    forward_euler!(stpm::SingleTPM, nₜ::Integer)
+    forward_euler!(stpm::SingleTPM, Δt)
 
 Predict the temperature at the next time step by the forward Euler method.
 - Explicit in time
 - First order in time
 In this function, the heat conduction equation is non-dimensionalized in time and length.
+
+# Arguments
+- `stpm` : Thermophysical model for a single asteroid
+- `Δt`   : Time step [sec]
 """
-function forward_euler!(stpm::SingleTPM)
+function forward_euler!(stpm::SingleTPM, Δt)
     T = stpm.temperature
     Nz = size(T, 1)
     Ns = size(T, 2)
 
     for nₛ in 1:Ns
-        λ = (stpm.thermo_params.λ isa Real ? stpm.thermo_params.λ : stpm.thermo_params.λ[nₛ])
+        P  = stpm.thermo_params.P
+        Δz = stpm.thermo_params.Δz
+        l  = (stpm.thermo_params.l isa Real ? stpm.thermo_params.l : stpm.thermo_params.l[nₛ])
+
+        λ = (Δt/P) / (Δz/l)^2 / 4π
+        λ ≥ 0.5 && error("The forward Euler method is unstable because λ = $λ. This should be less than 0.5.")
+
         for nz in 2:(Nz-1)
             stpm.SOLVER.T[nz] = (1-2λ)*T[nz, nₛ] + λ*(T[nz+1, nₛ] + T[nz-1, nₛ])  # Predict temperature at next time step
         end
@@ -73,15 +83,15 @@ end
 
 
 """
-    backward_euler!(stpm::SingleTPM)
+    backward_euler!(stpm::SingleTPM, Δt)
 
 Predict the temperature at the next time step by the backward Euler method.
 - Implicit in time (Unconditionally stable in the heat conduction equation)
 - First order in time
 - Second order in space
 In this function, the heat conduction equation is non-dimensionalized in time and length.
 """
-function backward_euler!(stpm::SingleTPM)
+function backward_euler!(stpm::SingleTPM, Δt)
     # T = stpm.temperature
     # Nz = size(T, 1)
     # Ns = size(T, 2)
@@ -114,15 +124,15 @@ end
 
 
 """
-    crank_nicolson!(stpm::SingleTPM)
+    crank_nicolson!(stpm::SingleTPM, Δt)
 
 Predict the temperature at the next time step by the Crank-Nicolson method.
 - Implicit in time (Unconditionally stable in the heat conduction equation)
 - Second order in time
 - Second order in space
 In this function, the heat conduction equation is non-dimensionalized in time and length.
 """
-function crank_nicolson!(stpm::SingleTPM)
+function crank_nicolson!(stpm::SingleTPM, Δt)
     # T = stpm.temperature
     # Nz = size(T, 1)
     # Ns = size(T, 2)

diff --git a/test/heat_conduction_1D.jl b/test/heat_conduction_1D.jl
@@ -65,10 +65,11 @@
     ##= Run TPM =##
     for nₜ in eachindex(ephem.time)
         nₜ == length(et_range) && break  # Stop to update the temperature at the final step
+        Δt = ephem.time[nₜ+1] - ephem.time[nₜ]
 
-        AsteroidThermoPhysicalModels.forward_euler!(stpm_FE)
-        AsteroidThermoPhysicalModels.backward_euler!(stpm_BE)
-        AsteroidThermoPhysicalModels.crank_nicolson!(stpm_CN)
+        AsteroidThermoPhysicalModels.forward_euler!(stpm_FE, Δt)
+        AsteroidThermoPhysicalModels.backward_euler!(stpm_BE, Δt)
+        AsteroidThermoPhysicalModels.crank_nicolson!(stpm_CN, Δt)
     end
 
     ##= Save data =##