updated the meaning of degree (#715)

one degree less

examples/continuity_2d.m

@@ -13,7 +13,7 @@ function example_continuity_2d()
 end
 
 % command to run the executable and generate the output file
-command = ['./example_continuity_2d -p continuity_2 -d 2 -l 6 -w 10 -n 10 -t 0.0001'];
+command = ['./example_continuity_2d -p continuity_2 -d 1 -l 6 -w 10 -n 10 -t 0.0001'];
 
 [status, cmdout] = system(command, '-echo');
 

examples/continuity_2d.py

@@ -24,7 +24,7 @@
         exit(1)
 
     print("asgard: running the continuity example")
-    os.system("./example_continuity_2d -p continuity_2 -d 2 -l 6 -w 10 -n 10 -t 0.0001")
+    os.system("./example_continuity_2d -p continuity_2 -d 1 -l 6 -w 10 -n 10 -t 0.0001")
 
     # the example above will run for 10 time steps and the -w 10 options
     # will tell the code to output on the final 10-th step

python/asgard.py

@@ -45,9 +45,8 @@ def __init__(self, filename, verbose = False):
             # print(fdata.keys())
 
             self.num_dimensions = fdata['ndims'][()]
-            self.pterms         = fdata['degree'][()]
-            self.porder         = self.pterms - 1
-            assert self.porder == 1, "only works with linear basis, others will be coming soon"
+            self.degree         = fdata['degree'][()]
+            assert self.degree == 1, "only works with linear basis, others will be coming soon"
 
             self.state = fdata['soln'][()]
             self.cells = fdata['elements'][()]
@@ -84,12 +83,12 @@ def __init__(self, filename, verbose = False):
             self.double_precision = True
             self.recsol = libasgard.asgard_make_dreconstruct_solution(
                 self.num_dimensions, self.num_cells, np.ctypeslib.as_ctypes(self.cells.reshape(-1,)),
-                self.pterms, np.ctypeslib.as_ctypes(self.state.reshape(-1,)))
+                self.degree, np.ctypeslib.as_ctypes(self.state.reshape(-1,)))
         else:
             self.double_precision = False
             self.recsol = libasgard.asgard_make_freconstruct_solution(
                 self.num_dimensions, self.num_cells, np.ctypeslib.as_ctypes(self.cells.reshape(-1,)),
-                self.pterms, np.ctypeslib.as_ctypes(self.state.reshape(-1,)))
+                self.degree, np.ctypeslib.as_ctypes(self.state.reshape(-1,)))
 
         libasgard.asgard_reconstruct_solution_setbounds(self.recsol,
                                                         np.ctypeslib.as_ctypes(self.dimension_min.reshape(-1,)),

python/pyasgard_test.py

@@ -60,7 +60,7 @@ def onedim_match_basis(self, ibasis, basis):
         libasgard = asgard.libasgard
         pntr = libasgard.asgard_make_dreconstruct_solution(
                 1, num_cells, np.ctypeslib.as_ctypes(cells),
-                2, np.ctypeslib.as_ctypes(state))
+                1, np.ctypeslib.as_ctypes(state))
 
         libasgard.asgard_reconstruct_solution_setbounds(
             pntr, np.ctypeslib.as_ctypes(dmin), np.ctypeslib.as_ctypes(dmax))
@@ -101,7 +101,7 @@ def twodim_match_basis(self, ibasis, basis):
         libasgard = asgard.libasgard
         pntr = libasgard.asgard_make_dreconstruct_solution(
                 2, num_cells, np.ctypeslib.as_ctypes(cells),
-                2, np.ctypeslib.as_ctypes(state))
+                1, np.ctypeslib.as_ctypes(state))
 
         libasgard.asgard_reconstruct_solution_setbounds(
             pntr, np.ctypeslib.as_ctypes(dmin), np.ctypeslib.as_ctypes(dmax))
@@ -151,7 +151,7 @@ def test2d_basis_reconstruct(self):
     # simple IO test
     def test_simple1d(self):
         print("\ntesting 1d plot")
-        os.system("./asgard -p continuity_1 -d 2 -l 6 -w 10 -t 0.01 1>/dev/null")
+        os.system("./asgard -p continuity_1 -d 1 -l 6 -w 10 -t 0.01 1>/dev/null")
 
         self.assertTrue(os.path.isfile("asgard_wavelet_10.h5"), "failed to run continuity_1")
 
@@ -178,7 +178,7 @@ def test_simple1d(self):
 
     def test_simple2d(self):
         print("\ntesting 2d plot")
