diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp index f380bf5cd..65684f064 100644 --- a/src/impl/particles_impl.ipp +++ b/src/impl/particles_impl.ipp @@ -48,6 +48,7 @@ namespace libcloudphxx n_part_to_init; // number of SDs to be initialized by source detail::rng rng; detail::config config; + as_t::as_t adve_scheme; // actual advection scheme used, might be different from opts_init.adve_scheme if courant>halo // pointer to collision kernel kernel_base *p_kernel; @@ -276,6 +277,7 @@ namespace libcloudphxx n_dims == 2 ? halo_size * (opts_init.nz + 1): // 2D halo_size * (opts_init.nz + 1) * opts_init.ny // 3D ), + adve_scheme(opts_init.adve_scheme), pure_const_multi (((opts_init.sd_conc) == 0) && (opts_init.sd_const_multi > 0 || opts_init.dry_sizes.size() > 0)) // coal prob can be greater than one only in sd_conc simulations { // note: there could be less tmp data spaces if _cell vectors diff --git a/src/impl/particles_impl_adve.ipp b/src/impl/particles_impl_adve.ipp index 80f906c21..2bab648e8 100644 --- a/src/impl/particles_impl_adve.ipp +++ b/src/impl/particles_impl_adve.ipp @@ -164,12 +164,12 @@ namespace libcloudphxx { if(n_dims==0) return; - if(opts_init.adve_scheme == as_t::euler) + if(adve_scheme == as_t::euler) { adve_calc >(true, halo_x); return; } - else if(opts_init.adve_scheme == as_t::implicit) + else if(adve_scheme == as_t::implicit) { adve_calc >(true, halo_x); return; diff --git a/src/impl/particles_impl_hskpng_ijk.ipp b/src/impl/particles_impl_hskpng_ijk.ipp index 8647c7bc8..a06b67a3c 100644 --- a/src/impl/particles_impl_hskpng_ijk.ipp +++ b/src/impl/particles_impl_hskpng_ijk.ipp @@ -38,11 +38,11 @@ namespace libcloudphxx thrust_device::vector &vi, const real_t &vd ) { - thrust::transform( - vx.begin(), vx.end(), // input - vi.begin(), // output - detail::divide_by_constant_and_cast(vd) // has to be done on doubles to avoid i==nx due to low precision of nvcc math - ); + thrust::transform( + vx.begin(), vx.end(), // input + vi.begin(), // output + detail::divide_by_constant_and_cast(vd) // has to be done on doubles to avoid i==nx due to low precision of nvcc math + ); } } helper; diff --git a/src/particles_init.ipp b/src/particles_init.ipp index 2398032ac..be4cf5bac 100644 --- a/src/particles_init.ipp +++ b/src/particles_init.ipp @@ -60,10 +60,6 @@ namespace libcloudphxx if (!courant_y.is_null()) pimpl->sync(courant_y, pimpl->courant_y); if (!courant_z.is_null()) pimpl->sync(courant_z, pimpl->courant_z); - // check if courants arent greater than 2 since it would break the predictor-corrector (halo of size 2 in the x direction) - assert(pimpl->opts_init.adve_scheme != as_t::pred_corr || (courant_x.is_null() || ((*(thrust::min_element(pimpl->courant_x.begin(), pimpl->courant_x.end()))) >= real_t(-2.) )) ); - assert(pimpl->opts_init.adve_scheme != as_t::pred_corr || (courant_x.is_null() || ((*(thrust::max_element(pimpl->courant_x.begin(), pimpl->courant_x.end()))) <= real_t( 2.) )) ); - if (pimpl->opts_init.chem_switch) for (int i = 0; i < chem_gas_n; ++i) pimpl->sync( diff --git a/src/particles_step.ipp b/src/particles_step.ipp index 7e4f07630..6919082e3 100644 --- a/src/particles_step.ipp +++ b/src/particles_step.ipp @@ -105,20 +105,20 @@ namespace libcloudphxx assert(*thrust::min_element(pimpl->rhod.begin(), pimpl->rhod.end()) >= 0); // check if courants are greater than 2 since it would break the predictor-corrector (halo of size 2 in the x direction) - #if !defined(NDEBUG) if(!(pimpl->opts_init.adve_scheme != as_t::pred_corr || (courant_x.is_null() || ((*(thrust::min_element(pimpl->courant_x.begin(), pimpl->courant_x.end()))) >= real_t(-2.) )) )) { - std::cerr << "Courant x less than -2 in pred_corr advection, minimum courant x: " << *(thrust::min_element(pimpl->courant_x.begin(), pimpl->courant_x.end())) << " at " << (thrust::min_element(pimpl->courant_x.begin(), pimpl->courant_x.end())) - pimpl->courant_x.begin() << std::endl; - debug::print(pimpl->courant_x); - assert(false); +#if !defined(NDEBUG) + std::cerr << "Courant x less than -2 in pred_corr advection, minimum courant x: " << *(thrust::min_element(pimpl->courant_x.begin(), pimpl->courant_x.end())) << " at " << (thrust::min_element(pimpl->courant_x.begin(), pimpl->courant_x.end())) - pimpl->courant_x.begin() << " using first order advection in this step" << std::endl; +#endif + pimpl->adve_scheme = as_t::euler; } if(!(pimpl->opts_init.adve_scheme != as_t::pred_corr || (courant_x.is_null() || ((*(thrust::max_element(pimpl->courant_x.begin(), pimpl->courant_x.end()))) <= real_t(2.) )) )) { - std::cerr << "Courant x more than 2 in pred_corr advection, maximum courant x: " << *(thrust::max_element(pimpl->courant_x.begin(), pimpl->courant_x.end())) << " at " << (thrust::max_element(pimpl->courant_x.begin(), pimpl->courant_x.end())) - pimpl->courant_x.begin() << std::endl; - debug::print(pimpl->courant_x); - assert(false); +#if !defined(NDEBUG) + std::cerr << "Courant x more than 2 in pred_corr advection, maximum courant x: " << *(thrust::max_element(pimpl->courant_x.begin(), pimpl->courant_x.end())) << " at " << (thrust::max_element(pimpl->courant_x.begin(), pimpl->courant_x.end())) - pimpl->courant_x.begin() << " using first order advection in this step" << std::endl; +#endif + pimpl->adve_scheme = as_t::euler; } - #endif if (pimpl->opts_init.chem_switch){ for (int i = 0; i < chem_gas_n; ++i){ @@ -327,6 +327,8 @@ namespace libcloudphxx // advection, it invalidates i,j,k and ijk! if (opts.adve) pimpl->adve(); + // revert to the desired adve scheme (in case we used eulerian this timestep for halo reasons) + pimpl->adve_scheme = pimpl->opts_init.adve_scheme; // sedimentation has to be done after advection, so that negative z doesnt crash hskpng_ijk in adve if (opts.sedi)