diff --git a/BUG_VERIFICATION.md b/BUG_VERIFICATION.md new file mode 100644 index 000000000..cff2bb021 --- /dev/null +++ b/BUG_VERIFICATION.md @@ -0,0 +1,243 @@ +# SPARTA AI-Found Bug Verification & Fix Audit + +This report independently verifies the list of ~111 bugs (Bugs **1–103, 105–112**; there is +no #104) reported by two AI tools, and audits the two fix branches against the canonical +baseline `origin/master` (`5aed836`): + +- **AB** = `origin/ai_bugfixes` (aborner, commit `d3f62a1`) — broad (128 src files, 33 KOKKOS). +- **CP** = `origin/copilot/fix-bugs-from-bugs-md` (`39c7a3a`) — narrower (78 src, 19 KOKKOS). + +For every bug we determined: (1) whether it is a *genuine* defect in `origin/master`, and +(2) whether each branch's fix is correct. The verified-correct fixes were then applied to this +branch (`claude/sparta-static-analysis-bugs-335wsw`, started fresh from `origin/master`); +false positives and incorrect/regressive fixes were rejected. + +## Headline results + +- **111** numbered entries audited. Of these, **3 are duplicates** of another entry + (23 = 67 = 98, same `fix_ablate` line) and **4 are NOT bugs** (40, 84, 93, 102). +- **~104 genuine defects** confirmed (REAL or PARTIAL-but-genuine). +- **Neither branch is complete or correct on its own:** + - **AB** has **2 incorrect fixes** (Bug 46 and Bug 40) and misses several pure-logic bugs + (3, 5, 9, 30, 47, 53, 55, 56, 57, 58, 81). + - **CP** has **2 incorrect/partial fixes** (Bug 40, and the `react_qk` half of Bug 72) and + misses most overflow/NaN/KOKKOS hardening (62, 64, 65, 73(½), 76, 77, 78, 79, 80, 87, 88, + 89, 90, 91, 92, 94, 96, 99, 100, 101, 102, 103, 110 …). + - Both branches **fix only the KOKKOS copies** of Bugs 62/77/78/79/80, leaving the identical + **CPU** defects in `update.cpp` / `geometry.cpp` unfixed. +- The deliverable branch applies the correct fix for every genuine bug (taking AB or CP + whichever is right, or a corrected/extended fix where both were wrong or incomplete), + **rejects** the 4 non-bugs, and **avoids AB's EPSZERO regression** in `collide_vss*`. +- **Build:** `make serial` links cleanly (`spa_serial`) with all non-KOKKOS fixes applied. + +## Legend + +- **Real?**: REAL · PARTIAL (genuine but narrow/defensive) · **NO** (not a bug) +- **AB / CP**: OK (correct) · DIFF (correct, alternative) · NO (not fixed) · **WRONG** · PARTIAL +- **Applied**: ✔ applied · ✔* applied (corrected/extended beyond both branches) · ✖ rejected + +| # | File(s) | Real? | AB | CP | Applied | Note | +|--|--|--|--|--|--|--| +| 1 | comm.cpp | REAL | OK | OK | ✔ | spurious double alloc of `rbuf` (bigint→int truncation) | +| 2 | react_tce.cpp | REAL | OK | OK | ✔ | unreachable react_prob warning moved out of `switch` | +| 3 | compute_reduce.cpp | REAL | NO | OK | ✔ | `narg`→`nargnew` bounds (replace & subset) — AB missed | +| 4 | react_qk.cpp | REAL | OK | WRONG | ✔ | see Bug 72; AB correct, CP pollutes `react_prob` | +| 5 | grid_custom.cpp | REAL | NO | OK | ✔ | `if(nnew-nold)`→`if(nnew>nold)` (huge memset on shrink) | +| 6 | compute_property_surf.cpp | REAL | OK | OK | ✔ | 3D pack_id `nsown`→`nchoose` OOB | +| 7 | utils.cpp | PARTIAL | OK | OK | ✔ | snprintf hardening (unbounded `cmd` string) | +| 8 | variable.cpp | REAL | OK | OK | ✔ | unchecked fopen of lock file | +| 9 | compute_count.cpp | REAL | NO | OK | ✔ | `if(imix<0)`→`if(igroup<0)` — AB missed | +| 10 | compute_gas_collision_tally.cpp | REAL | OK | OK | ✔ | `type2`→TYPE1 should be TYPE2 | +| 11 | compute_gas_reaction_tally.cpp | REAL | OK | OK | ✔ | vy2/vz2 pre-vel enum shift | +| 12 | compute_surf_reaction_tally.cpp | REAL | OK | OK | ✔ | `ID2POST||ID2POST`→`ID1POST||ID2POST` | +| 13 | fix_grid_check.cpp | REAL | OK | OK | ✔ | OOB `cells[icell]` in error msg + missing `continue` | +| 14 | timer.cpp | REAL | OK | OK | ✔ | uninit `timeout_start`; stray debug printf | +| 15 | surf_collide_specular.cpp | REAL | OK | OK | ✔ | wrapper() ignored noslip_flag | +| 16 | surf_collide_cll.cpp | REAL | OK | OK | ✔ | missing `if(copy) return` → double free | +| 17 | surf_collide_impulsive.cpp | REAL | OK | OK | ✔ | missing `if(copy) return` → double free | +| 18 | KOKKOS/fix_grid_check_kokkos.cpp | REAL | OK | OK | ✔ | OOB device access; add `return` | +| 19 | KOKKOS/compute_lambda_grid_kokkos.cpp | REAL | OK | OK | ✔ | missing `else` KNY/KNZ (subsumed by 38) | +| 20 | KOKKOS/fix_ave_histo_weight_kokkos.cpp | REAL | OK | OK | ✔ | realloc `>nmax`→`boltz | +| 49 | react_tce.cpp | REAL | OK | OK | ✔ | missing `break` inflates chem-rate tallies | +| 50 | fix_surf_temp.cpp | REAL | OK | OK | ✔ | uninit prefactor/threshold (no else) | +| 51 | fix_surf_temp.cpp | REAL | OK | OK | ✔ | stale cqw/fqw; re-resolve in init() | +| 52 | compute_lambda_grid.cpp | REAL | OK | OK | ✔ | CPU twin of 19/38 | +| 53 | update.cpp | REAL | NO | OK | ✔ | dangling *_active → double free — AB missed | +| 54 | react_bird.cpp | REAL | OK | OK | ✔ | uninit tally_reactions* in 1-arg ctor | +| 55 | particle.cpp | REAL | NO | OK | ✔ | size_restart int overflow — AB missed | +| 56 | grid_collate.cpp | REAL | NO | OK | ✔ | 32-bit memset/create — AB missed | +| 57 | grid_custom.cpp | REAL | NO | PARTIAL | ✔* | (size_t) on int **and** double memset (both missed double) | +| 58 | grid.cpp | REAL | NO | PARTIAL | ✔* | bigint on set1 **and** set2 (CP did set1 only) | +| 59 | surf_custom.cpp | REAL | OK | OK | ✔ | (size_t) memset | +| 60 | surf_comm.cpp | REAL | OK | PARTIAL | ✔ | bigint ×4 sites; CP missed dbuf → took AB | +| 61 | surf.cpp | PARTIAL | OK | PARTIAL | ✔ | effective fix = (size_t) on Tri/Line memset | +| 62 | update.h, geometry.cpp, KOKKOS update/geometry | REAL | PARTIAL | NO | ✔* | axi div-by-zero; **both left CPU geometry.cpp — fixed here** | +| 63 | fix_*/compute_* (sweep) | PARTIAL | OK | PARTIAL | ✔ partial | genuine %s sites applied via 28/32/94; pure-%d rejected | +| 64 | compute_gas_reaction_/collision_grid.cpp | PARTIAL | OK | NO | ✔ | (size_t) memset | +| 65 | fix_emit_face/face_file/surf.cpp | PARTIAL | OK | NO | ✔ | (size_t) memset maxactive/ntaskmax | +| 66 | fix_ablate_multi_inner.cpp | REAL | OK | OK | ✔ | /Ninterface==0 NaN (report's "SIGFPE" inexact) | +| 67 | fix_ablate.cpp | REAL | OK | OK | ✔ | duplicate of 23 | +| 68 | compute_tvib_grid.cpp | REAL | OK | OK | ✔ | groupspecies OOB → `index/maxmode` | +| 69 | compute_react_surf/boundary/isurf_grid.cpp | REAL | OK | OK | ✔ | strtok loop wipes prior matches | +| 70 | fix_temp_rescale.cpp | REAL | PARTIAL | OK | ✔ | t_current==0; CP guards both paths, AB only avg | +| 71 | compute_{eflux,grid,pflux,sonine,thermal,tvib,lambda}_grid.cpp | REAL | OK | OK | ✔ | memory_usage `=`→`+=` (under-report) | +| 72 | react_qk.cpp, react_tce_qk.cpp | REAL | OK | PARTIAL | ✔ | scratch `prob` vs `react_prob`; **CP wrong on react_qk** | +| 73 | KOKKOS/collide_vss_kokkos.cpp | REAL | OK | PARTIAL | ✔ | free_state before continue (PRNG race); CP missed ambipolar loop | +| 74 | KOKKOS/react_tce_kokkos.h | PARTIAL | OK | OK | ✔ | kb→boltz (premise "base uses boltz" was false; paired with 48) | +| 75 | KOKKOS/react_bird_kokkos.cpp | PARTIAL | OK | NO | ✔ | pool seed decorrelation (speculative; pool-seed only applied) | +| 76 | KOKKOS/collide_vss_kokkos.cpp, react_tce*.h | PARTIAL | OK | NO | ✔ | vr2>0 guard (the reachable one); ecc guard via 106 | +| 77 | update.cpp, KOKKOS/update_kokkos.cpp | REAL | PARTIAL | NO | ✔* | frac 0/0; **AB KOKKOS only — CPU fixed here** | +| 78 | update.cpp, KOKKOS/update_kokkos.cpp | PARTIAL | PARTIAL | NO | ✔* | clamp frac∈[0,1]; **CPU fixed here** | +| 79 | update.cpp, KOKKOS/update_kokkos.cpp | REAL | PARTIAL | NO | ✔* | stuck_iterate `==0`→`<=1e-14`; **CPU fixed here** | +| 80 | geometry.cpp, KOKKOS/geometry_kokkos.h | REAL | PARTIAL | NO | ✔* | catastrophic cancellation (Vieta); **CPU fixed here** | +| 81 | adapt_grid.cpp | REAL | NO | OK | ✔ | newcell→nglocalprev — AB missed | +| 82 | adapt_grid.cpp | PARTIAL | NO | NO | ✔* | (bigint) smalloc — neither fixed | +| 83 | adapt_grid.cpp | PARTIAL | NO | NO | ✔* | (bigint) clist/alist srealloc — neither fixed | +| **84** | adapt_grid.cpp/.h, grid_adapt.cpp | **NO** | NO | NO | **✖** | speculative int-widening; total size already bigint | +| 85 | KOKKOS/fix_ambipolar_kokkos.h | REAL | OK | OK | ✔ | -log(0)=+Inf → 1.0-drand() | +| 86 | KOKKOS/fix_vibmode_kokkos.h | REAL | OK | OK | ✔ | -log(0)→Inf cast to int (UB) | +| 87 | KOKKOS/fix_emit_*/diffuse/particle_kokkos | REAL | OK | NO | ✔ | 13 -log(drand) sites; CP missed all | +| 88 | create_isurf.cpp | REAL | OK | NO | ✔ | (bigint)/(size_t) nsurf*nbytes | +| 89 | create_*/dump_movie.cpp | PARTIAL | OK | NO | ✔ partial | dump_movie filename applied; create_* pure-%d rejected | +| 90 | create_isurf.cpp/.h | REAL | OK | NO | ✔ | maxsbuf int→bigint | +| 91 | create_isurf.cpp | REAL | OK | NO | ✔ | boxvol==0 div | +| 92 | create_isurf.cpp | REAL | OK | NO | ✔ | param==1 → Inf clamp | +| 93 | fix_move_surf/emit*.cpp | **NO** | DIFF | NO | **✖** | `int*sizeof` already promotes to 64-bit; casts inert | +| 94 | KOKKOS/fix_grid_check_kokkos.cpp | REAL | OK | NO | ✔ | sprintf→snprintf | +| 95 | fix_emit_face_file.cpp | REAL | OK | OK | ✔ | `for(m…;i++)` typo → infinite/OOB (HIGH) | +| 96 | fix_emit_face/face_file/surf.cpp | REAL | OK | NO | ✔ | subsonic /0 (massrho*soundspeed) | +| 97 | fix_temp_rescale.cpp | REAL | OK | OK | ✔ | global avg /0 (paired with 70) | +| 98 | fix_ablate.cpp | REAL | OK | OK | ✔ | duplicate of 23 | +| 99 | fix_ablate/ave_grid/histo/surf/time.cpp | REAL | OK | NO | ✔ | suffix leak on error path | +| 100 | compute_dt_grid.cpp + KOKKOS | REAL | OK | PARTIAL | ✔ | vrm_max==0 div; CP missed KOKKOS | +| 101 | compute_{thermal,eflux,pflux}_grid.cpp + KOKKOS | REAL | OK | NO | ✔ | volume==0 div in flux tallies | +| **102** | KOKKOS/compute_fft_grid_kokkos.cpp | **NO** | OK | NO | **✖** | sprintf into str[64] cannot overflow (cosmetic) | +| 103 | compute_lambda_grid.cpp, compute_reduce.cpp | REAL | OK | NO | ✔ | suffix leak on error | +| 105 | collide_vss.cpp | REAL | OK | OK | ✔ | volume==0 div in nattempt | +| 106 | react_qk/tce/tce_qk.cpp | PARTIAL | OK | OK | ✔ | ecc>0 guards (mostly already guarded; via 2/48/72) | +| 107 | surf_react_adsorb.cpp | PARTIAL | OK | OK | ✔ | CI vmag_sq>0 (ER dot==0 moot: dot hardcoded 2.0) | +| 108 | variable.cpp | REAL | OK | OK | ✔ | id leak on error | +| 109 | grid.cpp, input.cpp (sweep) | PARTIAL | DIFF | PARTIAL | ✔ partial | %s/filename sites applied (28/32); pure-%d rejected | +| 110 | dump.cpp, dump_grid.cpp, grid_id.cpp | REAL | OK | NO | ✔ | str[32] too small for deep cell ids | +| 111 | input.cpp | REAL | OK | OK | ✔ | commands[] leak on illegal `if` | +| 112 | grid.cpp | REAL | OK | OK | ✔ | list leak in grid group ops | + +## Notable findings + +**False positives (rejected — applying them would be wrong or pointless):** +- **Bug 40** (`fft2d_kokkos.cpp`): the report misread this function's *inverted* flag + convention — here `flag==-1` is the forward transform, so the existing `flag==1` scaling + already targets the inverse. **Both branches changed it and thereby moved scaling onto the + forward transform — an actual regression.** Rejected. +- **Bug 84** (`adapt_grid` int widening): the only genuinely size-scaling quantity + (`plevels[].nxyz` / total children) is already `bigint`; the `int` locals hold per-parent + subdivision factors that cannot overflow. Both branches correctly left it alone. +- **Bug 93** (`fix_move_surf`/emit `nsurf*sizeof(...)`): `sizeof` is `size_t`, so the + multiplication is already 64-bit. AB's added casts are inert; rejected to keep the diff honest. +- **Bug 102** (`compute_fft_grid_kokkos` str[64]): bounded macro text, no overflow possible. + +**Incorrect fixes in a branch (we took the other branch / corrected it):** +- **Bug 46** — AB not only failed to add the `vremax==0` guard but **deleted the pre-existing + `EPSZERO` division-by-zero guard** (and its macro) in both `collide_vss.cpp` and + `collide_vss_kokkos.cpp` — a regression. We took CP for the CPU file and applied only the + correct sub-edits (73, 76) to the KOKKOS file, **keeping EPSZERO**. +- **Bug 72 (react_qk.cpp)** — CP rewrote the first rejection-sampling loop to keep using + `react_prob` as the loop scratch, which leaves the sampling probability in `react_prob` and + pollutes the downstream reaction decision (`react_prob > random_prob`). AB correctly uses a + separate `prob` scratch (matching the file's own second loop). We took AB. + +**Completeness gaps both branches share (fixed here, marked ✔\*):** +- Bugs **62, 77, 78, 79, 80** were fixed by AB only in the KOKKOS copies; the identical CPU + defects in `update.cpp` and `geometry.cpp` were left unfixed by both. We ported the + verified-correct fixes to the CPU files (geometry.cpp Vieta rewrite was hand-verified: + product-of-roots `c/a` with `c = x[1]²−yhoriz²`). +- Bugs **57, 58, 60, 73, 100** were only partially fixed by the branch that touched them + (missing a second memset / set2 / dbuf / ambipolar loop / KOKKOS file); completed here. +- Bugs **82, 83** were fixed by neither branch; applied here. + +**Duplicates:** Bugs **23 = 67 = 98** are the same `fix_ablate.cpp` `idsource` shadowing line +(one fix). Bug **19** is a subset of **38**. Bugs **70/97** and **4/72/106** are related but +touch distinct sites. + +## Scope notes on hardening sweeps + +`sprintf`→`snprintf` and integer-cast "sweeps" (Bugs 63, 89, 109) were applied **only** where a +real overflow is possible — i.e. the format contains a `%s` of a filename / unbounded string. +Pure-numeric (`%d`, `BIGINT_FORMAT`) conversions into fixed buffers that cannot overflow were +treated as cosmetic non-bugs and not applied, to keep the deliverable scoped to genuine defects. + +## Verification performed + +- Each bug verified against `origin/master` (presence + genuineness) and against both branches + (fix correctness) by reading the actual source via `git show`. +- `make serial` builds and links `spa_serial` with all non-KOKKOS fixes applied. +- KOKKOS files are not part of a serial build (Kokkos not installed in `src/`); those edits are + either byte-identical to AB's (which compiles) or trivial targeted edits, and were diff-checked. + +## Appendix: `ai_bugfixes`-specific audit + +This appendix isolates what the **`ai_bugfixes` (AB)** branch got wrong or bundled, and where the +deliverable therefore departed from AB. + +### A. False positives in `ai_bugfixes` (AB changes that are not real bugs) + +| # | File | AB change | Why it is not a bug | +|--|--|--|--| +| **40** | `KOKKOS/fft2d_kokkos.cpp` | `if (flag == 1 …)` → `if (flag == -1 …)` | This function's flag convention is inverted (`-1` is the forward transform here), so the original `flag==1` scaling already targets the inverse. AB's change moves scaling onto the forward transform — a **regression**. | +| **93** | `fix_move_surf.cpp`, `fix_emit_face.cpp`, `fix_emit_face_file.cpp`, `fix_emit_surf.cpp` | added `(bigint)`/`(size_t)` casts to `nsurf*sizeof(...)`, `ntaskmax*sizeof(...)` | `sizeof` is `size_t`, so the products are already 64-bit. Casts are inert. | +| **102** | `KOKKOS/compute_fft_grid_kokkos.cpp` | `sprintf`→`snprintf` into `char str[64]` | Bounded macro text; cannot overflow. Cosmetic. | +| — | e.g. `collide_vss.cpp` "%d species…discrete model", `update_kokkos.cpp` BIGINT self-proc message, and similar | `sprintf`→`snprintf` on **bounded numeric** buffers | Defensive-only; no `%s`/unbounded input. Cosmetic non-bugs, not applied. | + +(AB correctly left **Bug 84** unchanged, so #84 is a report-level false positive, not an AB one.) + +### B. Bugs where the deliverable was *modified* from AB's version + +AB changed this code, but its version was wrong / regressive / carried unlisted extras, so AB was +**not** taken as-is: + +| # | File | AB's version | Deliverable | +|--|--|--|--| +| **46** | `collide_vss.cpp` | Did not add the `vremax==0` guard, and **deleted the pre-existing `EPSZERO` divide-by-zero guard** (+ its `#define`); also added `remain` outside the `volume>0` guard in the 2nd overload. | Took **copilot** (adds `vremax==0` guard, keeps `EPSZERO`). | +| **73 + 76** | `KOKKOS/collide_vss_kokkos.cpp` | Correct 73/76 edits **bundled with** the same `EPSZERO`-guard deletion (regression) + RNG-seed rewrites (unlisted). | Applied **only** the 73 (`free_state` before `continue`) and 76 (`vr2>0`) edits via targeted patches; **kept `EPSZERO`**; dropped the seed changes. | +| **75** | `KOKKOS/react_bird_kokkos.cpp` | Changed the `rand_pool` seed (listed) **and** rewrote the `random_backup` seed (unlisted). | Applied **only** the listed `rand_pool` seed change (`12345`→`54321`). | +| **54** | `react_bird.cpp` | `tally_reactions*` NULL init (listed) **plus** an `ncount>0` divide guard and a `sprintf`→`snprintf` (both unlisted). | Applied **only** Bug 54 (the two NULL inits). | +| **40 / 93 / 102** | see Part A | Wrong / inert changes. | **Rejected entirely.** | + +### C. Related context (AB fell short, but not an AB modification per se) + +- **AB fixed only the KOKKOS copies** of Bugs **62, 77, 78, 79, 80**; the identical CPU defects in + `update.cpp` / `geometry.cpp` were left unfixed by AB (and by copilot). The deliverable ports the + verified fixes to the CPU files. +- **AB never fixed** (deliverable took copilot or an original fix): Bugs **3, 5, 9, 30, 47, 53, 55, + 56, 81** (pure-logic bugs), and **82, 83** (fixed by neither branch). +- AB's unrelated `adapt_grid.cpp` warning rewrite dropped the message's `printf` arguments (a minor + regression). Avoided — the deliverable takes copilot's Bug 81 fix and does not import AB's + `adapt_grid.cpp` changes. diff --git a/src/KOKKOS/collide_vss_kokkos.cpp b/src/KOKKOS/collide_vss_kokkos.cpp index 6d4f5384e..12277101f 100644 --- a/src/KOKKOS/collide_vss_kokkos.cpp +++ b/src/KOKKOS/collide_vss_kokkos.cpp @@ -760,7 +760,6 @@ void CollideVSSKokkos::operator()(TagCollideCollisionsOne< NEARCP, GASTALLY, ATO else reduce.nreact_one++; } else { - rand_pool.free_state(rand_gen); continue; } @@ -1170,7 +1169,6 @@ void CollideVSSKokkos::operator()(TagCollideCollisionsOneAmbipolar< GASTALLY, AT else reduce.nreact_one++; } else { - rand_pool.free_state(rand_gen); continue; } @@ -1592,7 +1590,7 @@ void CollideVSSKokkos::SCATTER_TwoBodyScattering(Particle::OnePart *ip, vb = vr*sinX*cos(eps); wc = vr*sinX*sin(eps); } else { - double scale = sqrt((2.0 * postcoln.etrans) / (d_params(isp,jsp).mr * precoln.vr2)); + double scale = (precoln.vr2 > 0.0) ? sqrt((2.0 * postcoln.etrans) / (d_params(isp,jsp).mr * precoln.vr2)) : 0.0; double cosX = 2.0*pow(rand_gen.drand(),alpha_r) - 1.0; double sinX = sqrt(1.0 - cosX*cosX); vrc[0] = vi[0]-vj[0]; @@ -1799,7 +1797,7 @@ void CollideVSSKokkos::SCATTER_ThreeBodyScattering(Particle::OnePart *ip, vb = vr*sinX*cos(eps); wc = vr*sinX*sin(eps); } else { - double scale = sqrt((2.0*postcoln.etrans) / (mr*precoln.vr2)); + double scale = (precoln.vr2 > 0.0) ? sqrt((2.0*postcoln.etrans) / (mr*precoln.vr2)) : 0.0; vrc[0] = vi[0]-vj[0]; vrc[1] = vi[1]-vj[1]; vrc[2] = vi[2]-vj[2]; diff --git a/src/KOKKOS/compute_dt_grid_kokkos.cpp b/src/KOKKOS/compute_dt_grid_kokkos.cpp index 948d6c375..0c6b5c779 100644 --- a/src/KOKKOS/compute_dt_grid_kokkos.cpp +++ b/src/KOKKOS/compute_dt_grid_kokkos.cpp @@ -357,13 +357,15 @@ void ComputeDtGridKokkos::operator()(TagComputeDtGrid_ComputePerGrid, const int // cell dt based on transit time using maximum most probable speed double vrm_max = sqrt(2.0*boltz * d_temp_vector(i) / min_species_mass); - dt_candidate = transit_fraction*dx/vrm_max; - cell_dt_desired = MIN(dt_candidate,cell_dt_desired); - dt_candidate = transit_fraction*dy/vrm_max; - cell_dt_desired = MIN(dt_candidate,cell_dt_desired); - if (dimension == 3) { - dt_candidate = transit_fraction*dz/vrm_max; + if (vrm_max > 0.0) { + dt_candidate = transit_fraction*dx/vrm_max; + cell_dt_desired = MIN(dt_candidate,cell_dt_desired); + dt_candidate = transit_fraction*dy/vrm_max; cell_dt_desired = MIN(dt_candidate,cell_dt_desired); + if (dimension == 3) { + dt_candidate = transit_fraction*dz/vrm_max; + cell_dt_desired = MIN(dt_candidate,cell_dt_desired); + } } d_vector_grid(i) = cell_dt_desired; diff --git a/src/KOKKOS/compute_eflux_grid_kokkos.cpp b/src/KOKKOS/compute_eflux_grid_kokkos.cpp index 013220fcb..7f35f7e94 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; @@ -420,8 +420,12 @@ void ComputeEFluxGridKokkos::operator()(TagComputeEFluxGrid_post_process_grid, c h2 = d_etally(icell,mvv2v2) - 2.0*d_etally(icell,mvv2)*d_etally(icell,mv2)/summass - d_etally(icell,mv)*d_etally(icell,mv2v2)/summass + 2.0*d_etally(icell,mv)*d_etally(icell,mv2)*d_etally(icell,mv2)/summass/summass; - wt = 0.5 * fnum * d_cinfo[icell].weight / d_cinfo[icell].volume; - d_vec[icell] = wt/nsample * (h + h1 + h2); + if (d_cinfo[icell].volume > 0.0) { + wt = 0.5 * fnum * d_cinfo[icell].weight / d_cinfo[icell].volume; + d_vec[icell] = wt/nsample * (h + h1 + h2); + } else { + d_vec[icell] = 0.0; + } } } diff --git a/src/KOKKOS/compute_lambda_grid_kokkos.cpp b/src/KOKKOS/compute_lambda_grid_kokkos.cpp index f988a63fb..f2a062cd2 100644 --- a/src/KOKKOS/compute_lambda_grid_kokkos.cpp +++ b/src/KOKKOS/compute_lambda_grid_kokkos.cpp @@ -315,23 +315,43 @@ void ComputeLambdaGridKokkos::compute_per_grid_kokkos() sizeall += sizez; sizeall /= 3.0; } - if (l_noutputs == 1) l_vector_grid[i] = lambda / sizeall; - else l_array_grid(i,l_output_order[KNALL]) = lambda / sizeall; + if (sizeall > 0.0) { + if (l_noutputs == 1) l_vector_grid[i] = lambda / sizeall; + else l_array_grid(i,l_output_order[KNALL]) = lambda / sizeall; + } else { + if (l_noutputs == 1) l_vector_grid[i] = 0.0; + else l_array_grid(i,l_output_order[KNALL]) = 0.0; + } } if (l_knxflag) { - if (l_noutputs == 1) l_vector_grid[i] = lambda / sizex; - else l_array_grid(i,l_output_order[KNX]) = lambda / sizex; + if (sizex > 0.0) { + if (l_noutputs == 1) l_vector_grid[i] = lambda / sizex; + else l_array_grid(i,l_output_order[KNX]) = lambda / sizex; + } else { + if (l_noutputs == 1) l_vector_grid[i] = 0.0; + else l_array_grid(i,l_output_order[KNX]) = 0.0; + } } if (l_knyflag) { - if (l_noutputs == 1) l_vector_grid[i] = lambda / sizey; - l_array_grid(i,l_output_order[KNY]) = lambda / sizey; + if (sizey > 0.0) { + if (l_noutputs == 1) l_vector_grid[i] = lambda / sizey; + else l_array_grid(i,l_output_order[KNY]) = lambda / sizey; + } else { + if (l_noutputs == 1) l_vector_grid[i] = 0.0; + else l_array_grid(i,l_output_order[KNY]) = 0.0; + } } if (l_knzflag) { - if (l_noutputs == 1) l_vector_grid[i] = lambda / sizez; - l_array_grid(i,l_output_order[KNZ]) = lambda / sizez; + if (sizez > 0.0) { + if (l_noutputs == 1) l_vector_grid[i] = lambda / sizez; + else l_array_grid(i,l_output_order[KNZ]) = lambda / sizez; + } else { + if (l_noutputs == 1) l_vector_grid[i] = 0.0; + else l_array_grid(i,l_output_order[KNZ]) = 0.0; + } } }); } diff --git a/src/KOKKOS/compute_pflux_grid_kokkos.cpp b/src/KOKKOS/compute_pflux_grid_kokkos.cpp index a5ab56f6e..039cfd6bf 100644 --- a/src/KOKKOS/compute_pflux_grid_kokkos.cpp +++ b/src/KOKKOS/compute_pflux_grid_kokkos.cpp @@ -354,7 +354,7 @@ KOKKOS_INLINE_FUNCTION void ComputePFluxGridKokkos::operator()(TagComputePFluxGrid_post_process_grid_diag, const int &icell) const { double summass, summv, wt; summass = d_etally(icell,mass); - if (summass == 0.0) d_vec[icell] = 0.0; + if (summass == 0.0 || d_cinfo[icell].volume == 0.0) d_vec[icell] = 0.0; else{ wt = fnum * d_cinfo[icell].weight / d_cinfo[icell].volume; summv = d_etally(icell,mv); @@ -368,7 +368,7 @@ KOKKOS_INLINE_FUNCTION void ComputePFluxGridKokkos::operator()(TagComputePFluxGrid_post_process_grid_offdiag, const int &icell) const { double summass, wt; summass = d_etally(icell,mass); - if (summass == 0.0) d_vec[icell] = 0.0; + if (summass == 0.0 || d_cinfo[icell].volume == 0.0) d_vec[icell] = 0.0; else{ wt = fnum * d_cinfo[icell].weight / d_cinfo[icell].volume; d_vec[icell] = wt/nsample * (d_etally(icell,mvv) - 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 (istart < nsurf-2 && h_surf2tally[istart] != -1) istart++; + while (iend > 0 && h_surf2tally[iend] == -1) iend--; if (istart >= iend) { ntally = istart; break; diff --git a/src/KOKKOS/compute_thermal_grid_kokkos.cpp b/src/KOKKOS/compute_thermal_grid_kokkos.cpp index d022aba69..60bdfd076 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; @@ -281,7 +281,10 @@ void ComputeThermalGridKokkos::operator()(TagComputeThermalGrid_post_process_gri d_vec[icell] = mvsq - (mvx*mvx + mvy*mvy + mvz*mvz)/mass; d_vec[icell] *= prefactor; if (tflag) d_vec[icell] /= ncount; - else d_vec[icell] *= d_cinfo[icell].weight / d_cinfo[icell].volume / nsample; + else { + if (d_cinfo[icell].volume > 0.0) d_vec[icell] *= d_cinfo[icell].weight / d_cinfo[icell].volume / nsample; + else d_vec[icell] = 0.0; + } } } diff --git a/src/KOKKOS/create_particles_kokkos.cpp b/src/KOKKOS/create_particles_kokkos.cpp index bd0115ea2..c912003b4 100644 --- a/src/KOKKOS/create_particles_kokkos.cpp +++ b/src/KOKKOS/create_particles_kokkos.cpp @@ -423,7 +423,7 @@ void CreateParticlesKokkos::create_local(bigint np) // if using per-grid variables or per-grid custom attributes if (nfix_update_custom) - modify->update_custom(particle->nlocal-1,tempscale*temp_thermal, + modify->update_custom(inew,tempscale*temp_thermal, tempscale*temp_rot,tempscale*temp_vib, vstream_update_custom); } diff --git a/src/KOKKOS/fix_ambipolar_kokkos.h b/src/KOKKOS/fix_ambipolar_kokkos.h index f5eb10526..0cca1a8e8 100644 --- a/src/KOKKOS/fix_ambipolar_kokkos.h +++ b/src/KOKKOS/fix_ambipolar_kokkos.h @@ -93,8 +93,8 @@ class FixAmbipolarKokkos : public FixAmbipolar { rand_type rand_gen = rand_pool.get_state(); - const double vn = vscale * sqrt(-log(rand_gen.drand())); - const double vr = vscale * sqrt(-log(rand_gen.drand())); + const double vn = vscale * sqrt(-log(1.0 - rand_gen.drand())); + const double vr = vscale * sqrt(-log(1.0 - rand_gen.drand())); const double theta1 = MathConst::MY_2PI * rand_gen.drand(); const double theta2 = MathConst::MY_2PI * rand_gen.drand(); diff --git a/src/KOKKOS/fix_ave_histo_weight_kokkos.cpp b/src/KOKKOS/fix_ave_histo_weight_kokkos.cpp index f3275b124..b61d37e32 100644 --- a/src/KOKKOS/fix_ave_histo_weight_kokkos.cpp +++ b/src/KOKKOS/fix_ave_histo_weight_kokkos.cpp @@ -256,7 +256,7 @@ void FixAveHistoWeightKokkos::calculate_weights() if (grid->maxlocal > maxvectorwt) { memoryKK->destroy_kokkos(k_vectorwt,vectorwt); maxvectorwt = grid->maxlocal; - memory->create(vectorwt,maxvectorwt,"ave/histo/weight:vectorwt"); + memoryKK->create_kokkos(k_vectorwt,vectorwt,maxvectorwt,"ave/histo/weight:vectorwt"); } input->variable->compute_grid(m,vectorwt,1,0); k_vectorwt.modify_host(); @@ -268,7 +268,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"); } } @@ -376,7 +375,7 @@ void FixAveHistoWeightKokkos::bin_particles( KokkosBase* regionKKBase = dynamic_cast(region); - if (k_match.extent(0) > nmax) + if (k_match.extent(0) < nmax) MemKK::realloc_kokkos(k_match,"fix_ave_histo_weight:match",nmax); regionKKBase->match_all_kokkos(k_match); diff --git a/src/KOKKOS/fix_emit_face_kokkos.cpp b/src/KOKKOS/fix_emit_face_kokkos.cpp index 3b40812fd..c815e785b 100644 --- a/src/KOKKOS/fix_emit_face_kokkos.cpp +++ b/src/KOKKOS/fix_emit_face_kokkos.cpp @@ -486,7 +486,7 @@ void FixEmitFaceKokkos::operator()(TagFixEmitFace_perform_task, const int &i, in d_beta_un(cand) = beta_un; d_theta(cand) = MY_2PI * rand_gen.drand(); - d_vr(cand) = vscale_val * sqrt(-log(rand_gen.drand())); + d_vr(cand) = vscale_val * sqrt(-log(1.0 - rand_gen.drand())); d_erot(cand) = particle_kk_copy.obj.erot(ispecies,temp_rot,rand_gen); d_evib(cand) = particle_kk_copy.obj.evib(ispecies,temp_vib,rand_gen); d_id(cand) = MAXSMALLINT*rand_gen.drand(); @@ -545,7 +545,7 @@ void FixEmitFaceKokkos::operator()(TagFixEmitFace_perform_task, const int &i, in d_beta_un(cand) = beta_un; d_theta(cand) = MY_2PI * rand_gen.drand(); - d_vr(cand) = vscale_val * sqrt(-log(rand_gen.drand())); + d_vr(cand) = vscale_val * sqrt(-log(1.0 - rand_gen.drand())); d_erot(cand) = particle_kk_copy.obj.erot(ispecies,temp_rot,rand_gen); d_evib(cand) = particle_kk_copy.obj.evib(ispecies,temp_vib,rand_gen); d_id(cand) = MAXSMALLINT*rand_gen.drand(); diff --git a/src/KOKKOS/fix_emit_surf_kokkos.cpp b/src/KOKKOS/fix_emit_surf_kokkos.cpp index b38dbbebc..225f34cc5 100644 --- a/src/KOKKOS/fix_emit_surf_kokkos.cpp +++ b/src/KOKKOS/fix_emit_surf_kokkos.cpp @@ -608,7 +608,7 @@ void FixEmitSurfKokkos::operator()(TagFixEmitSurf_perform_task, const int &i, in else vnmag = beta_un*vscale + indot; const double theta = MY_2PI * rand_gen.drand(); - const double vr = vscale * sqrt(-log(rand_gen.drand())); + const double vr = vscale * sqrt(-log(1.0 - rand_gen.drand())); double vamag,vbmag; if (normalflag) { @@ -720,7 +720,7 @@ void FixEmitSurfKokkos::operator()(TagFixEmitSurf_perform_task, const int &i, in else vnmag = beta_un*vscale + indot; const double theta = MY_2PI * rand_gen.drand(); - const double vr = vscale * sqrt(-log(rand_gen.drand())); + const double vr = vscale * sqrt(-log(1.0 - rand_gen.drand())); double vamag,vbmag; if (normalflag) { diff --git a/src/KOKKOS/fix_grid_check_kokkos.cpp b/src/KOKKOS/fix_grid_check_kokkos.cpp index 587883a73..89975a44d 100644 --- a/src/KOKKOS/fix_grid_check_kokkos.cpp +++ b/src/KOKKOS/fix_grid_check_kokkos.cpp @@ -79,6 +79,7 @@ void FixGridCheckKokkos::end_of_step() if (icell < 0 || icell >= nglocal) { d_particle_problems(i) |= IS_IN_INVALID_CELL; local_nflag++; + return; } // does particle coord match icell bounds @@ -169,7 +170,7 @@ void FixGridCheckKokkos::end_of_step() MPI_Allreduce(&nflag,&all,1,MPI_INT,MPI_SUM,world); if (all && comm->me == 0) { char str[128]; - sprintf(str,"%d particles were in wrong cells on timestep " + snprintf(str, 128,"%d particles were in wrong cells on timestep " BIGINT_FORMAT,all,update->ntimestep); error->warning(FLERR,str); } @@ -183,14 +184,14 @@ void FixGridCheckKokkos::end_of_step() for (int i = 0; i < nlocal; ++i) { 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 + snprintf(str, 128, + "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) { - sprintf(str, + snprintf(str, 128, "Particle %d,%d on proc %d is outside cell " CELLINT_FORMAT " on timestep " BIGINT_FORMAT, i,particles[i].id,comm->me,cells[icell].id, @@ -198,7 +199,7 @@ void FixGridCheckKokkos::end_of_step() error->one(FLERR,str); } if (h_particle_problems(i) & IS_IN_SPLIT_CELL) { - sprintf(str, + snprintf(str, 128, "Particle %d,%d on proc %d is in split cell " CELLINT_FORMAT " on timestep " BIGINT_FORMAT, i,particles[i].id,comm->me,cells[icell].id, @@ -206,14 +207,14 @@ void FixGridCheckKokkos::end_of_step() error->one(FLERR,str); } if (h_particle_problems(i) & IS_IN_INTERIOR_CELL) { - sprintf(str, + snprintf(str, 128, "Particle %d,%d on proc %d is in interior cell " CELLINT_FORMAT " on timestep " BIGINT_FORMAT, i,particles[i].id,comm->me,cells[icell].id,update->ntimestep); error->one(FLERR,str); } if (h_particle_problems(i) & IS_IN_ZERO_VOLUME_CELL) { - sprintf(str, + snprintf(str, 128, "Particle %d,%d on proc %d is in volume=0 cell " CELLINT_FORMAT " on timestep " BIGINT_FORMAT, i,particles[i].id,comm->me,cells[icell].id,update->ntimestep); diff --git a/src/KOKKOS/fix_vibmode_kokkos.h b/src/KOKKOS/fix_vibmode_kokkos.h index 519416b89..d0379f35a 100644 --- a/src/KOKKOS/fix_vibmode_kokkos.h +++ b/src/KOKKOS/fix_vibmode_kokkos.h @@ -96,7 +96,7 @@ void FixVibmodeKokkos::update_custom_kokkos(int index, double temp_thermal, rand_type rand_gen = rand_pool.get_state(); for (int imode = 0; imode < nmode; imode++) { - const int ivib = static_cast (-log(rand_gen.drand()) * temp_vib / + const int ivib = static_cast (-log(1.0 - rand_gen.drand()) * temp_vib / d_species[isp].vibtemp[imode]); d_vibmode(index,imode) = ivib; evib += ivib * boltz * d_species[isp].vibtemp[imode]; 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/particle_kokkos.h b/src/KOKKOS/particle_kokkos.h index f7013f1ff..f5a69c200 100644 --- a/src/KOKKOS/particle_kokkos.h +++ b/src/KOKKOS/particle_kokkos.h @@ -259,7 +259,7 @@ double ParticleKokkos::erot(int isp, double temp_thermal, rand_type &erandom) co if (d_species[isp].rotdof < 2) return 0.0; if (d_species[isp].rotdof == 2) - eng = -log(erandom.drand()) * boltz * temp_thermal; + eng = -log(1.0 - erandom.drand()) * boltz * temp_thermal; else { a = 0.5*d_species[isp].rotdof-1.0; while (1) { @@ -289,12 +289,12 @@ double ParticleKokkos::evib(int isp, double temp_thermal, rand_type &erandom) co eng = 0.0; if (vibstyle == DISCRETE && d_species[isp].vibdof == 2) { - int ivib = static_cast (-log(erandom.drand()) * temp_thermal / + int ivib = static_cast (-log(1.0 - erandom.drand()) * temp_thermal / d_species[isp].vibtemp[0]); eng = ivib * boltz * d_species[isp].vibtemp[0]; } else if (vibstyle == SMOOTH || d_species[isp].vibdof >= 2) { if (d_species[isp].vibdof == 2) - eng = -log(erandom.drand()) * boltz * temp_thermal; + eng = -log(1.0 - erandom.drand()) * boltz * temp_thermal; else if (d_species[isp].vibdof > 2) { a = 0.5*d_species[isp].vibdof-1.; while (1) { diff --git a/src/KOKKOS/react_bird_kokkos.cpp b/src/KOKKOS/react_bird_kokkos.cpp index 2a46dc971..e7ddced1d 100644 --- a/src/KOKKOS/react_bird_kokkos.cpp +++ b/src/KOKKOS/react_bird_kokkos.cpp @@ -44,7 +44,7 @@ enum{ARRHENIUS,QUANTUM}; // other react files ReactBirdKokkos::ReactBirdKokkos(SPARTA *sparta, int narg, char **arg) : ReactBird(sparta, narg, arg), - rand_pool(12345 + comm->me + rand_pool(54321 + comm->me #ifdef SPARTA_KOKKOS_EXACT , sparta #endif diff --git a/src/KOKKOS/react_tce_kokkos.h b/src/KOKKOS/react_tce_kokkos.h index bc960f624..11322ac64 100644 --- a/src/KOKKOS/react_tce_kokkos.h +++ b/src/KOKKOS/react_tce_kokkos.h @@ -46,7 +46,7 @@ double bird_Evib(const int& nmode, const double& Tvib, // Comutes f for Newton's search method outlined in newtonTvib() double f = -Evib; - const double kb = 1.38064852e-23; + const double kb = boltz; for (int i = 0; i < nmode; i++) { const double vti = vibtemp[i]; @@ -64,7 +64,7 @@ double bird_dEvib(const int& nmode, const double& Tvib, const double vibtemp[]) // Comutes df for Newton's search method double df = 0.0; - const double kb = 1.38064852e-23; + const double kb = boltz; for (int i = 0; i < nmode; i++) { const double vti = vibtemp[i]; @@ -192,7 +192,7 @@ int attempt_kk(Particle::OnePart *ip, Particle::OnePart *jp, if (r->d_coeff[1]>((-1)*r->d_coeff[4])) e_excess = ecc - r->d_coeff[1]; else e_excess = ecc + r->d_coeff[4]; - if (e_excess <= 0.0) continue; + if (e_excess <= 0.0 || ecc <= 0.0) continue; if (!partialEnergy) { diff --git a/src/KOKKOS/remap3d_kokkos.cpp b/src/KOKKOS/remap3d_kokkos.cpp index eabffe36d..3d89063b6 100644 --- a/src/KOKKOS/remap3d_kokkos.cpp +++ b/src/KOKKOS/remap3d_kokkos.cpp @@ -324,10 +324,10 @@ struct remap_plan_3d_kokkos* RemapKokkos3d::remap_3d_cre out.ksize = out.khi - out.klo + 1; inarray = (struct extent_3d *) malloc(nprocs*sizeof(struct extent_3d)); - if (inarray == nullptr) return nullptr; + if (inarray == nullptr) { delete plan; return nullptr; } outarray = (struct extent_3d *) malloc(nprocs*sizeof(struct extent_3d)); - if (outarray == nullptr) return nullptr; + if (outarray == nullptr) { free(inarray); delete plan; return nullptr; } // combine input & output extents across all procs @@ -360,7 +360,10 @@ struct remap_plan_3d_kokkos* RemapKokkos3d::remap_3d_cre malloc(nsend*sizeof(struct pack_plan_3d)); if (plan->send_offset == nullptr || plan->send_size == nullptr || - plan->send_proc == nullptr || plan->packplan == nullptr) return nullptr; + plan->send_proc == nullptr || plan->packplan == nullptr) { + free(plan->send_offset); free(plan->send_size); free(plan->send_proc); free(plan->packplan); + free(outarray); free(inarray); delete plan; return nullptr; + } } if (nrecv) { @@ -393,7 +396,12 @@ struct remap_plan_3d_kokkos* RemapKokkos3d::remap_3d_cre if (plan->recv_offset == nullptr || plan->recv_size == nullptr || plan->recv_proc == nullptr || plan->recv_bufloc == nullptr || - plan->request == nullptr || plan->unpackplan == nullptr) return nullptr; + plan->request == nullptr || plan->unpackplan == nullptr) { + free(plan->recv_offset); free(plan->recv_size); free(plan->recv_proc); + free(plan->recv_bufloc); free(plan->request); free(plan->unpackplan); + if (nsend) { free(plan->send_offset); free(plan->send_size); free(plan->send_proc); free(plan->packplan); } + free(outarray); free(inarray); delete plan; return nullptr; + } } // store send info, with self as last entry @@ -529,14 +537,14 @@ struct remap_plan_3d_kokkos* RemapKokkos3d::remap_3d_cre // start with max length -- nprocs. Unused entries will be removed later int *local_cnts = (int*) malloc(2*nprocs*sizeof(int)); - if (local_cnts == nullptr) return nullptr; + if (local_cnts == nullptr) { free(outarray); free(inarray); delete plan; return nullptr; } int *local_sendcnts = local_cnts; int *local_recvcnts = (local_cnts + nprocs); // local arrays used to store the results of the allreduce int *global_cnts = (int*) malloc(2*nprocs*sizeof(int)); - if (global_cnts == nullptr) return nullptr; + if (global_cnts == nullptr) { free(local_cnts); free(outarray); free(inarray); delete plan; return nullptr; } int *global_sendcnts = global_cnts; int *global_recvcnts = (global_cnts + nprocs); @@ -615,7 +623,11 @@ struct remap_plan_3d_kokkos* RemapKokkos3d::remap_3d_cre if (plan->send_offset == nullptr || plan->send_size == nullptr || plan->sendcnts == nullptr || plan->sdispls == nullptr || - plan->packplan == nullptr) return nullptr; + plan->packplan == nullptr) { + free(plan->send_offset); free(plan->send_size); free(plan->sendcnts); + free(plan->sdispls); free(plan->packplan); + free(plan->commringlist); free(outarray); free(inarray); delete plan; return nullptr; + } // recv space @@ -653,7 +665,11 @@ struct remap_plan_3d_kokkos* RemapKokkos3d::remap_3d_cre if (plan->recv_offset == nullptr || plan->recv_size == nullptr || plan->rcvcnts == nullptr || plan->rdispls == nullptr || - plan->unpackplan == nullptr) return nullptr; + plan->unpackplan == nullptr) { + free(plan->recv_offset); free(plan->recv_size); free(plan->rcvcnts); free(plan->rdispls); free(plan->unpackplan); + if (nsend || nrecv) { free(plan->send_offset); free(plan->send_size); free(plan->sendcnts); free(plan->sdispls); free(plan->packplan); } + free(plan->commringlist); free(outarray); free(inarray); delete plan; return nullptr; + } } // store send info, with self as last entry diff --git a/src/KOKKOS/surf_collide_diffuse_kokkos.h b/src/KOKKOS/surf_collide_diffuse_kokkos.h index 25826164a..999cc1e59 100644 --- a/src/KOKKOS/surf_collide_diffuse_kokkos.h +++ b/src/KOKKOS/surf_collide_diffuse_kokkos.h @@ -228,10 +228,10 @@ class SurfCollideDiffuseKokkos : public SurfCollideDiffuse { int ispecies = p->ispecies; double vrm = sqrt(2.0*boltz * twall / d_species[ispecies].mass); - double vperp = vrm * sqrt(-log(rand_gen.drand())); + double vperp = vrm * sqrt(-log(1.0 - rand_gen.drand())); double theta = MathConst::MY_2PI * rand_gen.drand(); - double vtangent = vrm * sqrt(-log(rand_gen.drand())); + double vtangent = vrm * sqrt(-log(1.0 - rand_gen.drand())); double vtan1 = vtangent * sin(theta); double vtan2 = vtangent * cos(theta); @@ -323,11 +323,11 @@ class SurfCollideDiffuseKokkos : public SurfCollideDiffuse { if (d_species[isp].rotdof < 2) return 0.0; if (rotstyle == DISCRETE && d_species[isp].rotdof == 2) { - int irot = -log(rand_gen.drand()) * temp_thermal / + int irot = -log(1.0 - rand_gen.drand()) * temp_thermal / d_species[isp].rottemp[0]; eng = irot * boltz * d_species[isp].rottemp[0]; } else if (rotstyle == SMOOTH && d_species[isp].rotdof == 2) { - eng = -log(rand_gen.drand()) * boltz * temp_thermal; + eng = -log(1.0 - rand_gen.drand()) * boltz * temp_thermal; } else { a = 0.5*d_species[isp].rotdof-1.0; while (1) { @@ -362,12 +362,12 @@ class SurfCollideDiffuseKokkos : public SurfCollideDiffuse { eng = 0.0; if (vibstyle == DISCRETE && d_species[isp].vibdof == 2) { - int ivib = -log(rand_gen.drand()) * temp_thermal / + int ivib = -log(1.0 - rand_gen.drand()) * temp_thermal / d_species[isp].vibtemp[0]; eng = ivib * boltz * d_species[isp].vibtemp[0]; } else if (vibstyle == SMOOTH || d_species[isp].vibdof >= 2) { if (d_species[isp].vibdof == 2) - eng = -log(rand_gen.drand()) * boltz * temp_thermal; + eng = -log(1.0 - rand_gen.drand()) * boltz * temp_thermal; else if (d_species[isp].vibdof > 2) { a = 0.5*d_species[isp].vibdof-1.; while (1) { diff --git a/src/KOKKOS/surf_collide_specular_kokkos.cpp b/src/KOKKOS/surf_collide_specular_kokkos.cpp index 1191828ad..2ee7300b1 100644 --- a/src/KOKKOS/surf_collide_specular_kokkos.cpp +++ b/src/KOKKOS/surf_collide_specular_kokkos.cpp @@ -126,7 +126,7 @@ void SurfCollideSpecularKokkos::pre_collide() sr_kk_global_copy[nglob].copy((SurfReactGlobalKokkos*)(surf->sr[n])); sr_kk_global_copy[nglob].obj.pre_react(); sr_type_list[n] = 0; - sr_map[n] = nprob; + sr_map[n] = nglob; nglob++; } else if (strcmp(surf->sr[n]->style,"prob") == 0) { sr_kk_prob_copy[nprob].copy((SurfReactProbKokkos*)(surf->sr[n])); diff --git a/src/KOKKOS/update_kokkos.cpp b/src/KOKKOS/update_kokkos.cpp index 5f82f5a0f..23d32f320 100644 --- a/src/KOKKOS/update_kokkos.cpp +++ b/src/KOKKOS/update_kokkos.cpp @@ -705,7 +705,7 @@ template < int DIM, int SURF, int REACT, int OPT > void UpdateKokkos::move() if (error_flag) { char str[128]; - sprintf(str, + snprintf(str, sizeof(str), "Particle being sent to self proc " "on step " BIGINT_FORMAT, update->ntimestep); @@ -1053,22 +1053,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; @@ -1111,13 +1123,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; @@ -1419,7 +1437,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 7ac60fa88..b6e7ca8f5 100644 --- a/src/KOKKOS/update_kokkos.h +++ b/src/KOKKOS/update_kokkos.h @@ -209,12 +209,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/adapt_grid.cpp b/src/adapt_grid.cpp index 633c72bef..93d13484f 100644 --- a/src/adapt_grid.cpp +++ b/src/adapt_grid.cpp @@ -984,7 +984,7 @@ void AdaptGrid::candidates_coarsen() if (clhash->find(parentID) == clhash->end()) { if (cnum == cnummax) { cnummax += DELTA_LIST; - clist = (CList *) memory->srealloc(clist,cnummax*sizeof(CList), + clist = (CList *) memory->srealloc(clist,(bigint)cnummax*sizeof(CList), "adapt_grid:clist"); } (*clhash)[parentID] = cnum; @@ -1050,7 +1050,7 @@ void AdaptGrid::candidates_coarsen() (*clhash)[parentID] = cnum; if (cnum == cnummax) { cnummax += DELTA_LIST; - clist = (CList *) memory->srealloc(clist,cnummax*sizeof(CList), + clist = (CList *) memory->srealloc(clist,(bigint)cnummax*sizeof(CList), "adapt_grid:clist"); } clist[cnum].parentID = parentID; @@ -1376,7 +1376,7 @@ void AdaptGrid::particle_surf_comm() nsend = nreturn; memory->create(proclist,nsend,"adapt_grid:proclist"); - SendAdapt *sadapt = (SendAdapt *) memory->smalloc(nsend*sizeof(SendAdapt), + SendAdapt *sadapt = (SendAdapt *) memory->smalloc((bigint)nsend*sizeof(SendAdapt), "adapt_grid:sadapt"); for (int i = 0; i < nreturn; i++) { icell = outbuf[i].icell; @@ -1439,7 +1439,7 @@ void AdaptGrid::particle_surf_comm() if (alhash->find(parentID) == alhash->end()) { if (anum == anummax) { anummax += DELTA_LIST; - alist = (ActionList *) memory->srealloc(alist,anummax*sizeof(ActionList), + alist = (ActionList *) memory->srealloc(alist,(bigint)anummax*sizeof(ActionList), "adapt_grid:alist"); } (*alhash)[parentID] = anum; @@ -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 5e1567b53..f9569a426 100644 --- a/src/collide_vss.cpp +++ b/src/collide_vss.cpp @@ -139,15 +139,17 @@ double CollideVSS::attempt_collision(int icell, int np, double volume) double fnum = update->fnum; double dt = update->dt; - double nattempt; + double nattempt = 0.0; - if (remainflag) { - nattempt = 0.5 * np * (np-1) * - vremax[icell][0][0] * dt * fnum / volume + remain[icell][0][0]; - remain[icell][0][0] = nattempt - static_cast (nattempt); - } else { - nattempt = 0.5 * np * (np-1) * - vremax[icell][0][0] * dt * fnum / volume + random->uniform(); + if (volume > 0.0) { + if (remainflag) { + nattempt = 0.5 * np * (np-1) * + vremax[icell][0][0] * dt * fnum / volume + remain[icell][0][0]; + remain[icell][0][0] = nattempt - static_cast (nattempt); + } else { + nattempt = 0.5 * np * (np-1) * + vremax[icell][0][0] * dt * fnum / volume + random->uniform(); + } } return nattempt; @@ -161,21 +163,23 @@ double CollideVSS::attempt_collision(int icell, int igroup, int jgroup, double fnum = update->fnum; double dt = update->dt; - double nattempt; + double nattempt = 0.0; - // return 2x the value for igroup != jgroup, since no J,I pairing + if (volume > 0.0) { + // return 2x the value for igroup != jgroup, since no J,I pairing - double npairs; - if (igroup == jgroup) npairs = 0.5 * ngroup[igroup] * (ngroup[igroup]-1); - else npairs = ngroup[igroup] * (ngroup[jgroup]); - //else npairs = 0.5 * ngroup[igroup] * (ngroup[jgroup]); + double npairs; + if (igroup == jgroup) npairs = 0.5 * ngroup[igroup] * (ngroup[igroup]-1); + else npairs = ngroup[igroup] * (ngroup[jgroup]); + //else npairs = 0.5 * ngroup[igroup] * (ngroup[jgroup]); - nattempt = npairs * vremax[icell][igroup][jgroup] * dt * fnum / volume; + nattempt = npairs * vremax[icell][igroup][jgroup] * dt * fnum / volume; - if (remainflag) { - nattempt += remain[icell][igroup][jgroup]; - remain[icell][igroup][jgroup] = nattempt - static_cast (nattempt); - } else nattempt += random->uniform(); + if (remainflag) { + nattempt += remain[icell][igroup][jgroup]; + remain[icell][igroup][jgroup] = nattempt - static_cast (nattempt); + } else nattempt += random->uniform(); + } return nattempt; } @@ -210,6 +214,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; @@ -899,7 +904,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/comm.cpp b/src/comm.cpp index d0049791b..c535a5233 100644 --- a/src/comm.cpp +++ b/src/comm.cpp @@ -546,7 +546,6 @@ int Comm::send_cells_adapt(int nsend, int *procsend, char *inbuf, char **outbuf) memory->sfree(rbuf); maxrecvbuf = recvsize; rbuf = (char *) memory->smalloc(maxrecvbuf,"comm:rbuf"); - memory->create(rbuf,maxrecvbuf,"comm:rbuf"); memset(rbuf,0,maxrecvbuf); } 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_dt_grid.cpp b/src/compute_dt_grid.cpp index f6c236e19..6ee030374 100644 --- a/src/compute_dt_grid.cpp +++ b/src/compute_dt_grid.cpp @@ -643,13 +643,15 @@ void ComputeDtGrid::compute_per_grid() // cell dt based on transit time using maximum most probable speed vrm_max = sqrt(2.0*update->boltz * temp[i] / min_species_mass); - dt_candidate = transit_fraction*dx/vrm_max; - cell_dt_desired = MIN(dt_candidate,cell_dt_desired); - dt_candidate = transit_fraction*dy/vrm_max; - cell_dt_desired = MIN(dt_candidate,cell_dt_desired); - if (domain->dimension == 3) { - dt_candidate = transit_fraction*dz/vrm_max; + if (vrm_max > 0.0) { + dt_candidate = transit_fraction*dx/vrm_max; + cell_dt_desired = MIN(dt_candidate,cell_dt_desired); + dt_candidate = transit_fraction*dy/vrm_max; cell_dt_desired = MIN(dt_candidate,cell_dt_desired); + if (domain->dimension == 3) { + dt_candidate = transit_fraction*dz/vrm_max; + cell_dt_desired = MIN(dt_candidate,cell_dt_desired); + } } // per grid cell timestep = final cell_dt_desired for all criteria diff --git a/src/compute_eflux_grid.cpp b/src/compute_eflux_grid.cpp index 01de6b05b..9fe5a589e 100644 --- a/src/compute_eflux_grid.cpp +++ b/src/compute_eflux_grid.cpp @@ -365,8 +365,10 @@ void ComputeEFluxGrid::post_process_grid(int index, int nsample, 2.0*t[mv]*t[mv1]*t[mv1]/summass/summass; h2 = t[mvv2v2] - 2.0*t[mvv2]*t[mv2]/summass - t[mv]*t[mv2v2]/summass + 2.0*t[mv]*t[mv2]*t[mv2]/summass/summass; - wt = 0.5 * fnum * cinfo[icell].weight / cinfo[icell].volume; - vec[k] = wt/nsample * (h + h1 + h2); + if (cinfo[icell].volume > 0.0) { + wt = 0.5 * fnum * cinfo[icell].weight / cinfo[icell].volume; + vec[k] = wt/nsample * (h + h1 + h2); + } else vec[k] = 0.0; } k += nstride; } @@ -440,8 +442,8 @@ void ComputeEFluxGrid::reallocate() bigint ComputeEFluxGrid::memory_usage() { - bigint bytes; - bytes = nglocal * sizeof(double); - bytes = ntotal*nglocal * sizeof(double); + bigint bytes = 0; + bytes += nglocal * sizeof(double); + bytes += ntotal*nglocal * sizeof(double); return bytes; } diff --git a/src/compute_gas_collision_grid.cpp b/src/compute_gas_collision_grid.cpp index 659aee6a1..3e9971ed5 100644 --- a/src/compute_gas_collision_grid.cpp +++ b/src/compute_gas_collision_grid.cpp @@ -83,7 +83,7 @@ void ComputeGasCollisionGrid::compute_per_grid() void ComputeGasCollisionGrid::clear() { cinfo = grid->cinfo; - memset(vector_grid,0,nglocal*sizeof(double)); + memset(vector_grid,0,((size_t)nglocal)*sizeof(double)); } /* ---------------------------------------------------------------------- @@ -143,7 +143,7 @@ void ComputeGasCollisionGrid::reallocate() // also note if load-balancing is done, tallies will be lost // would need to implement (un)pack_grid_one() to avoid this - memset(vector_grid,0,nglocal*sizeof(double)); + memset(vector_grid,0,((size_t)nglocal)*sizeof(double)); } /* ---------------------------------------------------------------------- diff --git a/src/compute_gas_collision_tally.cpp b/src/compute_gas_collision_tally.cpp index 6c5f75e9e..d7afa2a19 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_grid.cpp b/src/compute_gas_reaction_grid.cpp index 216237991..d7241ce76 100644 --- a/src/compute_gas_reaction_grid.cpp +++ b/src/compute_gas_reaction_grid.cpp @@ -134,8 +134,8 @@ void ComputeGasReactionGrid::compute_per_grid() void ComputeGasReactionGrid::clear() { cinfo = grid->cinfo; - if (ncol == 0) memset(vector_grid,0,nglocal*sizeof(double)); - else if (nglocal) memset(&array_grid[0][0],0,nglocal*ncol*sizeof(double)); + if (ncol == 0) memset(vector_grid,0,((size_t)nglocal)*sizeof(double)); + else if (nglocal) memset(&array_grid[0][0],0,((size_t)nglocal)*ncol*sizeof(double)); } /* ---------------------------------------------------------------------- @@ -207,8 +207,8 @@ void ComputeGasReactionGrid::reallocate() // also note if load-balancing is done, tallies will be lost // would need to implement (un)pack_grid_one() to avoid this - if (ncol == 0) memset(vector_grid,0,nglocal*sizeof(double)); - else if (nglocal) memset(&array_grid[0][0],0,nglocal*ncol*sizeof(double)); + if (ncol == 0) memset(vector_grid,0,((size_t)nglocal)*sizeof(double)); + else if (nglocal) memset(&array_grid[0][0],0,((size_t)nglocal)*ncol*sizeof(double)); } /* ---------------------------------------------------------------------- diff --git a/src/compute_gas_reaction_tally.cpp b/src/compute_gas_reaction_tally.cpp index cf9303b29..cb8faad05 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_grid.cpp b/src/compute_grid.cpp index a7ac9a6bc..96807c46a 100644 --- a/src/compute_grid.cpp +++ b/src/compute_grid.cpp @@ -635,8 +635,8 @@ void ComputeGrid::reallocate() bigint ComputeGrid::memory_usage() { - bigint bytes; - bytes = nglocal * sizeof(double); - bytes = ntotal*nglocal * sizeof(double); + bigint bytes = 0; + bytes += nglocal * sizeof(double); + bytes += ntotal*nglocal * sizeof(double); return bytes; } diff --git a/src/compute_lambda_grid.cpp b/src/compute_lambda_grid.cpp index 728011d60..715e8d584 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; @@ -664,23 +666,43 @@ void ComputeLambdaGrid::compute_per_grid() sizeall += sizez; sizeall /= 3.0; } - if (noutputs == 1) vector_grid[i] = lambda / sizeall; - else array_grid[i][output_order[KNALL]] = lambda / sizeall; + if (sizeall > 0.0) { + if (noutputs == 1) vector_grid[i] = lambda / sizeall; + else array_grid[i][output_order[KNALL]] = lambda / sizeall; + } else { + if (noutputs == 1) vector_grid[i] = 0.0; + else array_grid[i][output_order[KNALL]] = 0.0; + } } if (knxflag) { - if (noutputs == 1) vector_grid[i] = lambda / sizex; - else array_grid[i][output_order[KNX]] = lambda / sizex; + if (sizex > 0.0) { + if (noutputs == 1) vector_grid[i] = lambda / sizex; + else array_grid[i][output_order[KNX]] = lambda / sizex; + } else { + if (noutputs == 1) vector_grid[i] = 0.0; + else array_grid[i][output_order[KNX]] = 0.0; + } } if (knyflag) { - if (noutputs == 1) vector_grid[i] = lambda / sizey; - array_grid[i][output_order[KNY]] = lambda / sizey; + if (sizey > 0.0) { + if (noutputs == 1) vector_grid[i] = lambda / sizey; + else array_grid[i][output_order[KNY]] = lambda / sizey; + } else { + if (noutputs == 1) vector_grid[i] = 0.0; + else array_grid[i][output_order[KNY]] = 0.0; + } } if (knzflag) { - if (noutputs == 1) vector_grid[i] = lambda / sizez; - array_grid[i][output_order[KNZ]] = lambda / sizez; + if (sizez > 0.0) { + if (noutputs == 1) vector_grid[i] = lambda / sizez; + else array_grid[i][output_order[KNZ]] = lambda / sizez; + } else { + if (noutputs == 1) vector_grid[i] = 0.0; + else array_grid[i][output_order[KNZ]] = 0.0; + } } } } @@ -736,10 +758,10 @@ void ComputeLambdaGrid::reallocate() bigint ComputeLambdaGrid::memory_usage() { - bigint bytes; - bytes = nglocal * sizeof(double); // vector_grid + bigint bytes = 0; + bytes += nglocal * sizeof(double); // vector_grid if (nrho_values > 1) - bytes = nglocal * nrho_values * sizeof(double); // array_grid1 + bytes += nglocal * nrho_values * sizeof(double); // array_grid1 bytes += nglocal * noutputs * sizeof(double); // array_grid bytes += nglocal * sizeof(double); // lambda_grid bytes += 2 * nglocal * ntotal * sizeof(double); // lambdainv + tauinv diff --git a/src/compute_pflux_grid.cpp b/src/compute_pflux_grid.cpp index fb7e79386..fc3d39e0d 100644 --- a/src/compute_pflux_grid.cpp +++ b/src/compute_pflux_grid.cpp @@ -303,7 +303,7 @@ void ComputePFluxGrid::post_process_grid(int index, int nsample, int mvv = emap[2]; for (int icell = lo; icell < hi; icell++) { summass = etally[icell][mass]; - if (summass == 0.0) vec[k] = 0.0; + if (summass == 0.0 || cinfo[icell].volume == 0.0) vec[k] = 0.0; else { wt = fnum * cinfo[icell].weight / cinfo[icell].volume; summv = etally[icell][mv]; @@ -328,7 +328,7 @@ void ComputePFluxGrid::post_process_grid(int index, int nsample, int mvv = emap[3]; for (int icell = lo; icell < hi; icell++) { summass = etally[icell][mass]; - if (summass == 0.0) vec[k] = 0.0; + if (summass == 0.0 || cinfo[icell].volume == 0.0) vec[k] = 0.0; else { wt = fnum * cinfo[icell].weight / cinfo[icell].volume; vec[k] = wt/nsample * (etally[icell][mvv] - @@ -409,8 +409,8 @@ void ComputePFluxGrid::reallocate() bigint ComputePFluxGrid::memory_usage() { - bigint bytes; - bytes = nglocal * sizeof(double); - bytes = ntotal*nglocal * sizeof(double); + bigint bytes = 0; + bytes += nglocal * sizeof(double); + bytes += ntotal*nglocal * sizeof(double); return bytes; } diff --git a/src/compute_property_surf.cpp b/src/compute_property_surf.cpp index 66a710d45..190731e3b 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 0bb4fd5c8..10b95c8bf 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 4b49c3e45..56ba57e49 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 c0d71adf0..f2272d3ca 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_sonine_grid.cpp b/src/compute_sonine_grid.cpp index 53f92a051..6425d88db 100644 --- a/src/compute_sonine_grid.cpp +++ b/src/compute_sonine_grid.cpp @@ -338,9 +338,9 @@ void ComputeSonineGrid::reallocate() bigint ComputeSonineGrid::memory_usage() { - bigint bytes; - bytes = nglocal * sizeof(double); // vector_grid - bytes = ntotal*nglocal * sizeof(double); // tally array + bigint bytes = 0; + bytes += nglocal * sizeof(double); // vector_grid + bytes += ntotal*nglocal * sizeof(double); // tally array bytes += nglocal*ngroup*4 * sizeof(double); // vcom return bytes; } diff --git a/src/compute_surf_reaction_tally.cpp b/src/compute_surf_reaction_tally.cpp index 87eb47ade..22da06f59 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_thermal_grid.cpp b/src/compute_thermal_grid.cpp index a1a0333c6..2ac8e6cfa 100644 --- a/src/compute_thermal_grid.cpp +++ b/src/compute_thermal_grid.cpp @@ -247,7 +247,11 @@ post_process_grid(int index, int nsample, vec[k] = mvsq - (mvx*mvx + mvy*mvy + mvz*mvz)/mass; vec[k] *= prefactor; if (tflag) vec[k] /= ncount; - else vec[k] *= cinfo[icell].weight / cinfo[icell].volume / nsample; + else { + if (cinfo[icell].volume > 0.0) + vec[k] *= cinfo[icell].weight / cinfo[icell].volume / nsample; + else vec[k] = 0.0; + } } k += nstride; } @@ -276,7 +280,7 @@ void ComputeThermalGrid::reallocate() bigint ComputeThermalGrid::memory_usage() { bigint bytes = 0; - bytes = nglocal * sizeof(double); - bytes = ntotal*nglocal * sizeof(double); + bytes += nglocal * sizeof(double); + bytes += ntotal*nglocal * sizeof(double); return bytes; } diff --git a/src/compute_tvib_grid.cpp b/src/compute_tvib_grid.cpp index aa4f5f37a..ed241fbeb 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; @@ -561,8 +561,8 @@ void ComputeTvibGrid::reallocate() bigint ComputeTvibGrid::memory_usage() { - bigint bytes; - bytes = nglocal * sizeof(double); - bytes = ntally*nglocal * sizeof(double); + bigint bytes = 0; + bytes += nglocal * sizeof(double); + bytes += ntally*nglocal * sizeof(double); return bytes; } diff --git a/src/create_isurf.cpp b/src/create_isurf.cpp index 9120ffcbf..cb336c3b2 100644 --- a/src/create_isurf.cpp +++ b/src/create_isurf.cpp @@ -988,9 +988,9 @@ void CreateISurf::comm_neigh_corners(int which) else if (which == CVAL) ncomm = 1 + ncorner; else if (which == INVAL) ncomm = 1 + ncorner*nmulti; - if (nsend*ncomm > maxsbuf) { + if ((bigint)nsend*ncomm > maxsbuf) { memory->destroy(sbuf); - maxsbuf = nsend*ncomm; + maxsbuf = (bigint)nsend*ncomm; memory->create(sbuf,maxsbuf,"createisurf:sbuf"); } @@ -1551,7 +1551,10 @@ void CreateISurf::set_cvalues_voxel() dx = cells[icell].hi[0] - cells[icell].lo[0]; dy = cells[icell].hi[1] - cells[icell].lo[1]; dz = cells[icell].hi[2] - cells[icell].lo[2]; - sfrac = (dx*dy*dz - cvol) / (dx*dy*dz); + + double boxvol = dx*dy*dz; + if (boxvol > 0.0) sfrac = (boxvol - cvol) / boxvol; + else sfrac = 0.0; if (sfrac < 0.0 || sfrac > 1.0) error->one(FLERR,"Calculated solid fraction above one or negative"); @@ -1848,6 +1851,7 @@ double CreateISurf::param2cval(double param, double v1) // param = (thresh - v0) / (v1 - v0) double v0; + if (param >= 1.0) param = 0.999999; v0 = (thresh - v1*param) / (1.0 - param); // bound by limits @@ -1887,11 +1891,11 @@ void CreateISurf::remove_old() else nbytes = sizeof(Surf::Tri); if (dim == 2) { - llines = (Surf::Line *) memory->smalloc(nsurf*nbytes,"createisurf:lines"); - memcpy(llines,surf->mylines,nsurf*nbytes); + llines = (Surf::Line *) memory->smalloc((bigint)nsurf*nbytes,"createisurf:lines"); + memcpy(llines,surf->mylines,(size_t)nsurf*nbytes); } else { - ltris = (Surf::Tri *) memory->smalloc(nsurf*nbytes,"createisurf:ltris"); - memcpy(ltris,surf->mytris,nsurf*nbytes); + ltris = (Surf::Tri *) memory->smalloc((bigint)nsurf*nbytes,"createisurf:ltris"); + memcpy(ltris,surf->mytris,(size_t)nsurf*nbytes); } surf->add_surfs(1,0,llines,ltris,ncustom,index_custom,cuvalues); diff --git a/src/create_isurf.h b/src/create_isurf.h index a23edc35a..63045f033 100644 --- a/src/create_isurf.h +++ b/src/create_isurf.h @@ -95,7 +95,7 @@ class CreateISurf : protected Pointers { int maxsend; double *sbuf; - int maxsbuf; + bigint maxsbuf; union ubuf { double d; diff --git a/src/custom.cpp b/src/custom.cpp index 42ad079c1..1c261ecf6 100644 --- a/src/custom.cpp +++ b/src/custom.cpp @@ -56,7 +56,7 @@ Custom::~Custom() for (int i = 0; i < naction; i++) { int action = actions[i].action; - if (action == FILESTYLE) { + if (action == FILESTYLE || action == FILECOARSE) { delete [] actions[i].fname; delete [] actions[i].cindex_file; delete [] actions[i].ctype_file; diff --git a/src/dump.cpp b/src/dump.cpp index 6c4deee15..1aeaa51ab 100644 --- a/src/dump.cpp +++ b/src/dump.cpp @@ -435,7 +435,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 6e513021e..2acea4d98 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 c507d8a4c..9b935efb9 100644 --- a/src/fix_ablate.cpp +++ b/src/fix_ablate.cpp @@ -98,8 +98,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; @@ -113,7 +115,7 @@ FixAblate::FixAblate(SPARTA *sparta, int narg, char **arg) : which = VARIABLE; int n = strlen(arg[5]); - char *idsource = new char[n]; + idsource = new char[n]; strcpy(idsource,&arg[5][2]); } else if (strcmp(arg[5],"random") == 0) { diff --git a/src/fix_ablate_multi_inner.cpp b/src/fix_ablate_multi_inner.cpp index 76945c2f6..2a591ce13 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 @@ -677,7 +678,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 8d719364c..9681dd803 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 2ede68c72..cb6d70016 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 d7640c26f..526202ab0 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 d9f60e06a..9d1882ad8 100644 --- a/src/fix_emit_face.cpp +++ b/src/fix_emit_face.cpp @@ -962,7 +962,7 @@ void FixEmitFace::subsonic_sort() maxactive = grid->nlocal; memory->create(activecell,maxactive,"emit/face:active"); } - memset(activecell,0,maxactive*sizeof(int)); + memset(activecell,0,((size_t)maxactive)*sizeof(int)); for (i = 0; i < ntask; i++) activecell[tasks[i].pcell] = 1; active_current = 1; } @@ -1072,7 +1072,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 * @@ -1111,7 +1111,7 @@ void FixEmitFace::grow_task() // set all new task bytes to 0 so valgrind won't complain // if bytes between fields are uninitialized - memset(&tasks[oldmax],0,(ntaskmax-oldmax)*sizeof(Task)); + memset(&tasks[oldmax],0,((size_t)ntaskmax-oldmax)*sizeof(Task)); // allocate vectors in each new task or set to NULL diff --git a/src/fix_emit_face_file.cpp b/src/fix_emit_face_file.cpp index dfc24c19f..864d1a634 100644 --- a/src/fix_emit_face_file.cpp +++ b/src/fix_emit_face_file.cpp @@ -239,6 +239,8 @@ void FixEmitFaceFile::init() // per-species vectors for mesh setting of species fractions // initialize to mixture settings + delete [] fflag; + delete [] fuser; fflag = new int[nspecies]; fuser = new double[nspecies]; for (isp = 0; isp < nspecies; isp++) { @@ -473,7 +475,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]; @@ -1116,7 +1118,7 @@ void FixEmitFaceFile::subsonic_sort() maxactive = grid->nlocal; memory->create(activecell,maxactive,"emit/face:active"); } - memset(activecell,0,maxactive*sizeof(int)); + memset(activecell,0,((size_t)maxactive)*sizeof(int)); for (i = 0; i < ntask; i++) activecell[tasks[i].pcell] = 1; active_current = 1; } @@ -1225,14 +1227,14 @@ 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); } vscale = tasks[i].vscale; - for (m = 0; m < nspecies; i++) { + for (m = 0; m < nspecies; m++) { ispecies = particle->mixture[imix]->species[m]; vscale[m] = sqrt(2.0 * update->boltz * temp_thermal_cell / species[ispecies].mass); @@ -1262,7 +1264,7 @@ void FixEmitFaceFile::grow_task() // set all new task bytes to 0 so valgrind won't complain // if bytes between fields are uninitialized - memset(&tasks[oldmax],0,(ntaskmax-oldmax)*sizeof(Task)); + memset(&tasks[oldmax],0,((size_t)ntaskmax-oldmax)*sizeof(Task)); // allocate vectors in each new task or set to NULL diff --git a/src/fix_emit_surf.cpp b/src/fix_emit_surf.cpp index 7edf802c1..40dc31492 100644 --- a/src/fix_emit_surf.cpp +++ b/src/fix_emit_surf.cpp @@ -1330,7 +1330,7 @@ void FixEmitSurf::subsonic_sort() maxactive = grid->nlocal; memory->create(activecell,maxactive,"emit/face:active"); } - memset(activecell,0,maxactive*sizeof(int)); + memset(activecell,0,((size_t)maxactive)*sizeof(int)); for (i = 0; i < ntask; i++) activecell[tasks[i].pcell] = 1; active_current = 1; } @@ -1455,7 +1455,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]; @@ -1494,7 +1494,7 @@ void FixEmitSurf::grow_task() // set all new task bytes to 0 so valgrind won't complain // if bytes between fields are uninitialized - memset(&tasks[oldmax],0,(ntaskmax-oldmax)*sizeof(Task)); + memset(&tasks[oldmax],0,((size_t)ntaskmax-oldmax)*sizeof(Task)); // allocate vectors in each new task or set to NULL // path and fracarea are allocated later to specific sizes diff --git a/src/fix_grid_check.cpp b/src/fix_grid_check.cpp index 49a7fc50e..f2660ffff 100644 --- a/src/fix_grid_check.cpp +++ b/src/fix_grid_check.cpp @@ -119,12 +119,13 @@ void FixGridCheck::end_of_step() if (outflag == ERROR) { char str[128]; 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); } nflag++; + continue; } // does particle coord match icell bounds @@ -218,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 b77c3924b..fbd53b260 100644 --- a/src/fix_halt.cpp +++ b/src/fix_halt.cpp @@ -214,7 +214,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_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 50e1aefc0..2bd8d25c2 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 849384c84..9638b6131 100644 --- a/src/grid.cpp +++ b/src/grid.cpp @@ -1473,8 +1473,8 @@ void Grid::set_inout() // create set1 and set2 lists so can swap between them int *set1,*set2; - memory->create(set1,nlocal*nface,"grid:set1"); - memory->create(set2,nlocal*nface,"grid:set2"); + memory->create(set1,(bigint)nlocal*nface,"grid:set1"); + memory->create(set2,(bigint)nlocal*nface,"grid:set2"); // initial set list = overlapped cells with corner values which are set @@ -2237,7 +2237,7 @@ 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; } @@ -2272,7 +2272,7 @@ 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; } @@ -2299,7 +2299,7 @@ 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_collate.cpp b/src/grid_collate.cpp index 2ee3489d8..affb42449 100644 --- a/src/grid_collate.cpp +++ b/src/grid_collate.cpp @@ -71,7 +71,7 @@ void Grid::collate_vector_implicit(int n, cellint *ids, int *proclist; memory->create(proclist,nsend,"grid:proclist"); double *in_rvous; - memory->create(in_rvous,2*nsend,"grid:in_rvous"); + memory->create(in_rvous,(bigint)2*nsend,"grid:in_rvous"); int m = 0; nsend = 0; @@ -146,7 +146,7 @@ void Grid::collate_array_implicit(int nrow, int ncol, cellint *ids, // zero output values for owned cells - if (nlocal) memset(&out[0][0],0,nlocal*ncol*sizeof(double)); + if (nlocal) memset(&out[0][0],0,(size_t)nlocal*ncol*sizeof(double)); // if I own grid cell, sum in values to out values directly // else nsend = # of tallies to contribute to rendezvous diff --git a/src/grid_custom.cpp b/src/grid_custom.cpp index 839a8d3e6..2f914670d 100644 --- a/src/grid_custom.cpp +++ b/src/grid_custom.cpp @@ -166,7 +166,7 @@ void Grid::reallocate_custom(int nold, int nnew) nnew,eicol[ewhich[ic]],"grid:eiarray"); if (nnew > nold) memset(&iarray[nold][0],0, - (nnew-nold)*eicol[ewhich[ic]]*sizeof(int)); + (size_t)(nnew-nold)*eicol[ewhich[ic]]*sizeof(int)); } } else { @@ -176,9 +176,9 @@ void Grid::reallocate_custom(int nold, int nnew) } else { double **darray = memory->grow(edarray[ewhich[ic]], nnew,edcol[ewhich[ic]],"grid:edarray"); - if (nnew - nold) + if (nnew > nold) memset(&darray[nold][0],0, - (nnew-nold)*edcol[ewhich[ic]]*sizeof(double)); + (size_t)(nnew-nold)*edcol[ewhich[ic]]*sizeof(double)); } } } 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 428226846..f26279e4c 100644 --- a/src/input.cpp +++ b/src/input.cpp @@ -239,8 +239,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); } } @@ -265,7 +265,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; @@ -306,8 +306,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); } @@ -1000,7 +1000,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++; @@ -1054,7 +1054,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++; @@ -1097,7 +1097,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; @@ -1122,7 +1122,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; @@ -1167,7 +1167,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); } } @@ -1245,7 +1245,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 12193f453..8131170d7 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/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/particle.cpp b/src/particle.cpp index 1ee07a278..a1e5ef623 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -707,7 +707,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); } } @@ -840,7 +840,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); } } @@ -889,7 +889,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); } } @@ -1411,12 +1411,14 @@ void Particle::read_restart_mixture(FILE *fp) int Particle::size_restart() { - int n = sizeof(int); - n = IROUNDUP(n); + bigint n = sizeof(int); + n = BIROUNDUP(n); n += nlocal * sizeof(OnePartRestart); n += nlocal * sizeof_custom(); - n = IROUNDUP(n); - return n; + n = BIROUNDUP(n); + if (n > MAXSMALLINT) + error->one(FLERR,"Per-processor particle count is too big for restart chunk"); + return static_cast(n); } /* ---------------------------------------------------------------------- diff --git a/src/react_bird.cpp b/src/react_bird.cpp index b64cb6a1a..4029ce20a 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..84858de75 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 && ecc > 0.0) 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; @@ -156,7 +157,8 @@ int ReactQK::attempt(Particle::OnePart *ip, Particle::OnePart *jp, iv = random->uniform()*(maxlev+0.99999999); evib = static_cast (iv * update->boltz*species[mspec].vibtemp[0]); - if (evib < ecc) prob = pow(1.0-evib/ecc,1.5 - r->coeff[6]); + if (evib < ecc && ecc > 0.0) 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 72fbc648f..67845cdae 100644 --- a/src/react_tce.cpp +++ b/src/react_tce.cpp @@ -103,7 +103,7 @@ int ReactTCE::attempt(Particle::OnePart *ip, Particle::OnePart *jp, if (r->coeff[1]>((-1)*r->coeff[4])) e_excess = ecc - r->coeff[1]; else e_excess = ecc + r->coeff[4]; - if (e_excess <= 0.0) continue; + if (e_excess <= 0.0 || ecc <= 0.0) continue; if (!partialEnergy) { @@ -152,9 +152,11 @@ int ReactTCE::attempt(Particle::OnePart *ip, Particle::OnePart *jp, case IONIZATION: case EXCHANGE: { - react_prob += r->coeff[2] * tgamma(z+2.5-r->coeff[5]) / MAX(1.0e-6,tgamma(z+r->coeff[3]+1.5)) * - pow(ecc-r->coeff[1],r->coeff[3]-1+r->coeff[5]) * - pow(1.0-r->coeff[1]/ecc,z+1.5-r->coeff[5]); + if (ecc > 0.0 && ecc > r->coeff[1]) { + react_prob += r->coeff[2] * tgamma(z+2.5-r->coeff[5]) / MAX(1.0e-6,tgamma(z+r->coeff[3]+1.5)) * + pow(ecc-r->coeff[1],r->coeff[3]-1+r->coeff[5]) * + pow(1.0-r->coeff[1]/ecc,z+1.5-r->coeff[5]); + } break; } @@ -172,21 +174,23 @@ int ReactTCE::attempt(Particle::OnePart *ip, Particle::OnePart *jp, int *sp2recomb = reactions[isp][jsp].sp2recomb; if (sp2recomb[recomb_species] != list[i]) continue; - react_prob += recomb_boost * recomb_density * r->coeff[2] * - tgamma(z+2.5-r->coeff[5]) / MAX(1.0e-6,tgamma(z+r->coeff[3]+1.5)) * - pow(ecc-r->coeff[1],r->coeff[3]-1+r->coeff[5]) * // extended to general recombination case with non-zero activation energy - pow(1.0-r->coeff[1]/ecc,z+1.5-r->coeff[5]); + if (ecc > 0.0 && ecc > r->coeff[1]) { + react_prob += recomb_boost * recomb_density * r->coeff[2] * + tgamma(z+2.5-r->coeff[5]) / MAX(1.0e-6,tgamma(z+r->coeff[3]+1.5)) * + pow(ecc-r->coeff[1],r->coeff[3]-1+r->coeff[5]) * // extended to general recombination case with non-zero activation energy + pow(1.0-r->coeff[1]/ecc,z+1.5-r->coeff[5]); + } 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 +240,8 @@ int ReactTCE::attempt(Particle::OnePart *ip, Particle::OnePart *jp, return list[i] + 1; } + + break; } } @@ -253,7 +259,7 @@ double ReactTCE::bird_Evib(int nmode, double Tvib, // Comutes f for Newton's search method outlined in newtonTvib() double f = -Evib; - double kb = 1.38064852e-23; + double kb = update->boltz; for (int i = 0; i < nmode; i++) { const double vti = vibtemp[i]; @@ -270,7 +276,7 @@ double ReactTCE::bird_dEvib(int nmode, double Tvib, double vibtemp[]) // Comutes df for Newton's search method double df = 0.0; - double kb = 1.38064852e-23; + double kb = update->boltz; for (int i = 0; i < nmode; i++) { const double vti = vibtemp[i]; diff --git a/src/react_tce_qk.cpp b/src/react_tce_qk.cpp index 36105ae9a..fa79b50a8 100644 --- a/src/react_tce_qk.cpp +++ b/src/react_tce_qk.cpp @@ -141,9 +141,11 @@ int ReactTCEQK::attempt_tce(Particle::OnePart *ip, Particle::OnePart *jp, case DISSOCIATION: case EXCHANGE: { - react_prob += r->coeff[2] * - pow(ecc-r->coeff[1],r->coeff[3]) * - pow(1.0-r->coeff[1]/ecc,r->coeff[5]); + if (ecc > 0.0) { + react_prob += r->coeff[2] * + pow(ecc-r->coeff[1],r->coeff[3]) * + pow(1.0-r->coeff[1]/ecc,r->coeff[5]); + } break; } @@ -229,8 +231,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 && ecc > 0.0) 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; @@ -259,7 +262,8 @@ int ReactTCEQK::attempt_qk(Particle::OnePart *ip, Particle::OnePart *jp, iv = random->uniform()*(maxlev+0.99999999); evib = static_cast (iv * update->boltz*species[mspec].vibtemp[0]); - if (evib < ecc) prob = pow(1.0-evib/ecc,1.5 - r->coeff[6]); + if (evib < ecc && ecc > 0.0) 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_isurf.cpp b/src/read_isurf.cpp index ffde42645..8707a0432 100644 --- a/src/read_isurf.cpp +++ b/src/read_isurf.cpp @@ -554,6 +554,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/surf.cpp b/src/surf.cpp index 670bb9605..070eebcf4 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; @@ -3078,11 +3078,11 @@ void Surf::grow(int old) if (domain->dimension == 2) { lines = (Surf::Line *) memory->srealloc(lines,nmax*sizeof(Line),"surf:lines"); - memset(&lines[old],0,(nmax-old)*sizeof(Line)); + memset(&lines[old],0,(size_t)(nmax-old)*sizeof(Line)); } else { tris = (Surf::Tri *) memory->srealloc(tris,nmax*sizeof(Tri),"surf:tris"); - memset(&tris[old],0,(nmax-old)*sizeof(Tri)); + memset(&tris[old],0,(size_t)(nmax-old)*sizeof(Tri)); } } diff --git a/src/surf_collide_cll.cpp b/src/surf_collide_cll.cpp index 2751c3fd3..836789bb3 100644 --- a/src/surf_collide_cll.cpp +++ b/src/surf_collide_cll.cpp @@ -136,6 +136,8 @@ SurfCollideCLL::SurfCollideCLL(SPARTA *sparta, int narg, char **arg) : SurfCollideCLL::~SurfCollideCLL() { + if (copy) return; + delete random; } diff --git a/src/surf_collide_impulsive.cpp b/src/surf_collide_impulsive.cpp index 314735245..263df8a92 100644 --- a/src/surf_collide_impulsive.cpp +++ b/src/surf_collide_impulsive.cpp @@ -141,6 +141,8 @@ SurfCollideImpulsive::SurfCollideImpulsive(SPARTA *sparta, int narg, char **arg) SurfCollideImpulsive::~SurfCollideImpulsive() { + if (copy) return; + delete random; } diff --git a/src/surf_collide_specular.cpp b/src/surf_collide_specular.cpp index 188bfdf0a..4cc32d3f1 100644 --- a/src/surf_collide_specular.cpp +++ b/src/surf_collide_specular.cpp @@ -160,7 +160,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_comm.cpp b/src/surf_comm.cpp index 89e2e1107..06aa1ece6 100644 --- a/src/surf_comm.cpp +++ b/src/surf_comm.cpp @@ -571,7 +571,7 @@ void Surf::spread_own2local_reduce(int n, int type, void *in, void *out) error->all(FLERR,"Overflow in spread_own2local_reduce"); bigint bbytes = (bigint) nlocal * n * sizeof(int); - memory->create(myvec,nlocal*n,"surf/spread:myvec"); + memory->create(myvec,(bigint)nlocal*n,"surf/spread:myvec"); memset(myvec,0,bbytes); if (n == 1) { @@ -602,7 +602,7 @@ void Surf::spread_own2local_reduce(int n, int type, void *in, void *out) error->all(FLERR,"Overflow in spread_own2local_reduce"); bigint bbytes = (bigint) nlocal * n * sizeof(double); - memory->create(myvec,nlocal*n,"surf/spread:myvec"); + memory->create(myvec,(bigint)nlocal*n,"surf/spread:myvec"); memset(myvec,0,bbytes); if (n == 1) { @@ -992,10 +992,10 @@ void Surf::spread_local2own(int n, int type, void *in, void *out) dbuf = NULL; if (type == INT) { iinput = (int *) in; - memory->create(ibuf,(n+1)*nunique,"spread/local2own:ibuf"); + memory->create(ibuf,(bigint)(n+1)*nunique,"spread/local2own:ibuf"); } else if (type == DOUBLE) { dinput = (double *) in; - memory->create(dbuf,(n+1)*nunique,"spread/local2own:dbuf"); + memory->create(dbuf,(bigint)(n+1)*nunique,"spread/local2own:dbuf"); } int isurf,index; diff --git a/src/surf_custom.cpp b/src/surf_custom.cpp index 0630b86d4..fe324d39d 100644 --- a/src/surf_custom.cpp +++ b/src/surf_custom.cpp @@ -150,24 +150,24 @@ void Surf::allocate_custom(int index) if (etype[index] == INT) { if (esize[index] == 0) { int *ivector = memory->create(eivec[ewhich[index]],n,"surf:eivec"); - if (ivector) memset(ivector,0,n*sizeof(int)); + if (ivector) memset(ivector,0,(size_t)n*sizeof(int)); eivec_local[ewhich[index]] = NULL; } else { int **iarray = memory->create(eiarray[ewhich[index]], n,eicol[ewhich[index]],"surf:eiarray"); - if (iarray) memset(&iarray[0][0],0,n*eicol[ewhich[index]]*sizeof(int)); + if (iarray) memset(&iarray[0][0],0,(size_t)n*eicol[ewhich[index]]*sizeof(int)); eiarray_local[ewhich[index]] = NULL; } } else { if (esize[index] == 0) { double *dvector = memory->create(edvec[ewhich[index]],n,"surf:edvec"); - if (dvector) memset(dvector,0,n*sizeof(double)); + if (dvector) memset(dvector,0,(size_t)n*sizeof(double)); edvec_local[ewhich[index]] = NULL; } else { double **darray = memory->create(edarray[ewhich[index]], n,edcol[ewhich[index]],"surf:eearray"); - if (darray) memset(&darray[0][0],0,n*edcol[ewhich[index]]*sizeof(double)); + if (darray) memset(&darray[0][0],0,(size_t)n*edcol[ewhich[index]]*sizeof(double)); edarray_local[ewhich[index]] = NULL; } } @@ -193,23 +193,23 @@ void Surf::reallocate_custom() if (etype[index] == INT) { if (esize[index] == 0) { int *ivector = memory->grow(eivec[ewhich[index]],nnew,"surf:eivec"); - if (nnew > nold) memset(&ivector[nold],0,(nnew-nold)*sizeof(int)); + if (nnew > nold) memset(&ivector[nold],0,(size_t)(nnew-nold)*sizeof(int)); } else { int **iarray = memory->grow(eiarray[ewhich[index]], nnew,eicol[ewhich[index]],"surf:eiarray"); if (nnew > nold) - memset(iarray[nold],0,(nnew-nold)*eicol[ewhich[index]]*sizeof(int)); + memset(iarray[nold],0,(size_t)(nnew-nold)*eicol[ewhich[index]]*sizeof(int)); } } else { if (esize[index] == 0) { double *dvector = memory->grow(edvec[ewhich[index]],nnew,"surf:edvec"); - if (nnew > nold) memset(&dvector[nold],0,(nnew-nold)*sizeof(double)); + if (nnew > nold) memset(&dvector[nold],0,(size_t)(nnew-nold)*sizeof(double)); } else { double **darray = memory->grow(edarray[ewhich[index]], nnew,edcol[ewhich[index]],"surf:eearray"); if (nnew > nold) - memset(darray[nold],0,(nnew-nold)*edcol[ewhich[index]]*sizeof(double)); + memset(darray[nold],0,(size_t)(nnew-nold)*edcol[ewhich[index]]*sizeof(double)); } } } diff --git a/src/surf_react_adsorb.cpp b/src/surf_react_adsorb.cpp index 98e6aaa23..582db61d5 100644 --- a/src/surf_react_adsorb.cpp +++ b/src/surf_react_adsorb.cpp @@ -824,9 +824,13 @@ int SurfReactAdsorb::react(Particle::OnePart *&ip, int isurf, double *norm, double dot = MathExtra::dot3(v,norm); double vmag_sq = MathExtra::lensq3(v); double E_i = 0.5 * species[ip->ispecies].mass * vmag_sq; - double cos_theta = abs(dot) / sqrt(vmag_sq); - prob_value[i] *= pow(E_i,r->energy_coeff[0]) * - pow(cos_theta,r->energy_coeff[1]); + if (vmag_sq > 0.0) { + double cos_theta = abs(dot) / sqrt(vmag_sq); + prob_value[i] *= pow(E_i,r->energy_coeff[0]) * + pow(cos_theta,r->energy_coeff[1]); + } else { + prob_value[i] = 0.0; + } } break; } diff --git a/src/timer.cpp b/src/timer.cpp index 01e5a2068..68c169e15 100644 --- a/src/timer.cpp +++ b/src/timer.cpp @@ -104,8 +104,10 @@ void Timer::init_timeout() _s_timeout = _timeout; if (_timeout < 0) _nextcheck = -1; - else + else { _nextcheck = _checkfreq; + timeout_start = MPI_Wtime(); + } } /* ---------------------------------------------------------------------- */ @@ -139,8 +141,6 @@ bool Timer::_check_timeout() // broadcast time to ensure all ranks act the same. MPI_Bcast(&walltime, 1, MPI_DOUBLE, 0, world); - printf("%g %g\n",walltime,_timeout); - if (walltime < _timeout) { _nextcheck += _checkfreq; return false; diff --git a/src/update.cpp b/src/update.cpp index ce2cb4d85..25a48457e 100644 --- a/src/update.cpp +++ b/src/update.cpp @@ -636,22 +636,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; @@ -694,13 +706,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; @@ -939,7 +957,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 @@ -1532,6 +1550,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 ecfa2dd1b..be5e9f72f 100644 --- a/src/update.h +++ b/src/update.h @@ -170,12 +170,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 e5b2dc05d..3b6fd6a75 100644 --- a/src/variable.cpp +++ b/src/variable.cpp @@ -653,10 +653,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); @@ -1283,8 +1286,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; diff --git a/src/write_grid.cpp b/src/write_grid.cpp index fb9eb6b8a..69e930c42 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 c15f8f8f0..d7a72b5fe 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 46377276c..60d5f5c7c 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 5a2e95c4f..79a30dffa 100644 --- a/src/write_surf.cpp +++ b/src/write_surf.cpp @@ -940,7 +940,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); } @@ -976,7 +976,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); }