More bc (#767)

*  moved the BC construction to the term_md

Author: mkstoyanov

python/asgardrun.sh

@@ -31,8 +31,6 @@ exename=$1
 
 shift
 
-echo $1
-
 ./$exename $@ -of _plt.h5
 
 @Python_EXECUTABLE@ -m asgard _plt.h5 $plt_opts

src/asgard_boundary_conditions.cpp

@@ -198,181 +198,6 @@ std::array<unscaled_bc_parts<P>, 2> make_unscaled_bc_parts(
   return {left_bc_parts, right_bc_parts};
 }
 
-// template<typename P>
-// bc_manager<P>::bc_manager(PDE<P> const &pde, elements::table const &table,
-//     hierarchy_manipulator<P> const &hier, coefficient_matrices<P> &cmats,
-//     connection_patterns const &conns,
-//     int const start_element, int const stop_element, P const t_init)
-//   : num_dims(pde.num_dims()), degree(hier.degree())
-// {
-//   tools::time_event timing("make unscaled bc");
-//   expect(start_element >= 0);
-//   expect(stop_element < table.size());
-//   expect(stop_element >= start_element);
-//
-//   // count the number of boundary condition vectors that we will need
-//   // the first vector (if present) is left unscaled, i.e., assumed 1
-//   // the others will use corresponding functions from the p-terms
-//   // so, we first count the number of basis that we will need
-//   int num_pt    = 0;
-//
-//   for (auto const &term : pde.get_terms()) {
-//     for (int d : indexof<int>(pde.num_dims())) {
-//       for (auto const &pt : term[d].get_partial_terms()) {
-//         if (not pt.left_bc_zero())
-//         {
-//           if (pt.left_bc_time_func())
-//             num_pt++;
-//           else
-//             num_const = 1;
-//         }
-//         if (not pt.right_bc_zero())
-//         {
-//           if (pt.right_bc_time_func())
-//             num_pt++;
-//           else
-//             num_const = 1;
-//         }
-//       }
-//     }
-//   }
-//
-//   std::vector<dimension<P>> const &dimensions = pde.get_dimensions();
-//   int const degree = hier.degree();
-//
-//   time_funcs.reserve(num_pt);
-//
-//   int64_t const vec_size = (stop_element - start_element + 1) * hier.block_size();
-//   unscaled_bc.resize(vec_size * (num_pt + num_const));
-//
-//   std::vector<P> temp;
-//
-//   P *unc = unscaled_bc.data();
-//   if (num_const > 0) {
-//     unc += vec_size;
-//     temp.resize(vec_size);
-//   }
-//
-//   for (int t : indexof<int>(pde.num_terms()))
-//   {
-//     std::vector<term<P>> const &term_md = pde.get_terms()[t];
-//
-//     for (int d : indexof<int>(pde.num_dims()))
-//     {
-//       dimension<P> const &dim = dimensions[d];
-//       term<P> const &term     = term_md[d];
-//
-//       std::vector<partial_term<P>> const &pterms = term.get_partial_terms();
-//       for (int pt : indexof<int>(pterms.size()))
-//       {
-//         auto const &ptt = pterms[pt];
-//
-//         if (not ptt.left_bc_zero())
-//         {
-//           std::vector<P> trace_bc = compute_left_boundary_condition(
-//               ptt.g_func(), ptt.dv_func(), t_init, dim,
-//               ptt.left_bc_funcs()[d]);
-//
-//           std::vector<std::vector<P>> p_term_left_bcs = generate_partial_bcs(
-//               dimensions, d, ptt.left_bc_funcs(), hier, cmats,
-//               conns, t_init, term_md, term.get_partial_terms(), t, pt, std::move(trace_bc));
-//
-//           if (ptt.left_bc_time_func())
-//             num_pt++;
-//           else
-//             num_const = 1;
-//         }
-//         if (not ptt.right_bc_zero())
-//         {
-//           if (ptt.right_bc_time_func())
-//             num_pt++;
-//           else
-//             num_const = 1;
-//         }
-//       }
-//     }
-//   }
-//
-//
-//
-//   unscaled_bc_parts<P> left_bc_parts;
-//   unscaled_bc_parts<P> right_bc_parts;
-//
-//   term_set<P> const &terms_vec_vec = pde.get_terms();
-//
-//   int const num_terms = static_cast<int>(terms_vec_vec.size());
-//   int const num_dims  = pde.num_dims();
-//
-//
-//   //int const degree = dimensions.front().get_degree();
-//
-//   for (int t : indexof<int>(num_terms))
-//   {
-//     std::vector<term<P>> const &term_md = terms_vec_vec[t];
-//
-//     std::vector<std::vector<std::vector<P>>> left_dim_pvecs;
-//     std::vector<std::vector<std::vector<P>>> right_dim_pvecs;
-//
-//     for (int d : indexof<int>(num_dims))
-//     {
-//       dimension<P> const &dim = dimensions[d];
-//
-//       term<P> const &term = term_md[d];
-//
-//       std::vector<std::vector<P>> left_pvecs;
-//       std::vector<std::vector<P>> right_pvecs;
-//
-//       std::vector<partial_term<P>> const &pterms = term.get_partial_terms();
-//       for (int pt : indexof<int>(pterms.size()))
-//       {
-//         partial_term<P> const &pterm = pterms[pt];
-//
-//         if (not pterm.left_bc_zero())
-//         {
-//           std::vector<P> trace_bc = compute_left_boundary_condition(
-//               pterm.g_func(), pterm.dv_func(), t_init, dim,
-//               pterm.left_bc_funcs()[d]);
-//
-//           std::vector<std::vector<P>> p_term_left_bcs = generate_partial_bcs(
-//               dimensions, d, pterm.left_bc_funcs(), hier, cmats,
-//               conns, t_init, term_md, pterms, t, pt, std::move(trace_bc));
-//
-//           // std::vector<P> combined =
-//           //     combine_dimensions(degree, table, start_element,
-//           //                        stop_element, p_term_left_bcs);
-//           //
-//           // left_pvecs.emplace_back(std::move(combined));
-//         }
-//
-//         if (not pterm.right_bc_zero())
-//         {
-//           std::vector<P> trace_bc = compute_right_boundary_condition(
-//               pterm.g_func(), pterm.dv_func(), t_init, dim,
-//               pterm.right_bc_funcs()[d]);
-//
-//           std::vector<std::vector<P>> p_term_right_bcs = generate_partial_bcs(
-//               dimensions, d, pterm.right_bc_funcs(), hier, cmats,
-//               conns, t_init, term_md, pterms, t, pt, std::move(trace_bc));
-//
-//           // std::vector<P> combined =
-//           //     combine_dimensions(degree, table, start_element,
-//           //                        stop_element, p_term_right_bcs);
-//           //
-//           // right_pvecs.emplace_back(std::move(combined));
-//         }
-//       }
-//
-//       left_dim_pvecs.emplace_back(std::move(left_pvecs));
-//       right_dim_pvecs.emplace_back(std::move(right_pvecs));
-//     }
-//
-//     left_bc_parts.emplace_back(std::move(left_dim_pvecs));
-//     right_bc_parts.emplace_back(std::move(right_dim_pvecs));
-//   }
-//
-//   //return {left_bc_parts, right_bc_parts};
-// }
-
 template<typename P>
 std::vector<P> generate_scaled_bc(unscaled_bc_parts<P> const &left_bc_parts,
                                   unscaled_bc_parts<P> const &right_bc_parts,
@@ -528,8 +353,6 @@ compute_right_boundary_condition(g_func_type<P> g_func, g_func_type<P> dv_func,
 
 /* explicit instantiations */
 #ifdef ASGARD_ENABLE_DOUBLE
-template class bc_manager<double>;
-
 template std::array<unscaled_bc_parts<double>, 2> make_unscaled_bc_parts(
     PDE<double> const &pde, elements::table const &table,
     hierarchy_manipulator<double> const &hier, coefficient_matrices<double> &cmats,
@@ -550,8 +373,6 @@ boundary_conditions::compute_right_boundary_condition(
 #endif
 
 #ifdef ASGARD_ENABLE_FLOAT
-template class bc_manager<float>;
-
 template std::array<unscaled_bc_parts<float>, 2>
 boundary_conditions::make_unscaled_bc_parts(
     PDE<float> const &pde, elements::table const &table,

src/asgard_boundary_conditions.hpp

@@ -8,35 +8,6 @@ template<typename P>
 using unscaled_bc_parts =
     std::vector<std::vector<std::vector<std::vector<P>>>>;
 
-template<typename P>
-class bc_manager {
-public:
-  // bc_manager(PDE<P> const &pde, elements::table const &table,
-  //            hierarchy_manipulator<P> const &hier, coefficient_matrices<P> &cmats,
-  //            connection_patterns const &conns,
-  //            int const start_element, int const stop_element, P const t_init = 0);
-
-private:
-  // struct partial_bc_data
-  // {
-  //   int dim = 0;
-  //   int id = -1;
-  //   scalar_func<P> fleft;
-  //   scalar_func<P> fright;
-  //   std::vector<P> vals;
-  // };
-  //
-  // // no product, maybe add it later at the level of terms_md
-  // // use either fixed coefficients or recompute for each t
-  //
-  // // useful local vars
-  // int num_dims = 0;
-  // int degree   = 0;
-  //
-  // std::vector<P> raw_bc;
-  // std::vector<partial_bc_data> bcs;
-};
-
 // generate all boundary conditions
 // add the new transformer here and the coefficient matrices
 template<typename P>

src/asgard_coefficients_mats.hpp

@@ -21,7 +21,7 @@ void gen_tri_cmat(legendre_basis<P> const &basis, P xleft, P xright, int level,
                   sfixed_func1d<P> const &rhs, P const rhs_const, flux_type flux,
                   boundary_type boundary, rhs_raw_data<P> &rhs_raw, block_tri_matrix<P> &coeff)
 {
-  static_assert(optype != operation_type::mass
+  static_assert(optype != operation_type::volume
                 and optype != operation_type::identity
                 and optype != operation_type::chain,
                 "identity, mass and chain operations yield diagonal matrices, "
@@ -33,19 +33,19 @@ void gen_tri_cmat(legendre_basis<P> const &basis, P xleft, P xright, int level,
                 "cannot use spatially dependant penalty term");
 
   if constexpr (optype == operation_type::grad) {
-    // the grad operation flips the dirichlet and free boundayr conditions
+    // the grad operation flips the fixed and free boundary conditions
     switch (boundary) {
-      case boundary_type::dirichlet:
-        boundary = boundary_type::free;
+      case boundary_type::bothsides:
+        boundary = boundary_type::none;
         break;
-      case boundary_type::free:
-        boundary = boundary_type::dirichlet;
+      case boundary_type::none:
+        boundary = boundary_type::bothsides;
         break;
-      case boundary_type::left_free:
-        boundary = boundary_type::right_free;
+      case boundary_type::right:
+        boundary = boundary_type::left;
         break;
-      case boundary_type::right_free:
-        boundary = boundary_type::left_free;
+      case boundary_type::left:
+        boundary = boundary_type::right;
         break;
       default: // periodic, do nothing since it is symmetric anyway
         break;
@@ -174,8 +174,8 @@ void gen_tri_cmat(legendre_basis<P> const &basis, P xleft, P xright, int level,
       if constexpr (optype == operation_type::penalty)
       {
         switch (boundary) {
-          case boundary_type::dirichlet:
-          case boundary_type::right_free: // dirichelt on the left
+          case boundary_type::bothsides:
+          case boundary_type::left: // dirichelt on the left
             smmat::axpy(nblock, escale * rhs_const, basis.to_left, coeff.diag(0));
             break;
           case boundary_type::periodic:
@@ -195,8 +195,8 @@ void gen_tri_cmat(legendre_basis<P> const &basis, P xleft, P xright, int level,
         }
 
         switch (boundary) {
-          case boundary_type::dirichlet:
-          case boundary_type::left_free: // dirichelt on the right
+          case boundary_type::bothsides:
+          case boundary_type::right: // dirichelt on the right
             smmat::axpy(nblock, escale * rhs_const, basis.to_right, coeff.diag(rmost));
             break;
           case boundary_type::periodic:
@@ -211,8 +211,8 @@ void gen_tri_cmat(legendre_basis<P> const &basis, P xleft, P xright, int level,
       {
         // look at the left-boundary
         switch (boundary) {
-          case boundary_type::free:
-          case boundary_type::left_free: // free on the left
+          case boundary_type::none:
+          case boundary_type::right: // free on the left
             smmat::axpy(nblock, -escale * ((rtype == rhs_type::is_const) ? rhs_const : rhs_vals[0][0]), basis.to_left, coeff.diag(0));
             break;
           case boundary_type::periodic: {
@@ -242,8 +242,8 @@ void gen_tri_cmat(legendre_basis<P> const &basis, P xleft, P xright, int level,
 
         // look at the right-boundary
         switch (boundary) {
-          case boundary_type::free:
-          case boundary_type::right_free: // free on the right
+          case boundary_type::none:
+          case boundary_type::left: // free on the right
             smmat::axpy(nblock, escale * ((rtype == rhs_type::is_const) ? rhs_const : rhs_raw.vals.back()), basis.to_right, coeff.diag(rmost));
             break;
           case boundary_type::periodic: {
@@ -278,7 +278,7 @@ template<typename P, operation_type optype>
 void gen_diag_cmat(legendre_basis<P> const &basis, int level,
                    P const rhs_const, block_diag_matrix<P> &coeff)
 {
-  static_assert(optype == operation_type::mass,
+  static_assert(optype == operation_type::volume,
                 "only mass matrices should be used to create mass terms");
 
   int const num_cells = fm::ipow2(level);
@@ -331,7 +331,7 @@ void gen_diag_cmat(legendre_basis<P> const &basis, P xleft, P xright, int level,
                    block_diag_matrix<P> &coeff)
 {
   ignore(rhs_f);
-  static_assert(optype == operation_type::mass,
+  static_assert(optype == operation_type::volume,
                 "only mass matrices should be used to create mass terms");
 
   int const num_cells = fm::ipow2(level);
@@ -351,7 +351,6 @@ void gen_diag_cmat(legendre_basis<P> const &basis, P xleft, P xright, int level,
 #pragma omp parallel for
     for (int i = 0; i < num_cells; i++) {
       P const l = xleft + i * dx; // left edge of cell i
-      rhs_pnts[i * basis.num_quad] = l;
       for (int k = 0; k < basis.num_quad; k++)
         rhs_pnts[i * basis.num_quad + k] = (0.5 * basis.qp[k] + 0.5) * dx + l;
     }

src/asgard_coefficients_mats_tests.cpp

@@ -129,7 +129,7 @@ TEMPLATE_TEST_CASE("simple div", "[div]", test_precs)
   }
 
   gen_tri_cmat<P, operation_type::div, rhs_type::is_const>(
-      basis, 0, 1, level, nullptr, 1, flux_type::central, boundary_type::free, rhs_raw, mat);
+      basis, 0, 1, level, nullptr, 1, flux_type::central, boundary_type::none, rhs_raw, mat);
 
   std::vector<P> const ref = {0, -4, 4, -4, 0, 4, -4, 0, 4, -4, 0, 4, -4, 0, 4,
                               -4, 0, 4, -4, 0, 4, -4, 4, 0};
@@ -153,7 +153,7 @@ TEMPLATE_TEST_CASE("simple mass", "[mass]", test_precs)
 
   block_diag_matrix<P> mat;
 
-  gen_diag_cmat<P, operation_type::mass>(basis, level, 1, mat);
+  gen_diag_cmat<P, operation_type::volume>(basis, level, 1, mat);
 
   for (int i = 0; i < 8; i++) {
     std::vector<P> ref = {1, 0, 0, 0, 1, 0, 0, 0, 1};
@@ -167,7 +167,7 @@ TEMPLATE_TEST_CASE("simple mass", "[mass]", test_precs)
         fx[i] = -3.5;
     };
 
-  gen_diag_cmat<P, operation_type::mass>(
+  gen_diag_cmat<P, operation_type::volume>(
       basis, 0, 1, level, cc, nullptr, mat);
 
   for (int i = 0; i < 8; i++) {

src/asgard_dimension.hpp

@@ -353,7 +353,7 @@ class separable_func
     std::copy(fdomain.begin(), fdomain.end(), consts_.begin());
   }
   //! set a function that is constant throughout the domain but has a time component
-  separable_func(std::vector<P> fdomain)
+  separable_func(std::vector<P> const &fdomain)
     : ignores_time_(true)
   {
     expect(static_cast<int>(fdomain.size()) <= max_num_dimensions);