Consolidate the verified static-analysis bug fixes from both review branches

src/KOKKOS/collide_vss_kokkos.cpp

@@ -982,7 +982,6 @@ void CollideVSSKokkos::operator()(TagCollideCollisionsOne< NEARCP, GASTALLY, ATO
       else
         reduce.nreact_one++;
     } else {
-      rand_pool.free_state(rand_gen);
       continue;
     }
 
@@ -2590,7 +2589,6 @@ void CollideVSSKokkos::operator()(TagCollideCollisionsOneAmbipolar< GASTALLY, AT
       else
         reduce.nreact_one++;
     } else {
-      rand_pool.free_state(rand_gen);
       continue;
     }
 

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;

src/KOKKOS/compute_lambda_grid_kokkos.cpp

@@ -326,12 +326,12 @@ void ComputeLambdaGridKokkos::compute_per_grid_kokkos()
 
       if (l_knyflag) {
         if (l_noutputs == 1) l_vector_grid[i] = lambda / sizey;
-        l_array_grid(i,l_output_order[KNY]) = lambda / sizey;
+        else l_array_grid(i,l_output_order[KNY]) = lambda / sizey;
       }
 
       if (l_knzflag) {
         if (l_noutputs == 1) l_vector_grid[i] = lambda / sizez;
-        l_array_grid(i,l_output_order[KNZ]) = lambda / sizez;
+        else l_array_grid(i,l_output_order[KNZ]) = lambda / sizez;
       }
     });
   }

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<ngroup; j++) {
     norm = d_vcom(icell,j,3);

src/KOKKOS/compute_surf_kokkos.cpp

@@ -226,7 +226,7 @@ int ComputeSurfKokkos::tallyinfo(surfint *&ptr)
 
   while (1) {
     while (istart < nsurf && h_surf2tally[istart] != -1) istart++;
-    while (h_surf2tally[iend] == -1 && iend > 0) iend--;
+    while (iend > 0 && h_surf2tally[iend] == -1) iend--;
     if (istart >= iend) {
       ntally = istart;
       break;

src/KOKKOS/compute_surf_kokkos.h

@@ -409,8 +409,9 @@ void surf_tally_kk(double /*dtremain*/, int isurf, int icell, int reaction,
         double r_coeff = 0.0;
         if (sr_type == 1)
           r_coeff = sr_kk_prob_copy[m].obj.d_coeffs(reaction-1,1);
-        a_array_surf_tally(itally,k++) += weight * r_coeff * fluxscale;
+        a_array_surf_tally(itally,k) += weight * r_coeff * fluxscale;
       }
+      k++;
       break;
     case ETOT:
       if (iorig) vsqpre = origmass * MathExtraKokkos::lensq3(vorig);

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;
 

src/KOKKOS/fft2d_kokkos.cpp

@@ -714,23 +714,23 @@ void FFT2dKokkos<DeviceType>::fft_2d_1d_only_kokkos(typename FFT_AT::t_FFT_DATA_
   // data is just an array of 0.0
 
 #if defined(FFT_KOKKOS_MKL_GPU)
-  if (flag == -1) {
+  if (flag == 1) {
     oneapi::mkl::dft::compute_forward(*(plan->desc_fast), (FFT_SCALAR*)d_data.data());
     oneapi::mkl::dft::compute_forward(*(plan->desc_slow), (FFT_SCALAR*)d_data.data());
   } else {
     oneapi::mkl::dft::compute_backward(*(plan->desc_fast), (FFT_SCALAR*)d_data.data());
     oneapi::mkl::dft::compute_backward(*(plan->desc_slow), (FFT_SCALAR*)d_data.data());
   }
 #elif defined(FFT_KOKKOS_MKL)
-  if (flag == -1) {
+  if (flag == 1) {
     DftiComputeForward(plan->handle_fast,d_data.data());
     DftiComputeForward(plan->handle_slow,d_data.data());
   } else {
     DftiComputeBackward(plan->handle_fast,d_data.data());
     DftiComputeBackward(plan->handle_slow,d_data.data());
   }
 #elif defined(FFT_KOKKOS_FFTW3) || defined(FFT_KOKKOS_NVPL)
-  if (flag == -1) {
+  if (flag == 1) {
     FFTW_API(execute_dft)(plan->plan_fast_forward,(FFT_KOKKOS_DATA*)d_data.data(),(FFT_KOKKOS_DATA*)d_data.data());
     FFTW_API(execute_dft)(plan->plan_slow_forward,(FFT_KOKKOS_DATA*)d_data.data(),(FFT_KOKKOS_DATA*)d_data.data());
   } else {
@@ -747,7 +747,7 @@ void FFT2dKokkos<DeviceType>::fft_2d_1d_only_kokkos(typename FFT_AT::t_FFT_DATA_
   kiss_fft_functor<DeviceType> f;
     typename FFT_AT::t_FFT_DATA_1d d_tmp =
      typename FFT_AT::t_FFT_DATA_1d(Kokkos::view_alloc("fft_2d:tmp",Kokkos::WithoutInitializing),d_data.extent(0));
-  if (flag == -1) {
+  if (flag == 1) {
     f = kiss_fft_functor<DeviceType>(d_data,d_tmp,plan->cfg_fast_forward,length1);
     Kokkos::parallel_for(total1/length1,f);
 
@@ -765,7 +765,7 @@ void FFT2dKokkos<DeviceType>::fft_2d_1d_only_kokkos(typename FFT_AT::t_FFT_DATA_
   // scaling if required
   // limit num to size of data
 
-  if (flag == 1 && plan->scaled) {
+  if (flag == -1 && plan->scaled) {
     FFT_SCALAR norm = plan->norm;
     int num = MIN(plan->normnum,nsize);
 

src/KOKKOS/fft3d_kokkos.cpp

@@ -915,7 +915,7 @@ void FFT3dKokkos<DeviceType>::fft_3d_1d_only_kokkos(typename FFT_AT::t_FFT_DATA_
   // data is just an array of 0.0
 
 #if defined(FFT_KOKKOS_MKL_GPU)
-  if (flag == -1) {
+  if (flag == 1) {
     oneapi::mkl::dft::compute_forward(*(plan->desc_fast), (FFT_SCALAR*)d_data.data());
     oneapi::mkl::dft::compute_forward(*(plan->desc_mid), (FFT_SCALAR*)d_data.data());
     oneapi::mkl::dft::compute_forward(*(plan->desc_slow), (FFT_SCALAR*)d_data.data());
@@ -925,7 +925,7 @@ void FFT3dKokkos<DeviceType>::fft_3d_1d_only_kokkos(typename FFT_AT::t_FFT_DATA_
     oneapi::mkl::dft::compute_backward(*(plan->desc_slow), (FFT_SCALAR*)d_data.data());
   }
 #elif defined(FFT_KOKKOS_MKL)
-  if (flag == -1) {
+  if (flag == 1) {
     DftiComputeForward(plan->handle_fast,d_data.data());
     DftiComputeForward(plan->handle_mid,d_data.data());
     DftiComputeForward(plan->handle_slow,d_data.data());
@@ -935,7 +935,7 @@ void FFT3dKokkos<DeviceType>::fft_3d_1d_only_kokkos(typename FFT_AT::t_FFT_DATA_
     DftiComputeBackward(plan->handle_slow,d_data.data());
   }
 #elif defined(FFT_KOKKOS_FFTW3) || defined(FFT_KOKKOS_NVPL)
-  if (flag == -1) {
+  if (flag == 1) {
     FFTW_API(execute_dft)(plan->plan_fast_forward,(FFT_KOKKOS_DATA*)d_data.data(),(FFT_KOKKOS_DATA*)d_data.data());
     FFTW_API(execute_dft)(plan->plan_mid_forward,(FFT_KOKKOS_DATA*)d_data.data(),(FFT_KOKKOS_DATA*)d_data.data());
     FFTW_API(execute_dft)(plan->plan_slow_forward,(FFT_KOKKOS_DATA*)d_data.data(),(FFT_KOKKOS_DATA*)d_data.data());
@@ -956,7 +956,7 @@ void FFT3dKokkos<DeviceType>::fft_3d_1d_only_kokkos(typename FFT_AT::t_FFT_DATA_
   kiss_fft_functor<DeviceType> f;
     typename FFT_AT::t_FFT_DATA_1d d_tmp =
      typename FFT_AT::t_FFT_DATA_1d(Kokkos::view_alloc("fft_3d:tmp",Kokkos::WithoutInitializing),d_data.extent(0));
-  if (flag == -1) {
+  if (flag == 1) {
     f = kiss_fft_functor<DeviceType>(d_data,d_tmp,plan->cfg_fast_forward,length1);
     Kokkos::parallel_for(total1/length1,f);
 
@@ -980,7 +980,7 @@ void FFT3dKokkos<DeviceType>::fft_3d_1d_only_kokkos(typename FFT_AT::t_FFT_DATA_
   // scaling if required
   // limit num to size of data
 
-  if (flag == 1 && plan->scaled) {
+  if (flag == -1 && plan->scaled) {
     FFT_SCALAR norm = plan->norm;
     int num = MIN(plan->normnum,nsize);
 

src/KOKKOS/fix_ave_histo_kokkos.cpp

@@ -505,19 +505,20 @@ void FixAveHistoKokkos::bin_particles(
   int n = particle->nlocal;
   int nmax = particle->maxlocal;
 
-  Region *region;
-  if (regionflag) region = domain->regions[iregion];
+  if (regionflag) {
+    Region *region = domain->regions[iregion];
 
-  if (!region->kokkos_flag)
-    error->all(FLERR,"KOKKOS package does not (yet) support chosen region style");
+    if (!region->kokkos_flag)
+      error->all(FLERR,"KOKKOS package does not (yet) support chosen region style");
 
-  KokkosBase* regionKKBase = dynamic_cast<KokkosBase*>(region);
+    KokkosBase* regionKKBase = dynamic_cast<KokkosBase*>(region);
 
-  if (k_match.extent(0) > nmax)
-    MemKK::realloc_kokkos(k_match,"fix_ave_histo_weight:match",nmax);
+    if (k_match.extent(0) < nmax)
+      MemKK::realloc_kokkos(k_match,"fix_ave_histo_weight:match",nmax);
 
-  regionKKBase->match_all_kokkos(k_match);
-  d_match = k_match.view_device();
+    regionKKBase->match_all_kokkos(k_match);
+    d_match = k_match.view_device();
+  }
 
   if (attribute == X) {
 
@@ -569,19 +570,20 @@ void FixAveHistoKokkos::bin_particles(
 
   d_values = mirror_view_from_raw_host_array<double,DeviceType>(values, n, stride);
 
-  Region *region;
-  if (regionflag) region = domain->regions[iregion];
+  if (regionflag) {
+    Region *region = domain->regions[iregion];
 
-  if (!region->kokkos_flag)
-    error->all(FLERR,"KOKKOS package does not (yet) support chosen region style");
+    if (!region->kokkos_flag)
+      error->all(FLERR,"KOKKOS package does not (yet) support chosen region style");
 
-  KokkosBase* regionKKBase = dynamic_cast<KokkosBase*>(region);
+    KokkosBase* regionKKBase = dynamic_cast<KokkosBase*>(region);
 
-  if (k_match.extent(0) < nmax)
-    MemKK::realloc_kokkos(k_match,"fix_ave_histo_weight:match",nmax);
+    if (k_match.extent(0) < nmax)
+      MemKK::realloc_kokkos(k_match,"fix_ave_histo_weight:match",nmax);
 
-  regionKKBase->match_all_kokkos(k_match);
-  d_match = k_match.view_device();
+    regionKKBase->match_all_kokkos(k_match);
+    d_match = k_match.view_device();
+  }
 
   if (regionflag && mixflag) {
     auto policy = RangePolicy<TagFixAveHisto_BinParticles1,DeviceType>(0, n);

src/KOKKOS/fix_ave_histo_weight_kokkos.cpp

@@ -293,7 +293,6 @@ void FixAveHistoWeightKokkos::calculate_weights()
   // explicit per-particle attributes
   // NOTE: need to allocate local storage
   } else {
-    printf("%d, %d\n", which[i] == VARIABLE, kind == PERGRID);
     error->all(FLERR,"Fix ave/histo/weight/kokkos option not yet supported");
   }
 }
@@ -329,19 +328,20 @@ void FixAveHistoWeightKokkos::bin_particles(
   int n = particle->nlocal;
   int nmax = particle->maxlocal;
 
-  Region *region;
-  if (regionflag) region = domain->regions[iregion];
+  if (regionflag) {
+    Region *region = domain->regions[iregion];
 
-  if (!region->kokkos_flag)
-    error->all(FLERR,"KOKKOS package does not (yet) support chosen region style");
+    if (!region->kokkos_flag)
+      error->all(FLERR,"KOKKOS package does not (yet) support chosen region style");
 
-  KokkosBase* regionKKBase = dynamic_cast<KokkosBase*>(region);
+    KokkosBase* regionKKBase = dynamic_cast<KokkosBase*>(region);
 
-  if (k_match.extent(0) < nmax)
-    MemKK::realloc_kokkos(k_match,"fix_ave_histo_weight:match",nmax);
+    if (k_match.extent(0) < nmax)
+      MemKK::realloc_kokkos(k_match,"fix_ave_histo_weight:match",nmax);
 
-  regionKKBase->match_all_kokkos(k_match);
-  d_match = k_match.view_device();
+    regionKKBase->match_all_kokkos(k_match);
+    d_match = k_match.view_device();
+  }
 
   if (attribute == X) {
 
@@ -393,19 +393,20 @@ void FixAveHistoWeightKokkos::bin_particles(
 
   d_values = mirror_view_from_raw_host_array<double,DeviceType>(values, n, stride);
 
-  Region *region;
-  if (regionflag) region = domain->regions[iregion];
+  if (regionflag) {
+    Region *region = domain->regions[iregion];
 
-  if (!region->kokkos_flag)
-    error->all(FLERR,"KOKKOS package does not (yet) support chosen region style");
+    if (!region->kokkos_flag)
+      error->all(FLERR,"KOKKOS package does not (yet) support chosen region style");
 
-  KokkosBase* regionKKBase = dynamic_cast<KokkosBase*>(region);
+    KokkosBase* regionKKBase = dynamic_cast<KokkosBase*>(region);
 
-  if (k_match.extent(0) > nmax)
-    MemKK::realloc_kokkos(k_match,"fix_ave_histo_weight:match",nmax);
+    if (k_match.extent(0) < nmax)
+      MemKK::realloc_kokkos(k_match,"fix_ave_histo_weight:match",nmax);
 
-  regionKKBase->match_all_kokkos(k_match);
-  d_match = k_match.view_device();
+    regionKKBase->match_all_kokkos(k_match);
+    d_match = k_match.view_device();
+  }
 
   if (regionflag && mixflag) {
     //auto policy = RangePolicy<TagFixAveHistoWeight_BinParticles1,DeviceType>(0, n);

src/KOKKOS/fix_grid_check_kokkos.cpp

@@ -191,9 +191,9 @@ void FixGridCheckKokkos::end_of_step()
       auto icell = particles[i].icell;
       if (h_particle_problems(i) & IS_IN_INVALID_CELL) {
         sprintf(str,
-                "Particle %d,%d on proc %d is in invalid cell " CELLINT_FORMAT
+                "Particle %d,%d on proc %d is in invalid cell index %d"
                 " on timestep " BIGINT_FORMAT,
-                i,particles[i].id,comm->me,cells[icell].id,update->ntimestep);
+                i,particles[i].id,comm->me,icell,update->ntimestep);
         error->one(FLERR,str);
       }
       if (h_particle_problems(i) & IS_OUTSIDE_CELL) {

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