Skip to content

Conversation

@tripatmn
Copy link

@tripatmn tripatmn commented Jan 26, 2026

User description

Description

Please include a summary of the changes and the related issue(s) if they exist.
Please also include relevant motivation and context.

Fixes #(issue) [optional]

Type of change

Please delete options that are not relevant.

  • Bug fix (non-breaking change which fixes an issue)
  • New feature (non-breaking change which adds functionality)
  • Something else

Scope

  • This PR comprises a set of related changes with a common goal

If you cannot check the above box, please split your PR into multiple PRs that each have a common goal.

How Has This Been Tested?

Please describe the tests that you ran to verify your changes.
Provide instructions so we can reproduce.
Please also list any relevant details for your test configuration

  • Test A
  • Test B

Test Configuration:

  • What computers and compilers did you use to test this:

Checklist

  • I have added comments for the new code
  • I added Doxygen docstrings to the new code
  • I have made corresponding changes to the documentation (docs/)
  • I have added regression tests to the test suite so that people can verify in the future that the feature is behaving as expected
  • I have added example cases in examples/ that demonstrate my new feature performing as expected.
    They run to completion and demonstrate "interesting physics"
  • I ran ./mfc.sh format before committing my code
  • New and existing tests pass locally with my changes, including with GPU capability enabled (both NVIDIA hardware with NVHPC compilers and AMD hardware with CRAY compilers) and disabled
  • This PR does not introduce any repeated code (it follows the DRY principle)
  • I cannot think of a way to condense this code and reduce any introduced additional line count

If your code changes any code source files (anything in src/simulation)

To make sure the code is performing as expected on GPU devices, I have:

  • Checked that the code compiles using NVHPC compilers
  • Checked that the code compiles using CRAY compilers
  • Ran the code on either V100, A100, or H100 GPUs and ensured the new feature performed as expected (the GPU results match the CPU results)
  • Ran the code on MI200+ GPUs and ensure the new features performed as expected (the GPU results match the CPU results)
  • Enclosed the new feature via nvtx ranges so that they can be identified in profiles
  • Ran a Nsight Systems profile using ./mfc.sh run XXXX --gpu -t simulation --nsys, and have attached the output file (.nsys-rep) and plain text results to this PR
  • Ran a Rocprof Systems profile using ./mfc.sh run XXXX --gpu -t simulation --rsys --hip-trace, and have attached the output file and plain text results to this PR.
  • Ran my code using various numbers of different GPUs (1, 2, and 8, for example) in parallel and made sure that the results scale similarly to what happens if you run without the new code/feature

PR Type

Enhancement, Tests


Description

  • Multiphase chemistry coupling framework: Integrated phase change with chemistry reactions by adding support for gas-phase-only reaction regions in liquid-dominated cells

  • New chemistry parameters: Extended chemistry_parameters derived type with multiphase, liquid_phase_idx, fuel_species_idx, and gas_phase_threshold fields for controlling multiphase chemistry behavior

  • Density safeguards: Implemented density floor enforcement (1.0e-12) and zero-density checks across chemistry, variable conversion, and phase change modules to prevent division errors

  • Evaporated mass tracking: Integrated automatic transfer of evaporated mass to fuel species in phase change module for proper chemistry coupling

  • Input validation: Added comprehensive validation subroutine s_check_inputs_multiphase_chemistry to verify parameter constraints and phase relaxation requirements

  • MPI synchronization: Extended MPI broadcast functionality in both pre-process and simulation modules to communicate multiphase chemistry parameters across all ranks

  • Comprehensive examples and documentation: Added 2D burning droplet example with multiple test cases, visualization utilities, validation framework, and detailed design/validation documentation

  • Test cases: Implemented Phase 1 validation tests covering chemistry skipping, evaporation mass transfer, input validation, and conservation checks


Diagram Walkthrough

flowchart LR
  A["Input Parameters<br/>multiphase, liquid_phase_idx,<br/>fuel_species_idx, gas_phase_threshold"]
  B["MPI Broadcast<br/>Synchronize across ranks"]
  C["Input Validation<br/>Check constraints"]
  D["Chemistry Module<br/>Skip in liquid cells<br/>Track evaporated mass"]
  E["Phase Change Module<br/>Transfer evaporated mass<br/>to fuel species"]
  F["Variable Conversion<br/>Apply density floor<br/>Handle multiphase EOS"]
  G["Simulation Results<br/>Coupled phase change<br/>+ chemistry reactions"]
  
  A --> B
  B --> C
  C --> D
  D --> E
  E --> F
  F --> G
Loading

File Walkthrough

Relevant files
Enhancement
7 files
m_chemistry.fpp
Multiphase chemistry coupling with phase change integration

src/common/m_chemistry.fpp

  • Added multiphase chemistry support with gas-phase-only reaction
    regions
  • Implemented density calculations handling both single-fluid and
    multiphase cases
  • Added safety checks for zero/near-zero densities to prevent division
    errors
  • Added liquid cell detection and skipping logic using volume fraction
    thresholds
  • Integrated evaporated mass tracking into species equations for phase
    change coupling
+110/-14
m_phase_change.fpp
Phase change and chemistry coupling integration                   

src/common/m_phase_change.fpp

  • Added tracking of initial liquid mass for evaporation calculations
  • Implemented automatic transfer of evaporated mass to fuel species
  • Added conditional logic for multiphase chemistry coupling in phase
    change module
+20/-1   
m_mpi_proxy.fpp
MPI communication for multiphase chemistry parameters       

src/pre_process/m_mpi_proxy.fpp

  • Added MPI broadcast for new multiphase chemistry parameters
  • Broadcasts multiphase, liquid_phase_idx, fuel_species_idx flags
  • Broadcasts gas_phase_threshold real parameter across all MPI ranks
+11/-0   
m_mpi_proxy.fpp
MPI synchronization of multiphase chemistry settings         

src/simulation/m_mpi_proxy.fpp

  • Extended MPI broadcast to include multiphase logical parameter
  • Added broadcasting of liquid_phase_idx and fuel_species_idx integers
  • Added broadcasting of gas_phase_threshold real parameter
+4/-2     
m_start_up.fpp
Chemistry parameters input file support                                   

src/pre_process/m_start_up.fpp

  • Added chem_params to the namelist read statement for preprocessing
  • Enables reading of multiphase chemistry parameters from input files
+1/-1     
m_derived_types.fpp
Chemistry parameters type extension for multiphase support

src/common/m_derived_types.fpp

  • Added four new fields to chemistry_parameters derived type
  • multiphase: logical flag to enable multiphase chemistry coupling
  • liquid_phase_idx: integer index of liquid phase fluid
  • fuel_species_idx: integer index of fuel species in mechanism
  • gas_phase_threshold: real parameter for minimum gas volume fraction
+7/-0     
viz.py
Visualization script for burning droplet results                 

examples/2D_burning_droplet/viz.py

  • New visualization script for burning droplet simulation results
  • Reads MFC binary data files and creates plots of temperature, species
    mass fractions, and flame structure
  • Provides command-line interface for selecting timesteps and output
    directories
  • Includes two main plotting functions: plot_burning_droplet() for 2D
    field visualization and plot_flame_structure() for radial profiles
+185/-0 
Error handling
2 files
m_variables_conversion.fpp
Density safeguards and multiphase EOS handling                     

src/common/m_variables_conversion.fpp

  • Added density floor enforcement (1.0e-12) for multiphase chemistry in
    liquid cells
  • Implemented conditional logic to use standard EOS in liquid-dominated
    cells
  • Added chemistry-specific energy calculations only for gas-phase
    regions
  • Ensured non-zero density before division operations in mixture
    calculations
+53/-19 
m_checker.fpp
Input validation for multiphase chemistry parameters         

src/simulation/m_checker.fpp

  • Added new validation subroutine s_check_inputs_multiphase_chemistry
  • Validates multiphase chemistry parameter constraints and ranges
  • Checks that phase relaxation is enabled when using multiphase
    chemistry
  • Verifies liquid phase index and fuel species index are within valid
    ranges
+17/-0   
Configuration changes
3 files
m_global_parameters.fpp
Default initialization of multiphase chemistry parameters

src/pre_process/m_global_parameters.fpp

  • Initialized new multiphase chemistry parameters with default values
  • multiphase = .false., liquid_phase_idx = 1, fuel_species_idx = 1
  • gas_phase_threshold = 0.01_wp (1% minimum gas fraction)
+4/-0     
m_global_parameters.fpp
Simulation module multiphase chemistry initialization       

src/simulation/m_global_parameters.fpp

  • Added initialization of multiphase chemistry parameters in simulation
    module
  • Sets same defaults as preprocessing module for consistency
