Better landau (#755)

* updated moments for the 1x1v case

CMakeLists.txt

@@ -382,6 +382,7 @@ set (components
      asgard_block_matrix
      asgard_boundary_conditions
      asgard_coefficients
+     asgard_coefficients_mats
      asgard_distribution
      asgard_discretization
      asgard_elements

src/asgard_coefficients_mats_tests.cpp

@@ -0,0 +1,103 @@
+#include "tests_general.hpp"
+
+#include "asgard_coefficients_mats.hpp"
+
+using namespace asgard;
+
+// functors have signature P(P)
+template<typename P, typename func_lhs, typename func_rhs>
+std::array<P, 2> test_mass_coeff_moments(func_lhs flhs, func_rhs frhs,
+                                         P xleft, P xright, int level, int degree)
+{
+  dimension<P> const dim(xleft, xright, level, degree, nullptr, nullptr, "dim");
+  hierarchy_manipulator<P> hier(degree, 1, {xleft, }, {xright, });
+
+  block_diag_matrix<P> ref; // reference matrix
+  gen_diag_cmat<P, coefficient_type::mass>(dim, level, 0, [&](int, P x, P)->P{ return frhs(x); }, ref);
+
+  partial_term<P> const pterm(coefficient_type::mass, nullptr,
+                              [&](P x, P)->P{ return flhs(x); });
+  level_mass_matrces<P> mass;
+  generate_partial_mass<P>(0, dim, pterm, hier, 0, mass);
+
+  if (mass.has_level(level))
+    invert_mass(degree + 1, mass[level], ref);
+
+  // project the lhs and rhs functions
+  std::vector<P> vlhs = hier.cell_project(
+    [&](std::vector<P> const &x, std::vector<P> &fx) {
+      for (auto i : indexof(x))
+        fx[i] = flhs(x[i]);
+    }, nullptr, level);
+  std::vector<P> vrhs = hier.cell_project(
+    [&](std::vector<P> const &x, std::vector<P> &fx) {
+      for (auto i : indexof(x))
+        fx[i] = frhs(x[i]);
+    }, nullptr, level);
+
+  // combine the two vectors into one, use two versions for the moments
+  int const pdof      = degree + 1;
+  int const num_cells = fm::ipow2(level);
+  std::vector<P> mom2(2 * pdof * num_cells), mom3(3 * pdof * num_cells);
+  for (int i : iindexof(num_cells))
+  {
+    std::copy_n(vlhs.data() + i * pdof, pdof, mom2.data() + 2 * i * pdof);
+    std::copy_n(vlhs.data() + i * pdof, pdof, mom3.data() + 3 * i * pdof);
+
+    std::copy_n(vrhs.data() + i * pdof, pdof, mom2.data() + 2 * i * pdof + pdof);
+    std::copy_n(vrhs.data() + i * pdof, pdof, mom3.data() + 3 * i * pdof + 2 * pdof);
+  }
+
+  block_diag_matrix<P> comp2, comp3;
+  gen_diag_mom_cmat<P, coefficient_type::mass, 1>(dim, level, 0, 2, mom2, comp2);
+  gen_diag_mom_cmat<P, coefficient_type::mass, 2>(dim, level, 0, 3, mom3, comp3);
+
+  block_diag_matrix<P> comp2_div, comp3_div;
+  gen_diag_mom_cmat_div<P, coefficient_type::mass, 1>(dim, level, 0, 2, mom2, comp2_div);
+  gen_diag_mom_cmat_div<P, coefficient_type::mass, 2>(dim, level, 0, 3, mom3, comp3_div);
+
+  // std::cout << " ref \n";
+  // ref.to_full().print(std::cout);
+  // std::cout << " comp2 \n";
+  // comp2.to_full().print(std::cout);
+  // std::cout << " comp3 \n";
+  // comp2_div.to_full().print(std::cout);
+
+  P const err2 = std::max(ref.to_full().max_diff(comp2.to_full()),
+                          ref.to_full().max_diff(comp2_div.to_full()));
+  P const err3 = std::max(ref.to_full().max_diff(comp3.to_full()),
+                          ref.to_full().max_diff(comp3_div.to_full()));
+
+  // std::cout << " err = " << err2 << "  " << err3 << "\n";
+
+  return {err2, err3};
+}
+
+TEMPLATE_TEST_CASE("projected mass", "[mass-coefficients]", test_precs)
+{
+  // compare the construction between a mass term constructed from simple function
+  // as opposed to the projection of the functions onto the Legendre basis
+  using P = TestType;
+  P constexpr tol = (std::is_same_v<P, double>) ? 1.E-11 : 1.E-4;
+
+  std::array<P, 2> err = {0, 0};
+
+  err = test_mass_coeff_moments<P>([](P)->P{ return 1; }, [](P)->P{ return -1; }, 0, 1, 2, 0);
+  REQUIRE(std::max(err[0], err[1]) < tol);
+
+  err = test_mass_coeff_moments<P>([](P)->P{ return 1; }, [](P)->P{ return -1; }, 0, 1, 5, 2);
+  REQUIRE(std::max(err[0], err[1]) < tol);
+
+  err = test_mass_coeff_moments<P>([](P)->P{ return 2; }, [](P)->P{ return 1; }, 0, 1, 5, 2);
+  REQUIRE(std::max(err[0], err[1]) < tol);
+
+  err = test_mass_coeff_moments<P>([](P)->P{ return 1; },
+                                   [](P x)->P{ return std::sin(x); },
+                                   -3.14, 3.14, 3, 0);
+  REQUIRE(std::max(err[0], err[1]) < tol);
+
+  err = test_mass_coeff_moments<P>([](P x)->P{ return 1 + 0.1 * std::sin(x); },
+                                   [](P x)->P{ return 1 + x + 0.1 *std::cos(x); },
+                                   -3.14, 3.14, 7, 1);
+  REQUIRE(std::max(err[0], err[1]) < 5.E-5);
+}

src/asgard_coefficients_tests.cpp

@@ -313,27 +313,6 @@ TEMPLATE_TEST_CASE("fokkerplanck2_complete_case4 terms", "[coefficients]",
   }
 }
 
-TEMPLATE_TEST_CASE("vlasov terms", "[coefficients]", test_precs)
-{
-  auto const gold_path =
-      coefficients_base_dir / "vlasov_lb_full_f_coefficients";
-
-  TestType const tol_factor = std::is_same_v<TestType, double> ? 1e-12 : 1e-3;
-
-  prog_opts opts;
-  opts.pde_choice     = PDE_opts::vlasov_lb_full_f;
-  opts.grid           = grid_type::dense;
-  opts.num_time_steps = 1;
-
-  SECTION("level [4,3], degree 2")
-  {
-    opts.start_levels = {4, 3};
-    opts.degree       = 2;
-
-    test_coefficients<TestType>(opts, gold_path, tol_factor);
-  }
-}
-
 template<typename P>
 class penalty_pde : public PDE<P>
 {

src/asgard_discretization.cpp

@@ -71,8 +71,6 @@ discretization_manager<precision>::discretization_manager(
 
     matrices.edata.electric_field.resize(fm::ipow2(dim.get_level()));
 
-    moms1d = moments1d<precision>(1, degree_, pde->max_level(), pde->get_dimensions());
-
     // if the inf-nrm is needed, initialize with a dummy value
     for (int d : indexof<int>(pde->num_dims()))
       for (int t : indexof<int>(pde->num_terms()))
@@ -82,8 +80,23 @@ discretization_manager<precision>::discretization_manager(
         }
   }
 
+  matrices.edata.num_moments = pde->required_moments();
+  if (matrices.edata.num_moments > 0) {
+    moms1d = moments1d<precision>(matrices.edata.num_moments, degree_,
+                                  pde->max_level(), pde->get_dimensions());
+    int const level = pde->get_dimensions().front().get_level();
+    moms1d->project_moments(level, state, grid.get_table(), matrices.edata.moments);
+    int const num_cells = fm::ipow2(level);
+    int const num_moms  = matrices.edata.num_moments;
+    hier.reconstruct1d(
+        num_moms, level, span2d<precision>((degree_ + 1), num_moms * num_cells,
+                                            matrices.edata.moments.data()));
+  }
+
   this->compute_coefficients();
 
+  comp_mats();
+
   auto const msg = grid.get_subgrid(get_rank());
   fixed_bc = boundary_conditions::make_unscaled_bc_parts(
         *pde, grid.get_table(), transformer, hier, matrices, conn, msg.row_start, msg.row_stop);

src/asgard_discretization.hpp

@@ -215,6 +215,24 @@ class discretization_manager
   std::vector<moment<precision>> &get_moments() const { return moments; }
   //! returns the coefficient matrices
   coefficient_matrices<precision> &get_cmatrices() const { return matrices; }
+  //! recomputes the moments given the state of interest
+  void compute_moments(std::vector<precision> const &f) {
+    if (moms1d) {
+      int const level = pde->get_dimensions().front().get_level();
+      moms1d->project_moments(level, f, grid.get_table(), matrices.edata.moments);
+      int const num_cells = fm::ipow2(level);
+      int const num_moms  = moms1d->num_mom();
+      hier.reconstruct1d(
+          num_moms, level, span2d<precision>((degree_ + 1), num_moms * num_cells,
+                                             matrices.edata.moments.data()));
+    }
+  }
+  //! (testing) recomputes the moments given the state of interest, keeps in hierarchical form
+  void compute_hmoments(std::vector<precision> const &f, std::vector<precision> &rmom) {
+    if (moms1d)
+      moms1d->project_moments(pde->get_dimensions().front().get_level(),
+                              f, grid.get_table(), rmom);
+  }
   //! recomputes the coefficients, can select sub
   void compute_coefficients(coeff_update_mode mode = coeff_update_mode::all) {
     generate_coefficients(*pde, matrices, conn, hier, time_, mode);
@@ -239,9 +257,15 @@ class discretization_manager
   void comp_mats() const { // two stream, compare matrices
     return;
     // reference check for operator construction, helps find problems
-    // auto ref = matrices.term_coeffs[4 + i].to_full(conn);
-    // auto com = matrices.term_coeffs[8 + i].to_full(conn);
+    auto ref = matrices.term_coeffs[8].to_full(conn);
+    auto com = matrices.term_coeffs[10].to_full(conn);
     // precision err = std::max(err, ref.max_diff(com));
+
+    std::cout << "  -- ref -- \n";
+    ref.print(std::cout);
+    std::cout << "  -- comp -- \n";
+    com.print(std::cout);
+    std::cout << "  -- ---- -- \n";
   }
 
 protected:

src/asgard_kronmult_matrix.hpp

@@ -29,6 +29,10 @@ struct coupled_term_data
   std::vector<P> electric_field;
   //! max-absolute value of the electric field
   std::optional<P> electric_field_infnrm;
+  //! number of computed moments
+  int num_moments = 0;
+  //! data for the computed moments
+  std::vector<P> moments;
 };
 
 /*!

src/asgard_moment.hpp

@@ -25,6 +25,8 @@ static constexpr resource sparse_resrc = resource::host;
 template<typename P>
 class moments1d {
 public:
+  //! empty constructor, no moments
+  moments1d() {}
   //! constructor, prepares the given number of momemnts, for dgree and up to the max_level
   moments1d(int num_mom, int degree, int max_level, std::vector<dimension<P>> const &dims);
 
@@ -46,6 +48,9 @@ class moments1d {
   void project_moment(int const mom, int const dim0_level, std::vector<P> const &state,
                       elements::table const &etable, std::vector<P> &moment) const;
 
+  //! \brief Returns the number of loaded moments
+  int num_mom() const { return num_mom_; }
+
 protected:
   /*!
    * \brief Computes the moment integrals over a sub-range of the domain
@@ -75,11 +80,11 @@ class moments1d {
 
 private:
   //! number of momemnts
-  int num_mom_;
+  int num_mom_ = 0;
   //! number of dimensions
-  int num_dims_;
+  int num_dims_ = 0;
   //! the degree of the basis
-  int degree_;
+  int degree_ = 0;
   //! ingeral of the canonical basis, each index holds num_mom_ * (degree_ + 1) entries
   std::array<vector2d<P>, max_num_dimensions> integ;
 };

src/asgard_moment_tests.cpp

@@ -233,48 +233,3 @@ struct distribution_test_init
 #ifdef ASGARD_USE_MPI
 static distribution_test_init const distrib_test_info;
 #endif
-
-TEMPLATE_TEST_CASE("CreateMomentReducedMatrix", "[moments]", test_precs)
-{
-  std::string const pde_choice = "vlasov";
-  fk::vector<int> const levels{4, 3};
-  auto constexpr tol_factor = get_tolerance<TestType>(100);
-
-  prog_opts opts;
-  opts.pde_choice     = PDE_opts::vlasov_lb_full_f;
-  opts.start_levels   = {4, 3};
-  opts.degree         = 2;
-  opts.grid           = grid_type::dense;
-  opts.num_time_steps = 1;
-
-  discretization_manager<TestType> disc(make_PDE<TestType>(opts),
-                                        verbosity_level::quiet);
-
-  auto &pde = disc.get_pde();
-
-  auto &moments = disc.get_moments();
-  REQUIRE(moments.size() > 0);
-
-  for (size_t i = 0; i < moments.size(); ++i)
-  {
-    moments[i].createFlist(pde);
-    moments[i].createMomentVector(pde, disc.get_grid().get_table());
-    moments[i].createMomentReducedMatrix(pde, disc.get_grid().get_table());
-
-    auto const gold_filename =
-        moment_base_dir /
-        ("moment_matrix_vlasov_d3_l4_3_m" + std::to_string(i + 1) + ".dat");
-    auto const gold_moment_matrix =
-        read_matrix_from_txt_file<TestType>(gold_filename);
-
-#ifdef ASGARD_USE_CUDA
-    rmse_comparison(
-        gold_moment_matrix,
-        moments[i].get_moment_matrix_dev().clone_onto_host().to_dense(),
-        tol_factor);
-#else
-    rmse_comparison(gold_moment_matrix,
-                    moments[i].get_moment_matrix_dev().to_dense(), tol_factor);
-#endif
-  }
-}

src/asgard_small_mats.hpp

@@ -84,6 +84,32 @@ void scal(int const &n, P alpha, P x[])
   for (int i = 0; i < n; i++)
     x[i] *= alpha;
 }
+//! scale x by alpha, n is the size of x, write the result in y
+template<typename P>
+void scal(int const &n, P alpha, P const x[], P y[])
+{
+  ASGARD_OMP_SIMD
+  for (int i = 0; i < n; i++)
+    y[i] = alpha * x[i];
+}
+//! B will be replaced by entry-wise mulitplication of each columns of A (nr x nc) by vector x
+template<typename P>
+void col_scal(int const &nr, int const &nc, P const x[], P const A[], P B[])
+{
+  ASGARD_PRAGMA_OMP_SIMD(collapse(2))
+  for (int c = 0; c < nc; c++)
+    for (int r = 0; r < nr; r++)
+      B[c * nr + r] = A[c * nr + r] * x[r];
+}
+//! entry-wise mulitplication of each columns of A (nr x nc) by vector alpha * x
+template<typename P>
+void col_scal(int const &nr, int const &nc, P alpha, P const x[], P A[])
+{
+  ASGARD_PRAGMA_OMP_SIMD(collapse(2))
+  for (int c = 0; c < nc; c++)
+    for (int r = 0; r < nr; r++)
+      A[c * nr + r] *= alpha * x[r];
+}
 //! matrix-vector multiplication y += A * x, A has size nr X nc
 template<typename P>
 void gemv1(int const &nr, int const &nc, P const A[], P const x[], P y[])