Skip to content

Silent blow-up when riemann_solver is unspecified: default dflt_int=-100 dispatches no solver, flux_n retains garbage #1410

@sbryngelson

Description

@sbryngelson

Summary

When a case file does not set riemann_solver, the Fortran variable is initialized to dflt_int = -100 (defined in m_global_parameters.fpp). The Fypp dispatch loop in m_riemann_solvers.fpp only handles solvers 1 (HLL), 2 (HLLC), 4 (HLLD), and 5 (LF):

#:for NAME, NUM in [('hll', 1), ('hllc', 2), ('hlld', 4), ('lf', 5)]
    if (riemann_solver == ${NUM}$) then
        call s_${NAME}$_riemann_solver(...)
    end if
#:endfor

With riemann_solver = -100, none of these branches execute. s_initialize_riemann_solver does not zero flux_n (it only conditionally zeros flux_src_vf for viscous/surface tension cases). In a release build (no -finit-real=snan), the flux_n allocatable array retains whatever uninitialized memory the allocator returns. On the very first RHS evaluation the flux divergence uses this garbage, causing immediate catastrophic blow-up (density → 10²⁶¹ after one step).

Reproduction

Add a 2D case that omits riemann_solver from the case dict. The simulation will blow up on step 1 with no diagnostic message.

Root cause

Three separate issues compound:

  1. No solver dispatch for dflt_int: riemann_solver = -100 is not in [1, 2, 4, 5], so flux_n is never written.
  2. flux_n not zeroed at startup: s_initialize_riemann_solver zeroes flux_src_vf only under viscous .or. surface_tension, never flux_n itself.
  3. Case validator silently skips the check: toolchain/mfc/case_validator.py line 663 has if riemann_solver is None: return, allowing the missing parameter to pass validation without a warning.

Impact

Any case file that forgets riemann_solver (including all cases using the default) blows up silently on step 1 with no useful error message. This was the root cause of a blow-up discovered while developing the 2D axisymmetric convergence test.

Suggested fixes

  • Short-term: zero flux_n in s_initialize_riemann_solver (safe, small cost), or add a @:PROHIBIT(riemann_solver == dflt_int, "riemann_solver must be specified") check at startup.
  • Longer-term: give riemann_solver a meaningful default (e.g. 2 for HLLC) in m_global_parameters.fpp, and remove the silent return in case_validator.py.

Metadata

Metadata

Assignees

No one assigned

    Labels

    bugSomething isn't working or doesn't seem right

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions