INTRO: Heat conduction + advection problem… 1D spherical coordinates.
For Dirichlet BC (bc_left_type = “Dirichlet” and bc_right_type = “Dirichlet”), the two codes (CN and MOL) return identical results up to machine precision. However, for Neuman BC (bc_left_type = “Neumann” and bc_right_type = “Neumann”), the temperature profiles differ rather significantly.
I’m pretty sure that there is no “obvious” bug in any code.
Perhaps someone can give me a hand here?
MOL Code:
# =========================================================================================
# Method-of-Lines Code
# =========================================================================================
#
# • Solves heat eq + grid advection: ∂T/∂t = α∇²T - v_grid ⋅ ∇T
# • Handles Dirichlet or Neumann BCs specified by constants below.
# • Uses Trapezoid() (Crank-Nicolson) with fixed dt.
# • Uses Double64 precision.
# • Standard MOL approach (solver integrates a, b, T).
#
# =========================================================================================
using DifferentialEquations, ComponentArrays, Printf
using OrdinaryDiffEq # For Trapezoid
using Plots, LaTeXStrings, Measures # Keep plotting essentials
using SciMLBase # For successful_retcode
using DoubleFloats # <<< Use Higher Precision
using Statistics # For mean
println("--- Minimal MOL Verification (Trapezoid, Double64) ---")
# ----------------------------------------------------------------------------
# Parameters (Double64)
# ----------------------------------------------------------------------------
const n_r = 101 # Spatial nodes (Int)
const tspan_end = Double64(0.05) # Simulation horizon
const alpha_target = Double64(1.5) # Target thermal diffusivity
const v_inner = Double64(3.0) # Prescribed inner boundary velocity
const v_outer = Double64(-3.0) # Prescribed outer boundary velocity
# --- Boundary Conditions --- <<< CHOOSE HERE >>>
#const bc_left_type = "Neumann" # "Dirichlet" or "Neumann"
#const bc_left_value = Double64(0.0) # Value (Temperature or dT/dr)
#const bc_right_type = "Neumann" # "Dirichlet" or "Neumann"
#const bc_right_value = Double64(-0.3) # Value (Temperature or dT/dr)
const bc_left_type = "Dirichlet"
const bc_left_value = Double64(100.0)
const bc_right_type = "Dirichlet"
const bc_right_value = Double64(20.0)
# --- Initial Conditions ---
const ra0 = Double64(1e-6) # Initial inner boundary position
const rf0 = Double64(1.0) # Initial outer radius
const T0_val = Double64(100.0) # Initial bulk temperature
# Time stepping to match CN
const cn_n_t = 30001
const cn_dt = tspan_end / Double64(cn_n_t - 1) # Will be Double64
println("BCs: $(bc_left_type) [$(bc_left_value)] / $(bc_right_type) [$(bc_right_value)]")
println("dt = ", cn_dt)
# ----------------------------------------------------------------------------
# Minimal Helper Functions (Return Double64)
# ----------------------------------------------------------------------------
radius_grid(a::Double64, b::Double64, nr::Int) = begin
h = (b - a) / Double64(nr - 1)
if !isfinite(h) || h <= 0.0 error("Invalid step h= $h in radius_grid") end
r = [a + Double64(i-1)*h for i in 1:nr]
dr = h
return (r, dr)
end
grid_velocity(r::Vector{Double64}, a::Double64, b::Double64, nr::Int, v_in::Double64, v_out::Double64)::Vector{Double64} = begin
ba_diff = b-a
if abs(ba_diff) < Double64(1e-15); return fill(Double64(0.0), nr);
else; return @. v_in + (r - a) / ba_diff * (v_out - v_in); end
end
# ----------------------------------------------------------------------------
# Simplified State Vector (Uses Double64)
# ----------------------------------------------------------------------------
make_state_simple(a::Double64, b::Double64, T::Vector{Double64}) = ComponentArray(a=a, b=b, T=T)
# ----------------------------------------------------------------------------
# rhs! Function (Spatial Operator L(T) including Advection)
# ----------------------------------------------------------------------------
function rhs_heat_final!(dX::ComponentArray{Double64}, X::ComponentArray{Double64}, _, t)
a::Double64 = X.a
b::Double64 = X.b
T::SubArray{Double64} = view(X.T, :)
# --- Grid Checks ---
if b <= a; dX .= 0.0; return nothing; end
local dr::Double64 = (b - a) / Double64(n_r - 1)
if dr <= Double64(1e-20); @warn "dr small ($dr)"; if dr <= 0.0; error("invalid dr."); end; end
r::Vector{Double64}, _ = radius_grid(a, b, n_r)
# --- Velocities ---
local dXa::Double64 = v_inner
local dXb::Double64 = v_outer
local v_grid::Vector{Double64} = grid_velocity(r, a, b, n_r, v_inner, v_outer)
# --- Temperature Equation: dT/dt = L(T) = α*∇²T - v_grid*∇T ---
dXT_calc = view(dX.T, :)
α::Double64 = alpha_target
# --- Precompute recurring terms ---
local dr2 = dr * dr
local inv_dr = Double64(1.0) / dr
local inv_2dr = inv_dr / Double64(2.0)
local two_alpha_dr2 = Double64(2.0) * α / dr2
@inbounds begin
# ----- Left Boundary (i=1) -----
if bc_left_type == "Dirichlet"
dXT_calc[1] = Double64(0.0) # Fixed value BC
elseif bc_left_type == "Neumann"
g_left = bc_left_value
v1 = v_grid[1]; inv_r1 = Double64(1.0) / r[1]
dXT_calc[1] = (T[2]-T[1])*two_alpha_dr2 + g_left * (α * (Double64(-2.0)*inv_dr + Double64(2.0)*inv_r1) - v1)
else; error("Unknown Left BC Type"); end
# ----- Interior Nodes (i=2 to n_r-1) -----
for i = 2:(n_r - 1)
local v_i::Double64 = v_grid[i]
local inv_ri::Double64 = Double64(1.0) / r[i]
local alpha_dr2::Double64 = α / dr2
local alpha_rdr::Double64 = α * inv_ri * inv_dr
local v_2dr::Double64 = v_i * inv_2dr
dXT_calc[i] = T[i-1]*(alpha_dr2 - alpha_rdr + v_2dr) +
T[i] *(-Double64(2.0) * alpha_dr2) +
T[i+1]*(alpha_dr2 + alpha_rdr - v_2dr)
end
# ----- Right Boundary (i=n_r) -----
if bc_right_type == "Dirichlet"
dXT_calc[n_r] = Double64(0.0) # Fixed value BC
elseif bc_right_type == "Neumann"
g_right = bc_right_value
vn = v_grid[n_r]; inv_rn = Double64(1.0) / r[n_r]
dXT_calc[n_r] = (T[n_r-1] - T[n_r]) * two_alpha_dr2 +
g_right * (α * (Double64(2.0)*inv_dr + Double64(2.0)*inv_rn) - vn)
else; error("Unknown Right BC Type"); end
end # end @inbounds
if any(!isfinite, dXT_calc); error("NaN/Inf dX.T"); end
dX.a = dXa; dX.b = dXb
return nothing
end
# ----------------------------------------------------------------------------
# Main run logic (outside function)
# ----------------------------------------------------------------------------
T₀_vec = fill(T0_val, n_r)
# Set initial state based on BC type
init_a = ra0
init_b = rf0
if bc_left_type == "Dirichlet"; T₀_vec[1] = bc_left_value; end
if bc_right_type == "Dirichlet"; T₀_vec[end] = bc_right_value; end
X₀ = make_state_simple(init_a, init_b, T₀_vec)
prob = ODEProblem(rhs_heat_final!, X₀, (Double64(0.0), tspan_end))
condition(u,t,i) = u.a >= u.b; affect!(i) = terminate!(i)
cb = DiscreteCallback(condition, affect!)
reltol=Double64(1e-12); abstol=Double64(1e-12) # Use tight tolerances
println("Running MOL simulation...")
sol = solve(prob, OrdinaryDiffEq.Trapezoid(autodiff=false);
reltol=reltol, abstol=abstol, adaptive=false, dt=cn_dt,
save_everystep=false, save_start=false, save_end=true,
callback=cb, maxiters = cn_n_t + 100)
println("MOL finished. Status: ", sol.retcode)
# ----------------------------------------------------------------------------
# Post-Processing
# ----------------------------------------------------------------------------
println("\n--- MOL Post-Processing ---")
final_state = isempty(sol.u) ? X₀ : sol.u[end]
final_t = isempty(sol.t) ? 0.0 : sol.t[end]
final_a = final_state.a; final_b = final_state.b
final_T_mol = final_state.T
final_dr = (final_b > final_a) ? (final_b - final_a) / Double64(n_r - 1) : NaN
@printf("MOL final t = %.6f\n", final_t)
@printf(" ra = %.17e rf = %.17e dr = %.3e\n", final_a, final_b, final_dr)
@printf(" T-min = %.17f T-max = %.17f T-mean = %.17f\n",
minimum(final_T_mol), maximum(final_T_mol), mean(final_T_mol))
# Plotting (Optional)
plot_valid = false
r_plot = Double64[]
if isfinite(final_a) && isfinite(final_b) && final_b > final_a + Double64(1e-15)
r_plot, _ = radius_grid(final_a, final_b, n_r)
plot_valid = true
end
plot_title = "MOL ($(bc_left_type)/$(bc_right_type), Double64)\n T(r) at t=$(round(Float64(final_t), sigdigits=3))s"
if plot_valid
plt = plot(Float64.(r_plot), Float64.(final_T_mol), xlabel="Radius (m)", ylabel="Temperature (K)", title=plot_title, size=(800, 600), legend=false)
display(plt)
else
println("Plotting skipped due to invalid final grid.")
end
println("\nFinal T vector (MOL):")
display(final_T_mol)
println("\nMOL Script finished.")
CN Code:
# =========================================================================================
# CN Code - Double64 Precision
# =========================================================================================
#
# • Solves conduction-advection problem with STANDARD Crank-Nicolson
# and 2nd ORDER ACCURATE Neumann Boundary Conditions (Dirichlet option handled).
# • Uses Double64 precision.
# • Corrected M matrix and RHS assembly for boundary conditions.
# • Uses Explicit Euler for grid update.
# • Wrapped main logic in a function to fix scope errors.
#
# =========================================================================================
# module CN_Final_Verified # Removed module
using LinearAlgebra
using Plots
using LaTeXStrings # For plot labels
using DoubleFloats # <<< Use Higher Precision
using Statistics # For mean
using Printf # For printing
# export run_cn # Not needed without module
println("--- Reference CN Verification (Double64) ---")
#---------------------------------------------------------------------
# Parameters (Double64)
#---------------------------------------------------------------------
const n_r = 101 # Spatial nodes (Int)
const total_time = Double64(0.05) # Simulation horizon
const n_t = 30001 # Number of time nodes
const α = Double64(1.5) # Target thermal diffusivity
const v_inner = Double64(3.0) # Prescribed inner boundary velocity
const v_outer = Double64(-3.0)# Prescribed outer boundary velocity
# --- Boundary Conditions --- <<< CHOOSE HERE >>>
const bc_left_type = "Neumann" # "Dirichlet" or "Neumann"
const bc_left_value = Double64(0.0) # Value (Temperature or dT/dr)
const bc_right_type = "Neumann" # "Dirichlet" or "Neumann"
const bc_right_value = Double64(-0.3) # Value (Temperature or dT/dr)
# const bc_left_type = "Dirichlet"
# const bc_left_value = Double64(100.0)
# const bc_right_type = "Dirichlet"
# const bc_right_value = Double64(20.0)
# --- Initial Conditions ---
const ra0 = Double64(1e-6) # Initial inner boundary position
const rf0 = Double64(1.0) # Initial outer radius
const T0_val = Double64(100.0) # Initial bulk temperature
println("BCs: $(bc_left_type) [$(bc_left_value)] / $(bc_right_type) [$(bc_right_value)]")
#---------------------------------------------------------------------
# Helper function to update boundary positions (Explicit Euler, Double64)
#---------------------------------------------------------------------
function update_boundaries!(ra::Double64, rf::Double64, dt::Double64, v_inner::Double64, v_outer::Double64)
new_ra = ra + v_inner * dt
new_rf = rf + v_outer * dt
return new_ra, new_rf
end
#---------------------------------------------------------------------
# Main simulation function.
#---------------------------------------------------------------------
function run_cn()
# Use constants defined outside the function
dt::Double64 = total_time / Double64(n_t - 1)
println("Reference CN using dt: ", dt)
# Initialize temperature (local to function)
T = fill(Double64(100.0), n_r)
# Set initial state based on BC type
if bc_left_type == "Dirichlet"; T[1] = bc_left_value; end
if bc_right_type == "Dirichlet"; T[end] = bc_right_value; end
t::Double64 = 0.0
# Set initial boundaries (local to function)
ra::Double64 = ra0
rf::Double64 = rf0
# Time-stepping loop using fixed number of steps
num_steps = n_t - 1
println("Reference CN running for $num_steps steps.")
final_r = Double64[] # Initialize local final_r
T_cn = similar(T) # Initialize local T_cn
for step in 1:num_steps
ra_prev::Double64 = ra; rf_prev::Double64 = rf
ra_new::Double64, rf_new::Double64 = update_boundaries!(ra, rf, dt, v_inner, v_outer)
if rf_new <= ra_new + Double64(1e-15); println("CN Warning: Domain collapse."); t=(step-1)*dt; final_r=[ra_prev + Double64(i-1)*(rf_prev-ra_prev)/Double64(n_r-1) for i in 1:n_r]; T_cn .= T; break; end
# Modify function-local variables (no `global` needed)
ra = ra_new; rf = rf_new
local dr::Double64 = (rf - ra) / Double64(n_r - 1)
if dr <= Double64(1e-20); @error "CN invalid dr."; t=(step-1)*dt; final_r=[ra_prev + Double64(i-1)*(rf_prev-ra_prev)/Double64(n_r-1) for i in 1:n_r]; T_cn .= T; break; end
r = [ra + Double64(i-1)*dr for i in 1:n_r]
A::Double64 = α * dt / (dr * dr)
B_dr::Double64 = α * dt / (Double64(2.0) * dr)
C_dr::Double64 = dt / (Double64(4.0) * dr)
inv_dr::Double64 = Double64(1.0) / dr
local velocity; local rf_ra_diff::Double64 = rf - ra
if abs(rf_ra_diff) < Double64(1e-15); velocity = r_val -> Double64(0.0)
else; velocity = r_val -> v_inner + (v_outer - v_inner) * (r_val - ra) / rf_ra_diff; end
M = zeros(Double64, n_r, n_r)
rhs = zeros(Double64, n_r)
#----- BOUNDARY CONDITIONS -----
dt_half = Double64(0.5) * dt
# --- LEFT BOUNDARY ---
if bc_left_type == "Dirichlet"
M[1, :] .= 0.0; M[1, 1] = Double64(1.0); rhs[1] = bc_left_value
elseif bc_left_type == "Neumann"
g_left = bc_left_value
v1 = velocity(r[1]); inv_r1 = Double64(1.0) / r[1]
L_T1_coeff_T1 = -2*α/(dr*dr); L_T1_coeff_T2 = 2*α/(dr*dr)
L_T1_const = g_left * (α * (Double64(-2.0)*inv_dr + Double64(2.0)*inv_r1) - v1)
M[1, 1] = Double64(1.0) - dt_half * L_T1_coeff_T1; M[1, 2] = -dt_half * L_T1_coeff_T2
rhs[1] = (Double64(1.0) + dt_half * L_T1_coeff_T1) * T[1] + (dt_half * L_T1_coeff_T2) * T[2] + dt_half * L_T1_const
else; error("Unknown left BC type"); end
# --- RIGHT BOUNDARY ---
if bc_right_type == "Dirichlet"
M[n_r, :] .= 0.0; M[n_r, n_r] = Double64(1.0); rhs[n_r] = bc_right_value
elseif bc_right_type == "Neumann"
g_right = bc_right_value
vn = velocity(r[end]); inv_rn = Double64(1.0) / r[end]
L_Tn_coeff_Tm1 = 2*α/(dr*dr); L_Tn_coeff_Tn = -2*α/(dr*dr)
L_Tn_const = g_right * (α*(Double64(2.0)*inv_dr + Double64(2.0)*inv_rn) - vn)
M[n_r, n_r-1] = -dt_half * L_Tn_coeff_Tm1; M[n_r, n_r] = Double64(1.0) - dt_half * L_Tn_coeff_Tn
rhs[n_r] = (dt_half * L_Tn_coeff_Tm1) * T[n_r-1] + (Double64(1.0) + dt_half * L_Tn_coeff_Tn) * T[n_r] + dt_half * L_Tn_const
else; error("Unknown right BC type"); end
#----- INTERIOR NODES -----
for j in 2:(n_r - 1)
local inv_rj::Double64 = Double64(1.0) / r[j]
local B_on_rj::Double64 = B_dr * inv_rj
local v_j::Double64 = velocity(r[j])
local C_adv_j::Double64 = C_dr * v_j
coeff_jm1 = α / (dr*dr) - α*inv_rj/dr + v_j/(Double64(2.0)*dr)
coeff_j = -Double64(2.0) * α / (dr*dr)
coeff_jp1 = α / (dr*dr) + α*inv_rj/dr - v_j/(Double64(2.0)*dr)
M[j, j-1] = - dt_half * coeff_jm1
M[j, j] = Double64(1.0) - dt_half * coeff_j
M[j, j+1] = - dt_half * coeff_jp1
rhs[j] = (dt_half * coeff_jm1) * T[j-1] + (Double64(1.0) + dt_half * coeff_j) * T[j] + (dt_half * coeff_jp1) * T[j+1]
end
# Solve the linear system
try T_new = M \ rhs; T = T_new # Update T
catch e
if isa(e, SingularException); @error "CN: Matrix M singular step $step, t=$t."; t=(step-1)*dt; final_r=[ra_prev + Double64(i-1)*(rf_prev-ra_prev)/Double64(n_r-1) for i in 1:n_r]; T_cn .= T; break;
else; rethrow(e); end
end
t = step * dt; final_r = r; T_cn .= T # Store last successful T and r
end # End time loop
println("Reference CN finished at t = ", t)
println("Reference CN final boundaries [ra, rf]: [", ra, ", ", rf, "]")
if isempty(final_r) && t >= total_time - dt/Double64(2.0); final_r=[ra + Double64(i-1)*(rf-ra)/Double64(n_r-1) for i in 1:n_r]; end
return final_r, T_cn, ra, rf # Return final calculated T
end
#---------------------------------------------------------------------
# Execution Logic
#---------------------------------------------------------------------
println("\n--- Running Reference CN (Corrected v3 Double64, $(bc_left_type)/$(bc_right_type) BCs) ---")
r_cn, T_cn, ra_cn, rf_cn = run_cn() # Call the function
println("Reference CN (Corrected v3 Double64) Final T range [T_min, T_max]: [", minimum(T_cn), ", ", maximum(T_cn), "] K")
# Plotting
plot_cn_title = "Standard CN v3 (Double64, $(bc_left_type)/$(bc_right_type) BCs)\n T(r) at t=0.05s"
if isempty(r_cn); println("\nReference CN plotting skipped.");
else; plot_cn = plot(Float64.(r_cn), Float64.(T_cn), xlabel="Radius (m)", ylabel="Temperature (K)", title=plot_cn_title, size=(800, 600), legend=false); display(plot_cn); println("\nReference CN Plot displayed."); end
println("\nFinal T vector (CN):")
display(T_cn)
println("\nReference CN Script finished.")