+4/-0     
case_dicts.py
Case dictionary parameter registration for multiphase chemistry

toolchain/mfc/run/case_dicts.py

  • Added parameter type definitions for new multiphase chemistry
    parameters
  • Registered multiphase, liquid_phase_idx, fuel_species_idx in both
    PRE_PROCESS and SIMULATION
  • Added gas_phase_threshold as REAL parameter type for both modules
+7/-2     
Documentation
8 files
case.py
2D burning droplet example with gas-phase chemistry           

examples/2D_burning_droplet/case.py

  • New example case for 2D burning droplet simulation with H2-O2
    chemistry
  • Implements smooth fuel-oxidizer transition using analytical
    expressions
  • Configures single-fluid model with species diffusion and reactions
  • Includes comprehensive documentation and parameter explanations
+352/-0 
visualize_phase1.py
Phase 1 validation data visualization script                         

examples/2D_burning_droplet/visualize_phase1.py

  • Visualization utility for Phase 1 validation restart data
  • Plots volume fractions, partial densities, and time evolution
  • Includes interface tracking and grid reading functionality
+318/-0 
case_liquid_droplet.py
Experimental liquid droplet combustion case template         

examples/2D_burning_droplet/case_liquid_droplet.py

  • Experimental case combining phase change with chemistry (3-fluid
    model)
  • Sets up liquid fuel, vapor, and oxidizer phases
  • Includes commented-out multiphase chemistry coupling configuration
  • Documents limitations and alternative approaches
+306/-0 
visualize_chemistry.py
Chemistry simulation results visualization utility             

examples/2D_burning_droplet/visualize_chemistry.py

  • Visualization script for chemistry test results
  • Plots species profiles, reactants, products, and energy evolution
  • Includes time evolution tracking for key species
+247/-0 
case_hydrocarbon.py
Hydrocarbon droplet combustion case template                         

examples/2D_burning_droplet/case_hydrocarbon.py

  • Template case for hydrocarbon (methane) droplet combustion using GRI30
    mechanism
  • Demonstrates setup for more complex chemistry mechanisms
  • Includes documentation on species configuration and mechanism file
    usage
  • Provides guidance for adapting to other fuel types
+253/-0 
COUPLING_DESIGN.md
Design document for phase change-chemistry coupling           

examples/2D_burning_droplet/COUPLING_DESIGN.md

  • Comprehensive design document for coupling phase change with chemistry
    in MFC
  • Describes current architecture limitations and proposed coupling
    approach
  • Details required code modifications across 8 major components (global
    parameters, pre-processor, chemistry, phase change, etc.)
  • Provides implementation priority with Phase 1 marked as completed
+313/-0 
README.md
Documentation for burning droplet simulation example         

examples/2D_burning_droplet/README.md

  • Comprehensive documentation for 2D burning droplet simulation example
  • Explains physics of droplet combustion and current limitations of MFC
    modules
  • Provides usage instructions, configuration details, and expected
    results
  • Documents workarounds for combining vaporization with combustion and
    references to multiphase chemistry coupling feature
+278/-0 
VALIDATION_RESULTS.md
Phase 1 validation results and implementation summary       

examples/2D_burning_droplet/VALIDATION_RESULTS.md

  • Summary of Phase 1 implementation completion and validation results
  • Documents completed components across chemistry, phase change, and
    variable conversion modules
  • Lists test results showing gas-only chemistry and phase change work,
    combined case has numerical challenges
  • Provides parameter reference, example usage, and list of modified
    files
+83/-0   
Tests
5 files
test_phase1_validation.py
Phase 1 multiphase chemistry validation test case               

examples/2D_burning_droplet/test_phase1_validation.py

  • Validation test case for Phase 1 multiphase chemistry implementation
  • 1D domain with liquid-gas interface for testing phase change coupling
  • Configurable chemistry and multiphase parameters for regression
    testing
  • Includes detailed comments on expected behavior
+299/-0 
test_phase_change_only.py
Phase change-only validation test case                                     

examples/2D_burning_droplet/test_phase_change_only.py

  • Test case for phase change without chemistry reactions
  • Validates basic multiphase flow with relaxation model
  • Provides baseline for comparing with chemistry-enabled runs
+299/-0 
validate_results.py
Validation checker for multiphase chemistry coupling         

examples/2D_burning_droplet/validate_results.py

  • New validation script for Phase 1 multiphase chemistry coupling
  • Provides framework for checking simulation outputs (NaN/Inf values,
    conservation laws, volume fractions)
  • Includes validation checklist and expected outcomes for six test cases
  • Generates summary report of validation status and next steps
+227/-0 
test_chemistry_only.py
Simple gas-phase chemistry test case                                         

examples/2D_burning_droplet/test_chemistry_only.py

  • New test case for gas-phase H2-O2 chemistry without phase change
  • Configures 1D domain with fuel (H2) on left and oxidizer (O2+N2) on
    right
  • Demonstrates pure combustion chemistry with diffusion and reactions
  • Includes command-line arguments for temperature and reaction control
+204/-0 
VALIDATION_PHASE1.md
Phase 1 validation test plan and specifications                   

examples/2D_burning_droplet/VALIDATION_PHASE1.md

  • Detailed validation test plan for Phase 1 multiphase chemistry
    coupling
  • Defines six test cases covering chemistry skipping, evaporation mass
    transfer, input validation, boundary conditions, conservation, and
    threshold sensitivity
  • Specifies quantitative metrics (mass conservation error < 1e-10) and
    qualitative checks
  • Provides success/failure criteria and test case parameters
+257/-0 

Summary by CodeRabbit

  • New Features

    • Multiphase chemistry coupling for droplet combustion simulations: chemistry now runs only in gas-phase regions while evaporated liquid mass transfers to fuel species.
    • New configuration parameters for multiphase modeling (liquid phase index, fuel species index, gas-phase threshold).
  • Documentation

    • Comprehensive 2D burning droplet example with design guide, validation tests, and usage instructions.
    • New visualization tools and test cases for chemistry-only, phase-change-only, and coupled multiphase scenarios.

✏️ Tip: You can customize this high-level summary in your review settings.

cursoragent and others added 13 commits January 24, 2026 22:18
This example demonstrates droplet combustion using MFC's chemistry module:
- Fuel vapor 'droplet' surrounded by oxidizer (O2 + AR)
- Species diffusion for fuel-oxidizer mixing
- Chemical reactions using h2o2.yaml mechanism
- Unity-Lewis transport model for stability

Includes:
- case.py: Main simulation case file
- README.md: Documentation and physics explanation
- viz.py: Visualization script for post-processing

This is a gas-phase combustion model where the droplet is
pre-vaporized, suitable for studying flame dynamics and
species transport in a droplet-like configuration.

Co-authored-by: tripatmn <tripatmn@mail.uc.edu>
Key improvements:
- Use analytical tanh profiles for smooth fuel-oxidizer transition
- Add --fast mode for quick testing (coarser grid, shorter time)
- Add --T_droplet option to control ignition temperature
- Better domain size (10mm) for realistic flame standoff
- Improved README with physics explanations and future work

The smooth transition profile mimics a diffusion flame structure
with realistic mixing layers between fuel and oxidizer regions.

Co-authored-by: tripatmn <tripatmn@mail.uc.edu>
Added:
- case_hydrocarbon.py: Template for CH4/GRI30 combustion
- Expanded README with:
  - File descriptions
  - Available mechanisms table
  - Guidance on combining phase change + chemistry
  - Troubleshooting section
  - Workarounds for current limitations

Co-authored-by: tripatmn <tripatmn@mail.uc.edu>
New files:
- case_liquid_droplet.py: True liquid droplet with phase change (3-fluid model)
  Uses relax_model=6 (pTg-equilibrium) for evaporation
  Based on 2D_phasechange_bubble parameters

Documentation improvements:
- Clear explanation of phase change vs chemistry limitation
- Table comparing the two approaches
- Guidance on which case file to use
- Workarounds for true burning droplet simulation
- Two-stage simulation approach explained

The current MFC architecture separates:
- Phase change: volume fractions (alpha), num_fluids > 1
- Chemistry: mass fractions (Y), num_fluids = 1

Coupling these requires future development.

Co-authored-by: tripatmn <tripatmn@mail.uc.edu>
COUPLING_DESIGN.md documents:
- Current architecture limitations
- Proposed equation system extension
- Required code modifications (9 files, ~880 lines)
- Implementation phases and testing strategy

Key insight: species should be tracked within gas phase only,
with evaporation providing source terms to fuel species.

