update energy equation with correct area expansion term

src/simulation/electronenergy.jl

@@ -1,7 +1,7 @@
 
 
 function update_electron_energy!(U, params, dt)
-    (;Δz_cell, Δz_edge, index, config, cache, Te_L, Te_R) = params
+    (;z_cell, Δz_cell, Δz_edge, index, config, cache, Te_L, Te_R) = params
     (;Aϵ, bϵ, nϵ, ue, ne, Tev, channel_area, dA_dz, κ) = cache
     implicit = params.config.implicit_energy
     explicit = 1 - implicit
@@ -97,30 +97,34 @@ function update_electron_energy!(U, params, dt)
             FR_factor_R = 5/3 * ueR - κR / ΔzR / neR
         end
 
-        # Fluxes at left and right boundary
-        FL = FL_factor_L * nϵL + FL_factor_C * nϵ0 + FL_factor_R * nϵR
-        FR = FR_factor_L * nϵL + FR_factor_C * nϵ0 + FR_factor_R * nϵR
-
         # Contribution to implicit part from fluxes
         Aϵ.d[i] = (FR_factor_C - FL_factor_C) / Δz
         Aϵ.dl[i-1] = (FR_factor_L - FL_factor_L) / Δz
         Aϵ.du[i] = (FR_factor_R - FL_factor_R) / Δz
 
+        # Contribution to implicit part from area expansion term
+        dlnA_dz = dA_dz[i] / channel_area[i]
+        dL, d0, dR = central_diff_coeffs(z_cell[i-1], z_cell[i], z_cell[i+1])
+        Aϵ.d[i] -= (5/3 * ue0 + d0 * κ0 / ne0) * dlnA_dz
+        Aϵ.dl[i-1] -= dL * κ0 / neL * dlnA_dz
+        Aϵ.du[i] -= dR * κ0 / neR * dlnA_dz
+
         # Contribution to implicit part from timestepping
         Aϵ.d[i]    = 1.0 + implicit * dt * Aϵ.d[i]
         Aϵ.dl[i-1] = implicit * dt * Aϵ.dl[i-1]
         Aϵ.du[i]   = implicit * dt * Aϵ.du[i]
 
         # Explicit flux
+        # dF / dz
+        FL = FL_factor_L * nϵL + FL_factor_C * nϵ0 + FL_factor_R * nϵR
+        FR = FR_factor_L * nϵL + FR_factor_C * nϵ0 + FR_factor_R * nϵR
         F_explicit = (FR - FL) / Δz
 
-        # Term to allow for changing area
-        dlnA_dz = dA_dz[i] / channel_area[i]
-        flux = 5/3 * nϵ0 * ue0
+        # (5/3 nϵ ue + κ dT/dz) d(lnA)/dz
+        Q_area = (5/3 * nϵ0 * ue0 + κ0 * (dL * Tev[i-1] + d0 * Tev[i] + dR * Tev[i+1])) * dlnA_dz
 
         # Explicit right-hand-side
-        bϵ[i] = nϵ[i] + dt * (Q - explicit * F_explicit)
-        #bϵ[i] -= dt * flux * dlnA_dz
+        bϵ[i] = nϵ[i] + dt * (Q - explicit * (F_explicit + Q_area))
     end
 
     # Solve equation system using Thomas algorithm

src/simulation/solution.jl

@@ -67,28 +67,37 @@ function solve(U, params, tspan; saveat)
 
         t += params.dt[]
 
-        # Update heavy species
-        ssprk22_step!(U, update_heavy_species!, params, t, params.dt[])
+        try
 
-        # Update electron quantities
-        update_electrons!(U, params, t)
+            # Update heavy species
+            ssprk22_step!(U, update_heavy_species!, params, t, params.dt[])
 
-        # Update plume geometry
-        update_plume_geometry!(U, params)
+            # Update electron quantities
+            update_electrons!(U, params, t)
 
-        # Check for NaNs and terminate if necessary
-        nandetected = false
-        infdetected = false
+            # Update plume geometry
+            update_plume_geometry!(U, params)
 
+        catch e
+            if (e isa InterruptException)
+                @warn("Simulation interrupted at time $t.")
+                retcode = :interrupt
+                break
+            else
+                @warn ("Error detected at time $t.")
+                retcode = :errordetected
+                break
+            end
+        end
+
+        # Check for NaNs and terminate if necessary
         @inbounds for j in 1:ncells, i in 1:nvars
             if isnan(U[i, j])
                 @warn("NaN detected in variable $i in cell $j at time $(t)")
-                nandetected = true
                 retcode = :NaNDetected
                 break
             elseif isinf(U[i, j])
                 @warn("Inf detected in variable $i in cell $j at time $(t)")
-                infdetected = true
                 retcode = :InfDetected
                 break
             end
@@ -121,10 +130,6 @@ function solve(U, params, tspan; saveat)
 
             save_ind += 1
         end
-
-        if nandetected || infdetected
-            break
-        end
     end
 
     ind = min(save_ind, length(savevals)+1)-1

src/simulation/sourceterms.jl

@@ -138,41 +138,17 @@ function electron_kinetic_energy(U, params, i)
 end
 
 function source_electron_energy(U, params, i)
-    (;cache) = params
-    ne = cache.ne[i]
-    ue = cache.ue[i]
-    ∇ϕ = cache.∇ϕ[i]
-    B = cache.B[i]
-    νe = cache.νe[i]
+    (;ne, ue, ∇ϕ) = params.cache
 
     # Compute ohmic heating term, which is the rate at which energy is transferred out of the electron
-    # drift (kinetic energy) into thermal inergy
-    if params.config.LANDMARK
-        # Neglect kinetic energy, so the rate of energy transfer into thermal energy is equivalent to
-        # the total input power into the electrons (j⃗ ⋅ E⃗ = -mₑnₑ|uₑ|²νₑ)
-        ohmic_heating  = ne * ue * ∇ϕ
-    else
-        # Do not neglect kinetic energy, so ohmic heating term is mₑnₑ|uₑ|²νₑ + ue ∇pe = 2nₑKνₑ + ue ∇pe
-        # where K is the electron bulk kinetic energy, 1/2 * mₑ|uₑ|²
-        K = params.cache.K[i]
-        νe = params.cache.νe[i]
-        ∇pe = params.cache.∇pe[i]
-        ∇Te = params.cache.∇Te[i]
-        ji = params.cache.ji[i]
-        Te = params.cache.Tev[i]
-        j_perp = ji - e * ne * ue
-        hall_param = e * B / (me * νe)
-        j_theta = hall_param * (j_perp - ji) + 1.5 * ne / B * ∇Te
-        σ = e^2 * ne / (me * νe)
-        ion_heating_rate = 3 * me / params.mi * νe * ne * Te
-        ohmic_heating = ne * ue * ∇ϕ - ion_heating_rate
-    end
+    # drift (kinetic energy) into thermal energy
+    ohmic_heating = ne[i] * ue[i] * ∇ϕ[i]
 
     # add excitation losses to total inelastic losses
     excitation_losses!(U, params, i)
 
     # compute wall losses
-    wall_loss = ne * wall_power_loss(params.config.wall_loss_model, U, params, i)
+    wall_loss = ne[i] * wall_power_loss(params.config.wall_loss_model, U, params, i)
 
     params.cache.wall_losses[i] = wall_loss
     params.cache.ohmic_heating[i] = ohmic_heating