Compare commits

...

10 Commits

85 changed files with 313 additions and 15334 deletions

9
.gitignore vendored
View File

@@ -69,8 +69,17 @@ subprojects/packagecache/
subprojects/hypre/
subprojects/qhull/
subprojects/libconstants/
subprojects/liblogging/
subprojects/libconfig/
subprojects/libcomposition/
subprojects/GridFire/
qhull.wrap
quill.wrap
yaml-cpp.wrap
cppad.wrap
subprojects/quill.wrap
.vscode/

View File

@@ -1 +0,0 @@
Use the utility `utils/atomic/convertWeightsToHeader.py` to generate atomicWeights.h

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

@@ -1,4 +0,0 @@
species_weight_dep = declare_dependency(
include_directories: include_directories('include'),
)
message('✅ SERiF species_weight dependency declared')

View File

@@ -1,469 +0,0 @@
CODATA 2022 + astrophysical constants
Most quantities have been converted to CGS
Generated on Feb 10, 2025
-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
Reference:
Mohr, P. , Tiesinga, E. , Newell, D. and Taylor, B. (2024), Codata Internationally Recommended 2022 Values of the Fundamental Physical Constants,
Codata Internationally Recommended 2022 Values of the Fundamental Physical Constants,
[online], https://tsapps.nist.gov/publication/get_pdf.cfm?pub_id=958002, https://physics.nist.gov/constants (Accessed February 10, 2025)
-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
Symbol Name Value Unit Uncertainty Source
-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
wienK Wien displacement law constant 2.89776850e-01 K cm 5.10000000e-07 CODATA2022
au1Hyper atomic unit of 1st hyperpolarizablity 3.20636151e-53 C^3 m^3 J^-2 2.80000000e-60 CODATA2022
au2Hyper atomic unit of 2nd hyperpolarizablity 6.23538080e-65 C^4 m^4 J^-3 1.10000000e-71 CODATA2022
auEDip atomic unit of electric dipole moment 8.47835309e-30 C m 7.30000000e-37 CODATA2022
auEPol atomic unit of electric polarizablity 1.64877727e-41 C^2 m^2 J^-1 1.60000000e-49 CODATA2022
auEQuad atomic unit of electric quadrupole moment 4.48655124e-40 C m^2 3.90000000e-47 CODATA2022
auMDip atomic unit of magn. dipole moment 1.85480190e-23 J T^-1 1.60000000e-30 CODATA2022
auMFlux atomic unit of magn. flux density 2.35051757e+09 G 7.10000000e-01 CODATA2022
muD deuteron magn. moment 4.33073482e-27 J T^-1 3.80000000e-34 CODATA2022
muD_Bhor deuteron magn. moment to Bohr magneton ratio 4.66975457e-04 5.00000000e-12 CODATA2022
muD_Nuc deuteron magn. moment to nuclear magneton ratio 8.57438233e-01 9.20000000e-09 CODATA2022
muD_e deuteron-electron magn. moment ratio -4.66434555e-04 5.00000000e-12 CODATA2022
muD_p deuteron-proton magn. moment ratio 3.07012208e-01 4.50000000e-09 CODATA2022
muD_n deuteron-neutron magn. moment ratio -4.48206520e-01 1.10000000e-07 CODATA2022
rgE electron gyromagn. ratio 1.76085963e+11 s^-1 T^-1 5.30000000e+01 CODATA2022
rgE_2pi electron gyromagn. ratio over 2 pi 2.80249532e+04 MHz T^-1 2.40000000e-03 CODATA2022
muE electron magn. moment -9.28476412e-24 J T^-1 8.00000000e-31 CODATA2022
muE_Bhor electron magn. moment to Bohr magneton ratio -1.00115965e+00 3.80000000e-12 CODATA2022
muE_Nuc electron magn. moment to nuclear magneton ratio -1.83828197e+03 8.50000000e-07 CODATA2022
muE_anom electron magn. moment anomaly 1.15965219e-03 3.80000000e-12 CODATA2022
muE_muP_shield electron to shielded proton magn. moment ratio -6.58227596e+02 7.10000000e-06 CODATA2022
muE_muH_shield electron to shielded helion magn. moment ratio 8.64058255e+02 1.00000000e-05 CODATA2022
muE_D electron-deuteron magn. moment ratio -2.14392349e+03 2.30000000e-05 CODATA2022
muE_mu electron-muon magn. moment ratio 2.06766989e+02 5.40000000e-06 CODATA2022
muE_n electron-neutron magn. moment ratio 9.60920500e+02 2.30000000e-04 CODATA2022
muE_p electron-proton magn. moment ratio -6.58210686e+02 6.60000000e-06 CODATA2022
mu0 magn. constant 1.25663706e-06 N A^-2 0.00000000e+00 CODATA2022
phi0 magn. flux quantum 2.06783385e-15 Wb 0.00000000e+00 CODATA2022
muMu muon magn. moment -4.49044799e-26 J T^-1 4.00000000e-33 CODATA2022
muMu_Bhor muon magn. moment to Bohr magneton ratio -4.84197045e-03 1.30000000e-10 CODATA2022
muMu_Nuc muon magn. moment to nuclear magneton ratio -8.89059698e+00 2.30000000e-07 CODATA2022
muMu_p muon-proton magn. moment ratio -3.18334512e+00 8.90000000e-08 CODATA2022
rgN neutron gyromagn. ratio 1.83247171e+08 s^-1 T^-1 4.30000000e+01 CODATA2022
rgN_2pi neutron gyromagn. ratio over 2 pi 2.91646950e+01 MHz T^-1 7.30000000e-06 CODATA2022
muN neutron magn. moment -9.66236450e-27 J T^-1 2.40000000e-33 CODATA2022
muN_Bhor neutron magn. moment to Bohr magneton ratio -1.04187563e-03 2.50000000e-10 CODATA2022
muN_Nuc neutron magn. moment to nuclear magneton ratio -1.91304273e+00 4.50000000e-07 CODATA2022
muN_p_shield neutron to shielded proton magn. moment ratio -6.84996940e-01 1.60000000e-07 CODATA2022
muN_e neutron-electron magn. moment ratio 1.04066882e-03 2.50000000e-10 CODATA2022
muN_p neutron-proton magn. moment ratio -6.84979340e-01 1.60000000e-07 CODATA2022
rgP proton gyromagn. ratio 2.67522187e+08 s^-1 T^-1 1.10000000e-01 CODATA2022
rgP_2pi proton gyromagn. ratio over 2 pi 4.25774813e+01 MHz T^-1 3.70000000e-06 CODATA2022
muP proton magn. moment 1.41060671e-26 J T^-1 1.20000000e-33 CODATA2022
muP_Bhor proton magn. moment to Bohr magneton ratio 1.52103221e-03 1.50000000e-11 CODATA2022
muP_Nuc proton magn. moment to nuclear magneton ratio 2.79284735e+00 2.80000000e-08 CODATA2022
muP_shield proton magn. shielding correction 2.56890000e-05 1.10000000e-08 CODATA2022
muP_n proton-neutron magn. moment ratio -1.45989805e+00 3.40000000e-07 CODATA2022
rgH shielded helion gyromagn. ratio 2.03789457e+08 s^-1 T^-1 2.40000000e+00 CODATA2022
rgH_2pi shielded helion gyromagn. ratio over 2 pi 3.24341015e+01 MHz T^-1 2.80000000e-06 CODATA2022
muH_shield shielded helion magn. moment -1.07455302e-26 J T^-1 9.30000000e-34 CODATA2022
muH_shield_Bhor shielded helion magn. moment to Bohr magneton rati -1.15867147e-03 1.40000000e-11 CODATA2022
muH_shield_Nuc shielded helion magn. moment to nuclear magneton r -2.12749772e+00 2.50000000e-08 CODATA2022
muH_shield_p shielded helion to proton magn. moment ratio -7.61766562e-01 1.20000000e-08 CODATA2022
muH_shield_p_shield shielded helion to shielded proton magn. moment ra -7.61786131e-01 3.30000000e-09 CODATA2022
muP_shield shielded proton magn. moment 1.41057047e-26 J T^-1 1.20000000e-33 CODATA2022
muP_shield_Bhor shielded proton magn. moment to Bohr magneton rati 1.52099313e-03 1.60000000e-11 CODATA2022
muP_shield_Nuc shielded proton magn. moment to nuclear magneton r 2.79277560e+00 3.00000000e-08 CODATA2022
aSi_220 {220} lattice spacing of silicon 1.92015571e-08 cm 3.20000000e-16 CODATA2022
aSi lattice spacing of silicon 1.92015576e-08 cm 5.00000000e-16 CODATA2022
mAlpha_e alpha particle-electron mass ratio 7.29429954e+03 2.40000000e-07 CODATA2022
mAlpha alpha particle mass 6.64465734e-24 g 2.00000000e-33 CODATA2022
EAlpha alpha particle mass energy equivalent 5.97192019e-03 erg 1.80000000e-12 CODATA2022
EAlpha_MeV alpha particle mass energy equivalent in MeV 5.97192019e-03 erg 1.76239430e-12 CODATA2022
mAlpha_u alpha particle mass in u 6.64465734e-24 g 1.04613961e-34 CODATA2022
MAlpha alpha particle molar mass 4.00150618e+00 g / mol 1.20000000e-09 CODATA2022
mAlpha_p alpha particle-proton mass ratio 3.97259969e+00 2.20000000e-10 CODATA2022
angstromStar Angstrom star 1.00001495e-08 cm 9.00000000e-15 CODATA2022
mU atomic mass constant 1.66053907e-24 g 5.00000000e-34 CODATA2022
mU_E atomic mass constant energy equivalent 1.49241809e-03 erg 4.50000000e-13 CODATA2022
mU_E_MeV atomic mass constant energy equivalent in MeV 1.49241809e-03 erg 4.48609458e-13 CODATA2022
mU_eV atomic mass unit-electron volt relationship 1.49241809e-03 erg 4.48609458e-13 CODATA2022
mU_hartree atomic mass unit-hartree relationship 3.42317769e+07 E_h 1.00000000e-02 CODATA2022
mU_Hz atomic mass unit-hertz relationship 2.25234272e+23 1 / s 6.80000000e+13 CODATA2022
mU_invm atomic mass unit-inverse meter relationship 7.51300661e+12 1 / cm 2.30000000e+03 CODATA2022
mU_J atomic mass unit-joule relationship 1.49241809e-03 erg 4.50000000e-13 CODATA2022
mU_K atomic mass unit-kelvin relationship 1.08095402e+13 K 3.30000000e+03 CODATA2022
mU_kg atomic mass unit-kilogram relationship 1.66053907e-24 g 5.00000000e-34 CODATA2022
au1Hyper atomic unit of 1st hyperpolarizability 3.20636131e-53 C^3 m^3 J^-2 1.50000000e-62 CODATA2022
au2Hyper atomic unit of 2nd hyperpolarizability 6.23537999e-65 C^4 m^4 J^-3 3.80000000e-74 CODATA2022
auAct atomic unit of action 1.05457182e-27 erg s 0.00000000e+00 CODATA2022
auC atomic unit of charge 1.60217663e-19 C 0.00000000e+00 CODATA2022
auCrho atomic unit of charge density 1.08120238e+12 C m^-3 4.90000000e+02 CODATA2022
auI atomic unit of current 6.62361824e-03 A 1.30000000e-14 CODATA2022
auMu atomic unit of electric dipole mom. 8.47835363e-30 C m 1.30000000e-39 CODATA2022
auE atomic unit of electric field 5.14220675e+11 V m^-1 7.80000000e+01 CODATA2022
auEFG atomic unit of electric field gradient 9.71736243e+21 V m^-2 2.90000000e+12 CODATA2022
auPol atomic unit of electric polarizability 1.64877727e-41 C^2 m^2 J^-1 5.00000000e-51 CODATA2022
auV atomic unit of electric potential 2.72113862e+01 V 5.30000000e-11 CODATA2022
auQuad atomic unit of electric quadrupole mom. 4.48655152e-40 C m^2 1.40000000e-49 CODATA2022
Eh atomic unit of energy 4.35974472e-11 erg 8.50000000e-23 CODATA2022
auF atomic unit of force 8.23872350e-03 dyn 1.20000000e-12 CODATA2022
auL atomic unit of length 5.29177211e-09 cm 8.00000000e-19 CODATA2022
auM atomic unit of mag. dipole mom. 1.85480202e-23 J T^-1 5.60000000e-33 CODATA2022
auB atomic unit of mag. flux density 2.35051757e+09 G 7.10000000e-01 CODATA2022
auChi atomic unit of magnetizability 7.89103660e-29 J T^-2 4.80000000e-38 CODATA2022
amu atomic unit of mass 9.10938370e-28 g 2.80000000e-37 CODATA2022
aup atomic unit of momentum 1.99285191e-19 dyn s 3.00000000e-29 CODATA2022
auEps atomic unit of permittivity 1.11265006e-10 F m^-1 1.70000000e-20 CODATA2022
aut atomic unit of time 2.41888433e-17 s 4.70000000e-29 CODATA2022
auv atomic unit of velocity 2.18769126e+08 cm / s 3.30000000e-02 CODATA2022
N_a Avogadro constant 6.02214076e+23 1 / mol 0.00000000e+00 CODATA2022
mu_B Bohr magneton 9.27401008e-24 J T^-1 2.80000000e-33 CODATA2022
mu_B_eV_T Bohr magneton in eV/T 5.78838181e-05 eV T^-1 1.70000000e-14 CODATA2022
mu_B_Hz_T Bohr magneton in Hz/T 1.39962449e+10 Hz T^-1 4.20000000e+00 CODATA2022
mu_B__mT Bohr magneton in inverse meters per tesla 4.66864481e+01 m^-1 T^-1 2.90000000e-07 CODATA2022
muB_K_T Bohr magneton in K/T 6.71713816e-01 K T^-1 2.00000000e-10 CODATA2022
rBhor Bohr radius 5.29177211e-09 cm 8.00000000e-19 CODATA2022
kB Boltzmann constant 1.38064900e-16 erg / K 0.00000000e+00 CODATA2022
kB_eV_K Boltzmann constant in eV/K 1.38064900e-16 erg / K 0.00000000e+00 CODATA2022
kB_Hz_K Boltzmann constant in Hz/K 2.08366191e+10 1 / (K s) 0.00000000e+00 CODATA2022
kB__MK Boltzmann constant in inverse meters per kelvin 6.95034570e-01 1 / (K cm) 4.00000000e-07 CODATA2022
Z0 characteristic impedance of vacuum 3.76730314e+02 ohm 5.61366546e-08 CODATA2022
re classical electron radius 2.81794033e-13 cm 1.30000000e-22 CODATA2022
lambda_compt Compton wavelength 2.42631024e-10 cm 7.30000000e-20 CODATA2022
lambda_compt_2pi Compton wavelength over 2 pi 3.86159268e-11 cm 1.80000000e-20 CODATA2022
G0 conductance quantum 7.74809173e-05 S 0.00000000e+00 CODATA2022
K_J-90 conventional value of Josephson constant 4.83597900e+14 Hz V^-1 0.00000000e+00 CODATA2022
R_K-90 conventional value of von Klitzing constant 2.58128070e+04 ohm 0.00000000e+00 CODATA2022
xu(CuKa1) Cu x unit 1.00207697e-11 cm 2.80000000e-18 CODATA2022
muD_e deuteron-electron mag. mom. ratio -4.66434555e-04 1.20000000e-12 CODATA2022
mD_e deuteron-electron mass ratio 3.67048297e+03 1.30000000e-07 CODATA2022
g_D deuteron g factor 8.57438234e-01 2.20000000e-09 CODATA2022
muD deuteron mag. mom. 4.33073509e-27 J T^-1 1.10000000e-35 CODATA2022
muD_Bhor deuteron mag. mom. to Bohr magneton ratio 4.66975457e-04 1.20000000e-12 CODATA2022
muD_Nuc deuteron mag. mom. to nuclear magneton ratio 8.57438234e-01 2.20000000e-09 CODATA2022
mD deuteron mass 3.34358377e-24 g 1.00000000e-33 CODATA2022
mD_E deuteron mass energy equivalent 3.00506323e-03 erg 9.10000000e-13 CODATA2022
mD_E_MeV deuteron mass energy equivalent in MeV 3.00506323e-03 erg 9.13240681e-13 CODATA2022
mD_u deuteron mass in u 3.34358377e-24 g 6.64215627e-35 CODATA2022
MD deuteron molar mass 2.01355321e+00 g / mol 6.10000000e-10 CODATA2022
muD_n deuteron-neutron mag. mom. ratio -4.48206530e-01 1.10000000e-07 CODATA2022
muD_p deuteron-proton mag. mom. ratio 3.07012209e-01 7.90000000e-10 CODATA2022
mD_p deuteron-proton mass ratio 1.99900750e+00 1.10000000e-10 CODATA2022
RrmsD deuteron rms charge radius 2.12799000e-13 cm 7.40000000e-17 CODATA2022
eps_0 electric constant 8.85418781e-12 F m^-1 1.30000000e-21 CODATA2022
qE_m electron charge to mass quotient -1.75882001e+11 C kg^-1 5.30000000e+01 CODATA2022
muE_D electron-deuteron mag. mom. ratio -2.14392349e+03 5.60000000e-06 CODATA2022
mE_D electron-deuteron mass ratio 2.72443711e-04 9.60000000e-15 CODATA2022
g_e electron g factor -2.00231930e+00 3.50000000e-13 CODATA2022
rg_e electron gyromag. ratio 1.76085963e+11 s^-1 T^-1 5.30000000e+01 CODATA2022
rg_e_2pi electron gyromag. ratio over 2 pi 2.80249516e+04 MHz T^-1 1.70000000e-04 CODATA2022
muE electron mag. mom. -9.28476470e-24 J T^-1 2.80000000e-33 CODATA2022
muE_anom electron mag. mom. anomaly 1.15965218e-03 1.80000000e-13 CODATA2022
muE_Bhor electron mag. mom. to Bohr magneton ratio -1.00115965e+00 1.80000000e-13 CODATA2022
muE_Nuc electron mag. mom. to nuclear magneton ratio -1.83828197e+03 1.10000000e-07 CODATA2022
mE electron mass 9.10938370e-28 g 2.80000000e-37 CODATA2022
mE_E electron mass energy equivalent 8.18710578e-07 erg 2.50000000e-16 CODATA2022
mE_MeV electron mass energy equivalent in MeV 8.18710578e-07 erg 2.40326495e-16 CODATA2022
mE_u electron mass in u 9.10938370e-28 g 2.65686251e-38 CODATA2022
ME electron molar mass 5.48579909e-04 g / mol 1.70000000e-13 CODATA2022
muE_mu electron-muon mag. mom. ratio 2.06766988e+02 4.60000000e-06 CODATA2022
mE_mu electron-muon mass ratio 4.83633169e-03 1.10000000e-10 CODATA2022
muE_n electron-neutron mag. mom. ratio 9.60920500e+02 2.30000000e-04 CODATA2022
mE_n electron-neutron mass ratio 5.43867344e-04 2.60000000e-13 CODATA2022
muE_p electron-proton mag. mom. ratio -6.58210688e+02 2.00000000e-07 CODATA2022
mE_p electron-proton mass ratio 5.44617021e-04 3.30000000e-14 CODATA2022
mE_tau electron-tau mass ratio 2.87585000e-04 1.90000000e-08 CODATA2022
mE_alpha electron to alpha particle mass ratio 1.37093355e-04 4.50000000e-15 CODATA2022
muE_H_shield electron to shielded helion mag. mom. ratio 8.64058257e+02 1.00000000e-05 CODATA2022
muE_p_shield electron to shielded proton mag. mom. ratio -6.58227597e+02 7.20000000e-06 CODATA2022
eV electron volt 1.60217663e-12 erg 0.00000000e+00 CODATA2022
eV_amu electron volt-atomic mass unit relationship 1.78266192e-33 g 5.31372501e-43 CODATA2022
eV_hartree electron volt-hartree relationship 3.67493222e-02 E_h 7.10000000e-14 CODATA2022
eV_Hz electron volt-hertz relationship 2.41798924e+14 1 / s 0.00000000e+00 CODATA2022
eV_invm electron volt-inverse meter relationship 8.06554394e+03 1 / cm 0.00000000e+00 CODATA2022
eV_J electron volt-joule relationship 1.60217663e-12 erg 0.00000000e+00 CODATA2022
eV_K electron volt-kelvin relationship 1.16045181e+04 K 0.00000000e+00 CODATA2022
eV_kg electron volt-kilogram relationship 1.78266192e-33 g 0.00000000e+00 CODATA2022
e elementary charge 1.60217663e-19 C 0.00000000e+00 CODATA2022
e_h elementary charge over h 2.41798926e+14 A J^-1 1.50000000e+06 CODATA2022
F Faraday constant 9.64853321e+04 C mol^-1 0.00000000e+00 CODATA2022
F_conv Faraday constant for conventional electric current 9.64853251e+04 C_90 mol^-1 1.20000000e-03 CODATA2022
G_F Fermi coupling constant 4.54379566e+00 s4 / (g2 cm4) 2.33738613e-06 CODATA2022
alpha fine-structure constant 7.29735257e-03 1.10000000e-12 CODATA2022
c1 first radiation constant 3.74177185e-05 St erg 0.00000000e+00 CODATA2022
c1L first radiation constant for spectral radiance 1.19104297e-05 St erg / rad2 0.00000000e+00 CODATA2022
Eh_amu hartree-atomic mass unit relationship 4.85087021e-32 g 1.46127438e-41 CODATA2022
Eh_eV hartree-electron volt relationship 4.35974472e-11 erg 8.49153616e-23 CODATA2022
Eh Hartree energy 4.35974472e-11 erg 8.50000000e-23 CODATA2022
Eh_E_eV Hartree energy in eV 4.35974472e-11 erg 8.49153616e-23 CODATA2022
Eh_Hz hartree-hertz relationship 6.57968392e+15 1 / s 1.30000000e+04 CODATA2022
Eh_invm hartree-inverse meter relationship 2.19474631e+05 1 / cm 4.30000000e-07 CODATA2022
Eh_J hartree-joule relationship 4.35974472e-11 erg 8.50000000e-23 CODATA2022
Eh_K hartree-kelvin relationship 3.15775025e+05 K 6.10000000e-07 CODATA2022
Eh_kg hartree-kilogram relationship 4.85087021e-32 g 9.40000000e-44 CODATA2022
mH_e helion-electron mass ratio 5.49588528e+03 2.40000000e-07 CODATA2022
mH helion mass 5.00641278e-24 g 1.50000000e-33 CODATA2022
mH_E helion mass energy equivalent 4.49953941e-03 erg 1.40000000e-12 CODATA2022
mH_E_eV helion mass energy equivalent in MeV 4.49953941e-03 erg 1.36185014e-12 CODATA2022
mH_u helion mass in u 5.00641278e-24 g 1.61072289e-34 CODATA2022
MH helion molar mass 3.01493225e+00 g / mol 9.10000000e-10 CODATA2022
mH_p helion-proton mass ratio 2.99315267e+00 1.30000000e-10 CODATA2022
Hz_amu hertz-atomic mass unit relationship 7.37249732e-48 g 2.15870079e-57 CODATA2022
Hz_eV hertz-electron volt relationship 6.62607015e-27 erg 0.00000000e+00 CODATA2022
Hz_Eh hertz-hartree relationship 1.51982985e-16 E_h 2.90000000e-28 CODATA2022
Hz_invm hertz-inverse meter relationship 3.33564095e-11 1 / cm 0.00000000e+00 CODATA2022
Hz_J hertz-joule relationship 6.62607015e-27 erg 0.00000000e+00 CODATA2022
Hz_K hertz-kelvin relationship 4.79924307e-11 K 0.00000000e+00 CODATA2022
Hz_kg hertz-kilogram relationship 7.37249732e-48 g 0.00000000e+00 CODATA2022
invalpha inverse fine-structure constant 1.37035999e+02 2.10000000e-08 CODATA2022
invm_amu inverse meter-atomic mass unit relationship 2.21021909e-39 g 6.64215627e-49 CODATA2022
invm_eV inverse meter-electron volt relationship 1.98644586e-18 erg 0.00000000e+00 CODATA2022
invm_Eh inverse meter-hartree relationship 4.55633525e-08 E_h 8.80000000e-20 CODATA2022
invm_Hz inverse meter-hertz relationship 2.99792458e+08 1 / s 0.00000000e+00 CODATA2022
invm_J inverse meter-joule relationship 1.98644586e-18 erg 0.00000000e+00 CODATA2022
invm_K inverse meter-kelvin relationship 1.43877688e-02 K 0.00000000e+00 CODATA2022
invm_kg inverse meter-kilogram relationship 2.21021909e-39 g 0.00000000e+00 CODATA2022
invG0 inverse of conductance quantum 1.29064037e+04 ohm 0.00000000e+00 CODATA2022
K_J Josephson constant 4.83597848e+14 Hz V^-1 0.00000000e+00 CODATA2022
J_amu joule-atomic mass unit relationship 1.11265006e-14 g 3.32107813e-24 CODATA2022
J_eV joule-electron volt relationship 1.00000000e+07 erg 0.00000000e+00 CODATA2022
J_Eh joule-hartree relationship 2.29371228e+17 E_h 4.50000000e+05 CODATA2022
J_Hz joule-hertz relationship 1.50919018e+33 1 / s 0.00000000e+00 CODATA2022
J_invm joule-inverse meter relationship 5.03411657e+22 1 / cm 0.00000000e+00 CODATA2022
J_K joule-kelvin relationship 7.24297052e+22 K 0.00000000e+00 CODATA2022
J_kg joule-kilogram relationship 1.11265006e-14 g 0.00000000e+00 CODATA2022
K_amu kelvin-atomic mass unit relationship 1.53617919e-37 g 4.64950939e-47 CODATA2022
K_eV kelvin-electron volt relationship 1.38064900e-16 erg 0.00000000e+00 CODATA2022
K_Eh kelvin-hartree relationship 3.16681156e-06 E_h 6.10000000e-18 CODATA2022
K_Hz kelvin-hertz relationship 2.08366191e+10 1 / s 0.00000000e+00 CODATA2022
K_invm kelvin-inverse meter relationship 6.95034800e-01 1 / cm 0.00000000e+00 CODATA2022
K_J kelvin-joule relationship 1.38064900e-16 erg 0.00000000e+00 CODATA2022
K_kg kelvin-kilogram relationship 1.53617919e-37 g 0.00000000e+00 CODATA2022
kg_amu kilogram-atomic mass unit relationship 1.00000000e+03 g 2.98897032e-07 CODATA2022
kg_eV kilogram-electron volt relationship 8.98755179e+23 erg 0.00000000e+00 CODATA2022
kg_Eh kilogram-hartree relationship 2.06148579e+34 E_h 4.00000000e+22 CODATA2022
kg_Hz kilogram-hertz relationship 1.35639249e+50 1 / s 0.00000000e+00 CODATA2022
kg_invm kilogram-inverse meter relationship 4.52443833e+39 1 / cm 0.00000000e+00 CODATA2022
kg_J kilogram-joule relationship 8.98755179e+23 erg 0.00000000e+00 CODATA2022
kg_K kilogram-kelvin relationship 6.50965726e+39 K 0.00000000e+00 CODATA2022
aSi lattice parameter of silicon 5.43102051e-08 cm 8.90000000e-16 CODATA2022
n0 Loschmidt constant (273.15 K, 101.325 kPa) 2.68678011e+19 1 / cm3 0.00000000e+00 CODATA2022
mu0 mag. constant 1.25663706e-06 N A^-2 1.90000000e-16 CODATA2022
Phi0 mag. flux quantum 2.06783385e-15 Wb 0.00000000e+00 CODATA2022
R molar gas constant 8.31446262e+07 erg / (K mol) 0.00000000e+00 CODATA2022
Mu molar mass constant 1.00000000e+00 g / mol 3.00000000e-10 CODATA2022
MC12 molar mass of carbon-12 1.20000000e+01 g / mol 3.60000000e-09 CODATA2022
h_mol molar Planck constant 3.99031271e-03 St g / mol 0.00000000e+00 CODATA2022
ch_mol molar Planck constant times c 1.19626566e+08 cm erg / mol 5.40000000e-02 CODATA2022
v_mol_ideal_100kPa molar volume of ideal gas (273.15 K, 100 kPa) 2.27109546e+04 cm3 / mol 0.00000000e+00 CODATA2022
v_mol_ideal_101.325kPa molar volume of ideal gas (273.15 K, 101.325 kPa) 2.24139695e+04 cm3 / mol 0.00000000e+00 CODATA2022
vSI_mol molar volume of silicon 1.20588320e+01 cm3 / mol 6.00000000e-07 CODATA2022
xu Mo x unit 1.00209952e-11 cm 5.30000000e-18 CODATA2022
lambda_compt_mu muon Compton wavelength 1.17344411e-12 cm 2.60000000e-20 CODATA2022
lambda_compt_mu_2pi muon Compton wavelength over 2 pi 1.86759431e-13 cm 4.20000000e-21 CODATA2022
mMu_e muon-electron mass ratio 2.06768283e+02 4.60000000e-06 CODATA2022
g_mu muon g factor -2.00233184e+00 1.30000000e-09 CODATA2022
muMu muon mag. mom. -4.49044830e-26 J T^-1 1.00000000e-33 CODATA2022
muMu_anom muon mag. mom. anomaly 1.16592089e-03 6.30000000e-10 CODATA2022
muMu_Bhor muon mag. mom. to Bohr magneton ratio -4.84197047e-03 1.10000000e-10 CODATA2022
muMu_Nuc muon mag. mom. to nuclear magneton ratio -8.89059703e+00 2.00000000e-07 CODATA2022
mMu muon mass 1.88353163e-25 g 4.20000000e-33 CODATA2022
mMu_E muon mass energy equivalent 1.69283380e-04 erg 3.80000000e-12 CODATA2022
mMu_E_MeV muon mass energy equivalent in MeV 1.69283380e-04 erg 3.68500626e-12 CODATA2022
mMu_E_u muon mass in u 1.88353163e-25 g 4.15134767e-33 CODATA2022
MMu muon molar mass 1.13428926e-01 g / mol 2.50000000e-09 CODATA2022
mMu_n muon-neutron mass ratio 1.12454517e-01 2.50000000e-09 CODATA2022
muMu_p muon-proton mag. mom. ratio -3.18334514e+00 7.10000000e-08 CODATA2022
mMu_p muon-proton mass ratio 1.12609526e-01 2.50000000e-09 CODATA2022
mMu_tau muon-tau mass ratio 5.94635000e-02 4.00000000e-06 CODATA2022
1 natural unit of action 1.05457182e-27 erg s 0.00000000e+00 CODATA2022
1_eV_s natural unit of action in eV s 1.05457182e-27 erg s 0.00000000e+00 CODATA2022
E natural unit of energy 8.18710578e-07 erg 2.50000000e-16 CODATA2022
E_MeV natural unit of energy in MeV 8.18710578e-07 erg 2.40326495e-16 CODATA2022
L natural unit of length 3.86159268e-11 cm 1.20000000e-20 CODATA2022
M natural unit of mass 9.10938370e-28 g 2.80000000e-37 CODATA2022
p natural unit of momentum 2.73092449e-17 dyn s 3.40000000e-25 CODATA2022
p_MeV_c natural unit of momentum in MeV/c 5.10998946e-01 MeV/c 3.10000000e-09 CODATA2022
t natural unit of time 1.28808867e-21 s 3.90000000e-31 CODATA2022
v natural unit of velocity 2.99792458e+10 cm / s 0.00000000e+00 CODATA2022
lambda_compt_n neutron Compton wavelength 1.31959091e-13 cm 7.50000000e-23 CODATA2022
lambda_compt_n_2pi neutron Compton wavelength over 2 pi 2.10019415e-14 cm 1.40000000e-23 CODATA2022
muN_e neutron-electron mag. mom. ratio 1.04066882e-03 2.50000000e-10 CODATA2022
mN_e neutron-electron mass ratio 1.83868366e+03 8.90000000e-07 CODATA2022
g_n neutron g factor -3.82608545e+00 9.00000000e-07 CODATA2022
rg_n neutron gyromag. ratio 1.83247171e+08 s^-1 T^-1 4.30000000e+01 CODATA2022
rg_n_2pi neutron gyromag. ratio over 2 pi 2.91646933e+01 MHz T^-1 6.90000000e-06 CODATA2022
muN neutron mag. mom. -9.66236510e-27 J T^-1 2.30000000e-33 CODATA2022
muN_Bhor neutron mag. mom. to Bohr magneton ratio -1.04187563e-03 2.50000000e-10 CODATA2022
muN_Nuc neutron mag. mom. to nuclear magneton ratio -1.91304273e+00 4.50000000e-07 CODATA2022
mN neutron mass 1.67492750e-24 g 9.50000000e-34 CODATA2022
mN_E neutron mass energy equivalent 1.50534976e-03 erg 8.60000000e-13 CODATA2022
mN_E_MeV neutron mass energy equivalent in MeV 1.50534976e-03 erg 8.65175382e-13 CODATA2022
mN_u neutron mass in u 1.67492750e-24 g 8.13664143e-34 CODATA2022
MN neutron molar mass 1.00866492e+00 g / mol 5.70000000e-10 CODATA2022
mN_mu neutron-muon mass ratio 8.89248406e+00 2.00000000e-07 CODATA2022
muN_p neutron-proton mag. mom. ratio -6.84979340e-01 1.60000000e-07 CODATA2022
mN_p neutron-proton mass ratio 1.00137842e+00 4.90000000e-10 CODATA2022
mN_tau neutron-tau mass ratio 5.28779000e-01 3.60000000e-05 CODATA2022
muN_p_shield neutron to shielded proton mag. mom. ratio -6.84996940e-01 1.60000000e-07 CODATA2022
G Newtonian constant of gravitation 6.67430000e-08 cm3 / (g s2) 1.50000000e-12 CODATA2022
G_hbar_c Newtonian constant of gravitation over h-bar c 6.70883000e-39 (GeV/c^2)^-2 1.50000000e-43 CODATA2022
mu_N nuclear magneton 5.05078375e-27 J T^-1 1.50000000e-36 CODATA2022
mu_N_eV_T nuclear magneton in eV/T 3.15245126e-08 eV T^-1 9.60000000e-18 CODATA2022
mu_N_invm_T nuclear magneton in inverse meters per tesla 2.54262343e-02 m^-1 T^-1 1.60000000e-10 CODATA2022
mu_N_K_T nuclear magneton in K/T 3.65826778e-04 K T^-1 1.10000000e-13 CODATA2022
mu_N_MHz_T nuclear magneton in MHz/T 7.62259323e+00 MHz T^-1 2.30000000e-09 CODATA2022
h Planck constant 6.62607015e-27 erg s 0.00000000e+00 CODATA2022
h_eV_s Planck constant in eV s 6.62607009e-27 erg s 4.00544159e-35 CODATA2022
hbar Planck constant over 2 pi 1.05457180e-27 erg s 1.30000000e-35 CODATA2022
hbar_eV_s Planck constant over 2 pi in eV s 1.05457181e-27 erg s 6.40870654e-36 CODATA2022
chbar_MeV Planck constant over 2 pi times c in MeV fm 3.16152675e-17 cm erg 1.92261196e-25 CODATA2022
planck_length Planck length 1.61625500e-33 cm 1.80000000e-38 CODATA2022
planck_mass Planck mass 2.17643400e-05 g 2.40000000e-10 CODATA2022
planck_mass_E_GeV Planck mass energy equivalent in GeV 1.95608143e+16 erg 2.24304729e+11 CODATA2022
planck_temp Planck temperature 1.41678400e+32 K 1.60000000e+27 CODATA2022
plank_time Planck time 5.39124700e-44 s 6.00000000e-49 CODATA2022
qP_mP proton charge to mass quotient 9.57883316e+07 C kg^-1 2.90000000e-02 CODATA2022
lambda_compt_p proton Compton wavelength 1.32140986e-13 cm 4.00000000e-23 CODATA2022
lambda_compt_p_2pi proton Compton wavelength over 2 pi 2.10308910e-14 cm 9.70000000e-24 CODATA2022
mP_e proton-electron mass ratio 1.83615267e+03 1.10000000e-07 CODATA2022
g_p proton g factor 5.58569469e+00 1.60000000e-09 CODATA2022
rg_p proton gyromag. ratio 2.67522187e+08 s^-1 T^-1 1.10000000e-01 CODATA2022
rg_p_2pi proton gyromag. ratio over 2 pi 4.25774789e+01 MHz T^-1 2.90000000e-07 CODATA2022
muP proton mag. mom. 1.41060680e-26 J T^-1 6.00000000e-36 CODATA2022
muP_Bhor proton mag. mom. to Bohr magneton ratio 1.52103220e-03 4.60000000e-13 CODATA2022
muP_Nuc proton mag. mom. to nuclear magneton ratio 2.79284734e+00 8.20000000e-10 CODATA2022
muP_shield proton mag. shielding correction 2.56890000e-05 1.10000000e-08 CODATA2022
mP proton mass 1.67262192e-24 g 5.10000000e-34 CODATA2022
mP_E proton mass energy equivalent 1.50327762e-03 erg 4.60000000e-13 CODATA2022
mP_E_MeV proton mass energy equivalent in MeV 1.50327762e-03 erg 4.64631224e-13 CODATA2022
mP_u proton mass in u 1.67262192e-24 g 8.80085705e-35 CODATA2022
MP proton molar mass 1.00727647e+00 g / mol 3.10000000e-10 CODATA2022
mP_mu proton-muon mass ratio 8.88024337e+00 2.00000000e-07 CODATA2022
muP_n proton-neutron mag. mom. ratio -1.45989805e+00 3.40000000e-07 CODATA2022
mP_n proton-neutron mass ratio 9.98623478e-01 4.90000000e-10 CODATA2022
RrmsP proton rms charge radius 8.41400000e-14 cm 1.90000000e-16 CODATA2022
mP_tau proton-tau mass ratio 5.28051000e-01 3.60000000e-05 CODATA2022
kappa quantum of circulation 3.63694755e+00 cm2 / s 1.10000000e-09 CODATA2022
2kappa quantum of circulation times 2 7.27389510e+00 cm2 / s 2.20000000e-09 CODATA2022
Rinf Rydberg constant 1.09737316e+05 1 / cm 2.10000000e-07 CODATA2022
cRinf_Hz Rydberg constant times c in Hz 3.28984196e+15 1 / s 6.40000000e+03 CODATA2022
hcRinf_eV Rydberg constant times hc in eV 2.17987236e-11 erg 4.16565925e-23 CODATA2022
hcRinf_J Rydberg constant times hc in J 2.17987236e-11 erg 4.20000000e-23 CODATA2022
S0_R_100kPa Sackur-Tetrode constant (1 K, 100 kPa) -1.15170754e+00 4.50000000e-10 CODATA2022
S0_R_101.325kPa Sackur-Tetrode constant (1 K, 101.325 kPa) -1.16487052e+00 4.50000000e-10 CODATA2022
c2 second radiation constant 1.43877688e+00 K cm 0.00000000e+00 CODATA2022
rg_h_shield shielded helion gyromag. ratio 2.03789457e+08 s^-1 T^-1 2.40000000e+00 CODATA2022
rg_h_shield_2pi shielded helion gyromag. ratio over 2 pi 3.24340997e+01 MHz T^-1 4.30000000e-07 CODATA2022
muH_shield shielded helion mag. mom. -1.07455309e-26 J T^-1 1.30000000e-34 CODATA2022
muH_shield_Bhor shielded helion mag. mom. to Bohr magneton ratio -1.15867147e-03 1.40000000e-11 CODATA2022
muH_shield_Nuc shielded helion mag. mom. to nuclear magneton rati -2.12749772e+00 2.50000000e-08 CODATA2022
muH_shield_p shielded helion to proton mag. mom. ratio -7.61766562e-01 8.90000000e-09 CODATA2022
muH_shield_p_shield shielded helion to shielded proton mag. mom. ratio -7.61786131e-01 3.30000000e-09 CODATA2022
rg_p_shield shielded proton gyromag. ratio 2.67515315e+08 s^-1 T^-1 2.90000000e+00 CODATA2022
rg_p_shield_2pi shielded proton gyromag. ratio over 2 pi 4.25763851e+01 MHz T^-1 5.30000000e-07 CODATA2022
muP_shield shielded proton mag. mom. 1.41057056e-26 J T^-1 1.50000000e-34 CODATA2022
muP_shield_Bhor shielded proton mag. mom. to Bohr magneton ratio 1.52099313e-03 1.70000000e-11 CODATA2022
muP_shield_Nuc shielded proton mag. mom. to nuclear magneton rati 2.79277560e+00 3.00000000e-08 CODATA2022
c speed of light in vacuum 2.99792458e+10 cm / s 0.00000000e+00 CODATA2022
g standard acceleration of gravity 9.80665000e+02 cm / s2 0.00000000e+00 CODATA2022
atm standard atmosphere 1.01325000e+06 P / s 0.00000000e+00 CODATA2022
sigma Stefan-Boltzmann constant 5.67037442e-05 g / (s3 K4) 0.00000000e+00 CODATA2022
lambda_compt_tau tau Compton wavelength 6.97771000e-14 cm 4.70000000e-18 CODATA2022
lambda_compt_tau_2pi tau Compton wavelength over 2 pi 1.11056000e-14 cm 1.00000000e-18 CODATA2022
mTau_e tau-electron mass ratio 3.47723000e+03 2.30000000e-01 CODATA2022
mTau tau mass 3.16754000e-24 g 2.10000000e-28 CODATA2022
mTau_E tau mass energy equivalent 2.84684000e-03 erg 1.90000000e-07 CODATA2022
mTau_E_MeV tau mass energy equivalent in MeV 2.84677949e-03 erg 2.56348261e-07 CODATA2022
mTau_u tau mass in u 3.16754469e-24 g 2.15870079e-28 CODATA2022
MTau tau molar mass 1.90754000e+00 g / mol 1.30000000e-04 CODATA2022
mTau_mu tau-muon mass ratio 1.68170000e+01 1.10000000e-03 CODATA2022
mTau_n tau-neutron mass ratio 1.89115000e+00 1.30000000e-04 CODATA2022
mTau_p tau-proton mass ratio 1.89376000e+00 1.30000000e-04 CODATA2022
sigma_T Thomson cross section 6.65245873e-25 cm2 6.00000000e-34 CODATA2022
muPsi_e triton-electron mag. mom. ratio -1.62051442e-03 2.10000000e-11 CODATA2022
mPsi_e triton-electron mass ratio 5.49692154e+03 2.70000000e-07 CODATA2022
g_psi triton g factor 5.95792493e+00 1.20000000e-08 CODATA2022
muPsi triton mag. mom. 1.50460952e-26 J T^-1 3.00000000e-35 CODATA2022
muPsi_Bhor triton mag. mom. to Bohr magneton ratio 1.62239367e-03 3.20000000e-12 CODATA2022
muPsi_Nuc triton mag. mom. to nuclear magneton ratio 2.97896247e+00 5.90000000e-09 CODATA2022
mPsi triton mass 5.00735674e-24 g 1.50000000e-33 CODATA2022
mPsi_E triton mass energy equivalent 4.50038781e-03 erg 1.40000000e-12 CODATA2022
mPsi_E_MeV triton mass energy equivalent in MeV 4.50038781e-03 erg 1.36185014e-12 CODATA2022
mPsi_u triton mass in u 5.00735674e-24 g 1.99264688e-34 CODATA2022
MPsi triton molar mass 3.01550072e+00 g / mol 9.20000000e-10 CODATA2022
muPsi_n triton-neutron mag. mom. ratio -1.55718553e+00 3.70000000e-07 CODATA2022
muPsi_p triton-proton mag. mom. ratio 1.06663991e+00 1.00000000e-08 CODATA2022
mPsi_p triton-proton mass ratio 2.99371703e+00 1.50000000e-10 CODATA2022
u unified atomic mass unit 1.66053907e-24 g 5.00000000e-34 CODATA2022
R_K von Klitzing constant 2.58128075e+04 ohm 0.00000000e+00 CODATA2022
theta_W weak mixing angle 2.22900000e-01 3.00000000e-04 CODATA2022
b_freq Wien frequency displacement law constant 5.87892576e+10 1 / (K s) 0.00000000e+00 CODATA2022
b_lambda Wien wavelength displacement law constant 2.89777196e-01 K cm 0.00000000e+00 CODATA2022
auMom_um atomic unit of mom.um 1.99285188e-19 dyn s 2.40000000e-27 CODATA2022
mE_h electron-helion mass ratio 1.81954307e-04 7.90000000e-15 CODATA2022
mE_psi electron-triton mass ratio 1.81920006e-04 9.00000000e-15 CODATA2022
g_h helion g factor -4.25525061e+00 5.00000000e-08 CODATA2022
muH helion mag. mom. -1.07461753e-26 J T^-1 1.30000000e-34 CODATA2022
muH_Bhor helion mag. mom. to Bohr magneton ratio -1.15874096e-03 1.40000000e-11 CODATA2022
muH_Nuc helion mag. mom. to nuclear magneton ratio -2.12762531e+00 2.50000000e-08 CODATA2022
n0 Loschmidt constant (273.15 K, 100 kPa) 2.65164580e+19 1 / cm3 0.00000000e+00 CODATA2022
p_um natural unit of mom.um 2.73092449e-17 dyn s 3.40000000e-25 CODATA2022
p_um_MeV_c natural unit of mom.um in MeV/c 5.10998946e-01 MeV/c 3.10000000e-09 CODATA2022
mN-mP neutron-proton mass difference 2.30557435e-27 g 8.20000000e-34 CODATA2022
mN-mP_E neutron-proton mass difference energy equivalent 2.07214689e-06 erg 7.40000000e-13 CODATA2022
mN-mP_E_MeV neutron-proton mass difference energy equivalent i 2.07214689e-06 erg 7.37001252e-13 CODATA2022
mN-mP_u neutron-proton mass difference in u 2.30557435e-27 g 8.13664143e-34 CODATA2022
stp standard-state pressure 1.00000000e+06 P / s 0.00000000e+00 CODATA2022
mAlpha alpha particle relative atomic mass 4.00150618e+00 6.30000000e-11 CODATA2022
muB_invm_T Bohr magneton in inverse meter per tesla 4.66864478e+01 m^-1 T^-1 1.40000000e-08 CODATA2022
kB_invm Boltzmann constant in inverse meter per kelvin 6.95034800e-01 1 / (K cm) 0.00000000e+00 CODATA2022
ampere-90 conventional value of ampere-90 1.00000009e+00 A 0.00000000e+00 CODATA2022
coulomb-90 conventional value of coulomb-90 1.00000009e+00 C 0.00000000e+00 CODATA2022
farad-90 conventional value of farad-90 9.99999982e-01 F 0.00000000e+00 CODATA2022
henry-90 conventional value of henry-90 1.00000002e+00 H 0.00000000e+00 CODATA2022
ohm-90 conventional value of ohm-90 1.00000002e+00 ohm 0.00000000e+00 CODATA2022
volt-90 conventional value of volt-90 1.00000011e+00 V 0.00000000e+00 CODATA2022
watt-90 conventional value of watt-90 1.00000020e+07 erg / s 0.00000000e+00 CODATA2022
mD_amu_rel deuteron relative atomic mass 2.01355321e+00 4.00000000e-11 CODATA2022
rg_e_MHz_T electron gyromag. ratio in MHz/T 2.80249514e+04 MHz T^-1 8.50000000e-06 CODATA2022
mE_amu_rel electron relative atomic mass 5.48579909e-04 1.60000000e-14 CODATA2022
e_hbar elementary charge over h-bar 1.51926745e+15 A J^-1 0.00000000e+00 CODATA2022
mH_amu_rel helion relative atomic mass 3.01493225e+00 9.70000000e-11 CODATA2022
h_shield_shift helion shielding shift 5.99674300e-05 1.00000000e-10 CODATA2022
F_hyperfine_Cs-133 hyperfine transition frequency of Cs-133 9.19263177e+09 1 / s 0.00000000e+00 CODATA2022
aSi_220_ideal lattice spacing of ideal Si (220) 1.92015572e-08 cm 3.20000000e-16 CODATA2022
eta luminous efficacy 6.83000000e-05 s3 rad2 cd/(g cm2) 0.00000000e+00 CODATA2022
rg_n_MHz_T neutron gyromag. ratio in MHz/T 2.91646931e+01 MHz T^-1 6.90000000e-06 CODATA2022
mN_amu_rel neutron relative atomic mass 1.00866492e+00 4.90000000e-10 CODATA2022
mu_N_invm nuclear magneton in inverse meter per tesla 2.54262341e-02 m^-1 T^-1 7.80000000e-12 CODATA2022
h_eV_Hz Planck constant in eV/Hz 6.62607015e-27 erg s 0.00000000e+00 CODATA2022
rg_p_MHz_T proton gyromag. ratio in MHz/T 4.25774785e+01 MHz T^-1 1.80000000e-08 CODATA2022
mP_amu_rel proton relative atomic mass 1.00727647e+00 5.30000000e-11 CODATA2022
lambda_compt_bar reduced Compton wavelength 3.86159268e-11 cm 1.20000000e-20 CODATA2022
lambda_compt_mu_bar reduced muon Compton wavelength 1.86759431e-13 cm 4.20000000e-21 CODATA2022
lambda_compt_n_bar reduced neutron Compton wavelength 2.10019416e-14 cm 1.20000000e-23 CODATA2022
hbar reduced Planck constant 1.05457182e-27 erg s 0.00000000e+00 CODATA2022
hbar_eV_s reduced Planck constant in eV s 1.05457182e-27 erg s 0.00000000e+00 CODATA2022
chbar_MeV reduced Planck constant times c in MeV fm 3.16152677e-17 cm erg 0.00000000e+00 CODATA2022
lambda_compt_p_bar reduced proton Compton wavelength 2.10308910e-14 cm 6.40000000e-24 CODATA2022
lambda_compt_tau_bar reduced tau Compton wavelength 1.11053800e-14 cm 7.50000000e-19 CODATA2022
rg_h_shield_MHz_T shielded helion gyromag. ratio in MHz/T 3.24340994e+01 MHz T^-1 3.80000000e-07 CODATA2022
rg_p_shield_MHz_T shielded proton gyromag. ratio in MHz/T 4.25763847e+01 MHz T^-1 4.60000000e-07 CODATA2022
d_shield-p_shield_HD shielding difference of d and p in HD 2.02000000e-08 2.00000000e-11 CODATA2022
t_shield-p_shield_HT shielding difference of t and p in HT 2.41400000e-08 2.00000000e-11 CODATA2022
tau_E tau energy equivalent 2.84684357e-03 erg 1.92261196e-07 CODATA2022
psi_amu_rel triton relative atomic mass 3.01550072e+00 1.20000000e-10 CODATA2022
muPsi_p triton to proton mag. mom. ratio 1.06663992e+00 2.10000000e-09 CODATA2022
epsilon0 vacuum electric permittivity 8.85418781e-12 F m^-1 1.30000000e-21 CODATA2022
mu0 vacuum mag. permeability 1.25663706e-06 N A^-2 1.90000000e-16 CODATA2022
mWZ W to Z mass ratio 8.81530000e-01 1.70000000e-04 CODATA2022
au Astronomical Unit 1.49597871e+13 cm 0.00000000e+00 IAU 2012 Resolution B2
pc Parsec 3.08567758e+18 cm 0.00000000e+00 Derived from au + IAU 2015 Resolution B 2 note [4]
kpc Kiloparsec 3.08567758e+21 cm 0.00000000e+00 Derived from au + IAU 2015 Resolution B 2 note [4]
L_bol0 Luminosity for absolute bolometric magnitude 0 3.01280000e+35 erg / s 0.00000000e+00 IAU 2015 Resolution B 2
L_sun Nominal solar luminosity 3.82800000e+33 erg / s 0.00000000e+00 IAU 2015 Resolution B 3
GM_sun Nominal solar mass parameter 1.32712440e+26 cm3 / s2 0.00000000e+00 IAU 2015 Resolution B 3
M_sun Solar mass 1.98840987e+33 g 4.46880543e+28 IAU 2015 Resolution B 3 + CODATA 2018
R_sun Nominal solar radius 6.95700000e+10 cm 0.00000000e+00 IAU 2015 Resolution B 3
GM_jup Nominal Jupiter mass parameter 1.26686530e+23 cm3 / s2 0.00000000e+00 IAU 2015 Resolution B 3
M_jup Jupiter mass 1.89812460e+30 g 4.26589589e+25 IAU 2015 Resolution B 3 + CODATA 2018
R_jup Nominal Jupiter equatorial radius 7.14920000e+09 cm 0.00000000e+00 IAU 2015 Resolution B 3
GM_earth Nominal Earth mass parameter 3.98600400e+20 cm3 / s2 0.00000000e+00 IAU 2015 Resolution B 3
M_earth Earth mass 5.97216787e+27 g 1.34220095e+23 IAU 2015 Resolution B 3 + CODATA 2018
R_earth Nominal Earth equatorial radius 6.37810000e+08 cm 0.00000000e+00 IAU 2015 Resolution B 3