Estimated modifications:
- m_global_parameters.fpp: Add multiphase_chemistry flag
- m_chemistry.fpp: Add gas-phase checks, evaporation sources
- m_phase_change.fpp: Add species source terms for evaporation
- m_riemann_solvers.fpp: Gas-phase weighted species fluxes
- m_variables_conversion.fpp: Gas-phase species conversion
- Pre-processor files: Species initialization in multiphase
- case_validator.py: New parameter validation

Co-authored-by: tripatmn <tripatmn@mail.uc.edu>
…droplet simulation

Add multiphase chemistry capability to enable chemistry calculations only in gas-phase
regions and couple evaporation with species transport. This is the first phase of
the burning droplet simulation feature.

Changes:
- Add multiphase chemistry parameters to chemistry_parameters type:
  - chem_params%multiphase: Enable/disable multiphase chemistry coupling
  - chem_params%liquid_phase_idx: Index of liquid phase fluid (default: 1)
  - chem_params%fuel_species_idx: Index of fuel species in mechanism (default: 1)
  - chem_params%gas_phase_threshold: Min gas volume fraction for chemistry (default: 0.01)

- Modify chemistry reaction flux to skip liquid-dominated cells when
  multiphase chemistry is enabled (m_chemistry.fpp)

- Add evaporation source term to species equations in phase change module:
  When liquid evaporates, the mass is added to the fuel species in the
  chemistry system (m_phase_change.fpp)

- Add MPI broadcast for new chemistry parameters (m_mpi_proxy.fpp)

- Add validation checks for multiphase chemistry configuration (m_checker.fpp):
  - Requires relax = T for phase change
  - Requires num_fluids >= 2
  - Validates liquid_phase_idx, fuel_species_idx, gas_phase_threshold ranges

- Update pre-processor global parameters with defaults (m_global_parameters.fpp)

Co-authored-by: tripatmn <tripatmn@mail.uc.edu>
- Mark Phase 1 as completed in COUPLING_DESIGN.md
- Add usage examples for new chem_params parameters
- Update case_liquid_droplet.py with multiphase chemistry configuration
- Update README.md with new multiphase chemistry coupling section

Co-authored-by: tripatmn <tripatmn@mail.uc.edu>
Add comprehensive validation suite for multiphase chemistry coupling:

- VALIDATION_PHASE1.md: Expected outcomes and test matrix
- test_phase1_validation.py: Minimal 1D test case
- validate_results.py: Results analysis script
- VALIDATION_RESULTS.md: Actual test results

Also update toolchain/mfc/run/case_dicts.py to recognize new
chem_params parameters in case files.

Results:
- Code compiles successfully
- Pre-processor works correctly
- Simulation starts with multiphase chemistry
- Test case needs numerical tuning (NaN after ~10 steps)
  due to stiff physics, not implementation bug

Co-authored-by: tripatmn <tripatmn@mail.uc.edu>
- Add test_phase_change_only.py for testing phase change without chemistry
- Add visualize_phase1.py for plotting simulation results
- Add generated figures showing:
  - Volume fractions at t=0, 50, 100
  - Partial densities at t=0, 50, 100
  - Time evolution of key quantities
  - Interface position tracking

Phase change simulation runs successfully for 100 time steps.
Chemistry case needs further tuning (NaN at t=0).

Co-authored-by: tripatmn <tripatmn@mail.uc.edu>
- test_chemistry_only.py: Single-fluid H2-O2 chemistry test
- visualize_chemistry.py: Visualization for chemistry results
- Generated figures showing species evolution

Gas-only chemistry runs successfully for 100 steps.
Shows H2O production at fuel-air interface.

Co-authored-by: tripatmn <tripatmn@mail.uc.edu>
… MPI broadcast

- Add chem_params to pre_process namelist (m_start_up.fpp)
- Add MPI broadcast for chem_params in pre_process (m_mpi_proxy.fpp)
- Add chem_params to PRE_PROCESS dict in case_dicts.py
- Fix primitive to conservative conversion for multiphase chemistry
- Fix temperature computation in chemistry module for multiphase
- Skip chemistry diffusion at liquid interfaces

Pre_process now generates NaN-free initial conditions for multiphase
chemistry cases. The simulation still crashes during initialization -
further investigation needed.

Co-authored-by: tripatmn <tripatmn@mail.uc.edu>
- Add multiphase chemistry check to primitive-to-flux conversion
- Skip chemistry EOS in liquid-dominated cells in flux computation
- Use standard stiffened gas EOS for liquid cells

Co-authored-by: tripatmn <tripatmn@mail.uc.edu>
- Restore test_phase1_validation.py settings
- Update VALIDATION_RESULTS.md with detailed findings:
  - Chemistry-only and phase-change-only tests pass
  - Combined test has numerical challenges at interface
  - Document known issues and recommended next steps
- List all modified files and parameter reference

Co-authored-by: tripatmn <tripatmn@mail.uc.edu>
@tripatmn tripatmn requested a review from sbryngelson as a code owner January 26, 2026 23:28
@codeant-ai
Copy link
Contributor

codeant-ai bot commented Jan 26, 2026

CodeAnt AI is reviewing your PR.


Thanks for using CodeAnt! 🎉

We're free for open-source projects. if you're enjoying it, help us grow by sharing.

Share on X ·
Reddit ·
LinkedIn

@coderabbitai
Copy link
Contributor

coderabbitai bot commented Jan 26, 2026

Caution

Review failed

The pull request is closed.

📝 Walkthrough

Walkthrough

Introduces Phase 1 multiphase chemistry coupling for a 2D burning droplet example in MFC. Adds design documents, configuration cases, test scripts, and visualization utilities alongside core logic changes: chemistry gating in liquid-dominated regions, evaporated-mass transfer to fuel species during phase-change, and liquid-cell handling in variable conversions. Infrastructure updates include MPI broadcasting, validation, and parameter initialization.

Changes

Cohort / File(s) Summary
Documentation & Design
examples/2D_burning_droplet/COUPLING_DESIGN.md, examples/2D_burning_droplet/README.md, examples/2D_burning_droplet/VALIDATION_PHASE1.md, examples/2D_burning_droplet/VALIDATION_RESULTS.md
Design document outlining Phase 1 implementation (multiphase chemistry gating, evaporated-mass transfer, parameter validation). README covering physics, configuration, running instructions, and phase-change + combustion workarounds. Validation specifications and implementation status for six test cases (chemistry skipping, mass transfer, input validation, boundary conditions, conservation, threshold sensitivity).
Example & Test Cases
examples/2D_burning_droplet/case.py, examples/2D_burning_droplet/case_hydrocarbon.py, examples/2D_burning_droplet/case_liquid_droplet.py, examples/2D_burning_droplet/test_chemistry_only.py, examples/2D_burning_droplet/test_phase1_validation.py, examples/2D_burning_droplet/test_phase_change_only.py
Python scripts defining H2-O2 burning droplet, GRI30 hydrocarbon, and 3-fluid liquid-droplet cases with configurable chemistry/diffusion flags. Test cases for chemistry-only, phase-change-only, and Phase 1 validation scenarios with JSON output. Support CLI options for temperature, scaling, and test modes.
Validation & Visualization
examples/2D_burning_droplet/validate_results.py, examples/2D_burning_droplet/visualize_chemistry.py, examples/2D_burning_droplet/visualize_phase1.py, examples/2D_burning_droplet/viz.py
Result validator scaffold with NaN/Inf and range checks, output directory discovery, and formatted validation report. Visualization utilities for chemistry data (species profiles, time-evolution), Phase 1 fields (volume/partial fractions, interface tracking), and 2D burning droplet (temperature/species fields, flame structure profiles).
Chemistry Module
src/common/m_chemistry.fpp
Adds multiphase gating to skip reaction and diffusion computation in liquid-dominated cells (alpha_gas < threshold). Computes rho_total from all partial densities, clamps species mass fractions to [0,1], and enforces non-zero density in diffusion calculations (+110/-14 lines).
Derived Types & Parameters
src/common/m_derived_types.fpp
Extends chemistry_parameters type with four new fields: multiphase (logical), liquid_phase_idx (integer), fuel_species_idx (integer), gas_phase_threshold (real).
Phase-Change Module
src/common/m_phase_change.fpp
Computes evaporated mass (dm_evap) and transfers it to fuel species when multiphase chemistry is enabled. Stores initial liquid mass (m1_old) and increments fuel species mass by evaporated amount.
Variable Conversions
src/common/m_variables_conversion.fpp
Adds liquid-dominated cell handling in conservative ↔ primitive and primitive → flux conversions. Skips chemistry EOS for liquid cells, uses standard EOS with zeroed species terms; otherwise retains chemistry path. Clamps density to 1e-12 minimum (+53/-19 lines).
Global Parameters & Initialization
src/pre_process/m_global_parameters.fpp, src/simulation/m_global_parameters.fpp
Sets default values for new multiphase chemistry parameters: multiphase = .false., liquid_phase_idx = 1, fuel_species_idx = 1, gas_phase_threshold = 0.01.
MPI Broadcasting
src/pre_process/m_mpi_proxy.fpp, src/simulation/m_mpi_proxy.fpp
Broadcasts new chemistry parameters (multiphase logical, liquid/fuel indices, gas-phase threshold) from rank 0 to all ranks during user-input phase.
Infrastructure & Validation
src/pre_process/m_start_up.fpp, src/simulation/m_checker.fpp, toolchain/mfc/run/case_dicts.py
Adds chem_params to pre_process namelist. New validation subroutine s_check_inputs_multiphase_chemistry enforces constraints (relax enabled, num_fluids ≥ 2, valid indices, 0 ≤ threshold ≤ 1). Toolchain parameter definitions for multiphase, liquid_phase_idx, fuel_species_idx, gas_phase_threshold.

