Skip to content

Fix chemistry‑related bugs #829

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 8 commits into
base: master
Choose a base branch
from
Open

Fix chemistry‑related bugs #829

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
bcaaa35
Chemistry Bug Fixes
DimAdam-01 May 6, 2025
4833d7b
Chemistry Bug Fixes
DimAdam-01 May 6, 2025
c8e7f9f
Format
DimAdam-01 May 6, 2025
652d0fa
Testing-Changes-In-CodeGenerator-For-Pyrometheus
DimAdam-01 May 8, 2025
e63683a
Merge branch 'MFlowCode:master' into my_new_branch
DimAdam-01 May 8, 2025
cc20ee9
Extend Frontier Run Time
DimAdam-01 May 8, 2025
86a5b77
Discard changes to .github/workflows/frontier/submit.sh
sbryngelson May 9, 2025
209dd94
Merge branch 'master' into my_new_branch
sbryngelson May 10, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions src/common/m_chemistry.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@
module m_chemistry

use m_thermochem, only: &
num_species, mol_weights, get_temperature, get_net_production_rates
num_species, molecular_weights, get_temperature, get_net_production_rates

use m_global_parameters

Expand Down Expand Up @@ -91,7 +91,7 @@
!$acc loop seq
do eqn = chemxb, chemxe

omega_m = mol_weights(eqn - chemxb + 1)*omega(eqn - chemxb + 1)
omega_m = molecular_weights(eqn - chemxb + 1)*omega(eqn - chemxb + 1)

Check warning on line 94 in src/common/m_chemistry.fpp

View check run for this annotation

Codecov / codecov/patch

src/common/m_chemistry.fpp#L94 

Added line #L94 was not covered by tests

rhs_vf(eqn)%sf(x, y, z) = rhs_vf(eqn)%sf(x, y, z) + omega_m

Expand Down
2 changes: 1 addition & 1 deletion src/post_process/m_mpi_proxy.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -165,7 +165,7 @@ contains
& 'alt_soundspeed', 'hypoelasticity', 'mhd', 'parallel_io', &
& 'rho_wrt', 'E_wrt', 'pres_wrt', 'gamma_wrt', 'sim_data', &
& 'heat_ratio_wrt', 'pi_inf_wrt', 'pres_inf_wrt', 'cons_vars_wrt', &
& 'prim_vars_wrt', 'c_wrt', 'qm_wrt','schlieren_wrt', &
& 'prim_vars_wrt', 'c_wrt', 'qm_wrt','schlieren_wrt','chem_wrt_T', &
& 'bubbles_euler', 'qbmm', 'polytropic', 'polydisperse', &
& 'file_per_process', 'relax', 'cf_wrt', &
& 'adv_n', 'ib', 'cfl_adap_dt', 'cfl_const_dt', 'cfl_dt', &
Expand Down
3 changes: 3 additions & 0 deletions src/pre_process/m_mpi_proxy.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -128,6 +128,9 @@ contains

call MPI_BCAST(patch_icpp(i)%model_spc, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)

