Persistent temperature oscillations behind shock waves in 1D Castro simulations (simplified‑SDC, Helmholtz EOS, no reactions)

[wangxu250]

​Hello everyone! I am using Castro to simulate 1D shock waves. The equation of state is Helmholtz EOS, and the composition is 50% C12 + 50% O16 (by mass). Reactions are disabled (castro.do_react = 0). The initial conditions are set by computing the pre-shock and post-shock states from the Rankine–Hugoniot relations together with the Helmholtz EOS, which serve as the initial condition for the simulation. I then iterate using the obtained numerical solution.

Behind the shock, the temperature profile exhibits significant unphysical oscillations (spikes). The amplitude of these oscillations decreases as I increase the grid resolution or reduce the CFL number, but due to limited computational resources I cannot continue refining the grid to see if they disappear entirely. Moreover, all flow variables (density, velocity, pressure, temperature) show some oscillations, but temperature is by far the most severely affected. At Mach numbers 1.5, 2.0, 2.5, and 3.0, using 32768 cells in a 1 cm domain (Δx ≈ 0.00003 cm) and CFL = 0.05, the temperature oscillations are still clearly visible. How can I resolve these oscillations with a coarser grid while preserving the physical nature of the solution?

# ------------------  INPUTS TO MAIN PROGRAM  -------------------
max_step = 500000

stop_time    = 1.e-5

# PROBLEM SIZE & GEOMETRY
geometry.is_periodic = 0
geometry.coord_sys   = 0        # 0 => cart, 1 => RZ  2=>spherical
geometry.prob_lo     = 0
geometry.prob_hi     = 1.0
amr.n_cell           = 32768

# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
# 0 = Interior           3 = Symmetry
# 1 = Inflow             4 = SlipWall
# 2 = Outflow            5 = NoSlipWall
# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
castro.lo_bc       =  2
castro.hi_bc       =  2

# WHICH PHYSICS
castro.do_hydro = 1
castro.do_react = 0

castro.small_temp = 1.e6
castro.small_dens = 1.e-5

castro.ppm_type = 1

castro.diffuse_temp = 1
castro.diffuse_cutoff_density = 1.e-2

castro.time_integration_method = 3

# TIME STEP CONTROL
castro.cfl            = 0.05     # cfl number for hyperbolic system
castro.init_shrink    = 0.1     # scale back initial timestep
castro.change_max     = 1.1     # max time step growth


# DIAGNOSTICS & VERBOSITY
castro.sum_interval   = 1       # timesteps between computing mass
castro.v              = 1       # verbosity in Castro.cpp
amr.v                 = 1       # verbosity in Amr.cpp

# REFINEMENT / REGRIDDING
amr.max_level       = 0       # maximum level number allowed
amr.ref_ratio       = 4 4  2 2 2 # refinement ratio
amr.regrid_int      = 2 2 2 2 # how often to regrid
amr.blocking_factor = 64      # block factor in grid generation
amr.max_grid_size   = 64
amr.n_error_buf     = 2 2 2 2 # number of buffer cells in error est

# CHECKPOINT FILES
#amr.restart=chk10000
amr.checkpoint_files_output = 1
amr.check_file      = chk        # root name of checkpoint file
amr.check_int       = 2500       # number of timesteps between checkpoints

# PLOTFILES
amr.plot_file       = plt        # root name of plotfile
amr.plot_per        = 1.e-11
amr.derive_plot_vars = ALL
castro.store_omegadot = 1

# problem initialization

problem.interface_pos=0.15

problem.T_left = 4.9893592120e8
problem.rho_left = 1.8430580379e7
problem.u_left = 2.3973962792e8

problem.T_right = 1.e8
problem.rho_right = 1.e7
problem.u_right = 3.5e3

# refinement

amr.refinement_indicators = tempgrad

amr.refine.tempgrad.max_level = 5
amr.refine.tempgrad.relative_gradient = 0.025
amr.refine.tempgrad.field_name = Temp

# Microphysics

integrator.rtol_spec = 1.e-6
integrator.atol_spec = 1.e-8
integrator.jacobian = 1