Skip to content

Commit e1ab56a

Browse files
claudestanmoore1
andcommitted
Consolidate the verified static-analysis bug fixes from both review branches
Consolidates claude/review-sparta-analysis-bugs-kyionn and claude/sparta-static-analysis-bugs-335wsw onto this branch, keeping only fixes for genuine defects that are still present here: - wrong-logic fixes: tally keyword/enum typos, compute count group check, tvib mode indexing, property/surf OOB, reaction2col wipe, rotc2 symmetry, vremax==0 NaN acceptance, react_qk/tce_qk react_prob scratch pollution, react_tce unreachable warnings + missing break, react_bird uninit tallies, adapt_grid coarsen newcell index, marching-cubes int corners, specular wrapper noslip, emit_face_file azimuth MY_2PI, tally_setup dangling *_active - KOKKOS: RNG state freed then reused in collide attempt loop, return-vs-continue in per-cell particle loops, lambda/grid missing else (also CPU), sonine dead OOB read, compute_surf tallyinfo OOB when nsurf==0, histo/weight inverted realloc + stray printf, grid-check OOB error message - crash/hang/NaN: mover frac 0/0 + clamp, stuck-particle epsilon, axi_remap axis guard, axisymmetric quadratic cancellation (Vieta) and a==0 linear case, temp rescale t_current==0, ablate Ninterface==0, subsonic emit 0-division, fix_surf_temp uninit units + stale cqw/fqw, variable lock-file fopen/fscanf, read_isurf unchecked fopen - hardening: snprintf on unbounded %s filename messages, str[32] cell-ID buffers, BIGINT_FORMAT in fix_halt - leaks on recoverable error paths (error->all throws): variable.cpp ids, suffix/list/commands leaks, write_isurf arg mutation + file leak - compute_reduce replace/subset bounds vs expanded args Fixes already present on this branch (from newer master work) were skipped, as were the audit entries both branches rejected or reverted: inert integer casts, unreachable guards, RNG seed/idiom changes, and cosmetic snprintf conversions. See BUG_VERIFICATION.md for the full consolidated table, skip rationale, and items flagged for domain sign-off (react_tce break, geometry Vieta/a==0, emit azimuth). Verified with a clean serial build. Co-authored-by: Stan Moore <stanmoore1@gmail.com> Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JAGuLGWVKDbVYbwNXjqZ6K
1 parent 7ee73d3 commit e1ab56a

62 files changed

Lines changed: 482 additions & 169 deletions

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

‎BUG_VERIFICATION.md‎