View File

@@ -1,10 +0,0 @@
#!/bin/sh
if [ "$#" -ne 2 ]; then
echo "Usage: $0 <input_file> <output_file>"
exit 1
fi
input_file="$1"
output_file="$2"
printf 'const char embeddedConstants[] = R"(' > "$output_file"
cat "$input_file" >> "$output_file"
printf ')";\n' >> "$output_file"

View File

@@ -1,17 +0,0 @@
data_file = files('const.dat')
command_file = files('format.sh')
output_file = meson.current_build_dir() + '/embedded_constants.h'
embedded_constants_h = custom_target('embed_constants',
input: data_file,
output: 'embedded_constants.h',
command: ['sh', '-c', command_file[0].full_path()+' @INPUT@ ' + output_file, '@INPUT@', '@OUTPUT@']
)
# Ensure the generated header is included
const_data_header = include_directories('.')
const_data_dep = declare_dependency(
include_directories: const_data_header,
sources: embedded_constants_h
)

View File

@@ -1 +0,0 @@
subdir('atomic')

View File

@@ -0,0 +1,3 @@
# Disabled for dev
#gridfire_p = subproject('GridFire')
#gridfire_dep = gridfire_p.get_variable('gridfire_dep')

View File

@@ -0,0 +1,7 @@
composition_p = subproject('libcomposition')
comp_dep = composition_p.get_variable('composition_dep')
spw_dep = composition_p.get_variable('species_weight_dep')
composition_dep = [
comp_dep,
spw_dep,
]

View File

@@ -0,0 +1,2 @@
config_p = subproject('libconfig')
config_dep = config_p.get_variable('config_dep')

View File

@@ -0,0 +1,4 @@
logging_p = subproject('liblogging')
logging_dep = logging_p.get_variable('logging_dep')
quill_dep = logging_p.get_variable('quill_dep')

View File

@@ -1 +1,4 @@
subdir('libconstants')
subdir('liblogging')
subdir('libconfig')
subdir('libcomposition')

View File

@@ -1,10 +1,9 @@
cmake = import('cmake')
subdir('fourdst')
subdir('GridFire')
subdir('mfem')
subdir('yaml-cpp')
subdir('quill')
subdir('boost')
subdir('opatIO')
subdir('mpi')

View File

@@ -1,10 +0,0 @@
quill_cmake_options = cmake.subproject_options()
quill_cmake_options.add_cmake_defines({
'BUILD_SHARED_LIBS': 'ON',
'CMAKE_SKIP_INSTALL_RULES': 'ON'
})
quill_sp = cmake.subproject(
'quill',
options: quill_cmake_options,
)
quill_dep = quill_sp.dependency('quill')

View File

@@ -1,12 +0,0 @@
yaml_cpp_cmake_options = cmake.subproject_options()
yaml_cpp_cmake_options.add_cmake_defines({
'CMAKE_POLICY_VERSION_MINIMUM': '3.5',
'BUILD_SHARED_LIBS': 'ON',
'CMAKE_SKIP_INSTALL_RULES': 'ON',
'YAML_CPP_BUILD_TESTS': 'OFF'
})
yaml_cpp_sp = cmake.subproject(
'yaml-cpp',
options: yaml_cpp_cmake_options,
)
yaml_cpp_dep = yaml_cpp_sp.dependency('yaml-cpp')

View File

@@ -22,11 +22,10 @@ py_mod = py_installation.extension_module(
config_dep,
composition_dep,
eos_dep,
species_weight_dep,
mfem_dep,
polysolver_dep,
trampoline_dep,
network_dep,
gridfire_dep,
],
cpp_args : ['-UNDEBUG'], # Example: Ensure assertions are enabled if needed
install : true,

View File

