Don't use pragma omp simd in Deposition (#1381)

Author: AlexanderSinn

src/fields/Fields.cpp

@@ -779,21 +779,23 @@ SetDirichletBoundaries (Array2<amrex::Real> RHS, const amrex::Box& solver_size,
         [=] AMREX_GPU_DEVICE (int i, int j) noexcept
         {
             const bool i_is_changing = (i < box_len0);
-            const bool i_lo_edge = (!i_is_changing) && (!j);
-            const bool i_hi_edge = (!i_is_changing) && j;
-            const bool j_lo_edge = i_is_changing && (!j);
-            const bool j_hi_edge = i_is_changing && j;
+            const int i_is_changing_i = static_cast<int>(i_is_changing);
+            const int i_not_changing_i = static_cast<int>(!i_is_changing);
+            const int i_lo_edge = static_cast<int>(!i_is_changing && (j == 0));
+            const int i_hi_edge = static_cast<int>(!i_is_changing && (j != 0));
+            const int j_lo_edge = static_cast<int>(i_is_changing && (j == 0));
+            const int j_hi_edge = static_cast<int>(i_is_changing && (j != 0));
 
-            const int i_idx = box_lo0 + i_hi_edge*(box_len0-1) + i_is_changing*i;
-            const int j_idx = box_lo1 + j_hi_edge*(box_len1-1) + (!i_is_changing)*(i-box_len0);
+            const int i_idx = box_lo0 + i_hi_edge*(box_len0-1) + i_is_changing_i*i;
+            const int j_idx = box_lo1 + j_hi_edge*(box_len1-1) + i_not_changing_i*(i-box_len0);
 
             const amrex::Real i_idx_offset = i_idx + (- i_lo_edge + i_hi_edge) * offset;
             const amrex::Real j_idx_offset = j_idx + (- j_lo_edge + j_hi_edge) * offset;
 
             const amrex::Real x = i_idx_offset * dx + offset0;
             const amrex::Real y = j_idx_offset * dy + offset1;
 
-            const amrex::Real dxdx = dx*dx*(!i_is_changing) + dy*dy*i_is_changing;
+            const amrex::Real dxdx = i_is_changing ? dy*dy : dx*dx;
 
             // atomic add because the corners of RHS get two values
             amrex::Gpu::Atomic::AddNoRet(&(RHS(i_idx, j_idx)),

src/laser/MultiLaser.cpp

@@ -44,9 +44,9 @@ MultiLaser::ReadParameters ()
     AMREX_ALWAYS_ASSERT(m_interp_order <= 3 && m_interp_order >= 0);
 
     bool mg_param_given = queryWithParser(pp, "MG_tolerance_rel", m_MG_tolerance_rel);
-    mg_param_given += queryWithParser(pp, "MG_tolerance_abs", m_MG_tolerance_abs);
-    mg_param_given += queryWithParser(pp, "MG_verbose", m_MG_verbose);
-    mg_param_given += queryWithParser(pp, "MG_average_rhs", m_MG_average_rhs);
+    mg_param_given = queryWithParser(pp, "MG_tolerance_abs", m_MG_tolerance_abs) || mg_param_given;
+    mg_param_given = queryWithParser(pp, "MG_verbose", m_MG_verbose) || mg_param_given;
+    mg_param_given = queryWithParser(pp, "MG_average_rhs", m_MG_average_rhs) || mg_param_given;
 
     // Raise warning if user specifies MG parameters without using the MG solver
     if (mg_param_given && (m_solver_type != "multigrid")) {

src/particles/deposition/DepositionUtil.H

@@ -240,10 +240,8 @@ SharedMemoryDeposition (int num_particles,
 
                         const int tile_id = (itile_x * ntile_y + itile_y);
 
-#ifdef AMREX_USE_OMP
-#pragma omp simd
-#endif
                         // deposit charge / current of all particles in this tile
+                        // cannot use SIMD here because particles might deposit into the same cell
                         for (int ip = a_offsets[tile_id]; ip < a_offsets[tile_id+1]; ++ip) {
                             do_deposit(a_indices[ip], ptd, field, idx_cache, idx_depos);
                         }
@@ -254,8 +252,8 @@ SharedMemoryDeposition (int num_particles,
     }
 #endif
     else {
-        // simple loop over particles, on CPU this only uses one thread
-        amrex::ParallelFor(num_particles,
+        // simple loop over particles, on CPU this only uses one thread and no SIMD
+        amrex::For(num_particles,
             [=] AMREX_GPU_DEVICE (int ip) {
                 if (is_valid(ip, ptd)) {
                     do_deposit(ip, ptd, field, idx_cache, idx_depos);