Skip to content

Conversation

DimAdam-01
Copy link
Contributor

@DimAdam-01 DimAdam-01 commented Jul 11, 2025

User description

Description

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

This update implements the diffusive fluxes for chemical species, includes thermal conduction in the energy equation, and computes the mixture’s viscosity.

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: This test examines one-dimensional multicomponent diffusion of four distinct species, as detailed in DOI 10.2514/6.2020-1751.
image image
  • Test B: Test B examines a two-dimensional premixed flame with an imposed perturbation, then tracks the ensuing thermo-diffusive instability growth.
image image

Test Configuration:

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

My Laptop, Delta (A100) and Delta AI (GH)

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


Description

  • Implement multicomponent diffusion fluxes for chemical species

  • Add thermal conduction to energy equation

  • Integrate mixture viscosity computation for chemistry

  • Add new test case for multicomponent diffusion validation


Diagram Walkthrough

flowchart LR
  A["Chemistry Module"] --> B["Diffusion Flux Computation"]
  A --> C["Viscosity Calculation"]
  B --> D["Species Mass Fluxes"]
  B --> E["Thermal Conduction"]
  D --> F["Energy Equation Update"]
  E --> F
  C --> G["Riemann Solver Integration"]
  F --> H["RHS Updates"]
  G --> H
Loading

File Walkthrough

Relevant files
Enhancement
3 files
m_rhs.fpp
Add diffusion flux computation and RHS updates                     
+131/-37
m_chemistry.fpp
Implement diffusion flux and viscosity calculations           
+166/-1 
m_riemann_solvers.fpp
Integrate viscosity computation in Riemann solvers             
+51/-0   
Tests
10 files
1dHardcodedIC.fpp
Add hardcoded initial conditions for diffusion test           
+30/-0   
cases.py
Add multicomponent diffusion test case                                     
+16/-2   
golden-metadata.txt
Update test metadata for new environment                                 
+10/-144
golden-metadata.txt
Update test metadata for new environment                                 
+10/-144
golden-metadata.txt
Update test metadata for new environment                                 
+10/-144
golden-metadata.txt
Update test metadata for new environment                                 
+10/-144
golden-metadata.txt
Update test metadata for new environment                                 
+10/-144
golden-metadata.txt
Update test metadata for new environment                                 
+10/-144
golden-metadata.txt
Update test metadata for new environment                                 
+10/-144
golden-metadata.txt
Add new test metadata file                                                             
+18/-0   
Bug fix
1 files
m_variables_conversion.fpp
Fix tolerance variable reference in chemistry                       
+1/-2     
Configuration changes
4 files
m_start_up.fpp
Add GPU memory updates for chemistry parameters                   
+3/-0     
m_global_parameters.fpp
Update GPU memory for chemistry parameters                             
+2/-0     
case.py
Update wave speeds and average state parameters                   
+2/-2     
input.py
Add OpenACC directive to Pyrometheus code generation         
+1/-1     
Dependencies
1 files
pyproject.toml
Update Pyrometheus dependency version                                       
+2/-2     
Additional files
5 files
golden.txt +24/-24 
golden.txt +56/-0   
golden.txt +24/-24 
golden.txt +24/-24 
golden.txt +24/-24 

@sbryngelson
Copy link
Member

@DimAdam-01 can you fix your merge conflicts

@sbryngelson
Copy link
Member

status? also please use descriptive commit messages. also, spamming commits causes the CI to run dozens of times and clogs up the self hosted runners.

@DimAdam-01 DimAdam-01 marked this pull request as draft September 2, 2025 17:32
@sbryngelson
Copy link
Member

this seems to be looking pretty good. is there anything missing at the moment? @DimAdam-01

@DimAdam-01
Copy link
Contributor Author

DimAdam-01 commented Sep 24, 2025

this seems to be looking pretty good. is there anything missing at the moment? @DimAdam-01

No everything is there, just to mention, in order to get the code to run on the Frontier , I had to modify the Pyrometheus code generation script (in toolchain/mfc/run/input.py). I added the directive_offload="acc", which was necessary for the code to compile and pass all tests. @sbryngelson

@sbryngelson
Copy link
Member

this seems to be looking pretty good. is there anything missing at the moment? @DimAdam-01

No everything is there, just to mention, in order to get the code to run on the Frontier , I had to modify the Pyrometheus code generation script (in toolchain/mfc/run/input.py). I added the directive_offload="acc", which was necessary for the code to compile and pass all tests. @sbryngelson

Got it. That's fine. @prathi-wind can fix that, oacc is default for now anyway. I saw his latest push to his draft PR and the only failing OMP test was chemistry, so somethings there need to be fixed anyway.