Sequence Diagram(s)

sequenceDiagram
    participant PC as Phase-Change<br/>(s_infinite_<br/>relaxation)
    participant VConv as Variables<br/>Conversion
    participant Chem as Chemistry<br/>(Reaction/Diffusion)
    
    Note over PC,Chem: Time Step Advancement
    
    PC->>PC: Compute liquid mass (m1)
    PC->>PC: Store m1_old
    PC->>PC: Relax to equilibrium<br/>(evaporate liquid)
    PC->>PC: Compute evaporated mass<br/>dm_evap = m1_old - m1
    
    alt Multiphase Chemistry Enabled
        PC->>PC: Transfer dm_evap<br/>to fuel_species_idx
    end
    
    PC->>VConv: Updated conserved<br/>variables (q_cons)
    
    VConv->>VConv: Compute alpha_liquid<br/>from volume fractions
    
    alt Liquid-Dominated Cell
        VConv->>VConv: Use standard EOS<br/>(skip chemistry)
        VConv->>VConv: Zero species energy
    else Gas Cell
        VConv->>VConv: Use chemistry EOS<br/>compute Ys, T, e_mix
    end
    
    VConv->>Chem: Converted primitives<br/>(q_prim)
    
    Chem->>Chem: Compute alpha_gas<br/>from volume fractions
    
    alt alpha_gas > threshold
        Chem->>Chem: Compute diffusion flux
        Chem->>Chem: Compute reaction flux
    else alpha_gas ≤ threshold
        Chem->>Chem: Skip chemistry<br/>(return zero flux)
    end
Loading

Estimated code review effort

🎯 4 (Complex) | ⏱️ ~50 minutes

Possibly related PRs

  • PR #1084: Modifies src/common/m_chemistry.fpp and chemistry_parameters type; overlaps with transport_model parameter additions in diffusion/reaction loop logic and derived-type extensions.

Suggested labels

Review effort 4/5, size:L

Suggested reviewers

  • sbryngelson

Poem

🐰 A droplet burns with vapor's dance,
Chemistry and phase-change advance,
Liquid gates now guard the flame,
Multiphase coupling—what a game!
Evaporation feeds the fire,
As fuel species rise ever higher! 🔥

✨ Finishing touches
  • 📝 Generate docstrings

Thanks for using CodeRabbit! It's free for OSS, and your support helps us grow. If you like it, consider giving us a shout-out.

❤️ Share

Comment @coderabbitai help to get the list of available commands and usage tips.

@codeant-ai codeant-ai bot added the size:XXL This PR changes 1000+ lines, ignoring generated files label Jan 26, 2026
@qodo-code-review
Copy link
Contributor

PR Reviewer Guide 🔍

Here are some key observations to aid the review process:

⏱️ Estimated effort to review: 4 🔵🔵🔵🔵⚪
🧪 PR contains tests
🔒 No security concerns identified
⚡ Recommended focus areas for review

Possible Issue

The MPI broadcast type used for chem_params%gas_phase_threshold differs between pre-process and simulation proxies. Pre-process uses mpi_p while simulation uses MPI_DOUBLE_PRECISION. If mpi_p is not guaranteed to be double precision across builds, this can lead to corrupted gas_phase_threshold on non-root ranks and hard-to-debug behavior.