-        os.system("./asgard -p continuity_2 -d 2 -l 6 -w 10 -t 0.00001 1>/dev/null")
+        os.system("./asgard -p continuity_2 -d 1 -l 6 -w 10 -t 0.00001 1>/dev/null")
 
         self.assertTrue(os.path.isfile("asgard_wavelet_10.h5"), "failed to run continuity_2")
 
@@ -211,7 +211,7 @@ def test_simple2d(self):
 
     def test_simple3d(self):
         print("\ntesting 3d plot (longer test)")
-        os.system("./asgard -p continuity_3 -d 2 -l 8 -w 10 -t 0.0001 1>/dev/null")
+        os.system("./asgard -p continuity_3 -d 1 -l 8 -w 10 -t 0.0001 1>/dev/null")
 
         self.assertTrue(os.path.isfile("asgard_wavelet_10.h5"), "failed to run continuity_3")
 

src/adapt.cpp

@@ -132,7 +132,6 @@ distributed_grid<P>::distributed_grid(options const &cli_opts,
   plan_ = get_plan(get_num_ranks(), table_);
 }
 
-// FIXME assumes uniform degree across levels
 template<typename P>
 fk::vector<P> distributed_grid<P>::get_initial_condition(
     std::vector<dimension<P>> &dims, P const mult, int const num_terms,
@@ -167,7 +166,7 @@ fk::vector<P> distributed_grid<P>::get_initial_condition(
                               mult]() {
     auto const subgrid     = this->get_subgrid(get_rank());
     auto const vector_size = (subgrid.col_stop - subgrid.col_start + 1) *
-                             std::pow(dims[0].get_degree(), dims.size());
+                             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++)
     {

src/adapt.hpp

@@ -21,8 +21,8 @@ namespace asgard::adapt
 // the purpose of this class is to adapt the set of the active elements given
 // the coefficients in the initial condition/solution vector x. during gemv that
 // drives explicit time advance, for each term, each element connection reads
-// from a deg^dim segment x beginning at x[i*deg^dim] and writes deg^dim
-// coefficients to y[j*deg^dim].
+// from a n^dim segment x beginning at x[i*n^dim] and writes n^dim
+// coefficients to y[j*n^dim], where n = degree + 1.
 
 // elements responsible for coefficients with low absolute value may be removed
 // from the grid (grid coarsening) elements responsible for coefficients with

src/adapt_tests.cpp

@@ -119,7 +119,7 @@ TEMPLATE_TEST_CASE("adapt - 1d, scattered coarsen/refine", "[adapt]",
     return;
   }
 
-  auto const degree     = 3;
+  auto const degree     = 2;
   auto const level      = 4;
   auto const pde_choice = PDE_opts::continuity_1;
   auto const num_dims   = 1;
@@ -137,7 +137,7 @@ TEMPLATE_TEST_CASE("adapt - 2d, all zero", "[adapt]", test_precs)
   {
     return;
   }
-  auto const degree     = 2;
+  auto const degree     = 1;
   auto const level      = 5;
   auto const pde_choice = PDE_opts::continuity_2;
   auto const num_dims   = 2;
@@ -161,7 +161,7 @@ TEMPLATE_TEST_CASE("adapt - 3d, scattered, contiguous refine/adapt", "[adapt]",
   {
     return;
   }
-  auto const degree     = 4;
+  auto const degree     = 3;
   auto const level      = 4;
   auto const pde_choice = PDE_opts::continuity_3;
   auto const num_dims   = 3;
@@ -209,7 +209,7 @@ void test_initial(parser const &problem, std::string const &gold_filepath)
 
 TEMPLATE_TEST_CASE("initial - diffusion 1d", "[adapt]", test_precs)
 {
-  auto const degree     = 4;
+  auto const degree     = 3;
   auto const level      = 3;
   auto const pde_choice = PDE_opts::diffusion_1;
   auto const num_dims   = 1;
@@ -234,7 +234,7 @@ TEMPLATE_TEST_CASE("initial - diffusion 2d", "[adapt]", test_precs)
   {
     return;
   }
-  auto const degree     = 3;
+  auto const degree     = 2;
   auto const level      = 2;
   auto const pde_choice = PDE_opts::diffusion_2;
   auto const num_dims   = 2;

src/asgard.cpp

@@ -39,7 +39,23 @@ void simulate(parser const &cli_input, std::unique_ptr<PDE<precision>> &pde)
   node_out() << "ASGarD problem configuration:" << '\n';
   node_out() << "  selected PDE: " << cli_input.get_pde_string()
              << '\n';
-  node_out() << "  degree: " << degree << '\n';
+  switch (degree)
+  {
+  case 0:
+    node_out() << "  degree: constant (0) \n";
+    break;
+  case 1:
+    node_out() << "  degree: linear (1) \n";
+    break;
+  case 2:
+    node_out() << "  degree: quadratic (2) \n";
+    break;
+  case 3:
+    node_out() << "  degree: cubic (3) \n";
+    break;
+  default:
+    node_out() << "  degree: " << degree << '\n';
+  };
   node_out() << "  N steps: " << opts.num_time_steps << '\n';
   node_out() << "  write freq: " << opts.wavelet_output_freq << '\n';
   node_out() << "  realspace freq: " << opts.realspace_output_freq
@@ -68,8 +84,7 @@ void simulate(parser const &cli_input, std::unique_ptr<PDE<precision>> &pde)
 
   adapt::distributed_grid adaptive_grid(*pde, opts);
   node_out() << "  degrees of freedom: "
-             << adaptive_grid.size() * static_cast<uint64_t>(std::pow(
-                                           degree, pde->num_dims()))
+             << adaptive_grid.size() * fm::ipow(degree + 1, pde->num_dims())
              << '\n';
 
   node_out() << "  generating: basis operator..." << '\n';
@@ -86,8 +101,7 @@ void simulate(parser const &cli_input, std::unique_ptr<PDE<precision>> &pde)
   auto initial_condition =
       adaptive_grid.get_initial_condition(*pde, transformer, opts);
   node_out() << "  degrees of freedom (post initial adapt): "
-                     << adaptive_grid.size() * static_cast<uint64_t>(std::pow(
-                                                   degree, pde->num_dims()))
+                     << adaptive_grid.size() * fm::ipow(degree + 1, pde->num_dims())
                      << '\n';
 
   // -- regen mass mats after init conditions - TODO: check dims/rechaining?
@@ -173,7 +187,7 @@ void simulate(parser const &cli_input, std::unique_ptr<PDE<precision>> &pde)
   std::vector<size_t> sizes(pde->num_dims);
   for (int i = 0; i < pde->num_dims; i++)
   {
-    sizes[i] = pde->get_dimensions()[i].get_degree() *
+    sizes[i] = (pde->get_dimensions()[i].get_degree() + 1) *
                fm::two_raised_to(pde->get_dimensions()[i].get_level());
   }
   ml_plot.set_var("initial_condition",

src/asgard_dimension.hpp

@@ -56,7 +56,7 @@ struct dimension
 
     for (int i = 0; i <= level_; ++i)
     {
-      auto const max_dof = fm::two_raised_to(i) * degree_;
+      auto const max_dof = fm::two_raised_to(i) * (degree + 1);
       expect(max_dof < INT_MAX);
       this->mass_.push_back(eye<P>(max_dof));
     }
@@ -84,7 +84,7 @@ struct dimension
 
   void set_degree(int const degree)
   {
-    expect(degree > 0);
+    expect(degree >= 0);
     degree_ = degree;
   }
 

src/asgard_grid_1d.hpp

@@ -46,13 +46,13 @@ class connect_1d
   }
   /*!
    * \brief Creates a new connectivity matrix by expanding each element with
-   *        a block of size (porder+1) by (porder+1)
+   *        a block of size (degree+1) by (degree+1)
    */
-  connect_1d(connect_1d const &elem_connect, int porder)
-      : levels(-1), rows(elem_connect.rows * (porder + 1)),
+  connect_1d(connect_1d const &elem_connect, int degree)
+      : levels(-1), rows(elem_connect.rows * (degree + 1)),
         pntr(rows + 1, 0), diag(rows)
   {
-    int const block_rows = porder + 1;
+    int const block_rows = degree + 1;
 
     pntr[0] = 0;
     for (int row = 0; row < elem_connect.num_rows(); row++)

src/asgard_indexset.cpp

@@ -324,13 +324,13 @@ indexset compute_ancestry_completion(indexset const &iset,
  * \brief Helper method, fills the indexes with the polynomial degree of freedom
  *
  * The cells are the current set of cells to process,
- * pterm is the number of polynomial terms,
+ * pdof is the number of polynomial terms,
  * e.g., 2 for linear and 3 for quadratic.
  * tsize is the size of the tensor within a cell,
- * i.e., tsize = pterms to power num_dimensions
+ * i.e., tsize = pdof to power num_dimensions
  */
 template<typename itype>
-void complete_poly_order(span2d<itype> const &cells, int64_t pterms,
+void complete_poly_order(span2d<itype> const &cells, int64_t pdof,
                          int64_t tsize, span2d<int> indexes)
 {
   int num_dimensions = cells.stride();
@@ -348,59 +348,55 @@ void complete_poly_order(span2d<itype> const &cells, int64_t pterms,
 
       for (int d = num_dimensions - 1; d >= 0; d--)
       {
-        idx[d] = cell[d] * pterms + static_cast<int>(t % pterms);
-        t /= pterms;
+        idx[d] = cell[d] * pdof + static_cast<int>(t % pdof);
+        t /= pdof;
       }
     }
   }
 }
 
-vector2d<int> complete_poly_order(vector2d<int> const &cells, int porder)
+vector2d<int> complete_poly_order(vector2d<int> const &cells, int degree)
 {
-  int num_dimensions = cells.stride();
+  int const num_dimensions = cells.stride();
 
-  int64_t num_cells = cells.num_strips();
+  int64_t const num_cells = cells.num_strips();
 
-  int64_t pterms = porder + 1;
+  int64_t const pdof = degree + 1;
 
-  int64_t tsize = pterms;
-  for (int64_t d = 1; d < num_dimensions; d++)
-    tsize *= pterms;
+  int64_t const tsize = fm::ipow(pdof, num_dimensions);
 
   vector2d<int> indexes(num_dimensions, tsize * num_cells);
 
   complete_poly_order(
-      span2d(num_dimensions, num_cells, cells[0]), pterms, tsize,
+      span2d(num_dimensions, num_cells, cells[0]), pdof, tsize,
       span2d(num_dimensions, tsize * num_cells, indexes[0]));
 
   return indexes;
 }
 
 vector2d<int> complete_poly_order(vector2d<int> const &cells,
-                                  indexset const &padded, int porder)
+                                  indexset const &padded, int degree)
 {
   expect(padded.num_indexes() == 0 or padded.num_dimensions() == cells.stride());
 
   int num_dimensions = cells.stride();
 
-  int64_t num_cells  = cells.num_strips();
-  int64_t num_padded = padded.num_indexes();
+  int64_t const num_cells  = cells.num_strips();
+  int64_t const num_padded = padded.num_indexes();
 
-  int64_t pterms = porder + 1;
+  int64_t const pdof = degree + 1;
 
-  int64_t tsize = pterms;
-  for (int64_t d = 1; d < num_dimensions; d++)
-    tsize *= pterms;
+  int64_t const tsize = fm::ipow(pdof, num_dimensions);
 
   vector2d<int> indexes(num_dimensions, tsize * (num_cells + num_padded));
 
   complete_poly_order(
-      span2d(num_dimensions, num_cells, cells[0]), pterms, tsize,
+      span2d(num_dimensions, num_cells, cells[0]), pdof, tsize,
       span2d(num_dimensions, tsize * num_cells, indexes[0]));
 
   if (num_padded > 0)
     complete_poly_order(
-        span2d(num_dimensions, num_padded, padded[0]), pterms, tsize,
+        span2d(num_dimensions, num_padded, padded[0]), pdof, tsize,
         span2d(num_dimensions, tsize * num_padded, indexes[tsize * num_cells]));
 
   return indexes;

src/asgard_indexset.hpp

@@ -405,20 +405,16 @@ indexset compute_ancestry_completion(indexset const &iset,
  *
  * Given the cells, the returned list of indexes
  * will hold all indexes of the corresponding degrees of freedom.
- *
- * porder counts 1 for linears, 2 for quadratics, and so on.
  */
-vector2d<int> complete_poly_order(vector2d<int> const &cells, int porder);
+vector2d<int> complete_poly_order(vector2d<int> const &cells, int degree);
 
 /*!
  * \brief Completes the cells to indexes of degrees of freedom
  *
  * Given the active cells and padded cells, the returned list of indexes
  * will hold all indexes of the corresponding degrees of freedom.
- *
- * porder counts 1 for linears, 2 for quadratics, and so on.
  */
 vector2d<int> complete_poly_order(vector2d<int> const &cells,
-                                  indexset const &padded, int porder);
+                                  indexset const &padded, int degree);
 
 } // namespace asgard

src/asgard_indexset_tests.cpp

@@ -125,11 +125,11 @@ TEST_CASE("connectivity full and expanded to dof", "[connectivity]")
 
   // expand the cells by adding the degrees of freedom for quadratic basis
   // i.e., each entry in the sparse matrix is replaced with a 3x3 block
-  int const porder = 2;
-  connect_1d expanded(cells, porder);
-  REQUIRE(expanded.num_rows() == (porder + 1) * 8);
+  int const degree = 2;
+  connect_1d expanded(cells, degree);
+  REQUIRE(expanded.num_rows() == (degree + 1) * 8);
   // there are fewer connection since we removed the self-connection
-  REQUIRE(expanded.num_connections() == 60 * (porder + 1) * (porder + 1));
+  REQUIRE(expanded.num_connections() == 60 * (degree + 1) * (degree + 1));
 
   // compare the connectivity to the 12-th element (first in cell 4)
   REQUIRE(expanded.row_end(12) - expanded.row_begin(12) == 21);