Lines changed: 127 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,127 @@
1+
# Consolidated bug list — static-analysis review branches
2+
3+
This is the consolidation of the two audit branches
4+
`claude/review-sparta-analysis-bugs-kyionn` and `claude/sparta-static-analysis-bugs-335wsw`
5+
(both audited the ~111 AI-reported bugs against `origin/master` @ `5aed836`), re-verified
6+
against **this** branch, which already carries ~182 newer master commits plus the KOKKOS
7+
porting work. Only defects that are still present here and are genuine were applied.
8+
Bug numbers (#N) refer to the original numbered report used by both audit branches.
9+
10+
The two branches share one fix commit; `kyionn` then reverted inert casts, unreachable
11+
guards, an RNG-seed change, and the `-log(1.0-x)` RNG hardening. Those reverts were all
12+
verified correct, so the consolidation started from `kyionn` and was pruned further.
13+
14+
## Applied fixes (grouped)
15+
16+
### Wrong results / wrong logic
17+
| File(s) | Bug | Defect |
18+
|--|--|--|
19+
| `compute_gas_collision_tally.cpp` | 10 | `type2` keyword mapped to `TYPE1` |
20+
| `compute_gas_reaction_tally.cpp` | 11 | `vy2/pre`→`VX2PRE`, `vz2/pre`→`VY2PRE` enum shift |
21+
| `compute_surf_reaction_tally.cpp` | 12 | `ID2POST\|\|ID2POST` should be `ID1POST\|\|ID2POST` (datatype) |
22+
| `compute_count.cpp` | 9 | group-not-found check tested `imix<0` instead of `igroup<0` |
23+
| `compute_tvib_grid.cpp` | 68 | mode output used `groupspecies[index]` instead of `groupspecies[index/maxmode]` → OOB read |
24+
| `compute_property_surf.cpp` | 6 | 3D `pack_id` looped `nsown` over `cglobal[nchoose]` → OOB read |
25+
| `compute_react_surf/boundary/isurf_grid.cpp` | 69 | `reaction2col` re-zeroed inside the strtok loop, wiping earlier matches (array is pre-zeroed at create) |
26+
| `collide_vss.cpp` | 47 | `rotc2` not assigned symmetrically (`[jsp][isp]` left stale) |
27+
| `collide_vss.cpp` | 46 | first-attempt cell with `vremax==0` → `0/0` NaN compare accepts a zero-relative-velocity collision → NaN velocities |
28+
| `react_qk.cpp`, `react_tce_qk.cpp` | 4/72 | vibrational-level rejection loop used `react_prob` as scratch, polluting the downstream `react_prob > random_prob` test; also `prob` left stale/uninitialized when `evib >= ecc` |
29+
| `react_tce.cpp` | 2 | react_prob warnings sat unreachable inside the `switch`; moved after it |
30+
| `react_tce.cpp` | 49 | missing `break` after a reaction is selected: in `computeChemRates` mode every later reaction in the list is also tallied with the same cumulative probability — **changes chem-rate output; flagged for domain sign-off** |
31+
| `react_bird.cpp` | 54 | `tally_reactions{,_all}` not NULLed in the 1-arg ctor → uninitialized delete |
32+
| `adapt_grid.cpp` | 81 | `newcell = grid->nlocal-1` after `coarsen_cell`; correct index is `nlocal` before the call |
33+
| `marching_cubes.h` | 30 | corner values `v000..v111` declared `int`, assigned/averaged as `double` → truncation in interpolation |
34+
| `surf_collide_specular.cpp` | 15 | `wrapper()` ignored `noslip_flag` and always did specular reflection |
35+
| `fix_emit_face_file.cpp` | 24 | `perform_task` azimuth drawn from `MY_PI` instead of `MY_2PI` (subsonic twin already fixed upstream) — biases tangential velocities positive |
36+
| `update.cpp` | 53 | `tally_setup()` deleted `*_active` arrays without NULLing → dangling/double delete when counts drop to 0 |
37+
38+
### KOKKOS-specific wrong results / races
39+
| File(s) | Bug | Defect |
40+
|--|--|--|
41+
| `collide_vss_kokkos.cpp` | 73 | `rand_pool.free_state()` before `continue` inside the attempt loop: the state is reused after being freed (race/correlated streams) and freed again at loop exit |
42+
| `compute_eflux_grid_kokkos.cpp`, `compute_thermal_grid_kokkos.cpp` | 44 | `return` instead of `continue` in the per-cell particle loop: one out-of-mixture particle skips the rest of the cell's particles |
43+
| `compute_lambda_grid_kokkos.cpp` + CPU `compute_lambda_grid.cpp` | 19/38/52 | missing `else` on the KNY/KNZ outputs: with a single output, writes the unallocated `array_grid` |
44+
| `compute_sonine_grid_kokkos.cpp` | 21 | dead `d_particles[icell]` read in `normalize_vcom` (indexes particles by cell id) — OOB device read |
45+
| `compute_surf_kokkos.cpp` | 39 | `h_surf2tally[iend]` evaluated before `iend > 0` — OOB read when `nsurf==0` (istart side already fixed upstream) |
46+
| `fix_ave_histo_weight_kokkos.cpp` | 20 | second `bin_particles` overload reallocs `k_match` when it is too **big** (`>` vs `<`) → OOB write when too small |
47+
| `fix_ave_histo_weight_kokkos.cpp` | 22 | stray debug `printf` |
48+
| `fix_grid_check_kokkos.cpp` | 37 | invalid-cell error message read `cells[icell].id` with the invalid `icell` → OOB host read (device-side `return` already fixed upstream) |
49+
50+
### Crashes / hangs / NaN poisoning
51+
| File(s) | Bug | Defect |
52+
|--|--|--|
53+
| `update.cpp`, `update_kokkos.cpp` | 77/78 | box-exit fraction `0/0` when `xnew==x` on the crossed face; clamped to `[0,1]` (guards only fire in degenerate states) |
54+
| `update.cpp`, `update_kokkos.cpp` | 79 | stuck-particle detection compared `minparam == 0.0`; a bounce advancing by ~1e-17 loops forever → `<= 1e-14` |
55+
| `update.h`, `update_kokkos.h` | 62 | `axi_remap` divides by `x[1]==0` when a particle sits exactly on the axis → NaN velocities |
56+
| `geometry.cpp`, `geometry_kokkos.h` | 80 | axisymmetric quadratics: catastrophic cancellation in `(-b ± sqrt)/a` (Vieta form used for the cancelling root), and `a==0` (trajectory parallel to cone slope) returned "no hit" instead of solving the linear case — **numerically verified, flagged for domain sign-off**; also `xc[1]==0` axis guard in the collision-point velocity rotation |
57+
| `fix_temp_rescale.cpp` | 70/97 | `sqrt(t_target/t_current)` with `t_current==0` (perfectly cold cell) → inf/NaN velocities; both paths guarded |
58+
| `fix_ablate_multi_inner.cpp` | 66 | `total/Ninterface` with `Ninterface==0` → NaN into ablation values (both call sites) |
59+
| `fix_emit_face.cpp`, `fix_emit_face_file.cpp`, `fix_emit_surf.cpp` | 96 | subsonic pressure correction divides by `massrho_cell*soundspeed_cell == 0` (cold/empty cell) → NaN vstream |
60+
| `fix_surf_temp.cpp` | 50 | neither `si` nor `cgs` left `prefactor`/`threshold` uninitialized → now errors out |
61+
| `fix_surf_temp.cpp` | 51 | `cqw`/`fqw` resolved once at construction; stale/dangling if the compute/fix list changes between runs → re-resolved in `init()` |
62+
| `variable.cpp` | 8/34 | universe-variable lock file: unchecked `fopen` (NULL deref) and unchecked `fscanf` (uninitialized `nextindex`) |
63+
| `read_isurf.cpp` | 31 | unchecked binary `fopen` in `read_corners_parallel` → NULL `fseek/fread` segfault |
64+
| `fix_grid_check.cpp` | 27 | inside-surfs error printed `icell` as the third coordinate; now prints `x[2]` |
65+
66+
### Buffer overflows (unbounded `%s` / undersized buffers)
67+
| File(s) | Bug | Defect |
68+
|--|--|--|
69+
| `dump.cpp`, `dump_grid.cpp`, `grid_id.cpp` | 110 | `char str[32]` too small for deep multi-level cell-ID strings and wide user formats; `id_num2str` now bounds its writes (128) |
70+
| `input.cpp`, `move_surf.cpp`, `particle.cpp`, `write_grid.cpp`, `write_isurf.cpp`, `write_restart.cpp`, `write_surf.cpp`, `dump_movie.cpp` | 28/32/89/109 | `sprintf(str128, "... %s", filename)` with unbounded filenames → `snprintf` (only `%s` sites; pure-numeric formats left alone) |
71+
| `surf.cpp` | 33 | suffix style name `sprintf` into `estyle[256]` from unbounded `arg[1]` → `snprintf` (×2) |
72+
| `utils.cpp` | 7 | `missing_cmd_args` formats unbounded command name into `msg[128]` → `snprintf` |
73+
| `input.cpp` | 29 | "Unknown command" built in a leaked heap buffer → stack buffer with truncation |
74+
| `fix_halt.cpp` | 26 | `%ld` for a `bigint` → `BIGINT_FORMAT` |
75+
76+
### Memory leaks on recoverable error paths (`error->all` throws; library callers recover)
77+
| File(s) | Bug | Defect |
78+
|--|--|--|
79+
| `variable.cpp` | 108 | `id` leaked at every `error->all` between allocation and `delete [] id` (compute/fix/surf-collide/surf-react/custom/`v_` blocks) |
80+
| `compute_reduce.cpp`, `compute_lambda_grid.cpp`, `fix_ave_{grid,histo,surf,time}.cpp` | 3(part)/103/99 | `suffix` leaked on the malformed-bracket error path |
81+
| `grid.cpp` | 112 | `list` leaked when a group ID doesn't exist (3 sites) |
82+
| `input.cpp` | 111 | `commands[]` leaked on illegal `if` command |
83+
| `write_isurf.cpp` | 35 | command mutated `arg[4]` in place (truncates the filename for any repeat invocation) and leaked `file` |
84+
85+
### Argument-parsing bounds
86+
| File(s) | Bug | Defect |
87+
|--|--|--|
88+
| `compute_reduce.cpp` | 3 | `replace`/`subset` bounds checked against `narg` instead of `nargnew` (expanded args) |
89+
90+
## Already fixed on this branch (skipped — no change needed)
91+
Bugs 1 (`comm.cpp` double alloc), 5 (`grid_custom` shrink memset), 13 (CPU grid-check invalid cell),
92+
14 (`timer.cpp`), 16/17 (CLL/impulsive dtor double free — RNG moved to base class),
93+
18 (device-side grid-check `return`), 23/67/98 (`fix_ablate` `idsource` shadow),
94+
25 (`fflag/fuser` re-init leak), 36 (`custom.cpp` FILECOARSE leak),
95+
41 (`sr_map` nglob), 42 (DualView host alloc), 43 (`create_particles_kokkos` `inew`),
96+
48/74 (hardcoded kb, CPU+KOKKOS), 55 (`size_restart` → `size_restart_big`),
97+
56 (`grid_collate` memset), 57 (`grid_custom` array memsets), 59-array (`surf_custom` memsets),
98+
71 (`memory_usage` `=` vs `+=`), 88 (`remove_old` smalloc casts), 95 (`for(m…;i++)` typo),
99+
20-first-site, 39-istart-side, and one of the two `MY_PI` azimuth sites.
100+
101+
## Rejected / intentionally not applied
102+
- **40** (`fft2d_kokkos` flag flip): not a bug — both original fix branches regressed it; rejected by both audits.
103+
- **84, 93, 102**: not bugs (already-64-bit products, bounded numeric buffer).
104+
- **59-vector, 60-own2local, 61, 82, 83** and the emit-file `memset`/`create` casts (58, 60-local2own, 65, part of 56): inert — `int*size_t` already promotes, and `Memory::create(..., int n, ...)` truncates a `bigint` count argument right back, so the casts change nothing.
105+
- **75** (KOKKOS RNG pool seed change): speculative, changes reproducibility.
106+
- **85/86/87** (`-log(1.0-drand())` hardening): guards a ~2⁻⁵³-rare `-log(0)` at the cost of shifting every RNG baseline; reverted by the second audit per maintainer preference.
107+
- **100** (`vrm_max>0` in compute dt/grid), **105** (`volume>0` in `attempt_collision`): unreachable — upstream checks already exclude the zero cases.
108+
- **101/38-guards** (volume/size zero-guards in eflux/pflux/thermal/lambda outputs): cell edge lengths are always positive, and `volume==0` cells cannot legally hold particles (fix grid/check errors on them) — output-side guards would only mask an already-invalid state.
109+
- **76** (`vr2>0` scatter guards): unreachable once the bug-46 `vremax==0` guard exists — a `vr2==0` pair can no longer be accepted for collision; skipping also keeps CPU/KOKKOS scatter code identical.
110+
- **106-remnant** (`ecc>0` guards in TCE/QK attempt): the `e_excess <= 0` checks upstream already imply `ecc > coeff[1]`, making the pow bases positive.
111+
- **107** (`vmag_sq>0` in surf_react_adsorb): a zero-velocity particle cannot reach a surface.
112+
- **45** (remap3d_kokkos malloc-failure cleanup): leaks only on an OOM path that aborts the run.
113+
- **cosmetic `sprintf`→`snprintf`** on bounded pure-numeric formats.
114+
115+
## Flagged for domain sign-off (applied, but they change physics output)
116+
- **49** `react_tce.cpp` missing `break`: chem-rate tallies (`computeChemRates`) will decrease to one reaction per collision.
117+
- **80** `geometry.cpp`/`geometry_kokkos.h` Vieta rewrite + `a==0` linear branch: intersection times change at the roundoff level, and trajectories exactly parallel to a cone's slope now hit instead of pass through.
118+
- **24** `fix_emit_face_file.cpp` `MY_2PI`: emitted tangential velocity distribution changes (was biased to one half-plane).
119+
- **46/79** and the frac clamps alter behavior only in previously-NaN/hung states.
120+
121+
## Verification
122+
- Every applied fix was checked against the current code of this branch (many of the
123+
original 111 were already fixed by the ~182 master commits merged since the audit
124+
baseline; those were skipped rather than re-applied).
125+
- `make serial` compiles and links `spa_serial` cleanly with all changes.
126+
- KOKKOS edits mirror their verified CPU twins or are minimal targeted patches
127+
(reviewed by hand; not part of the serial build).

‎src/KOKKOS/collide_vss_kokkos.cpp‎

Lines changed: 0 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -982,7 +982,6 @@ void CollideVSSKokkos::operator()(TagCollideCollisionsOne< NEARCP, GASTALLY, ATO
982982
else
983983
reduce.nreact_one++;
984984
} else {
985-
rand_pool.free_state(rand_gen);
986985
continue;
987986
}
988987

@@ -2590,7 +2589,6 @@ void CollideVSSKokkos::operator()(TagCollideCollisionsOneAmbipolar< GASTALLY, AT
25902589
else
25912590
reduce.nreact_one++;
25922591
} else {
2593-
rand_pool.free_state(rand_gen);
25942592
continue;
25952593
}
25962594

‎src/KOKKOS/compute_eflux_grid_kokkos.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -241,7 +241,7 @@ void ComputeEFluxGridKokkos::operator()(TagComputeEFluxGrid_compute_per_grid, co
241241

242242
const int ispecies = d_particles[i].ispecies;
243243
const int igroup = d_s2g(imix,ispecies);
244-
if (igroup < 0) return;
244+
if (igroup < 0) continue;
245245

246246
const double mass = d_species[ispecies].mass;
247247
double *v = d_particles[i].v;

‎src/KOKKOS/compute_lambda_grid_kokkos.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -326,12 +326,12 @@ void ComputeLambdaGridKokkos::compute_per_grid_kokkos()
326326

327327
if (l_knyflag) {
328328
if (l_noutputs == 1) l_vector_grid[i] = lambda / sizey;
329-
l_array_grid(i,l_output_order[KNY]) = lambda / sizey;
329+
else l_array_grid(i,l_output_order[KNY]) = lambda / sizey;
330330
}
331331

332332
if (l_knzflag) {
333333
if (l_noutputs == 1) l_vector_grid[i] = lambda / sizez;
334-
l_array_grid(i,l_output_order[KNZ]) = lambda / sizez;
334+
else l_array_grid(i,l_output_order[KNZ]) = lambda / sizez;
335335
}
336336
});
337337
}

