better-faster transformer (#733)

* better-faster hierarchical transformer

CMakeLists.txt

@@ -554,13 +554,13 @@ if (ASGARD_BUILD_TESTS)
 
   # components with MPI-enabled testing
   set (mpi_test_components
-    adapt
-    distribution
-    time_advance
+    asgard_adapt
+    asgard_distribution
+    asgard_time_advance
   )
 
   if (ASGARD_USE_SCALAPACK)
-    list(APPEND mpi_test_components cblacs_grid fast_math scalapack_matrix_info scalapack_vector_info)
+    list(APPEND mpi_test_components asgard_cblacs_grid asgard_fast_math asgard_scalapack_matrix_info asgard_scalapack_vector_info)
   endif()
 
   foreach (component IN LISTS components)

src/asgard_adapt.cpp

@@ -67,26 +67,6 @@ remap_for_delete(std::vector<int64_t> const &deleted_indices,
   return new_to_old;
 }
 
-// helper to find new levels for each dimension after adapting table
-static std::vector<int>
-get_levels(elements::table const &adapted_table, int const num_dims)
-{
-  assert(num_dims > 0);
-  auto const flat_table = adapted_table.get_active_table();
-  auto const coord_size = num_dims * 2;
-  std::vector<int> max_levels(num_dims, 0);
-  for (int64_t element = 0; element < adapted_table.size(); ++element)
-  {
-    fk::vector<int, mem_type::const_view> coords(
-        flat_table, element * coord_size, (element + 1) * coord_size - 1);
-    for (auto i = 0; i < num_dims; ++i)
-    {
-      max_levels[i] = std::max(coords(i), max_levels[i]);
-    }
-  }
-  return max_levels;
-}
-
 template<typename P>
 static void update_levels(elements::table const &adapted_table, PDE<P> &pde,
                           bool const rechain = false)
@@ -129,98 +109,6 @@ distributed_grid<P>::distributed_grid(int max_level, prog_opts const &options)
   plan_ = get_plan(get_num_ranks(), table_);
 }
 
-template<typename P>
-fk::vector<P> distributed_grid<P>::get_initial_condition(
-    std::vector<dimension<P>> &dims, P const mult, int const num_terms,
-    std::vector<std::vector<term<P>>> &terms,
-    basis::wavelet_transform<P, resource::host> const &transformer,
-    prog_opts const &options)
-{
-  // get unrefined condition
-
-  // TODO: this needs to be refactored to allow dimensions to have different
-  // number of md_funcs
-  auto const num_md_funcs = dims.back().initial_condition.size();
-  std::vector<std::vector<vector_func<P>>> v_functions;
-  for (auto const &dim : dims)
-  {
-    // every dimension should have the same number of functions defined
-    expect(dim.initial_condition.size() == num_md_funcs);
-  }
-
-  for (size_t i = 0; i < num_md_funcs; i++)
-  {
-    v_functions.push_back(std::vector<vector_func<P>>());
-    for (auto const &dim : dims)
-    {
-      // add the ith function for this dimension
-      v_functions[i].push_back(dim.initial_condition[i]);
-    }
-  }
-
-  P const time             = 0;
-  auto const initial_unref = [this, &v_functions, &dims, &transformer, time,
-                              mult]() {
-    auto const subgrid     = this->get_subgrid(get_rank());
-    auto const vector_size = (subgrid.col_stop - subgrid.col_start + 1) *
-                             fm::ipow(dims[0].get_degree() + 1, dims.size());
-    fk::vector<P> initial(vector_size);
-    for (size_t i = 0; i < v_functions.size(); i++)
-    {
-      // TODO temp add scalar time func to initial conditions with multi-D func
-      // PR
-      auto const combined = transform_and_combine_dimensions(
-          dims, v_functions[i], this->get_table(), transformer,
-          subgrid.col_start, subgrid.col_stop, dims[0].get_degree(), time,
-          mult);
-      initial = initial + combined;
-    }
-    return initial;
-  };
-
-  if (not options.adapt_threshold)
-  {
-    return initial_unref();
-  }
-
-  // refine
-  fk::vector<P> refine_y(initial_unref());
-  auto refining = true;
-  while (refining)
-  {
-    auto const old_y   = fk::vector<P>(refine_y);
-    auto const refined = this->refine(old_y, options);
-    refining           = old_y.size() != refined.size();
-    // update_levels(this->get_table(), pde);
-    update_levels(this->get_table(), dims, num_terms, terms);
-
-    // reproject
-    refine_y = initial_unref();
-  }
-
-  // coarsen
-  auto const coarse_y = this->coarsen(refine_y, options);
-  update_levels(this->get_table(), dims, num_terms, terms);
-
-  return initial_unref();
-}
-
-template<typename P>
-void distributed_grid<P>::get_initial_condition(
-    std::vector<dimension<P>> const &dims,
-    std::vector<vector_func<P>> const &v_functions, P const mult,
-    basis::wavelet_transform<P, resource::host> const &transformer,
-    fk::vector<P, mem_type::view> result)
-{
-  // get unrefined condition
-  P const time       = 0;
-  auto const subgrid = this->get_subgrid(get_rank());
-  // TODO temp add scalar time func to initial conditions with multi-D func PR
-  transform_and_combine_dimensions(
-      dims, v_functions, this->get_table(), transformer, subgrid.col_start,
-      subgrid.col_stop, dims[0].get_degree(), time, mult, result);
-}
-
 template<typename P>
 fk::vector<P>
 distributed_grid<P>::coarsen_solution(PDE<P> &pde, fk::vector<P> const &x)

