Improved testing (#763)

* fixed some bugs, most notably in restarting

src/asgard_block_matrix.cpp

@@ -214,12 +214,18 @@ void gemm_block_tri(int const n, block_tri_matrix<P> const &A, block_tri_matrix<
                     block_tri_matrix<P> &C)
 {
   int const M = A.nrows();
+  expect(M >= 1);
   expect(A.nblock() == B.nblock());
   expect(A.nblock() == C.nblock());
   expect(A.nblock() == n * n);
   expect(B.nrows() == M);
   expect(C.nrows() == M);
 
+  if (M == 1) {
+    smmat::gemm<0>(n, A.diag(0), B.diag(0), C.diag(0));
+    return;
+  }
+
   smmat::gemm<0>(n, A.diag(0), B.lower(0), C.lower(0));
   smmat::gemm<1>(n, A.lower(0), B.diag(M - 1), C.lower(0));
 

src/asgard_discretization.cpp

@@ -193,6 +193,10 @@ void discretization_manager<precision>::start_cold()
     std::cout << '\n';
 
     std::cout << sgrid;
+    if (options.adapt_threshold)
+      std::cout << "  adaptive tolerance: " << options.adapt_threshold.value() << '\n';
+    else
+      std::cout << "  non-adaptive\n";
   }
 
   { // setting up the time-step approach
@@ -299,6 +303,8 @@ template<typename precision>
 void discretization_manager<precision>::restart_from_file()
 {
 #ifdef ASGARD_USE_HIGHFIVE
+  tools::time_event timing_("restart from file");
+
   time_data<precision> dtime;
   h5manager<precision>::read(pde2.options().restart_file, high_verbosity(), pde2, sgrid,
                              dtime, state);
@@ -316,7 +322,19 @@ void discretization_manager<precision>::restart_from_file()
   terms = term_manager<precision>(pde2);
 
   terms.prapare_workspace(sgrid);
-  // initialize the moments here, we already have the the state
+
+  // moments, identical to the start_cold() case
+  mom_deps deps = terms.find_deps();
+  if (deps.num_moments > 0) {
+    moms1d = moments1d(deps.num_moments, degree_, pde2.max_level(), pde2.domain());
+    if (deps.poisson) {
+      poisson = solvers::poisson(degree_, pde2.domain().xleft(0), pde2.domain().xright(0),
+                                 sgrid.current_level(0));
+
+      terms.cdata.electric_field.resize(fm::ipow2(sgrid.current_level(0)));
+    }
+  }
+
   if (stepper.needed_precon() == preconditioner_opts::adi) {
     precision const substep
         = (options.step_method.value() == time_advance::method::cn) ? 0.5 : 1;
@@ -331,6 +349,10 @@ void discretization_manager<precision>::restart_from_file()
     if (not options.subtitle.empty())
       std::cout << "subtitle: " << options.subtitle << '\n';
     std::cout << sgrid;
+    if (options.adapt_threshold)
+      std::cout << "  adaptive tolerance: " << options.adapt_threshold.value() << '\n';
+    else
+      std::cout << "  non-adaptive\n";
     std::cout << stepper;
     if (high_verbosity())
       progress_report();

src/asgard_io.cpp

@@ -411,9 +411,7 @@ void h5manager<P>::read(std::string const &filename, bool silent, PDEv2<P> &pde,
       }
     } else if (stop >= 0) { // overriding the stop time, dt is not set
       if (n >= 0) {
-        P const fdt = H5Easy::load<P>(file, "dtime_dt");
-        dtime = time_data<P>(sm, typename time_data<P>::input_dt{fdt}, n);
-        fstop += dtime.stop_time_;
+        dtime = time_data<P>(sm, n, typename time_data<P>::input_stop_time{stop - curr_time});
       } else {
         P const fdt = H5Easy::load<P>(file, "dtime_dt");
         dtime = time_data<P>(sm,
@@ -455,12 +453,12 @@ void h5manager<P>::read(std::string const &filename, bool silent, PDEv2<P> &pde,
 
     grid.iset_.indexes_ = H5Easy::load<std::vector<int>>(file, "grid_indexes");
 
-    grid.dsort_ = dimension_sort(grid.iset_);
-
     if (grid.iset_.indexes_.size() != static_cast<size_t>(num_dims * num_indexes))
       throw std::runtime_error("file corruption detected: wrong number of sparse grid "
                                "indexes found in the file");
 
+    grid.dsort_ = dimension_sort(grid.iset_);
+
     // checking the max levels, we can reset the max level for the simulation
     // first we follow the same logic for specifying either all dims or a single int
     // then we do not allow the max level to be reduced below the current level
@@ -494,7 +492,7 @@ void h5manager<P>::read(std::string const &filename, bool silent, PDEv2<P> &pde,
     if (not pde.options_.adapt_threshold) { // no adapt is loaded
       if (not pde.options_.set_no_adapt) { // adaptivity wasn't explicitly canceled
         double const adapt = H5Easy::load<double>(file, "grid_adapt_threshold");
-        if (adapt > 0) // if negative, then adaptivity was never set to begind with
+        if (adapt > 0) // if negative, then adaptivity was never set to begin with
           pde.options_.adapt_threshold = adapt;
       }
     }

src/asgard_io_tests.cpp

@@ -1,6 +1,6 @@
 #include "tests_general.hpp"
 
-// workaround for missing include issue with highfive
+// reintroduce private headers
 #include <highfive/H5Easy.hpp>
 #include <highfive/H5DataType.hpp>
 #include <highfive/H5DataSpace.hpp>
@@ -96,7 +96,7 @@ TEMPLATE_TEST_CASE("save/restart logic", "[io]", test_precs)
   std::string const filename = (std::is_same_v<TestType, double>)
       ? "_asgard_dsave_test1.h5" : "_asgard_fsave_test1.h5";
 
-  SECTION("simple resrat")
+  SECTION("simple restart")
   {
     int const num_dims = 4;
     std::string const title    = "basic io test";
@@ -228,3 +228,113 @@ TEMPLATE_TEST_CASE("save/restart logic", "[io]", test_precs)
                         "cannot reset the final time to an instance before the current time");
   }
 }
+
+TEMPLATE_TEST_CASE("save/restart logic (longer)", "[io]", test_precs)
+{
+  using P = TestType;
+  SECTION("basic restart")
+  {
+    using pde = pde_contcos;
+
+    // 1. make a pde and set 4 time-steps, advance in time and save the state
+    //    - check initial and final error, and file-existing
+    // 2. restart and set 4 additional time-steps, advance in time
+    //    - make sure restarted matches the saved and new end-time is set
+    // 3. compare against a one-shot integration using 4 steps
+
+    auto options = make_opts("-l 5 -d 2 -dt 0.01 -n 4 -of _asg_testfile.h5");
+    discretization_manager<P> disc(make_testpde<pde, P>(2, options));
+    REQUIRE(get_qoi_indicator<pde, P>(disc) < 1.E-2);
+    advance_time(disc);
+    REQUIRE(get_qoi_indicator<pde, P>(disc) < 1.E-2);
+
+    disc.save_final_snapshot();
+    REQUIRE(std::filesystem::exists("_asg_testfile.h5"));
+
+    auto ropts = make_opts("-n 4 -dt 0.01 -restart _asg_testfile.h5");
+    discretization_manager<P> rdisc(make_testpde<pde, P>(2, ropts));
+
+    REQUIRE(std::abs(get_qoi_indicator<pde, P>(disc) - get_qoi_indicator<pde, P>(rdisc)) < 1.E-10);
+
+    REQUIRE(std::abs(rdisc.time_params().stop_time() - 0.08) < 2.E-9); // updated the stop time
+    advance_time(rdisc);
+
+    REQUIRE(std::abs(rdisc.time_params().time() - 0.08) < 1.E-8);
+
+    options = make_opts("-l 5 -d 2 -dt 0.01 -n 8");
+    discretization_manager<P> reff(make_testpde<pde, P>(2, options));
+    advance_time(reff);
+
+    REQUIRE(std::abs(reff.time_params().time() - 0.08) < 1.E-8);
+
+    REQUIRE(std::abs(get_qoi_indicator<pde, P>(reff) - get_qoi_indicator<pde, P>(rdisc)) < 1.E-10);
+  }
+
+  SECTION("adaptive restart")
+  {
+    using pde = pde_contcos;
+
+    // similar to above but uses adaptivity
+
+    auto options = make_opts("-l 8 -d 2 -dt 0.01 -n 8 -a 1.E-2 -of _asg_testfile.h5");
+    discretization_manager<P> disc(make_testpde<pde, P>(2, options));
+    REQUIRE(get_qoi_indicator<pde, P>(disc) < 1.E-2);
+    advance_time(disc, 4);
+    REQUIRE(get_qoi_indicator<pde, P>(disc) < 1.E-2);
+
+    disc.save_final_snapshot();
+    REQUIRE(std::filesystem::exists("_asg_testfile.h5"));
+
+    auto ropts = make_opts("-restart _asg_testfile.h5");
+    discretization_manager<P> rdisc(make_testpde<pde, P>(2, ropts));
+
+    REQUIRE(rdisc.get_sgrid().num_indexes() == disc.get_sgrid().num_indexes());
+    REQUIRE(std::abs(get_qoi_indicator<pde, P>(disc) - get_qoi_indicator<pde, P>(rdisc)) < 1.E-10);
+
+    REQUIRE(std::abs(rdisc.time_params().stop_time() - 0.08) < 2.E-9); // updated the stop time
+    advance_time(rdisc);
+
+    REQUIRE(std::abs(rdisc.time_params().time() - 0.08) < 1.E-8);
+
+    options = make_opts("-l 8 -d 2 -dt 0.01 -n 8 -a 1.E-2");
+    discretization_manager<P> reff(make_testpde<pde, P>(2, options));
+    advance_time(reff);
+
+    REQUIRE(rdisc.get_sgrid().num_indexes() == reff.get_sgrid().num_indexes());
+
+    REQUIRE(std::abs(reff.time_params().time() - 0.08) < 1.E-8);
+
+    REQUIRE(std::abs(get_qoi_indicator<pde, P>(reff) - get_qoi_indicator<pde, P>(rdisc)) < 1.E-10);
+  }
+
+  SECTION("moments restart")
+  {
+    using pde = pde_twostream;
+
+    // similar to above but uses adaptivity
+
+    auto options = make_opts("-l 7 -d 2 -dt 1.953125E-3 -n 8 -a 1.E-6 -of _asg_testfile.h5");
+    discretization_manager<P> disc(make_testpde<pde, P>(2, options));
+    double const ienergy = get_qoi_indicator<pde, P>(disc);
+    advance_time(disc, 4);
+    double tol = (std::is_same_v<P, double>) ? 1.E-8 : 1.E-5;
+    REQUIRE(std::abs(ienergy - get_qoi_indicator<pde, P>(disc)) < tol);
+
+    disc.save_final_snapshot();
+    REQUIRE(std::filesystem::exists("_asg_testfile.h5"));
+
+    auto ropts = make_opts("-restart _asg_testfile.h5");
+    discretization_manager<P> rdisc(make_testpde<pde, P>(2, ropts));
+
+    REQUIRE(rdisc.get_sgrid().num_indexes() == disc.get_sgrid().num_indexes());
+    REQUIRE(std::abs(get_qoi_indicator<pde, P>(disc) - get_qoi_indicator<pde, P>(rdisc)) < 1.E-10);
+
+    REQUIRE(std::abs(rdisc.time_params().stop_time() - 1.5625E-2) < 1.E-10); // updated the stop time
+    advance_time(rdisc);
+
+    REQUIRE(std::abs(rdisc.time_params().time() - 1.5625E-2) < 1.E-10);
+
+    advance_time(disc);
+    REQUIRE(std::abs(get_qoi_indicator<pde, P>(rdisc) - get_qoi_indicator<pde, P>(disc)) < 1.E-8);
+  }
+}

src/asgard_pde.cpp

@@ -0,0 +1,49 @@
+#include "asgard_pde.hpp"
+
+#include "device/asgard_kronmult_common.hpp"
+
+namespace asgard
+{
+
+template<typename P>
+void builtin_v<P>::positive(std::vector<P> const &x, std::vector<P> &y)
+{
+#pragma omp parallel for
+  for (size_t i = 0; i < x.size(); i++)
+    y[i] = std::max(P{0}, x[i]);
+}
+template<typename P>
+void builtin_v<P>::negative(std::vector<P> const &x, std::vector<P> &y)
+{
+#pragma omp parallel for
+  for (size_t i = 0; i < x.size(); i++)
+    y[i] = std::min(P{0}, x[i]);
+}
+
+template<typename P>
+void builtin_v<P>::sin(std::vector<P> const &x, std::vector<P> &y) {
+  ASGARD_OMP_PARFOR_SIMD
+  for (size_t i = 0; i < x.size(); i++)
+    y[i] = std::sin(x[i]);
+}
+template<typename P>
+void builtin_v<P>::cos(std::vector<P> const &x, std::vector<P> &y) {
+  ASGARD_OMP_PARFOR_SIMD
+  for (size_t i = 0; i < x.size(); i++)
+    y[i] = std::cos(x[i]);
+}
+template<typename P>
+void builtin_v<P>::dcos(std::vector<P> const &x, std::vector<P> &y) {
+  ASGARD_OMP_PARFOR_SIMD
+  for (size_t i = 0; i < x.size(); i++)
+    y[i] = -std::sin(x[i]);
+}
+
+#ifdef ASGARD_ENABLE_DOUBLE
+  template struct builtin_v<double>;
+#endif
+
+#ifdef ASGARD_ENABLE_FLOAT
+  template struct builtin_v<float>;
+#endif
+} // namespace asgard

src/asgard_pde.hpp

@@ -112,4 +112,65 @@ std::unique_ptr<PDE<P>> make_PDE(std::string const &opts)
   return make_PDE<P>(make_opts(opts));
 }
 
+/*!
+ * \internal
+ * \brief Wraps around commonly used vector functions
+ *
+ * \endinternal
+ */
+template<typename P>
+struct builtin_v {
+  //! y is equal to x with all negative values replaced by zero
+  static void positive(std::vector<P> const &x, std::vector<P> &y);
+  //! y is equal to x with all positive values replaced by zero
+  static void negative(std::vector<P> const &x, std::vector<P> &y);
+
+  //! vector version of std::sin()
+  static void sin(std::vector<P> const &x, std::vector<P> &y);
+  //! vector version of std::cos()
+  static void cos(std::vector<P> const &x, std::vector<P> &y);
+  //! vector version of derivative of std::cos(), i.e., -std::sin()
+  static void dcos(std::vector<P> const &x, std::vector<P> &y);
+
+  //! wrapper for sin
+};
+
+/*!
+ * \internal
+ * \brief Wraps around commonly used functions, with time parameter
+ *
+ * \endinternal
+ */
+template<typename P>
+struct builtin_t {
+  //! overloads with dummy time parameter
+  static void sin(std::vector<P> const &x, P, std::vector<P> &y) {
+    builtin_v<P>::sin(x, y);
+  }
+  //! overloads with dummy time parameter
+  static void cos(std::vector<P> const &x, P, std::vector<P> &y) {
+    builtin_v<P>::cos(x, y);
+  }
+  //! overloads with dummy time parameter
+  static void dcos(std::vector<P> const &x, P, std::vector<P> &y) {
+    builtin_v<P>::dcos(x, y);
+  }
+};
+
+/*!
+ * \internal
+ * \brief Wraps around commonly used functions, scalar variant
+ *
+ * \endinternal
+ */
+template<typename P>
+struct builtin_s {
+  //! std::sin(x)
+  static P sin(P x) { return std::sin(x); }
+  //! std::sin(x)
+  static P cos(P x) { return std::cos(x); }
+  //! d/dx std::cos(x) = -std::sin(x)
+  static P dcos(P x) { return -std::sin(x); }
+};
+
 } // namespace asgard