Merge pull request #367 from pdziekan/fix_pred_corr_large_C

pred_corr: if C>2 use euler this step

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<real_t, device> rng;
       detail::config<real_t> 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<real_t, n_t> *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

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<detail::adve_helper_expl<real_t> >(true, halo_x);
         return;
       }
-      else if(opts_init.adve_scheme == as_t::implicit)
+      else if(adve_scheme == as_t::implicit)
       {
         adve_calc<detail::adve_helper_impl<real_t> >(true, halo_x);
         return;

src/impl/particles_impl_hskpng_ijk.ipp

@@ -38,11 +38,11 @@ namespace libcloudphxx
           thrust_device::vector<thrust_size_t> &vi,
           const real_t &vd
         ) {
-	  thrust::transform(
-	    vx.begin(), vx.end(),                                // input
-	    vi.begin(),                                          // output
-	    detail::divide_by_constant_and_cast<double, thrust_size_t>(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<double, thrust_size_t>(vd) // has to be done on doubles to avoid i==nx due to low precision of nvcc math
+          );
         }
       } helper;
 

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(

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)