src/asgard_adapt.hpp

@@ -4,6 +4,26 @@
 
 namespace asgard::adapt
 {
+// helper to find new levels for each dimension after adapting table
+inline std::vector<int>
+get_levels(elements::table const &adapted_table, int const num_dims)
+{
+  assert(num_dims > 0);
+  auto const flat_table = adapted_table.get_active_table();
+  auto const coord_size = num_dims * 2;
+  std::vector<int> max_levels(num_dims, 0);
+  for (int64_t element = 0; element < adapted_table.size(); ++element)
+  {
+    fk::vector<int, mem_type::const_view> coords(
+        flat_table, element * coord_size, (element + 1) * coord_size - 1);
+    for (auto i = 0; i < num_dims; ++i)
+    {
+      max_levels[i] = std::max(coords(i), max_levels[i]);
+    }
+  }
+  return max_levels;
+}
+
 // this class bundles
 // 1) the element table (set of active elements and their coordinates) and
 // 2) the distribution plan that maps ranks to the active elements whose
@@ -40,29 +60,6 @@ class distributed_grid
       : distributed_grid(pde.max_level(), pde.options())
   {}
 
-  // driver routines
-  fk::vector<P> get_initial_condition(
-      PDE<P> &pde,
-      basis::wavelet_transform<P, resource::host> const &transformer)
-  {
-    return this->get_initial_condition(
-        pde.get_dimensions(),
-        pde.has_exact_time() ? pde.exact_time(0.0) : static_cast<P>(1.0),
-        pde.num_terms(), pde.get_terms(), transformer, pde.options());
-  }
-
-  fk::vector<P> get_initial_condition(
-      std::vector<dimension<P>> &dims, P const mult, int const num_terms,
-      std::vector<std::vector<term<P>>> &terms,
-      basis::wavelet_transform<P, resource::host> const &transformer,
-      prog_opts const &cli_opts);
-
-  void get_initial_condition(
-      std::vector<dimension<P>> const &dims,
-      std::vector<vector_func<P>> const &v_functions, P const mult,
-      basis::wavelet_transform<P, resource::host> const &transformer,
-      fk::vector<P, mem_type::view> result);
-
   fk::vector<P> coarsen_solution(PDE<P> &pde, fk::vector<P> const &x);
   fk::vector<P>
   refine_solution(PDE<P> &pde, fk::vector<P> const &x);

src/asgard_adapt_tests.cpp

@@ -170,21 +170,16 @@ void test_initial(prog_opts const &opts, std::string const &gold_filepath)
         [](P old_value) { return std::abs(old_value) < 1.e-14; }, 0.0);
     return raw;
   }();