@@ -18,7 +18,7 @@
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
#
# *********************************************************************** #
project('SERiF', 'cpp', version: '0.0.1a', default_options: ['cpp_std=c++23'], meson_version: '>=1.6.0')
project('SERiF', 'cpp', version: 'v0.0.1a', default_options: ['cpp_std=c++23'], meson_version: '>=1.5.0')
# Add default visibility for all C++ targets
add_project_arguments('-fvisibility=default', language: 'cpp')
@@ -37,12 +37,12 @@ endif
# Pass the DATA_DIR definition to the compiler
add_project_arguments('-DDATA_DIR=' + data_dir, language : 'cpp')
# We disable these because of CppAD
add_project_arguments('-Wno-bitwise-instead-of-logical', language: 'cpp')
# Build external dependencies first so that all the embedded resources are available to the other targets
subdir('build-config')
subdir('assets/static')
if get_option('build_debug_utils')
subdir('utils/debugUtils')
endif
@@ -54,9 +54,9 @@ if get_option('build_tests')
subdir('tests')
endif
if get_option('build_python')
subdir('build-python')
endif
# if get_option('build_python')
# subdir('build-python')
# endif
if get_option('build_post_run_utils')
subdir('utils')

View File

@@ -1,32 +0,0 @@
# Define the library
composition_sources = files(
'private/composition.cpp',
)
composition_headers = files(
'public/composition.h'
)
dependencies = [
probe_dep,
quill_dep,
species_weight_dep
]
# Define the libcomposition library so it can be linked against by other parts of the build system
libcomposition = library('composition',
composition_sources,
include_directories: include_directories('public'),
cpp_args: ['-fvisibility=default'],
dependencies: dependencies,
install : true)
composition_dep = declare_dependency(
include_directories: include_directories('public'),
link_with: libcomposition,
dependencies: dependencies,
sources: composition_sources,
)
# Make headers accessible
install_headers(composition_headers, subdir : 'SERiF/composition')

View File

@@ -1,668 +0,0 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: March 26, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#include "composition.h"
#include "quill/LogMacros.h"
#include <stdexcept>
#include <unordered_map>
#include <vector>
#include <array>
#include <ranges>
#include <utility>
#include "atomicSpecies.h"
namespace serif::composition {
CompositionEntry::CompositionEntry() : m_symbol("H-1"), m_isotope(chemSpecies::species.at("H-1")),
m_initialized(false) {
}
CompositionEntry::CompositionEntry(const std::string& symbol, const bool massFracMode) : m_symbol(symbol), m_isotope(chemSpecies::species.at(symbol)), m_massFracMode(massFracMode) {
setSpecies(symbol);
}
CompositionEntry::CompositionEntry(const CompositionEntry& entry) :
m_symbol(entry.m_symbol),
m_isotope(entry.m_isotope),
m_massFracMode(entry.m_massFracMode),
m_massFraction(entry.m_massFraction),
m_numberFraction(entry.m_numberFraction),
m_relAbundance(entry.m_relAbundance),
m_initialized(entry.m_initialized) {}
void CompositionEntry::setSpecies(const std::string& symbol) {
if (m_initialized) {
throw std::runtime_error("Composition entry is already initialized.");
}
if (chemSpecies::species.count(symbol) == 0) {
throw std::runtime_error("Invalid symbol.");
}
m_symbol = symbol;
m_isotope = chemSpecies::species.at(symbol);
m_initialized = true;
}
std::string CompositionEntry::symbol() const {
return m_symbol;
}
double CompositionEntry::mass_fraction() const {
if (!m_massFracMode) {
throw std::runtime_error("Composition entry is in number fraction mode.");
}
return m_massFraction;
}
double CompositionEntry::mass_fraction(double meanMolarMass) const {
if (m_massFracMode) {
return m_massFraction;
}
return m_relAbundance / meanMolarMass;
}
double CompositionEntry::number_fraction() const {
if (m_massFracMode) {
throw std::runtime_error("Composition entry is in mass fraction mode.");
}
return m_numberFraction;
}
double CompositionEntry::number_fraction(double totalMoles) const {
if (m_massFracMode) {
return m_relAbundance / totalMoles;
}
return m_numberFraction;
}
double CompositionEntry::rel_abundance() const {
return m_relAbundance;
}
chemSpecies::Species CompositionEntry::isotope() const {
return m_isotope;
}
void CompositionEntry::setMassFraction(double mass_fraction) {
if (!m_massFracMode) {
throw std::runtime_error("Composition entry is in number fraction mode.");
}
m_massFraction = mass_fraction;
m_relAbundance = m_massFraction / m_isotope.mass();
}
void CompositionEntry::setNumberFraction(double number_fraction) {
if (m_massFracMode) {
throw std::runtime_error("Composition entry is in mass fraction mode.");
}
m_numberFraction = number_fraction;
m_relAbundance = m_numberFraction * m_isotope.mass();
}
bool CompositionEntry::setMassFracMode(double meanParticleMass) {
if (m_massFracMode) {
return false;
}
m_massFracMode = true;
m_massFraction = m_relAbundance / meanParticleMass;
return true;
}
bool CompositionEntry::setNumberFracMode(double specificNumberDensity) {
if (!m_massFracMode) {
return false;
}
m_massFracMode = false;
m_numberFraction = m_relAbundance / specificNumberDensity;
return true;
}
bool CompositionEntry::getMassFracMode() const {
return m_massFracMode;
}
Composition::Composition(const std::vector<std::string>& symbols) {
for (const auto& symbol : symbols) {
registerSymbol(symbol);
}
}
Composition::Composition(const std::set<std::string>& symbols) {
for (const auto& symbol : symbols) {
registerSymbol(symbol);
}
}
Composition::Composition(const std::vector<std::string>& symbols, const std::vector<double>& fractions, bool massFracMode) : m_massFracMode(massFracMode) {
if (symbols.size() != fractions.size()) {
LOG_ERROR(m_logger, "The number of symbols and fractions must be equal.");
throw std::runtime_error("The number of symbols and fractions must be equal.");
}
validateComposition(fractions);
for (const auto &symbol : symbols) {
registerSymbol(symbol);
}
for (size_t i = 0; i < symbols.size(); ++i) {
if (m_massFracMode) {
setMassFraction(symbols[i], fractions[i]);
} else {
setNumberFraction(symbols[i], fractions[i]);
}
}
finalize();
}
Composition::Composition(const Composition &composition) {
m_finalized = composition.m_finalized;
m_specificNumberDensity = composition.m_specificNumberDensity;
m_meanParticleMass = composition.m_meanParticleMass;
m_massFracMode = composition.m_massFracMode;
m_registeredSymbols = composition.m_registeredSymbols;
m_compositions = composition.m_compositions;
}
Composition& Composition::operator=(const Composition &other) {
if (this != &other) {
m_finalized = other.m_finalized;
m_specificNumberDensity = other.m_specificNumberDensity;
m_meanParticleMass = other.m_meanParticleMass;
m_massFracMode = other.m_massFracMode;
m_registeredSymbols = other.m_registeredSymbols;
m_compositions = other.m_compositions;
// note: m_config remains bound to the same singleton, so we skip it
}
return *this;
}
void Composition::registerSymbol(const std::string& symbol, bool massFracMode) {
if (!isValidSymbol(symbol)) {
LOG_ERROR(m_logger, "Invalid symbol: {}", symbol);
throw std::runtime_error("Invalid symbol.");
}
// If no symbols have been registered allow mode to be set
if (m_registeredSymbols.size() == 0) {
m_massFracMode = massFracMode;
} else {
if (m_massFracMode != massFracMode) {
LOG_ERROR(m_logger, "Composition is in mass fraction mode. Cannot register symbol in number fraction mode.");
throw std::runtime_error("Composition is in mass fraction mode. Cannot register symbol in number fraction mode.");
}
}
if (m_registeredSymbols.find(symbol) != m_registeredSymbols.end()) {
LOG_WARNING(m_logger, "Symbol {} is already registered.", symbol);
return;
}
m_registeredSymbols.insert(symbol);
CompositionEntry entry(symbol, m_massFracMode);
m_compositions[symbol] = entry;
LOG_INFO(m_logger, "Registered symbol: {}", symbol);
}
void Composition::registerSymbol(const std::vector<std::string>& symbols, bool massFracMode) {
for (const auto& symbol : symbols) {
registerSymbol(symbol, massFracMode);
}
}
std::set<std::string> Composition::getRegisteredSymbols() const {
return m_registeredSymbols;
}
void Composition::validateComposition(const std::vector<double>& fractions) const {
if (!isValidComposition(fractions)) {
LOG_ERROR(m_logger, "Invalid composition.");
throw std::runtime_error("Invalid composition.");
}
}
bool Composition::isValidComposition(const std::vector<double>& fractions) const {
double sum = 0.0;
for (const auto& fraction : fractions) {
sum += fraction;
}
if (sum < 0.999999 || sum > 1.000001) {
LOG_ERROR(m_logger, "The sum of fractions must be equal to 1.");
return false;
}
return true;
}
bool Composition::isValidSymbol(const std::string& symbol) {
return chemSpecies::species.contains(symbol);
}
double Composition::setMassFraction(const std::string& symbol, const double& mass_fraction) {
if (!m_registeredSymbols.contains(symbol)) {
LOG_ERROR(m_logger, "Symbol {} is not registered.", symbol);
throw std::runtime_error("Symbol is not registered.");
}
if (!m_massFracMode) {
LOG_ERROR(m_logger, "Composition is in number fraction mode.");
throw std::runtime_error("Composition is in number fraction mode.");
}
if (mass_fraction < 0.0 || mass_fraction > 1.0) {
LOG_ERROR(m_logger, "Mass fraction must be between 0 and 1 for symbol {}. Currently it is {}.", symbol, mass_fraction);
throw std::runtime_error("Mass fraction must be between 0 and 1.");
}
m_finalized = false;
const double old_mass_fraction = m_compositions.at(symbol).mass_fraction();
m_compositions.at(symbol).setMassFraction(mass_fraction);
return old_mass_fraction;
}
std::vector<double> Composition::setMassFraction(const std::vector<std::string>& symbols, const std::vector<double>& mass_fractions) {
if (symbols.size() != mass_fractions.size()) {
LOG_ERROR(m_logger, "The number of symbols and mass fractions must be equal.");
throw std::runtime_error("The number of symbols and mass fractions must be equal.");
}
std::vector<double> old_mass_fractions;
old_mass_fractions.reserve(symbols.size());
for (size_t i = 0; i < symbols.size(); ++i) {
old_mass_fractions.push_back(setMassFraction(symbols[i], mass_fractions[i]));
}
return old_mass_fractions;
}
double Composition::setNumberFraction(const std::string& symbol, const double& number_fraction) {
if (m_registeredSymbols.find(symbol) == m_registeredSymbols.end()) {
LOG_ERROR(m_logger, "Symbol {} is not registered.", symbol);
throw std::runtime_error("Symbol is not registered.");
}
if (m_massFracMode) {
LOG_ERROR(m_logger, "Composition is in mass fraction mode.");
throw std::runtime_error("Composition is in mass fraction mode.");
}
if (number_fraction < 0.0 || number_fraction > 1.0) {
LOG_ERROR(m_logger, "Number fraction must be between 0 and 1 for symbol {}. Currently it is {}.", symbol, number_fraction);
throw std::runtime_error("Number fraction must be between 0 and 1.");
}
m_finalized = false;
double old_number_fraction = m_compositions.at(symbol).number_fraction();
m_compositions.at(symbol).setNumberFraction(number_fraction);
return old_number_fraction;
}
std::vector<double> Composition::setNumberFraction(const std::vector<std::string>& symbols, const std::vector<double>& number_fractions) {
if (symbols.size() != number_fractions.size()) {
LOG_ERROR(m_logger, "The number of symbols and number fractions must be equal.");
throw std::runtime_error("The number of symbols and number fractions must be equal.");
}
std::vector<double> old_number_fractions;
old_number_fractions.reserve(symbols.size());
for (size_t i = 0; i < symbols.size(); ++i) {
old_number_fractions.push_back(setNumberFraction(symbols[i], number_fractions[i]));
}
return old_number_fractions;
}
bool Composition::finalize(const bool norm) {
bool finalized = false;
if (m_massFracMode) {
finalized = finalizeMassFracMode(norm);
} else {
finalized = finalizeNumberFracMode(norm);
}
if (finalized) {
m_finalized = true;
}
return finalized;
}
bool Composition::finalizeMassFracMode(bool norm) {
std::vector<double> mass_fractions;
mass_fractions.reserve(m_compositions.size());
for (const auto& [_, entry] : m_compositions) {
mass_fractions.push_back(entry.mass_fraction());
}
if (norm) {
double sum = 0.0;
for (const auto& mass_fraction : mass_fractions) {
sum += mass_fraction;
}
for (int i = 0; i < static_cast<int>(mass_fractions.size()); ++i) {
mass_fractions[i] /= sum;
}
for (auto& [symbol, entry] : m_compositions) {
setMassFraction(symbol, entry.mass_fraction() / sum);
}
}
try {
validateComposition(mass_fractions);
} catch (const std::runtime_error& e) {
double massSum = 0.0;
for (const auto& [_, entry] : m_compositions) {
massSum += entry.mass_fraction();
}
LOG_ERROR(m_logger, "Composition is invalid (Total mass {}).", massSum);
m_finalized = false;
return false;
}
for (const auto& [_, entry] : m_compositions) {
m_specificNumberDensity += entry.rel_abundance();
}
m_meanParticleMass = 1.0/m_specificNumberDensity;
return true;
}
bool Composition::finalizeNumberFracMode(bool norm) {
std::vector<double> number_fractions;
number_fractions.reserve(m_compositions.size());
for (const auto& [_, entry] : m_compositions) {
number_fractions.push_back(entry.number_fraction());
}
if (norm) {
double sum = 0.0;
for (const auto& number_fraction : number_fractions) {
sum += number_fraction;
}
for (auto& [symbol, entry] : m_compositions) {
setNumberFraction(symbol, entry.number_fraction() / sum);
}
}
try {
validateComposition(number_fractions);
} catch (const std::runtime_error& e) {
double numberSum = 0.0;
for (const auto& [_, entry] : m_compositions) {
numberSum += entry.number_fraction();
}
LOG_ERROR(m_logger, "Composition is invalid (Total number {}).", numberSum);
m_finalized = false;
return false;
}
for (const auto& [_, entry] : m_compositions) {
m_meanParticleMass += entry.rel_abundance();
}
m_specificNumberDensity = 1.0/m_meanParticleMass;
return true;
}
Composition Composition::mix(const Composition& other, double fraction) const {
if (!m_finalized || !other.m_finalized) {
LOG_ERROR(m_logger, "Compositions have not both been finalized.");
throw std::runtime_error("Compositions have not been finalized (Consider running .finalize()).");
}
if (fraction < 0.0 || fraction > 1.0) {
LOG_ERROR(m_logger, "Fraction must be between 0 and 1.");
throw std::runtime_error("Fraction must be between 0 and 1.");
}
std::set<std::string> mixedSymbols = other.getRegisteredSymbols();
// Get the union of the two sets
mixedSymbols.insert(m_registeredSymbols.begin(), m_registeredSymbols.end());
Composition mixedComposition(mixedSymbols);
for (const auto& symbol : mixedSymbols) {
double otherMassFrac = 0.0;
const double thisMassFrac = hasSymbol(symbol) ? getMassFraction(symbol) : 0.0;
otherMassFrac = other.hasSymbol(symbol) ? other.getMassFraction(symbol) : 0.0;
double massFraction = fraction * thisMassFrac + otherMassFrac * (1-fraction);
mixedComposition.setMassFraction(symbol, massFraction);
}
mixedComposition.finalize();
return mixedComposition;
}
double Composition::getMassFraction(const std::string& symbol) const {
if (!m_finalized) {
LOG_ERROR(m_logger, "Composition has not been finalized.");
throw std::runtime_error("Composition has not been finalized (Consider running .finalize()).");
}
if (!m_compositions.contains(symbol)) {
LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol);
throw std::runtime_error("Symbol is not in the composition.");
}
if (m_massFracMode) {
return m_compositions.at(symbol).mass_fraction();
} else {
return m_compositions.at(symbol).mass_fraction(m_meanParticleMass);
}
}
std::unordered_map<std::string, double> Composition::getMassFraction() const {
std::unordered_map<std::string, double> mass_fractions;
for (const auto &symbol: m_compositions | std::views::keys) {
mass_fractions[symbol] = getMassFraction(symbol);
}
return mass_fractions;
}
double Composition::getNumberFraction(const std::string& symbol) const {
if (!m_finalized) {
LOG_ERROR(m_logger, "Composition has not been finalized.");
throw std::runtime_error("Composition has not been finalized (Consider running .finalize()).");
}
if (!m_compositions.contains(symbol)) {
LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol);
throw std::runtime_error("Symbol is not in the composition.");
}
if (!m_massFracMode) {
return m_compositions.at(symbol).number_fraction();
} else {
return m_compositions.at(symbol).number_fraction(m_specificNumberDensity);
}
}
std::unordered_map<std::string, double> Composition::getNumberFraction() const {
std::unordered_map<std::string, double> number_fractions;
for (const auto &symbol: m_compositions | std::views::keys) {
number_fractions[symbol] = getNumberFraction(symbol);
}
return number_fractions;
}
std::pair<CompositionEntry, GlobalComposition> Composition::getComposition(const std::string& symbol) const {
if (!m_finalized) {
LOG_ERROR(m_logger, "Composition has not been finalized.");
throw std::runtime_error("Composition has not been finalized (Consider running .finalize()).");
}
if (!m_compositions.contains(symbol)) {
LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol);
throw std::runtime_error("Symbol is not in the composition.");
}
return {m_compositions.at(symbol), {m_specificNumberDensity, m_meanParticleMass}};
}
std::pair<std::unordered_map<std::string, CompositionEntry>, GlobalComposition> Composition::getComposition() const {
if (!m_finalized) {
LOG_ERROR(m_logger, "Composition has not been finalized.");
throw std::runtime_error("Composition has not been finalized (Consider running .finalize()).");
}
return {m_compositions, {m_specificNumberDensity, m_meanParticleMass}};
}
double Composition::getMeanParticleMass() const {
if (!m_finalized) {
LOG_ERROR(m_logger, "Composition has not been finalized.");
throw std::runtime_error("Composition has not been finalized (Consider running .finalize()).");
}
return m_meanParticleMass;
}
double Composition::getMeanAtomicNumber() const {
if (!m_finalized) {
LOG_ERROR(m_logger, "Composition must be finalized before getting the mean atomic mass number.");
throw std::runtime_error("Composition not finalized. Cannot retrieve mean atomic mass number.");
}
double zSum = 0.0;
// Loop through all registered species in the composition.
for (const auto &val: m_compositions | std::views::values) {
zSum += (val.mass_fraction() * val.m_isotope.z())/val.m_isotope.a();
}
const double mean_A = m_meanParticleMass * zSum;
return mean_A;
}
Composition Composition::subset(const std::vector<std::string>& symbols, std::string method) const {
const std::array<std::string, 2> methods = {"norm", "none"};
if (std::ranges::find(methods, method) == methods.end()) {
const std::string errorMessage = "Invalid method: " + method + ". Valid methods are 'norm' and 'none'.";
LOG_ERROR(m_logger, "Invalid method: {}. Valid methods are norm and none.", method);
throw std::runtime_error(errorMessage);
}
Composition subsetComposition;
for (const auto& symbol : symbols) {
if (!m_compositions.contains(symbol)) {
LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol);
throw std::runtime_error("Symbol is not in the composition.");
} else {
subsetComposition.registerSymbol(symbol);
}
subsetComposition.setMassFraction(symbol, m_compositions.at(symbol).mass_fraction());
}
if (method == "norm") {
const bool isNorm = subsetComposition.finalize(true);
if (!isNorm) {
LOG_ERROR(m_logger, "Subset composition is invalid.");
throw std::runtime_error("Subset composition is invalid.");
}
}
return subsetComposition;
}
void Composition::setCompositionMode(const bool massFracMode) {
if (!m_finalized) {
LOG_ERROR(m_logger, "Composition has not been finalized. Mode cannot be set unless composition is finalized.");
throw std::runtime_error("Composition has not been finalized (Consider running .finalize()). The mode cannot be set unless the composition is finalized.");
}
bool okay = true;
for (auto &entry: m_compositions | std::views::values) {
if (massFracMode) {
okay = entry.setMassFracMode(m_meanParticleMass);
} else {
okay = entry.setNumberFracMode(m_specificNumberDensity);
}
if (!okay) {
LOG_ERROR(m_logger, "Composition mode could not be set.");
throw std::runtime_error("Composition mode could not be set due to an unknown error.");
}
}
m_massFracMode = massFracMode;
}
CanonicalComposition Composition::getCanonicalComposition(bool harsh) const {
if (!m_finalized) {
LOG_ERROR(m_logger, "Composition has not been finalized.");
throw std::runtime_error("Composition has not been finalized (Consider running .finalize()).");
}
CanonicalComposition canonicalComposition;
const std::array<std::string, 7> canonicalH = {
"H-1", "H-2", "H-3", "H-4", "H-5", "H-6", "H-7"
};
const std::array<std::string, 8> canonicalHe = {
"He-3", "He-4", "He-5", "He-6", "He-7", "He-8", "He-9", "He-10"
};
for (const auto& symbol : canonicalH) {
if (hasSymbol(symbol)) {
canonicalComposition.X += getMassFraction(symbol);
}
}
for (const auto& symbol : canonicalHe) {
if (hasSymbol(symbol)) {
canonicalComposition.Y += getMassFraction(symbol);
}
}
for (const auto& symbol : getRegisteredSymbols()) {
const bool isHSymbol = std::ranges::find(canonicalH, symbol) != std::end(canonicalH);
const bool isHeSymbol = std::ranges::find(canonicalHe, symbol) != std::end(canonicalHe);
if (isHSymbol || isHeSymbol) {
continue; // Skip canonical H and He symbols
}
canonicalComposition.Z += getMassFraction(symbol);
}
const double Z = 1.0 - (canonicalComposition.X + canonicalComposition.Y);
if (std::abs(Z - canonicalComposition.Z) > 1e-6) {
if (!harsh) {
LOG_WARNING(m_logger, "Validation composition Z (X-Y = {}) is different than canonical composition Z ({}) (∑a_i where a_i != H/He).", Z, canonicalComposition.Z);
}
else {
LOG_ERROR(m_logger, "Validation composition Z (X-Y = {}) is different than canonical composition Z ({}) (∑a_i where a_i != H/He).", Z, canonicalComposition.Z);
throw std::runtime_error("Validation composition Z (X-Y = " + std::to_string(Z) + ") is different than canonical composition Z (" + std::to_string(canonicalComposition.Z) + ") (∑a_i where a_i != H/He).");
}
}
return canonicalComposition;
}
bool Composition::hasSymbol(const std::string& symbol) const {
return m_compositions.count(symbol) > 0;
}
/// OVERLOADS
Composition Composition::operator+(const Composition& other) const {
return mix(other, 0.5);
}
std::ostream& operator<<(std::ostream& os, const GlobalComposition& comp) {
os << "Global Composition: \n";
os << "\tSpecific Number Density: " << comp.specificNumberDensity << "\n";
os << "\tMean Particle Mass: " << comp.meanParticleMass << "\n";
return os;
}
std::ostream& operator<<(std::ostream& os, const CompositionEntry& entry) {
os << "<" << entry.m_symbol << " : m_frac = " << entry.mass_fraction() << ">";
return os;
}
std::ostream& operator<<(std::ostream& os, const Composition& composition) {
os << "Composition: \n";
for (const auto& [symbol, entry] : composition.m_compositions) {
os << entry << "\n";
}
return os;
}
} // namespace serif::composition

View File