if (chemistry) then
call MPI_BCAST(patch_icpp(i)%Y, size(patch_icpp(i)%Y), mpi_p, 0, MPI_COMM_WORLD, ierr)
end if
! Broadcast IB variables
call MPI_BCAST(patch_ib(i)%geometry, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
call MPI_BCAST(patch_ib(i)%model_filepath, len(patch_ib(i)%model_filepath), MPI_CHARACTER, 0, MPI_COMM_WORLD, ierr)
Expand Down
12 changes: 6 additions & 6 deletions src/simulation/include/inline_riemann.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,8 +37,8 @@
call get_species_enthalpies_rt(T_L, h_iL)
call get_species_enthalpies_rt(T_R, h_iR)

h_iL = h_iL*gas_constant/mol_weights*T_L
h_iR = h_iR*gas_constant/mol_weights*T_R
h_iL = h_iL*gas_constant/molecular_weights*T_L
h_iR = h_iR*gas_constant/molecular_weights*T_R
call get_species_specific_heats_r(T_L, Cp_iL)
call get_species_specific_heats_r(T_R, Cp_iR)

Expand All @@ -48,17 +48,17 @@

if (abs(T_L - T_R) < eps) then
! Case when T_L and T_R are very close
Cp_avg = sum(Yi_avg(:)*(0.5_wp*Cp_iL(i) + 0.5_wp*Cp_iR(:))*gas_constant/mol_weights(:))
Cv_avg = sum(Yi_avg(:)*((0.5_wp*Cp_iL(i) + 0.5_wp*Cp_iR(:))*gas_constant/mol_weights(:) - gas_constant/mol_weights(:)))
Cp_avg = sum(Yi_avg(:)*(0.5_wp*Cp_iL(:) + 0.5_wp*Cp_iR(:))*gas_constant/molecular_weights(:))
Cv_avg = sum(Yi_avg(:)*((0.5_wp*Cp_iL(:) + 0.5_wp*Cp_iR(:))*gas_constant/molecular_weights(:) - gas_constant/molecular_weights(:)))
else
! Normal calculation when T_L and T_R are sufficiently different
Cp_avg = sum(Yi_avg(:)*(h_iR(:) - h_iL(:))/(T_R - T_L))
Cv_avg = sum(Yi_avg(:)*((h_iR(:) - h_iL(:))/(T_R - T_L) - gas_constant/mol_weights(:)))
Cv_avg = sum(Yi_avg(:)*((h_iR(:) - h_iL(:))/(T_R - T_L) - gas_constant/molecular_weights(:)))
end if

gamma_avg = Cp_avg/Cv_avg

Phi_avg(:) = (gamma_avg - 1._wp)*(vel_avg_rms/2.0_wp - h_avg_2(:)) + gamma_avg*gas_constant/mol_weights(:)*T_avg
Phi_avg(:) = (gamma_avg - 1._wp)*(vel_avg_rms/2.0_wp - h_avg_2(:)) + gamma_avg*gas_constant/molecular_weights(:)*T_avg
c_sum_Yi_Phi = sum(Yi_avg(:)*Phi_avg(:))
end if

Expand Down
12 changes: 6 additions & 6 deletions src/simulation/m_riemann_solvers.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -510,8 +510,8 @@ contains
call get_mixture_molecular_weight(Ys_L, MW_L)
call get_mixture_molecular_weight(Ys_R, MW_R)

Xs_L(:) = Ys_L(:)*MW_L/mol_weights(:)
Xs_R(:) = Ys_R(:)*MW_R/mol_weights(:)
Xs_L(:) = Ys_L(:)*MW_L/molecular_weights(:)
Xs_R(:) = Ys_R(:)*MW_R/molecular_weights(:)

R_gas_L = gas_constant/MW_L
R_gas_R = gas_constant/MW_R
Expand Down Expand Up @@ -1015,8 +1015,8 @@ contains
Y_L = qL_prim_rs${XYZ}$_vf(j, k, l, i)
Y_R = qR_prim_rs${XYZ}$_vf(j + 1, k, l, i)

flux_rs${XYZ}$_vf(j, k, l, i) = (s_M*Y_R*rho_R*vel_R(dir_idx(norm_dir)) &
- s_P*Y_L*rho_L*vel_L(dir_idx(norm_dir)) &
flux_rs${XYZ}$_vf(j, k, l, i) = (s_M*Y_R*rho_R*vel_R(dir_idx(1)) &
- s_P*Y_L*rho_L*vel_L(dir_idx(1)) &
+ s_M*s_P*(Y_L*rho_L - Y_R*rho_R)) &
/(s_M - s_P)
flux_src_rs${XYZ}$_vf(j, k, l, i) = 0._wp
Expand Down Expand Up @@ -2540,8 +2540,8 @@ contains
call get_mixture_molecular_weight(Ys_L, MW_L)
call get_mixture_molecular_weight(Ys_R, MW_R)

Xs_L(:) = Ys_L(:)*MW_L/mol_weights(:)
Xs_R(:) = Ys_R(:)*MW_R/mol_weights(:)
Xs_L(:) = Ys_L(:)*MW_L/molecular_weights(:)
Xs_R(:) = Ys_R(:)*MW_R/molecular_weights(:)

R_gas_L = gas_constant/MW_L
R_gas_R = gas_constant/MW_R
Expand Down
6 changes: 3 additions & 3 deletions toolchain/mfc/run/input.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -76,10 +76,10 @@ def generate_fpp(self, target) -> None:
# Write the generated Fortran code to the m_thermochem.f90 file with the chosen precision
common.file_write(
os.path.join(modules_dir, "m_thermochem.f90"),
pyro.codegen.fortran90.gen_thermochem_code(
pyro.FortranCodeGenerator().generate(
"m_thermochem",
self.get_cantera_solution(),
module_name="m_thermochem",
real_type=real_type
pyro.CodeGenerationOptions(scalar_type = real_type)
),
True
)
Expand Down
4 changes: 2 additions & 2 deletions toolchain/pyproject.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,8 +38,8 @@ dependencies = [
"matplotlib",

# Chemistry
"cantera==3.0.1",
"pyrometheus==1.0.2"
"cantera==3.1.0",
"pyrometheus@git+https://github.com/DimAdam-01/pyrometheus.git@main"
Copy link
Member

Choose a reason for hiding this comment

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

you can't make master MFC use your fork of pyrometheus....

Copy link
Contributor Author

Choose a reason for hiding this comment

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

The original version encountered build errors on Frontier, so I introduced a few modifications to verify whether it would compile successfully there, and it did.

Copy link
Member

Choose a reason for hiding this comment

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

then we need to create a PR for the Pyrometheus fix first

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yes, I’m speaking with Esteban at the moment.

]

[tool.hatch.metadata]
Expand Down
Loading