From bc95f93ca9d16ba85f514c1d65661c8d9f799e19 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 26 Jun 2026 18:51:42 +0000 Subject: [PATCH 1/3] Fix vibrational/rotational energy sampling cutoff for many DOF Particle::evib() and Particle::erot() sampled internal energy from the equilibrium distribution x^a*exp(-x) (x = E/kT, a = dof/2 - 1) using acceptance-rejection with candidate energies drawn uniformly on [0, 10 kT]. The fixed 10 kT cutoff is below the mean of the distribution (mean = a+1) once the number of degrees of freedom is large. For highly polyatomic species such as SF6, this truncates the high-energy tail of the distribution, systematically biasing the sampled vibrational energy low and distorting its shape (violating detailed balance). Scale the candidate range with the number of degrees of freedom using mean + ~9 standard deviations (a + 1 + 9*sqrt(a+1)) so the cutoff always covers the tail while keeping rejection efficiency roughly constant. The rotational path is bounded (rotdof <= 3) so it was latent there, but is fixed for consistency. Fixes #64. Co-Authored-By: stanmoore1 Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01H4ebZSujHcJJhY6bNMnmZ6 --- src/particle.cpp | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/src/particle.cpp b/src/particle.cpp index 1ee07a278..cc41a4cec 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -1036,9 +1036,13 @@ double Particle::erot(int isp, double temp_thermal, RanKnuth *erandom) eng = -log(erandom->uniform()) * update->boltz * temp_thermal; } else { a = 0.5*particle->species[isp].rotdof-1.0; + // sample E/kT from the equilibrium distribution x^a * exp(-x) + // candidate range must cover the high-energy tail: the mode is at + // x = a, the mean is a+1, the std dev is sqrt(a+1); use mean + ~9 std + // devs so the cut-off scales with rotdof rather than a fixed 10 kT + double xmax = a + 1.0 + 9.0*sqrt(a+1.0); while (1) { - // energy cut-off at 10 kT - erm = 10.0*erandom->uniform(); + erm = xmax*erandom->uniform(); b = pow(erm/a,a) * exp(a-erm); if (b > erandom->uniform()) break; } @@ -1077,9 +1081,14 @@ double Particle::evib(int isp, double temp_thermal, RanKnuth *erandom) eng = -log(erandom->uniform()) * update->boltz * temp_thermal; else if (species[isp].vibdof > 2) { a = 0.5*particle->species[isp].vibdof-1.; + // sample E/kT from the equilibrium distribution x^a * exp(-x) + // candidate range must cover the high-energy tail: the mode is at + // x = a, the mean is a+1, the std dev is sqrt(a+1); use mean + ~9 std + // devs so the cut-off scales with vibdof (a fixed 10 kT cut-off is + // below the mean for molecules with many vibrational modes) + double xmax = a + 1.0 + 9.0*sqrt(a+1.0); while (1) { - // energy cut-off at 10 kT - erm = 10.0*erandom->uniform(); + erm = xmax*erandom->uniform(); b = pow(erm/a,a) * exp(a-erm); if (b > erandom->uniform()) break; } From 37ca2bd3034781749474da0e699912ed26de8622 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 26 Jun 2026 21:31:25 +0000 Subject: [PATCH 2/3] Apply evib/erot sampling cutoff fix to KOKKOS package The KOKKOS GPU paths mirror the CPU samplers and had the same fixed 10 kT candidate cutoff in ParticleKokkos::erot/evib and the inlined erot/evib in SurfCollideDiffuseKokkos. Scale the cutoff with the number of degrees of freedom (a + 1 + 9*sqrt(a+1)) to match the CPU fix, so CPU and GPU runs sample the same vibrational/rotational energy distribution. Co-Authored-By: stanmoore1 Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01H4ebZSujHcJJhY6bNMnmZ6 --- src/KOKKOS/particle_kokkos.h | 14 ++++++++++---- src/KOKKOS/surf_collide_diffuse_kokkos.h | 14 ++++++++++---- 2 files changed, 20 insertions(+), 8 deletions(-) diff --git a/src/KOKKOS/particle_kokkos.h b/src/KOKKOS/particle_kokkos.h index f7013f1ff..76af20d13 100644 --- a/src/KOKKOS/particle_kokkos.h +++ b/src/KOKKOS/particle_kokkos.h @@ -262,9 +262,12 @@ double ParticleKokkos::erot(int isp, double temp_thermal, rand_type &erandom) co eng = -log(erandom.drand()) * boltz * temp_thermal; else { a = 0.5*d_species[isp].rotdof-1.0; + // candidate range must cover the tail of x^a*exp(-x) (mode a, mean a+1, + // std dev sqrt(a+1)); scale the cut-off with dof rather than fixing it + // at 10 kT, which is below the mean for large dof + double xmax = a + 1.0 + 9.0*sqrt(a+1.0); while (1) { - // energy cut-off at 10 kT - erm = 10.0*erandom.drand(); + erm = xmax*erandom.drand(); b = pow(erm/a,a) * exp(a-erm); if (b > erandom.drand()) break; } @@ -297,9 +300,12 @@ double ParticleKokkos::evib(int isp, double temp_thermal, rand_type &erandom) co eng = -log(erandom.drand()) * boltz * temp_thermal; else if (d_species[isp].vibdof > 2) { a = 0.5*d_species[isp].vibdof-1.; + // candidate range must cover the tail of x^a*exp(-x) (mode a, mean a+1, + // std dev sqrt(a+1)); scale the cut-off with dof rather than fixing it + // at 10 kT, which is below the mean for large dof + double xmax = a + 1.0 + 9.0*sqrt(a+1.0); while (1) { - // energy cut-off at 10 kT - erm = 10.0*erandom.drand(); + erm = xmax*erandom.drand(); b = pow(erm/a,a) * exp(a-erm); if (b > erandom.drand()) break; } diff --git a/src/KOKKOS/surf_collide_diffuse_kokkos.h b/src/KOKKOS/surf_collide_diffuse_kokkos.h index 25826164a..4af451c0c 100644 --- a/src/KOKKOS/surf_collide_diffuse_kokkos.h +++ b/src/KOKKOS/surf_collide_diffuse_kokkos.h @@ -330,9 +330,12 @@ class SurfCollideDiffuseKokkos : public SurfCollideDiffuse { eng = -log(rand_gen.drand()) * boltz * temp_thermal; } else { a = 0.5*d_species[isp].rotdof-1.0; + // candidate range must cover the tail of x^a*exp(-x) (mode a, mean a+1, + // std dev sqrt(a+1)); scale the cut-off with dof rather than fixing it + // at 10 kT, which is below the mean for large dof + double xmax = a + 1.0 + 9.0*sqrt(a+1.0); while (1) { - // energy cut-off at 10 kT - erm = 10.0*rand_gen.drand(); + erm = xmax*rand_gen.drand(); b = pow(erm/a,a) * exp(a-erm); if (b > rand_gen.drand()) break; } @@ -370,9 +373,12 @@ class SurfCollideDiffuseKokkos : public SurfCollideDiffuse { eng = -log(rand_gen.drand()) * boltz * temp_thermal; else if (d_species[isp].vibdof > 2) { a = 0.5*d_species[isp].vibdof-1.; + // candidate range must cover the tail of x^a*exp(-x) (mode a, mean a+1, + // std dev sqrt(a+1)); scale the cut-off with dof rather than fixing it + // at 10 kT, which is below the mean for large dof + double xmax = a + 1.0 + 9.0*sqrt(a+1.0); while (1) { - // energy cut-off at 10 kT - erm = 10.0*rand_gen.drand(); + erm = xmax*rand_gen.drand(); b = pow(erm/a,a) * exp(a-erm); if (b > rand_gen.drand()) break; } From e5071ed63a2ae146e55480208e0ae3020b8b5366 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 7 Jul 2026 22:25:03 +0000 Subject: [PATCH 3/3] Rebless in.vibrate logs for the energy sampling cutoff fix The vibrational energy sampling cutoff change shifts the RNG stream in the vibrate example (CO2 has vibdof=8, which is sampled via the rejection branch during create_particles), so the mpi_1 and mpi_4 gold logs no longer match. Regenerate both reference logs with the fixed code. Verified reproducible: an unchanged example (circle.diffuse) still passes against its committed gold at 1e-7 in this environment, and the new logs match independently-generated runs exactly. Co-Authored-By: stanmoore1 Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01H4ebZSujHcJJhY6bNMnmZ6 --- ...pi_1.vibrate => log.07Jul26.mpi_1.vibrate} | 69 +++++++++-------- ...pi_4.vibrate => log.07Jul26.mpi_4.vibrate} | 77 ++++++++++--------- 2 files changed, 74 insertions(+), 72 deletions(-) rename examples/vibrate/{log.12Jul24.mpi_1.vibrate => log.07Jul26.mpi_1.vibrate} (56%) rename examples/vibrate/{log.12Jul24.mpi_4.vibrate => log.07Jul26.mpi_4.vibrate} (54%) diff --git a/examples/vibrate/log.12Jul24.mpi_1.vibrate b/examples/vibrate/log.07Jul26.mpi_1.vibrate similarity index 56% rename from examples/vibrate/log.12Jul24.mpi_1.vibrate rename to examples/vibrate/log.07Jul26.mpi_1.vibrate index 5a47d9b7f..263d74f9f 100644 --- a/examples/vibrate/log.12Jul24.mpi_1.vibrate +++ b/examples/vibrate/log.07Jul26.mpi_1.vibrate @@ -1,4 +1,4 @@ -SPARTA (7 Mar 2024) +SPARTA (24 Sep 2025) Running on 1 MPI task(s) ################################################################################ # test of vibrational energy modes @@ -22,8 +22,8 @@ create_box 0 1e-5 0 1e-5 0 1e-5 Created orthogonal box = (0 0 0) to (1e-05 1e-05 1e-05) create_grid 2 2 2 Created 8 child grid cells - CPU time = 0.000811451 secs - create/ghost percent = 98.6133 1.38665 + CPU time = 0.000814469 secs + create/ghost percent = 97.2833 2.71674 species co2.species N2 CO2 vibfile co2.species.vib @@ -39,7 +39,7 @@ fix 1 vibmode create_particles mix n 0 twopass Created 1000 particles - CPU time = 0.00164516 secs + CPU time = 0.00144881 secs variable collrate equal "ncoll*step*2/np" @@ -61,63 +61,64 @@ stats_style step cpu np nattempt ncoll c_St c_Sr c_Srv c_Sv stats 10 run 100 -WARNING: Using compute grid tvib with fix vibmode may give incorrect temperature, use compute tvib/grid instead (/home/runner/work/sparta/sparta/src/compute_grid.cpp:203) +WARNING: Using compute grid tvib with fix vibmode may give incorrect temperature, use compute tvib/grid instead (/home/user/sparta/src/compute_grid.cpp:203) Memory usage per proc in Mbytes: - particles (ave,min,max) = 1.9375 1.9375 1.9375 + particles (ave,min,max) = 1.8125 1.8125 1.8125 grid (ave,min,max) = 1.51379 1.51379 1.51379 surf (ave,min,max) = 0 0 0 - total (ave,min,max) = 3.45288 3.45288 3.45288 + total (ave,min,max) = 3.32788 3.32788 3.32788 Step CPU Np Natt Ncoll c_St c_Sr c_Srv c_Sv - 0 0 1000 0 0 20527.039 20857.18 0 0 - 10 0.002129962 1000 39 26 20205.949 20691.401 204.40929 661.99501 - 20 0.004260896 1000 41 32 19568.703 20782.232 415.73498 964.11516 - 30 0.006338579 1000 41 20 19231.273 20644.583 622.01036 1289.9613 - 40 0.008423717 1000 45 30 18809.903 20632.587 774.40678 1458.8993 - 50 0.010497413 1000 44 28 18453.821 20329.389 988.81379 1699.0558 - 60 0.012657381 1000 42 30 18346.728 20057.611 1138.0057 1876.238 - 70 0.014718163 1000 43 25 18087.449 19796.238 1276.5398 2033.2845 - 80 0.016783203 1000 43 28 18095.118 19426.564 1378.0915 2145.6484 - 90 0.018804391 1000 43 26 17644.746 19502.423 1565.1087 2347.3082 - 100 0.020832432 1000 43 31 17636.431 19349.211 1662.1014 2456.254 -Loop time of 0.0208741 on 1 procs for 100 steps with 1000 particles + 0 0 1000 0 0 20307.245 20991.439 0 0 + 10 0.001701039 1000 41 31 20136.929 20715.968 153.70002 553.49515 + 20 0.003346812 1000 42 27 19684.277 20809.904 300.58986 827.67411 + 30 0.005084473 1000 44 30 19189.196 20798.65 489.49675 1104.0232 + 40 0.006898945 1000 44 28 18700.169 20948.487 659.87707 1301.6076 + 50 0.008736291 1000 42 25 18590.564 20466.038 878.81439 1563.6939 + 60 0.010530764 1000 42 29 17949.78 20613.8 1046.2097 1750.0843 + 70 0.012353185 1000 43 28 17761.724 20111.423 1200.1754 1922.2281 + 80 0.014190234 1000 44 26 17863.5 20005.609 1289.8841 2023.4138 + 90 0.015823542 1000 43 25 17620.526 19846.638 1382.3751 2124.3913 + 100 0.017426353 1000 45 27 17533.167 19573.719 1453.4421 2205.8191 +Loop time of 0.0174433 on 1 procs for 100 steps with 1000 particles +Performance: 5732.876 timesteps/s, 5.733 Mparticle-step/s MPI task timing breakdown: Section | min time | avg time | max time |%varavg| %total --------------------------------------------------------------- -Move | 0.018738 | 0.018738 | 0.018738 | 0.0 | 89.77 -Coll | 0.001148 | 0.001148 | 0.001148 | 0.0 | 5.50 -Sort | 0.00020238 | 0.00020238 | 0.00020238 | 0.0 | 0.97 -Comm | 1.3567e-05 | 1.3567e-05 | 1.3567e-05 | 0.0 | 0.06 +Move | 0.015794 | 0.015794 | 0.015794 | 0.0 | 90.54 +Coll | 0.00097996 | 0.00097996 | 0.00097996 | 0.0 | 5.62 +Sort | 0.00013298 | 0.00013298 | 0.00013298 | 0.0 | 0.76 +Comm | 2.2031e-05 | 2.2031e-05 | 2.2031e-05 | 0.0 | 0.13 Modify | 0 | 0 | 0 | 0.0 | 0.00 -Output | 0.0007668 | 0.0007668 | 0.0007668 | 0.0 | 3.67 -Other | | 4.933e-06 | | | 0.02 +Output | 0.00050734 | 0.00050734 | 0.00050734 | 0.0 | 2.91 +Other | | 6.989e-06 | | | 0.04 Particle moves = 100000 (0.1M) -Cells touched = 332469 (0.332M) +Cells touched = 328596 (0.329M) Particle comms = 0 (0K) -Boundary collides = 232448 (0.232M) +Boundary collides = 228741 (0.229M) Boundary exits = 0 (0K) SurfColl checks = 0 (0K) SurfColl occurs = 0 (0K) Surf reactions = 0 (0K) -Collide attempts = 4239 (4.24K) -Collide occurs = 2938 (2.94K) +Collide attempts = 4218 (4.22K) +Collide occurs = 2876 (2.88K) Reactions = 0 (0K) Particles stuck = 0 Axisymm bad moves = 0 -Particle-moves/CPUsec/proc: 4.79062e+06 +Particle-moves/CPUsec/proc: 5.73288e+06 Particle-moves/step: 1000 -Cell-touches/particle/step: 3.32469 +Cell-touches/particle/step: 3.28596 Particle comm iterations/step: 1 Particle fraction communicated: 0 -Particle fraction colliding with boundary: 2.32448 +Particle fraction colliding with boundary: 2.28741 Particle fraction exiting boundary: 0 Surface-checks/particle/step: 0 Surface-collisions/particle/step: 0 Surf-reactions/particle/step: 0 -Collision-attempts/particle/step: 0.04239 -Collisions/particle/step: 0.02938 +Collision-attempts/particle/step: 0.04218 +Collisions/particle/step: 0.02876 Reactions/particle/step: 0 Particles: 1000 ave 1000 max 1000 min diff --git a/examples/vibrate/log.12Jul24.mpi_4.vibrate b/examples/vibrate/log.07Jul26.mpi_4.vibrate similarity index 54% rename from examples/vibrate/log.12Jul24.mpi_4.vibrate rename to examples/vibrate/log.07Jul26.mpi_4.vibrate index dcc8c5ff5..358792f73 100644 --- a/examples/vibrate/log.12Jul24.mpi_4.vibrate +++ b/examples/vibrate/log.07Jul26.mpi_4.vibrate @@ -1,4 +1,4 @@ -SPARTA (7 Mar 2024) +SPARTA (24 Sep 2025) Running on 4 MPI task(s) ################################################################################ # test of vibrational energy modes @@ -22,8 +22,8 @@ create_box 0 1e-5 0 1e-5 0 1e-5 Created orthogonal box = (0 0 0) to (1e-05 1e-05 1e-05) create_grid 2 2 2 Created 8 child grid cells - CPU time = 0.000860103 secs - create/ghost percent = 93.2626 6.73745 + CPU time = 0.000879556 secs + create/ghost percent = 91.6968 8.30317 species co2.species N2 CO2 vibfile co2.species.vib @@ -39,7 +39,7 @@ fix 1 vibmode create_particles mix n 0 twopass Created 1000 particles - CPU time = 0.0012523 secs + CPU time = 0.00118294 secs variable collrate equal "ncoll*step*2/np" @@ -61,67 +61,68 @@ stats_style step cpu np nattempt ncoll c_St c_Sr c_Srv c_Sv stats 10 run 100 -WARNING: Using compute grid tvib with fix vibmode may give incorrect temperature, use compute tvib/grid instead (/home/runner/work/sparta/sparta/src/compute_grid.cpp:203) +WARNING: Using compute grid tvib with fix vibmode may give incorrect temperature, use compute tvib/grid instead (/home/user/sparta/src/compute_grid.cpp:203) Memory usage per proc in Mbytes: - particles (ave,min,max) = 1.9375 1.9375 1.9375 + particles (ave,min,max) = 1.8125 1.8125 1.8125 grid (ave,min,max) = 1.51379 1.51379 1.51379 surf (ave,min,max) = 0 0 0 - total (ave,min,max) = 3.45169 3.45169 3.45169 + total (ave,min,max) = 3.32669 3.32669 3.32669 Step CPU Np Natt Ncoll c_St c_Sr c_Srv c_Sv - 0 0 1000 0 0 19750.722 21250.869 0 0 - 10 0.001099912 1000 40 29 19584.106 20569.578 179.94805 581.31383 - 20 0.00204286 1000 45 35 19517.758 20359.032 333.16715 872.99551 - 30 0.003018749 1000 45 35 19023.886 20314.638 495.38299 1085.8964 - 40 0.00399013 1000 43 29 18478.529 20230.419 775.61845 1437.0413 - 50 0.004915165 1000 43 33 18189.411 20164.428 907.17702 1580.6091 - 60 0.005841441 1000 44 28 17823.085 20197.993 1059.5758 1775.4883 - 70 0.006737521 1000 44 32 17559.06 19893.231 1235.0641 1970.1942 - 80 0.007771981 1000 44 29 16907.053 20287.099 1438.5541 2191.7101 - 90 0.008669063 1000 43 25 16632.298 19978.181 1598.7441 2367.6259 - 100 0.009600359 1000 45 27 16568.79 19614.025 1741.1069 2529.6526 -Loop time of 0.0096493 on 4 procs for 100 steps with 1000 particles + 0 0 1000 0 0 20604.382 21301.272 0 0 + 10 0.001017648 1000 43 32 20499.65 20858.013 150.15022 522.48317 + 20 0.001881933 1000 41 32 20509.003 20359.192 322.99105 868.56297 + 30 0.002708559 1000 41 27 20109.663 20119.263 499.14626 1083.5941 + 40 0.003511703 1000 43 27 19613.153 20167.821 667.20383 1317.1638 + 50 0.004320494 1000 44 31 19352.856 20248.73 795.75074 1469.0043 + 60 0.005139096 1000 46 33 18972.622 20395.593 953.88246 1651.0698 + 70 0.005959587 1000 43 31 18602.262 20130.448 1107.9058 1823.2868 + 80 0.006810711 1000 43 31 18193.647 20309.316 1280.61 2018.1609 + 90 0.007636682 1000 45 30 17604.824 20072.144 1467.603 2220.0972 + 100 0.008497566 1000 43 25 17390.064 19762.296 1727.1455 2495.3107 +Loop time of 0.00853214 on 4 procs for 100 steps with 1000 particles +Performance: 11720.384 timesteps/s, 11.720 Mparticle-step/s MPI task timing breakdown: Section | min time | avg time | max time |%varavg| %total --------------------------------------------------------------- -Move | 0.0045529 | 0.0045854 | 0.0045998 | 0.0 | 47.52 -Coll | 0.00032341 | 0.00032793 | 0.00033003 | 0.0 | 3.40 -Sort | 7.2893e-05 | 7.7937e-05 | 8.4618e-05 | 0.0 | 0.81 -Comm | 0.0038536 | 0.0041274 | 0.0042462 | 0.2 | 42.77 +Move | 0.0036794 | 0.0038398 | 0.0040027 | 0.2 | 45.00 +Coll | 0.00023008 | 0.000236 | 0.00024157 | 0.0 | 2.77 +Sort | 4.0392e-05 | 4.2848e-05 | 4.7076e-05 | 0.0 | 0.50 +Comm | 0.0039121 | 0.0040888 | 0.0043192 | 0.2 | 47.92 Modify | 0 | 0 | 0 | 0.0 | 0.00 -Output | 0.0003919 | 0.00049521 | 0.00078913 | 0.0 | 5.13 -Other | | 3.536e-05 | | | 0.37 +Output | 0.00022122 | 0.00029706 | 0.00051278 | 0.0 | 3.48 +Other | | 2.768e-05 | | | 0.32 Particle moves = 100000 (0.1M) -Cells touched = 326495 (0.326M) -Particle comms = 74165 (74.2K) -Boundary collides = 226537 (0.227M) +Cells touched = 333151 (0.333M) +Particle comms = 74058 (74.1K) +Boundary collides = 233026 (0.233M) Boundary exits = 0 (0K) SurfColl checks = 0 (0K) SurfColl occurs = 0 (0K) Surf reactions = 0 (0K) -Collide attempts = 4218 (4.22K) -Collide occurs = 2938 (2.94K) +Collide attempts = 4287 (4.29K) +Collide occurs = 2954 (2.95K) Reactions = 0 (0K) Particles stuck = 0 Axisymm bad moves = 0 -Particle-moves/CPUsec/proc: 2.59086e+06 +Particle-moves/CPUsec/proc: 2.9301e+06 Particle-moves/step: 1000 -Cell-touches/particle/step: 3.26495 +Cell-touches/particle/step: 3.33151 Particle comm iterations/step: 1 -Particle fraction communicated: 0.74165 -Particle fraction colliding with boundary: 2.26537 +Particle fraction communicated: 0.74058 +Particle fraction colliding with boundary: 2.33026 Particle fraction exiting boundary: 0 Surface-checks/particle/step: 0 Surface-collisions/particle/step: 0 Surf-reactions/particle/step: 0 -Collision-attempts/particle/step: 0.04218 -Collisions/particle/step: 0.02938 +Collision-attempts/particle/step: 0.04287 +Collisions/particle/step: 0.02954 Reactions/particle/step: 0 -Particles: 250 ave 267 max 230 min -Histogram: 1 0 1 0 0 0 0 0 0 2 +Particles: 250 ave 272 max 206 min +Histogram: 1 0 0 0 0 0 0 1 1 1 Cells: 2 ave 2 max 2 min Histogram: 4 0 0 0 0 0 0 0 0 0 GhostCell: 6 ave 6 max 6 min