-  auto const pde = make_PDE<P>(opts);
 
-  basis::wavelet_transform<P, resource::host> const transformer(
-      *pde, verbosity_level::quiet);
-  adapt::distributed_grid<P> adaptive_grid(*pde);
-  generate_dimension_mass_mat<P>(*pde, transformer);
+  discretization_manager<P> disc(make_PDE<P>(opts));
 
-  auto const test =
-      adaptive_grid.get_initial_condition(*pde, transformer);
+  fk::vector<P> const test = disc.current_state();
 
   REQUIRE(gold.size() >= test.size());
 
   auto constexpr tol_factor = get_tolerance<P>(100);
-  auto const my_subgrid     = adaptive_grid.get_subgrid(get_rank());
-  auto const segment_size   = element_segment_size(*pde);
+  auto const my_subgrid     = disc.get_grid().get_subgrid(get_rank());
+  auto const segment_size   = disc.get_hiermanip().block_size();
   fk::vector<P, mem_type::const_view> const my_gold(
       gold, my_subgrid.col_start * segment_size,
       (my_subgrid.col_stop + 1) * segment_size - 1);

src/asgard_coefficients.cpp

@@ -568,6 +568,8 @@ fk::matrix<P> generate_coefficients(
 }
 
 #ifdef ASGARD_ENABLE_DOUBLE
+template class coefficient_matrix_manager<double>;
+
 template fk::matrix<double> generate_coefficients<double>(
     dimension<double> const &dim, partial_term<double> const &pterm,
     basis::wavelet_transform<double, resource::host> const &transformer,
@@ -589,6 +591,8 @@ template void generate_dimension_mass_mat<double>(
 #endif
 
 #ifdef ASGARD_ENABLE_FLOAT
+template class coefficient_matrix_manager<float>;
+
 template fk::matrix<float> generate_coefficients<float>(
     dimension<float> const &dim, partial_term<float> const &pterm,
     basis::wavelet_transform<float, resource::host> const &transformer,

src/asgard_coefficients.hpp

@@ -23,4 +23,46 @@ fk::matrix<P> generate_coefficients(
     dimension<P> const &dim, partial_term<P> const &pterm,
     basis::wavelet_transform<P, resource::host> const &transformer,
     int const level, P const time = 0.0, bool const rotate = true);
+
+/*!
+ * \internal
+ * \brief Stores matrices associated with a pde
+ *
+ * Stores the matrices of a PDE discretization. Starting with the mass matrices,
+ * which are stored in block-diagonal format.
+ *
+ * \endinternal
+ */
+template<typename P>
+class coefficient_matrix_manager
+{
+public:
+  coefficient_matrix_manager(PDE<P> *pde)
+      : pde_(pde), num_dims_(pde_->num_dims())
+  {
+    for (int d : indexof<int>(num_dims_))
+      if (pde_->get_dimensions()[d].volume_jacobian_dV)
+      {
+        mass_[d] = std::make_unique<mass_matrix<P>>();
+        if (dv_funcs_.empty())
+          dv_funcs_.resize(num_dims_, nullptr); // nullptr means default to 1
+        dv_funcs_[d] = [=](std::vector<P> const &x, std::vector<P> &fx)
+            -> void {
+          for (auto i : indexof(x))
+            fx[i] = pde_->get_dimensions()[d].volume_jacobian_dV(x[i], 0);
+        };
+      }
+  }
+
+  auto &mass() const { return mass_; }
+  auto const &dv() const { return dv_funcs_; }
+
+private:
+  PDE<P> *pde_ = nullptr;
+
+  int num_dims_;
+  mutable std::array<std::unique_ptr<mass_matrix<P>>, max_num_dimensions> mass_;
+  std::vector<function_1d<P>> dv_funcs_;
+};
+
 } // namespace asgard