‎src/KOKKOS/compute_sonine_grid_kokkos.cpp‎

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -225,9 +225,6 @@ void ComputeSonineGridKokkos::operator()(TagComputeSonineGrid_compute_vcom, cons
225225

226226
KOKKOS_INLINE_FUNCTION
227227
void ComputeSonineGridKokkos::operator()(TagComputeSonineGrid_normalize_vcom, const int &icell) const {
228-
const int ispecies = d_particles[icell].ispecies;
229-
const int igroup = d_s2g(imix,ispecies);
230-
231228
double norm;
232229
for (int j=0; j<ngroup; j++) {
233230
norm = d_vcom(icell,j,3);

‎src/KOKKOS/compute_surf_kokkos.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -226,7 +226,7 @@ int ComputeSurfKokkos::tallyinfo(surfint *&ptr)
226226

227227
while (1) {
228228
while (istart < nsurf && h_surf2tally[istart] != -1) istart++;
229-
while (h_surf2tally[iend] == -1 && iend > 0) iend--;
229+
while (iend > 0 && h_surf2tally[iend] == -1) iend--;
230230
if (istart >= iend) {
231231
ntally = istart;
232232
break;

‎src/KOKKOS/compute_thermal_grid_kokkos.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -166,7 +166,7 @@ void ComputeThermalGridKokkos::operator()(TagComputeThermalGrid_compute_per_grid
166166

167167
const int ispecies = d_particles[i].ispecies;
168168
const int igroup = d_s2g(imix,ispecies);
169-
if (igroup < 0) return;
169+
if (igroup < 0) continue;
170170

171171
const int icell = d_particles[i].icell;
172172

‎src/KOKKOS/fix_ave_histo_weight_kokkos.cpp‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -293,7 +293,6 @@ void FixAveHistoWeightKokkos::calculate_weights()
293293
// explicit per-particle attributes
294294
// NOTE: need to allocate local storage
295295
} else {
296-
printf("%d, %d\n", which[i] == VARIABLE, kind == PERGRID);
297296
error->all(FLERR,"Fix ave/histo/weight/kokkos option not yet supported");
298297
}
299298
}
@@ -401,7 +400,7 @@ void FixAveHistoWeightKokkos::bin_particles(
401400

402401
KokkosBase* regionKKBase = dynamic_cast<KokkosBase*>(region);
403402

404-
if (k_match.extent(0) > nmax)
403+
if (k_match.extent(0) < nmax)
405404
MemKK::realloc_kokkos(k_match,"fix_ave_histo_weight:match",nmax);
406405

407406
regionKKBase->match_all_kokkos(k_match);

‎src/KOKKOS/fix_grid_check_kokkos.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -191,9 +191,9 @@ void FixGridCheckKokkos::end_of_step()
191191
auto icell = particles[i].icell;
192192
if (h_particle_problems(i) & IS_IN_INVALID_CELL) {
193193
sprintf(str,
194-
"Particle %d,%d on proc %d is in invalid cell " CELLINT_FORMAT
194+
"Particle %d,%d on proc %d is in invalid cell index %d"
195195
" on timestep " BIGINT_FORMAT,
196-
i,particles[i].id,comm->me,cells[icell].id,update->ntimestep);
196+
i,particles[i].id,comm->me,icell,update->ntimestep);
197197
error->one(FLERR,str);
198198
}
199199
if (h_particle_problems(i) & IS_OUTSIDE_CELL) {

‎src/KOKKOS/geometry_kokkos.h‎

Lines changed: 44 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -156,10 +156,20 @@ bool axi_horizontal_line(double tdelta, double *x, double *v,
156156
double arg = yhoriz*yhoriz*a - v[2]*v[2]*x[1]*x[1];
157157
if (arg < 0.0) return false;
158158
double sarg = sqrt(arg);
159+
double c = x[1]*x[1] - yhoriz*yhoriz;
159160

160161
nc = 2;
161-
double tone = (b - sarg) / a;
162-
double ttwo = (b + sarg) / a;
162+
double tone, ttwo;
163+
if (b > 0.0) {
164+
ttwo = (b + sarg) / a;
165+
tone = c / (b + sarg);
166+
} else if (b < 0.0) {
167+
tone = (b - sarg) / a;
168+
ttwo = c / (b - sarg);
169+
} else {
170+
tone = -sarg / a;
171+
ttwo = sarg / a;
172+
}
163173
t1 = MIN(tone,ttwo);
164174
t2 = MAX(tone,ttwo);
165175

@@ -265,20 +275,33 @@ bool axi_line_intersect(double tdelta, double *x, double *v,
265275
double dconst = x21*v1[1] - y21*v1[0];
266276

267277
double a = x21sq*(v[1]*v[1] + v[2]*v[2]) - y21sq*v[0]*v[0];
268-
if (a == 0.0) return false;
269278
double b = x21sq*x[1]*v[1] - y21sq*x[0]*v[0] - y21*v[0]*dconst;
270279
double c = x21sq*x[1]*x[1] - y21sq*x[0]*x[0] -
271280
2.0*y21*x[0]*dconst - dconst*dconst;
272281

273-
double arg = b*b - a*c;
274-
if (arg < 0.0) return false;
275-
double sarg = sqrt(arg);
276-
277-
nc = 2;
278-
double tone = (-b - sarg) / a;
279-
double ttwo = (-b + sarg) / a;
280-
t1 = MIN(tone,ttwo);
281-
t2 = MAX(tone,ttwo);
282+
if (a == 0.0) {
283+
if (b == 0.0) return false;
284+
nc = 1;
285+
t1 = t2 = -0.5 * c / b;
286+
} else {
287+
double arg = b*b - a*c;
288+
if (arg < 0.0) return false;
289+
double sarg = sqrt(arg);
290+
nc = 2;
291+
double tone, ttwo;
292+
if (b > 0.0) {
293+
tone = (-b - sarg) / a;
294+
ttwo = c / (-b - sarg);
295+
} else if (b < 0.0) {
296+
tone = c / (-b + sarg);
297+
ttwo = (-b + sarg) / a;
298+
} else {
299+
tone = -sarg / a;
300+
ttwo = sarg / a;
301+
}
302+
t1 = MIN(tone,ttwo);
303+
t2 = MAX(tone,ttwo);
304+
}
282305
}
283306

284307
// if selfflag, particle starts on surf line segment
@@ -341,11 +364,16 @@ bool axi_line_intersect(double tdelta, double *x, double *v,
341364
if (v1[1] == v2[1]) xc[1] = v1[1];
342365
xc[2] = 0.0;
343366

344-
double rn = ynew / xc[1];
345-
double wn = znew / xc[1];
346367
vc[0] = v[0];
347-
vc[1] = v[1]*rn + v[2]*wn;
348-
vc[2] = -v[1]*wn + v[2]*rn;
368+
if (xc[1] > 0.0) {
369+
double rn = ynew / xc[1];
370+
double wn = znew / xc[1];
371+
vc[1] = v[1]*rn + v[2]*wn;
372+
vc[2] = -v[1]*wn + v[2]*rn;
373+
} else {
374+
vc[1] = v[1];
375+
vc[2] = v[2];
376+
}
349377

350378
// test that xc is within line segment bounds
351379
// y-test for vertical line, else x-test

0 commit comments

Comments
 (0)