Multifluid acoustic source adds the full mass source to every partial density

[adob]

Summary

In src/simulation/m_acoustic_src.fpp, the acoustic-source routine computes one scalar mass source, mass_src, for each affected cell. In multifluid 5-equation and 6-equation cases, that same full scalar source is then added to every partial-density continuity equation.

For a cell with partial densities (m_k = \alpha_k \rho_k), the current update is effectively:

$$ \dot m_k \leftarrow \dot m_k + S $$

for every fluid (k), where (S = \texttt{mass_src}).

For (N) fluids this makes the total mixture-density source

$$ \sum_k \dot m_k = N S $$

rather than (S), and it changes the local composition unless all partial densities are equal.

For two fluids, for example,

$$ m_1' = m_1 + S,\Delta t, \qquad m_2' = m_2 + S,\Delta t, $$

so the total density increment is $2S\Delta t$. The mass fractions also change spuriously because the same absolute increment is applied to both phases.

A composition-preserving update is to distribute the mixture mass source according to the local partial-density fractions:

$$ \dot m_k \leftarrow \dot m_k + S\frac{m_k}{\sum_j m_j}. $$

Then the partial-density increments sum to the intended mixture source (S) and the source does not create or destroy composition.

Affected revisions checked

This behavior is present in:

• MFC v5.7.0, commit 0b3d56c8e1c9a3aa8d112681b72c87ab69fe5f13
• Current upstream master checked at commit fef3c276a4e7deaa39d8736deaa36c9ce9ad5bf5

Current code

On upstream master at fef3c276 in src/simulation/m_acoustic_src.fpp, the relevant code is:

! Update the rhs variables
$:GPU_PARALLEL_LOOP(private='[j, k, l]',collapse=3)
do l = 0, p
    do k = 0, n
        do j = 0, m
            $:GPU_LOOP(parallelism='[seq]')
            do q = eqn_idx%cont%beg, eqn_idx%cont%end
                rhs_vf(q)%sf(j, k, l) = rhs_vf(q)%sf(j, k, l) + mass_src(j, k, l)
            end do
            $:GPU_LOOP(parallelism='[seq]')
            do q = eqn_idx%mom%beg, eqn_idx%mom%end
                rhs_vf(q)%sf(j, k, l) = rhs_vf(q)%sf(j, k, l) + &
                    mom_src(q - eqn_idx%cont%end, j, k, l)
            end do
            rhs_vf(eqn_idx%E)%sf(j, k, l) = &
                rhs_vf(eqn_idx%E)%sf(j, k, l) + e_src(j, k, l)
        end do
    end do
end do
$:END_GPU_PARALLEL_LOOP()

The same routine computes mass_src once per cell:

mass_src(j, k, l) = mass_src(j, k, l) + mass_src_diff

and computes the acoustic energy source once from the same mass_src_diff:

E_src(j, k, l) = E_src(j, k, l) + &
    mass_src_diff*c**2._wp*gamma_mix

This is consistent with mass_src being a mixture-level acoustic source, not a separate per-fluid source.

Why this is problematic

For a multifluid conservative state,

$$ \rho = \sum_k \alpha_k \rho_k. $$

The acoustic source computes a single density/mass source for the mixture. Adding that same value to every alpha_rho_k equation multiplies the total density forcing by the number of fluids.

It also changes composition. With two unequal partial densities (m_1 \ne m_2), the current update gives

$$ Y_1' = \frac{m_1 + \delta} {m_1 + m_2 + 2\delta}, $$

which is generally not equal to the original

$$ Y_1 = \frac{m_1}{m_1+m_2}. $$

Thus an acoustic wave source can create minority-phase mass even when the intended acoustic forcing should preserve the local mixture composition.

This is especially important for diffuse-interface gas/liquid cases where one phase may be present only at a trace volume fraction. The source can artificially increase the trace phase and alter phasic thermodynamic quantities.

Expected behavior

For 5eq/6eq cases with more than one fluid, the scalar mixture mass source should be partitioned among partial densities so that

$$ \sum_k \Delta(\alpha_k\rho_k) = S,\Delta t $$

and, absent an explicitly phase-selective source, the local mass fractions should remain unchanged.

One natural implementation is

$$ \Delta(\alpha_k\rho_k) = S,\Delta t \frac{\alpha_k\rho_k}{\rho}. $$

Single-fluid behavior should remain unchanged.

