[obsolete] grid refinement compatibility

CMakeLists.txt

@@ -116,7 +116,7 @@ target_include_directories(cloudphxx_lgrngn
 # enabling c++14
 target_compile_features(cloudphxx_lgrngn PUBLIC cxx_std_14)
 # config-specific flags
-target_compile_options(cloudphxx_lgrngn PRIVATE $<$<AND:$<COMPILE_LANGUAGE:CXX>,$<CONFIG:DEBUG>>:                  -g -Og -DTHRUST_DEBUG>)
+target_compile_options(cloudphxx_lgrngn PRIVATE $<$<AND:$<COMPILE_LANGUAGE:CXX>,$<CONFIG:DEBUG>>:                  -g -Og -DTHRUST_DEBUG -Wfatal-errors>)
 target_compile_options(cloudphxx_lgrngn PRIVATE $<$<AND:$<COMPILE_LANGUAGE:CXX>,$<CONFIG:RELEASE>>:                -Ofast -march=native -Winline -DNDEBUG >)
 target_compile_options(cloudphxx_lgrngn PRIVATE $<$<AND:$<COMPILE_LANGUAGE:CXX>,$<CONFIG:RELWITHDEBINFO>>:         -Ofast -march=native>)
 target_compile_options(cloudphxx_lgrngn PRIVATE $<$<AND:$<COMPILE_LANGUAGE:CXX>,$<CONFIG:RELWITHDEBINFOPORTABLE>>: -Ofast>)

include/libcloudph++/lgrngn/opts_init.hpp

@@ -46,6 +46,14 @@ namespace libcloudphxx
       // Eulerian component parameters
       int nx, ny, nz;
       real_t dx, dy, dz, dt;
+      unsigned int n_ref; // number of grid refinements per 'normal' cell; refined cell size is (dx, dy, dz) / n_ref; 
+      // condensation (and sedimentation?) is resolved on refined cells, rest (coalescence, advection, etc.) on normal cells
+      // input arrays on the refined grid: pressure p (if supplied), rv, tht
+      // internal arrays on the refined grid: RH, T, eta
+      // TODO:
+      // dynamic viscosity eta on both grids?
+      // what about cell volume dv?
+      // what about density rhod?
 
       // no. of substeps 
       int sstp_cond, sstp_coal; 
@@ -65,7 +73,7 @@ namespace libcloudphxx
       // is it allowed to change dt during simulation through opts.dt
       bool variable_dt_switch;
 
-      // or, alternatively to sd_conc_mean, multiplicity of all SDs = const
+      // alternatively to sd_conc, multiplicity of all SDs = const
       unsigned long long  sd_const_multi;
 
       // max no. of super-droplets in the system
@@ -236,7 +244,8 @@ namespace libcloudphxx
         open_side_walls(false),
         periodic_topbot_walls(false),
         variable_dt_switch(false),
-        rng_seed_init_switch(false)
+        rng_seed_init_switch(false),
+        n_ref(1)
       {}
 
       // dtor (just to silence -Winline warnings)

src/detail/debug.hpp

@@ -0,0 +1,44 @@
+#pragma once
+
+namespace libcloudphxx
+{
+  namespace lgrngn
+  {
+//#if !defined(NDEBUG) // TODO (CMake defaults)
+    namespace debug
+    {
+      template <typename vec_t>
+      void print(const vec_t &v)
+      {
+        thrust::copy(v.begin(), v.end(), std::ostream_iterator<typename vec_t::value_type>(std::cerr, "\t"));
+        std::cerr << std::endl;
+      }
+
+      template <typename vec_bgn_t, typename vec_end_t>
+      void print(const vec_bgn_t &bgn, const vec_end_t &end)
+      {
+        thrust::copy(bgn, end, std::ostream_iterator<typename vec_bgn_t::value_type>(std::cerr, "\t"));
+        std::cerr << std::endl;
+      }
+
+      template<class T>
+      void print(ref_grid<T> v)
+      {
+        std::cerr << "normal grid values:" << std::endl;
+        print(v.cbegin(), v.cend());
+        std::cerr << "refined grid values:" << std::endl;
+        print(v.cbegin_ref(), v.cend_ref());
+      }
+
+      template<class T>
+      void print(ref_part<T> v)
+      {
+        std::cerr << "particle parameters on normal grid:" << std::endl;
+        print(v.cbegin(), v.cend());
+        std::cerr << "particle parameters on refined grid:" << std::endl;
+        print(v.cbegin_ref(), v.cend_ref());
+      }
+    };
+//#endif
+  };
+};

src/detail/grid_refinement_helpers.hpp

@@ -0,0 +1,230 @@
+// helper classes to store arrays associated with grid refinement
+#include <thrust/host_vector.h>
+
+#pragma once
+
+namespace libcloudphxx
+{
+  namespace lgrngn
+  {
+    template<class T>
+    class ref_common
+    {
+      public:
+      const bool ref_flag; // true if refinement is actually done
+
+      protected:
+      // ctor
+      ref_common(const bool ref_flag):
+        ref_flag(ref_flag)
+      {}
+
+      ref_common() = delete;
+    };
+
+
+    // for single values of type T
+    template<class T>
+    class ref_val : ref_common<T>
+    {
+      private:
+      using parent_t = ref_common<T>;
+      T val, val_ref;
+
+      public:
+      ref_val(T val, T val_ref):
+        val(val),
+        val_ref(val_ref),
+        parent_t(true)
+       {}
+      ref_val() = delete;
+
+      T& get()
+      {
+        return val;
+      }
+      /*
+      const T& get()
+      {
+        return val;
+      }
+      */
+      T& get_ref()
+      {
+        //return ref_flag ? val_ref : val;
+        return val_ref;
+      }
+      /*
+      const T& get_ref()
+      {
+        //return ref_flag ? val_ref : val;
+        return val_ref;
+      }
+      */
+    };
+
+    template<class T, bool host>
+    class ref_arr_common : public ref_common<T>
+    {
+      private:
+      using parent_t = ref_common<T>;
+      using vec_t = typename std::conditional<host, thrust::host_vector<T>, thrust_device::vector<T>>::type;
+
+      protected:
+      using parent_t::parent_t;
+      vec_t arr, arr_ref; // actual data, arr_ref initialized only if refinement is actually done
+
+      void resize(thrust_size_t size, thrust_size_t size_ref)
+      {
+        arr.resize(size);
+        if(this->ref_flag) arr_ref.resize(size_ref);
+      }
+      void resize(ref_val<thrust_size_t> size)
+      {
+        arr.resize(size.get());
+        if(this->ref_flag) arr_ref.resize(size.get_ref());
+      }
+      void reserve(thrust_size_t size, thrust_size_t size_ref)
+      {
+        arr.reserve(size);
+        if(this->ref_flag) arr_ref.reserve(size_ref);
+      }
+      void reserve(ref_val<thrust_size_t> size)
+      {
+        arr.reserve(size.get());
+        if(this->ref_flag) arr_ref.reserve(size.get_ref());
+      }
+
+      public:
+
+      auto begin()
+      {
+        return arr.begin();
+      }
+
+      const auto cbegin()
+      {
+        return arr.begin();
+      }
+
+      auto begin_ref()
+      {
+        return this->ref_flag ? arr_ref.begin() : arr.begin();
+      }
+
+      auto cbegin_ref()
+      {
+        return this->ref_flag ? arr_ref.begin() : arr.begin();
+      }
+
+      auto end()
+      {
+        return arr.end();
+      }
+
+      const auto cend()
+      {
+        return arr.end();
+      }
+
+      auto end_ref()
+      {
+        return this->ref_flag ? arr_ref.end() : arr.end();
+      }
+
+      auto cend_ref()
+      {
+        return this->ref_flag ? arr_ref.end() : arr.end();
+      }
+
+      auto ptr()
+      {
+        return &arr;
+      }
+
+      auto ptr_ref()
+      {
+        return &arr_ref;
+      }
+
+      vec_t& get()
+      {
+        return arr;
+      }
+
+      vec_t& get_ref()
+      {
+        return arr_ref;
+      }
+    };
+
+    // for arrays of the size of the grid that need values
+    // in the normal grid and in the refined grid (or in only one of these)
+    template <class T, bool host = false> // value type and container type
+    class ref_grid : public ref_arr_common<T, host>
+    {
+      using parent_t = ref_arr_common<T, host>;
+
+      public:
+      // ctor
+      ref_grid(const int &n_cell, const int &n_cell_ref):
+        parent_t(n_cell != n_cell_ref)
+      {
+        parent_t::resize(n_cell, n_cell_ref);
+      }
+
+      ref_grid(ref_val<thrust_size_t> n):
+        ref_grid(n.get(), n.get_ref())
+        {}
+
+      ref_grid() = delete;
+
+      // grid size cant be 0, ensure that correct arrays are pointed to
+      auto begin()
+      {
+        assert(this->arr.size() > 0);
+        return parent_t::begin();
+      }
+      auto begin_ref()
+      {
+        assert((this->ref_flag && this->arr_ref.size() > 0) || (!this->ref_flag && this->arr.size() > 0));
+        return parent_t::begin_ref();
+      }
+      auto end()
+      {
+        assert(this->arr.size() > 0);
+        return parent_t::end();
+      }
+      auto end_ref()
+      {
+        assert((this->ref_flag && this->arr_ref.size() > 0) || (!this->ref_flag && this->arr.size() > 0));
+        return parent_t::end_ref();
+      }
+    };
+
+    // for arrays of the size of the number of particles that need
+    // two independent values: for the normal and for the refined grid (e.g. ijk)
+    template <class T, bool host = false> // value type and container type
+    class ref_part : public ref_arr_common<T, host>
+    {
+      using parent_t = ref_arr_common<T, host>;
+
+      public:
+      // ctor
+      ref_part(const int &n_ref):
+        parent_t(n_ref > 1)
+      {assert(n_ref > 0);}
+
+      ref_part() = delete;
+
+      void resize(int n)
+      {
+        parent_t::resize(n, n);
+      }
+      void reserve(int n)
+      {
+        parent_t::reserve(n, n);
+      }
+    };
+  };
+};

src/detail/thrust.hpp

@@ -9,24 +9,5 @@ namespace libcloudphxx
   namespace lgrngn
   {
     typedef thrust_device::vector<int>::size_type thrust_size_t;
-
-//#if !defined(NDEBUG) // TODO (CMake defaults)
-    namespace debug
-    {
-      template <typename vec_t>
-      void print(const vec_t &v)
-      {
-	thrust::copy(v.begin(), v.end(), std::ostream_iterator<typename vec_t::value_type>(std::cerr, "\t"));
-	std::cerr << std::endl;
-      }
-
-      template <typename vec_bgn_t, typename vec_end_t>
-      void print(const vec_bgn_t &bgn, const vec_end_t &end)
-      {
-	thrust::copy(bgn, end, std::ostream_iterator<typename vec_bgn_t::value_type>(std::cerr, "\t"));
-	std::cerr << std::endl;
-      }
-    };
-//#endif
   };
 };