@DimAdam-01 DimAdam-01 marked this pull request as ready for review September 25, 2025 00:56
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 diffusion flux energy correction accumulates within the species loop, subtracting h_kYsrho_Vic repeatedly; this looks like it should be applied once after computing rho_Vic and the raw species fluxes, not accumulated per species. Please verify conservation and that Mass_Diffu_Energy is not over-subtracted.

$:GPU_LOOP(parallelism='[seq]')
do eqn = chemxb, chemxe
    Mass_Diffu_Flux(eqn - chemxb + 1) = rho_cell*mass_diffusivities_mixavg_Cell(eqn - chemxb + 1)* &
                                        molecular_weights(eqn - chemxb + 1)/MW_cell*dXk_dxi(eqn - chemxb + 1)
    rho_Vic = rho_Vic + Mass_Diffu_Flux(eqn - chemxb + 1)
    Mass_Diffu_Energy = Mass_Diffu_Energy + h_k(eqn - chemxb + 1)*Mass_Diffu_Flux(eqn - chemxb + 1)
end do

! Apply corrections for mass conservation
$:GPU_LOOP(parallelism='[seq]')
do eqn = chemxb, chemxe
    Mass_Diffu_Energy = Mass_Diffu_Energy - h_k(eqn - chemxb + 1)*Ys_cell(eqn - chemxb + 1)*rho_Vic
    Mass_Diffu_Flux(eqn - chemxb + 1) = Mass_Diffu_Flux(eqn - chemxb + 1) - rho_Vic*Ys_cell(eqn - chemxb + 1)
end do

! Add thermal conduction contribution
Mass_Diffu_Energy = lambda_Cell*dT_dxi + Mass_Diffu_Energy

! Update flux arrays
Duplicate Logic

The per-direction RHS updates for viscous/surface_tension vs diffusion duplicate the same divergence pattern three times with only dx/dy/dz and index strides differing; consider factoring into a small helper to reduce maintenance risk and keep diffusion and viscous paths consistent.

    if (surface_tension .or. viscous) then
        $:GPU_PARALLEL_LOOP(collapse=3)
        do l = 0, p
            do k = 0, n
                do j = 0, m
                    $:GPU_LOOP(parallelism='[seq]')
                    do i = momxb, E_idx
                        rhs_vf(i)%sf(j, k, l) = &
                            rhs_vf(i)%sf(j, k, l) + 1._wp/dx(j)* &
                            (flux_src_n_in(i)%sf(j - 1, k, l) &
                             - flux_src_n_in(i)%sf(j, k, l))
                    end do
                end do
            end do
        end do
    end if

    if (chem_params%diffusion) then
        $:GPU_PARALLEL_LOOP(collapse=3)
        do l = 0, p
            do k = 0, n
                do j = 0, m
                    $:GPU_LOOP(parallelism='[seq]')
                    do i = chemxb, chemxe
                        rhs_vf(i)%sf(j, k, l) = &
                            rhs_vf(i)%sf(j, k, l) + 1._wp/dx(j)* &
                            (flux_src_n_in(i)%sf(j - 1, k, l) &
                             - flux_src_n_in(i)%sf(j, k, l))
                    end do

                    if (.not. viscous) then
                        rhs_vf(E_idx)%sf(j, k, l) = &
                            rhs_vf(E_idx)%sf(j, k, l) + 1._wp/dx(j)* &
                            (flux_src_n_in(E_idx)%sf(j - 1, k, l) &
                             - flux_src_n_in(E_idx)%sf(j, k, l))
                    end if
                end do
            end do
        end do
    end if