! Chemistry parameters (for multiphase chemistry coupling)
if (chemistry) then
    #:for VAR in [ 'diffusion', 'reactions', 'multiphase' ]
        call MPI_BCAST(chem_params%${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
    #:endfor
    #:for VAR in [ 'gamma_method', 'transport_model', 'liquid_phase_idx', 'fuel_species_idx' ]
        call MPI_BCAST(chem_params%${VAR}$, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
    #:endfor
    call MPI_BCAST(chem_params%gas_phase_threshold, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)
end if
Possible Issue

Evaporated mass is added directly to the fuel species conservative variable, but there is no visible corresponding adjustment to ensure overall mass/species consistency (e.g., whether this double-counts with existing vapor-phase continuity updates, or whether the fuel species variable is expected to represent the evaporating fluid). This needs validation for conservation and for correct coupling semantics between fluid mass transfer and species mass transfer.

! Multiphase chemistry coupling: Add evaporated mass to fuel species
! When liquid evaporates, the mass becomes fuel vapor which should be
! tracked by the species equations
if (chemistry .and. chem_params%multiphase) then
    ! Compute evaporated mass (positive when liquid evaporates)
    dm_evap = m1_old - q_cons_vf(lp + contxb - 1)%sf(j, k, l)

    ! Add evaporated mass to fuel species
    ! The fuel species is at index chem_params%fuel_species_idx
    if (dm_evap > 0.0_wp) then
        q_cons_vf(chemxb + chem_params%fuel_species_idx - 1)%sf(j, k, l) = &
            q_cons_vf(chemxb + chem_params%fuel_species_idx - 1)%sf(j, k, l) + dm_evap
    end if
end if
Duplication / Consistency

There are multiple duplicated “skip liquid-dominated cells” checks and threshold comparisons across temperature evaluation, reactions, and diffusion paths, using slightly different forms (e.g., alpha_gas < gas_phase_threshold vs liquid fraction > (1 - gas_phase_threshold)). This increases the risk of inconsistent behavior at the interface and makes future changes error-prone; consider centralizing the gas/liquid eligibility logic (or at least making the condition identical everywhere) and confirming the advection index (advxb + liquid_phase_idx - 1) is consistent for both conservative and primitive representations.

! For multiphase chemistry, handle liquid cells specially
if (chem_params%multiphase) then
    ! Get gas volume fraction from advection variables
    alpha_gas = 1.0_wp - q_cons_vf(advxb + chem_params%liquid_phase_idx - 1)%sf(x, y, z)

    ! Skip liquid cells - set default temperature
    if (alpha_gas < chem_params%gas_phase_threshold) then
        q_T_sf%sf(x, y, z) = 300.0_wp
        cycle
    end if

    ! Compute total density from all partial densities
    rho_total = 0.0_wp
    do i = contxb, contxe
        rho_total = rho_total + q_cons_vf(i)%sf(x, y, z)
    end do

    ! Ensure non-zero density
    if (rho_total < 1.0e-10_wp) then
        q_T_sf%sf(x, y, z) = 300.0_wp
        cycle
    end if

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

High-level Suggestion

To improve numerical stability at the liquid-gas interface, replace the sharp gas_phase_threshold with a smooth blending function. This will create a gradual transition between the multiphase and chemistry models. [High-level, importance: 9]

Solution Walkthrough:

Before:

subroutine s_compute_chemistry_reaction_flux(...)
    ...
    do x = bounds(1)%beg, bounds(1)%end
        if (chem_params%multiphase) then
            alpha_gas = 1.0_wp - q_prim_qp(advxb + chem_params%liquid_phase_idx - 1)%sf(x, y, z)

            ! Sharp cutoff for chemistry
            if (alpha_gas < chem_params%gas_phase_threshold) then
                cycle ! Skip chemistry for this cell
            end if
        end if

        ! Compute reaction rates (omega_m)
        ...
        rhs_vf(eqn)%sf(x, y, z) = rhs_vf(eqn)%sf(x, y, z) + omega_m
    end do
    ...
end subroutine

After:

subroutine s_compute_chemistry_reaction_flux(...)
    ...
    do x = bounds(1)%beg, bounds(1)%end
        blending_factor = 1.0_wp
        if (chem_params%multiphase) then
            alpha_gas = 1.0_wp - q_prim_qp(advxb + chem_params%liquid_phase_idx - 1)%sf(x, y, z)

            ! Smooth blending function (e.g., smoothstep)
            ! It transitions from 0 to 1 around the threshold.
            x_norm = (alpha_gas - threshold_min) / (threshold_max - threshold_min)
            x_norm = max(0.0_wp, min(1.0_wp, x_norm))
            blending_factor = x_norm*x_norm * (3.0_wp - 2.0_wp * x_norm)
        end if

        ! Compute reaction rates (omega_m)
        ...
        ! Scale the reaction source term by the blending factor
        rhs_vf(eqn)%sf(x, y, z) = rhs_vf(eqn)%sf(x, y, z) + blending_factor * omega_m
    end do
    ...
end subroutine

Comment on lines +323 to +324
rho_L = max(q_prim_qp(1)%sf(x, y, z), 1.0e-12_wp)
rho_R = max(q_prim_qp(1)%sf(x + offsets(1), y + offsets(2), z + offsets(3)), 1.0e-12_wp)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggestion: In multiphase chemistry mode, calculate total gas-phase density by summing all non-liquid fluid densities instead of assuming fluid 1 is the gas phase. [possible issue, importance: 9]

Suggested change
rho_L = max(q_prim_qp(1)%sf(x, y, z), 1.0e-12_wp)
rho_R = max(q_prim_qp(1)%sf(x + offsets(1), y + offsets(2), z + offsets(3)), 1.0e-12_wp)
if (chem_params%multiphase) then
rho_L = 0.0_wp
rho_R = 0.0_wp
do i = 1, num_fluids
if (i /= chem_params%liquid_phase_idx) then
rho_L = rho_L + q_prim_qp(i)%sf(x, y, z)
rho_R = rho_R + q_prim_qp(i)%sf(x + offsets(1), y + offsets(2), z + offsets(3))
end if
end do
else
rho_L = q_prim_qp(1)%sf(x, y, z)
rho_R = q_prim_qp(1)%sf(x + offsets(1), y + offsets(2), z + offsets(3))
end if
rho_L = max(rho_L, 1.0e-12_wp)
rho_R = max(rho_R, 1.0e-12_wp)

Comment on lines +276 to +289
! Multiphase chemistry coupling: Add evaporated mass to fuel species
! When liquid evaporates, the mass becomes fuel vapor which should be
! tracked by the species equations
if (chemistry .and. chem_params%multiphase) then
! Compute evaporated mass (positive when liquid evaporates)
dm_evap = m1_old - q_cons_vf(lp + contxb - 1)%sf(j, k, l)

! Add evaporated mass to fuel species
! The fuel species is at index chem_params%fuel_species_idx
if (dm_evap > 0.0_wp) then
q_cons_vf(chemxb + chem_params%fuel_species_idx - 1)%sf(j, k, l) = &
q_cons_vf(chemxb + chem_params%fuel_species_idx - 1)%sf(j, k, l) + dm_evap
end if
end if
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggestion: Modify the phase change logic to handle both evaporation and condensation by adding or subtracting mass from the fuel species, ensuring mass conservation. [possible issue, importance: 8]

Suggested change
! Multiphase chemistry coupling: Add evaporated mass to fuel species
! When liquid evaporates, the mass becomes fuel vapor which should be
! tracked by the species equations
if (chemistry .and. chem_params%multiphase) then
! Compute evaporated mass (positive when liquid evaporates)
dm_evap = m1_old - q_cons_vf(lp + contxb - 1)%sf(j, k, l)
! Add evaporated mass to fuel species
! The fuel species is at index chem_params%fuel_species_idx
if (dm_evap > 0.0_wp) then
q_cons_vf(chemxb + chem_params%fuel_species_idx - 1)%sf(j, k, l) = &
q_cons_vf(chemxb + chem_params%fuel_species_idx - 1)%sf(j, k, l) + dm_evap
end if
end if
if (chemistry .and. chem_params%multiphase) then
! Compute evaporated/condensed mass
! dm_evap > 0 for evaporation, dm_evap < 0 for condensation
dm_evap = m1_old - q_cons_vf(lp + contxb - 1)%sf(j, k, l)
! Add/remove mass from the fuel species to conserve mass
! The fuel species is at index chem_params%fuel_species_idx
q_cons_vf(chemxb + chem_params%fuel_species_idx - 1)%sf(j, k, l) = &
q_cons_vf(chemxb + chem_params%fuel_species_idx - 1)%sf(j, k, l) + dm_evap
end if

Comment on lines +169 to +177
if (chemistry) then
#:for VAR in [ 'diffusion', 'reactions', 'multiphase' ]
call MPI_BCAST(chem_params%${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
#:endfor
#:for VAR in [ 'gamma_method', 'transport_model', 'liquid_phase_idx', 'fuel_species_idx' ]
call MPI_BCAST(chem_params%${VAR}$, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
#:endfor
call MPI_BCAST(chem_params%gas_phase_threshold, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)
end if
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggestion: Use MPI_DOUBLE_PRECISION instead of mpi_p when broadcasting chem_params%gas_phase_threshold to ensure type consistency across MPI ranks. [general, importance: 7]

Suggested change
if (chemistry) then
#:for VAR in [ 'diffusion', 'reactions', 'multiphase' ]
call MPI_BCAST(chem_params%${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
#:endfor
#:for VAR in [ 'gamma_method', 'transport_model', 'liquid_phase_idx', 'fuel_species_idx' ]
call MPI_BCAST(chem_params%${VAR}$, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
#:endfor
call MPI_BCAST(chem_params%gas_phase_threshold, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)
end if
if (chemistry) then
...
call MPI_BCAST(chem_params%gas_phase_threshold, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
...
end if

Comment on lines +171 to +173
matching = [d for d in data_dirs if str(args.timestep) in d]
if matching:
data_dir = matching[0]
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggestion: Fix ambiguous timestep directory matching by parsing the numeric part of the directory name and performing an exact integer comparison instead of using a substring search. [possible issue, importance: 8]

Suggested change
matching = [d for d in data_dirs if str(args.timestep) in d]
if matching:
data_dir = matching[0]
def _get_ts_from_dir(d):
try:
return int(os.path.basename(d).lstrip('D'))
except (ValueError, IndexError):
return -1
matching = [d for d in data_dirs if _get_ts_from_dir(d) == args.timestep]
if matching:
data_dir = matching[0]

Comment on lines +31 to +35
for name in os.listdir(case_dir):
path = os.path.join(case_dir, name)
if os.path.isdir(path) and name.startswith('D'):
dirs.append(path)
return sorted(dirs)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggestion: Implement numerical sorting for timestep directories instead of lexicographical sorting to ensure correct ordering (e.g., D2 before D10). [possible issue, importance: 8]

Suggested change
for name in os.listdir(case_dir):
path = os.path.join(case_dir, name)
if os.path.isdir(path) and name.startswith('D'):
dirs.append(path)
return sorted(dirs)
for name in os.listdir(case_dir):
path = os.path.join(case_dir, name)
if os.path.isdir(path) and name.startswith('D'):
try:
# Extract number to sort numerically
num = int(name.lstrip('D'))
dirs.append((num, path))
except ValueError:
continue # Ignore directories that don't end in a number
# Sort by the numeric part of the directory name
dirs.sort()
return [path for num, path in dirs]

Comment on lines +21 to +26
def read_binary_data(filepath, nx, ny, precision=2):
"""Read MFC binary data file."""
dtype = np.float64 if precision == 2 else np.float32
with open(filepath, 'rb') as f:
data = np.fromfile(f, dtype=dtype)
return data.reshape((ny + 1, nx + 1))
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggestion: Add a try-except block to the read_binary_data function to handle ValueError during data reshaping, preventing crashes from malformed files. [possible issue, importance: 7]

Suggested change
def read_binary_data(filepath, nx, ny, precision=2):
"""Read MFC binary data file."""
dtype = np.float64 if precision == 2 else np.float32
with open(filepath, 'rb') as f:
data = np.fromfile(f, dtype=dtype)
return data.reshape((ny + 1, nx + 1))
def read_binary_data(filepath, nx, ny, precision=2):
"""Read MFC binary data file."""
dtype = np.float64 if precision == 2 else np.float32
with open(filepath, 'rb') as f:
data = np.fromfile(f, dtype=dtype)
try:
return data.reshape((ny + 1, nx + 1))
except ValueError as e:
print(f"Error reshaping {filepath}: {e}")
return None

@codeant-ai
Copy link
Contributor

codeant-ai bot commented Jan 26, 2026

Nitpicks 🔍

🔒 No security issues identified
⚡ Recommended areas for review

  • Mass Conservation
    Evaporated mass (dm_evap) is added to a fuel species variable but there is no explicit update/consistency check of other conserved quantities (total cell mass/density, species normalization, or energy/entropy). Confirm the species slot is a conserved mass (not a mass fraction). If species are stored as mass fractions you'd need to convert mass->fraction (divide by local density) and/or update global conserved totals so the solver remains consistent.

  • Suspicious assignment
    In s_convert_conservative_to_primitive_variables the code sets every continuity primitive slot to the same total density (rho_K) unconditionally. This is likely incorrect: continuity variables are per-fluid partial densities or specific continuity slots and should not all be set to the total mixture density. Verify intended target indices and variable meanings before keeping this assignment.

  • Possible Index OOB
    The code unconditionally writes to q_cons_vf at index (chemxb + chem_params%fuel_species_idx - 1). Validate that chem_params%fuel_species_idx is within the valid species range and that chemxb offset plus this index maps to the expected conserved-species slot. An out-of-range or wrong-slot index can corrupt memory or update the wrong conserved field.

  • Reaction density not multiphase-aware
    s_compute_chemistry_reaction_flux sets rho = q_cons_qp(contxe)%sf(...) without accounting for multiphase logic used elsewhere (gas-phase density vs. total). If multiphase is enabled the reaction rates need gas-phase density; using the raw continuity primitive may be incorrect. Confirm and adapt to use gas-phase density consistent with other routines.

  • Inconsistent volume-fraction source
    Different routines read the liquid/gas volume fraction from different field arrays:

  • s_compute_q_T_sf reads it from q_cons_vf(...)

  • Other routines (s_compute_T_from_primitives, diffusion, reaction) read it from q_prim_vf or q_prim_qp.
    This mismatch can produce inconsistent decisions about whether a cell is "liquid" or "gas" and lead to divergent behavior (skipping chemistry or using defaults incorrectly). Verify which field actually stores the advected liquid fraction and unify reads across the module.

  • Data shape mismatch
    read_binary_data assumes the binary contains exactly (ny+1)*(nx+1) elements and calls reshape directly. If the file size/precision doesn't match, reshape will raise an error or produce corrupted data. Validate file size / element count before reshaping and provide a clear error message or fallback.

  • Binary read robustness
    The code assumes the restart binary layout is exactly float64 contiguous with size num_vars * num_cells and silently infers a new num_vars when sizes match. This can hide endianness, header/footer bytes, wrong dtype, or different ordering (Fortran vs C) and lead to incorrect reshapes or silently wrong data.

  • Binary read robustness
    read_mfc_data opens the file and passes the file object into np.fromfile, which can be unreliable across platforms and may behave differently for file objects vs. filenames. Large files should also consider memory mapping. Confirm this reliably reads the MFC binary layout, handles endianness, and fails with informative errors when format differs.

  • Binary-read / reshape assumptions
    The code uses np.fromfile(..., dtype=np.float64) and then computes num_vars = len(data) // num_cells and reshapes into (num_vars, num_cells). This assumes the restart files are raw float64 arrays with no header and that the length is an exact multiple of num_cells. If the file uses a different dtype, has metadata, or timestep filenames are zero-padded, the reshape will fail or produce incorrect data. Consider validating file format, dtype, and handling unexpected sizes more robustly.

  • Directory selection fragility
    find_data_files sorts directories lexicographically and main selects timesteps by substring matching. This can lead to wrong selection when directories are named like D1, D10, D2 (lexicographic order) or when a substring matches multiple entries. Consider numerical sorting and stricter timestep matching (e.g. exact numeric parse).

@tripatmn tripatmn closed this Jan 26, 2026
@tripatmn tripatmn deleted the muscl-usage-eec40 branch January 26, 2026 23:32

$:GPU_LOOP(parallelism='[seq]')
do i = 1, contxe
qK_prim_vf(i)%sf(j, k, l) = rho_K
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggestion: Logic bug: the new loop overwrites all primitive continuity variables with the total density (rho_K) which corrupts primitive state (per-species partial densities / continuity variables). This will make downstream code use incorrect primitive continuity values and break conversions (mass fractions, fluxes, etc.). Restore copying from the conservative state (or the previous correct mapping) instead of writing the total density to every continuity primitive. [logic error]

Severity Level: Critical 🚨
- ❌ Conservative-to-primitive conversion corrupted for chem cells.
- ❌ Flux/species advection calculations use wrong densities.
- ⚠️ Energy and mass diagnostics may be inaccurate.
Suggested change
qK_prim_vf(i)%sf(j, k, l) = rho_K
qK_prim_vf(i)%sf(j, k, l) = qK_cons_vf(i)%sf(j, k, l)
Steps of Reproduction ✅
1. Open src/common/m_variables_conversion.fpp and locate subroutine
s_convert_conservative_to_primitive_variables. The added block appears in the chemistry
branch where rho_K is computed (see lines around 719-730 in the PR hunk).

2. Run a simulation with chemistry enabled such that chem_params%multiphase is true and
you have a liquid-dominated cell where species conservative masses sum to zero (so rho_K
computed by the loop at line ~719 is near zero). This situation occurs inside
s_convert_conservative_to_primitive_variables when chemistry is true (the code path that
computes rho_K from qK_cons_vf at lines ~719-724).

3. In that case the code at lines 722-726 forces rho_K to 1.0e-12_wp, then the loop at
lines 728-730 writes rho_K into every primitive continuity slot qK_prim_vf(i)%sf(j,k,l).
This overwrites per-species primitive continuity values with the total density instead of
preserving per-species continuity.

4. Downstream routines that immediately read qK_prim_vf continuity fields (for example
s_convert_primitive_to_flux_variables in the same file which reads qK_prim_vf into
alpha_rho_K at lines ~1263-1276) will therefore receive incorrect
continuity/partial-density values. The observed effect is corrupted primitive state and
incorrect species/flux calculations. The code pattern (writing rho_K into all contxe
slots) appears unintended compared to the surrounding conversions (which normally copy
from qK_cons_vf or compute per-species values).
Prompt for AI Agent 🤖
This is a comment left during a code review.

**Path:** src/common/m_variables_conversion.fpp
**Line:** 730:730
**Comment:**
	*Logic Error: Logic bug: the new loop overwrites all primitive continuity variables with the total density (`rho_K`) which corrupts primitive state (per-species partial densities / continuity variables). This will make downstream code use incorrect primitive continuity values and break conversions (mass fractions, fluxes, etc.). Restore copying from the conservative state (or the previous correct mapping) instead of writing the total density to every continuity primitive.

Validate the correctness of the flagged issue. If correct, How can I resolve this? If you propose a fix, implement it and please make it concise.

formatter_class=argparse.ArgumentDefaultsHelpFormatter
)

parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT",
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggestion: The argparse --mfc argument sets default="{}" (a string) while type=json.loads parses CLI values to a dict; argparse does not apply type to the default value, so when the option is omitted args.mfc will be the string "{}" and later args.mfc.get(...) will raise an AttributeError. Change the default to an actual dict (e.g. {}) so args.mfc is a mapping when not provided. [type error]

Severity Level: Critical 🚨
- ❌ Example crashes at startup with AttributeError.
- ⚠️ Example unusable without providing --mfc manually.
- ⚠️ CI checks running examples may fail.
Suggested change
parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT",
parser.add_argument("--mfc", type=json.loads, default={}, metavar="DICT",
Steps of Reproduction ✅
1. Run the example script without supplying --mfc: `python3
examples/2D_burning_droplet/case_hydrocarbon.py`. The parser option is declared at
`examples/2D_burning_droplet/case_hydrocarbon.py:36-37` and the default value for `--mfc`
is the string `"{}"`.

2. Argument parsing happens next at `examples/2D_burning_droplet/case_hydrocarbon.py:45`
where `args = parser.parse_args()` is executed; because argparse applies `type=json.loads`
only to provided CLI values, the default remains the string `"{}"`.

3. The code later constructs the `case` dictionary and uses `args.mfc.get("mpi", True)` at
`examples/2D_burning_droplet/case_hydrocarbon.py:199-200` (the `"parallel_io"` entry). At
that moment `args.mfc` is the string `"{}"`, so `args.mfc.get(...)` raises AttributeError:
"'str' object has no attribute 'get'".

4. Reproduced outcome: running the script with no `--mfc` raises AttributeError during
module execution and the example fails to produce the case dictionary. This is a real
runtime crash on the normal example execution path.
Prompt for AI Agent 🤖
This is a comment left during a code review.

**Path:** examples/2D_burning_droplet/case_hydrocarbon.py
**Line:** 36:36
**Comment:**
	*Type Error: The argparse `--mfc` argument sets `default="{}"` (a string) while `type=json.loads` parses CLI values to a dict; argparse does not apply `type` to the `default` value, so when the option is omitted `args.mfc` will be the string `"{}"` and later `args.mfc.get(...)` will raise an AttributeError. Change the default to an actual dict (e.g. `{}`) so `args.mfc` is a mapping when not provided.

Validate the correctness of the flagged issue. If correct, How can I resolve this? If you propose a fix, implement it and please make it concise.

formatter_class=argparse.ArgumentDefaultsHelpFormatter
)

parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT")
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggestion: The argparse --mfc argument sets default="{}" (a string) while the code later expects args.mfc to be a dict and calls .get(...), which will raise AttributeError if the flag is not provided; set the default to an actual dict so args.mfc.get(...) works when --mfc is omitted. [possible bug]

Severity Level: Critical 🚨
- ❌ Example script crashes at startup when --mfc omitted.
- ❌ Case generation aborted before printing case JSON.
- ⚠️ Output "parallel_io" value never computed.
Suggested change
parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT")
parser.add_argument("--mfc", type=json.loads, default={}, metavar="DICT")
Steps of Reproduction ✅
1. Run the example script without the --mfc flag:

   - Command: python examples/2D_burning_droplet/case_liquid_droplet.py

   - parse_args() invoked at examples/2D_burning_droplet/case_liquid_droplet.py:43.

2. Because add_argument used default="{}" (a string) at line 38, args.mfc becomes the
literal string "{}" (defaults are not converted by the type function).

3. The code later evaluates args.mfc.get(...) when building the case dict at:

   - examples/2D_burning_droplet/case_liquid_droplet.py:209 ("parallel_io": "T" if
   args.mfc.get("mpi", True) else "F").

4. At runtime this raises AttributeError: 'str' object has no attribute 'get', aborting
the script before outputting the case JSON.
Prompt for AI Agent 🤖
This is a comment left during a code review.

**Path:** examples/2D_burning_droplet/case_liquid_droplet.py
**Line:** 38:38
**Comment:**
	*Possible Bug: The argparse `--mfc` argument sets `default="{}"` (a string) while the code later expects `args.mfc` to be a dict and calls `.get(...)`, which will raise AttributeError if the flag is not provided; set the default to an actual dict so `args.mfc.get(...)` works when `--mfc` is omitted.

Validate the correctness of the flagged issue. If correct, How can I resolve this? If you propose a fix, implement it and please make it concise.

# -------------------------------------------------------------------------
# Fluid Properties (ideal gas)
# -------------------------------------------------------------------------
"fluid_pp(1)%gamma": 1.0 / (gamma - 1),
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggestion: Incorrect thermodynamic value: fluid_pp(1)%gamma is set to 1.0 / (gamma - 1) which computes cp/(R) or similar, not the ratio of specific heats; if gamma is the ratio cp/cv (e.g. 1.4) the solver likely expects that value directly. Assign the gamma variable instead of computing 1.0/(gamma-1). [logic error]

Severity Level: Critical 🚨
- ❌ Example inputs incorrect thermodynamic ratio.
- ⚠️ Simulation results from example may be physically invalid.
Suggested change
"fluid_pp(1)%gamma": 1.0 / (gamma - 1),
"fluid_pp(1)%gamma": gamma,
Steps of Reproduction ✅
1. Inspect the constant definition of `gamma` in the file:

   - File: examples/2D_burning_droplet/test_chemistry_only.py, line 46: `gamma = 1.4`.

2. Run the example script to print the case JSON:

   - Command: python3 examples/2D_burning_droplet/test_chemistry_only.py

   - The script prints the `case` dictionary at `if __name__ == "__main__":` (line 203).

3. In the printed JSON, locate the key `"fluid_pp(1)%gamma"` (set at line 162).

   - The produced value is 1.0 / (1.4 - 1) = 2.5 due to the current expression on line
   162.

4. Observe that this value (2.5) does not match the local `gamma` variable (1.4 at line
46) and so the example emits an inconsistent thermodynamic parameter.

Note: Reproduction requires only running this example and inspecting its output; no
additional modules are needed to confirm the incorrect value.
Prompt for AI Agent 🤖
This is a comment left during a code review.

**Path:** examples/2D_burning_droplet/test_chemistry_only.py
**Line:** 162:162
**Comment:**
	*Logic Error: Incorrect thermodynamic value: `fluid_pp(1)%gamma` is set to `1.0 / (gamma - 1)` which computes cp/(R) or similar, not the ratio of specific heats; if `gamma` is the ratio cp/cv (e.g. 1.4) the solver likely expects that value directly. Assign the `gamma` variable instead of computing `1.0/(gamma-1)`.

Validate the correctness of the flagged issue. If correct, How can I resolve this? If you propose a fix, implement it and please make it concise.

formatter_class=argparse.ArgumentDefaultsHelpFormatter
)

parser.add_argument("--no-chemistry", action="store_true", default=True,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggestion: The argparse flag --no-chemistry is declared with action="store_true" but also sets default=True, which means args.no_chemistry will always be True (whether the flag is passed or not) and the chemistry configuration block will never run; remove or set the default to False so the flag behaves as intended. [logic error]

Severity Level: Critical 🚨
- ❌ Chemistry block skipped in example case output.
- ⚠️ Multiphase chemistry parameters never get set.
- ⚠️ Example fails to validate phase-change chemistry logic.
Suggested change
parser.add_argument("--no-chemistry", action="store_true", default=True,
parser.add_argument("--no-chemistry", action="store_true", default=False,
Steps of Reproduction ✅
1. Run the example script: python3 examples/2D_burning_droplet/test_phase_change_only.py
(argument parsing occurs at lines 22-36; parser defined at lines 22-26 and args parsed at
line 36).

2. Inspect parsed value: due to the definition on lines 28-29
(parser.add_argument("--no-chemistry", action="store_true", default=True,...)),
args.no_chemistry is True even when no flag is provided.

3. Execution reaches the chemistry configuration guard at line 267: "if not
args.no_chemistry:" (lines 267-274 contain the chemistry block). Because args.no_chemistry
is always True, the condition is False and the entire chemistry block is skipped.

4. Observe output: printed JSON from lines 298-299 (print(json.dumps(case))) does not
contain any chemistry-related keys (e.g., "chemistry", "cantera_file"), verifying
chemistry was never enabled in the produced case dictionary.
Prompt for AI Agent 🤖
This is a comment left during a code review.

**Path:** examples/2D_burning_droplet/test_phase_change_only.py
**Line:** 28:28
**Comment:**
	*Logic Error: The argparse flag `--no-chemistry` is declared with `action="store_true"` but also sets `default=True`, which means `args.no_chemistry` will always be True (whether the flag is passed or not) and the chemistry configuration block will never run; remove or set the default to False so the flag behaves as intended.

Validate the correctness of the flagged issue. If correct, How can I resolve this? If you propose a fix, implement it and please make it concise.

@codeant-ai
Copy link
Contributor

codeant-ai bot commented Jan 26, 2026

CodeAnt AI finished reviewing your PR.

Copy link
Contributor

@cubic-dev-ai cubic-dev-ai bot left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

7 issues found across 37 files

Prompt for AI agents (all issues)

Check if these issues are valid — if so, understand the root cause of each and fix them.


<file name="examples/2D_burning_droplet/README.md">

<violation number="1" location="examples/2D_burning_droplet/README.md:35">
P2: README contradicts itself about whether phase change and chemistry are coupled, which can mislead users on supported workflows.</violation>
</file>

<file name="examples/2D_burning_droplet/VALIDATION_PHASE1.md">

<violation number="1" location="examples/2D_burning_droplet/VALIDATION_PHASE1.md:127">
P2: Validation doc expects chemistry inactive at threshold, but implementation treats `alpha_gas == gas_phase_threshold` as active. Update the test matrix to match the `<` comparison in code.</violation>
</file>

<file name="src/simulation/m_mpi_proxy.fpp">

<violation number="1" location="src/simulation/m_mpi_proxy.fpp:132">
P2: MPI_BCAST uses MPI_DOUBLE_PRECISION for chem_params%gas_phase_threshold even though the field is real(wp) and mpi_p is the project’s MPI type for wp. This can break single-precision builds or cause MPI datatype mismatches; use mpi_p instead.</violation>
</file>

<file name="examples/2D_burning_droplet/validate_results.py">

<violation number="1" location="examples/2D_burning_droplet/validate_results.py:222">
P2: run_validation always returns success after printing a manual checklist and never executes the declared validation checks, so invalid outputs would still pass when this script is used programmatically.</violation>
</file>

<file name="examples/2D_burning_droplet/viz.py">

<violation number="1" location="examples/2D_burning_droplet/viz.py:42">
P2: Visualization hard-codes nx/ny=200 even though the case supports different resolutions, which can lead to reshape errors or incorrect plots for non-default runs.</violation>

<violation number="2" location="examples/2D_burning_droplet/viz.py:88">
P2: Saving to output_dir assumes the directory already exists; passing a new directory will raise FileNotFoundError. Create the directory before calling savefig.</violation>

<violation number="3" location="examples/2D_burning_droplet/viz.py:171">
P2: Timestep selection uses substring matching on directory paths, which can pick the wrong timestep when multiple directories contain the requested digits (e.g., 1 matches D0010).</violation>
</file>

Since this is your first cubic review, here's how it works:

  • cubic automatically reviews your code and comments on bugs and improvements
  • Teach cubic by replying to its comments. cubic learns from your replies and gets better over time
  • Ask questions if you need clarification on any suggestion

Reply with feedback, questions, or to request a fix. Tag @cubic-dev-ai to re-run a review.

| Model | 5/6-equation + relaxation | 5-equation + species |
| Physics | Evaporation/condensation | Reactions/diffusion |

**The modules are not yet coupled.** This means:
Copy link
Contributor

@cubic-dev-ai cubic-dev-ai bot Jan 26, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

P2: README contradicts itself about whether phase change and chemistry are coupled, which can mislead users on supported workflows.

Prompt for AI agents
Check if this issue is valid — if so, understand the root cause and fix it. At examples/2D_burning_droplet/README.md, line 35:

<comment>README contradicts itself about whether phase change and chemistry are coupled, which can mislead users on supported workflows.</comment>

<file context>
@@ -0,0 +1,278 @@
+| Model | 5/6-equation + relaxation | 5-equation + species |
+| Physics | Evaporation/condensation | Reactions/diffusion |
+
+**The modules are not yet coupled.** This means:
+- Phase change can vaporize a droplet but won't react the vapor
+- Chemistry can burn fuel vapor but needs pre-vaporized initial conditions
</file context>
Fix with Cubic

|-----------|----------------|------------|-------|
| 0.01 | 0.005 | NO | Below threshold |
| 0.01 | 0.015 | YES | Above threshold |
| 0.01 | 0.01 | NO | At threshold (< not <=) |
Copy link
Contributor

@cubic-dev-ai cubic-dev-ai bot Jan 26, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

P2: Validation doc expects chemistry inactive at threshold, but implementation treats alpha_gas == gas_phase_threshold as active. Update the test matrix to match the < comparison in code.

Prompt for AI agents
Check if this issue is valid — if so, understand the root cause and fix it. At examples/2D_burning_droplet/VALIDATION_PHASE1.md, line 127:

<comment>Validation doc expects chemistry inactive at threshold, but implementation treats `alpha_gas == gas_phase_threshold` as active. Update the test matrix to match the `<` comparison in code.</comment>

<file context>
@@ -0,0 +1,257 @@
+|-----------|----------------|------------|-------|
+| 0.01 | 0.005 | NO | Below threshold |
+| 0.01 | 0.015 | YES | Above threshold |
+| 0.01 | 0.01 | NO | At threshold (< not <=) |
+| 0.10 | 0.05 | NO | Higher threshold |
+| 0.10 | 0.15 | YES | Above higher threshold |
</file context>
Fix with Cubic

call MPI_BCAST(chem_params%${VAR}$, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
#:endfor
call MPI_BCAST(chem_params%gas_phase_threshold, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
Copy link
Contributor

@cubic-dev-ai cubic-dev-ai bot Jan 26, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

P2: MPI_BCAST uses MPI_DOUBLE_PRECISION for chem_params%gas_phase_threshold even though the field is real(wp) and mpi_p is the project’s MPI type for wp. This can break single-precision builds or cause MPI datatype mismatches; use mpi_p instead.

Prompt for AI agents
Check if this issue is valid — if so, understand the root cause and fix it. At src/simulation/m_mpi_proxy.fpp, line 132:

<comment>MPI_BCAST uses MPI_DOUBLE_PRECISION for chem_params%gas_phase_threshold even though the field is real(wp) and mpi_p is the project’s MPI type for wp. This can break single-precision builds or cause MPI datatype mismatches; use mpi_p instead.</comment>

<file context>
@@ -121,13 +121,15 @@ contains
                 call MPI_BCAST(chem_params%${VAR}$, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
             #:endfor
+
+            call MPI_BCAST(chem_params%gas_phase_threshold, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
         end if
 
</file context>
Suggested change
call MPI_BCAST(chem_params%gas_phase_threshold, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
call MPI_BCAST(chem_params%gas_phase_threshold, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)
Fix with Cubic

print("Mark each item as [X] when verified.")
print()

return True
Copy link
Contributor

@cubic-dev-ai cubic-dev-ai bot Jan 26, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

P2: run_validation always returns success after printing a manual checklist and never executes the declared validation checks, so invalid outputs would still pass when this script is used programmatically.

Prompt for AI agents
Check if this issue is valid — if so, understand the root cause and fix it. At examples/2D_burning_droplet/validate_results.py, line 222:

<comment>run_validation always returns success after printing a manual checklist and never executes the declared validation checks, so invalid outputs would still pass when this script is used programmatically.</comment>

<file context>
@@ -0,0 +1,227 @@
+    print("Mark each item as [X] when verified.")
+    print()
+    
+    return True
+
+
</file context>
Fix with Cubic


# Read grid dimensions from case (you may need to parse case.py)
# For now, assuming standard dimensions
nx = 200
Copy link
Contributor

@cubic-dev-ai cubic-dev-ai bot Jan 26, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

P2: Visualization hard-codes nx/ny=200 even though the case supports different resolutions, which can lead to reshape errors or incorrect plots for non-default runs.

Prompt for AI agents
Check if this issue is valid — if so, understand the root cause and fix it. At examples/2D_burning_droplet/viz.py, line 42:

<comment>Visualization hard-codes nx/ny=200 even though the case supports different resolutions, which can lead to reshape errors or incorrect plots for non-default runs.</comment>

<file context>
@@ -0,0 +1,185 @@
+    
+    # Read grid dimensions from case (you may need to parse case.py)
+    # For now, assuming standard dimensions
+    nx = 200
+    ny = 200
+    
</file context>
Fix with Cubic

plt.tight_layout()

if output_dir:
output_file = os.path.join(output_dir, 'burning_droplet.png')
Copy link
Contributor

@cubic-dev-ai cubic-dev-ai bot Jan 26, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

P2: Saving to output_dir assumes the directory already exists; passing a new directory will raise FileNotFoundError. Create the directory before calling savefig.

Prompt for AI agents
Check if this issue is valid — if so, understand the root cause and fix it. At examples/2D_burning_droplet/viz.py, line 88:

<comment>Saving to output_dir assumes the directory already exists; passing a new directory will raise FileNotFoundError. Create the directory before calling savefig.</comment>

<file context>
@@ -0,0 +1,185 @@
+    plt.tight_layout()
+    
+    if output_dir:
+        output_file = os.path.join(output_dir, 'burning_droplet.png')
+        plt.savefig(output_file, dpi=200, bbox_inches='tight')
+        print(f"Saved: {output_file}")
</file context>
Fix with Cubic

if args.timestep == -1:
data_dir = data_dirs[-1]
else:
matching = [d for d in data_dirs if str(args.timestep) in d]
Copy link
Contributor

@cubic-dev-ai cubic-dev-ai bot Jan 26, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

P2: Timestep selection uses substring matching on directory paths, which can pick the wrong timestep when multiple directories contain the requested digits (e.g., 1 matches D0010).

Prompt for AI agents
Check if this issue is valid — if so, understand the root cause and fix it. At examples/2D_burning_droplet/viz.py, line 171:

<comment>Timestep selection uses substring matching on directory paths, which can pick the wrong timestep when multiple directories contain the requested digits (e.g., 1 matches D0010).</comment>

<file context>
@@ -0,0 +1,185 @@
+    if args.timestep == -1:
+        data_dir = data_dirs[-1]
+    else:
+        matching = [d for d in data_dirs if str(args.timestep) in d]
+        if matching:
+            data_dir = matching[0]
</file context>
Fix with Cubic

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

Review effort 4/5 size:XXL This PR changes 1000+ lines, ignoring generated files

Development

Successfully merging this pull request may close these issues.

2 participants