Materials starter - Heat equation convergence

MODEL "Materials starter - Heat equation convergence"
DESCRIPTION "Separate spatial discretization error from time-integration error for a homogeneous slab's sine temperature mode."
VERSION "1.0.0"
TAGS "materials", "starter", "heat", "PDE", "convergence"

strict

input FourierNumber = 0.1
  VALIDATE BETWEEN 0.001 AND 0.5 MESSAGE "Use a positive dimensionless time alpha*t/L^2 from 0.001 to 0.5 in this example."
input SlabLength_mm = 10
  VALIDATE BETWEEN 0.001 AND 10000 MESSAGE "Use a positive full slab length from 0.001 to 10000 millimeters."
input Diffusivity_mm2_s = 1
  VALIDATE BETWEEN 0.000001 AND 100 MESSAGE "Use a positive constant diffusivity from 0.000001 to 100 square millimeters per second."

Length_m = CONVERT(SlabLength_mm, "mm", "m")
Diffusivity_m2_s = CONVERT(Diffusivity_mm2_s, "mm^2/s", "m^2/s")
# Solve in xi=x/L and Fo=alpha*t/L^2, with theta=(T-T_boundary)/amplitude.
# The normalized domain and field are dimensionless, avoiding hidden units
# in the PDE lambda. Both faces have theta=0; initially theta=sin(pi*xi).
# At Fo=0.1: continuum center 0.3727078389; 21/41-node centers about
# 0.3734643407/0.3728969372; spatial error ratio about 4.0006.
CoarseField = PDE_SOLVE_1D(LAMBDA(time, position, value, gradient, curvature, curvature), LAMBDA(position, SIN(PI() * position)), 0, 0, 0, 1, 21, FourierNumber, {samples: 2, rtol: 1e-9, atol: 1e-12})
FineField = PDE_SOLVE_1D(LAMBDA(time, position, value, gradient, curvature, curvature), LAMBDA(position, SIN(PI() * position)), 0, 0, 0, 1, 41, FourierNumber, {samples: 2, rtol: 1e-9, atol: 1e-12})
CoarseCenter = INDEX(CoarseField, 2, 11)
FineCenter = INDEX(FineField, 2, 21)
ContinuumCenter = EXP(-((PI()^2) * FourierNumber))
CoarseSpacing = 1 / 20
FineSpacing = 1 / 40
CoarseSemidiscrete = EXP(-((4 / CoarseSpacing^2) * SIN(PI() * CoarseSpacing / 2)^2 * FourierNumber))
FineSemidiscrete = EXP(-((4 / FineSpacing^2) * SIN(PI() * FineSpacing / 2)^2 * FourierNumber))
CoarseSpatialError = ABS(CoarseCenter - ContinuumCenter)
FineSpatialError = ABS(FineCenter - ContinuumCenter)

output ElapsedTime_s = FourierNumber * Length_m^2 / Diffusivity_m2_s
output Center_nx21 = CoarseCenter
output Center_nx41 = FineCenter
output ExactContinuumCenter = ContinuumCenter
output SpatialError_nx21 = CoarseSpatialError
output SpatialError_nx41 = FineSpatialError
output SpatialErrorRatio = CoarseSpatialError / FineSpatialError
output ObservedSpatialOrder = LN(SpatialErrorRatio) / LN(2)
output TemporalError_nx21 = ABS(CoarseCenter - CoarseSemidiscrete)
output TemporalError_nx41 = ABS(FineCenter - FineSemidiscrete)
output Assumptions = "Homogeneous 1D slab, constant positive thermal diffusivity, equal fixed face temperatures, no source term, and a sine initial excess-temperature profile. nx counts endpoints. The exact semidiscrete mode isolates time-integration error; solver rtol/atol are local controls, not global error bounds."
output ScenarioGuide = "Save Baseline at FourierNumber=0.1. Set it to 0.2 and save Later cooling, then restore Baseline or reset. At fixed Fourier number, changing length or diffusivity changes elapsed physical time, not the normalized field."

Profile = MAP(SEQUENCE(41, 1, 0, 1), (sample_index) => {
  position: sample_index / 40,
  numerical: INDEX(FineField, 2, sample_index + 1),
  continuum: ContinuumCenter * SIN(PI() * sample_index / 40)
})
A1 = "Elapsed physical time (s)"
B1 = ElapsedTime_s
A2 = "Center, nx=21"
B2 = Center_nx21
A3 = "Center, nx=41"
B3 = Center_nx41
A4 = "Exact continuum center"
B4 = ExactContinuumCenter
A5 = "Absolute spatial error, nx=21"
B5 = SpatialError_nx21
A6 = "Absolute spatial error, nx=41"
B6 = SpatialError_nx41
A7 = "Spatial error ratio (approximately 4)"
B7 = SpatialErrorRatio
A8 = "Observed spatial order (approximately 2)"
B8 = ObservedSpatialOrder
A9 = "Time-integration error against semidiscrete mode, nx=21"
B9 = TemporalError_nx21
A10 = "Time-integration error against semidiscrete mode, nx=41"
B10 = TemporalError_nx41
A12 = Assumptions
A13 = ScenarioGuide

Charts_1!type = "charts"
Charts_1!config = <toml>
  version = 1
  kind = "charts"
  title = "Heat convergence"
  [layout]
  title = "Normalized final slab temperature"
  x_axis = "Position x/L"
  y_axis = "Excess temperature / initial amplitude"
  [[charts]]
  type = "line"
  smooth = false
  bind = "Profile"
  x = "position"
  y = "numerical"
  label = "41 nodes"
  [[charts]]
  type = "line"
  smooth = false
  bind = "Profile"
  x = "position"
  y = "continuum"
  label = "Continuum reference"
</toml>

END MODEL