Proposed fix

A local patch partitions the mixture mass source according to the current partial-density fractions:

! Update the rhs variables
$:GPU_PARALLEL_LOOP(private='[j, k, l, q, myRho]',collapse=3)
do l = 0, p
    do k = 0, n
        do j = 0, m
            if ((model_eqns == model_eqns_5eq .or. &
                 model_eqns == model_eqns_6eq) .and. &
                num_fluids > 1) then

                ! mass_src is a mixture-density source. Distribute it
                ! across partial densities by the local mass fractions
                ! so the source does not create or destroy composition.
                myRho = 0._wp
                $:GPU_LOOP(parallelism='[seq]')
                do q = eqn_idx%cont%beg, eqn_idx%cont%end
                    myRho = myRho + q_cons_vf(q)%sf(j, k, l)
                end do

                if (myRho > 0._wp) then
                    $:GPU_LOOP(parallelism='[seq]')
                    do q = eqn_idx%cont%beg, eqn_idx%cont%end
                        rhs_vf(q)%sf(j, k, l) = &
                            rhs_vf(q)%sf(j, k, l) + &
                            mass_src(j, k, l) * &
                            q_cons_vf(q)%sf(j, k, l) / myRho
                    end do
                end if
            else
                $:GPU_LOOP(parallelism='[seq]')
                do q = eqn_idx%cont%beg, eqn_idx%cont%end
                    rhs_vf(q)%sf(j, k, l) = &
                        rhs_vf(q)%sf(j, k, l) + mass_src(j, k, l)
                end do
            end if

            $:GPU_LOOP(parallelism='[seq]')
            do q = eqn_idx%mom%beg, eqn_idx%mom%end
                rhs_vf(q)%sf(j, k, l) = &
                    rhs_vf(q)%sf(j, k, l) + &
                    mom_src(q - eqn_idx%cont%end, j, k, l)
            end do

            rhs_vf(eqn_idx%E)%sf(j, k, l) = &
                rhs_vf(eqn_idx%E)%sf(j, k, l) + e_src(j, k, l)
        end do
    end do
end do
$:END_GPU_PARALLEL_LOOP()

The patch applies cleanly as one commit on top of current upstream master:

• Upstream base: fef3c276
• Local patch: f7f7aac0 — Preserve composition in multiphase acoustic source

Minimal regression test

A small regression can isolate the issue without bubble dynamics or pressure relaxation.

Initial condition

Use either 5eq or 6eq with two fluids in a uniform cell/domain:

• uniform pressure
• uniform velocity
• unequal partial densities, for example (m_1 \gg m_2)
• both fluids present
• one acoustic source with nonzero mass forcing

Choose a source/time at which mass_src != 0.

Quantities to check

For one RHS evaluation or sufficiently short timestep, let

$$ \Delta m_k = m_k(t+\Delta t)-m_k(t). $$

The regression should verify:

1. Total source conservation

$$ \sum_k \Delta m_k \approx \Delta t,S. $$

2. Composition preservation

$$ \frac{m_k(t+\Delta t)} {\sum_j m_j(t+\Delta t)} \approx \frac{m_k(t)} {\sum_j m_j(t)}. $$

With the current upstream implementation, the first condition becomes approximately (N\Delta t S), and the second condition fails for unequal partial densities.

With the proposed partitioning, both conditions should be satisfied up to time-integration/numerical error.

It would be useful to run this regression for both model_eqns = 5eq and model_eqns = 6eq, plus a single-fluid case to confirm unchanged behavior.

Scope

This report concerns only the acoustic-source distribution into the multifluid continuity equations.

A separate 6eq moving heterogeneous-EOS contact problem observed in MFC v5.7.0 does not reproduce on current upstream master and is not part of this report.

Environment used during investigation

The issue was identified while using MFC for a water/argon weak-acoustic multiphase case.

Relevant local build environment:

• NVIDIA RTX 4090
• NVHPC 25.9
• CUDA 12.9
• OpenACC GPU build
• double precision

The core issue described above is visible directly in the RHS source update and is not expected to depend on the GPU/compiler configuration.

Suggested resolution

1. Partition mass_src among multifluid partial-density equations according to local mass fractions, or otherwise document and implement the intended phase-specific source semantics.
2. Add a regression checking both total mixture mass injection and composition preservation.
3. Retain the existing behavior for single-fluid cases.