@@ -1,509 +0,0 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: March 26, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#pragma once
#include <iostream>
#include <string>
#include <unordered_map>
#include <set>
#include <utility>
#include "probe.h"
#include "config.h"
#include "atomicSpecies.h"
namespace serif::composition {
struct CanonicalComposition {
double X = 0.0; ///< Mass fraction of Hydrogen.
double Y = 0.0; ///< Mass fraction of Helium.
double Z = 0.0; ///< Mass fraction of Metals.
friend std::ostream& operator<<(std::ostream& os, const CanonicalComposition& composition) {
os << "<CanonicalComposition: "
<< "X = " << composition.X << ", "
<< "Y = " << composition.Y << ", "
<< "Z = " << composition.Z << ">";
return os;
}
};
/**
* @brief Represents the global composition of a system. This tends to be used after finalize and is primarily for internal use.
*/
struct GlobalComposition {
double specificNumberDensity; ///< The specific number density of the composition (\sum_{i} X_i m_i. Where X_i is the number fraction of the ith species and m_i is the mass of the ith species).
double meanParticleMass; ///< The mean particle mass of the composition (\sum_{i} \frac{n_i}{m_i}. where n_i is the number fraction of the ith species and m_i is the mass of the ith species).
// Overload the output stream operator for GlobalComposition
friend std::ostream& operator<<(std::ostream& os, const GlobalComposition& comp);
};
/**
* @brief Represents an entry in the composition with a symbol and mass fraction.
*/
struct CompositionEntry {
std::string m_symbol; ///< The chemical symbol of the species.
chemSpecies::Species m_isotope; ///< The isotope of the species.
bool m_massFracMode = true; ///< The mode of the composition entry. True if mass fraction, false if number fraction.
double m_massFraction = 0.0; ///< The mass fraction of the species.
double m_numberFraction = 0.0; ///< The number fraction of the species.
double m_relAbundance = 0.0; ///< The relative abundance of the species for converting between mass and number fractions.
bool m_initialized = false; ///< True if the composition entry has been initialized.
/**
* @brief Default constructor.
*/
CompositionEntry();
/**
* @brief Constructs a CompositionEntry with the given symbol and mode.
* @param symbol The chemical symbol of the species.
* @param massFracMode True if mass fraction mode, false if number fraction mode.
* *Example Usage:*
* @code
* CompositionEntry entry("H", true);
* @endcode
*/
CompositionEntry(const std::string& symbol, bool massFracMode=true);
/**
* @brief Copy constructor.
* @param entry The CompositionEntry to copy.
*/
CompositionEntry(const CompositionEntry& entry);
/**
* @brief Sets the species for the composition entry.
* @param symbol The chemical symbol of the species.
*/
void setSpecies(const std::string& symbol);
/**
* @brief Gets the chemical symbol of the species.
* @return The chemical symbol of the species.
*/
std::string symbol() const;
/**
* @brief Gets the mass fraction of the species.
* @return The mass fraction of the species.
*/
double mass_fraction() const;
/**
* @brief Gets the mass fraction of the species given the mean molar mass.
* @param meanMolarMass The mean molar mass.
* @return The mass fraction of the species.
*/
double mass_fraction(double meanMolarMass) const;
/**
* @brief Gets the number fraction of the species.
* @return The number fraction of the species.
*/
double number_fraction() const;
/**
* @brief Gets the number fraction of the species given the total moles.
* @param totalMoles The total moles.
* @return The number fraction of the species.
*/
double number_fraction(double totalMoles) const;
/**
* @brief Gets the relative abundance of the species.
* @return The relative abundance of the species.
*/
double rel_abundance() const;
/**
* @brief Gets the isotope of the species.
* @return The isotope of the species.
*/
chemSpecies::Species isotope() const;
/**
* @brief Gets the mode of the composition entry.
* @return True if mass fraction mode, false if number fraction mode.
*/
bool getMassFracMode() const;
/**
* @brief Sets the mass fraction of the species.
* @param mass_fraction The mass fraction to set.
*/
void setMassFraction(double mass_fraction);
/**
* @brief Sets the number fraction of the species.
* @param number_fraction The number fraction to set.
*/
void setNumberFraction(double number_fraction);
/**
* @brief Sets the mode to mass fraction mode.
* @param meanMolarMass The mean molar mass.
* @return True if the mode was successfully set, false otherwise.
*/
bool setMassFracMode(double meanMolarMass);
/**
* @brief Sets the mode to number fraction mode.
* @param totalMoles The total moles.
* @return True if the mode was successfully set, false otherwise.
*/
bool setNumberFracMode(double totalMoles);
/**
* @brief Overloaded output stream operator for CompositionEntry.
* @param os The output stream.
* @param entry The CompositionEntry to output.
* @return The output stream.
*/
friend std::ostream& operator<<(std::ostream& os, const CompositionEntry& entry);
};
/**
* @brief Manages the composition of elements.
* @details The composition is a collection of elements with their respective mass fractions.
* The general purpose of this class is to provide a standardized interface for managing the composition of
* any part of 4DSSE. There are a few rules when using this class.
* - Only species in the atomicSpecies.h database can be used. There are 1000s (All species from AME2020) in there so it should not be a problem.
* - Before a mass fraction can be set with a particular instance of Composition, the symbol must be registered. (i.e. register He-3 before setting its mass fraction)
* - Before any composition information can be retrived (e.g. getComposition), the composition must be finalized (call to .finalize()). This checks if the total mass fraction sums to approximatly 1 (within 1 part in 10^8)
* - Any changes made to the composition after finalization will "unfinalize" the composition. This means that the composition must be finalized again before any information can be retrived.
* - The mass fraction of any individual species must be no more than 1 and no less than 0.
* - The only exception to the finalize rule is if the compositon was constructed with symbols and mass fractions at instantiation time. In this case, the composition is automatically finalized.
* however, this means that the composition passed to the constructor must be valid.
*
* *Example Usage:* Constructing a finalized composition with symbols and mass fractions:
* @code
* std::vector<std::string> symbols = {"H", "He"};
* std::vector<double> mass_fractions = {0.7, 0.3};
* Composition comp(symbols, mass_fractions);
* @endcode
* *Example Usage:* Constructing a composition with symbols and finalizing it later:
* @code
* std::vector<std::string> symbols = {"H", "He"};
* Composition comp(symbols);
* comp.setComposition("H", 0.7);
* comp.setComposition("He", 0.3);
* comp.finalize();
* @endcode
*/
class Composition {
private:
serif::config::Config& m_config = serif::config::Config::getInstance();
serif::probe::LogManager& m_logManager = serif::probe::LogManager::getInstance();
quill::Logger* m_logger = m_logManager.getLogger("log");
bool m_finalized = false; ///< True if the composition is finalized.
double m_specificNumberDensity = 0.0; ///< The specific number density of the composition (\sum_{i} X_i m_i. Where X_i is the number fraction of the ith species and m_i is the mass of the ith species).
double m_meanParticleMass = 0.0; ///< The mean particle mass of the composition (\sum_{i} \frac{n_i}{m_i}. where n_i is the number fraction of the ith species and m_i is the mass of the ith species).
bool m_massFracMode = true; ///< True if mass fraction mode, false if number fraction mode.
std::set<std::string> m_registeredSymbols; ///< The registered symbols.
std::unordered_map<std::string, CompositionEntry> m_compositions; ///< The compositions.
/**
* @brief Checks if the given symbol is valid.
* @details A symbol is valid if it is in the atomic species database (species in atomicSpecies.h). These include all the isotopes from AME2020.
* @param symbol The symbol to check.
* @return True if the symbol is valid, false otherwise.
*/
static bool isValidSymbol(const std::string& symbol);
/**
* @brief Checks if the given mass fractions are valid.
* @param mass_fractions The mass fractions to check.
* @return True if the mass fractions are valid, false otherwise.
*/
bool isValidComposition(const std::vector<double>& fractions) const;
/**
* @brief Validates the given mass fractions.
* @param mass_fractions The mass fractions to validate.
* @throws std::invalid_argument if the mass fractions are invalid.
*/
void validateComposition(const std::vector<double>& fractions) const;
/**
* @brief Finalizes the composition in mass fraction mode.
* @param norm If true, the composition will be normalized to sum to 1.
* @return True if the composition is successfully finalized, false otherwise.
*/
bool finalizeMassFracMode(bool norm);
/**
* @brief Finalizes the composition in number fraction mode.
* @param norm If true, the composition will be normalized to sum to 1.
* @return True if the composition is successfully finalized, false otherwise.
*/
bool finalizeNumberFracMode(bool norm);
public:
/**
* @brief Default constructor.
*/
Composition() = default;
/**
* @brief Default destructor.
*/
~Composition() = default;
/**
* @brief Finalizes the composition.
* @param norm If true, the composition will be normalized to sum to 1 [Default False]
* @return True if the composition is successfully finalized, false otherwise.
*/
bool finalize(bool norm=false);
/**
* @brief Constructs a Composition with the given symbols.
* @param symbols The symbols to initialize the composition with.
* *Example Usage:*
* @code
* std::vector<std::string> symbols = {"H", "O"};
* Composition comp(symbols);
* @endcode
*/
explicit Composition(const std::vector<std::string>& symbols);
/**
* @brief Constructs a Composition with the given symbols as a set.
* @param symbols The symbols to initialize the composition with.
* *Example Usage:*
* @code
* std::set<std::string> symbols = {"H", "O"};
* Composition comp(symbols);
* @endcode
*/
explicit Composition(const std::set<std::string>& symbols);
/**
* @brief Constructs a Composition with the given symbols and mass fractions.
* @param symbols The symbols to initialize the composition with.
* @param mass_fractions The mass fractions corresponding to the symbols.
* @param massFracMode True if mass fraction mode, false if number fraction mode.
* *Example Usage:*
* @code
* std::vector<std::string> symbols = {"H", "O"};
* std::vector<double> mass_fractions = {0.1, 0.9};
* Composition comp(symbols, mass_fractions);
* @endcode
*/
Composition(const std::vector<std::string>& symbols, const std::vector<double>& mass_fractions, bool massFracMode=true);
/**
* @brief Constructs a Composition from another Composition.
* @param composition The Composition to copy.
*/
Composition(const Composition& composition);
Composition& operator=(Composition const& other);
/**
* @brief Registers a new symbol.
* @param symbol The symbol to register.
* @param massFracMode True if mass fraction mode, false if number fraction mode.
* *Example Usage:*
* @code
* Composition comp;
* comp.registerSymbol("H");
* @endcode
*/
void registerSymbol(const std::string& symbol, bool massFracMode=true);
/**
* @brief Registers multiple new symbols.
* @param symbols The symbols to register.
* @param massFracMode True if mass fraction mode, false if number fraction mode.
* *Example Usage:*
* @code
* std::vector<std::string> symbols = {"H", "O"};
* Composition comp;
* comp.registerSymbol(symbols);
* @endcode
*/
void registerSymbol(const std::vector<std::string>& symbols, bool massFracMode=true);
/**
* @brief Gets the registered symbols.
* @return A set of registered symbols.
*/
[[nodiscard]] std::set<std::string> getRegisteredSymbols() const;
/**
* @brief Sets the mass fraction for a given symbol.
* @param symbol The symbol to set the mass fraction for.
* @param mass_fraction The mass fraction to set.
* @return The mass fraction that was set.
* *Example Usage:*
* @code
* Composition comp;
* comp.setMassFraction("H", 0.1);
* @endcode
*/
double setMassFraction(const std::string& symbol, const double& mass_fraction);
/**
* @brief Sets the mass fraction for multiple symbols.
* @param symbols The symbols to set the mass fraction for.
* @param mass_fractions The mass fractions corresponding to the symbols.
* @return A vector of mass fractions that were set.
* *Example Usage:*
* @code
* std::vector<std::string> symbols = {"H", "O"};
* std::vector<double> mass_fractions = {0.1, 0.9};
* Composition comp;
* comp.setMassFraction(symbols, mass_fractions);
* @endcode
*/
std::vector<double> setMassFraction(const std::vector<std::string>& symbols, const std::vector<double>& mass_fractions);
/**
* @brief Sets the number fraction for a given symbol.
* @param symbol The symbol to set the number fraction for.
* @param number_fraction The number fraction to set.
* @return The number fraction that was set.
*/
double setNumberFraction(const std::string& symbol, const double& number_fraction);
/**
* @brief Sets the number fraction for multiple symbols.
* @param symbols The symbols to set the number fraction for.
* @param number_fractions The number fractions corresponding to the symbols.
* @return A vector of number fractions that were set.
*/
std::vector<double> setNumberFraction(const std::vector<std::string>& symbols, const std::vector<double>& number_fractions);
/**
* @brief Mix two compositions together with a given fraction.
* @param other The other composition to mix with.
* @param fraction The fraction of the other composition to mix with. This is the fraction of the other composition wrt. to the current. i.e. fraction=1 would mean that 50% of the new composition is from the other and 50% from the current).
*/
Composition mix(const Composition& other, double fraction) const;
/**
* @brief Gets the mass fractions of all compositions.
* @return An unordered map of compositions with their mass fractions.
*/
[[nodiscard]] std::unordered_map<std::string, double> getMassFraction() const;
/**
* @brief Gets the mass fraction for a given symbol.
* @param symbol The symbol to get the mass fraction for.
* @return The mass fraction for the given symbol.
*/
[[nodiscard]] double getMassFraction(const std::string& symbol) const;
/**
* @brief Gets the number fraction for a given symbol.
* @param symbol The symbol to get the number fraction for.
* @return The number fraction for the given symbol.
*/
[[nodiscard]] double getNumberFraction(const std::string& symbol) const;
/**
* @brief Gets the number fractions of all compositions.
* @return An unordered map of compositions with their number fractions.
*/
[[nodiscard]] std::unordered_map<std::string, double> getNumberFraction() const;
/**
* @brief Gets the composition entry and global composition for a given symbol.
* @param symbol The symbol to get the composition for.
* @return A pair containing the CompositionEntry and GlobalComposition for the given symbol.
*/
[[nodiscard]] std::pair<CompositionEntry, GlobalComposition> getComposition(const std::string& symbol) const;
/**
* @brief Gets all composition entries and the global composition.
* @return A pair containing an unordered map of CompositionEntries and the GlobalComposition.
*/
[[nodiscard]] std::pair<std::unordered_map<std::string, CompositionEntry>, GlobalComposition> getComposition() const;
/**
* @brief Compute the mean particle mass of the composition.
* @return Mean particle mass in g.
*/
[[nodiscard]] double getMeanParticleMass() const;
/**
* @brief Compute the mean atomic mass number of the composition.
* @return Mean atomic mass number.
*/
[[nodiscard]] double getMeanAtomicNumber() const;
/**
* @brief Gets a subset of the composition.
* @param symbols The symbols to include in the subset.
* @param method The method to use for the subset (default is "norm").
* @return A Composition object containing the subset.
*/
Composition subset(const std::vector<std::string>& symbols, std::string method="norm") const;
/**
* @brief Check if a symbol is registered.
* @param symbol The symbol to check.
* @return True if the symbol is registered, false otherwise.
*/
bool hasSymbol(const std::string& symbol) const;
/**
* @brief Sets the composition mode.
* @param massFracMode True if mass fraction mode, false if number fraction mode.
*/
void setCompositionMode(bool massFracMode);
/**
* @brief Gets the current canonical composition (X, Y, Z).
* @param harsh If true, this will throw an error if X-Y != Z where Z is computed as the sum of all other elements.
* @return True if mass fraction mode, false if number fraction mode.
*
* @throws std::runtime_error if the composition is not finalized or if the canonical composition cannot be computed.
* @throws std::runtime_error if harsh is true and the canonical composition is not valid.
*/
[[nodiscard]] CanonicalComposition getCanonicalComposition(bool harsh=false) const;
/**
* @brief Overloaded output stream operator for Composition.
* @param os The output stream.
* @param composition The Composition to output.
* @return The output stream.
*/
friend std::ostream& operator<<(std::ostream& os, const Composition& composition);
// Overload the + operator to call mix with a fraction of 0.5
/**
* @brief Overloads the + operator to mix two compositions together with a fraction of 0.5.
* @param other The other composition to mix with.
* @return The mixed composition.
*/
Composition operator+(const Composition& other) const;
};
}; // namespace serif::composition

14
src/config/config.h Normal file
View File

@@ -0,0 +1,14 @@
#pragma once
#include <string>
namespace serif::config {
struct SERiFConfig {
struct EOS {
struct Helm {
std::string logFile = "log";
};
};
};
}

View File

@@ -1,26 +1 @@
# Define the library
config_sources = files(
'private/config.cpp',
)
config_headers = files(
'public/config.h'
)
# Define the libconfig library so it can be linked against by other parts of the build system
libconfig = library('config',
config_sources,
include_directories: include_directories('public'),
cpp_args: ['-fvisibility=default'],
dependencies: [yaml_cpp_dep],
install : true)
config_dep = declare_dependency(
include_directories: include_directories('public'),
link_with: libconfig,
sources: config_sources,
dependencies: [yaml_cpp_dep],
)
# Make headers accessible
install_headers(config_headers, subdir : '4DSSE/config')
include_directories('.', install_dir: )

View File

@@ -1,109 +0,0 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: March 20, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#include <string>
#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
#include <map>
#include "yaml-cpp/yaml.h"
#include "config.h"
namespace serif::config {
Config::Config() {}
Config::~Config() {}
Config& Config::getInstance() {
static Config instance;
return instance;
}
bool Config::loadConfig(const std::string& configFile) {
configFilePath = configFile;
try {
yamlRoot = YAML::LoadFile(configFile);
} catch (YAML::BadFile& e) {
std::cerr << "Error: " << e.what() << std::endl;
return false;
}
m_loaded = true;
return true;
}
bool Config::isKeyInCache(const std::string &key) {
return configMap.find(key) != configMap.end();
}
void Config::addToCache(const std::string &key, const YAML::Node &node) {
configMap[key] = node;
}
void Config::registerUnknownKey(const std::string &key) {
unknownKeys.push_back(key);
}
bool Config::has(const std::string &key) {
if (!m_loaded) {
throw std::runtime_error("Error! Config file not loaded");
}
if (isKeyInCache(key)) { return true; }
YAML::Node node = YAML::Clone(yamlRoot);
std::istringstream keyStream(key);
std::string subKey;
while (std::getline(keyStream, subKey, ':')) {
if (!node[subKey]) {
registerUnknownKey(key);
return false;
}
node = node[subKey]; // go deeper
}
// Key exists and is of the requested type
addToCache(key, node);
return true;
}
void recurse_keys(const YAML::Node& node, std::vector<std::string>& keyList, const std::string& path = "") {
if (node.IsMap()) {
for (const auto& it : node) {
auto key = it.first.as<std::string>();
auto new_path = path.empty() ? key : path + ":" + key;
recurse_keys(it.second, keyList, new_path);
}
} else {
keyList.push_back(path);
}
}
std::vector<std::string> Config::keys() const {
std::vector<std::string> keyList;
YAML::Node node = YAML::Clone(yamlRoot);
recurse_keys(node, keyList);
return keyList;
}
}

View File

@@ -1,245 +0,0 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: March 26, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#pragma once
#include <string>
#include <iostream>
#include <sstream>
#include <vector>
#include <map>
#include <algorithm>
#include <stdexcept>
// Required for YAML parsing
#include "yaml-cpp/yaml.h"
// -- Forward Def of Resource manager to let it act as a friend of Config --
namespace serif::resource { class ResourceManager; }
class configTestPrivateAccessor; // Forward declaration for test utility
namespace serif::config {
/**
* @class Config
* @brief Singleton class to manage configuration settings loaded from a YAML file.
*/
class Config {
private:
/**
* @brief Private constructor to prevent instantiation.
*/
Config();
/**
* @brief Destructor.
*/
~Config();
YAML::Node yamlRoot; ///< Root node of the YAML configuration.
std::string configFilePath; ///< Path to the configuration file.
bool debug = false; ///< Flag to enable debug output.
bool loaded = false; ///< Flag to indicate if the configuration has been loaded.
std::map<std::string, YAML::Node> configMap; ///< Cache for the location of configuration settings.
std::vector<std::string> unknownKeys; ///< Cache for the existence of configuration settings.
/**
* @brief Get a value from the configuration cache.
* @tparam T Type of the value to retrieve.
* @param key Key of the configuration value.
* @param defaultValue Default value to return if the key does not exist.
* @return Configuration value of type T.
*/
template <typename T>
T getFromCache(const std::string &key, T defaultValue) {
if (configMap.find(key) != configMap.end()) {
try {
return configMap[key].as<T>();
} catch (const YAML::Exception& e) {
return defaultValue;
}
}
return defaultValue;
}
/**
* @brief Check if a key exists in the configuration cache.
* @param key Key to check.
* @return True if the key exists in the cache, false otherwise.
*/
bool isKeyInCache(const std::string &key);
/**
* @brief Add a key-value pair to the configuration cache.
* @param key Key of the configuration value.
* @param node YAML node containing the configuration value.
*/
void addToCache(const std::string &key, const YAML::Node &node);
/**
* @brief Register a key as not found in the configuration.
* @param key Key that was not found.
*/
void registerUnknownKey(const std::string &key);
bool m_loaded = false;
// Only friends can access get without a default value
template <typename T>
T get(const std::string &key) {
if (!m_loaded) {
throw std::runtime_error("Error! Config file not loaded");
}
if (has(key)) {
return getFromCache<T>(key, T());
} else {
throw std::runtime_error("Error! Key not found in config file");
}
}
public:
/**
* @brief Get the singleton instance of the Config class.
* @return Reference to the Config instance.
*/
static Config& getInstance();
Config (const Config&) = delete;
Config& operator= (const Config&) = delete;
Config (Config&&) = delete;
Config& operator= (Config&&) = delete;
void setDebug(bool debug) { this->debug = debug; }
/**
* @brief Load configuration from a YAML file.
* @param configFilePath Path to the YAML configuration file.
* @return True if the configuration was loaded successfully, false otherwise.
*/
bool loadConfig(const std::string& configFilePath);
/**
* @brief Get the input table from the configuration.
* @return Input table as a string.
*/
std::string getInputTable() const;
/**
* @brief Get a configuration value by key.
* @tparam T Type of the value to retrieve.
* @param key Key of the configuration value.
* @param defaultValue Default value to return if the key does not exist.
* @return Configuration value of type T.
*
* @example
* @code
* Config& config = Config::getInstance();
* config.loadConfig("example.yaml");
* int maxIter = config.get<int>("opac:lowTemp:numeric:maxIter", 10);
*/
template <typename T>
T get(const std::string &key, T defaultValue) {
if (!m_loaded) {
// ONLY THROW ERROR IF HARSH OR WARN CONFIGURATION
#if defined(CONFIG_HARSH)
throw std::runtime_error("Error! Config file not loaded. To disable this error, recompile with CONFIG_HARSH=0");
#elif defined(CONFIG_WARN)
std::cerr << "Warning! Config file not loaded. This instance of 4DSSE was compiled with CONFIG_WARN so the code will continue using only default values" << std::endl;
#endif
}
// --- Check if the key has already been checked for existence
if (std::find(unknownKeys.begin(), unknownKeys.end(), key) != unknownKeys.end()) {
return defaultValue; // If the key has already been added to the unknown cache do not traverse the YAML tree or hit the cache
}
// --- Check if the key is already in the cache (avoid traversing YAML nodes)
if (isKeyInCache(key)) {
return getFromCache<T>(key, defaultValue);
}
// --- If the key is not in the cache, check the YAML file
else {
YAML::Node node = YAML::Clone(yamlRoot);
std::istringstream keyStream(key);
std::string subKey;
while (std::getline(keyStream, subKey, ':')) {
if (!node[subKey]) {
// Key does not exist
registerUnknownKey(key);
return defaultValue;
}
node = node[subKey]; // go deeper
}
try {
// Key exists and is of the requested type
addToCache(key, node);
return node.as<T>();
} catch (const YAML::Exception& e) {
// Key is not of the requested type
registerUnknownKey(key);
return defaultValue; // return default value if the key does not exist
}
}
}
/**
* @brief Check if the key exists in the given config file
* @param key Key to check;
* @return boolean true or false
*/
bool has(const std::string &key);
/**
* @brief Get all keys defined in the configuration file.
* @return Vector of all keys in the configuration file.
*/
std::vector<std::string> keys() const;
/**
* @brief Print the configuration file path and the YAML root node.
* @param os Output stream.
* @param config Config object to print.
* @return Output stream.
*/
friend std::ostream& operator<<(std::ostream& os, const Config& config) {
if (!config.m_loaded) {
os << "Config file not loaded" << std::endl;
return os;
}
if (!config.debug) {
os << "Config file: " << config.configFilePath << std::endl;
} else{
// Print entire YAML file from root
os << "Config file: " << config.configFilePath << std::endl;
os << config.yamlRoot << std::endl;
}
return os;
}
// Setup gTest class as a friend
friend class ::configTestPrivateAccessor; // Friend declaration for global test accessor
// -- Resource Manager is a friend of config so it can create a seperate instance
friend class serif::resource::ResourceManager; // Adjusted friend declaration
};
}

View File

@@ -1,24 +0,0 @@
# Define the library
const_sources = files(
'private/const.cpp',
)
const_headers = files(
'public/const.h'
)
# Define the libconst library so it can be linked against by other parts of the build system
libconst = library('const',
const_sources,
include_directories: include_directories('public'),
cpp_args: ['-fvisibility=default'],
dependencies: [const_data_dep],
install : true)
const_dep = declare_dependency(
include_directories: include_directories('public'),
link_with: libconst,
)
# Make headers accessible
install_headers(const_headers, subdir : '4DSSE/const')

View File

@@ -1,126 +0,0 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: March 17, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#include <iostream>
#include <map>
#include <vector>
#include <string>
#include <fstream>
#include <sstream>
#include <set>
#include "const.h"
#include "embedded_constants.h" // Generated at build time by meson
namespace serif::constant {
Constants::Constants() {
loaded_ = initialize();
}
bool Constants::initialize() {
return load();
}
Constant Constants::get(const std::string& name) const {
auto it = constants_.find(name);
if (it != constants_.end()) {
return it->second;
} else {
throw std::out_of_range("Constant '" + name + "' not found.");
}
}
Constant Constants::operator[](const std::string& name) const {
return this->get(name);
}
bool Constants::has(const std::string& name) const {
return constants_.find(name) != constants_.end();
}
std::set<std::string> Constants::keys() const {
std::set<std::string> keys;
for (const auto& pair : constants_) {
keys.insert(pair.first);
}
return keys;
}
std::string Constants::trim(const std::string& str) {
size_t first = str.find_first_not_of(" \t");
if (first == std::string::npos) return "";
size_t last = str.find_last_not_of(" \t");
return str.substr(first, last - first + 1);
}
bool Constants::load() {
std::istringstream fileStream(embeddedConstants);
std::string line;
bool data_section = false;
int line_count = 0;
while (std::getline(fileStream, line)) {
line_count++;
// Detect start of data section (double divider line)
if (!data_section) {
if (line.find("Symbol") != std::string::npos) { // Find header row
std::getline(fileStream, line); // Skip dashed divider
std::getline(fileStream, line); // Skip second dashed divider
data_section = true;
}
continue;
}
// Ensure the line is long enough to contain all fields
if (line.length() < 125) continue;
// Define exact column widths from Python script
int start = 0;
const std::string symbol = trim(line.substr(start, col_widths_[0])); start += col_widths_[0];
const std::string name = trim(line.substr(start, col_widths_[1])); start += col_widths_[1];
const std::string valueStr = line.substr(start, col_widths_[2]); start += col_widths_[2];
const std::string unit = trim(line.substr(start, col_widths_[3])); start += col_widths_[3]; // Only trim the unit
const std::string uncertaintyStr = line.substr(start, col_widths_[4]); start += col_widths_[4];
const std::string reference = trim(line.substr(start, col_widths_[5])); // Only trim reference
// Convert numerical fields safely
double value = 0.0, uncertainty = 0.0;
try {
value = std::stod(valueStr);
} catch (...) {
std::cerr << "Warning: Invalid value in line " << line_count << ": " << valueStr << std::endl;
}
try {
uncertainty = std::stod(uncertaintyStr);
} catch (...) {
std::cerr << "Warning: Invalid uncertainty in line " << line_count << ": " << uncertaintyStr << std::endl;
}
// Store in map
constants_.emplace(symbol, Constant{name, value, uncertainty, unit, reference});
}
loaded_ = true;
return true;
}
}

View File

@@ -1,141 +0,0 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: March 17, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#pragma once
#include <string>
#include <iostream>
#include <set>
#include <map>
namespace serif::constant {
/**
* @brief Structure to hold a constant's details.
*/
struct Constant {
const std::string name; ///< Name of the constant
const double value; ///< Value of the constant
const double uncertainty; ///< Uncertainty in the constant's value
const std::string unit; ///< Unit of the constant
const std::string reference; ///< Reference for the constant's value
/**
* @brief Parameterized constructor.
* @param name The name of the constant.
* @param value The value of the constant.
* @param uncertainty The uncertainty in the constant's value.
* @param unit The unit of the constant.
* @param reference The reference for the constant's value.
*/
Constant(const std::string& name, const double value, const double uncertainty, const std::string& unit, const std::string& reference)
: name(name), value(value), uncertainty(uncertainty), unit(unit), reference(reference) {}
/**
* @brief overload the << operator for pretty printing
*/
friend std::ostream& operator<<(std::ostream& os, const Constant& c) {
os << "<" << c.name << ": ";
os << c.value << "±" << c.uncertainty << " ";
os << c.unit << " (" << c.reference << ")>\n";
return os;
}
};
/**
* @brief Class to manage a collection of constants.
*/
class Constants {
private:
bool loaded_ = false; ///< Flag to indicate if constants are loaded
const int col_widths_[6] = {25, 52, 20, 20, 17, 45}; // From the python script used to generate the constants file
std::map<std::string, Constant> constants_; ///< Map to store constants by name
/**
* @brief Default constructor. Private to avoid direct instantiation
*/
Constants();
/**
* @brief Load constants from the embedded header file.
* @return True if loading was successful, false otherwise.
*/
bool load();
/**
* @brief Initialize constants.
* @return True if initialization was successful, false otherwise.
*/
bool initialize();
/**
* @brief Trim leading and trailing whitespace from a string.
* @param str The string to trim.
* @return The trimmed string.
*/
std::string trim(const std::string& str);
public:
/**
* @brief get instance of constants singleton
* @return instance of constants
*/
static Constants& getInstance() {
static Constants instance;
return instance;
}
/**
* @brief Check if constants are loaded.
* @return True if constants are loaded, false otherwise.
*/
bool isLoaded() const { return loaded_; }
/**
* @brief Get a constant by key.
* @param key The name of the constant to retrieve.
* @return The constant associated with the given key.
*/
Constant get(const std::string& key) const;
/**
* @brief Overloaded subscript operator to access constants by key.
* @param key The name of the constant to retrieve.
* @return The constant associated with the given key.
* @throws std::out_of_range if the key is not found.
*/
Constant operator[](const std::string& key) const;
/**
* @brief Check if a constant exists by key.
* @param key The name of the constant to check.
* @return True if the constant exists, false otherwise.
*/
bool has(const std::string& key) const;
/**
* @brief Get a list of all constant keys.
* @return A vector of all constant keys.
*/
std::set<std::string> keys() const;
};
} // namespace serif::const

View File

