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 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; } 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; }