| 3308 | } |
| 3309 | |
| 3310 | advance_status |
| 3311 | Castro::check_for_negative_density () |
| 3312 | { |
| 3313 | BL_PROFILE("check_for_negative_density()"); |
| 3314 | |
| 3315 | MultiFab& S_old = get_old_data(State_Type); |
| 3316 | MultiFab& S_new = get_new_data(State_Type); |
| 3317 | |
| 3318 | ReduceOps<ReduceOpMax, ReduceOpMax> reduce_op; |
| 3319 | ReduceData<int, int> reduce_data(reduce_op); |
| 3320 | using ReduceTuple = typename decltype(reduce_data)::Type; |
| 3321 | |
| 3322 | #ifdef AMREX_USE_OMP |
| 3323 | #pragma omp parallel |
| 3324 | #endif |
| 3325 | for (MFIter mfi(S_new, TilingIfNotGPU()); mfi.isValid(); ++mfi) { |
| 3326 | const Box& bx = mfi.tilebox(); |
| 3327 | |
| 3328 | const auto S_old_arr = S_old.array(mfi); |
| 3329 | const auto S_new_arr = S_new.array(mfi); |
| 3330 | |
| 3331 | reduce_op.eval(bx, reduce_data, |
| 3332 | [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept -> ReduceTuple |
| 3333 | { |
| 3334 | int rho_check_failed = 0; |
| 3335 | int X_check_failed = 0; |
| 3336 | |
| 3337 | const Real rho = S_new_arr(i,j,k,URHO); |
| 3338 | const Real rhoInv = 1.0_rt / rho; |
| 3339 | |
| 3340 | // Optionally, the user can ignore this if the starting |
| 3341 | // density is lower than a certain threshold. This is useful |
| 3342 | // if the minimum density occurs in material that is not |
| 3343 | // dynamically important; in that case, a density reset suffices. |
| 3344 | |
| 3345 | if (S_old_arr(i,j,k,URHO) >= retry_small_density_cutoff && rho < small_dens) { |
| 3346 | #ifndef AMREX_USE_GPU |
| 3347 | std::cout << "Invalid density = " << rho << " at index " << i << ", " << j << ", " << k << "\n"; |
| 3348 | #endif |
| 3349 | rho_check_failed = 1; |
| 3350 | } |
| 3351 | |
| 3352 | if (rho >= castro::abundance_failure_rho_cutoff) { |
| 3353 | |
| 3354 | for (int n = 0; n < NumSpec; ++n) { |
| 3355 | Real X = S_new_arr(i,j,k,UFS+n) * rhoInv; |
| 3356 | |
| 3357 | if (X < -castro::abundance_failure_tolerance || |
| 3358 | X > 1.0_rt + castro::abundance_failure_tolerance) { |
| 3359 | #ifndef AMREX_USE_GPU |
| 3360 | std::cout << "Invalid X[" << n << "] = " << X << " in zone " |
| 3361 | << i << ", " << j << ", " << k |
| 3362 | << " with density = " << rho << std::endl; |
| 3363 | #elif defined(ALLOW_GPU_PRINTF) |
| 3364 | AMREX_DEVICE_PRINTF("Invalid X[%d] = %g in zone (%d,%d,%d) with density = %g\n", |
| 3365 | n, X, i, j, k, rho); |
| 3366 | #endif |
| 3367 | X_check_failed = 1; |