@@ -3,6 +3,8 @@
#include "helm.h"
#include <string>
#include <ranges>
namespace serif::eos {
EOS::EOS(const EOSio& reader) : m_reader(reader) {}
EOS::EOS(const std::string& filename, const EOSFormat format) : m_reader(EOSio(filename, format)) {}
@@ -15,7 +17,13 @@ namespace serif::eos {
q.rho = in.density; // Density in g/cm^3
q.abar = in.composition.getMeanParticleMass(); // Mean atomic mass in g
q.zbar = in.composition.getMeanAtomicNumber(); // Mean atomic number (dimensionless)
//TODO: Impliment actual get Mean atomic Number Method
q.zbar = 0;
for (const auto& sp : in.composition | std::views::keys) {
q.zbar += sp.z();
}
q.zbar /= static_cast<double>(in.composition.size());
helmholtz::HELMEOSOutput tempOutput;
tempOutput = helmholtz::get_helm_EOS(q, *std::get<std::unique_ptr<helmholtz::HELMTable>>(m_reader.getTable()));
@@ -61,4 +69,4 @@ namespace serif::eos {
const EOSio& EOS::getReader() const {
return m_reader;
}
}
}

View File

@@ -36,10 +36,12 @@
#include "helm.h"
#include "const.h"
#include "probe.h"
#include "config.h"
#include "fourdst/constants/const.h"
#include "fourdst/config/config.h"
#include "quill/LogMacros.h"
#include "quill/Logger.h"
#include "fourdst/logging/logging.h"
using namespace std;
@@ -126,10 +128,8 @@ namespace serif::eos::helmholtz {
// this function reads in the HELM table and stores in the above arrays
std::unique_ptr<HELMTable> read_helm_table(const std::string &filename) {
serif::config::Config& config = serif::config::Config::getInstance();
auto logFile = config.get<std::string>("EOS:Helm:LogFile", "log");
serif::probe::LogManager& logManager = serif::probe::LogManager::getInstance();
quill::Logger* logger = logManager.getLogger(logFile);
fourdst::logging::LogManager& logManager = fourdst::logging::LogManager::getInstance();
quill::Logger* logger = logManager.getLogger("log");
LOG_INFO(logger, "read_helm_table : Reading HELM table from file {}", filename);
// Make a unique pointer to the HELMTable
@@ -237,10 +237,8 @@ namespace serif::eos::helmholtz {
and returns the calculated quantities in the input
***/
serif::eos::helmholtz::HELMEOSOutput get_helm_EOS(serif::eos::helmholtz::HELMEOSInput &q, const serif::eos::helmholtz::HELMTable &table) {
serif::config::Config& config = serif::config::Config::getInstance();
auto logFile = config.get<std::string>("EOS:Helm:LogFile", "log");
serif::probe::LogManager& logManager = serif::probe::LogManager::getInstance();
quill::Logger* logger = logManager.getLogger(logFile);
fourdst::logging::LogManager& logManager = fourdst::logging::LogManager::getInstance();
quill::Logger* logger = logManager.getLogger("log");
fourdst::constant::Constants& constants = fourdst::constant::Constants::getInstance();
const double pi = std::numbers::pi;

View File

@@ -3,7 +3,7 @@
#include "EOSio.h"
#include "helm.h"
#include <string>
#include "composition.h"
#include "fourdst/composition/composition.h"
#include <iomanip>
namespace serif::eos {
@@ -15,7 +15,7 @@ namespace serif::eos {
* required to query the Equation of State.
*/
struct EOSInput {
serif::composition::Composition composition; ///< The composition of the system.
fourdst::composition::Composition composition; ///< The composition of the system.
double density; ///< The density of the system in cgs (g/cm^3).
double temperature; ///< The temperature of the system in cgs (K).
friend std::ostream& operator<<(std::ostream& os, const EOSInput& input) {

View File

@@ -5,22 +5,16 @@
# Utility Libraries
subdir('types')
subdir('misc')
subdir('config')
subdir('probe')
# Physically Informed Libraries
# subdir('constants')
subdir('composition')
# Asset Libraries
subdir('eos')
subdir('meshIO')
# Resouce Manager Libraries
subdir('resource')
# Resouce Manager Libraries Disabled for dev
#subdir('resource')
# Physics Libraries
subdir('network')
subdir('polytrope')
# Python Bindings

View File

@@ -1,38 +0,0 @@
# Define the library
network_sources = files(
'private/network.cpp',
'private/approx8.cpp'
)
network_headers = files(
'public/network.h',
'public/approx8.h'
)
dependencies = [
boost_dep,
const_dep,
quill_dep,
mfem_dep,
config_dep,
probe_dep,
species_weight_dep,
composition_dep,
]
# Define the libnetwork library so it can be linked against by other parts of the build system
libnetwork = library('network',
network_sources,
include_directories: include_directories('public'),
dependencies: dependencies,
install : true)
network_dep = declare_dependency(
include_directories: include_directories('public'),
link_with: libnetwork,
sources: network_sources,
dependencies: dependencies,
)
# Make headers accessible
install_headers(network_headers, subdir : '4DSSE/network')

View File

@@ -1,542 +0,0 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: March 21, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#include <cmath>
#include <stdexcept>
#include <array>
#include <boost/numeric/odeint.hpp>
#include "const.h"
#include "config.h"
#include "quill/LogMacros.h"
#include "approx8.h"
#include "network.h"
/* Nuclear reaction network in cgs units based on Frank Timmes' "approx8".
At this time it does neither screening nor neutrino losses. It includes
the following 8 isotopes:
h1
he3
he4
c12
n14
o16
ne20
mg24
and the following nuclear reactions:
---pp chain---
p(p,e+)d
d(p,g)he3
he3(he3,2p)he4
---CNO cycle---
c12(p,g)n13 - n13 -> c13 + p -> n14
n14(p,g)o15 - o15 + p -> c12 + he4
n14(a,g)f18 - proceeds to ne20
n15(p,a)c12 - / these two n15 reactions are \ CNO I
n15(p,g)o16 - \ used to calculate a fraction / CNO II
o16(p,g)f17 - f17 + e -> o17 and then o17 + p -> n14 + he4
---alpha captures---
c12(a,g)o16
triple alpha
o16(a,g)ne20
ne20(a,g)mg24
c12(c12,a)ne20
c12(o16,a)mg24
At present the rates are all evaluated using a fitting function.
The coefficients to the fit are from reaclib.jinaweb.org .
*/
namespace serif::network::approx8{
// using namespace std;
using namespace boost::numeric::odeint;
//helper functions
// a function to multiply two arrays and then sum the resulting elements: sum(a*b)
double sum_product( const vec7 &a, const vec7 &b){
if (a.size() != b.size()) {
throw std::runtime_error("Error: array size mismatch in sum_product");
}
double sum=0;
for (size_t i=0; i < a.size(); i++) {
sum += a[i] * b[i];
}
return sum;
}
// the fit to the nuclear reaction rates is of the form:
// exp( a0 + a1/T9 + a2/T9^(1/3) + a3*T9^(1/3) + a4*T9 + a5*T9^(5/3) + log(T9) )
// this function returns an array of the T9 terms in that order, where T9 is the temperatures in GigaKelvin
vec7 get_T9_array(const double &T) {
vec7 arr;
const double T9=1e-9*T;
const double T913=pow(T9,1./3.);
arr[0]=1;
arr[1]=1/T9;
arr[2]=1/T913;
arr[3]=T913;
arr[4]=T9;
arr[5]=pow(T9,5./3.);
arr[6]=log(T9);
return arr;
}
// this function uses the two preceding functions to evaluate the rate given the T9 array and the coefficients
double rate_fit(const vec7 &T9, const vec7 &coef){
return exp(sum_product(T9,coef));
}
// p + p -> d; this, like some of the other rates, this is a composite of multiple fits
double pp_rate(const vec7 &T9) {
constexpr vec7 a1 = {-34.78630, 0,-3.511930, 3.100860, -0.1983140, 1.262510e-2, -1.025170};
constexpr vec7 a2 = { -4.364990e+1,-2.460640e-3,-2.750700,-4.248770e-1,1.598700e-2,-6.908750e-4,-2.076250e-1};
return rate_fit(T9,a1) + rate_fit(T9,a2);
}
// p + d -> he3
double dp_rate(const vec7 &T9) {
constexpr vec7 a1 = {7.528980, 0, -3.720800, 0.8717820, 0, 0,-0.6666670};
constexpr vec7 a2 = {8.935250, 0, -3.720800, 0.1986540, 0, 0, 0.3333330};
return rate_fit(T9,a1) + rate_fit(T9,a2);
}
// he3 + he3 -> he4 + 2p
double he3he3_rate(const vec7 &T9){
constexpr vec7 a = {2.477880e+01,0,-12.27700,-0.1036990,-6.499670e-02,1.681910e-02,-6.666670e-01};
return rate_fit(T9,a);
}
// he3(he3,2p)he4
double he3he4_rate(const vec7 &T9){
constexpr vec7 a1 = {1.560990e+01,0.000000e+00,-1.282710e+01,-3.082250e-02,-6.546850e-01,8.963310e-02,-6.666670e-01};
constexpr vec7 a2 = {1.770750e+01,0.000000e+00,-1.282710e+01,-3.812600e+00,9.422850e-02,-3.010180e-03,1.333330e+00};
return rate_fit(T9,a1) + rate_fit(T9,a2);
}
// he4 + he4 + he4 -> c12
double triple_alpha_rate(const vec7 &T9){
constexpr vec7 a1 = {-9.710520e-01,0.000000e+00,-3.706000e+01,2.934930e+01,-1.155070e+02,-1.000000e+01,-1.333330e+00};
constexpr vec7 a2 = {-1.178840e+01,-1.024460e+00,-2.357000e+01,2.048860e+01,-1.298820e+01,-2.000000e+01,-2.166670e+00};
constexpr vec7 a3 = {-2.435050e+01,-4.126560e+00,-1.349000e+01,2.142590e+01,-1.347690e+00,8.798160e-02,-1.316530e+01};
return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3);
}
// c12 + p -> n13
double c12p_rate(const vec7 &T9){
constexpr vec7 a1={1.714820e+01,0.000000e+00,-1.369200e+01,-2.308810e-01,4.443620e+00,-3.158980e+00,-6.666670e-01};
constexpr vec7 a2={1.754280e+01,-3.778490e+00,-5.107350e+00,-2.241110e+00,1.488830e-01,0.000000e+00,-1.500000e+00};
return rate_fit(T9,a1) + rate_fit(T9,a2);
}
// c12 + he4 -> o16
double c12a_rate(const vec7 &T9){
constexpr vec7 a1={6.965260e+01,-1.392540e+00,5.891280e+01,-1.482730e+02,9.083240e+00,-5.410410e-01,7.035540e+01};
constexpr vec7 a2={2.546340e+02,-1.840970e+00,1.034110e+02,-4.205670e+02,6.408740e+01,-1.246240e+01,1.373030e+02};
return rate_fit(T9,a1) + rate_fit(T9,a2);
}
// n14(p,g)o15 - o15 + p -> c12 + he4
double n14p_rate(const vec7 &T9){
constexpr vec7 a1={1.701000e+01,0.000000e+00,-1.519300e+01,-1.619540e-01,-7.521230e+00,-9.875650e-01,-6.666670e-01};
constexpr vec7 a2={2.011690e+01,0.000000e+00,-1.519300e+01,-4.639750e+00,9.734580e+00,-9.550510e+00,3.333330e-01};
constexpr vec7 a3={7.654440e+00,-2.998000e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
constexpr vec7 a4={6.735780e+00,-4.891000e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,6.820000e-02};
return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3) + rate_fit(T9,a4);
}
// n14(a,g)f18 assumed to go on to ne20
double n14a_rate(const vec7 &T9){
constexpr vec7 a1={2.153390e+01,0.000000e+00,-3.625040e+01,0.000000e+00,0.000000e+00,-5.000000e+00,-6.666670e-01};
constexpr vec7 a2={1.968380e-01,-5.160340e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
constexpr vec7 a3={1.389950e+01,-1.096560e+01,-5.622700e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3);
}
// n15(p,a)c12 (CNO I)
double n15pa_rate(const vec7 &T9){
constexpr vec7 a1 = {2.747640e+01,0.000000e+00,-1.525300e+01,1.593180e+00,2.447900e+00,-2.197080e+00,-6.666670e-01};
constexpr vec7 a2 = {-4.873470e+00,-2.021170e+00,0.000000e+00,3.084970e+01,-8.504330e+00,-1.544260e+00,-1.500000e+00};
constexpr vec7 a3 = {2.089720e+01,-7.406000e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
constexpr vec7 a4 = {-6.575220e+00,-1.163800e+00,0.000000e+00,2.271050e+01,-2.907070e+00,2.057540e-01,-1.500000e+00};
return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3) + rate_fit(T9,a4);
}
// n15(p,g)o16 (CNO II)
double n15pg_rate(const vec7 &T9){
constexpr vec7 a1 = {2.001760e+01,0.000000e+00,-1.524000e+01,3.349260e-01,4.590880e+00,-4.784680e+00,-6.666670e-01};
constexpr vec7 a2 = {6.590560e+00,-2.923150e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
constexpr vec7 a3 = {1.454440e+01,-1.022950e+01,0.000000e+00,0.000000e+00,4.590370e-02,0.000000e+00,-1.500000e+00};
return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3);
}
double n15pg_frac(const vec7 &T9){
const double f1=n15pg_rate(T9);
const double f2=n15pa_rate(T9);
return f1/(f1+f2);
}
// o16(p,g)f17 then f17 -> o17(p,a)n14
double o16p_rate(const vec7 &T9){
constexpr vec7 a={1.909040e+01,0.000000e+00,-1.669600e+01,-1.162520e+00,2.677030e-01,-3.384110e-02,-6.666670e-01};
return rate_fit(T9,a);
}
// o16(a,g)ne20
double o16a_rate(const vec7 &T9){
constexpr vec7 a1={2.390300e+01,0.000000e+00,-3.972620e+01,-2.107990e-01,4.428790e-01,-7.977530e-02,-6.666670e-01};
constexpr vec7 a2={3.885710e+00,-1.035850e+01,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
constexpr vec7 a3={9.508480e+00,-1.276430e+01,0.000000e+00,-3.659250e+00,7.142240e-01,-1.075080e-03,-1.500000e+00};
return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3);
}
// ne20(a,g)mg24
double ne20a_rate(const vec7 &T9){
constexpr vec7 a1={2.450580e+01,0.000000e+00,-4.625250e+01,5.589010e+00,7.618430e+00,-3.683000e+00,-6.666670e-01};
constexpr vec7 a2={-3.870550e+01,-2.506050e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
constexpr vec7 a3={1.983070e+00,-9.220260e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00};
constexpr vec7 a4={-8.798270e+00,-1.278090e+01,0.000000e+00,1.692290e+01,-2.573250e+00,2.089970e-01,-1.500000e+00};
return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3) + rate_fit(T9,a4);
}
// c12(c12,a)ne20
double c12c12_rate(const vec7 &T9){
constexpr vec7 a={6.128630e+01,0.000000e+00,-8.416500e+01,-1.566270e+00,-7.360840e-02,-7.279700e-02,-6.666670e-01};
return rate_fit(T9,a);
}
// c12(o16,a)mg24
double c12o16_rate(const vec7 &T9){
constexpr vec7 a={4.853410e+01,3.720400e-01,-1.334130e+02,5.015720e+01,-3.159870e+00,1.782510e-02,-2.370270e+01};
return rate_fit(T9,a);
}
// for Boost ODE solvers either a struct or a class is required
// a Jacobian matrix for implicit solvers
void Jacobian::operator() ( const vector_type &y, matrix_type &J, double /* t */, vector_type &dfdt ) const {
fourdst::constant::Constants& constants = fourdst::constant::Constants::getInstance();
const double avo = constants.get("N_a").value;
const double clight = constants.get("c").value;
// EOS
const vec7 T9=get_T9_array(y[Approx8Net::iTemp]);
// evaluate rates once per call
const double rpp=pp_rate(T9);
const double r33=he3he3_rate(T9);
const double r34=he3he4_rate(T9);
const double r3a=triple_alpha_rate(T9);
const double rc12p=c12p_rate(T9);
const double rc12a=c12a_rate(T9);
const double rn14p=n14p_rate(T9);
const double rn14a=n14a_rate(T9);
const double ro16p=o16p_rate(T9);
const double ro16a=o16a_rate(T9);
const double rne20a=ne20a_rate(T9);
const double r1212=c12c12_rate(T9);
const double r1216=c12o16_rate(T9);
const double pFrac=n15pg_frac(T9);
const double aFrac=1-pFrac;
const double yh1 = y[Approx8Net::ih1];
const double yhe3 = y[Approx8Net::ihe3];
const double yhe4 = y[Approx8Net::ihe4];
const double yc12 = y[Approx8Net::ic12];
const double yn14 = y[Approx8Net::in14];
const double yo16 = y[Approx8Net::io16];
const double yne20 = y[Approx8Net::ine20];
// zero all elements to begin
for (int i=0; i < Approx8Net::nVar; i++) {
for (int j=0; j < Approx8Net::nVar; j++) {
J(i,j)=0.0;
}
}
// h1 jacobian elements
J(Approx8Net::ih1,Approx8Net::ih1) = -3*yh1*rpp - 2*yc12*rc12p -2*yn14*rn14p -2*yo16*ro16p;
J(Approx8Net::ih1,Approx8Net::ihe3) = 2*yhe3*r33 - yhe4*r34;
J(Approx8Net::ih1,Approx8Net::ihe4) = -yhe3*r34;
J(Approx8Net::ih1,Approx8Net::ic12) = -2*yh1*rc12p;
J(Approx8Net::ih1,Approx8Net::in14) = -2*yh1*rn14p;
J(Approx8Net::ih1,Approx8Net::io16) = -2*yh1*ro16p;
// he3 jacobian elements
J(Approx8Net::ihe3,Approx8Net::ih1) = yh1*rpp;
J(Approx8Net::ihe3,Approx8Net::ihe3) = -2*yhe3*r33 - yhe4*r34;
J(Approx8Net::ihe3,Approx8Net::ihe4) = -yhe3*r34;
// he4 jacobian elements
J(Approx8Net::ihe4,Approx8Net::ih1) = yn14*aFrac*rn14p + yo16*ro16p;
J(Approx8Net::ihe4,Approx8Net::ihe3) = yhe3*r33 - yhe4*r34;
J(Approx8Net::ihe4,Approx8Net::ihe4) = yhe3*r34 - 1.5*yhe4*yhe4*r3a - yc12*rc12a - 1.5*yn14*rn14a - yo16*ro16a - yne20*rne20a;
J(Approx8Net::ihe4,Approx8Net::ic12) = -yhe4*rc12a + yc12*r1212 + yo16*r1216;
J(Approx8Net::ihe4,Approx8Net::in14) = yh1*aFrac*rn14p - 1.5*yhe4*rn14a;
J(Approx8Net::ihe4,Approx8Net::io16) = yh1*ro16p - yhe4*ro16a + yc12*r1216;
J(Approx8Net::ihe4,Approx8Net::ine20) = -yhe4*rne20a;
// c12 jacobian elements
J(Approx8Net::ic12,Approx8Net::ih1) = -yc12*rc12p + yn14*aFrac*rn14p;
J(Approx8Net::ic12,Approx8Net::ihe4) = 0.5*yhe4*yhe4*r3a - yhe4*rc12a;
J(Approx8Net::ic12,Approx8Net::ic12) = -yh1*rc12p - yhe4*rc12a - yo16*r1216 - 2*yc12*r1212;
J(Approx8Net::ic12,Approx8Net::in14) = yh1*yn14*aFrac*rn14p;
J(Approx8Net::ic12,Approx8Net::io16) = -yc12*r1216;
// n14 jacobian elements
J(Approx8Net::in14,Approx8Net::ih1) = yc12*rc12p - yn14*rn14p + yo16*ro16p;
J(Approx8Net::in14,Approx8Net::ihe4) = -yn14*rn14a;
J(Approx8Net::in14,Approx8Net::ic12) = yh1*rc12p;
J(Approx8Net::in14,Approx8Net::in14) = -yh1*rn14p - yhe4*rn14a;
J(Approx8Net::in14,Approx8Net::io16) = yo16*ro16p;
// o16 jacobian elements
J(Approx8Net::io16,Approx8Net::ih1) = yn14*pFrac*rn14p - yo16*ro16p;
J(Approx8Net::io16,Approx8Net::ihe4) = yc12*rc12a - yo16*ro16a;
J(Approx8Net::io16,Approx8Net::ic12) = yhe4*rc12a - yo16*r1216;
J(Approx8Net::io16,Approx8Net::in14) = yh1*pFrac*rn14p;
J(Approx8Net::io16,Approx8Net::io16) = yh1*ro16p - yc12*r1216 -yhe4*ro16a;
// ne20 jacobian elements
J(Approx8Net::ine20,Approx8Net::ihe4) = yn14*rn14a + yo16*ro16a - yne20*rne20a;
J(Approx8Net::ine20,Approx8Net::ic12) = yc12*r1212;
J(Approx8Net::ine20,Approx8Net::in14) = yhe4*rn14a;
J(Approx8Net::ine20,Approx8Net::io16) = yo16*ro16a;
J(Approx8Net::ine20,Approx8Net::ine20) = -yhe4*rne20a;
// mg24 jacobian elements
J(Approx8Net::img24,Approx8Net::ihe4) = yne20*rne20a;
J(Approx8Net::img24,Approx8Net::ic12) = yo16*r1216;
J(Approx8Net::img24,Approx8Net::io16) = yc12*r1216;
J(Approx8Net::img24,Approx8Net::ine20) = yhe4*rne20a;
// energy accounting
for (int j=0; j<Approx8Net::nIso; j++) {
for (int i=0; i<Approx8Net::nIso; i++) {
J(Approx8Net::iEnergy,j) += J(i,j)*Approx8Net::mIon[i];
}
J(Approx8Net::iEnergy,j) *= -avo*clight*clight;
}
}
void ODE::operator() ( const vector_type &y, vector_type &dydt, double /* t */) const {
const fourdst::constant::Constants& constants = fourdst::constant::Constants::getInstance();
const double avo = constants.get("N_a").value;
const double clight = constants.get("c").value;
// EOS
const double T = y[Approx8Net::iTemp];
const double den = y[Approx8Net::iDensity];
const vec7 T9=get_T9_array(T);
// rates
const double rpp=den*pp_rate(T9);
const double r33=den*he3he3_rate(T9);
const double r34=den*he3he4_rate(T9);
const double r3a=den*den*triple_alpha_rate(T9);
const double rc12p=den*c12p_rate(T9);
const double rc12a=den*c12a_rate(T9);
const double rn14p=den*n14p_rate(T9);
const double rn14a=n14a_rate(T9);
const double ro16p=den*o16p_rate(T9);
const double ro16a=den*o16a_rate(T9);
const double rne20a=den*ne20a_rate(T9);
const double r1212=den*c12c12_rate(T9);
const double r1216=den*c12o16_rate(T9);
const double pFrac=n15pg_frac(T9);
const double aFrac=1-pFrac;
const double yh1 = y[Approx8Net:: ih1];
const double yhe3 = y[Approx8Net:: ihe3];
const double yhe4 = y[Approx8Net:: ihe4];
const double yc12 = y[Approx8Net:: ic12];
const double yn14 = y[Approx8Net:: in14];
const double yo16 = y[Approx8Net:: io16];
const double yne20 = y[Approx8Net::ine20];
dydt[Approx8Net::ih1] = -1.5*yh1*yh1*rpp;
dydt[Approx8Net::ih1] += yhe3*yhe3*r33;
dydt[Approx8Net::ih1] += -yhe3*yhe4*r34;
dydt[Approx8Net::ih1] += -2*yh1*yc12*rc12p;
dydt[Approx8Net::ih1] += -2*yh1*yn14*rn14p;
dydt[Approx8Net::ih1] += -2*yh1*yo16*ro16p;
dydt[Approx8Net::ihe3] = 0.5*yh1*yh1*rpp;
dydt[Approx8Net::ihe3] += -yhe3*yhe3*r33;
dydt[Approx8Net::ihe3] += -yhe3*yhe4*r34;
dydt[Approx8Net::ihe4] = 0.5*yhe3*yhe3*r33;
dydt[Approx8Net::ihe4] += yhe3*yhe4*r34;
dydt[Approx8Net::ihe4] += -yhe4*yc12*rc12a;
dydt[Approx8Net::ihe4] += yh1*yn14*aFrac*rn14p;
dydt[Approx8Net::ihe4] += yh1*yo16*ro16p;
dydt[Approx8Net::ihe4] += -0.5*yhe4*yhe4*yhe4*r3a;
dydt[Approx8Net::ihe4] += -yhe4*yo16*ro16a;
dydt[Approx8Net::ihe4] += 0.5*yc12*yc12*r1212;
dydt[Approx8Net::ihe4] += yc12*yo16*r1216;
dydt[Approx8Net::ihe4] += -yhe4*yne20*rne20a;
dydt[Approx8Net::ic12] = (1./6.)*yhe4*yhe4*yhe4*r3a;
dydt[Approx8Net::ic12] += -yhe4*yc12*rc12a;
dydt[Approx8Net::ic12] += -yh1*yc12*rc12p;
dydt[Approx8Net::ic12] += yh1*yn14*aFrac*rn14p;
dydt[Approx8Net::ic12] += -yc12*yc12*r1212;
dydt[Approx8Net::ic12] += -yc12*yo16*r1216;
dydt[Approx8Net::in14] = yh1*yc12*rc12p;
dydt[Approx8Net::in14] += -yh1*yn14*rn14p;
dydt[Approx8Net::in14] += yh1*yo16*ro16p;
dydt[Approx8Net::in14] += -yhe4*yn14*rn14a;
dydt[Approx8Net::io16] = yhe4*yc12*rc12a;
dydt[Approx8Net::io16] += yh1*yn14*pFrac*rn14p;
dydt[Approx8Net::io16] += -yh1*yo16*ro16p;
dydt[Approx8Net::io16] += -yc12*yo16*r1216;
dydt[Approx8Net::io16] += -yhe4*yo16*ro16a;
dydt[Approx8Net::ine20] = 0.5*yc12*yc12*r1212;
dydt[Approx8Net::ine20] += yhe4*yn14*rn14a;
dydt[Approx8Net::ine20] += yhe4*yo16*ro16a;
dydt[Approx8Net::ine20] += -yhe4*yne20*rne20a;
dydt[Approx8Net::img24] = yc12*yo16*r1216;
dydt[Approx8Net::img24] += yhe4*yne20*rne20a;
dydt[Approx8Net::iTemp] = 0.;
dydt[Approx8Net::iDensity] = 0.;
// energy accounting
double eNuc = 0.;
for (int i=0; i<Approx8Net::nIso; i++) {
eNuc += Approx8Net::mIon[i]*dydt[i];
}
dydt[Approx8Net::iEnergy] = -eNuc*avo*clight*clight;
}
Approx8Network::Approx8Network() : Network(APPROX8) {}
NetOut Approx8Network::evaluate(const NetIn &netIn) {
m_y = convert_netIn(netIn);
m_tMax = netIn.tMax;
m_dt0 = netIn.dt0;
const double stiff_abs_tol = m_config.get<double>("Network:Approx8:Stiff:AbsTol", 1.0e-6);
const double stiff_rel_tol = m_config.get<double>("Network:Approx8:Stiff:RelTol", 1.0e-6);
const double nonstiff_abs_tol = m_config.get<double>("Network:Approx8:NonStiff:AbsTol", 1.0e-6);
const double nonstiff_rel_tol = m_config.get<double>("Network:Approx8:NonStiff:RelTol", 1.0e-6);
int num_steps = -1;
if (m_stiff) {
LOG_DEBUG(m_logger, "Using stiff solver for Approx8Network");
num_steps = integrate_const(
make_dense_output<rosenbrock4<double>>(stiff_abs_tol, stiff_rel_tol),
std::make_pair(ODE(), Jacobian()),
m_y,
0.0,
m_tMax,
m_dt0
);
} else {
LOG_DEBUG(m_logger, "Using non stiff solver for Approx8Network");
num_steps = integrate_const (
make_dense_output<runge_kutta_dopri5<vector_type>>(nonstiff_abs_tol, nonstiff_rel_tol),
ODE(),
m_y,
0.0,
m_tMax,
m_dt0
);
}
double ySum = 0.0;
for (int i = 0; i < Approx8Net::nIso; i++) {
m_y[i] *= Approx8Net::aIon[i];
ySum += m_y[i];
}
for (int i = 0; i < Approx8Net::nIso; i++) {
m_y[i] /= ySum;
}
NetOut netOut;
std::vector<double> outComposition;
outComposition.reserve(Approx8Net::nVar);
for (int i = 0; i < Approx8Net::nIso; i++) {
outComposition.push_back(m_y[i]);
}
netOut.energy = m_y[Approx8Net::iEnergy];
netOut.num_steps = num_steps;
const std::vector<std::string> symbols = {"H-1", "He-3", "He-4", "C-12", "N-14", "O-16", "Ne-20", "Mg-24"};
netOut.composition = serif::composition::Composition(symbols, outComposition);
return netOut;
}
void Approx8Network::setStiff(bool stiff) {
m_stiff = stiff;
}
vector_type Approx8Network::convert_netIn(const NetIn &netIn) {
vector_type y(Approx8Net::nVar, 0.0);
y[Approx8Net::ih1] = netIn.composition.getMassFraction("H-1");
y[Approx8Net::ihe3] = netIn.composition.getMassFraction("He-3");
y[Approx8Net::ihe4] = netIn.composition.getMassFraction("He-4");
y[Approx8Net::ic12] = netIn.composition.getMassFraction("C-12");
y[Approx8Net::in14] = netIn.composition.getMassFraction("N-14");
y[Approx8Net::io16] = netIn.composition.getMassFraction("O-16");
y[Approx8Net::ine20] = netIn.composition.getMassFraction("Ne-20");
y[Approx8Net::img24] = netIn.composition.getMassFraction("Mg-24");
y[Approx8Net::iTemp] = netIn.temperature;
y[Approx8Net::iDensity] = netIn.density;
y[Approx8Net::iEnergy] = netIn.energy;
double ySum = 0.0;
for (int i = 0; i < Approx8Net::nIso; i++) {
y[i] /= Approx8Net::aIon[i];
ySum += y[i];
}
for (int i = 0; i < Approx8Net::nIso; i++) {
y[i] /= ySum;
}
return y;
}
};
// main program

View File

@@ -1,68 +0,0 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Authors: Aaron Dotter, Emily Boudreaux
// Last Modified: March 21, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#include "network.h"
#include "approx8.h"
#include "probe.h"
#include "quill/LogMacros.h"
namespace serif::network {
Network::Network(const NetworkFormat format) :
m_config(serif::config::Config::getInstance()),
m_logManager(serif::probe::LogManager::getInstance()),
m_logger(m_logManager.getLogger("log")),
m_format(format) {
if (format == NetworkFormat::UNKNOWN) {
LOG_ERROR(m_logger, "nuclearNetwork::Network::Network() called with UNKNOWN format");
throw std::runtime_error("nuclearNetwork::Network::Network() called with UNKNOWN format");
}
}
NetworkFormat Network::getFormat() const {
return m_format;
}
NetworkFormat Network::setFormat(const NetworkFormat format) {
const NetworkFormat oldFormat = m_format;
m_format = format;
return oldFormat;
}
NetOut Network::evaluate(const NetIn &netIn) {
NetOut netOut;
switch (m_format) {
case APPROX8: {
approx8::Approx8Network network;
netOut = network.evaluate(netIn);
break;
}
case UNKNOWN: {
LOG_ERROR(m_logger, "Network format {} is not implemented.", FormatStringLookup.at(m_format));
throw std::runtime_error("Network format not implemented.");
}
default: {
LOG_ERROR(m_logger, "Unknown network format.");
throw std::runtime_error("Unknown network format.");
}
}
return netOut;
}
}

View File

@@ -1,333 +0,0 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: March 21, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#pragma once
#include <array>
#include <boost/numeric/odeint.hpp>
#include "network.h"
/**
* @file approx8.h
* @brief Header file for the Approx8 nuclear reaction network.
*
* This file contains the definitions and declarations for the Approx8 nuclear reaction network.
* The network is based on Frank Timmes' "approx8" and includes 8 isotopes and various nuclear reactions.
* The rates are evaluated using a fitting function with coefficients from reaclib.jinaweb.org.
*/
namespace serif::network::approx8{
/**
* @typedef vector_type
* @brief Alias for a vector of doubles using Boost uBLAS.
*/
typedef boost::numeric::ublas::vector< double > vector_type;
/**
* @typedef matrix_type
* @brief Alias for a matrix of doubles using Boost uBLAS.
*/
typedef boost::numeric::ublas::matrix< double > matrix_type;
/**
* @typedef vec7
* @brief Alias for a std::array of 7 doubles.
*/
typedef std::array<double,7> vec7;
using namespace boost::numeric::odeint;
/**
* @struct Approx8Net
* @brief Contains constants and arrays related to the nuclear network.
*/
struct Approx8Net{
static constexpr int ih1=0;
static constexpr int ihe3=1;
static constexpr int ihe4=2;
static constexpr int ic12=3;
static constexpr int in14=4;
static constexpr int io16=5;
static constexpr int ine20=6;
static constexpr int img24=7;
static constexpr int iTemp=img24+1;
static constexpr int iDensity =iTemp+1;
static constexpr int iEnergy=iDensity+1;
static constexpr int nIso=img24+1; // number of isotopes
static constexpr int nVar=iEnergy+1; // number of variables
static constexpr std::array<int,nIso> aIon = {
1,
3,
4,
12,
14,
16,
20,
24
};
static constexpr std::array<double,nIso> mIon = {
1.67262164e-24,
5.00641157e-24,
6.64465545e-24,
1.99209977e-23,
2.32462686e-23,
2.65528858e-23,
3.31891077e-23,
3.98171594e-23
};
};
/**
* @brief Multiplies two arrays and sums the resulting elements.
* @param a First array.
* @param b Second array.
* @return Sum of the product of the arrays.
* @example
* @code
* vec7 a = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0};
* vec7 b = {0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5};
* double result = sum_product(a, b);
* @endcode
*/
double sum_product( const vec7 &a, const vec7 &b);
/**
* @brief Returns an array of T9 terms for the nuclear reaction rate fit.
* @param T Temperature in GigaKelvin.
* @return Array of T9 terms.
* @example
* @code
* double T = 1.5;
* vec7 T9_array = get_T9_array(T);
* @endcode
*/
vec7 get_T9_array(const double &T);
/**
* @brief Evaluates the nuclear reaction rate given the T9 array and coefficients.
* @param T9 Array of T9 terms.
* @param coef Array of coefficients.
* @return Evaluated rate.
* @example
* @code
* vec7 T9 = get_T9_array(1.5);
* vec7 coef = {1.0, 0.1, 0.01, 0.001, 0.0001, 0.00001, 0.000001};
* double rate = rate_fit(T9, coef);
* @endcode
*/
double rate_fit(const vec7 &T9, const vec7 &coef);
/**
* @brief Calculates the rate for the reaction p + p -> d.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double pp_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction p + d -> he3.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double dp_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction he3 + he3 -> he4 + 2p.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double he3he3_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction he3(he3,2p)he4.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double he3he4_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction he4 + he4 + he4 -> c12.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double triple_alpha_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction c12 + p -> n13.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double c12p_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction c12 + he4 -> o16.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double c12a_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction n14(p,g)o15 - o15 + p -> c12 + he4.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double n14p_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction n14(a,g)f18 assumed to go on to ne20.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double n14a_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction n15(p,a)c12 (CNO I).
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double n15pa_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction n15(p,g)o16 (CNO II).
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double n15pg_rate(const vec7 &T9);
/**
* @brief Calculates the fraction for the reaction n15(p,g)o16.
* @param T9 Array of T9 terms.
* @return Fraction of the reaction.
*/
double n15pg_frac(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction o16(p,g)f17 then f17 -> o17(p,a)n14.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double o16p_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction o16(a,g)ne20.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double o16a_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction ne20(a,g)mg24.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double ne20a_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction c12(c12,a)ne20.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double c12c12_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction c12(o16,a)mg24.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double c12o16_rate(const vec7 &T9);
/**
* @struct Jacobian
* @brief Functor to calculate the Jacobian matrix for implicit solvers.
*/
struct Jacobian {
/**
* @brief Calculates the Jacobian matrix.
* @param y State vector.
* @param J Jacobian matrix.
* @param dfdt Derivative of the state vector.
*/
void operator() ( const vector_type &y, matrix_type &J, double /* t */, vector_type &dfdt ) const;
};
/**
* @struct ODE
* @brief Functor to calculate the derivatives for the ODE solver.
*/
struct ODE {
/**
* @brief Calculates the derivatives of the state vector.
* @param y State vector.
* @param dydt Derivative of the state vector.
*/
void operator() ( const vector_type &y, vector_type &dydt, double /* t */) const;
};
/**
* @class Approx8Network
* @brief Class for the Approx8 nuclear reaction network.
*/
class Approx8Network final : public Network {
public:
Approx8Network();
/**
* @brief Evaluates the nuclear network.
* @param netIn Input parameters for the network.
* @return Output results from the network.
*/
NetOut evaluate(const NetIn &netIn) override;
/**
* @brief Sets whether the solver should use a stiff method.
* @param stiff Boolean indicating if a stiff method should be used.
*/
void setStiff(bool stiff);
/**
* @brief Checks if the solver is using a stiff method.
* @return Boolean indicating if a stiff method is being used.
*/
bool isStiff() const { return m_stiff; }
private:
vector_type m_y;
double m_tMax = 0;
double m_dt0 = 0;
bool m_stiff = false;
/**
* @brief Converts the input parameters to the internal state vector.
* @param netIn Input parameters for the network.
* @return Internal state vector.
*/
static vector_type convert_netIn(const NetIn &netIn);
};
} // namespace nnApprox8

View File

@@ -1,129 +0,0 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Authors: Emily Boudreaux, Aaron Dotter
// Last Modified: March 21, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#pragma once
#include <vector>
#include "probe.h"
#include "config.h"
#include "quill/Logger.h"
#include "composition.h"
#include <unordered_map>
namespace serif::network {
enum NetworkFormat {
APPROX8, ///< Approx8 nuclear reaction network format.
UNKNOWN,
};
static inline std::unordered_map<NetworkFormat, std::string> FormatStringLookup = {
{APPROX8, "Approx8"},
{UNKNOWN, "Unknown"}
};
/**
* @struct NetIn
* @brief Input structure for the network evaluation.
*
* This structure holds the input parameters required for the network evaluation.
*
* Example usage:
* @code
* nuclearNetwork::NetIn netIn;
* netIn.composition = {1.0, 0.0, 0.0};
* netIn.tmax = 1.0e6;
* netIn.dt0 = 1.0e-3;
* netIn.temperature = 1.0e8;
* netIn.density = 1.0e5;
* netIn.energy = 1.0e12;
* @endcode
*/
struct NetIn {
serif::composition::Composition composition; ///< Composition of the network
double tMax; ///< Maximum time
double dt0; ///< Initial time step
double temperature; ///< Temperature in Kelvin
double density; ///< Density in g/cm^3
double energy; ///< Energy in ergs
};
/**
* @struct NetOut
* @brief Output structure for the network evaluation.
*
* This structure holds the output results from the network evaluation.
*
* Example usage:
* @code
* nuclearNetwork::NetOut netOut = network.evaluate(netIn);
* std::vector<double> composition = netOut.composition;
* int steps = netOut.num_steps;
* double energy = netOut.energy;
* @endcode
*/
struct NetOut {
serif::composition::Composition composition; ///< Composition of the network after evaluation
int num_steps; ///< Number of steps taken in the evaluation
double energy; ///< Energy in ergs after evaluation
};
/**
* @class Network
* @brief Class for network evaluation.
*
* This class provides methods to evaluate the network based on the input parameters.
*
* Example usage:
* @code
* nuclearNetwork::Network network;
* nuclearNetwork::NetIn netIn;
* // Set netIn parameters...
* nuclearNetwork::NetOut netOut = network.evaluate(netIn);
* @endcode
*/
class Network {
public:
explicit Network(const NetworkFormat format = NetworkFormat::APPROX8);
virtual ~Network() = default;
NetworkFormat getFormat() const;
NetworkFormat setFormat(const NetworkFormat format);
/**
* @brief Evaluate the network based on the input parameters.
*
* @param netIn Input parameters for the network evaluation.
* @return NetOut Output results from the network evaluation.
*/
virtual NetOut evaluate(const NetIn &netIn);
protected:
serif::config::Config& m_config; ///< Configuration instance
serif::probe::LogManager& m_logManager; ///< Log manager instance
quill::Logger* m_logger; ///< Logger instance
NetworkFormat m_format; ///< Format of the network
};
} // namespace nuclearNetwork

View File

@@ -35,7 +35,7 @@ dependencies = [
probe_dep,
quill_dep,
config_dep,
resourceManager_dep,
# resourceManager_dep,
types_dep,
]

View File

@@ -28,20 +28,24 @@
#include "mfem.hpp"
#include "4DSTARTypes.h"
#include "config.h"
#include "fourdst/config/config.h"
// #include "fourdst/constants/const.h"
#include "integrators.h"
#include "mfem.hpp"
#include "polytropeOperator.h"
#include "polyCoeff.h"
#include "probe.h"
#include "resourceManager.h"
#include "resourceManagerTypes.h"
// #include "resourceManager.h"
// #include "resourceManagerTypes.h"
#include "utilities.h"
#include "quill/LogMacros.h"
#include "fourdst/logging/logging.h"
namespace serif::polytrope {
static mesh::MeshIO loadedMesh("assets/dynamic/mesh/core_midres.msh");
namespace laneEmden {
double a (const int k, const double n) { // NOLINT(*-no-recursion)
@@ -74,8 +78,8 @@ namespace laneEmden {
PolySolver::PolySolver(mfem::Mesh& mesh, const double n, const double order)
: m_config(serif::config::Config::getInstance()), // Updated
m_logManager(serif::probe::LogManager::getInstance()),
: // m_config(fourdst::config::Config::getInstance()), // Updated
m_logManager(fourdst::logging::LogManager::getInstance()),
m_logger(m_logManager.getLogger("log")),
m_polytropicIndex(n),
m_feOrder(order),
@@ -102,11 +106,12 @@ mfem::Mesh& PolySolver::prepareMesh(const double n) {
if (n > 4.99 || n < 0.0) {
throw std::runtime_error("The polytropic index n must be less than 5.0 and greater than 0.0. Currently it is " + std::to_string(n));
}
const serif::resource::ResourceManager& rm = serif::resource::ResourceManager::getInstance();
const serif::resource::types::Resource& genericResource = rm.getResource("mesh:polySphere");
const auto &meshIO = std::get<std::unique_ptr<serif::mesh::MeshIO>>(genericResource);
meshIO->LinearRescale(polycoeff::x1(n)); // Assumes polycoeff is now serif::polytrope::polycoeff
return meshIO->GetMesh();
// Disabled for dev
// const serif::resource::ResourceManager& rm = serif::resource::ResourceManager::getInstance();
// const serif::resource::types::Resource& genericResource = rm.getResource("mesh:polySphere");
return loadedMesh.GetMesh();
}
PolySolver::~PolySolver() = default;
@@ -378,10 +383,10 @@ void PolySolver::setInitialGuess() const {
(*m_phi)(phiCenterDofs[i]) = 0.0;
}
if (m_config.get<bool>("Poly:Solver:ViewInitialGuess", false)) {
serif::probe::glVisView(*m_theta, m_mesh, "θ init");
serif::probe::glVisView(*m_phi, m_mesh, "φ init");
}
// if (m_config.get<bool>("Poly:Solver:ViewInitialGuess", false)) {
// serif::probe::glVisView(*m_theta, m_mesh, "θ init");
// serif::probe::glVisView(*m_phi, m_mesh, "φ init");
// }
}
@@ -390,25 +395,25 @@ void PolySolver::saveAndViewSolution(const mfem::BlockVector& state_vector) cons
mfem::Vector& x_theta = x_block.GetBlock(0);
mfem::Vector& x_phi = x_block.GetBlock(1);
if (m_config.get<bool>("Poly:Output:View", false)) {
serif::probe::glVisView(x_theta, *m_feTheta, "θ Solution");
serif::probe::glVisView(x_phi, *m_fePhi, "ɸ Solution");
}
// if (m_config.get<bool>("Poly:Output:View", false)) {
// serif::probe::glVisView(x_theta, *m_feTheta, "θ Solution");
// serif::probe::glVisView(x_phi, *m_fePhi, "ɸ Solution");
// }
// --- Extract the Solution ---
if (m_config.get<bool>("Poly:Output:1D:Save", true)) {
const auto solutionPath = m_config.get<std::string>("Poly:Output:1D:Path", "polytropeSolution_1D.csv");
auto derivSolPath = "d" + solutionPath;
const auto rayCoLatitude = m_config.get<double>("Poly:Output:1D:RayCoLatitude", 0.0);
const auto rayLongitude = m_config.get<double>("Poly:Output:1D:RayLongitude", 0.0);
const auto raySamples = m_config.get<int>("Poly:Output:1D:RaySamples", 100);
const std::vector rayDirection = {rayCoLatitude, rayLongitude};
serif::probe::getRaySolution(x_theta, *m_feTheta, rayDirection, raySamples, solutionPath);
// Probe::getRaySolution(x_phi, *m_fePhi, rayDirection, raySamples, derivSolPath);
}
// if (m_config.get<bool>("Poly:Output:1D:Save", true)) {
// const auto solutionPath = m_config.get<std::string>("Poly:Output:1D:Path", "polytropeSolution_1D.csv");
// auto derivSolPath = "d" + solutionPath;
//
// const auto rayCoLatitude = m_config.get<double>("Poly:Output:1D:RayCoLatitude", 0.0);
// const auto rayLongitude = m_config.get<double>("Poly:Output:1D:RayLongitude", 0.0);
// const auto raySamples = m_config.get<int>("Poly:Output:1D:RaySamples", 100);
//
// const std::vector rayDirection = {rayCoLatitude, rayLongitude};
//
// serif::probe::getRaySolution(x_theta, *m_feTheta, rayDirection, raySamples, solutionPath);
// // Probe::getRaySolution(x_phi, *m_fePhi, rayDirection, raySamples, derivSolPath);
// }
}
void PolySolver::setOperatorEssentialTrueDofs() const {
@@ -418,15 +423,24 @@ void PolySolver::setOperatorEssentialTrueDofs() const {
void PolySolver::LoadSolverUserParams(double &newtonRelTol, double &newtonAbsTol, int &newtonMaxIter, int &newtonPrintLevel,
double &gmresRelTol, double &gmresAbsTol, int &gmresMaxIter, int &gmresPrintLevel) const {
newtonRelTol = m_config.get<double>("Poly:Solver:Newton:RelTol", 1.e-4);
newtonAbsTol = m_config.get<double>("Poly:Solver:Newton:AbsTol", 1.e-6);
newtonMaxIter = m_config.get<int>("Poly:Solver:Newton:MaxIter", 10);
newtonPrintLevel = m_config.get<int>("Poly:Solver:Newton:PrintLevel", 3);
gmresRelTol = m_config.get<double>("Poly:Solver:GMRES:RelTol", 1.e-12);
gmresAbsTol = m_config.get<double>("Poly:Solver:GMRES:AbsTol", 1.e-12);
gmresMaxIter = m_config.get<int>("Poly:Solver:GMRES:MaxIter", 200);
gmresPrintLevel = m_config.get<int>("Poly:Solver:GMRES:PrintLevel", -1);
// newtonAbsTol = m_config.get<double>("Poly:Solver:Newton:AbsTol", 1.e-6);
// newtonMaxIter = m_config.get<int>("Poly:Solver:Newton:MaxIter", 10);
// newtonPrintLevel = m_config.get<int>("Poly:Solver:Newton:PrintLevel", 3);
//
// gmresRelTol = m_config.get<double>("Poly:Solver:GMRES:RelTol", 1.e-12);
// gmresAbsTol = m_config.get<double>("Poly:Solver:GMRES:AbsTol", 1.e-12);
// gmresMaxIter = m_config.get<int>("Poly:Solver:GMRES:MaxIter", 200);
// gmresPrintLevel = m_config.get<int>("Poly:Solver:GMRES:PrintLevel", -1);
// The config system has changed, for now just hard code these.
newtonAbsTol = 1.0e-6;
newtonMaxIter = 10;
newtonPrintLevel = 3;
gmresRelTol = 1.0e-12;
gmresAbsTol = 1.0e-12;
gmresMaxIter = 200;
gmresPrintLevel = -1;
LOG_DEBUG(m_logger, "Newton Solver (relTol: {:0.2E}, absTol: {:0.2E}, maxIter: {}, printLevel: {})", newtonRelTol, newtonAbsTol, newtonMaxIter, newtonPrintLevel);
LOG_DEBUG(m_logger, "GMRES Solver (relTol: {:0.2E}, absTol: {:0.2E}, maxIter: {}, printLevel: {})", gmresRelTol, gmresAbsTol, gmresMaxIter, gmresPrintLevel);

View File

@@ -27,10 +27,10 @@
#include "integrators.h"
#include "4DSTARTypes.h"
#include "polytropeOperator.h"
#include "config.h"
#include "fourdst/config/config.h"
#include "meshIO.h"
#include "probe.h"
#include "quill/Logger.h"
#include "fourdst/logging/logging.h"
namespace serif {
namespace polytrope {
@@ -282,8 +282,8 @@ public: // Public methods
private: // Private Attributes
// --- Configuration and Logging ---
serif::config::Config& m_config; ///< Reference to the global configuration manager instance.
serif::probe::LogManager& m_logManager; ///< Reference to the global log manager instance.
// fourdst::config::Config& m_config; ///< Reference to the global configuration manager instance.
fourdst::logging::LogManager& m_logManager; ///< Reference to the global log manager instance.
quill::Logger* m_logger; ///< Pointer to the specific logger instance for this class.
// --- Physical and Discretization Parameters ---
@@ -619,4 +619,4 @@ private: // Private methods
};
} // namespace polytrope
} // namespace serif
} // namespace serif

View File

@@ -22,13 +22,13 @@
#include <cmath>
#include "integrators.h"
#include "config.h"
#include "fourdst/config/config.h"
#include <string>
namespace serif::polytrope::polyMFEMUtils {
NonlinearPowerIntegrator::NonlinearPowerIntegrator(const double n) :
m_polytropicIndex(n),
m_epsilon(serif::config::Config::getInstance().get<double>("Poly:Solver:Epsilon", 1.0e-8)) {
m_epsilon(1e-8) {
if (m_polytropicIndex < 0.0) {
throw std::invalid_argument("Polytropic index must be non-negative.");

View File

@@ -26,7 +26,7 @@
#include "mfem_smout.h"
#include <memory>
#include "config.h"
#include "fourdst/config/config.h"
namespace serif {
namespace polytrope {

View File

@@ -22,8 +22,10 @@
#include "mfem.hpp"
#include <string>
#include "config.h"
#include "fourdst/config/config.h"
#include "probe.h"
#include "quill/Logger.h"
#include "fourdst/logging/logging.h"
/**
@@ -69,8 +71,8 @@ namespace serif::polytrope {
*/
virtual void AssembleElementGrad (const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::DenseMatrix &elmat) override;
private:
serif::config::Config& m_config = serif::config::Config::getInstance();
serif::probe::LogManager& m_logManager = serif::probe::LogManager::getInstance();
// fourdst::config::Config& m_config = fourdst::config::Config::getInstance();
fourdst::logging::LogManager& m_logManager = fourdst::logging::LogManager::getInstance();
quill::Logger* m_logger = m_logManager.getLogger("log");
double m_polytropicIndex;
double m_epsilon;

View File

@@ -25,6 +25,8 @@
#include <memory>
#include "probe.h"
#include "fourdst/logging/logging.h"
#include "quill/Logger.h"
namespace serif::polytrope {
@@ -291,7 +293,7 @@ public:
private:
// --- Logging ---
serif::probe::LogManager& m_logManager = serif::probe::LogManager::getInstance(); ///< Reference to the global log manager.
fourdst::logging::LogManager& m_logManager = fourdst::logging::LogManager::getInstance(); ///< Reference to the global log manager.
quill::Logger* m_logger = m_logManager.getLogger("log"); ///< Pointer to the specific logger instance.
// --- Input Bilinear/Nonlinear Forms (owned) ---
@@ -394,4 +396,4 @@ private:
void update_preconditioner(const mfem::Operator &grad) const;
};
} // namespace serif::polytrope
} // namespace serif::polytrope

View File

@@ -31,7 +31,8 @@ dependencies = [
config_dep,
mfem_dep,
quill_dep,
macros_dep
macros_dep,
logging_dep,
]
# Define the liblogger library so it can be linked against by other parts of the build system

View File

@@ -18,14 +18,8 @@
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#include "quill/Backend.h"
#include "quill/Frontend.h"
#include "quill/Logger.h"
#include "quill/sinks/ConsoleSink.h"
#include "quill/sinks/FileSink.h"
#include "quill/LogMacros.h"
#include <stdexcept>
#include <stdexcept>
#include <string>
#include <iostream>
#include <chrono>
@@ -36,11 +30,18 @@
#include "mfem.hpp"
#include "config.h"
#include "fourdst/config/config.h"
#include "probe.h"
#include <thread>
#include "warning_control.h"
#include "quill/LogMacros.h"
#include "quill/Logger.h"
#include "fourdst/logging/logging.h"
namespace serif::probe {
@@ -56,27 +57,29 @@ void wait(int seconds) {
void glVisView(mfem::GridFunction& u, mfem::Mesh& mesh,
const std::string& windowTitle, const std::string& keyset) {
serif::config::Config& config = serif::config::Config::getInstance();
quill::Logger* logger = LogManager::getInstance().getLogger("log");
if (config.get<bool>("Probe:GLVis:Visualization", true)) {
quill::Logger* logger = fourdst::logging::LogManager::getInstance().getLogger("log");
// if (config.get<bool>("Probe:GLVis:Visualization", true)) {
std::string usedKeyset;
LOG_INFO(logger, "Visualizing solution using GLVis...");
LOG_INFO(logger, "Window title: {}", windowTitle);
if (keyset.empty()) {
usedKeyset = config.get<std::string>("Probe:GLVis:DefaultKeyset", "");
// usedKeyset = config.get<std::string>("Probe:GLVis:DefaultKeyset", "");
usedKeyset = "";
} else {
usedKeyset = keyset;
}
LOG_INFO(logger, "Keyset: {}", usedKeyset);
const auto vishost = config.get<std::string>("Probe:GLVis:Host", "localhost");
const auto visport = config.get<int>("Probe:GLVis:Port", 19916);
// const auto vishost = config.get<std::string>("Probe:GLVis:Host", "localhost");
// const auto visport = config.get<int>("Probe:GLVis:Port", 19916);
const std::string vishost = "localhost";
const int visport = 19916;
mfem::socketstream sol_sock(vishost.c_str(), visport);
sol_sock.precision(8);
sol_sock << "solution\n" << mesh << u
<< "window_title '" << windowTitle <<
"'\n" << "keys " << usedKeyset << "\n";
sol_sock << std::flush;
}
// }
}
void glVisView(mfem::Vector &vec, mfem::FiniteElementSpace &fes, const std::string &windowTitle, const std::string& keyset) {
@@ -108,14 +111,15 @@ double getMeshRadius(mfem::Mesh& mesh) {
std::pair<std::vector<double>, std::vector<double>> getRaySolution(mfem::GridFunction& u, mfem::Mesh& mesh,
const std::vector<double>& rayDirection,
int numSamples, std::string filename) {
serif::config::Config& config = serif::config::Config::getInstance();
serif::probe::LogManager& logManager = serif::probe::LogManager::getInstance();
// fourdst::config::Config& config = fourdst::config::Config::getInstance();
fourdst::logging::LogManager& logManager = fourdst::logging::LogManager::getInstance();
quill::Logger* logger = logManager.getLogger("log");
LOG_INFO(logger, "Getting ray solution...");
// Check if the directory to write to exists
// If it does not exist and MakeDir is true create it
// Otherwise throw an exception
bool makeDir = config.get<bool>("Probe:GetRaySolution:MakeDir", true);
// bool makeDir = config.get<bool>("Probe:GetRaySolution:MakeDir", true);
bool makeDir = true;
std::filesystem::path path = filename;
if (makeDir) {
@@ -202,62 +206,4 @@ std::pair<std::vector<double>, std::vector<double>> getRaySolution(mfem::Vector
DEPRECATION_WARNING_ON
return getRaySolution(gf, *fes.GetMesh(), rayDirection, numSamples, filename);
}
LogManager::LogManager() {
serif::config::Config& config = serif::config::Config::getInstance();
quill::Backend::start();
auto CLILogger = quill::Frontend::create_or_get_logger(
"root",
quill::Frontend::create_or_get_sink<quill::ConsoleSink>("sink_id_1"));
newFileLogger(config.get<std::string>("Probe:LogManager:DefaultLogName", "4DSSE.log"), "log");
loggerMap.emplace("stdout", CLILogger);
}
LogManager::~LogManager() = default;
quill::Logger* LogManager::getLogger(const std::string& loggerName) {
auto it = loggerMap.find(loggerName); // Find *once*
if (it == loggerMap.end()) {
throw std::runtime_error("Cannot find logger " + loggerName);
}
return it->second; // Return the raw pointer from the shared_ptr
}
std::vector<std::string> LogManager::getLoggerNames() {
std::vector<std::string> loggerNames;
loggerNames.reserve(loggerMap.size());
for (const auto& pair : loggerMap) { // Use range-based for loop and const auto&
loggerNames.push_back(pair.first);
}
return loggerNames;
}
std::vector<quill::Logger*> LogManager::getLoggers() {
std::vector<quill::Logger*> loggers;
loggers.reserve(loggerMap.size());
for (const auto& pair : loggerMap) {
loggers.push_back(pair.second); // Get the raw pointer
}
return loggers;
}
quill::Logger* LogManager::newFileLogger(const std::string& filename,
const std::string& loggerName) {
auto file_sink = quill::Frontend::create_or_get_sink<quill::FileSink>(
filename,
[]() {
quill::FileSinkConfig cfg;
cfg.set_open_mode('w');
return cfg;
}(),
quill::FileEventNotifier{});
// Get the raw pointer.
quill::Logger* rawLogger = quill::Frontend::create_or_get_logger(loggerName, std::move(file_sink));
// Create a shared_ptr from the raw pointer.
loggerMap.emplace(loggerName, rawLogger);
return rawLogger;
}
} // namespace Probe
} // namespace Probe

View File

@@ -22,15 +22,13 @@
#pragma once
#include <string>
#include <map>
#include <vector>
#include <utility>
#include "mfem.hpp"
#include "quill/Logger.h"
/**
* @brief The Probe namespace contains utility functions for debugging and logging.
* @brief The Probe namespace contains utility functions for debugging.
*/
namespace serif::probe {
/**
@@ -57,7 +55,7 @@ namespace serif::probe {
/**
* @brief Visualize a vector using GLVis.
* @param vec The vector to visualize.
* @param mesh The mesh associated with the vector.
* @param fes The mesh associated with the vector.
* @param windowTitle The title of the visualization window.
* @param keyset The keyset to use for visualization.
*/
@@ -73,65 +71,4 @@ namespace serif::probe {
const std::vector<double>& rayDirection, int numSamples, std::string filename="");
/**
* @brief Class to manage logging operations.
*/
class LogManager {
private:
/**
* @brief Private constructor for singleton pattern.
*/
LogManager();
/**
* @brief Destructor.
*/
~LogManager();
// Map to store pointers to quill loggers (raw pointers as quill deals with its own memory managment in a seperated, detatched, thread)
std::map<std::string, quill::Logger*> loggerMap;
// Prevent copying and assignment (Rule of Zero)
LogManager(const LogManager&) = delete;
LogManager& operator=(const LogManager&) = delete;
public:
/**
* @brief Get the singleton instance of LogManager.
* @return The singleton instance of LogManager.
*/
static LogManager& getInstance() {
static LogManager instance;
return instance;
}
/**
* @brief Get a logger by name.
* @param loggerName The name of the logger.
* @return A pointer to the logger.
*/
quill::Logger* getLogger(const std::string& loggerName);
/**
* @brief Get the names of all loggers.
* @return A vector of logger names.
*/
std::vector<std::string> getLoggerNames();
/**
* @brief Get all loggers.
* @return A vector of pointers to the loggers.
*/
std::vector<quill::Logger*> getLoggers();
/**
* @brief Create a new file logger.
* @param filename The name of the log file.
* @param loggerName The name of the logger.
* @return A pointer to the new logger.
*/
quill::Logger* newFileLogger(const std::string& filename,
const std::string& loggerName);
};
} // namespace Probe

View File

@@ -9,13 +9,15 @@
#include "bindings.h"
#include "species.h"
namespace py = pybind11;
std::string sv_to_string(std::string_view sv) {
return std::string(sv);
}
std::string get_ostream_str(const serif::composition::Composition& comp) {
std::string get_ostream_str(const fourdst::composition::Composition& comp) {
std::ostringstream oss;
oss << comp;
return oss.str();
@@ -24,38 +26,38 @@ std::string get_ostream_str(const serif::composition::Composition& comp) {
void register_comp_bindings(pybind11::module &comp_submodule) {
// --- Bindings for composition and species module ---
py::class_<serif::composition::GlobalComposition>(comp_submodule, "GlobalComposition")
.def_readonly("specificNumberDensity", &serif::composition::GlobalComposition::specificNumberDensity)
.def_readonly("meanParticleMass", &serif::composition::GlobalComposition::meanParticleMass)
py::class_<fourdst::composition::GlobalComposition>(comp_submodule, "GlobalComposition")
.def_readonly("specificNumberDensity", &fourdst::composition::GlobalComposition::specificNumberDensity)
.def_readonly("meanParticleMass", &fourdst::composition::GlobalComposition::meanParticleMass)
.def("__repr__", // Add a string representation for easy printing in Python
[](const serif::composition::GlobalComposition &gc) {
[](const fourdst::composition::GlobalComposition &gc) {
return "<GlobalComposition(specNumDens=" + std::to_string(gc.specificNumberDensity) +
", meanMass=" + std::to_string(gc.meanParticleMass) + ")>";
});
py::class_<serif::composition::CompositionEntry>(comp_submodule, "CompositionEntry")
.def("symbol", &serif::composition::CompositionEntry::symbol)
py::class_<fourdst::composition::CompositionEntry>(comp_submodule, "CompositionEntry")
.def("symbol", &fourdst::composition::CompositionEntry::symbol)
.def("mass_fraction",
py::overload_cast<>(&serif::composition::CompositionEntry::mass_fraction, py::const_),
py::overload_cast<>(&fourdst::composition::CompositionEntry::mass_fraction, py::const_),
"Gets the mass fraction of the species.")
.def("mass_fraction",
py::overload_cast<double>(&serif::composition::CompositionEntry::mass_fraction, py::const_),
py::overload_cast<double>(&fourdst::composition::CompositionEntry::mass_fraction, py::const_),
py::arg("meanMolarMass"), // Name the argument in Python
"Gets the mass fraction of the species given the mean molar mass.")
.def("number_fraction",
py::overload_cast<>(&serif::composition::CompositionEntry::number_fraction, py::const_),
py::overload_cast<>(&fourdst::composition::CompositionEntry::number_fraction, py::const_),
"Gets the number fraction of the species.")
.def("number_fraction",
py::overload_cast<double>(&serif::composition::CompositionEntry::number_fraction, py::const_),
py::overload_cast<double>(&fourdst::composition::CompositionEntry::number_fraction, py::const_),
py::arg("totalMoles"),
"Gets the number fraction of the species given the total moles.")
.def("rel_abundance", &serif::composition::CompositionEntry::rel_abundance)
.def("isotope", &serif::composition::CompositionEntry::isotope) // Assuming Species is bound or convertible
.def("getMassFracMode", &serif::composition::CompositionEntry::getMassFracMode)
.def("rel_abundance", &fourdst::composition::CompositionEntry::rel_abundance)
.def("isotope", &fourdst::composition::CompositionEntry::isotope) // Assuming Species is bound or convertible
.def("getMassFracMode", &fourdst::composition::CompositionEntry::getMassFracMode)
.def("__repr__", // Optional: nice string representation
[](const serif::composition::CompositionEntry &ce) {
[](const fourdst::composition::CompositionEntry &ce) {
// You might want to include more info here now
return "<CompositionEntry(symbol='" + ce.symbol() + "', " +
"mass_frac=" + std::to_string(ce.mass_fraction()) + ", " +
@@ -63,7 +65,7 @@ void register_comp_bindings(pybind11::module &comp_submodule) {
});
// --- Binding for the main Composition class ---
py::class_<serif::composition::Composition>(comp_submodule, "Composition")
py::class_<fourdst::composition::Composition>(comp_submodule, "Composition")
// Constructors
.def(py::init<>(), "Default constructor")
.def(py::init<const std::vector<std::string>&>(),
@@ -75,61 +77,61 @@ void register_comp_bindings(pybind11::module &comp_submodule) {
"Constructor taking symbols, fractions, and mode (True=Mass, False=Number)")
// Methods
.def("finalize", &serif::composition::Composition::finalize, py::arg("norm") = false,
.def("finalize", &fourdst::composition::Composition::finalize, py::arg("norm") = false,
"Finalize the composition, optionally normalizing fractions to sum to 1.")
.def("registerSymbol", py::overload_cast<const std::string&, bool>(&serif::composition::Composition::registerSymbol),
.def("registerSymbol", py::overload_cast<const std::string&, bool>(&fourdst::composition::Composition::registerSymbol),
py::arg("symbol"), py::arg("massFracMode") = true, "Register a single symbol.")
.def("registerSymbol", py::overload_cast<const std::vector<std::string>&, bool>(&serif::composition::Composition::registerSymbol),
.def("registerSymbol", py::overload_cast<const std::vector<std::string>&, bool>(&fourdst::composition::Composition::registerSymbol),
py::arg("symbols"), py::arg("massFracMode") = true, "Register multiple symbols.")
.def("getRegisteredSymbols", &serif::composition::Composition::getRegisteredSymbols,
.def("getRegisteredSymbols", &fourdst::composition::Composition::getRegisteredSymbols,
"Get the set of registered symbols.")
.def("setMassFraction", py::overload_cast<const std::string&, const double&>(&serif::composition::Composition::setMassFraction),
.def("setMassFraction", py::overload_cast<const std::string&, const double&>(&fourdst::composition::Composition::setMassFraction),
py::arg("symbol"), py::arg("mass_fraction"), "Set mass fraction for a single symbol (requires massFracMode). Returns old value.")
.def("setMassFraction", py::overload_cast<const std::vector<std::string>&, const std::vector<double>&>(&serif::composition::Composition::setMassFraction),
.def("setMassFraction", py::overload_cast<const std::vector<std::string>&, const std::vector<double>&>(&fourdst::composition::Composition::setMassFraction),
py::arg("symbols"), py::arg("mass_fractions"), "Set mass fractions for multiple symbols (requires massFracMode). Returns list of old values.")
.def("setNumberFraction", py::overload_cast<const std::string&, const double&>(&serif::composition::Composition::setNumberFraction),
.def("setNumberFraction", py::overload_cast<const std::string&, const double&>(&fourdst::composition::Composition::setNumberFraction),
py::arg("symbol"), py::arg("number_fraction"), "Set number fraction for a single symbol (requires !massFracMode). Returns old value.")
.def("setNumberFraction", py::overload_cast<const std::vector<std::string>&, const std::vector<double>&>(&serif::composition::Composition::setNumberFraction),
.def("setNumberFraction", py::overload_cast<const std::vector<std::string>&, const std::vector<double>&>(&fourdst::composition::Composition::setNumberFraction),
py::arg("symbols"), py::arg("number_fractions"), "Set number fractions for multiple symbols (requires !massFracMode). Returns list of old values.")
.def("mix", &serif::composition::Composition::mix, py::arg("other"), py::arg("fraction"),
.def("mix", &fourdst::composition::Composition::mix, py::arg("other"), py::arg("fraction"),
"Mix with another composition. Returns new Composition.")
.def("getMassFraction", py::overload_cast<const std::string&>(&serif::composition::Composition::getMassFraction, py::const_),
.def("getMassFraction", py::overload_cast<const std::string&>(&fourdst::composition::Composition::getMassFraction, py::const_),
py::arg("symbol"), "Get mass fraction for a symbol (calculates if needed). Requires finalization.")
.def("getMassFraction", py::overload_cast<>(&serif::composition::Composition::getMassFraction, py::const_),
.def("getMassFraction", py::overload_cast<>(&fourdst::composition::Composition::getMassFraction, py::const_),
"Get dictionary of all mass fractions. Requires finalization.")
.def("getNumberFraction", py::overload_cast<const std::string&>(&serif::composition::Composition::getNumberFraction, py::const_),
.def("getNumberFraction", py::overload_cast<const std::string&>(&fourdst::composition::Composition::getNumberFraction, py::const_),
py::arg("symbol"), "Get number fraction for a symbol (calculates if needed). Requires finalization.")
.def("getNumberFraction", py::overload_cast<>(&serif::composition::Composition::getNumberFraction, py::const_),
.def("getNumberFraction", py::overload_cast<>(&fourdst::composition::Composition::getNumberFraction, py::const_),
"Get dictionary of all number fractions. Requires finalization.")
// Note: pybind11 automatically converts std::pair to a Python tuple
.def("getComposition", py::overload_cast<const std::string&>(&serif::composition::Composition::getComposition, py::const_),
.def("getComposition", py::overload_cast<const std::string&>(&fourdst::composition::Composition::getComposition, py::const_),
py::arg("symbol"), "Returns a tuple (CompositionEntry, GlobalComposition) for the symbol. Requires finalization.")
// Binding the version returning map<string, Entry> requires a bit more care or helper function
// to convert the map to a Python dict if needed directly. Let's bind the pair version for now.
.def("getComposition", py::overload_cast<>(&serif::composition::Composition::getComposition, py::const_),
.def("getComposition", py::overload_cast<>(&fourdst::composition::Composition::getComposition, py::const_),
"Returns a tuple (dict[str, CompositionEntry], GlobalComposition) for all symbols. Requires finalization.")
.def("subset", &serif::composition::Composition::subset, py::arg("symbols"), py::arg("method") = "norm",
.def("subset", &fourdst::composition::Composition::subset, py::arg("symbols"), py::arg("method") = "norm",
"Create a new Composition containing only the specified symbols.")
.def("hasSymbol", &serif::composition::Composition::hasSymbol, py::arg("symbol"),
.def("hasSymbol", &fourdst::composition::Composition::hasSymbol, py::arg("symbol"),
"Check if a symbol is registered.")
.def("setCompositionMode", &serif::composition::Composition::setCompositionMode, py::arg("massFracMode"),
.def("setCompositionMode", &fourdst::composition::Composition::setCompositionMode, py::arg("massFracMode"),
"Set the mode (True=Mass, False=Number). Requires finalization before switching.")
// Operator overload
.def(py::self + py::self, "Mix equally with another composition.") // Binds operator+
// Add __repr__ or __str__
.def("__repr__", [](const serif::composition::Composition &comp) {
.def("__repr__", [](const fourdst::composition::Composition &comp) {
return get_ostream_str(comp); // Use helper for C++ operator<<
});
@@ -139,25 +141,25 @@ void register_comp_bindings(pybind11::module &comp_submodule) {
void register_species_bindings(pybind11::module &chem_submodule) {
// --- Bindings for species module ---
py::class_<chemSpecies::Species>(chem_submodule, "Species")
.def("mass", &chemSpecies::Species::mass, "Get atomic mass (amu)")
.def("massUnc", &chemSpecies::Species::massUnc, "Get atomic mass uncertainty (amu)")
.def("bindingEnergy", &chemSpecies::Species::bindingEnergy, "Get binding energy (keV/nucleon?)") // Check units
.def("betaDecayEnergy", &chemSpecies::Species::betaDecayEnergy, "Get beta decay energy (keV?)") // Check units
.def("betaCode", [](const chemSpecies::Species& s){ return sv_to_string(s.betaCode()); }, "Get beta decay code") // Convert string_view
.def("name", [](const chemSpecies::Species& s){ return sv_to_string(s.name()); }, "Get species name (e.g., 'H-1')") // Convert string_view
.def("el", [](const chemSpecies::Species& s){ return sv_to_string(s.el()); }, "Get element symbol (e.g., 'H')") // Convert string_view
.def("nz", &chemSpecies::Species::nz, "Get NZ value")
.def("n", &chemSpecies::Species::n, "Get neutron number N")
.def("z", &chemSpecies::Species::z, "Get proton number Z")
.def("a", &chemSpecies::Species::a, "Get mass number A")
py::class_<fourdst::atomic::Species>(chem_submodule, "Species")
.def("mass", &fourdst::atomic::Species::mass, "Get atomic mass (amu)")
.def("massUnc", &fourdst::atomic::Species::massUnc, "Get atomic mass uncertainty (amu)")
.def("bindingEnergy", &fourdst::atomic::Species::bindingEnergy, "Get binding energy (keV/nucleon?)") // Check units
.def("betaDecayEnergy", &fourdst::atomic::Species::betaDecayEnergy, "Get beta decay energy (keV?)") // Check units
.def("betaCode", [](const fourdst::atomic::Species& s){ return sv_to_string(s.betaCode()); }, "Get beta decay code") // Convert string_view
.def("name", [](const fourdst::atomic::Species& s){ return sv_to_string(s.name()); }, "Get species name (e.g., 'H-1')") // Convert string_view
.def("el", [](const fourdst::atomic::Species& s){ return sv_to_string(s.el()); }, "Get element symbol (e.g., 'H')") // Convert string_view
.def("nz", &fourdst::atomic::Species::nz, "Get NZ value")
.def("n", &fourdst::atomic::Species::n, "Get neutron number N")
.def("z", &fourdst::atomic::Species::z, "Get proton number Z")
.def("a", &fourdst::atomic::Species::a, "Get mass number A")
.def("__repr__",
[](const chemSpecies::Species &s) {
[](const fourdst::atomic::Species &s) {
std::ostringstream oss;
oss << s;
return oss.str();
});
chem_submodule.attr("species") = py::cast(chemSpecies::species); // Expose the species map
chem_submodule.attr("species") = py::cast(fourdst::atomic::species); // Expose the species map
}

View File

@@ -4,7 +4,6 @@ bindings_headers = files('bindings.h')
dependencies = [
composition_dep,
species_weight_dep,
python3_dep,
pybind11_dep,
]

View File

@@ -14,7 +14,7 @@ template <typename T>
void def_config_get(py::module &m) {
m.def("get",
[](const std::string &key, T defaultValue) {
return serif::config::Config::getInstance().get<T>(key, defaultValue);
return fourdst::config::Config::getInstance().get<T>(key, defaultValue);
},
py::arg("key"), py::arg("defaultValue"),
"Get configuration value (type inferred from default)");
@@ -28,28 +28,28 @@ void register_config_bindings(pybind11::module &config_submodule) {
config_submodule.def("loadConfig",
[](const std::string& configFilePath) {
return serif::config::Config::getInstance().loadConfig(configFilePath);
return fourdst::config::Config::getInstance().loadConfig(configFilePath);
},
py::arg("configFilePath"),
"Load configuration from a YAML file.");
config_submodule.def("has",
[](const std::string &key) {
return serif::config::Config::getInstance().has(key);
return fourdst::config::Config::getInstance().has(key);
},
py::arg("key"),
"Check if a key exists in the configuration.");
config_submodule.def("keys",
[]() {
return py::cast(serif::config::Config::getInstance().keys());
return py::cast(fourdst::config::Config::getInstance().keys());
},
"Get a list of all configuration keys.");
config_submodule.def("__repr__",
[]() {
std::ostringstream oss;
oss << serif::config::Config::getInstance(); // Use the existing operator<<
oss << fourdst::config::Config::getInstance(); // Use the existing operator<<
return std::string("<fourdsse_bindings.config module accessing C++ Singleton>\n") + oss.str();
});
}

View File

@@ -5,7 +5,7 @@ bindings_headers = files('bindings.h')
dependencies = [
eos_dep,
config_dep,
resourceManager_dep,
# resourceManager_dep,
python3_dep,
pybind11_dep,
]

View File

@@ -1,8 +1 @@
subdir('composition')
subdir('const')
subdir('config')
subdir('mfem')
subdir('eos')
subdir('polytrope')
subdir('network')
message('Python bindings currently disabled for development')

View File

@@ -10,7 +10,7 @@ bindings_headers = files(
dependencies = [
config_dep,
resourceManager_dep,
# resourceManager_dep,
python3_dep,
pybind11_dep,
mpi_dep,

View File

@@ -13,22 +13,22 @@ namespace py = pybind11;
void register_network_bindings(pybind11::module &network_submodule) {
py::enum_<serif::network::NetworkFormat>(network_submodule, "NetworkFormat")
.value("APPROX8", serif::network::NetworkFormat::APPROX8)
.value("UNKNOWN", serif::network::NetworkFormat::UNKNOWN)
.def("__str__", [](const serif::network::NetworkFormat format) {
return serif::network::FormatStringLookup[format];
py::enum_<gridfire::NetworkFormat>(network_submodule, "NetworkFormat")
.value("APPROX8", gridfire::NetworkFormat::APPROX8)
.value("UNKNOWN", gridfire::NetworkFormat::UNKNOWN)
.def("__str__", [](const gridfire::NetworkFormat format) {
return gridfire::FormatStringLookup[format];
});
py::class_<serif::network::NetIn>(network_submodule, "NetIn")
py::class_<gridfire::NetIn>(network_submodule, "NetIn")
.def(py::init<>())
.def_readwrite("composition", &serif::network::NetIn::composition)
.def_readwrite("tMax", &serif::network::NetIn::tMax)
.def_readwrite("dt0", &serif::network::NetIn::dt0)
.def_readwrite("temperature", &serif::network::NetIn::temperature)
.def_readwrite("density", &serif::network::NetIn::density)
.def_readwrite("energy", &serif::network::NetIn::energy)
.def("__repr__", [](const serif::network::NetIn &netIn) {
.def_readwrite("composition", &gridfire::NetIn::composition)
.def_readwrite("tMax", &gridfire::NetIn::tMax)
.def_readwrite("dt0", &gridfire::NetIn::dt0)
.def_readwrite("temperature", &gridfire::NetIn::temperature)
.def_readwrite("density", &gridfire::NetIn::density)
.def_readwrite("energy", &gridfire::NetIn::energy)
.def("__repr__", [](const gridfire::NetIn &netIn) {
std::stringstream ss;
ss << "NetIn(composition=" << netIn.composition
<< ", tMax=" << netIn.tMax
@@ -39,11 +39,11 @@ void register_network_bindings(pybind11::module &network_submodule) {
return ss.str();
});
py::class_<serif::network::NetOut>(network_submodule, "NetOut")
.def_readonly("composition", &serif::network::NetOut::composition)
.def_readonly("num_steps", &serif::network::NetOut::num_steps)
.def_readonly("energy", &serif::network::NetOut::energy)
.def("__repr__", [](const serif::network::NetOut &netOut) {
py::class_<gridfire::NetOut>(network_submodule, "NetOut")
.def_readonly("composition", &gridfire::NetOut::composition)
.def_readonly("num_steps", &gridfire::NetOut::num_steps)
.def_readonly("energy", &gridfire::NetOut::energy)
.def("__repr__", [](const gridfire::NetOut &netOut) {
std::stringstream ss;
ss << "NetOut(composition=" << netOut.composition
<< ", num_steps=" << netOut.num_steps
@@ -51,14 +51,14 @@ void register_network_bindings(pybind11::module &network_submodule) {
return ss.str();
});
py::class_<serif::network::Network>(network_submodule, "Network")
.def(py::init<serif::network::NetworkFormat>(), py::arg("format"))
.def("evaluate", &serif::network::Network::evaluate, py::arg("netIn"))
.def("getFormat", &serif::network::Network::getFormat)
.def("setFormat", &serif::network::Network::setFormat, py::arg("format"))
.def("__repr__", [](const serif::network::Network &network) {
py::class_<gridfire::Network>(network_submodule, "Network")
.def(py::init<gridfire::NetworkFormat>(), py::arg("format"))
.def("evaluate", &gridfire::Network::evaluate, py::arg("netIn"))
.def("getFormat", &gridfire::Network::getFormat)
.def("setFormat", &gridfire::Network::setFormat, py::arg("format"))
.def("__repr__", [](const gridfire::Network &network) {
std::stringstream ss;
ss << "Network(format=" << serif::network::FormatStringLookup[network.getFormat()] << ")";
ss << "Network(format=" << gridfire::FormatStringLookup[network.getFormat()] << ")";
return ss.str();
});
@@ -68,7 +68,7 @@ void register_network_bindings(pybind11::module &network_submodule) {
}
void register_approx8_bindings(pybind11::module &network_submodule) {
using namespace serif::network::approx8;
using namespace gridfire::approx8;
py::class_<vec7>(network_submodule, "vec7")
.def("__getitem__", [](const vec7 &v, const size_t i) {

View File

@@ -5,7 +5,7 @@ bindings_headers = files('bindings.h')
dependencies = [
python3_dep,
pybind11_dep,
network_dep,
gridfire_dep,
]
shared_module('py_network',

View File

@@ -5,7 +5,7 @@ bindings_headers = files('bindings.h')
dependencies = [
polysolver_dep,
config_dep,
resourceManager_dep,
# resourceManager_dep,
python3_dep,
pybind11_dep,
]

View File

@@ -10,7 +10,6 @@ resourceManager_headers = files(
)
dependencies = [
yaml_cpp_dep,
opatio_dep,
eos_dep,
quill_dep,
@@ -37,4 +36,4 @@ resourceManager_dep = declare_dependency(
)
# Make headers accessible
install_headers(resourceManager_headers, subdir : '4DSSE/resource')
install_headers(resourceManager_headers, subdir : '4DSSE/resource')

View File

@@ -29,7 +29,7 @@
#include "resourceManagerTypes.h"
#include "debug.h"
#include "config.h"
#include "fourdst/config/config.h"
#define STRINGIFY(x) #x
#define TOSTRING(x) STRINGIFY(x)
@@ -37,7 +37,8 @@
namespace serif::resource {
ResourceManager::ResourceManager() {
const std::string defaultDataDir = TOSTRING(DATA_DIR);
m_dataDir = m_config.get<std::string>("Data:Dir", defaultDataDir);
// m_dataDir = m_config.get<std::string>("Data:Dir", defaultDataDir);
m_dataDir = defaultDataDir;
// -- Get the index file path using filesystem to make it a system safe path
const std::string indexFilePath = m_dataDir + "/index.yaml";
// TODO Add checks to make sure data dir exists and index.yaml exists

View File

@@ -22,13 +22,12 @@
#include <vector>
#include <string>
#include <stdexcept>
#include <unordered_map>
#include "resourceManagerTypes.h"
#include "config.h"
#include "probe.h"
#include "quill/LogMacros.h"
#include "fourdst/config/config.h"
#include "quill/Logger.h"
#include "fourdst/logging/logging.h"
/**
* @class ResourceManager
@@ -55,11 +54,11 @@ namespace serif::resource {
*/
ResourceManager& operator=(const ResourceManager&) = delete;
serif::config::Config& m_config = serif::config::Config::getInstance();
serif::probe::LogManager& m_logManager = serif::probe::LogManager::getInstance();
// fourdst::config::Config& m_config = fourdst::config::Config::getInstance();
fourdst::logging::LogManager& m_logManager = fourdst::logging::LogManager::getInstance();
quill::Logger* m_logger = m_logManager.getLogger("log");
serif::config::Config m_resourceConfig;
// fourdst::config::Config m_resourceConfig;
std::string m_dataDir;
std::unordered_map<std::string, types::Resource> m_resources;

View File

@@ -0,0 +1,4 @@
[wrap-git]
url = https://github.com/4D-STAR/GridFire.git
revision = v0.7.4rc2
depth = 1

View File

@@ -0,0 +1,7 @@
[wrap-git]
url = https://github.com/4D-STAR/libcomposition.git
revision = v2.3.0
depth = 1
[provide]
libcomposition = composition_dep

View File

@@ -0,0 +1,7 @@
[wrap-git]
url = https://github.com/4D-STAR/libconfig.git
revision = v2.0.2
depth = 1
[provide]
libconfig = config_dep

View File

@@ -1,6 +1,6 @@
[wrap-git]
url = git@github.com:4D-STAR/libconstants.git
revision = v1.1
url = https://github.com/4D-STAR/libconstants.git
revision = v1.1.1
depth = 1
[provide]

View File

@@ -0,0 +1,7 @@
[wrap-git]
url = https://github.com/4D-STAR/liblogging.git
revision = v1.1.1
depth = 1
[provide]
liblogging = logging_dep

View File

@@ -1,5 +0,0 @@
[wrap-git]
url = https://github.com/odygrd/quill
revision = v8.1.1
[cmake]

View File

@@ -1,5 +0,0 @@
[wrap-git]
url = https://github.com/jbeder/yaml-cpp.git
revision = yaml-cpp-0.7.0
[cmake]

View File

@@ -1,206 +0,0 @@
#include <gtest/gtest.h>
#include <stdexcept>
#include <string>
#include <algorithm>
#include "atomicSpecies.h"
#include "composition.h"
#include "config.h"
std::string EXAMPLE_FILENAME = std::string(getenv("MESON_SOURCE_ROOT")) + "/tests/config/example.yaml";
/**
* @brief Test suite for the composition class.
*/
class compositionTest : public ::testing::Test {};
/**
* @brief Test the constructor of the composition class.
*/
TEST_F(compositionTest, isotopeMasses) {
EXPECT_NO_THROW(chemSpecies::species.at("H-1"));
EXPECT_DOUBLE_EQ(chemSpecies::species.at("H-1").mass(), 1.007825031898);
EXPECT_DOUBLE_EQ(chemSpecies::species.at("He-3").mass(), 3.0160293219700001);
EXPECT_DOUBLE_EQ(chemSpecies::species.at("He-4").mass(),4.0026032541300003);
}
TEST_F(compositionTest, constructor) {
serif::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME);
EXPECT_NO_THROW(serif::composition::Composition comp);
}
TEST_F(compositionTest, registerSymbol) {
serif::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME);
serif::composition::Composition comp;
EXPECT_NO_THROW(comp.registerSymbol("H-1"));
EXPECT_NO_THROW(comp.registerSymbol("He-4"));
EXPECT_THROW(comp.registerSymbol("H-19"), std::runtime_error);
EXPECT_THROW(comp.registerSymbol("He-21"), std::runtime_error);
std::set<std::string> registeredSymbols = comp.getRegisteredSymbols();
EXPECT_TRUE(registeredSymbols.find("H-1") != registeredSymbols.end());
EXPECT_TRUE(registeredSymbols.find("He-4") != registeredSymbols.end());
EXPECT_TRUE(registeredSymbols.find("H-19") == registeredSymbols.end());
EXPECT_TRUE(registeredSymbols.find("He-21") == registeredSymbols.end());
}
TEST_F(compositionTest, setGetComposition) {
serif::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME);
serif::composition::Composition comp;
comp.registerSymbol("H-1");
comp.registerSymbol("He-4");
EXPECT_DOUBLE_EQ(comp.setMassFraction("H-1", 0.5), 0.0);
EXPECT_DOUBLE_EQ(comp.setMassFraction("He-4", 0.5), 0.0);
EXPECT_DOUBLE_EQ(comp.setMassFraction("H-1", 0.6), 0.5);
EXPECT_DOUBLE_EQ(comp.setMassFraction("He-4", 0.4), 0.5);
EXPECT_NO_THROW(comp.finalize());
EXPECT_DOUBLE_EQ(comp.getMassFraction("H-1"), 0.6);
EXPECT_THROW(comp.setMassFraction("He-3", 0.3), std::runtime_error);
EXPECT_NO_THROW(comp.setMassFraction({"H-1", "He-4"}, {0.5, 0.5}));
EXPECT_THROW(auto r = comp.getComposition("H-1"), std::runtime_error);
EXPECT_TRUE(comp.finalize());
EXPECT_DOUBLE_EQ(comp.getComposition("H-1").first.mass_fraction(), 0.5);
EXPECT_NO_THROW(comp.setMassFraction({"H-1", "He-4"}, {0.6, 0.6}));
EXPECT_FALSE(comp.finalize());
EXPECT_THROW(auto r = comp.getComposition("H-1"), std::runtime_error);
}
TEST_F(compositionTest, setGetNumberFraction) {
serif::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME);
serif::composition::Composition comp;
comp.registerSymbol("H-1", false);
comp.registerSymbol("He-4", false);
EXPECT_DOUBLE_EQ(comp.setNumberFraction("H-1", 0.5), 0.0);
EXPECT_DOUBLE_EQ(comp.setNumberFraction("He-4", 0.5), 0.0);
EXPECT_DOUBLE_EQ(comp.setNumberFraction("H-1", 0.6), 0.5);
EXPECT_DOUBLE_EQ(comp.setNumberFraction("He-4", 0.4), 0.5);
EXPECT_NO_THROW(comp.finalize());
EXPECT_DOUBLE_EQ(comp.getNumberFraction("H-1"), 0.6);
EXPECT_THROW(comp.setNumberFraction("He-3", 0.3), std::runtime_error);
}
TEST_F(compositionTest, subset) {
serif::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME);
serif::composition::Composition comp;
comp.registerSymbol("H-1");
comp.registerSymbol("He-4");
comp.setMassFraction("H-1", 0.6);
comp.setMassFraction("He-4", 0.4);
EXPECT_NO_THROW(comp.finalize());
std::vector<std::string> symbols = {"H-1"};
serif::composition::Composition subsetComp = comp.subset(symbols, "norm");
EXPECT_TRUE(subsetComp.finalize());
EXPECT_DOUBLE_EQ(subsetComp.getMassFraction("H-1"), 1.0);
}
TEST_F(compositionTest, finalizeWithNormalization) {
serif::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME);
serif::composition::Composition comp;
comp.registerSymbol("H-1");
comp.registerSymbol("He-4");
comp.setMassFraction("H-1", 0.3);
comp.setMassFraction("He-4", 0.3);
EXPECT_TRUE(comp.finalize(true));
EXPECT_DOUBLE_EQ(comp.getMassFraction("H-1"), 0.5);
EXPECT_DOUBLE_EQ(comp.getMassFraction("He-4"), 0.5);
}
TEST_F(compositionTest, finalizeWithoutNormalization) {
serif::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME);
serif::composition::Composition comp;
comp.registerSymbol("H-1");
comp.registerSymbol("He-4");
comp.setMassFraction("H-1", 0.5);
comp.setMassFraction("He-4", 0.5);
EXPECT_TRUE(comp.finalize(false));
EXPECT_DOUBLE_EQ(comp.getMassFraction("H-1"), 0.5);
EXPECT_DOUBLE_EQ(comp.getMassFraction("He-4"), 0.5);
}
TEST_F(compositionTest, getComposition) {
serif::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME);
serif::composition::Composition comp;
comp.registerSymbol("H-1");
comp.registerSymbol("He-4");
comp.setMassFraction("H-1", 0.6);
comp.setMassFraction("He-4", 0.4);
EXPECT_NO_THROW(comp.finalize());
auto compositionEntry = comp.getComposition("H-1");
EXPECT_DOUBLE_EQ(compositionEntry.first.mass_fraction(), 0.6);
EXPECT_DOUBLE_EQ(compositionEntry.second.meanParticleMass, 1.4382769310381101);
EXPECT_DOUBLE_EQ(compositionEntry.second.specificNumberDensity, 1.0/1.4382769310381101);
}
TEST_F(compositionTest, setCompositionMode) {
serif::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME);
serif::composition::Composition comp;
comp.registerSymbol("H-1");
comp.registerSymbol("He-4");
comp.setMassFraction("H-1", 0.6);
comp.setMassFraction("He-4", 0.4);
EXPECT_NO_THROW(comp.finalize());
EXPECT_DOUBLE_EQ(comp.getMassFraction("H-1"), 0.6);
EXPECT_DOUBLE_EQ(comp.getMassFraction("He-4"), 0.4);
EXPECT_NO_THROW(comp.setCompositionMode(false));
EXPECT_NO_THROW(comp.setNumberFraction("H-1", 0.9));
EXPECT_NO_THROW(comp.setNumberFraction("He-4", 0.1));
EXPECT_THROW(comp.setCompositionMode(true), std::runtime_error);
EXPECT_NO_THROW(comp.finalize());
EXPECT_NO_THROW(comp.setCompositionMode(true));
}
TEST_F(compositionTest, hasSymbol) {
serif::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME);
serif::composition::Composition comp;
comp.registerSymbol("H-1");
comp.registerSymbol("He-4");
comp.setMassFraction("H-1", 0.6);
comp.setMassFraction("He-4", 0.4);
EXPECT_NO_THROW(comp.finalize());
EXPECT_TRUE(comp.hasSymbol("H-1"));
EXPECT_TRUE(comp.hasSymbol("He-4"));
EXPECT_FALSE(comp.hasSymbol("H-2"));
EXPECT_FALSE(comp.hasSymbol("He-3"));
}
TEST_F(compositionTest, mix) {
serif::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME);
serif::composition::Composition comp1;
comp1.registerSymbol("H-1");
comp1.registerSymbol("He-4");
comp1.setMassFraction("H-1", 0.6);
comp1.setMassFraction("He-4", 0.4);
EXPECT_NO_THROW(comp1.finalize());
serif::composition::Composition comp2;
comp2.registerSymbol("H-1");
comp2.registerSymbol("He-4");
comp2.setMassFraction("H-1", 0.4);
comp2.setMassFraction("He-4", 0.6);
EXPECT_NO_THROW(comp2.finalize());
serif::composition::Composition mixedComp = comp1 + comp2;
EXPECT_TRUE(mixedComp.finalize());
EXPECT_DOUBLE_EQ(mixedComp.getMassFraction("H-1"), 0.5);
EXPECT_DOUBLE_EQ(mixedComp.getMassFraction("He-4"), 0.5);
serif::composition::Composition mixedComp2 = comp1.mix(comp2, 0.25);
EXPECT_TRUE(mixedComp2.finalize());
EXPECT_DOUBLE_EQ(mixedComp2.getMassFraction("H-1"), 0.45);
EXPECT_DOUBLE_EQ(mixedComp2.getMassFraction("He-4"), 0.55);
}

View File

@@ -1,23 +0,0 @@
# Test files for const
test_sources = [
'compositionTest.cpp',
]
foreach test_file : test_sources
exe_name = test_file.split('.')[0]
message('Building test: ' + exe_name)
# Create an executable target for each test
test_exe = executable(
exe_name,
test_file,
dependencies: [gtest_dep, species_weight_dep, gtest_main, composition_dep],
install_rpath: '@loader_path/../../src' # Ensure runtime library path resolves correctly
)
# Add the executable as a test
test(
exe_name,
test_exe,
env: ['MESON_SOURCE_ROOT=' + meson.project_source_root(), 'MESON_BUILD_ROOT=' + meson.project_build_root()])
endforeach

View File

@@ -1,91 +0,0 @@
#include "composition.h"
#include "config.h"
#include <string>
int main(int argv, char* argc[]) {
std::string pathToConfigFile;
if (argv == 2) {
pathToConfigFile = argc[1];
} else {
pathToConfigFile = "config.json";
}
serif::config::Config::getInstance().loadConfig(pathToConfigFile);
serif::composition::Composition comp;
std::vector<std::string> symbols = {"H-1", "He-4"};
comp.registerSymbol(symbols);
comp.setMassFraction("H-1", 0.7);
comp.setMassFraction("He-4", 0.3);
comp.finalize();
std::cout << "=== Mass Fraction Mode ===" << std::endl;
std::cout << "\tH-1 Mass Frac: " << comp.getMassFraction("H-1") << std::endl;
std::cout << "\tH-1 Number Frac: " << comp.getNumberFraction("H-1") << std::endl;
std::cout << "\tHe-4 Mass Frac: " << comp.getMassFraction("He-4") << std::endl;
std::cout << "\tHe-4 Number Frac: " << comp.getNumberFraction("He-4") << std::endl;
std::cout << "\tMass Frac Sum: " << comp.getMassFraction("H-1") + comp.getMassFraction("He-4") << std::endl;
std::cout << "\tNumber Frac Sum: " << comp.getNumberFraction("H-1") + comp.getNumberFraction("He-4") << std::endl;
serif::composition::Composition comp2;
comp2.registerSymbol(symbols, false);
comp2.setNumberFraction("H-1", comp.getNumberFraction("H-1"));
comp2.setNumberFraction("He-4", comp.getNumberFraction("He-4"));
comp2.finalize();
std::cout << "=== Number Fraction Mode ===" << std::endl;
std::cout << "\tH-1 Mass Frac: " << comp2.getMassFraction("H-1") << std::endl;
std::cout << "\tH-1 Number Frac: " << comp2.getNumberFraction("H-1") << std::endl;
std::cout << "\tHe-4 Mass Frac: " << comp2.getMassFraction("He-4") << std::endl;
std::cout << "\tHe-4 Number Frac: " << comp2.getNumberFraction("He-4") << std::endl;
std::cout << "\tMass Frac Sum: " << comp2.getMassFraction("H-1") + comp2.getMassFraction("He-4") << std::endl;
std::cout << "\tNumber Frac Sum: " << comp2.getNumberFraction("H-1") + comp2.getNumberFraction("He-4") << std::endl;
std::vector<std::string> symbols3 = {
"H-1",
"He-3",
"He-4",
"Li-6",
"Li-7",
"O-16",
"C-12",
"Fe-56",
"Ne-20",
"N-14",
"Mg-24",
"Si-28",
"S-32",
"Ca-40",
};
std::vector<double> mass_fractions = {
0.71243,
1e-5,
0.27,
1e-11,
1e-9,
0.009,
0.003,
0.0014,
0.0015,
0.001,
0.0006,
0.0006,
0.0004,
0.00006
};
serif::composition::Composition comp3(symbols3, mass_fractions, true);
std::cout << "=== Mass Fraction Mode ===" << std::endl;
double massFracSum = 0.0;
double numberFracSum = 0.0;
for (const auto& symbol : symbols3) {
std::cout << "\t" << symbol << " Mass Frac: " << comp3.getMassFraction(symbol) << std::endl;
std::cout << "\t" << symbol << " Number Frac: " << comp3.getNumberFraction(symbol) << std::endl;
massFracSum += comp3.getMassFraction(symbol);
numberFracSum += comp3.getNumberFraction(symbol);
}
std::cout << "\tMass Frac Sum: " << massFracSum << std::endl;
std::cout << "\tNumber Frac Sum: " << numberFracSum << std::endl;
return 0;
}

View File

@@ -1,2 +0,0 @@
Composition:
Tracked: "H-1"

View File

@@ -1 +0,0 @@
executable('composition_sandbox', 'comp.cpp', dependencies: [composition_dep, config_dep])

View File

@@ -1,111 +0,0 @@
#include <gtest/gtest.h>
#include "config.h"
#include <iostream>
#include <string>
#include <vector>
#include <set>
#include <sstream>
#include <algorithm>
std::string EXAMPLE_FILENAME = std::string(getenv("MESON_SOURCE_ROOT")) + "/tests/config/example.yaml";
/**
* @file configTest.cpp
* @brief Unit tests for the Config class.
*/
class configTestPrivateAccessor {
public:
static bool callIsKeyInCache(serif::config::Config& config, const std::string& key) {
return config.isKeyInCache(key);
}
static int callCacheSize(serif::config::Config& config) {
return config.configMap.size();
}
static void callAddToCache(serif::config::Config& config, const std::string& key, const YAML::Node& node) {
config.addToCache(key, node);
}
static void callRegisterKeyNotFound(serif::config::Config& config, const std::string& key) {
config.registerUnknownKey(key);
}
static bool CheckIfKeyUnknown(serif::config::Config& config, const std::string& key) {
if (std::find(config.unknownKeys.begin(), config.unknownKeys.end(), key) == config.unknownKeys.end()) {
return false;
}
return true;
}
};
/**
* @brief Test suite for the Config class.
*/
class configTest : public ::testing::Test {};
/**
* @brief Test the constructor of the Config class.
*/
TEST_F(configTest, constructor) {
EXPECT_NO_THROW(serif::config::Config::getInstance());
}
TEST_F(configTest, loadConfig) {
serif::config::Config& config = serif::config::Config::getInstance();
EXPECT_TRUE(config.loadConfig(EXAMPLE_FILENAME));
}
TEST_F(configTest, singletonTest) {
serif::config::Config& config1 = serif::config::Config::getInstance();
serif::config::Config& config2 = serif::config::Config::getInstance();
EXPECT_EQ(&config1, &config2);
}
TEST_F(configTest, getTest) {
serif::config::Config& config = serif::config::Config::getInstance();
config.loadConfig(EXAMPLE_FILENAME);
int maxIter = config.get<int>("opac:lowTemp:numeric:maxIter", 10);
EXPECT_EQ(maxIter, 100);
EXPECT_NE(maxIter, 10);
std::string logLevel = config.get<std::string>("logLevel", "DEBUG");
EXPECT_EQ(logLevel, "INFO");
EXPECT_NE(logLevel, "DEBUG");
float polytropicIndex = config.get<float>("poly:physics:index", 2);
EXPECT_EQ(polytropicIndex, 1.5);
EXPECT_NE(polytropicIndex, 2);
float polytropicIndex2 = config.get<float>("poly:physics:index2", 2.0);
EXPECT_EQ(polytropicIndex2, 2.0);
}
TEST_F(configTest, secondSingletonTest) {
serif::config::Config& config = serif::config::Config::getInstance();
EXPECT_EQ(config.get<int>("opac:lowTemp:numeric:maxIter", 10), 100);
}
TEST_F(configTest, isKeyInCacheTest) {
serif::config::Config& config = serif::config::Config::getInstance();
config.loadConfig(EXAMPLE_FILENAME);
EXPECT_TRUE(configTestPrivateAccessor::callIsKeyInCache(config, "opac:lowTemp:numeric:maxIter"));
EXPECT_FALSE(configTestPrivateAccessor::callIsKeyInCache(config, "opac:lowTemp:numeric:maxIter2"));
}
TEST_F(configTest, cacheSize) {
serif::config::Config& config = serif::config::Config::getInstance();
config.loadConfig(EXAMPLE_FILENAME);
EXPECT_EQ(configTestPrivateAccessor::callCacheSize(config), 3);
EXPECT_NE(configTestPrivateAccessor::callCacheSize(config), 4);
config.get<std::string>("outputDir", "DEBUG");
EXPECT_EQ(configTestPrivateAccessor::callCacheSize(config), 4);
}
TEST_F(configTest, unknownKeyTest) {
serif::config::Config& config = serif::config::Config::getInstance();
config.loadConfig(EXAMPLE_FILENAME);
config.get<int>("opac:lowTemp:numeric:random", 10);
EXPECT_FALSE(configTestPrivateAccessor::CheckIfKeyUnknown(config, "opac:lowTemp:numeric:maxIter"));
EXPECT_TRUE(configTestPrivateAccessor::CheckIfKeyUnknown(config, "opac:lowTemp:numeric:random"));
}

View File

@@ -1,32 +0,0 @@
# High level options
logLevel: "INFO"
outputDir: output
# Module options
poly:
numeric:
newtonTol: 1e-6
newtonMaxIter: 100
gmresTol: 1e-6
gmresMaxIter: 100
physics:
index: 1.5
# Module options
opac:
highTemp:
physics:
table: "/path/to/highTempTable.dat"
numeric:
tol: 1e-6
maxIter: 100
lowTemp:
physics:
table: "/path/to/lowTempTable.dat"
numeric:
tol: 1e-6
maxIter: 100
mesh:
structure:
refine: 2

View File

@@ -1,24 +0,0 @@
# Test files for const
test_sources = [
'configTest.cpp',
]
foreach test_file : test_sources
exe_name = test_file.split('.')[0]
message('Building test: ' + exe_name)
# Create an executable target for each test
test_exe = executable(
exe_name,
test_file,
dependencies: [gtest_dep, config_dep, gtest_main],
include_directories: include_directories('../../src/config/public'),
install_rpath: '@loader_path/../../src' # Ensure runtime library path resolves correctly
)
# Add the executable as a test
test(
exe_name,
test_exe,
env: ['MESON_SOURCE_ROOT=' + meson.project_source_root(), 'MESON_BUILD_ROOT=' + meson.project_build_root()])
endforeach

View File

@@ -26,7 +26,7 @@ std::string TEST_CONFIG = std::string(getenv("MESON_SOURCE_ROOT")) + "/tests/tes
*/
TEST_F(eosTest, read_helm_table) {
serif::config::Config::getInstance().loadConfig(TEST_CONFIG);
fourdst::config::Config::getInstance().loadConfig(TEST_CONFIG);
const serif::resource::ResourceManager& rm = serif::resource::ResourceManager::getInstance();
auto& eos = std::get<std::unique_ptr<serif::eos::EOSio>>(rm.getResource("eos:helm"));
const auto& table = eos->getTable();
@@ -87,7 +87,7 @@ TEST_F(eosTest, get_helm_EOS) {
}
TEST_F(eosTest, eos_using_composition) {
serif::composition::Composition composition;
fourdst::composition::Composition composition;
composition.registerSymbol({
"H-1",
"He-4",

View File

@@ -4,15 +4,10 @@ gtest_main = dependency('gtest_main', required: true)
gtest_nomain_dep = dependency('gtest', main: false, required : true)
# Subdirectories for unit and integration tests
subdir('meshIO')
subdir('probe')
subdir('config')
subdir('eos')
subdir('resource')
subdir('network')
subdir('composition')
#subdir('meshIO')
#subdir('probe')
#subdir('eos')
#subdir('resource')
subdir('poly')
# Subdirectories for sandbox tests
subdir('composition_sandbox')

View File

@@ -1,63 +0,0 @@
#include <gtest/gtest.h>
#include <string>
#include "approx8.h"
#include "config.h"
#include "network.h"
#include "composition.h"
#include <vector>
std::string TEST_CONFIG = std::string(getenv("MESON_SOURCE_ROOT")) + "/tests/testsConfig.yaml";
class approx8Test : public ::testing::Test {};
/**
* @brief Test the constructor of the Config class.
*/
TEST_F(approx8Test, constructor) {
serif::config::Config& config = serif::config::Config::getInstance();
config.loadConfig(TEST_CONFIG);
EXPECT_NO_THROW(serif::network::approx8::Approx8Network());
}
TEST_F(approx8Test, setStiff) {
serif::network::approx8::Approx8Network network;
EXPECT_NO_THROW(network.setStiff(true));
EXPECT_TRUE(network.isStiff());
EXPECT_NO_THROW(network.setStiff(false));
EXPECT_FALSE(network.isStiff());
}
TEST_F(approx8Test, evaluate) {
serif::network::approx8::Approx8Network network;
serif::network::NetIn netIn;
std::vector<double> comp = {0.708, 2.94e-5, 0.276, 0.003, 0.0011, 9.62e-3, 1.62e-3, 5.16e-4};
std::vector<std::string> symbols = {"H-1", "He-3", "He-4", "C-12", "N-14", "O-16", "Ne-20", "Mg-24"};
serif::composition::Composition composition;
composition.registerSymbol(symbols, true);
composition.setMassFraction(symbols, comp);
bool isFinalized = composition.finalize(true);
EXPECT_TRUE(isFinalized);
netIn.composition = composition;
netIn.temperature = 1e7;
netIn.density = 1e2;
netIn.energy = 0.0;
netIn.tMax = 3.15e17;
netIn.dt0 = 1e12;
serif::network::NetOut netOut;
EXPECT_NO_THROW(netOut = network.evaluate(netIn));
double energyFraction = netOut.energy / 1.6433051127589775E+18;
double H1MassFraction = netOut.composition.getMassFraction("H-1")/ 0.50166262445895604;
double He4MassFraction = netOut.composition.getMassFraction("He-4") / 0.48172273720971226;
double relError = 1e-6;
EXPECT_NEAR(H1MassFraction, 1.0, relError);
EXPECT_NEAR(He4MassFraction, 1.0, relError);
EXPECT_NEAR(energyFraction, 1.0, relError);
}

View File

@@ -1,25 +0,0 @@
# Test files for network
test_sources = [
'approx8Test.cpp',
]
foreach test_file : test_sources
exe_name = test_file.split('.')[0]
message('Building test: ' + exe_name)
# Create an executable target for each test
test_exe = executable(
exe_name,
test_file,
dependencies: [gtest_dep, gtest_main, network_dep, species_weight_dep, composition_dep],
include_directories: include_directories('../../src/network/public'),
link_with: libnetwork, # Link the dobj library
install_rpath: '@loader_path/../../src' # Ensure runtime library path resolves correctly
)
# Add the executable as a test
test(
exe_name,
test_exe,
env: ['MESON_SOURCE_ROOT=' + meson.project_source_root(), 'MESON_BUILD_ROOT=' + meson.project_build_root()])
endforeach

View File

@@ -11,7 +11,7 @@ foreach test_file : test_sources
test_exe = executable(
exe_name,
test_file,
dependencies: [gtest_dep, polysolver_dep, probe_dep, quill_dep, config_dep, gtest_main, resourceManager_dep],
dependencies: [gtest_dep, polysolver_dep, probe_dep, quill_dep, config_dep, gtest_main],
install_rpath: '@loader_path/../../src' # Ensure runtime library path resolves correctly
)

View File

@@ -24,8 +24,8 @@
#include "quill/LogMacros.h"
#include "mfem.hpp"
#include "polySolver.h"
#include "probe.h"
#include "config.h"
#include "fourdst/logging/logging.h"
#include "fourdst/config/config.h"
std::string CONFIG_FILENAME = std::string(getenv("MESON_SOURCE_ROOT")) + "/tests/testsConfig.yaml";
@@ -35,20 +35,16 @@ class polyTest : public ::testing::Test {};
TEST_F(polyTest, Solve) {
using namespace serif::polytrope;
serif::config::Config& config = serif::config::Config::getInstance();
config.loadConfig(CONFIG_FILENAME);
serif::probe::LogManager& logManager = serif::probe::LogManager::getInstance();
fourdst::logging::LogManager& logManager = fourdst::logging::LogManager::getInstance();
quill::Logger* logger = logManager.getLogger("log");
LOG_INFO(logger, "Starting polytrope solve test 1...");
config.loadConfig(CONFIG_FILENAME);
double polytropicIndex = config.get<double>("Tests:Poly:Index", 1);
LOG_INFO(logger, "Solving polytrope with n = {:0.2f}", polytropicIndex);
// LOG_INFO(logger, "Solving polytrope with n = {:0.2f}", 1);
PolySolver polytrope(polytropicIndex, 1);
PolySolver polytrope(1, 1);
LOG_INFO(logger, "Solving polytrope...");
EXPECT_NO_THROW(polytrope.solve());
LOG_INFO(logger, "Polytrope solved.");

View File

@@ -13,38 +13,9 @@
std::string TEST_CONFIG = std::string(getenv("MESON_SOURCE_ROOT")) + "/tests/testsConfig.yaml";
std::string getLastLine(const std::string& filename) {
std::ifstream file(filename);
std::string line, lastLine;
if (!file.is_open()) {
throw std::runtime_error("Could not open file");
}
while (std::getline(file, line)) {
lastLine = line;
}
return lastLine; // Returns the last non-empty line
}
std::string stripTimestamps(const std::string& logLine) {
std::regex logPattern(R"(\d+:\d+:\d+\.\d+\s+\[\d+\]\s+probeTest\.cpp:\d+\s+LOG_INFO\s+[A-Za-z]*\s+(.*))");
std::smatch match;
if (std::regex_match(logLine, match, logPattern) && match.size() > 1) {
return match[1].str(); // Extract log message after timestamp
}
return logLine; // Return as-is if pattern doesn't match
}
class probeTest : public ::testing::Test {};
TEST_F(probeTest, DefaultConstructorTest) {
serif::config::Config::getInstance().loadConfig(TEST_CONFIG);
EXPECT_NO_THROW(serif::probe::LogManager::getInstance());
}
TEST_F(probeTest, waitTest) {
auto start = std::chrono::high_resolution_clock::now();
serif::probe::wait(1);
@@ -53,45 +24,3 @@ TEST_F(probeTest, waitTest) {
EXPECT_LE(elapsed.count(), 1.1);
}
TEST_F(probeTest, getLoggerTest) {
serif::config::Config::getInstance().loadConfig(TEST_CONFIG);
serif::probe::LogManager& logManager = serif::probe::LogManager::getInstance();
const std::string loggerName = "testLog";
const std::string filename = "test.log";
quill::Logger* logger = logManager.newFileLogger(filename, loggerName);
EXPECT_NE(logger, nullptr);
LOG_INFO(logger, "This is a test message");
// Wait for the log to be written by calling getLastLine until it is non empty
std::string lastLine;
while (lastLine.empty()) {
lastLine = getLastLine("test.log");
}
EXPECT_EQ(stripTimestamps(lastLine), "This is a test message");
}
TEST_F(probeTest, newFileLoggerTest) {
serif::config::Config::getInstance().loadConfig(TEST_CONFIG);
serif::probe::LogManager& logManager = serif::probe::LogManager::getInstance();
const std::string loggerName = "newLog";
const std::string filename = "newLog.log";
quill::Logger* logger = logManager.newFileLogger(filename, loggerName);
EXPECT_NE(logger, nullptr);
LOG_INFO(logger, "This is a new test message");
// Wait for the log to be written by calling getLastLine until it is non empty
std::string lastLine;
while (lastLine.empty()) {
lastLine = getLastLine(filename);
}
EXPECT_EQ(stripTimestamps(lastLine), "This is a new test message");
}
TEST_F(probeTest, getLoggerNames) {
serif::config::Config::getInstance().loadConfig(TEST_CONFIG);
serif::probe::LogManager& logManager = serif::probe::LogManager::getInstance();
std::vector<std::string> loggerNames = logManager.getLoggerNames();
EXPECT_EQ(loggerNames.size(), 4);
EXPECT_EQ(loggerNames.at(0), "log");
EXPECT_EQ(loggerNames.at(1), "newLog");
EXPECT_EQ(loggerNames.at(2), "stdout");
EXPECT_EQ(loggerNames.at(3), "testLog");
}

View File

@@ -28,12 +28,12 @@ class resourceManagerTest : public ::testing::Test {};
* @brief Test the constructor of the resourceManager class.
*/
TEST_F(resourceManagerTest, constructor) {
serif::config::Config::getInstance().loadConfig(TEST_CONFIG);
fourdst::config::Config::getInstance().loadConfig(TEST_CONFIG);
EXPECT_NO_THROW(serif::resource::ResourceManager::getInstance());
}
TEST_F(resourceManagerTest, getAvaliableResources) {
serif::config::Config::getInstance().loadConfig(TEST_CONFIG);
fourdst::config::Config::getInstance().loadConfig(TEST_CONFIG);
serif::resource::ResourceManager& rm = serif::resource::ResourceManager::getInstance();
std::vector<std::string> resources = rm.getAvailableResources();
std::set<std::string> expected = {"eos:helm", "mesh:polySphere"};
@@ -42,7 +42,7 @@ TEST_F(resourceManagerTest, getAvaliableResources) {
}
TEST_F(resourceManagerTest, getResource) {
serif::config::Config::getInstance().loadConfig(TEST_CONFIG);
fourdst::config::Config::getInstance().loadConfig(TEST_CONFIG);
serif::resource::ResourceManager& rm = serif::resource::ResourceManager::getInstance();
std::string name = "eos:helm";
const serif::resource::types::Resource &r = rm.getResource(name);

View File

@@ -9,7 +9,7 @@
int main(int argv, char **argc) {
std::string CONFIG_FILENAME = std::string(getenv("MESON_SOURCE_ROOT")) + "/tests/testsConfig.yaml";
serif::config::Config& config = serif::config::Config::getInstance();
fourdst::config::Config& config = fourdst::config::Config::getInstance();
config.loadConfig(CONFIG_FILENAME);
// Read the mesh from the given mesh file
std::string meshFile = argc[1];

View File

@@ -1 +1 @@
subdir('meshView')
#subdir('meshView')