From 456f0be998aa791719c6d402057a984859de423e Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 20 Aug 2026 15:04:52 +0000 Subject: [PATCH 1/4] Consolidate the verified static-analysis bug fixes from both review branches Consolidates claude/review-sparta-analysis-bugs-kyionn and claude/sparta-static-analysis-bugs-335wsw onto this branch, keeping only fixes for genuine defects still present here, re-verified line-by-line against the full 111-entry bug list (74 automated site checks, plus reachability proofs for every rejected entry): - wrong-logic fixes: tally keyword/enum typos, compute count group check, tvib mode indexing, property/surf OOB, reaction2col wipe, rotc2 symmetry, vremax==0 NaN acceptance, react_qk/tce_qk react_prob scratch pollution, react_tce unreachable warnings + missing break, react_bird uninit tallies, adapt_grid coarsen newcell index, marching-cubes int corners, specular wrapper noslip, emit_face_file azimuth MY_2PI, tally_setup dangling *_active - KOKKOS: RNG state freed then reused in collide attempt loop, return-vs-continue in per-cell particle loops, lambda/grid missing else (also CPU), sonine dead OOB read, compute_surf tallyinfo OOB when nsurf==0, ave/histo(/weight) inverted realloc conditions + stray printf, grid-check OOB error message, compute_surf ECHEM conditional column advance (also CPU), bin_particles uninitialized Region pointer when no region is defined - crash/hang/NaN: mover frac 0/0 + clamp, stuck-particle epsilon, axi_remap axis guard, axisymmetric quadratic cancellation (Vieta) and a==0 linear case, temp rescale t_current==0, ablate Ninterface==0, subsonic emit 0-division, fix_surf_temp uninit units + stale cqw/fqw, variable lock-file fopen/fscanf, read_isurf unchecked fopen - hardening: snprintf for unbounded %s filename/name messages across the code base (input, write_*, read_*, sparta, output restart names with BIGINT timesteps, domain/modify/surf suffix styles, mixture, fix_print, surf_react_prob, PYTHON), truncation-safe stats.cpp line/format building, str[32] cell-ID buffers, BIGINT_FORMAT in fix_halt - leaks on recoverable error paths (errors throw for library callers): variable.cpp ids, suffix/list/commands leaks, write_isurf arg mutation + file leak; compute_reduce replace/subset bounds vs expanded args - KOKKOS FFT convention: fft_2d_1d_only_kokkos and fft_3d_1d_only_kokkos ran forward plans for flag==-1 and normalized on flag==1, opposite of the host fft2d.cpp/fft3d.cpp and inconsistent with their own cuFFT/hipFFT branches; both now match the CPU convention (timing-only) Fixes already present on this branch (from newer master work) were skipped, as were entries rejected with verified rationale: inert integer casts, unreachable guards (collide volume, dt/grid vrm_max, ecc, adsorb zero-velocity, piston uprime==0), RNG seed/idiom changes that shift reproducibility baselines, the upstream fft2d one-line scaling flip, and cosmetic snprintf conversions of provably bounded formats. Items flagged for domain sign-off: react_tce break, geometry Vieta/a==0, emit azimuth. Verified with a clean serial build. Co-authored-by: aborner1 <42391399+aborner1@users.noreply.github.com> Co-authored-by: Stan Moore Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01JAGuLGWVKDbVYbwNXjqZ6K --- src/KOKKOS/collide_vss_kokkos.cpp | 2 - src/KOKKOS/compute_eflux_grid_kokkos.cpp | 2 +- src/KOKKOS/compute_lambda_grid_kokkos.cpp | 4 +- src/KOKKOS/compute_sonine_grid_kokkos.cpp | 3 - src/KOKKOS/compute_surf_kokkos.cpp | 2 +- src/KOKKOS/compute_surf_kokkos.h | 3 +- src/KOKKOS/compute_thermal_grid_kokkos.cpp | 2 +- src/KOKKOS/fft2d_kokkos.cpp | 10 +-- src/KOKKOS/fft3d_kokkos.cpp | 10 +-- src/KOKKOS/fix_ave_histo_kokkos.cpp | 38 +++++------ src/KOKKOS/fix_ave_histo_weight_kokkos.cpp | 39 ++++++------ src/KOKKOS/fix_grid_check_kokkos.cpp | 4 +- src/KOKKOS/geometry_kokkos.h | 60 +++++++++++++----- src/KOKKOS/update_kokkos.cpp | 32 +++++++--- src/KOKKOS/update_kokkos.h | 14 +++-- src/PYTHON/python_impl.cpp | 38 +++++------ src/adapt_grid.cpp | 3 +- src/collide_vss.cpp | 3 +- src/compute_count.cpp | 2 +- src/compute_gas_collision_tally.cpp | 2 +- src/compute_gas_reaction_tally.cpp | 4 +- src/compute_lambda_grid.cpp | 8 ++- src/compute_property_surf.cpp | 2 +- src/compute_react_boundary.cpp | 1 - src/compute_react_isurf_grid.cpp | 1 - src/compute_react_surf.cpp | 1 - src/compute_reduce.cpp | 8 ++- src/compute_surf.cpp | 3 +- src/compute_surf_reaction_tally.cpp | 2 +- src/compute_tvib_grid.cpp | 2 +- src/domain.cpp | 2 +- src/dump.cpp | 2 +- src/dump_grid.cpp | 2 +- src/dump_movie.cpp | 4 +- src/fix_ablate.cpp | 4 +- src/fix_ablate_multi_inner.cpp | 6 +- src/fix_ave_grid.cpp | 4 +- src/fix_ave_histo.cpp | 4 +- src/fix_ave_surf.cpp | 4 +- src/fix_ave_time.cpp | 4 +- src/fix_emit_face.cpp | 2 +- src/fix_emit_face_file.cpp | 4 +- src/fix_emit_surf.cpp | 2 +- src/fix_grid_check.cpp | 4 +- src/fix_halt.cpp | 2 +- src/fix_print.cpp | 2 +- src/fix_surf_temp.cpp | 12 +++- src/fix_temp_rescale.cpp | 7 ++- src/geometry.cpp | 60 +++++++++++++----- src/grid.cpp | 15 ++++- src/grid_id.cpp | 3 +- src/input.cpp | 22 +++---- src/marching_cubes.h | 2 +- src/mixture.cpp | 2 +- src/modify.cpp | 4 +- src/move_surf.cpp | 2 +- src/output.cpp | 4 +- src/particle.cpp | 6 +- src/react_bird.cpp | 2 + src/react_qk.cpp | 6 +- src/react_tce.cpp | 8 ++- src/react_tce_qk.cpp | 6 +- src/read_grid.cpp | 4 +- src/read_isurf.cpp | 6 ++ src/read_restart.cpp | 16 ++--- src/read_surf.cpp | 8 +-- src/sparta.cpp | 12 ++-- src/stats.cpp | 22 ++++--- src/surf.cpp | 4 +- src/surf_collide_specular.cpp | 3 +- src/surf_react_prob.cpp | 2 +- src/update.cpp | 33 +++++++--- src/update.h | 14 +++-- src/utils.cpp | 2 +- src/variable.cpp | 73 +++++++++++++++++----- src/write_grid.cpp | 2 +- src/write_isurf.cpp | 13 ++-- src/write_restart.cpp | 8 +-- src/write_surf.cpp | 4 +- 79 files changed, 472 insertions(+), 272 deletions(-) diff --git a/src/KOKKOS/collide_vss_kokkos.cpp b/src/KOKKOS/collide_vss_kokkos.cpp index fcaf48986..06ab9bab8 100644 --- a/src/KOKKOS/collide_vss_kokkos.cpp +++ b/src/KOKKOS/collide_vss_kokkos.cpp @@ -982,7 +982,6 @@ void CollideVSSKokkos::operator()(TagCollideCollisionsOne< NEARCP, GASTALLY, ATO else reduce.nreact_one++; } else { - rand_pool.free_state(rand_gen); continue; } @@ -2590,7 +2589,6 @@ void CollideVSSKokkos::operator()(TagCollideCollisionsOneAmbipolar< GASTALLY, AT else reduce.nreact_one++; } else { - rand_pool.free_state(rand_gen); continue; } diff --git a/src/KOKKOS/compute_eflux_grid_kokkos.cpp b/src/KOKKOS/compute_eflux_grid_kokkos.cpp index d0436d9d5..d0d02659a 100644 --- a/src/KOKKOS/compute_eflux_grid_kokkos.cpp +++ b/src/KOKKOS/compute_eflux_grid_kokkos.cpp @@ -241,7 +241,7 @@ void ComputeEFluxGridKokkos::operator()(TagComputeEFluxGrid_compute_per_grid, co const int ispecies = d_particles[i].ispecies; const int igroup = d_s2g(imix,ispecies); - if (igroup < 0) return; + if (igroup < 0) continue; const double mass = d_species[ispecies].mass; double *v = d_particles[i].v; diff --git a/src/KOKKOS/compute_lambda_grid_kokkos.cpp b/src/KOKKOS/compute_lambda_grid_kokkos.cpp index f988a63fb..2f2243743 100644 --- a/src/KOKKOS/compute_lambda_grid_kokkos.cpp +++ b/src/KOKKOS/compute_lambda_grid_kokkos.cpp @@ -326,12 +326,12 @@ void ComputeLambdaGridKokkos::compute_per_grid_kokkos() if (l_knyflag) { if (l_noutputs == 1) l_vector_grid[i] = lambda / sizey; - l_array_grid(i,l_output_order[KNY]) = lambda / sizey; + else l_array_grid(i,l_output_order[KNY]) = lambda / sizey; } if (l_knzflag) { if (l_noutputs == 1) l_vector_grid[i] = lambda / sizez; - l_array_grid(i,l_output_order[KNZ]) = lambda / sizez; + else l_array_grid(i,l_output_order[KNZ]) = lambda / sizez; } }); } diff --git a/src/KOKKOS/compute_sonine_grid_kokkos.cpp b/src/KOKKOS/compute_sonine_grid_kokkos.cpp index ce786ee9c..0ffa1bc21 100644 --- a/src/KOKKOS/compute_sonine_grid_kokkos.cpp +++ b/src/KOKKOS/compute_sonine_grid_kokkos.cpp @@ -225,9 +225,6 @@ void ComputeSonineGridKokkos::operator()(TagComputeSonineGrid_compute_vcom, cons KOKKOS_INLINE_FUNCTION void ComputeSonineGridKokkos::operator()(TagComputeSonineGrid_normalize_vcom, const int &icell) const { - const int ispecies = d_particles[icell].ispecies; - const int igroup = d_s2g(imix,ispecies); - double norm; for (int j=0; j 0) iend--; + while (iend > 0 && h_surf2tally[iend] == -1) iend--; if (istart >= iend) { ntally = istart; break; diff --git a/src/KOKKOS/compute_surf_kokkos.h b/src/KOKKOS/compute_surf_kokkos.h index 9f2355ad7..25021fe44 100644 --- a/src/KOKKOS/compute_surf_kokkos.h +++ b/src/KOKKOS/compute_surf_kokkos.h @@ -409,8 +409,9 @@ void surf_tally_kk(double /*dtremain*/, int isurf, int icell, int reaction, double r_coeff = 0.0; if (sr_type == 1) r_coeff = sr_kk_prob_copy[m].obj.d_coeffs(reaction-1,1); - a_array_surf_tally(itally,k++) += weight * r_coeff * fluxscale; + a_array_surf_tally(itally,k) += weight * r_coeff * fluxscale; } + k++; break; case ETOT: if (iorig) vsqpre = origmass * MathExtraKokkos::lensq3(vorig); diff --git a/src/KOKKOS/compute_thermal_grid_kokkos.cpp b/src/KOKKOS/compute_thermal_grid_kokkos.cpp index 0a2db7ca2..0aa7a4aa5 100644 --- a/src/KOKKOS/compute_thermal_grid_kokkos.cpp +++ b/src/KOKKOS/compute_thermal_grid_kokkos.cpp @@ -166,7 +166,7 @@ void ComputeThermalGridKokkos::operator()(TagComputeThermalGrid_compute_per_grid const int ispecies = d_particles[i].ispecies; const int igroup = d_s2g(imix,ispecies); - if (igroup < 0) return; + if (igroup < 0) continue; const int icell = d_particles[i].icell; diff --git a/src/KOKKOS/fft2d_kokkos.cpp b/src/KOKKOS/fft2d_kokkos.cpp index 038947b77..fd274a5a1 100644 --- a/src/KOKKOS/fft2d_kokkos.cpp +++ b/src/KOKKOS/fft2d_kokkos.cpp @@ -714,7 +714,7 @@ void FFT2dKokkos::fft_2d_1d_only_kokkos(typename FFT_AT::t_FFT_DATA_ // data is just an array of 0.0 #if defined(FFT_KOKKOS_MKL_GPU) - if (flag == -1) { + if (flag == 1) { oneapi::mkl::dft::compute_forward(*(plan->desc_fast), (FFT_SCALAR*)d_data.data()); oneapi::mkl::dft::compute_forward(*(plan->desc_slow), (FFT_SCALAR*)d_data.data()); } else { @@ -722,7 +722,7 @@ void FFT2dKokkos::fft_2d_1d_only_kokkos(typename FFT_AT::t_FFT_DATA_ oneapi::mkl::dft::compute_backward(*(plan->desc_slow), (FFT_SCALAR*)d_data.data()); } #elif defined(FFT_KOKKOS_MKL) - if (flag == -1) { + if (flag == 1) { DftiComputeForward(plan->handle_fast,d_data.data()); DftiComputeForward(plan->handle_slow,d_data.data()); } else { @@ -730,7 +730,7 @@ void FFT2dKokkos::fft_2d_1d_only_kokkos(typename FFT_AT::t_FFT_DATA_ DftiComputeBackward(plan->handle_slow,d_data.data()); } #elif defined(FFT_KOKKOS_FFTW3) || defined(FFT_KOKKOS_NVPL) - if (flag == -1) { + if (flag == 1) { FFTW_API(execute_dft)(plan->plan_fast_forward,(FFT_KOKKOS_DATA*)d_data.data(),(FFT_KOKKOS_DATA*)d_data.data()); FFTW_API(execute_dft)(plan->plan_slow_forward,(FFT_KOKKOS_DATA*)d_data.data(),(FFT_KOKKOS_DATA*)d_data.data()); } else { @@ -747,7 +747,7 @@ void FFT2dKokkos::fft_2d_1d_only_kokkos(typename FFT_AT::t_FFT_DATA_ kiss_fft_functor f; typename FFT_AT::t_FFT_DATA_1d d_tmp = typename FFT_AT::t_FFT_DATA_1d(Kokkos::view_alloc("fft_2d:tmp",Kokkos::WithoutInitializing),d_data.extent(0)); - if (flag == -1) { + if (flag == 1) { f = kiss_fft_functor(d_data,d_tmp,plan->cfg_fast_forward,length1); Kokkos::parallel_for(total1/length1,f); @@ -765,7 +765,7 @@ void FFT2dKokkos::fft_2d_1d_only_kokkos(typename FFT_AT::t_FFT_DATA_ // scaling if required // limit num to size of data - if (flag == 1 && plan->scaled) { + if (flag == -1 && plan->scaled) { FFT_SCALAR norm = plan->norm; int num = MIN(plan->normnum,nsize); diff --git a/src/KOKKOS/fft3d_kokkos.cpp b/src/KOKKOS/fft3d_kokkos.cpp index 74418f2c2..cfa854912 100644 --- a/src/KOKKOS/fft3d_kokkos.cpp +++ b/src/KOKKOS/fft3d_kokkos.cpp @@ -915,7 +915,7 @@ void FFT3dKokkos::fft_3d_1d_only_kokkos(typename FFT_AT::t_FFT_DATA_ // data is just an array of 0.0 #if defined(FFT_KOKKOS_MKL_GPU) - if (flag == -1) { + if (flag == 1) { oneapi::mkl::dft::compute_forward(*(plan->desc_fast), (FFT_SCALAR*)d_data.data()); oneapi::mkl::dft::compute_forward(*(plan->desc_mid), (FFT_SCALAR*)d_data.data()); oneapi::mkl::dft::compute_forward(*(plan->desc_slow), (FFT_SCALAR*)d_data.data()); @@ -925,7 +925,7 @@ void FFT3dKokkos::fft_3d_1d_only_kokkos(typename FFT_AT::t_FFT_DATA_ oneapi::mkl::dft::compute_backward(*(plan->desc_slow), (FFT_SCALAR*)d_data.data()); } #elif defined(FFT_KOKKOS_MKL) - if (flag == -1) { + if (flag == 1) { DftiComputeForward(plan->handle_fast,d_data.data()); DftiComputeForward(plan->handle_mid,d_data.data()); DftiComputeForward(plan->handle_slow,d_data.data()); @@ -935,7 +935,7 @@ void FFT3dKokkos::fft_3d_1d_only_kokkos(typename FFT_AT::t_FFT_DATA_ DftiComputeBackward(plan->handle_slow,d_data.data()); } #elif defined(FFT_KOKKOS_FFTW3) || defined(FFT_KOKKOS_NVPL) - if (flag == -1) { + if (flag == 1) { FFTW_API(execute_dft)(plan->plan_fast_forward,(FFT_KOKKOS_DATA*)d_data.data(),(FFT_KOKKOS_DATA*)d_data.data()); FFTW_API(execute_dft)(plan->plan_mid_forward,(FFT_KOKKOS_DATA*)d_data.data(),(FFT_KOKKOS_DATA*)d_data.data()); FFTW_API(execute_dft)(plan->plan_slow_forward,(FFT_KOKKOS_DATA*)d_data.data(),(FFT_KOKKOS_DATA*)d_data.data()); @@ -956,7 +956,7 @@ void FFT3dKokkos::fft_3d_1d_only_kokkos(typename FFT_AT::t_FFT_DATA_ kiss_fft_functor f; typename FFT_AT::t_FFT_DATA_1d d_tmp = typename FFT_AT::t_FFT_DATA_1d(Kokkos::view_alloc("fft_3d:tmp",Kokkos::WithoutInitializing),d_data.extent(0)); - if (flag == -1) { + if (flag == 1) { f = kiss_fft_functor(d_data,d_tmp,plan->cfg_fast_forward,length1); Kokkos::parallel_for(total1/length1,f); @@ -980,7 +980,7 @@ void FFT3dKokkos::fft_3d_1d_only_kokkos(typename FFT_AT::t_FFT_DATA_ // scaling if required // limit num to size of data - if (flag == 1 && plan->scaled) { + if (flag == -1 && plan->scaled) { FFT_SCALAR norm = plan->norm; int num = MIN(plan->normnum,nsize); diff --git a/src/KOKKOS/fix_ave_histo_kokkos.cpp b/src/KOKKOS/fix_ave_histo_kokkos.cpp index 89bab835a..2322f9a29 100644 --- a/src/KOKKOS/fix_ave_histo_kokkos.cpp +++ b/src/KOKKOS/fix_ave_histo_kokkos.cpp @@ -505,19 +505,20 @@ void FixAveHistoKokkos::bin_particles( int n = particle->nlocal; int nmax = particle->maxlocal; - Region *region; - if (regionflag) region = domain->regions[iregion]; + if (regionflag) { + Region *region = domain->regions[iregion]; - if (!region->kokkos_flag) - error->all(FLERR,"KOKKOS package does not (yet) support chosen region style"); + if (!region->kokkos_flag) + error->all(FLERR,"KOKKOS package does not (yet) support chosen region style"); - KokkosBase* regionKKBase = dynamic_cast(region); + KokkosBase* regionKKBase = dynamic_cast(region); - if (k_match.extent(0) > nmax) - MemKK::realloc_kokkos(k_match,"fix_ave_histo_weight:match",nmax); + if (k_match.extent(0) < nmax) + MemKK::realloc_kokkos(k_match,"fix_ave_histo_weight:match",nmax); - regionKKBase->match_all_kokkos(k_match); - d_match = k_match.view_device(); + regionKKBase->match_all_kokkos(k_match); + d_match = k_match.view_device(); + } if (attribute == X) { @@ -569,19 +570,20 @@ void FixAveHistoKokkos::bin_particles( d_values = mirror_view_from_raw_host_array(values, n, stride); - Region *region; - if (regionflag) region = domain->regions[iregion]; + if (regionflag) { + Region *region = domain->regions[iregion]; - if (!region->kokkos_flag) - error->all(FLERR,"KOKKOS package does not (yet) support chosen region style"); + if (!region->kokkos_flag) + error->all(FLERR,"KOKKOS package does not (yet) support chosen region style"); - KokkosBase* regionKKBase = dynamic_cast(region); + KokkosBase* regionKKBase = dynamic_cast(region); - if (k_match.extent(0) < nmax) - MemKK::realloc_kokkos(k_match,"fix_ave_histo_weight:match",nmax); + if (k_match.extent(0) < nmax) + MemKK::realloc_kokkos(k_match,"fix_ave_histo_weight:match",nmax); - regionKKBase->match_all_kokkos(k_match); - d_match = k_match.view_device(); + regionKKBase->match_all_kokkos(k_match); + d_match = k_match.view_device(); + } if (regionflag && mixflag) { auto policy = RangePolicy(0, n); diff --git a/src/KOKKOS/fix_ave_histo_weight_kokkos.cpp b/src/KOKKOS/fix_ave_histo_weight_kokkos.cpp index ee0f6e638..7dde85694 100644 --- a/src/KOKKOS/fix_ave_histo_weight_kokkos.cpp +++ b/src/KOKKOS/fix_ave_histo_weight_kokkos.cpp @@ -293,7 +293,6 @@ void FixAveHistoWeightKokkos::calculate_weights() // explicit per-particle attributes // NOTE: need to allocate local storage } else { - printf("%d, %d\n", which[i] == VARIABLE, kind == PERGRID); error->all(FLERR,"Fix ave/histo/weight/kokkos option not yet supported"); } } @@ -329,19 +328,20 @@ void FixAveHistoWeightKokkos::bin_particles( int n = particle->nlocal; int nmax = particle->maxlocal; - Region *region; - if (regionflag) region = domain->regions[iregion]; + if (regionflag) { + Region *region = domain->regions[iregion]; - if (!region->kokkos_flag) - error->all(FLERR,"KOKKOS package does not (yet) support chosen region style"); + if (!region->kokkos_flag) + error->all(FLERR,"KOKKOS package does not (yet) support chosen region style"); - KokkosBase* regionKKBase = dynamic_cast(region); + KokkosBase* regionKKBase = dynamic_cast(region); - if (k_match.extent(0) < nmax) - MemKK::realloc_kokkos(k_match,"fix_ave_histo_weight:match",nmax); + if (k_match.extent(0) < nmax) + MemKK::realloc_kokkos(k_match,"fix_ave_histo_weight:match",nmax); - regionKKBase->match_all_kokkos(k_match); - d_match = k_match.view_device(); + regionKKBase->match_all_kokkos(k_match); + d_match = k_match.view_device(); + } if (attribute == X) { @@ -393,19 +393,20 @@ void FixAveHistoWeightKokkos::bin_particles( d_values = mirror_view_from_raw_host_array(values, n, stride); - Region *region; - if (regionflag) region = domain->regions[iregion]; + if (regionflag) { + Region *region = domain->regions[iregion]; - if (!region->kokkos_flag) - error->all(FLERR,"KOKKOS package does not (yet) support chosen region style"); + if (!region->kokkos_flag) + error->all(FLERR,"KOKKOS package does not (yet) support chosen region style"); - KokkosBase* regionKKBase = dynamic_cast(region); + KokkosBase* regionKKBase = dynamic_cast(region); - if (k_match.extent(0) > nmax) - MemKK::realloc_kokkos(k_match,"fix_ave_histo_weight:match",nmax); + if (k_match.extent(0) < nmax) + MemKK::realloc_kokkos(k_match,"fix_ave_histo_weight:match",nmax); - regionKKBase->match_all_kokkos(k_match); - d_match = k_match.view_device(); + regionKKBase->match_all_kokkos(k_match); + d_match = k_match.view_device(); + } if (regionflag && mixflag) { //auto policy = RangePolicy(0, n); diff --git a/src/KOKKOS/fix_grid_check_kokkos.cpp b/src/KOKKOS/fix_grid_check_kokkos.cpp index ca9e0f4d4..2e2babb73 100644 --- a/src/KOKKOS/fix_grid_check_kokkos.cpp +++ b/src/KOKKOS/fix_grid_check_kokkos.cpp @@ -191,9 +191,9 @@ void FixGridCheckKokkos::end_of_step() auto icell = particles[i].icell; if (h_particle_problems(i) & IS_IN_INVALID_CELL) { sprintf(str, - "Particle %d,%d on proc %d is in invalid cell " CELLINT_FORMAT + "Particle %d,%d on proc %d is in invalid cell index %d" " on timestep " BIGINT_FORMAT, - i,particles[i].id,comm->me,cells[icell].id,update->ntimestep); + i,particles[i].id,comm->me,icell,update->ntimestep); error->one(FLERR,str); } if (h_particle_problems(i) & IS_OUTSIDE_CELL) { diff --git a/src/KOKKOS/geometry_kokkos.h b/src/KOKKOS/geometry_kokkos.h index 319fa8d7f..48f15a7ef 100644 --- a/src/KOKKOS/geometry_kokkos.h +++ b/src/KOKKOS/geometry_kokkos.h @@ -156,10 +156,20 @@ bool axi_horizontal_line(double tdelta, double *x, double *v, double arg = yhoriz*yhoriz*a - v[2]*v[2]*x[1]*x[1]; if (arg < 0.0) return false; double sarg = sqrt(arg); + double c = x[1]*x[1] - yhoriz*yhoriz; nc = 2; - double tone = (b - sarg) / a; - double ttwo = (b + sarg) / a; + double tone, ttwo; + if (b > 0.0) { + ttwo = (b + sarg) / a; + tone = c / (b + sarg); + } else if (b < 0.0) { + tone = (b - sarg) / a; + ttwo = c / (b - sarg); + } else { + tone = -sarg / a; + ttwo = sarg / a; + } t1 = MIN(tone,ttwo); t2 = MAX(tone,ttwo); @@ -265,20 +275,33 @@ bool axi_line_intersect(double tdelta, double *x, double *v, double dconst = x21*v1[1] - y21*v1[0]; double a = x21sq*(v[1]*v[1] + v[2]*v[2]) - y21sq*v[0]*v[0]; - if (a == 0.0) return false; double b = x21sq*x[1]*v[1] - y21sq*x[0]*v[0] - y21*v[0]*dconst; double c = x21sq*x[1]*x[1] - y21sq*x[0]*x[0] - 2.0*y21*x[0]*dconst - dconst*dconst; - double arg = b*b - a*c; - if (arg < 0.0) return false; - double sarg = sqrt(arg); - - nc = 2; - double tone = (-b - sarg) / a; - double ttwo = (-b + sarg) / a; - t1 = MIN(tone,ttwo); - t2 = MAX(tone,ttwo); + if (a == 0.0) { + if (b == 0.0) return false; + nc = 1; + t1 = t2 = -0.5 * c / b; + } else { + double arg = b*b - a*c; + if (arg < 0.0) return false; + double sarg = sqrt(arg); + nc = 2; + double tone, ttwo; + if (b > 0.0) { + tone = (-b - sarg) / a; + ttwo = c / (-b - sarg); + } else if (b < 0.0) { + tone = c / (-b + sarg); + ttwo = (-b + sarg) / a; + } else { + tone = -sarg / a; + ttwo = sarg / a; + } + t1 = MIN(tone,ttwo); + t2 = MAX(tone,ttwo); + } } // if selfflag, particle starts on surf line segment @@ -341,11 +364,16 @@ bool axi_line_intersect(double tdelta, double *x, double *v, if (v1[1] == v2[1]) xc[1] = v1[1]; xc[2] = 0.0; - double rn = ynew / xc[1]; - double wn = znew / xc[1]; vc[0] = v[0]; - vc[1] = v[1]*rn + v[2]*wn; - vc[2] = -v[1]*wn + v[2]*rn; + if (xc[1] > 0.0) { + double rn = ynew / xc[1]; + double wn = znew / xc[1]; + vc[1] = v[1]*rn + v[2]*wn; + vc[2] = -v[1]*wn + v[2]*rn; + } else { + vc[1] = v[1]; + vc[2] = v[2]; + } // test that xc is within line segment bounds // y-test for vertical line, else x-test diff --git a/src/KOKKOS/update_kokkos.cpp b/src/KOKKOS/update_kokkos.cpp index e216976ae..1fbecfb16 100644 --- a/src/KOKKOS/update_kokkos.cpp +++ b/src/KOKKOS/update_kokkos.cpp @@ -1515,22 +1515,34 @@ void UpdateKokkos::operator()(TagUpdateMove frac = 1.0; if (xnew[0] < lo[0]) { - frac = (lo[0]-x[0]) / (xnew[0]-x[0]); + if (xnew[0] != x[0]) frac = (lo[0]-x[0]) / (xnew[0]-x[0]); + else frac = 0.0; + if (frac < 0.0) frac = 0.0; + else if (frac > 1.0) frac = 1.0; outface = XLO; } else if (xnew[0] >= hi[0]) { - frac = (hi[0]-x[0]) / (xnew[0]-x[0]); + if (xnew[0] != x[0]) frac = (hi[0]-x[0]) / (xnew[0]-x[0]); + else frac = 0.0; + if (frac < 0.0) frac = 0.0; + else if (frac > 1.0) frac = 1.0; outface = XHI; } if (DIM != 1) { if (xnew[1] < lo[1]) { - newfrac = (lo[1]-x[1]) / (xnew[1]-x[1]); + if (xnew[1] != x[1]) newfrac = (lo[1]-x[1]) / (xnew[1]-x[1]); + else newfrac = 0.0; + if (newfrac < 0.0) newfrac = 0.0; + else if (newfrac > 1.0) newfrac = 1.0; if (newfrac < frac) { frac = newfrac; outface = YLO; } } else if (xnew[1] >= hi[1]) { - newfrac = (hi[1]-x[1]) / (xnew[1]-x[1]); + if (xnew[1] != x[1]) newfrac = (hi[1]-x[1]) / (xnew[1]-x[1]); + else newfrac = 0.0; + if (newfrac < 0.0) newfrac = 0.0; + else if (newfrac > 1.0) newfrac = 1.0; if (newfrac < frac) { frac = newfrac; outface = YHI; @@ -1573,13 +1585,19 @@ void UpdateKokkos::operator()(TagUpdateMove if (DIM == 3) { if (xnew[2] < lo[2]) { - newfrac = (lo[2]-x[2]) / (xnew[2]-x[2]); + if (xnew[2] != x[2]) newfrac = (lo[2]-x[2]) / (xnew[2]-x[2]); + else newfrac = 0.0; + if (newfrac < 0.0) newfrac = 0.0; + else if (newfrac > 1.0) newfrac = 1.0; if (newfrac < frac) { frac = newfrac; outface = ZLO; } } else if (xnew[2] >= hi[2]) { - newfrac = (hi[2]-x[2]) / (xnew[2]-x[2]); + if (xnew[2] != x[2]) newfrac = (hi[2]-x[2]) / (xnew[2]-x[2]); + else newfrac = 0.0; + if (newfrac < 0.0) newfrac = 0.0; + else if (newfrac > 1.0) newfrac = 1.0; if (newfrac < frac) { frac = newfrac; outface = ZHI; @@ -1915,7 +1933,7 @@ void UpdateKokkos::operator()(TagUpdateMove // stuck_iterate = consecutive iterations particle is immobile - if (minparam == 0.0) stuck_iterate++; + if (minparam <= 1.0e-14) stuck_iterate++; else stuck_iterate = 0; // reset post-bounce xnew diff --git a/src/KOKKOS/update_kokkos.h b/src/KOKKOS/update_kokkos.h index 91438eaae..e8f1b9391 100644 --- a/src/KOKKOS/update_kokkos.h +++ b/src/KOKKOS/update_kokkos.h @@ -276,12 +276,14 @@ class UpdateKokkos : public Update { double znew = x[2]; x[1] = sqrt(ynew*ynew + znew*znew); x[2] = 0.0; - double rn = ynew / x[1]; - double wn = znew / x[1]; - double vy = v[1]; - double vz = v[2]; - v[1] = vy*rn + vz*wn; - v[2] = -vy*wn + vz*rn; + if (x[1] > 0.0) { + double rn = ynew / x[1]; + double wn = znew / x[1]; + double vy = v[1]; + double vz = v[2]; + v[1] = vy*rn + vz*wn; + v[2] = -vy*wn + vz*rn; + } }; typedef void (UpdateKokkos::*FnPtr)(); diff --git a/src/PYTHON/python_impl.cpp b/src/PYTHON/python_impl.cpp index 849ad2032..22fd18354 100644 --- a/src/PYTHON/python_impl.cpp +++ b/src/PYTHON/python_impl.cpp @@ -112,7 +112,7 @@ void PythonImpl::command(int narg, char **arg) int ifunc = find(arg[0]); if (ifunc < 0) { char msg[128]; - sprintf(msg, "Python invoke of unknown function: %s", arg[0]); + snprintf(msg, 128, "Python invoke of unknown function: %s", arg[0]); error->all(FLERR, msg); } @@ -121,7 +121,7 @@ void PythonImpl::command(int narg, char **arg) str = input->variable->python_style(pfuncs[ifunc].ovarname, pfuncs[ifunc].name); if (!str) { char msg[128]; - sprintf(msg, + snprintf(msg, 128, "Python variable %s does not match variable %s " "registered with Python function %s", arg[0], pfuncs[ifunc].ovarname, pfuncs[ifunc].name); @@ -160,7 +160,7 @@ void PythonImpl::command(int narg, char **arg) err = execute_file(arg[1]); else { char msg[128]; - sprintf(msg, "Could not open python source file %s for processing", arg[1]); + snprintf(msg, 128, "Could not open python source file %s for processing", arg[1]); error->all(FLERR, msg); } } @@ -188,7 +188,7 @@ void PythonImpl::command(int narg, char **arg) ninput = utils::inumeric(FLERR, arg[iarg + 1], false, sparta); if (ninput < 0) { char msg[128]; - sprintf(msg, "Invalid number of python input arguments: %i", ninput); + snprintf(msg, 128, "Invalid number of python input arguments: %i", ninput); error->all(FLERR, msg); } iarg += 2; @@ -226,7 +226,7 @@ void PythonImpl::command(int narg, char **arg) iarg++; } else { char msg[128]; - sprintf(msg, "Unknown python command keyword: %s", arg[iarg]); + snprintf(msg, 128, "Unknown python command keyword: %s", arg[iarg]); error->all(FLERR, msg); } } @@ -255,7 +255,7 @@ void PythonImpl::command(int narg, char **arg) if (fp == nullptr) { PyUtils::Print_Errors(); char msg[128]; - sprintf(msg, "Could not open Python file: %s", pyfile); + snprintf(msg, 128, "Could not open Python file: %s", pyfile); error->all(FLERR, msg); } @@ -263,7 +263,7 @@ void PythonImpl::command(int narg, char **arg) if (err) { PyUtils::Print_Errors(); char msg[128]; - sprintf(msg, "Could not open Python file: %s", pyfile); + snprintf(msg, 128, "Could not open Python file: %s", pyfile); error->all(FLERR, msg); } fclose(fp); @@ -273,7 +273,7 @@ void PythonImpl::command(int narg, char **arg) if (err) { PyUtils::Print_Errors(); char msg[128]; - sprintf(msg, "Could not process Python string: %s", herestr); + snprintf(msg, 128, "Could not process Python string: %s", herestr); error->all(FLERR, msg); } } @@ -286,14 +286,14 @@ void PythonImpl::command(int narg, char **arg) if (!pFunc) { PyUtils::Print_Errors(); char msg[128]; - sprintf(msg, "Could not find Python function %s", pfuncs[ifunc].name); + snprintf(msg, 128, "Could not find Python function %s", pfuncs[ifunc].name); error->all(FLERR, msg); } if (!PyCallable_Check(pFunc)) { PyUtils::Print_Errors(); char msg[128]; - sprintf(msg, "Python function %s is not callable", pfuncs[ifunc].name); + snprintf(msg, 128, "Python function %s is not callable", pfuncs[ifunc].name); error->all(FLERR, msg); } @@ -323,7 +323,7 @@ void PythonImpl::invoke_function(int ifunc, char *result, double *dvalue) if (!pArgs) { char msg[128]; - sprintf(msg, "Could not prepare arguments for Python function %s", pfuncs[ifunc].name); + snprintf(msg, 128, "Could not prepare arguments for Python function %s", pfuncs[ifunc].name); error->all(FLERR, msg); } @@ -334,7 +334,7 @@ void PythonImpl::invoke_function(int ifunc, char *result, double *dvalue) str = input->variable->retrieve(pfuncs[ifunc].svalue[i]); if (!str) { char msg[128]; - sprintf(msg, "Could not evaluate Python function %s input variable: %s", + snprintf(msg, 128, "Could not evaluate Python function %s input variable: %s", pfuncs[ifunc].name, pfuncs[ifunc].svalue[i]); error->all(FLERR, msg); } @@ -350,7 +350,7 @@ void PythonImpl::invoke_function(int ifunc, char *result, double *dvalue) str = input->variable->retrieve(pfuncs[ifunc].svalue[i]); if (!str) { char msg[128]; - sprintf(msg, "Could not evaluate Python function %s input variable: %s", + snprintf(msg, 128, "Could not evaluate Python function %s input variable: %s", pfuncs[ifunc].name, pfuncs[ifunc].svalue[i]); error->all(FLERR, msg); } @@ -366,7 +366,7 @@ void PythonImpl::invoke_function(int ifunc, char *result, double *dvalue) str = input->variable->retrieve(pfuncs[ifunc].svalue[i]); if (!str) { char msg[128]; - sprintf(msg, "Could not evaluate Python function %s input variable: %s", + snprintf(msg, 128, "Could not evaluate Python function %s input variable: %s", pfuncs[ifunc].name, pfuncs[ifunc].svalue[i]); error->all(FLERR, msg); } @@ -378,7 +378,7 @@ void PythonImpl::invoke_function(int ifunc, char *result, double *dvalue) pValue = PyCapsule_New((void *)sparta, nullptr, nullptr); } else { char msg[128]; - sprintf(msg, "Unsupported variable type: %i", itype); + snprintf(msg, 128, "Unsupported variable type: %i", itype); error->all(FLERR, msg); } PyTuple_SetItem(pArgs, i, pValue); @@ -393,7 +393,7 @@ void PythonImpl::invoke_function(int ifunc, char *result, double *dvalue) if (!pValue) { PyUtils::Print_Errors(); char msg[128]; - sprintf(msg, "Python evaluation of function %s failed", pfuncs[ifunc].name); + snprintf(msg, 128, "Python evaluation of function %s failed", pfuncs[ifunc].name); error->all(FLERR, msg); } @@ -543,7 +543,7 @@ int PythonImpl::create_entry(char *name, int ninput, int noutput, error->all(FLERR, "Missing python format keyword"); else if (format && ((int) strlen(format) != ninput + noutput)) { char msg[128]; - sprintf(msg, "Input/output arguments (%i) and format characters (%zu) are inconsistent", + snprintf(msg, 128, "Input/output arguments (%i) and format characters (%zu) are inconsistent", (ninput + noutput), strlen(format)); error->all(FLERR, msg); } @@ -626,7 +626,7 @@ int PythonImpl::create_entry(char *name, int ninput, int noutput, } else { char msg[128]; - sprintf(msg, "Invalid python format character: %i", type); + snprintf(msg, 128, "Invalid python format character: %i", type); error->all(FLERR, msg); } } @@ -646,7 +646,7 @@ int PythonImpl::create_entry(char *name, int ninput, int noutput, pfuncs[ifunc].otype = STRING; else { char msg[128]; - sprintf(msg, "Invalid python return format character: %i", type); + snprintf(msg, 128, "Invalid python return format character: %i", type); error->all(FLERR, msg); } diff --git a/src/adapt_grid.cpp b/src/adapt_grid.cpp index 63859edbc..d6df6bd36 100644 --- a/src/adapt_grid.cpp +++ b/src/adapt_grid.cpp @@ -1526,6 +1526,7 @@ int AdaptGrid::perform_coarsen() // coarsen parentID to become a new child cell + int nglocalprev = grid->nlocal; grid->coarsen_cell(parentID,plevel,plo,phi,nchild, alist[i].index,alist[i].nsurf,alist[i].np, alist[i].surfs,alist[i].particles,cut2d,cut3d); @@ -1533,7 +1534,7 @@ int AdaptGrid::perform_coarsen() cells = grid->cells; cinfo = grid->cinfo; sinfo = grid->sinfo; - newcell = grid->nlocal - 1; + newcell = nglocalprev; // if new child has no surfs and any of its children was INSIDE // then type of new child cell = INSIDE diff --git a/src/collide_vss.cpp b/src/collide_vss.cpp index 41cc91ec2..735096d16 100644 --- a/src/collide_vss.cpp +++ b/src/collide_vss.cpp @@ -225,6 +225,7 @@ int CollideVSS::test_collision(int icell, int igroup, int jgroup, double vre = vro*prefactor[ispecies][jspecies]; vremax[icell][igroup][jgroup] = MAX(vre,vremax[icell][igroup][jgroup]); + if (vremax[icell][igroup][jgroup] == 0.0) return 0; if (vre/vremax[icell][igroup][jgroup] < random->uniform()) return 0; precoln.vr2 = vr2; return 1; @@ -1049,7 +1050,7 @@ void CollideVSS::read_param_file(char *fname) params[isp][jsp].alpha = params[jsp][isp].alpha = atof(words[5]); if (relaxflag == VARIABLE) { params[isp][jsp].rotc1 = params[jsp][isp].rotc1 = atof(words[6]); - params[isp][jsp].rotc2 = atof(words[7]); + params[isp][jsp].rotc2 = params[jsp][isp].rotc2 = atof(words[7]); params[isp][jsp].rotc3 = params[jsp][isp].rotc3 = (MY_PI+MY_PI2*MY_PI2)*params[isp][jsp].rotc2; if(params[isp][jsp].rotc2 > 0) diff --git a/src/compute_count.cpp b/src/compute_count.cpp index 36a0ba4c8..8880862da 100644 --- a/src/compute_count.cpp +++ b/src/compute_count.cpp @@ -47,7 +47,7 @@ ComputeCount::ComputeCount(SPARTA *sparta, int narg, char **arg) : if (imix < 0) error->all(FLERR,"Unknown species/mixture in compute count command"); int igroup = particle->mixture[imix]->find_group(ptr+1); - if (imix < 0) + if (igroup < 0) error->all(FLERR,"Unknown mixture group in compute count command"); *ptr = '/'; if (nvalues == maxvalues) allocate(maxvalues+DELTAVALUES); diff --git a/src/compute_gas_collision_tally.cpp b/src/compute_gas_collision_tally.cpp index 9d4bc661d..a3dee2bd3 100644 --- a/src/compute_gas_collision_tally.cpp +++ b/src/compute_gas_collision_tally.cpp @@ -55,7 +55,7 @@ ComputeGasCollisionTally::ComputeGasCollisionTally(SPARTA *sparta, int narg, cha else if (strcmp(arg[iarg],"id1") == 0) which[nvalue++] = ID1; else if (strcmp(arg[iarg],"id2") == 0) which[nvalue++] = ID2; else if (strcmp(arg[iarg],"type1") == 0) which[nvalue++] = TYPE1; - else if (strcmp(arg[iarg],"type2") == 0) which[nvalue++] = TYPE1; + else if (strcmp(arg[iarg],"type2") == 0) which[nvalue++] = TYPE2; else if (strcmp(arg[iarg],"vx1/pre") == 0) which[nvalue++] = VX1PRE; else if (strcmp(arg[iarg],"vy1/pre") == 0) which[nvalue++] = VY1PRE; else if (strcmp(arg[iarg],"vz1/pre") == 0) which[nvalue++] = VZ1PRE; diff --git a/src/compute_gas_reaction_tally.cpp b/src/compute_gas_reaction_tally.cpp index 058b6737e..3302c033b 100644 --- a/src/compute_gas_reaction_tally.cpp +++ b/src/compute_gas_reaction_tally.cpp @@ -68,8 +68,8 @@ ComputeGasReactionTally::ComputeGasReactionTally(SPARTA *sparta, int narg, char else if (strcmp(arg[iarg],"vy1/pre") == 0) which[nvalue++] = VY1PRE; else if (strcmp(arg[iarg],"vz1/pre") == 0) which[nvalue++] = VZ1PRE; else if (strcmp(arg[iarg],"vx2/pre") == 0) which[nvalue++] = VX2PRE; - else if (strcmp(arg[iarg],"vy2/pre") == 0) which[nvalue++] = VX2PRE; - else if (strcmp(arg[iarg],"vz2/pre") == 0) which[nvalue++] = VY2PRE; + else if (strcmp(arg[iarg],"vy2/pre") == 0) which[nvalue++] = VY2PRE; + else if (strcmp(arg[iarg],"vz2/pre") == 0) which[nvalue++] = VZ2PRE; else if (strcmp(arg[iarg],"vx1/post") == 0) which[nvalue++] = VX1POST; else if (strcmp(arg[iarg],"vy1/post") == 0) which[nvalue++] = VY1POST; else if (strcmp(arg[iarg],"vz1/post") == 0) which[nvalue++] = VZ1POST; diff --git a/src/compute_lambda_grid.cpp b/src/compute_lambda_grid.cpp index 2016249b7..d93a43abb 100644 --- a/src/compute_lambda_grid.cpp +++ b/src/compute_lambda_grid.cpp @@ -194,8 +194,10 @@ ComputeLambdaGrid::ComputeLambdaGrid(SPARTA *sparta, int narg, char **arg) : char *ptr = strchr(suffix,'['); if (ptr) { - if (suffix[strlen(suffix)-1] != ']') + if (suffix[strlen(suffix)-1] != ']') { + delete [] suffix; error->all(FLERR,"Invalid nrho in compute lambda/grid command"); + } nrhoindex[i] = atoi(ptr+1); *ptr = '\0'; } else nrhoindex[i] = 0; @@ -675,12 +677,12 @@ void ComputeLambdaGrid::compute_per_grid() if (knyflag) { if (noutputs == 1) vector_grid[i] = lambda / sizey; - array_grid[i][output_order[KNY]] = lambda / sizey; + else array_grid[i][output_order[KNY]] = lambda / sizey; } if (knzflag) { if (noutputs == 1) vector_grid[i] = lambda / sizez; - array_grid[i][output_order[KNZ]] = lambda / sizez; + else array_grid[i][output_order[KNZ]] = lambda / sizez; } } } diff --git a/src/compute_property_surf.cpp b/src/compute_property_surf.cpp index 85f0569b7..738056a90 100644 --- a/src/compute_property_surf.cpp +++ b/src/compute_property_surf.cpp @@ -247,7 +247,7 @@ void ComputePropertySurf::pack_id(int n) Surf::Tri *tris; if (distributed) tris = surf->mytris; else tris = surf->tris; - for (int i = 0; i < nsown; i++) { + for (int i = 0; i < nchoose; i++) { m = cglobal[i]; if (tris[m].mask & groupbit) buf[n] = tris[m].id; else buf[n] = 0.0; diff --git a/src/compute_react_boundary.cpp b/src/compute_react_boundary.cpp index ce41acce5..a0bb7d6da 100644 --- a/src/compute_react_boundary.cpp +++ b/src/compute_react_boundary.cpp @@ -66,7 +66,6 @@ ComputeReactBoundary(SPARTA *sparta, int narg, char **arg) : char *ptr = copy; while ((ptr = strtok(ptr,"/")) != (char *) NULL) { for (int ireaction = 0; ireaction < ntotal; ireaction++) { - reaction2col[ireaction][icol] = 0; if (which == REACTANT) { if (surf->sr[isr]->match_reactant(ptr,ireaction)) reaction2col[ireaction][icol] = 1; diff --git a/src/compute_react_isurf_grid.cpp b/src/compute_react_isurf_grid.cpp index 84aaa69a3..a7eac68a4 100644 --- a/src/compute_react_isurf_grid.cpp +++ b/src/compute_react_isurf_grid.cpp @@ -75,7 +75,6 @@ ComputeReactISurfGrid(SPARTA *sparta, int narg, char **arg) : char *ptr = copy; while ((ptr = strtok(ptr,"/")) != (char *) NULL) { for (int ireaction = 0; ireaction < ntotal; ireaction++) { - reaction2col[ireaction][icol] = 0; if (which == REACTANT) { if (surf->sr[isr]->match_reactant(ptr,ireaction)) reaction2col[ireaction][icol] = 1; diff --git a/src/compute_react_surf.cpp b/src/compute_react_surf.cpp index f218a4c47..8a9219c14 100644 --- a/src/compute_react_surf.cpp +++ b/src/compute_react_surf.cpp @@ -67,7 +67,6 @@ ComputeReactSurf::ComputeReactSurf(SPARTA *sparta, int narg, char **arg) : char *ptr = copy; while ((ptr = strtok(ptr,"/")) != (char *) NULL) { for (int ireaction = 0; ireaction < ntotal; ireaction++) { - reaction2col[ireaction][icol] = 0; if (which == REACTANT) { if (surf->sr[isr]->match_reactant(ptr,ireaction)) reaction2col[ireaction][icol] = 1; diff --git a/src/compute_reduce.cpp b/src/compute_reduce.cpp index c925ab439..e366a0bdf 100644 --- a/src/compute_reduce.cpp +++ b/src/compute_reduce.cpp @@ -133,8 +133,10 @@ ComputeReduce::ComputeReduce(SPARTA *spa, int narg, char **arg) : char *ptr = strchr(suffix,'['); if (ptr) { - if (suffix[strlen(suffix)-1] != ']') + if (suffix[strlen(suffix)-1] != ']') { + delete [] suffix; error->all(FLERR,"Illegal compute reduce command"); + } argindex[nvalues] = atoi(ptr+1); *ptr = '\0'; } else argindex[nvalues] = 0; @@ -158,7 +160,7 @@ ComputeReduce::ComputeReduce(SPARTA *spa, int narg, char **arg) : while (iarg < nargnew) { if (strcmp(arg[iarg],"replace") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal compute reduce command"); + if (iarg+3 > nargnew) error->all(FLERR,"Illegal compute reduce command"); if (mode != MINN && mode != MAXX) error->all(FLERR,"Compute reduce replace requires min or max mode"); int col1 = atoi(arg[iarg+1]) - 1; @@ -171,7 +173,7 @@ ComputeReduce::ComputeReduce(SPARTA *spa, int narg, char **arg) : replace[col1] = col2; iarg += 3; } else if (strcmp(arg[iarg],"subset") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal compute reduce command"); + if (iarg+2 > nargnew) error->all(FLERR,"Illegal compute reduce command"); int n = strlen(arg[iarg+1]) + 1; subsetID = new char[n]; strcpy(subsetID,arg[iarg+1]); diff --git a/src/compute_surf.cpp b/src/compute_surf.cpp index 3a0b5222b..c829fe229 100644 --- a/src/compute_surf.cpp +++ b/src/compute_surf.cpp @@ -614,8 +614,9 @@ void ComputeSurf::surf_tally(double /*dtremain*/, int isurf, int icell, int reac if (reaction && !transparent) { sr = surf->sr[isr]; r_coeff = sr->reaction_coeff(reaction-1); - vec[k++] += weight * r_coeff * fluxscale; + vec[k] += weight * r_coeff * fluxscale; } + k++; break; case ETOT: if (iorig) vsqpre = origmass * MathExtra::lensq3(vorig); diff --git a/src/compute_surf_reaction_tally.cpp b/src/compute_surf_reaction_tally.cpp index ef2768ceb..73713a1d4 100644 --- a/src/compute_surf_reaction_tally.cpp +++ b/src/compute_surf_reaction_tally.cpp @@ -288,7 +288,7 @@ int ComputeSurfReactionTally::datatype(int icol) if (sizeof(surfint) == sizeof(bigint)) return BIGINT; } if (which[icol-1] == IDPRE) return INT; - if (which[icol-1] == ID2POST || which[icol-1] == ID2POST) return INT; + if (which[icol-1] == ID1POST || which[icol-1] == ID2POST) return INT; if (which[icol-1] == TYPEPRE) return INT; if (which[icol-1] == TYPE1POST || which[icol-1] == TYPE2POST) return INT; diff --git a/src/compute_tvib_grid.cpp b/src/compute_tvib_grid.cpp index ab7063c39..3cad375ee 100644 --- a/src/compute_tvib_grid.cpp +++ b/src/compute_tvib_grid.cpp @@ -499,7 +499,7 @@ void ComputeTvibGrid::post_process_grid(int index, int /*nsample*/, evib = emap[2*imode]; count = evib+1; for (isp = 0; isp < nsp; isp++) { - ispecies = groupspecies[index][isp]; + ispecies = groupspecies[index / maxmode][isp]; theta = species[ispecies].vibtemp[imode]; if (theta == 0.0 || etally[icell][count] == 0.0) { tspecies_mode[isp][imode] = 0.0; diff --git a/src/domain.cpp b/src/domain.cpp index 2c7f84c06..a07fe3fe6 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -440,7 +440,7 @@ void Domain::add_region(int narg, char **arg) if (sparta->suffix_enable) { if (sparta->suffix) { char estyle[256]; - sprintf(estyle,"%s/%s",arg[1],sparta->suffix); + snprintf(estyle,256,"%s/%s",arg[1],sparta->suffix); if (0) return; diff --git a/src/dump.cpp b/src/dump.cpp index 2c99d383a..2b27c34a9 100644 --- a/src/dump.cpp +++ b/src/dump.cpp @@ -459,7 +459,7 @@ void Dump::openfile() int Dump::convert_string(int n, double *mybuf) { int i,j; - char str[32]; + char str[128]; int offset = 0; int m = 0; diff --git a/src/dump_grid.cpp b/src/dump_grid.cpp index aa267f0e3..5c6643173 100644 --- a/src/dump_grid.cpp +++ b/src/dump_grid.cpp @@ -365,7 +365,7 @@ void DumpGrid::write_string(int n, double *mybuf) void DumpGrid::write_text(int n, double *mybuf) { int i,j; - char str[32]; + char str[128]; int m = 0; for (i = 0; i < n; i++) { diff --git a/src/dump_movie.cpp b/src/dump_movie.cpp index 59113a77f..a9af140c9 100644 --- a/src/dump_movie.cpp +++ b/src/dump_movie.cpp @@ -46,7 +46,7 @@ void DumpMovie::openfile() if ((comm->me == 0) && (fp == NULL)) { #ifdef SPARTA_FFMPEG - sprintf(moviecmd,"ffmpeg -v error -y -r %.2f -f image2pipe -c:v ppm -i - " + snprintf(moviecmd,1024,"ffmpeg -v error -y -r %.2f -f image2pipe -c:v ppm -i - " "-r 24.0 -b:v %dk %s ", framerate, bitrate, filename); #else error->one(FLERR,"Support for writing movies not included"); @@ -60,7 +60,7 @@ void DumpMovie::openfile() if (fp == NULL) { char str[128]; - sprintf(str,"Failed to open FFmpeg pipeline to file %s",filename); + snprintf(str,128,"Failed to open FFmpeg pipeline to file %s",filename); error->one(FLERR,str); } } diff --git a/src/fix_ablate.cpp b/src/fix_ablate.cpp index d723cb1a2..a5510af3a 100644 --- a/src/fix_ablate.cpp +++ b/src/fix_ablate.cpp @@ -100,8 +100,10 @@ FixAblate::FixAblate(SPARTA *sparta, int narg, char **arg) : char *ptr = strchr(suffix,'['); if (ptr) { - if (suffix[strlen(suffix)-1] != ']') + if (suffix[strlen(suffix)-1] != ']') { + delete [] suffix; error->all(FLERR,"Illegal fix ablate command"); + } argindex = atoi(ptr+1); *ptr = '\0'; } else argindex = 0; diff --git a/src/fix_ablate_multi_inner.cpp b/src/fix_ablate_multi_inner.cpp index e3523cfef..b377a9258 100644 --- a/src/fix_ablate_multi_inner.cpp +++ b/src/fix_ablate_multi_inner.cpp @@ -89,7 +89,8 @@ void FixAblate::decrement_multid_outside() Ninterface = find_ninter(); total = celldelta[icell]; - perout = total/Ninterface; + if (Ninterface > 0) perout = total / Ninterface; + else perout = 0.0; // iterate to find the number of vertices around each corner // also assign perout to the interface points @@ -700,7 +701,8 @@ void FixAblate::decrement_multiv_multid_outside() Ninterface = find_ninter(); total = celldelta[icell]; - perout = total/Ninterface; + if (Ninterface > 0) perout = total / Ninterface; + else perout = 0.0; for (i = 0; i < ncorner; i++) { diff --git a/src/fix_ave_grid.cpp b/src/fix_ave_grid.cpp index f232d1377..9d5427332 100644 --- a/src/fix_ave_grid.cpp +++ b/src/fix_ave_grid.cpp @@ -103,8 +103,10 @@ FixAveGrid::FixAveGrid(SPARTA *sparta, int narg, char **arg) : char *ptr = strchr(suffix,'['); if (ptr) { - if (suffix[strlen(suffix)-1] != ']') + if (suffix[strlen(suffix)-1] != ']') { + delete [] suffix; error->all(FLERR,"Illegal fix ave/grid command"); + } argindex[i] = atoi(ptr+1); *ptr = '\0'; } else argindex[i] = 0; diff --git a/src/fix_ave_histo.cpp b/src/fix_ave_histo.cpp index debae8e6b..7d774be0e 100644 --- a/src/fix_ave_histo.cpp +++ b/src/fix_ave_histo.cpp @@ -152,8 +152,10 @@ FixAveHisto::FixAveHisto(SPARTA *spa, int narg, char **arg) : char *ptr = strchr(suffix,'['); if (ptr) { - if (suffix[strlen(suffix)-1] != ']') + if (suffix[strlen(suffix)-1] != ']') { + delete [] suffix; error->all(FLERR,"Illegal fix ave/histo command"); + } argindex[i] = atoi(ptr+1); *ptr = '\0'; } else argindex[i] = 0; diff --git a/src/fix_ave_surf.cpp b/src/fix_ave_surf.cpp index 4b681a19b..4661edaa2 100644 --- a/src/fix_ave_surf.cpp +++ b/src/fix_ave_surf.cpp @@ -107,8 +107,10 @@ FixAveSurf::FixAveSurf(SPARTA *sparta, int narg, char **arg) : char *ptr = strchr(suffix,'['); if (ptr) { - if (suffix[strlen(suffix)-1] != ']') + if (suffix[strlen(suffix)-1] != ']') { + delete [] suffix; error->all(FLERR,"Illegal fix ave/surf command"); + } argindex[i] = atoi(ptr+1); *ptr = '\0'; } else argindex[i] = 0; diff --git a/src/fix_ave_time.cpp b/src/fix_ave_time.cpp index 77ff4eb8a..0abea6fb0 100644 --- a/src/fix_ave_time.cpp +++ b/src/fix_ave_time.cpp @@ -99,8 +99,10 @@ FixAveTime::FixAveTime(SPARTA *sparta, int narg, char **arg) : char *ptr = strchr(suffix,'['); if (ptr) { - if (suffix[strlen(suffix)-1] != ']') + if (suffix[strlen(suffix)-1] != ']') { + delete [] suffix; error->all(FLERR,"Illegal fix ave/time command"); + } argindex[i] = atoi(ptr+1); *ptr = '\0'; } else argindex[i] = 0; diff --git a/src/fix_emit_face.cpp b/src/fix_emit_face.cpp index 199316a0c..12b6349d8 100644 --- a/src/fix_emit_face.cpp +++ b/src/fix_emit_face.cpp @@ -1106,7 +1106,7 @@ void FixEmitFace::subsonic_grid() tempmax = MAX(tempmax,temp_thermal_cell); } - if (np) { + if (np && massrho_cell*soundspeed_cell > 0.0) { ndim = tasks[i].ndim; sign = tasks[i].normal[ndim]; vstream[ndim] += sign * diff --git a/src/fix_emit_face_file.cpp b/src/fix_emit_face_file.cpp index 7409308ca..cb9d2af18 100644 --- a/src/fix_emit_face_file.cpp +++ b/src/fix_emit_face_file.cpp @@ -479,7 +479,7 @@ void FixEmitFaceFile::perform_task() v[ndim] = beta_un*vscale[isp]*normal[ndim] + vstream[ndim]; - theta = MY_PI * random->uniform(); + theta = MY_2PI * random->uniform(); vr = vscale[isp] * sqrt(-log(random->uniform())); v[pdim] = vr * sin(theta) + vstream[pdim]; v[qdim] = vr * cos(theta) + vstream[qdim]; @@ -1233,7 +1233,7 @@ void FixEmitFaceFile::subsonic_grid() tempmax = MAX(tempmax,temp_thermal_cell); } - if (np) { + if (np && massrho_cell*soundspeed_cell > 0.0) { sign = normal[ndim]; vstream[ndim] += sign * (tasks[i].press - press_cell) / (massrho_cell*soundspeed_cell); diff --git a/src/fix_emit_surf.cpp b/src/fix_emit_surf.cpp index a6923ff18..d64d84419 100644 --- a/src/fix_emit_surf.cpp +++ b/src/fix_emit_surf.cpp @@ -1530,7 +1530,7 @@ void FixEmitSurf::subsonic_grid() else normal = norm_vstream; } - if (np) { + if (np && massrho_cell*soundspeed_cell > 0.0) { vsmag = (psubsonic - press_cell) / (massrho_cell*soundspeed_cell); vstream[0] += vsmag*normal[0]; vstream[1] += vsmag*normal[1]; diff --git a/src/fix_grid_check.cpp b/src/fix_grid_check.cpp index f206b5f36..e4c29ce96 100644 --- a/src/fix_grid_check.cpp +++ b/src/fix_grid_check.cpp @@ -219,9 +219,9 @@ void FixGridCheck::end_of_step() if (outflag == ERROR) { char str[128]; sprintf(str, - "Particle %d,%d on proc %d at %g %g %d is inside surfs in cell " + "Particle %d,%d on proc %d at %g %g %g is inside surfs in cell " CELLINT_FORMAT " on timestep " BIGINT_FORMAT, - i,particles[i].id,comm->me,x[0],x[1],icell,cells[icell].id, + i,particles[i].id,comm->me,x[0],x[1],x[2],cells[icell].id, update->ntimestep); error->one(FLERR,str); } diff --git a/src/fix_halt.cpp b/src/fix_halt.cpp index 5c304e851..767672887 100644 --- a/src/fix_halt.cpp +++ b/src/fix_halt.cpp @@ -237,7 +237,7 @@ void FixHalt::end_of_step() // print message with ID of fix halt in case multiple instances char message[128]; - sprintf(message, "Fix halt condition for fix-id %s met on step %ld with value %g", + sprintf(message, "Fix halt condition for fix-id %s met on step " BIGINT_FORMAT " with value %g", id, update->ntimestep, attvalue); if (eflag == HARD) { error->all(FLERR, message); diff --git a/src/fix_print.cpp b/src/fix_print.cpp index 8f1fdeacd..19d0febf0 100644 --- a/src/fix_print.cpp +++ b/src/fix_print.cpp @@ -60,7 +60,7 @@ FixPrint::FixPrint(SPARTA *sparta, int narg, char **arg) : else fp = fopen(arg[iarg+1],"a"); if (fp == NULL) { char str[128]; - sprintf(str,"Cannot open fix print file %s",arg[iarg+1]); + snprintf(str,128,"Cannot open fix print file %s",arg[iarg+1]); error->one(FLERR,str); } } diff --git a/src/fix_surf_temp.cpp b/src/fix_surf_temp.cpp index cb753198e..015446782 100644 --- a/src/fix_surf_temp.cpp +++ b/src/fix_surf_temp.cpp @@ -145,7 +145,7 @@ FixSurfTemp::FixSurfTemp(SPARTA *sparta, int narg, char **arg) : } else if (strcmp(update->unit_style,"cgs") == 0) { prefactor = 1.0 / (emi * SB_CGS); threshold = 1.0e-3; - } + } else error->all(FLERR,"Fix surf/temp requires si or cgs units"); // trigger one-time initialization of custom per-surf temperatures @@ -173,6 +173,16 @@ int FixSurfTemp::setmask() void FixSurfTemp::init() { + if (source == COMPUTE) { + icompute = modify->find_compute(id_qw); + if (icompute < 0) error->all(FLERR,"Could not find fix surf/temp compute ID"); + cqw = modify->compute[icompute]; + } else if (source == FIX) { + ifix = modify->find_fix(id_qw); + if (ifix < 0) error->all(FLERR,"Could not find fix surf/temp fix ID"); + fqw = modify->fix[ifix]; + } + if (!firstflag) return; firstflag = 0; diff --git a/src/fix_temp_rescale.cpp b/src/fix_temp_rescale.cpp index 692e2315f..13b242781 100644 --- a/src/fix_temp_rescale.cpp +++ b/src/fix_temp_rescale.cpp @@ -172,7 +172,8 @@ void FixTempRescale::end_of_step_no_average(double t_target) // vscale = scale factor for thermal velocity components - vscale = sqrt(t_target/t_current); + if (t_current <= 0.0) vscale = 1.0; + else vscale = sqrt(t_target/t_current); // 2nd pass: loop over particles in cell // rescale thermal velocity components @@ -294,7 +295,9 @@ void FixTempRescale::end_of_step_average(double t_target) // scale all particles in all cells by vscale t_current /= n_current; - double vscale = sqrt(t_target/t_current); + double vscale; + if (t_current <= 0.0) vscale = 1.0; + else vscale = sqrt(t_target/t_current); // loop over grid cells to rescale velocity of particles in each // single-particle cells are also rescaled, their vcom = 0.0 diff --git a/src/geometry.cpp b/src/geometry.cpp index 329c3d7cf..7b538d036 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -763,20 +763,33 @@ bool axi_line_intersect(double tdelta, double *x, double *v, double dconst = x21*v1[1] - y21*v1[0]; double a = x21sq*(v[1]*v[1] + v[2]*v[2]) - y21sq*v[0]*v[0]; - if (a == 0.0) return false; double b = x21sq*x[1]*v[1] - y21sq*x[0]*v[0] - y21*v[0]*dconst; double c = x21sq*x[1]*x[1] - y21sq*x[0]*x[0] - 2.0*y21*x[0]*dconst - dconst*dconst; - double arg = b*b - a*c; - if (arg < 0.0) return false; - double sarg = sqrt(arg); - - nc = 2; - double tone = (-b - sarg) / a; - double ttwo = (-b + sarg) / a; - t1 = MIN(tone,ttwo); - t2 = MAX(tone,ttwo); + if (a == 0.0) { + if (b == 0.0) return false; + nc = 1; + t1 = t2 = -0.5 * c / b; + } else { + double arg = b*b - a*c; + if (arg < 0.0) return false; + double sarg = sqrt(arg); + nc = 2; + double tone, ttwo; + if (b > 0.0) { + tone = (-b - sarg) / a; + ttwo = c / (-b - sarg); + } else if (b < 0.0) { + tone = c / (-b + sarg); + ttwo = (-b + sarg) / a; + } else { + tone = -sarg / a; + ttwo = sarg / a; + } + t1 = MIN(tone,ttwo); + t2 = MAX(tone,ttwo); + } } // if selfflag, particle starts on surf line segment @@ -839,11 +852,16 @@ bool axi_line_intersect(double tdelta, double *x, double *v, if (v1[1] == v2[1]) xc[1] = v1[1]; xc[2] = 0.0; - double rn = ynew / xc[1]; - double wn = znew / xc[1]; vc[0] = v[0]; - vc[1] = v[1]*rn + v[2]*wn; - vc[2] = -v[1]*wn + v[2]*rn; + if (xc[1] > 0.0) { + double rn = ynew / xc[1]; + double wn = znew / xc[1]; + vc[1] = v[1]*rn + v[2]*wn; + vc[2] = -v[1]*wn + v[2]*rn; + } else { + vc[1] = v[1]; + vc[2] = v[2]; + } // test that xc is within line segment bounds // y-test for vertical line, else x-test @@ -929,10 +947,20 @@ bool axi_horizontal_line(double tdelta, double *x, double *v, double arg = yhoriz*yhoriz*a - v[2]*v[2]*x[1]*x[1]; if (arg < 0.0) return false; double sarg = sqrt(arg); + double c = x[1]*x[1] - yhoriz*yhoriz; nc = 2; - double tone = (b - sarg) / a; - double ttwo = (b + sarg) / a; + double tone, ttwo; + if (b > 0.0) { + ttwo = (b + sarg) / a; + tone = c / (b + sarg); + } else if (b < 0.0) { + tone = (b - sarg) / a; + ttwo = c / (b - sarg); + } else { + tone = -sarg / a; + ttwo = sarg / a; + } t1 = MIN(tone,ttwo); t2 = MAX(tone,ttwo); diff --git a/src/grid.cpp b/src/grid.cpp index 435fe3a8c..0c3c722fb 100644 --- a/src/grid.cpp +++ b/src/grid.cpp @@ -2459,7 +2459,10 @@ void Grid::group(int narg, char **arg) int jgroup; for (int iarg = 3; iarg < narg; iarg++) { jgroup = find_group(arg[iarg]); - if (jgroup == -1) error->all(FLERR,"Group ID does not exist"); + if (jgroup == -1) { + delete [] list; + error->all(FLERR,"Group ID does not exist"); + } list[iarg-3] = jgroup; } @@ -2494,7 +2497,10 @@ void Grid::group(int narg, char **arg) int jgroup; for (int iarg = 3; iarg < narg; iarg++) { jgroup = find_group(arg[iarg]); - if (jgroup == -1) error->all(FLERR,"Group ID does not exist"); + if (jgroup == -1) { + delete [] list; + error->all(FLERR,"Group ID does not exist"); + } list[iarg-3] = jgroup; } @@ -2521,7 +2527,10 @@ void Grid::group(int narg, char **arg) int jgroup; for (int iarg = 3; iarg < narg; iarg++) { jgroup = find_group(arg[iarg]); - if (jgroup == -1) error->all(FLERR,"Group ID does not exist"); + if (jgroup == -1) { + delete [] list; + error->all(FLERR,"Group ID does not exist"); + } list[iarg-3] = jgroup; } diff --git a/src/grid_id.cpp b/src/grid_id.cpp index 48ff66203..55a7f03a5 100644 --- a/src/grid_id.cpp +++ b/src/grid_id.cpp @@ -541,7 +541,8 @@ void Grid::id_num2str(cellint id, char *str) newbits = plevels[level].newbits; mask = (1L << newbits) - 1; ichild = id & mask; - sprintf(&str[offset],CELLINT_FORMAT,ichild); + if (offset >= 110) break; // prevent buffer overflow + snprintf(&str[offset], 128 - offset, CELLINT_FORMAT, ichild); offset = strlen(str); id = id >> newbits; if (!id) return; diff --git a/src/input.cpp b/src/input.cpp index d8ed9295a..a1fd48c1f 100644 --- a/src/input.cpp +++ b/src/input.cpp @@ -242,8 +242,8 @@ void Input::file() // execute the command if (execute_command()) { - char *str = new char[maxline+32]; - sprintf(str,"Unknown command: %s",line); + char str[128]; + snprintf(str,128,"Unknown command: %.90s",line); error->all(FLERR,str); } } @@ -268,7 +268,7 @@ void Input::file(const char *filename) infile = fopen(filename,"r"); if (infile == NULL) { char str[128]; - sprintf(str,"Cannot open input script %s",filename); + snprintf(str,128,"Cannot open input script %s",filename); error->one(FLERR,str); } infiles[0] = infile; @@ -309,8 +309,8 @@ char *Input::one(const char *single) // execute the command and return its name if (execute_command()) { - char *str = new char[maxline+32]; - sprintf(str,"Unknown command: %s",line); + char str[128]; + snprintf(str,128,"Unknown command: %.90s",line); error->all(FLERR,str); } @@ -1010,7 +1010,7 @@ void Input::ifthenelse() ncommands = 0; for (int i = first; i <= last; i++) { int n = strlen(arg[i]) + 1; - if (n == 1) error->all(FLERR,"Illegal if command"); + if (n == 1) { for (int j = 0; j < ncommands; j++) delete [] commands[j]; delete [] commands; error->all(FLERR,"Illegal if command"); } commands[ncommands] = new char[n]; strcpy(commands[ncommands],arg[i]); ncommands++; @@ -1064,7 +1064,7 @@ void Input::ifthenelse() ncommands = 0; for (int i = first; i <= last; i++) { int n = strlen(arg[i]) + 1; - if (n == 1) error->all(FLERR,"Illegal if command"); + if (n == 1) { for (int j = 0; j < ncommands; j++) delete [] commands[j]; delete [] commands; error->all(FLERR,"Illegal if command"); } commands[ncommands] = new char[n]; strcpy(commands[ncommands],arg[i]); ncommands++; @@ -1107,7 +1107,7 @@ void Input::include() infile = fopen(arg[0],"r"); if (infile == NULL) { char str[128]; - sprintf(str,"Cannot open input script %s",arg[0]); + snprintf(str,128,"Cannot open input script %s",arg[0]); error->one(FLERR,str); } infiles[nfile++] = infile; @@ -1132,7 +1132,7 @@ void Input::jump() infile = fopen(arg[0],"r"); if (infile == NULL) { char str[128]; - sprintf(str,"Cannot open input script %s",arg[0]); + snprintf(str,128,"Cannot open input script %s",arg[0]); error->one(FLERR,str); } infiles[nfile-1] = infile; @@ -1177,7 +1177,7 @@ void Input::log() if (logfile == NULL) { char str[128]; - sprintf(str,"Cannot open logfile %s",arg[0]); + snprintf(str,128,"Cannot open logfile %s",arg[0]); error->one(FLERR,str); } } @@ -1255,7 +1255,7 @@ void Input::print() else fp = fopen(arg[iarg+1],"a"); if (fp == NULL) { char str[128]; - sprintf(str,"Cannot open print file %s",arg[iarg+1]); + snprintf(str,128,"Cannot open print file %s",arg[iarg+1]); error->one(FLERR,str); } } diff --git a/src/marching_cubes.h b/src/marching_cubes.h index 2957d4c3a..388c41f3e 100644 --- a/src/marching_cubes.h +++ b/src/marching_cubes.h @@ -33,7 +33,7 @@ class MarchingCubes : protected Pointers { double thresh; double *lo,*hi; - int v000,v001,v010,v011,v100,v101,v110,v111; + double v000,v001,v010,v011,v100,v101,v110,v111; double v000iso,v001iso,v010iso,v011iso,v100iso,v101iso,v110iso,v111iso; double inval[8][6]; double i0,i1,i2,i3,i4,i5,i6,i7,i8,i9,i10,i11; diff --git a/src/mixture.cpp b/src/mixture.cpp index e45b963f9..b2d0a7e1e 100644 --- a/src/mixture.cpp +++ b/src/mixture.cpp @@ -222,7 +222,7 @@ void Mixture::init() if (err) { char str[128]; - sprintf(str,"Mixture %s fractions exceed 1.0",id); + snprintf(str,128,"Mixture %s fractions exceed 1.0",id); error->all(FLERR,str); } diff --git a/src/modify.cpp b/src/modify.cpp index 78daa6040..ba5062398 100644 --- a/src/modify.cpp +++ b/src/modify.cpp @@ -342,7 +342,7 @@ void Modify::add_fix(int narg, char **arg) if (sparta->suffix_enable) { if (sparta->suffix) { char estyle[256]; - sprintf(estyle,"%s/%s",arg[1],sparta->suffix); + snprintf(estyle,256,"%s/%s",arg[1],sparta->suffix); if (0) return; @@ -438,7 +438,7 @@ void Modify::add_compute(int narg, char **arg) if (sparta->suffix_enable) { if (sparta->suffix) { char estyle[256]; - sprintf(estyle,"%s/%s",arg[1],sparta->suffix); + snprintf(estyle,256,"%s/%s",arg[1],sparta->suffix); if (0) return; diff --git a/src/move_surf.cpp b/src/move_surf.cpp index 742107050..3cb1433d0 100644 --- a/src/move_surf.cpp +++ b/src/move_surf.cpp @@ -327,7 +327,7 @@ void MoveSurf::readfile() fp = fopen(file,"r"); if (fp == NULL) { char str[128]; - sprintf(str,"Cannot open move surf file %s",file); + snprintf(str,128,"Cannot open move surf file %s",file); error->one(FLERR,str); } } diff --git a/src/output.cpp b/src/output.cpp index 08176260b..88410ea1e 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -300,7 +300,7 @@ void Output::write(bigint ntimestep) char *file = new char[strlen(restart1) + 16]; char *ptr = strchr(restart1,'*'); *ptr = '\0'; - sprintf(file,"%s" BIGINT_FORMAT "%s",restart1,ntimestep,ptr+1); + snprintf(file,strlen(restart1)+16,"%s" BIGINT_FORMAT "%s",restart1,ntimestep,ptr+1); *ptr = '*'; if (last_restart != ntimestep) restart->write(file); delete [] file; @@ -386,7 +386,7 @@ void Output::write_restart(bigint ntimestep) char *file = new char[strlen(restart1) + 16]; char *ptr = strchr(restart1,'*'); *ptr = '\0'; - sprintf(file,"%s" BIGINT_FORMAT "%s",restart1,ntimestep,ptr+1); + snprintf(file,strlen(restart1)+16,"%s" BIGINT_FORMAT "%s",restart1,ntimestep,ptr+1); *ptr = '*'; restart->write(file); delete [] file; diff --git a/src/particle.cpp b/src/particle.cpp index 66229eff0..3556d5afc 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -787,7 +787,7 @@ void Particle::add_species(int narg, char **arg) fp = fopen(arg[0],"r"); if (fp == NULL) { char str[128]; - sprintf(str,"Cannot open species file %s",arg[0]); + snprintf(str,128,"Cannot open species file %s",arg[0]); error->one(FLERR,str); } } @@ -920,7 +920,7 @@ void Particle::add_species(int narg, char **arg) fp = fopen(arg[rotindex],"r"); if (fp == NULL) { char str[128]; - sprintf(str,"Cannot open rotation file %s",arg[rotindex]); + snprintf(str,128,"Cannot open rotation file %s",arg[rotindex]); error->one(FLERR,str); } } @@ -969,7 +969,7 @@ void Particle::add_species(int narg, char **arg) fp = fopen(arg[vibindex],"r"); if (fp == NULL) { char str[128]; - sprintf(str,"Cannot open vibration file %s",arg[vibindex]); + snprintf(str,128,"Cannot open vibration file %s",arg[vibindex]); error->one(FLERR,str); } } diff --git a/src/react_bird.cpp b/src/react_bird.cpp index 63dd5a881..b2728600c 100644 --- a/src/react_bird.cpp +++ b/src/react_bird.cpp @@ -71,6 +71,8 @@ ReactBird::ReactBird(SPARTA *sparta) : React(sparta) reactions = NULL; list_ij = NULL; sp2recomb_ij = NULL; + tally_reactions = NULL; + tally_reactions_all = NULL; } /* ---------------------------------------------------------------------- */ diff --git a/src/react_qk.cpp b/src/react_qk.cpp index 5c22554e4..442af603a 100644 --- a/src/react_qk.cpp +++ b/src/react_qk.cpp @@ -126,8 +126,9 @@ int ReactQK::attempt(Particle::OnePart *ip, Particle::OnePart *jp, do { iv = static_cast (random->uniform()*(maxlev+0.99999999)); evib = static_cast (iv / inverse_kT); - if (evib < ecc) react_prob = pow(1.0-evib/ecc,1.5-omega); - } while (random->uniform() < react_prob); + if (evib < ecc) prob = pow(1.0-evib/ecc,1.5-omega); + else prob = 0.0; + } while (random->uniform() < prob); ilevel = static_cast (fabs(r->coeff[4]) * inverse_kT); if (iv >= ilevel) react_prob = 1.0; @@ -157,6 +158,7 @@ int ReactQK::attempt(Particle::OnePart *ip, Particle::OnePart *jp, evib = static_cast (iv * update->boltz*species[mspec].vibtemp[0]); if (evib < ecc) prob = pow(1.0-evib/ecc,1.5 - r->coeff[6]); + else prob = 0.0; } while (random->uniform() < prob); ilevel = static_cast diff --git a/src/react_tce.cpp b/src/react_tce.cpp index 07391996b..51b42059e 100644 --- a/src/react_tce.cpp +++ b/src/react_tce.cpp @@ -179,14 +179,14 @@ int ReactTCE::attempt(Particle::OnePart *ip, Particle::OnePart *jp, break; } - if (react_prob < 0) error->warning(FLERR,"Negative reaction probability"); - else if (react_prob > 1) error->warning(FLERR,"Reaction probability greater than 1"); - default: error->one(FLERR,"Unknown outcome in reaction"); break; } + if (react_prob < 0) error->warning(FLERR,"Negative reaction probability"); + else if (react_prob > 1) error->warning(FLERR,"Reaction probability greater than 1"); + // test against random number to see if this reaction occurs // if it does, reset species of I,J and optional K to product species // J particle is destroyed in recombination reaction, set species = -1 @@ -236,6 +236,8 @@ int ReactTCE::attempt(Particle::OnePart *ip, Particle::OnePart *jp, return list[i] + 1; } + + break; } } diff --git a/src/react_tce_qk.cpp b/src/react_tce_qk.cpp index 36105ae9a..c088d41be 100644 --- a/src/react_tce_qk.cpp +++ b/src/react_tce_qk.cpp @@ -229,8 +229,9 @@ int ReactTCEQK::attempt_qk(Particle::OnePart *ip, Particle::OnePart *jp, do { iv = static_cast (random->uniform()*(maxlev+0.99999999)); evib = static_cast (iv / inverse_kT); - if (evib < ecc) react_prob = pow(1.0-evib/ecc,1.5-omega); - } while (random->uniform() < react_prob); + if (evib < ecc) prob = pow(1.0-evib/ecc,1.5-omega); + else prob = 0.0; + } while (random->uniform() < prob); ilevel = static_cast (fabs(r->coeff[4]) * inverse_kT); if (iv >= ilevel) react_prob = 1.0; @@ -260,6 +261,7 @@ int ReactTCEQK::attempt_qk(Particle::OnePart *ip, Particle::OnePart *jp, evib = static_cast (iv * update->boltz*species[mspec].vibtemp[0]); if (evib < ecc) prob = pow(1.0-evib/ecc,1.5 - r->coeff[6]); + else prob = 0.0; } while (random->uniform() < prob); ilevel = static_cast diff --git a/src/read_grid.cpp b/src/read_grid.cpp index 5e8c751b5..b7d1be74a 100644 --- a/src/read_grid.cpp +++ b/src/read_grid.cpp @@ -403,7 +403,7 @@ void ReadGrid::open(char *file) else { #ifdef SPARTA_GZIP char gunzip[128]; - sprintf(gunzip,"gunzip -c %s",file); + snprintf(gunzip,128,"gunzip -c %s",file); fp = popen(gunzip,"r"); #else error->one(FLERR,"Cannot open gzipped file"); @@ -412,7 +412,7 @@ void ReadGrid::open(char *file) if (fp == NULL) { char str[128]; - sprintf(str,"Cannot open file %s",file); + snprintf(str,128,"Cannot open file %s",file); error->one(FLERR,str); } } diff --git a/src/read_isurf.cpp b/src/read_isurf.cpp index 770ef5ead..05d778bda 100644 --- a/src/read_isurf.cpp +++ b/src/read_isurf.cpp @@ -611,6 +611,12 @@ void ReadISurf::read_corners_parallel(char *gridfile) else if (precision == DOUBLE) memory->create(dbuf,nvalues,"readisurf:dbuf"); fp = fopen(gridfile,"rb"); + if (fp == NULL) { + char str[128]; + snprintf(str,128,"Cannot open read_isurf grid corner point file %s", + gridfile); + error->one(FLERR,str); + } if (precision == INT) { fseek(fp,offset*sizeof(uint8_t)+dim*sizeof(int),SEEK_SET); tmp = fread(ibuf,sizeof(uint8_t),nvalues,fp); diff --git a/src/read_restart.cpp b/src/read_restart.cpp index 782ee0f73..c77c73458 100644 --- a/src/read_restart.cpp +++ b/src/read_restart.cpp @@ -110,7 +110,7 @@ void ReadRestart::command(int narg, char **arg) fp = fopen(hfile,"rb"); if (fp == NULL) { char str[128]; - sprintf(str,"Cannot open restart file %s",hfile); + snprintf(str,128,"Cannot open restart file %s",hfile); error->one(FLERR,str); } if (multiproc) delete [] hfile; @@ -464,7 +464,7 @@ void ReadRestart::file_search(char *infile, char *outfile) ptr = strchr(infile,'*'); *ptr = '\0'; - sprintf(outfile,"%s" BIGINT_FORMAT "%s",infile,maxnum,ptr+1); + snprintf(outfile,strlen(infile)+16,"%s" BIGINT_FORMAT "%s",infile,maxnum,ptr+1); *ptr = '*'; // clean up @@ -905,7 +905,7 @@ void ReadRestart::read_gp_multi_file_less_procs(char *file) fp = fopen(procfile,"rb"); if (fp == NULL) { char str[128]; - sprintf(str,"Cannot open restart file %s",procfile); + snprintf(str,128,"Cannot open restart file %s",procfile); error->one(FLERR,str); } @@ -991,7 +991,7 @@ void ReadRestart::read_gp_multi_file_more_procs(char *file) fp = fopen(procfile,"rb"); if (fp == NULL) { char str[128]; - sprintf(str,"Cannot open restart file %s",procfile); + snprintf(str,128,"Cannot open restart file %s",procfile); error->one(FLERR,str); } delete [] procfile; @@ -1093,7 +1093,7 @@ void ReadRestart::read_gp_multi_file_less_procs_memlimit(char *file) fp = fopen(procfile,"rb"); if (fp == NULL) { char str[128]; - sprintf(str,"Cannot open restart file %s",procfile); + snprintf(str,128,"Cannot open restart file %s",procfile); error->one(FLERR,str); } @@ -1242,7 +1242,7 @@ void ReadRestart::read_gp_multi_file_more_procs_memlimit(char *file) fp = fopen(procfile,"rb"); if (fp == NULL) { char str[128]; - sprintf(str,"Cannot open restart file %s",procfile); + snprintf(str,128,"Cannot open restart file %s",procfile); error->one(FLERR,str); } delete [] procfile; @@ -1641,7 +1641,7 @@ void ReadRestart::read_surfs_multi_file_less_procs(char *file) fp = fopen(procfile,"rb"); if (fp == NULL) { char str[128]; - sprintf(str,"Cannot open restart file %s",procfile); + snprintf(str,128,"Cannot open restart file %s",procfile); error->one(FLERR,str); } @@ -1734,7 +1734,7 @@ void ReadRestart::read_surfs_multi_file_more_procs(char *file) fp = fopen(procfile,"rb"); if (fp == NULL) { char str[128]; - sprintf(str,"Cannot open restart file %s",procfile); + snprintf(str,128,"Cannot open restart file %s",procfile); error->one(FLERR,str); } delete [] procfile; diff --git a/src/read_surf.cpp b/src/read_surf.cpp index c9d8e59b9..405a9dfa5 100644 --- a/src/read_surf.cpp +++ b/src/read_surf.cpp @@ -668,7 +668,7 @@ void ReadSurf::base(char *file) fp = fopen(hfile,"r"); if (fp == NULL) { char str[128]; - sprintf(str,"Cannot open surface base file %s",hfile); + snprintf(str,128,"Cannot open surface base file %s",hfile); error->one(FLERR,str); } delete [] hfile; @@ -2251,7 +2251,7 @@ void ReadSurf::open(char *file) else { #ifdef SPARTA_GZIP char gunzip[128]; - sprintf(gunzip,"gunzip -c %s",file); + snprintf(gunzip,128,"gunzip -c %s",file); fp = popen(gunzip,"r"); #else error->one(FLERR,"Cannot open gzipped file"); @@ -2260,7 +2260,7 @@ void ReadSurf::open(char *file) if (fp == NULL) { char str[128]; - sprintf(str,"Cannot open file %s",file); + snprintf(str,128,"Cannot open file %s",file); error->one(FLERR,str); } } @@ -2341,7 +2341,7 @@ void ReadSurf::file_search(char *infile, char *outfile) ptr = strchr(infile,'*'); *ptr = '\0'; - sprintf(outfile,"%s" BIGINT_FORMAT "%s",infile,maxnum,ptr+1); + snprintf(outfile,strlen(infile)+16,"%s" BIGINT_FORMAT "%s",infile,maxnum,ptr+1); *ptr = '*'; // clean up diff --git a/src/sparta.cpp b/src/sparta.cpp index 8027ab462..3521bfbd3 100644 --- a/src/sparta.cpp +++ b/src/sparta.cpp @@ -258,7 +258,7 @@ SPARTA::SPARTA(int narg, char **arg, MPI_Comm communicator) else infile = fopen(arg[inflag],"r"); if (infile == NULL) { char str[128]; - sprintf(str,"Cannot open input script %s",arg[inflag]); + snprintf(str,128,"Cannot open input script %s",arg[inflag]); error->one(FLERR,str); } } @@ -289,7 +289,7 @@ SPARTA::SPARTA(int narg, char **arg, MPI_Comm communicator) screen = NULL; else { char str[128]; - sprintf(str,"%s.%d",arg[screenflag],universe->iworld); + snprintf(str,128,"%s.%d",arg[screenflag],universe->iworld); screen = fopen(str,"w"); if (screen == NULL) error->one(FLERR,"Cannot open screen file"); } @@ -297,7 +297,7 @@ SPARTA::SPARTA(int narg, char **arg, MPI_Comm communicator) screen = NULL; else { char str[128]; - sprintf(str,"%s.%d",arg[partscreenflag],universe->iworld); + snprintf(str,128,"%s.%d",arg[partscreenflag],universe->iworld); screen = fopen(str,"w"); if (screen == NULL) error->one(FLERR,"Cannot open screen file"); } else screen = NULL; @@ -313,7 +313,7 @@ SPARTA::SPARTA(int narg, char **arg, MPI_Comm communicator) logfile = NULL; else { char str[128]; - sprintf(str,"%s.%d",arg[logflag],universe->iworld); + snprintf(str,128,"%s.%d",arg[logflag],universe->iworld); logfile = fopen(str,"w"); if (logfile == NULL) error->one(FLERR,"Cannot open logfile"); } @@ -321,7 +321,7 @@ SPARTA::SPARTA(int narg, char **arg, MPI_Comm communicator) logfile = NULL; else { char str[128]; - sprintf(str,"%s.%d",arg[partlogflag],universe->iworld); + snprintf(str,128,"%s.%d",arg[partlogflag],universe->iworld); logfile = fopen(str,"w"); if (logfile == NULL) error->one(FLERR,"Cannot open logfile"); } else logfile = NULL; @@ -330,7 +330,7 @@ SPARTA::SPARTA(int narg, char **arg, MPI_Comm communicator) infile = fopen(arg[inflag],"r"); if (infile == NULL) { char str[128]; - sprintf(str,"Cannot open input script %s",arg[inflag]); + snprintf(str,128,"Cannot open input script %s",arg[inflag]); error->one(FLERR,str); } } else infile = NULL; diff --git a/src/stats.cpp b/src/stats.cpp index 385000ac0..79e8df09c 100644 --- a/src/stats.cpp +++ b/src/stats.cpp @@ -168,9 +168,10 @@ void Stats::init() } n = strlen(format[i]); - sprintf(&format[i][n],"%s ",ptr); + snprintf(&format[i][n],32-n,"%s ",ptr); } - strcat(format[nfield-1],"\n"); + n = strlen(format[nfield-1]); + snprintf(&format[nfield-1][n],32-n,"\n"); delete [] format_line; @@ -222,9 +223,11 @@ void Stats::init() void Stats::header() { int loc = 0; - for (int i = 0; i < nfield; i++) - loc += sprintf(&line[loc],"%s ",keyword[i]); - sprintf(&line[loc],"\n"); + for (int i = 0; i < nfield; i++) { + loc += snprintf(&line[loc],MAXLINE-loc,"%s ",keyword[i]); + if (loc > MAXLINE-2) loc = MAXLINE-2; + } + snprintf(&line[loc],MAXLINE-loc,"\n"); if (me == 0) { if (screen) fprintf(screen,"%s",line); @@ -280,15 +283,18 @@ void Stats::compute(int flag) (this->*vfunc[ifield])(); cache_keyword[ifield] = keyword[ifield]; if (vtype[ifield] == FLOAT) { - loc += sprintf(&line[loc],format[ifield],dvalue); + loc += snprintf(&line[loc],MAXLINE-loc,format[ifield],dvalue); + if (loc > MAXLINE-1) loc = MAXLINE-1; cache_type[ifield] = 2; // SPARTA_DOUBLE cache_data[ifield].d = dvalue; } else if (vtype[ifield] == INT) { - loc += sprintf(&line[loc],format[ifield],ivalue); + loc += snprintf(&line[loc],MAXLINE-loc,format[ifield],ivalue); + if (loc > MAXLINE-1) loc = MAXLINE-1; cache_type[ifield] = 0; // SPARTA_INT cache_data[ifield].i = ivalue; } else if (vtype[ifield] == BIGINT) { - loc += sprintf(&line[loc],format[ifield],bivalue); + loc += snprintf(&line[loc],MAXLINE-loc,format[ifield],bivalue); + if (loc > MAXLINE-1) loc = MAXLINE-1; cache_type[ifield] = (sizeof(bigint) == 8) ? 4 : 0; // SPARTA_INT64 cache_data[ifield].b = bivalue; } diff --git a/src/surf.cpp b/src/surf.cpp index 5a2c1cf35..fed157ca7 100644 --- a/src/surf.cpp +++ b/src/surf.cpp @@ -2030,7 +2030,7 @@ void Surf::add_collide(int narg, char **arg) if (sparta->suffix_enable) { if (sparta->suffix) { char estyle[256]; - sprintf(estyle,"%s/%s",arg[1],sparta->suffix); + snprintf(estyle,256,"%s/%s",arg[1],sparta->suffix); if (0) return; @@ -2103,7 +2103,7 @@ void Surf::add_react(int narg, char **arg) if (sparta->suffix_enable) { if (sparta->suffix) { char estyle[256]; - sprintf(estyle,"%s/%s",arg[1],sparta->suffix); + snprintf(estyle,256,"%s/%s",arg[1],sparta->suffix); if (0) return; diff --git a/src/surf_collide_specular.cpp b/src/surf_collide_specular.cpp index 559ee5e07..602968d7c 100644 --- a/src/surf_collide_specular.cpp +++ b/src/surf_collide_specular.cpp @@ -166,7 +166,8 @@ void SurfCollideSpecular::wrapper(Particle::OnePart *p, double *norm, if (flags) noslip_flag = flags[0]; - MathExtra::reflect3(p->v,norm); + if (noslip_flag) MathExtra::negate3(p->v); + else MathExtra::reflect3(p->v,norm); } /* ---------------------------------------------------------------------- diff --git a/src/surf_react_prob.cpp b/src/surf_react_prob.cpp index 38f7e4cd2..7f78a48fe 100644 --- a/src/surf_react_prob.cpp +++ b/src/surf_react_prob.cpp @@ -308,7 +308,7 @@ void SurfReactProb::readfile(char *fname) fp = fopen(fname,"r"); if (fp == NULL) { char str[128]; - sprintf(str,"Cannot open reaction file %s",fname); + snprintf(str,128,"Cannot open reaction file %s",fname); error->one(FLERR,str); } } diff --git a/src/update.cpp b/src/update.cpp index b3f443045..f210546d3 100644 --- a/src/update.cpp +++ b/src/update.cpp @@ -938,22 +938,34 @@ template < int DIM, int SURF, int OPT > void Update::move() frac = 1.0; if (xnew[0] < lo[0]) { - frac = (lo[0]-x[0]) / (xnew[0]-x[0]); + if (xnew[0] != x[0]) frac = (lo[0]-x[0]) / (xnew[0]-x[0]); + else frac = 0.0; + if (frac < 0.0) frac = 0.0; + else if (frac > 1.0) frac = 1.0; outface = XLO; } else if (xnew[0] >= hi[0]) { - frac = (hi[0]-x[0]) / (xnew[0]-x[0]); + if (xnew[0] != x[0]) frac = (hi[0]-x[0]) / (xnew[0]-x[0]); + else frac = 0.0; + if (frac < 0.0) frac = 0.0; + else if (frac > 1.0) frac = 1.0; outface = XHI; } if (DIM != 1) { if (xnew[1] < lo[1]) { - newfrac = (lo[1]-x[1]) / (xnew[1]-x[1]); + if (xnew[1] != x[1]) newfrac = (lo[1]-x[1]) / (xnew[1]-x[1]); + else newfrac = 0.0; + if (newfrac < 0.0) newfrac = 0.0; + else if (newfrac > 1.0) newfrac = 1.0; if (newfrac < frac) { frac = newfrac; outface = YLO; } } else if (xnew[1] >= hi[1]) { - newfrac = (hi[1]-x[1]) / (xnew[1]-x[1]); + if (xnew[1] != x[1]) newfrac = (hi[1]-x[1]) / (xnew[1]-x[1]); + else newfrac = 0.0; + if (newfrac < 0.0) newfrac = 0.0; + else if (newfrac > 1.0) newfrac = 1.0; if (newfrac < frac) { frac = newfrac; outface = YHI; @@ -996,13 +1008,19 @@ template < int DIM, int SURF, int OPT > void Update::move() if (DIM == 3) { if (xnew[2] < lo[2]) { - newfrac = (lo[2]-x[2]) / (xnew[2]-x[2]); + if (xnew[2] != x[2]) newfrac = (lo[2]-x[2]) / (xnew[2]-x[2]); + else newfrac = 0.0; + if (newfrac < 0.0) newfrac = 0.0; + else if (newfrac > 1.0) newfrac = 1.0; if (newfrac < frac) { frac = newfrac; outface = ZLO; } } else if (xnew[2] >= hi[2]) { - newfrac = (hi[2]-x[2]) / (xnew[2]-x[2]); + if (xnew[2] != x[2]) newfrac = (hi[2]-x[2]) / (xnew[2]-x[2]); + else newfrac = 0.0; + if (newfrac < 0.0) newfrac = 0.0; + else if (newfrac > 1.0) newfrac = 1.0; if (newfrac < frac) { frac = newfrac; outface = ZHI; @@ -1241,7 +1259,7 @@ template < int DIM, int SURF, int OPT > void Update::move() // stuck_iterate = consecutive iterations particle is immobile - if (minparam == 0.0) stuck_iterate++; + if (minparam <= 1.0e-14) stuck_iterate++; else stuck_iterate = 0; // reset post-bounce xnew @@ -1895,6 +1913,7 @@ int Update::tally_setup() delete [] blist_active; glist_compute = slist_compute = blist_compute = NULL; + glist_active = slist_active = blist_active = NULL; nglist_compute = nslist_compute = nblist_compute = 0; for (int i = 0; i < modify->ncompute; i++) { diff --git a/src/update.h b/src/update.h index e1bf529fb..3fff20675 100644 --- a/src/update.h +++ b/src/update.h @@ -191,12 +191,14 @@ class Update : protected Pointers { double znew = x[2]; x[1] = sqrt(ynew*ynew + znew*znew); x[2] = 0.0; - double rn = ynew / x[1]; - double wn = znew / x[1]; - double vy = v[1]; - double vz = v[2]; - v[1] = vy*rn + vz*wn; - v[2] = -vy*wn + vz*rn; + if (x[1] > 0.0) { + double rn = ynew / x[1]; + double wn = znew / x[1]; + double vy = v[1]; + double vz = v[2]; + v[1] = vy*rn + vz*wn; + v[2] = -vy*wn + vz*rn; + } }; typedef void (Update::*FnPtr)(); diff --git a/src/utils.cpp b/src/utils.cpp index 3e2f62ce5..a1a156388 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -123,7 +123,7 @@ void utils::missing_cmd_args(const std::string &file, int line, const std::strin Error *error) { char msg[128]; - sprintf(msg,"Illegal %s command: missing argument(s)",cmd.c_str()); + snprintf(msg,sizeof(msg),"Illegal %s command: missing argument(s)",cmd.c_str()); if (error) error->all(file.c_str(), line, msg); } diff --git a/src/variable.cpp b/src/variable.cpp index 4bf1a00f9..6024ef7c0 100644 --- a/src/variable.cpp +++ b/src/variable.cpp @@ -654,10 +654,13 @@ int Variable::next(int narg, char **arg) delete random; FILE *fp = fopen("tmp.sparta.variable.lock","r"); + if (fp == NULL) error->one(FLERR,"Could not open variable lock file for reading"); int tmp = fscanf(fp,"%d",&nextindex); + if (tmp != 1) error->one(FLERR,"Failed to read index from variable lock file"); //printf("READ %d %d\n",universe->me,nextindex); fclose(fp); fp = fopen("tmp.sparta.variable.lock","w"); + if (fp == NULL) error->one(FLERR,"Could not open variable lock file for writing"); fprintf(fp,"%d\n",nextindex+1); //printf("WRITE %d %d\n",universe->me,nextindex+1); fclose(fp); @@ -1325,8 +1328,10 @@ double Variable::evaluate(char *str, Tree **tree) strcpy(id,&word[2]); int icompute = modify->find_compute(id); - if (icompute < 0) + if (icompute < 0) { + delete [] id; error->all(FLERR,"Invalid compute ID in variable formula"); + } Compute *compute = modify->compute[icompute]; delete [] id; @@ -1602,7 +1607,10 @@ double Variable::evaluate(char *str, Tree **tree) strcpy(id,&word[2]); int ifix = modify->find_fix(id); - if (ifix < 0) error->all(FLERR,"Invalid fix ID in variable formula"); + if (ifix < 0) { + delete [] id; + error->all(FLERR,"Invalid fix ID in variable formula"); + } Fix *fix = modify->fix[ifix]; delete [] id; @@ -1825,30 +1833,42 @@ double Variable::evaluate(char *str, Tree **tree) int icustom,size,type; if (cwhich == PARTICLE_CUSTOM) { - if (tree == NULL || treestyle != PARTICLE) + if (tree == NULL || treestyle != PARTICLE) { + delete [] id; error->all(FLERR,"Per-particle custom attribute in " "non particle-style variable formula"); + } icustom = particle->find_custom(id); - if (icustom < 0) + if (icustom < 0) { + delete [] id; error->all(FLERR,"Invalid custom attribute ID in variable formula"); + } size = particle->esize[icustom]; type = particle->etype[icustom]; } else if (cwhich == GRID_CUSTOM) { - if (tree == NULL || treestyle != GRID) + if (tree == NULL || treestyle != GRID) { + delete [] id; error->all(FLERR,"Per-grid custom attribute in " "non grid-style variable formula"); + } icustom = grid->find_custom(id); - if (icustom < 0) + if (icustom < 0) { + delete [] id; error->all(FLERR,"Invalid custom attribute ID in variable formula"); + } size = grid->esize[icustom]; type = grid->etype[icustom]; } else if (cwhich == SURF_CUSTOM) { - if (tree == NULL || treestyle != SURF) + if (tree == NULL || treestyle != SURF) { + delete [] id; error->all(FLERR,"Per-surf custom attribute in " "non surf-style variable formula"); + } icustom = surf->find_custom(id); - if (icustom < 0) + if (icustom < 0) { + delete [] id; error->all(FLERR,"Invalid custom attribute ID in variable formula"); + } size = surf->esize[icustom]; type = surf->etype[icustom]; } @@ -1947,8 +1967,10 @@ double Variable::evaluate(char *str, Tree **tree) strcpy(id,&word[3]); int isc = surf->find_collide(id); - if (isc < 0) + if (isc < 0) { + delete [] id; error->all(FLERR,"Invalid surf collide ID in variable formula"); + } SurfCollide *sc = surf->sc[isc]; delete [] id; @@ -2002,8 +2024,10 @@ double Variable::evaluate(char *str, Tree **tree) strcpy(id,&word[3]); int isr = surf->find_react(id); - if (isr < 0) + if (isr < 0) { + delete [] id; error->all(FLERR,"Invalid surf reaction ID in variable formula"); + } SurfReact *sr = surf->sr[isr]; delete [] id; @@ -2054,10 +2078,14 @@ double Variable::evaluate(char *str, Tree **tree) strcpy(id,&word[2]); int ivar = find(id); - if (ivar < 0) + if (ivar < 0) { + delete [] id; error->all(FLERR,"Invalid variable name in variable formula"); - if (eval_in_progress[ivar]) + } + if (eval_in_progress[ivar]) { + delete [] id; error->all(FLERR,"Variable has circular dependency"); + } // parse zero or one trailing brackets // point i beyond last bracket @@ -2092,8 +2120,10 @@ double Variable::evaluate(char *str, Tree **tree) style[ivar] != GRID && style[ivar] != SURF) { char *var = retrieve(id); - if (var == NULL) + if (var == NULL) { + delete [] id; error->all(FLERR,"Invalid variable evaluation in variable formula"); + } if (tree) { Tree *newtree = new Tree(); newtree->type = VALUE; @@ -2106,9 +2136,11 @@ double Variable::evaluate(char *str, Tree **tree) } else if (nbracket == 0 && style[ivar] == PARTICLE) { - if (tree == NULL || treestyle != PARTICLE) + if (tree == NULL || treestyle != PARTICLE) { + delete [] id; error->all(FLERR,"Per-particle variable in " "non particle-style variable formula"); + } Tree *newtree; evaluate(data[ivar][0],&newtree); treestack[ntreestack++] = newtree; @@ -2118,9 +2150,11 @@ double Variable::evaluate(char *str, Tree **tree) } else if (nbracket == 0 && style[ivar] == GRID) { - if (tree == NULL || treestyle != GRID) + if (tree == NULL || treestyle != GRID) { + delete [] id; error->all(FLERR,"Per-grid variable in " "non grid-style variable formula"); + } Tree *newtree; evaluate(data[ivar][0],&newtree); treestack[ntreestack++] = newtree; @@ -2130,16 +2164,21 @@ double Variable::evaluate(char *str, Tree **tree) } else if (nbracket == 0 && style[ivar] == SURF) { - if (tree == NULL || treestyle != SURF) + if (tree == NULL || treestyle != SURF) { + delete [] id; error->all(FLERR,"Per-surf variable in " "non surf-style variable formula"); + } Tree *newtree; evaluate(data[ivar][0],&newtree); treestack[ntreestack++] = newtree; // unrecognized variable - } else error->all(FLERR,"Mismatched variable in variable formula"); + } else { + delete [] id; + error->all(FLERR,"Mismatched variable in variable formula"); + } delete [] id; diff --git a/src/write_grid.cpp b/src/write_grid.cpp index c26aed20a..418b80c1a 100644 --- a/src/write_grid.cpp +++ b/src/write_grid.cpp @@ -80,7 +80,7 @@ void WriteGrid::command(int narg, char **arg) fp = fopen(arg[0],"w"); if (!fp) { char str[128]; - sprintf(str,"Cannot open file %s",arg[0]); + snprintf(str,128,"Cannot open file %s",arg[0]); error->one(FLERR,str); } } diff --git a/src/write_isurf.cpp b/src/write_isurf.cpp index 3dfb558ee..112602bf8 100644 --- a/src/write_isurf.cpp +++ b/src/write_isurf.cpp @@ -63,11 +63,14 @@ void WriteISurf::command(int narg, char **arg) char *ptr; int n = strlen(arg[4]) + 16; char *file = new char[n]; + char *arg4_copy = new char[strlen(arg[4]) + 1]; + strcpy(arg4_copy, arg[4]); - if ((ptr = strchr(arg[4],'*'))) { + if ((ptr = strchr(arg4_copy,'*'))) { *ptr = '\0'; - sprintf(file,"%s" BIGINT_FORMAT "%s",arg[4],update->ntimestep,ptr+1); - } else strcpy(file,arg[4]); + sprintf(file,"%s" BIGINT_FORMAT "%s",arg4_copy,update->ntimestep,ptr+1); + } else strcpy(file,arg4_copy); + delete [] arg4_copy; // ablation fix ID @@ -124,7 +127,7 @@ void WriteISurf::command(int narg, char **arg) fp = fopen(file,"wb"); if (!fp) { char str[128]; - sprintf(str,"Cannot open grid corner point file %s",file); + snprintf(str,128,"Cannot open grid corner point file %s",file); error->one(FLERR,str); } } @@ -160,6 +163,8 @@ void WriteISurf::command(int narg, char **arg) 100.0*(time2-time1)/time_total,100.0*(time3-time2)/time_total); } } + + delete [] file; } /* ---------------------------------------------------------------------- diff --git a/src/write_restart.cpp b/src/write_restart.cpp index e52a1fb80..8b1377829 100644 --- a/src/write_restart.cpp +++ b/src/write_restart.cpp @@ -202,7 +202,7 @@ void WriteRestart::write(char *file) fp = fopen(hfile,"wb"); if (fp == NULL) { char str[128]; - sprintf(str,"Cannot open restart file %s",hfile); + snprintf(str,128,"Cannot open restart file %s",hfile); error->one(FLERR,str); } if (multiproc) delete [] hfile; @@ -255,7 +255,7 @@ void WriteRestart::write(char *file) fp = fopen(multiname,"wb"); if (fp == NULL) { char str[128]; - sprintf(str,"Cannot open restart file %s",multiname); + snprintf(str,128,"Cannot open restart file %s",multiname); error->one(FLERR,str); } write_int(PROCSPERFILE,nclusterprocs); @@ -387,7 +387,7 @@ void WriteRestart::write_less_memory(char *file) fp = fopen(hfile,"wb"); if (fp == NULL) { char str[128]; - sprintf(str,"Cannot open restart file %s",hfile); + snprintf(str,128,"Cannot open restart file %s",hfile); error->one(FLERR,str); } if (multiproc) delete [] hfile; @@ -440,7 +440,7 @@ void WriteRestart::write_less_memory(char *file) fp = fopen(multiname,"wb"); if (fp == NULL) { char str[128]; - sprintf(str,"Cannot open restart file %s",multiname); + snprintf(str,128,"Cannot open restart file %s",multiname); error->one(FLERR,str); } write_int(PROCSPERFILE,nclusterprocs); diff --git a/src/write_surf.cpp b/src/write_surf.cpp index 08acc0e8b..65d29387a 100644 --- a/src/write_surf.cpp +++ b/src/write_surf.cpp @@ -949,7 +949,7 @@ void WriteSurf::write_base(char *file) fp = fopen(hfile,"w"); if (fp == NULL) { char str[128]; - sprintf(str,"Cannot open surface base file %s",hfile); + snprintf(str,128,"Cannot open surface base file %s",hfile); error->one(FLERR,str); } @@ -985,7 +985,7 @@ void WriteSurf::open(char *file) fp = fopen(onefile,"w"); if (fp == NULL) { char str[128]; - sprintf(str,"Cannot open surface file %s",onefile); + snprintf(str,128,"Cannot open surface file %s",onefile); error->one(FLERR,str); } From 0f645cc4e66e45aeafcd8f412b4d79e3f4acedf8 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 20 Aug 2026 15:49:14 +0000 Subject: [PATCH 2/4] Revert two fixes that example-deck regression testing rejected Reproduced the 8 failing CTest decks locally and bisected against a master-built binary; two applied entries from the bug list were the causes and are reverted here: - marching_cubes.h int-to-double corner values (reported bug 30): in.ablation.3d hard-errors at step ~50 with "Some cell faces do not have zero or 2 triangles" and in.exp2imp.sphere.3d results shift. The integer truncation of ablated corner values is load-bearing for the marching-cubes shared-face topology consistency. The file is byte-identical to master again. - react_tce.cpp warning relocation (reported bug 2): numerics are unchanged (in.chem / in.chem.gastally stats are bit-identical to master), but the resurrected per-collision "Reaction probability greater than 1" warnings flood logs and break reference-log comparison; the check also fires on the cumulative roulette probability, which legitimately exceeds 1. Master placement restored; only the computeChemRates missing-break fix remains in this file. After these reverts all four decks (ablation.3d, chem, chem.gastally, exp2imp.sphere.3d) complete against the same inputs with bit-identical stats to a master-built binary; serial build is clean. Co-authored-by: Stan Moore Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01JAGuLGWVKDbVYbwNXjqZ6K --- src/marching_cubes.h | 2 +- src/react_tce.cpp | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/marching_cubes.h b/src/marching_cubes.h index 388c41f3e..2957d4c3a 100644 --- a/src/marching_cubes.h +++ b/src/marching_cubes.h @@ -33,7 +33,7 @@ class MarchingCubes : protected Pointers { double thresh; double *lo,*hi; - double v000,v001,v010,v011,v100,v101,v110,v111; + int v000,v001,v010,v011,v100,v101,v110,v111; double v000iso,v001iso,v010iso,v011iso,v100iso,v101iso,v110iso,v111iso; double inval[8][6]; double i0,i1,i2,i3,i4,i5,i6,i7,i8,i9,i10,i11; diff --git a/src/react_tce.cpp b/src/react_tce.cpp index 51b42059e..4d1e89f5f 100644 --- a/src/react_tce.cpp +++ b/src/react_tce.cpp @@ -179,14 +179,14 @@ int ReactTCE::attempt(Particle::OnePart *ip, Particle::OnePart *jp, break; } + if (react_prob < 0) error->warning(FLERR,"Negative reaction probability"); + else if (react_prob > 1) error->warning(FLERR,"Reaction probability greater than 1"); + default: error->one(FLERR,"Unknown outcome in reaction"); break; } - if (react_prob < 0) error->warning(FLERR,"Negative reaction probability"); - else if (react_prob > 1) error->warning(FLERR,"Reaction probability greater than 1"); - // test against random number to see if this reaction occurs // if it does, reset species of I,J and optional K to product species // J particle is destroyed in recombination reaction, set species = -1 From df449c6c02a0662a891839d0ee48d8affca17e6c Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 20 Aug 2026 16:10:48 +0000 Subject: [PATCH 3/4] react_tce: replace dead probability warnings with a screen-only warn-once The react_prob validity warnings sat between the last case's break and the default: label of the reaction-type switch, i.e. inside the switch with no path to them - dead code that has never fired. The check itself is meaningful: react_prob accumulates the roulette-wheel sum of TCE reaction probabilities, and a sum above 1.0 means a reaction is guaranteed and the excess probability mass is clipped, biasing reaction rates (the classic under-resolved-chemistry symptom); a negative value means bad reaction-file coefficients. Restore the diagnostic without the two failure modes that a naive relocation has: - warn once per run (flag reset in init), not once per clipped collision per rank, which would flood output on production decks - write it to the screen only (logflag 0): a warning line inside the log's stats block breaks log parsers, including the regression harness in tools/testing, whose olog reader float()s every line between the Step header and "Loop time of" Verified on in.chem / in.chem.gastally: the log file is byte-identical in stats to a master-built binary (existing gold-standard logs still pass, confirmed by parsing the new log with tools/testing olog.py), and stdout carries exactly one warning per run. Co-authored-by: Stan Moore Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01JAGuLGWVKDbVYbwNXjqZ6K --- src/react_tce.cpp | 31 +++++++++++++++++++++++++++---- src/react_tce.h | 3 +++ 2 files changed, 30 insertions(+), 4 deletions(-) diff --git a/src/react_tce.cpp b/src/react_tce.cpp index 4d1e89f5f..6cba03862 100644 --- a/src/react_tce.cpp +++ b/src/react_tce.cpp @@ -30,7 +30,10 @@ enum{DISSOCIATION,EXCHANGE,IONIZATION,RECOMBINATION}; // other files /* ---------------------------------------------------------------------- */ ReactTCE::ReactTCE(SPARTA *sparta, int narg, char **arg) : - ReactBird(sparta, narg, arg) {} + ReactBird(sparta, narg, arg) +{ + prob_warn_flag = 0; +} /* ---------------------------------------------------------------------- */ @@ -39,6 +42,8 @@ void ReactTCE::init() if (!collide || strcmp(collide->style,"vss") != 0) error->all(FLERR,"React tce can only be used with collide vss"); + prob_warn_flag = 0; + ReactBird::init(); } @@ -179,14 +184,32 @@ int ReactTCE::attempt(Particle::OnePart *ip, Particle::OnePart *jp, break; } - if (react_prob < 0) error->warning(FLERR,"Negative reaction probability"); - else if (react_prob > 1) error->warning(FLERR,"Reaction probability greater than 1"); - default: error->one(FLERR,"Unknown outcome in reaction"); break; } + // sum of reaction probabilities should be < 1 for a valid TCE scheme, + // else reaction rates are biased by clipping + // warn only once per run to avoid flooding output, and only to the + // screen: a warning inside the log's stats block would break + // log-parsing tools (e.g. tools/testing regression harness) + + if (!prob_warn_flag) { + if (react_prob < 0.0) { + prob_warn_flag = 1; + error->warning(FLERR,"Negative TCE reaction probability, " + "check reaction file coefficients " + "(further warnings suppressed)",0); + } else if (react_prob > 1.0) { + prob_warn_flag = 1; + error->warning(FLERR,"TCE reaction probability exceeded 1.0, " + "chemistry may be under-resolved, " + "consider reducing timestep or fnum " + "(further warnings suppressed)",0); + } + } + // test against random number to see if this reaction occurs // if it does, reset species of I,J and optional K to product species // J particle is destroyed in recombination reaction, set species = -1 diff --git a/src/react_tce.h b/src/react_tce.h index 77f9da441..b36e30b5d 100644 --- a/src/react_tce.h +++ b/src/react_tce.h @@ -45,6 +45,9 @@ class ReactTCE : public ReactBird { double bird_Evib(int nmode, double Tvib, double VibTemp[], double Evib); + + protected: + int prob_warn_flag; // 1 after warning once about an invalid react_prob }; } From 72a843a5fcde180a1ba079d8e88371d0495de6e3 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 20 Aug 2026 16:22:48 +0000 Subject: [PATCH 4/4] Harden olog against messages in the stats block; log the TCE warning again tools/testing/olog.py float()ed every whitespace token of every line between the Step header and "Loop time of", so any message emitted mid-run (a WARNING, a fix grid/check report) crashed the parse and failed the regression test regardless of the gold standard. Skip blank and non-numeric lines instead, and guard the last-timestep print for the empty case. With the parser tolerant, drop the screen-only restriction on the react_tce warn-once diagnostic so it also lands in the log file, where users doing log-only runs will see it. Verified with the harness's own code paths: extract_data over the warning-bearing in.chem log passes and yields 11 rows with data identical to a master-built binary (CPU column excluded, as the harness does); the shipped gold standard log.11Sep23.mpi_1.chem still parses identically under the hardened reader, so no reference logs need regeneration. Co-authored-by: Stan Moore Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01JAGuLGWVKDbVYbwNXjqZ6K --- src/react_tce.cpp | 8 +++----- tools/pizza/olog.py | 13 ++++++++++--- 2 files changed, 13 insertions(+), 8 deletions(-) diff --git a/src/react_tce.cpp b/src/react_tce.cpp index 6cba03862..146013d61 100644 --- a/src/react_tce.cpp +++ b/src/react_tce.cpp @@ -191,22 +191,20 @@ int ReactTCE::attempt(Particle::OnePart *ip, Particle::OnePart *jp, // sum of reaction probabilities should be < 1 for a valid TCE scheme, // else reaction rates are biased by clipping - // warn only once per run to avoid flooding output, and only to the - // screen: a warning inside the log's stats block would break - // log-parsing tools (e.g. tools/testing regression harness) + // warn only once per run to avoid flooding output if (!prob_warn_flag) { if (react_prob < 0.0) { prob_warn_flag = 1; error->warning(FLERR,"Negative TCE reaction probability, " "check reaction file coefficients " - "(further warnings suppressed)",0); + "(further warnings suppressed)"); } else if (react_prob > 1.0) { prob_warn_flag = 1; error->warning(FLERR,"TCE reaction probability exceeded 1.0, " "chemistry may be under-resolved, " "consider reducing timestep or fnum " - "(further warnings suppressed)",0); + "(further warnings suppressed)"); } } diff --git a/tools/pizza/olog.py b/tools/pizza/olog.py index 47ee3d9d3..4923b0c00 100644 --- a/tools/pizza/olog.py +++ b/tools/pizza/olog.py @@ -289,15 +289,22 @@ def read_one(self,*list): # split chunk into entries # parse each entry for numeric fields, append to data + # skip blank lines and any non-numeric line (e.g. a WARNING message + # emitted mid-run) so embedded messages do not break the parse lines = chunk.split("\n") for line in lines: words = line.split() - self.data.append([float(i) for i in words]) + if not words: continue + try: + self.data.append([float(i) for i in words]) + except ValueError: + continue # print last timestep of chunk - print(int(self.data[len(self.data)-1][0]), end=' ') - sys.stdout.flush() + if self.data: + print(int(self.data[len(self.data)-1][0]), end=' ') + sys.stdout.flush() return eof