elseif (idir == 2) then ! y-direction

    if (surface_tension) then
        $:GPU_PARALLEL_LOOP(collapse=3)
        do l = 0, p
            do k = 0, n
                do j = 0, m
                    rhs_vf(c_idx)%sf(j, k, l) = &
                        rhs_vf(c_idx)%sf(j, k, l) + 1._wp/dy(k)* &
                        q_prim_vf(c_idx)%sf(j, k, l)* &
                        (flux_src_n_in(advxb)%sf(j, k, l) - &
                         flux_src_n_in(advxb)%sf(j, k - 1, l))
                end do
            end do
        end do
    end if

    if (cyl_coord .and. ((bc_y%beg == -2) .or. (bc_y%beg == -14))) then
        if (viscous) then
            if (p > 0) then
                call s_compute_viscous_stress_tensor(q_prim_vf, &
                                                     dq_prim_dx_vf(mom_idx%beg:mom_idx%end), &
                                                     dq_prim_dy_vf(mom_idx%beg:mom_idx%end), &
                                                     dq_prim_dz_vf(mom_idx%beg:mom_idx%end), &
                                                     tau_Re_vf, &
                                                     idwbuff(1), idwbuff(2), idwbuff(3))
            else
                call s_compute_viscous_stress_tensor(q_prim_vf, &
                                                     dq_prim_dx_vf(mom_idx%beg:mom_idx%end), &
                                                     dq_prim_dy_vf(mom_idx%beg:mom_idx%end), &
                                                     dq_prim_dy_vf(mom_idx%beg:mom_idx%end), &
                                                     tau_Re_vf, &
                                                     idwbuff(1), idwbuff(2), idwbuff(3))
            end if

            $:GPU_PARALLEL_LOOP(collapse=2)
            do l = 0, p
                do j = 0, m
                    $:GPU_LOOP(parallelism='[seq]')
                    do i = momxb, E_idx
                        rhs_vf(i)%sf(j, 0, l) = &
                            rhs_vf(i)%sf(j, 0, l) + 1._wp/(y_cc(1) - y_cc(-1))* &
                            (tau_Re_vf(i)%sf(j, -1, l) &
                             - tau_Re_vf(i)%sf(j, 1, l))
                    end do
                end do
            end do

        end if

        $:GPU_PARALLEL_LOOP(collapse=3)
        do l = 0, p
            do k = 1, n
                do j = 0, m
                    $:GPU_LOOP(parallelism='[seq]')
                    do i = momxb, E_idx
                        rhs_vf(i)%sf(j, k, l) = &
                            rhs_vf(i)%sf(j, k, l) + 1._wp/dy(k)* &
                            (flux_src_n_in(i)%sf(j, k - 1, l) &
                             - flux_src_n_in(i)%sf(j, k, l))
                    end do
                end do
            end do
        end do

    else

        if (viscous .or. surface_tension) then
            $:GPU_PARALLEL_LOOP(collapse=3)
            do l = 0, p
                do k = 0, n
                    do j = 0, m
                        $:GPU_LOOP(parallelism='[seq]')
                        do i = momxb, E_idx
                            rhs_vf(i)%sf(j, k, l) = &
                                rhs_vf(i)%sf(j, k, l) + 1._wp/dy(k)* &
                                (flux_src_n_in(i)%sf(j, k - 1, l) &
                                 - flux_src_n_in(i)%sf(j, k, l))
                        end do
                    end do
                end do
            end do
        end if

        if (chem_params%diffusion) then
            $:GPU_PARALLEL_LOOP(collapse=3)
            do l = 0, p
                do k = 0, n
                    do j = 0, m
                        $:GPU_LOOP(parallelism='[seq]')
                        do i = chemxb, chemxe
                            rhs_vf(i)%sf(j, k, l) = &
                                rhs_vf(i)%sf(j, k, l) + 1._wp/dy(k)* &
                                (flux_src_n_in(i)%sf(j, k - 1, l) &
                                 - flux_src_n_in(i)%sf(j, k, l))
                        end do
                        if (.not. viscous) then
                            rhs_vf(E_idx)%sf(j, k, l) = &
                                rhs_vf(E_idx)%sf(j, k, l) + 1._wp/dy(k)* &
                                (flux_src_n_in(E_idx)%sf(j, k - 1, l) &
                                 - flux_src_n_in(E_idx)%sf(j, k, l))
                        end if
                    end do
                end do
            end do
        end if
    end if

    ! Applying the geometrical viscous Riemann source fluxes calculated as average
    ! of values at cell boundaries
    if (cyl_coord) then
        if ((bc_y%beg == -2) .or. (bc_y%beg == -14)) then

            $:GPU_PARALLEL_LOOP(collapse=3)
            do l = 0, p
                do k = 1, n
                    do j = 0, m
                        $:GPU_LOOP(parallelism='[seq]')
                        do i = momxb, E_idx
                            rhs_vf(i)%sf(j, k, l) = &
                                rhs_vf(i)%sf(j, k, l) - 5.e-1_wp/y_cc(k)* &
                                (flux_src_n_in(i)%sf(j, k - 1, l) &
                                 + flux_src_n_in(i)%sf(j, k, l))
                        end do
                    end do
                end do
            end do

            if (viscous) then
                $:GPU_PARALLEL_LOOP(collapse=2)
                do l = 0, p
                    do j = 0, m
                        $:GPU_LOOP(parallelism='[seq]')
                        do i = momxb, E_idx
                            rhs_vf(i)%sf(j, 0, l) = &
                                rhs_vf(i)%sf(j, 0, l) - 1._wp/y_cc(0)* &
                                tau_Re_vf(i)%sf(j, 0, l)
                        end do
                    end do
                end do
            end if
        else

            $:GPU_PARALLEL_LOOP(collapse=3)
            do l = 0, p
                do k = 0, n
                    do j = 0, m
                        $:GPU_LOOP(parallelism='[seq]')
                        do i = momxb, E_idx
                            rhs_vf(i)%sf(j, k, l) = &
                                rhs_vf(i)%sf(j, k, l) - 5.e-1_wp/y_cc(k)* &
                                (flux_src_n_in(i)%sf(j, k - 1, l) &
                                 + flux_src_n_in(i)%sf(j, k, l))
                        end do
                    end do
                end do
            end do

        end if
    end if

elseif (idir == 3) then ! z-direction

    if (surface_tension) then
        $:GPU_PARALLEL_LOOP(collapse=3)
        do l = 0, p
            do k = 0, n
                do j = 0, m
                    rhs_vf(c_idx)%sf(j, k, l) = &
                        rhs_vf(c_idx)%sf(j, k, l) + 1._wp/dz(l)* &
                        q_prim_vf(c_idx)%sf(j, k, l)* &
                        (flux_src_n_in(advxb)%sf(j, k, l) - &
                         flux_src_n_in(advxb)%sf(j, k, l - 1))
                end do
            end do
        end do
    end if

    if (viscous .or. surface_tension) then
        $:GPU_PARALLEL_LOOP(collapse=3)
        do l = 0, p
            do k = 0, n
                do j = 0, m
                    $:GPU_LOOP(parallelism='[seq]')
                    do i = momxb, E_idx
                        rhs_vf(i)%sf(j, k, l) = &
                            rhs_vf(i)%sf(j, k, l) + 1._wp/dz(l)* &
                            (flux_src_n_in(i)%sf(j, k, l - 1) &
                             - flux_src_n_in(i)%sf(j, k, l))
                    end do
                end do
            end do
        end do
    end if

    if (chem_params%diffusion) then
        $:GPU_PARALLEL_LOOP(collapse=3)
        do l = 0, p
            do k = 0, n
                do j = 0, m
                    $:GPU_LOOP(parallelism='[seq]')
                    do i = chemxb, chemxe
                        rhs_vf(i)%sf(j, k, l) = &
                            rhs_vf(i)%sf(j, k, l) + 1._wp/dz(l)* &
                            (flux_src_n_in(i)%sf(j, k, l - 1) &
                             - flux_src_n_in(i)%sf(j, k, l))
                    end do
                    if (.not. viscous) then
                        rhs_vf(E_idx)%sf(j, k, l) = &
                            rhs_vf(E_idx)%sf(j, k, l) + 1._wp/dz(l)* &
                            (flux_src_n_in(E_idx)%sf(j, k, l - 1) &
                             - flux_src_n_in(E_idx)%sf(j, k, l))
                    end if
                end do
            end do
        end do
    end if
Likely Typo

The assignment to q_prim_vf at advxb is performed twice consecutively with the same value, which appears redundant and could mask a missing intended variable initialization.

q_prim_vf(advxb)%sf(i, 0, 0) = 1.0_wp
q_prim_vf(advxb)%sf(i, 0, 0) = 1.0_wp

@sbryngelson sbryngelson requested a review from Copilot September 25, 2025 18:03
Copy link
Contributor

@Copilot Copilot AI left a comment

Choose a reason for hiding this comment

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

Pull Request Overview

This PR implements multicomponent diffusion fluxes, thermal conduction, and mixture viscosity for chemical species in the MFC solver. It adds new chemistry functionality to enable more sophisticated multispecies flow simulations with transport phenomena.

Key changes include:

  • Implementation of species mass diffusion flux computations using mixture-averaged approaches
  • Addition of thermal conduction in the energy equation for chemistry cases
  • Integration of mixture viscosity calculations for chemical species

Reviewed Changes

Copilot reviewed 20 out of 24 changed files in this pull request and generated 4 comments.

Show a summary per file
File Description
src/common/m_chemistry.fpp Implements core diffusion flux computation and viscosity calculation routines
src/simulation/m_rhs.fpp Adds diffusion flux RHS updates and memory allocation for chemistry diffusion
src/simulation/m_riemann_solvers.fpp Integrates viscosity computation and adds flux source initialization for chemistry
src/pre_process/include/1dHardcodedIC.fpp Adds hardcoded initial conditions for multicomponent diffusion test case
toolchain/mfc/test/cases.py Adds new test case for multicomponent diffusion validation
src/simulation/m_start_up.fpp Updates GPU memory management for chemistry parameters
src/simulation/m_global_parameters.fpp Updates GPU memory handling for chemistry parameters
src/common/m_variables_conversion.fpp Fixes tolerance variable reference
toolchain/pyproject.toml Updates Pyrometheus dependency version
toolchain/mfc/run/input.py Adds OpenACC directive to Pyrometheus code generation
Comments suppressed due to low confidence (1)

src/simulation/m_rhs.fpp:1

  • Code duplication detected. The RHS flux update pattern is repeated for viscous/surface_tension and diffusion cases. Consider extracting this into a helper subroutine to follow DRY principles.
!>

@sbryngelson sbryngelson merged commit eb152c5 into MFlowCode:master Sep 28, 2025
38 of 39 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Development

Successfully merging this pull request may close these issues.

2 participants