From c2fcb0f8b727f31c811f2a71dfbcd1f66ac50754 Mon Sep 17 00:00:00 2001
From: Miroslav Stoyanov <30537612+mkstoyanov@users.noreply.github.com>
Date: Tue, 4 Feb 2025 16:52:31 -0500
Subject: [PATCH] Improved testing (#763)
* fixed some bugs, most notably in restarting
---
src/asgard_block_matrix.cpp | 6 +
src/asgard_discretization.cpp | 24 ++-
src/asgard_io.cpp | 10 +-
src/asgard_io_tests.cpp | 114 +++++++++++++-
src/asgard_pde.cpp | 49 ++++++
src/asgard_pde.hpp | 61 ++++++++
src/asgard_testpdes.hpp | 278 ++++++++++++++++++++++++++++++++++
src/continuity.cpp | 2 +-
src/pde/asgard_pde_base.hpp | 2 +
testing/tests_general.hpp | 4 +-
10 files changed, 539 insertions(+), 11 deletions(-)
create mode 100644 src/asgard_pde.cpp
create mode 100644 src/asgard_testpdes.hpp
diff --git a/src/asgard_block_matrix.cpp b/src/asgard_block_matrix.cpp
index 0ef427486..a50c905e3 100644
--- a/src/asgard_block_matrix.cpp
+++ b/src/asgard_block_matrix.cpp
@@ -214,12 +214,18 @@ void gemm_block_tri(int const n, block_tri_matrix
const &A, block_tri_matrix<
block_tri_matrix
&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));
diff --git a/src/asgard_discretization.cpp b/src/asgard_discretization.cpp
index 1c95c4149..4d2ee6778 100644
--- a/src/asgard_discretization.cpp
+++ b/src/asgard_discretization.cpp
@@ -193,6 +193,10 @@ void discretization_manager::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
void discretization_manager::restart_from_file()
{
#ifdef ASGARD_USE_HIGHFIVE
+ tools::time_event timing_("restart from file");
+
time_data dtime;
h5manager::read(pde2.options().restart_file, high_verbosity(), pde2, sgrid,
dtime, state);
@@ -316,7 +322,19 @@ void discretization_manager::restart_from_file()
terms = term_manager(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::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();
diff --git a/src/asgard_io.cpp b/src/asgard_io.cpp
index 1afe11db8..4000ea8dc 100644
--- a/src/asgard_io.cpp
+++ b/src/asgard_io.cpp
@@ -411,9 +411,7 @@ void h5manager::read(std::string const &filename, bool silent, PDEv2
&pde,
}
} else if (stop >= 0) { // overriding the stop time, dt is not set
if (n >= 0) {
- P const fdt = H5Easy::load
(file, "dtime_dt");
- dtime = time_data
(sm, typename time_data
::input_dt{fdt}, n);
- fstop += dtime.stop_time_;
+ dtime = time_data
(sm, n, typename time_data
::input_stop_time{stop - curr_time});
} else {
P const fdt = H5Easy::load
(file, "dtime_dt");
dtime = time_data
(sm,
@@ -455,12 +453,12 @@ void h5manager
::read(std::string const &filename, bool silent, PDEv2
&pde,
grid.iset_.indexes_ = H5Easy::load>(file, "grid_indexes");
- grid.dsort_ = dimension_sort(grid.iset_);
-
if (grid.iset_.indexes_.size() != static_cast(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::read(std::string const &filename, bool silent, PDEv2
&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(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;
}
}
diff --git a/src/asgard_io_tests.cpp b/src/asgard_io_tests.cpp
index c1c52263d..ce20e7334 100644
--- a/src/asgard_io_tests.cpp
+++ b/src/asgard_io_tests.cpp
@@ -1,6 +1,6 @@
#include "tests_general.hpp"
-// workaround for missing include issue with highfive
+// reintroduce private headers
#include
#include
#include
@@ -96,7 +96,7 @@ TEMPLATE_TEST_CASE("save/restart logic", "[io]", test_precs)
std::string const filename = (std::is_same_v)
? "_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 disc(make_testpde(2, options));
+ REQUIRE(get_qoi_indicator(disc) < 1.E-2);
+ advance_time(disc);
+ REQUIRE(get_qoi_indicator(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 rdisc(make_testpde(2, ropts));
+
+ REQUIRE(std::abs(get_qoi_indicator(disc) - get_qoi_indicator(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 reff(make_testpde(2, options));
+ advance_time(reff);
+
+ REQUIRE(std::abs(reff.time_params().time() - 0.08) < 1.E-8);
+
+ REQUIRE(std::abs(get_qoi_indicator(reff) - get_qoi_indicator(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 disc(make_testpde(2, options));
+ REQUIRE(get_qoi_indicator(disc) < 1.E-2);
+ advance_time(disc, 4);
+ REQUIRE(get_qoi_indicator(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 rdisc(make_testpde(2, ropts));
+
+ REQUIRE(rdisc.get_sgrid().num_indexes() == disc.get_sgrid().num_indexes());
+ REQUIRE(std::abs(get_qoi_indicator(disc) - get_qoi_indicator(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 reff(make_testpde(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(reff) - get_qoi_indicator(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 disc(make_testpde(2, options));
+ double const ienergy = get_qoi_indicator(disc);
+ advance_time(disc, 4);
+ double tol = (std::is_same_v) ? 1.E-8 : 1.E-5;
+ REQUIRE(std::abs(ienergy - get_qoi_indicator(disc)) < tol);
+
+ disc.save_final_snapshot();
+ REQUIRE(std::filesystem::exists("_asg_testfile.h5"));
+
+ auto ropts = make_opts("-restart _asg_testfile.h5");
+ discretization_manager rdisc(make_testpde(2, ropts));
+
+ REQUIRE(rdisc.get_sgrid().num_indexes() == disc.get_sgrid().num_indexes());
+ REQUIRE(std::abs(get_qoi_indicator(disc) - get_qoi_indicator(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(rdisc) - get_qoi_indicator(disc)) < 1.E-8);
+ }
+}
diff --git a/src/asgard_pde.cpp b/src/asgard_pde.cpp
new file mode 100644
index 000000000..d9b417e15
--- /dev/null
+++ b/src/asgard_pde.cpp
@@ -0,0 +1,49 @@
+#include "asgard_pde.hpp"
+
+#include "device/asgard_kronmult_common.hpp"
+
+namespace asgard
+{
+
+template
+void builtin_v::positive(std::vector
const &x, std::vector
&y)
+{
+#pragma omp parallel for
+ for (size_t i = 0; i < x.size(); i++)
+ y[i] = std::max(P{0}, x[i]);
+}
+template
+void builtin_v::negative(std::vector
const &x, std::vector
&y)
+{
+#pragma omp parallel for
+ for (size_t i = 0; i < x.size(); i++)
+ y[i] = std::min(P{0}, x[i]);
+}
+
+template
+void builtin_v::sin(std::vector
const &x, std::vector
&y) {
+ ASGARD_OMP_PARFOR_SIMD
+ for (size_t i = 0; i < x.size(); i++)
+ y[i] = std::sin(x[i]);
+}
+template
+void builtin_v::cos(std::vector
const &x, std::vector
&y) {
+ ASGARD_OMP_PARFOR_SIMD
+ for (size_t i = 0; i < x.size(); i++)
+ y[i] = std::cos(x[i]);
+}
+template
+void builtin_v::dcos(std::vector
const &x, std::vector
&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;
+#endif
+
+#ifdef ASGARD_ENABLE_FLOAT
+ template struct builtin_v;
+#endif
+} // namespace asgard
diff --git a/src/asgard_pde.hpp b/src/asgard_pde.hpp
index 43220ebe1..28feca5e7 100644
--- a/src/asgard_pde.hpp
+++ b/src/asgard_pde.hpp
@@ -112,4 +112,65 @@ std::unique_ptr> make_PDE(std::string const &opts)
return make_PDE(make_opts(opts));
}
+/*!
+ * \internal
+ * \brief Wraps around commonly used vector functions
+ *
+ * \endinternal
+ */
+template
+struct builtin_v {
+ //! y is equal to x with all negative values replaced by zero
+ static void positive(std::vector const &x, std::vector
&y);
+ //! y is equal to x with all positive values replaced by zero
+ static void negative(std::vector
const &x, std::vector
&y);
+
+ //! vector version of std::sin()
+ static void sin(std::vector
const &x, std::vector
&y);
+ //! vector version of std::cos()
+ static void cos(std::vector
const &x, std::vector
&y);
+ //! vector version of derivative of std::cos(), i.e., -std::sin()
+ static void dcos(std::vector
const &x, std::vector
&y);
+
+ //! wrapper for sin
+};
+
+/*!
+ * \internal
+ * \brief Wraps around commonly used functions, with time parameter
+ *
+ * \endinternal
+ */
+template
+struct builtin_t {
+ //! overloads with dummy time parameter
+ static void sin(std::vector const &x, P, std::vector
&y) {
+ builtin_v
::sin(x, y);
+ }
+ //! overloads with dummy time parameter
+ static void cos(std::vector
const &x, P, std::vector
&y) {
+ builtin_v
::cos(x, y);
+ }
+ //! overloads with dummy time parameter
+ static void dcos(std::vector
const &x, P, std::vector
&y) {
+ builtin_v
::dcos(x, y);
+ }
+};
+
+/*!
+ * \internal
+ * \brief Wraps around commonly used functions, scalar variant
+ *
+ * \endinternal
+ */
+template
+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
diff --git a/src/asgard_testpdes.hpp b/src/asgard_testpdes.hpp
new file mode 100644
index 000000000..b3b0cd11f
--- /dev/null
+++ b/src/asgard_testpdes.hpp
@@ -0,0 +1,278 @@
+#pragma once
+
+#include "asgard.hpp"
+
+/*!
+ * \internal
+ * \file asgard_testpde.cpp
+ * \brief Simple PDEs used for testing, not included in the library
+ * \author The ASGarD Team
+ * \ingroup asgard_testing
+ *
+ * \endinternal
+ */
+
+namespace asgard
+{
+
+#ifndef __ASGARD_DOXYGEN_SKIP_INTERNAL
+
+/*!
+ * \internal
+ * \ingroup asgard_testing
+ * \brief Simpler version of the continuity-md example using cos-waves in +/- 1.5 PI
+ *
+ * \endinternal
+ */
+struct pde_contcos {};
+/*!
+ * \internal
+ * \ingroup asgard_testing
+ * \brief Most simple PDE using Poisson solver, but no analytic solution
+ *
+ * \endinternal
+ */
+struct pde_twostream {};
+
+#endif
+
+/*!
+ * \internal
+ * \ingroup asgard_testing
+ * \brief Creates a simple test pde for the given dimensions and options
+ *
+ * \endinternal
+ */
+template
+PDEv2 make_testpde(int num_dims, prog_opts options) {
+
+ if constexpr (std::is_same_v)
+ {
+ options.title = "Simplified Continuity " + std::to_string(num_dims) + "D";
+
+ std::vector> ranges(num_dims, {-P{3} * PI / P{2}, P{3} * PI / P{2}});
+
+ pde_domain domain(ranges);
+
+ int const max_level = options.max_level();
+
+ P const dx = domain.min_cell_size(max_level);
+
+ options.default_dt = 0.5 * 0.1 * dx;
+
+ PDEv2
pde(std::move(options), std::move(domain));
+
+ term_1d
div = term_div
(1, flux_type::upwind, boundary_type::periodic);
+
+ // the multi-dimensional divergence, initially set to identity in md
+ std::vector> ops(num_dims);
+ for (int d = 0; d < num_dims; d++)
+ {
+ ops[d] = div;
+ pde += term_md(ops);
+ ops[d] = term_identity{};
+ }
+
+ // put the time-parameter inside one of the cos-functions
+ // tests the non-separable in time capabilities
+ auto cos_1t = [](std::vector
const &x, P t, std::vector
&fx) ->
+ void {
+ for (size_t i = 0; i < x.size(); i++)
+ fx[i] = std::cos(t) * std::cos(x[i]);
+ };
+
+ std::vector> func_md(num_dims, builtin_t::cos);
+
+ func_md[0] = cos_1t;
+ pde.add_initial(func_md);
+ func_md[0] = builtin_t
::cos;
+
+ pde.add_source({func_md, builtin_s
::dcos}); // derivative in time
+
+ func_md[0] = builtin_t
::dcos;
+ pde.add_source({func_md, builtin_s
::cos});
+ func_md[0] = cos_1t;
+
+ // compute the spacial derivatives
+ for (int d = 1; d < num_dims; d++)
+ {
+ func_md[d] = builtin_t
::dcos;
+ pde.add_source(func_md);
+ func_md[d] = builtin_t
::cos;
+ }
+
+ return pde;
+ }
+ else if constexpr (std::is_same_v)
+ {
+ options.title = "Test Two Stream Instability";
+
+ // the domain has one position and one velocity dimension: 1x1v
+ pde_domain domain(position_dims{1}, velocity_dims{1},
+ {{-2 * PI, 2 * PI}, {-2 * PI, 2 * PI}});
+
+ options.default_degree = 2;
+
+ // the CFL is more complicated
+ int const k = options.degree.value_or(options.default_degree.value());
+ int const n = (1 << options.max_level());
+ options.default_dt = 3.0 / (2 * (2 * k + 1) * n);
+
+ options.default_step_method = time_advance::method::rk2;
+
+ PDEv2
pde(options, domain);
+
+ pde += term_md
(std::vector>{
+ term_div(1, flux_type::upwind, boundary_type::periodic),
+ term_mass
(builtin_v
::positive)
+ });
+
+ pde += term_md
(std::vector>{
+ term_div(1, flux_type::downwind, boundary_type::periodic),
+ term_mass
(builtin_v
::negative),
+ });
+
+ pde += term_md
(std::vector>{
+ mass_electric(builtin_v
::positive),
+ term_div
(1, flux_type::upwind, boundary_type::dirichlet)
+ });
+
+ pde += term_md
(std::vector>{
+ mass_electric(builtin_v
::negative),
+ term_div
(1, flux_type::downwind, boundary_type::dirichlet)
+ });
+
+ // initial conditions in x and v
+ auto ic_x = [](std::vector
const &x, P /* time */, std::vector
&fx) ->
+ void {
+ for (size_t i = 0; i < x.size(); i++)
+ fx[i] = 1.0 - 0.5 * std::cos(0.5 * x[i]);
+ };
+
+ auto ic_v = [](std::vector
const &v, P /* time */, std::vector
&fv) ->
+ void {
+ P const c = P{1} / std::sqrt(PI);
+
+ for (size_t i = 0; i < v.size(); i++)
+ fv[i] = c * v[i] * v[i] * std::exp(-v[i] * v[i]);
+ };
+
+ pde.add_initial(asgard::separable_func
({ic_x, ic_v}));
+
+ return pde;
+
+ } else {
+ rassert(false, "Incorrect pde type for make_testpde");
+ }
+}
+
+/*!
+ * \internal
+ * \ingroup asgard_testing
+ * \brief Returns an indicator of the "health" of the PDE
+ *
+ * In cases when the PDE has a known analytic solution, this will simply return
+ * the L^2 error between the current state and the known solution.
+ *
+ * In other cases, the indicator can be different, e.g., measuring energy conservation.
+ *
+ * \endinternal
+ */
+template
+double get_qoi_indicator(asgard::discretization_manager const &disc) {
+
+ if constexpr (std::is_same_v)
+ {
+ // there is no analytic solution, using the sum of particle potential and kinetic
+ // energy as the indicator, it is not zero but must be near constant
+ int const num_moms = 3;
+ int const pdof = disc.degree() + 1;
+ moments1d moms(num_moms, pdof - 1, disc.get_pde2().max_level(),
+ disc.get_pde2().domain());
+ std::vector mom_vec;
+
+ moms.project_moments(disc.get_sgrid(), disc.current_state(), mom_vec);
+
+ disc.do_poisson_update(disc.current_state()); // update the electric field
+
+ auto const &efield = disc.get_terms().cdata.electric_field;
+
+ int const level0 = disc.get_sgrid().current_level(0);
+ int const num_cell = fm::ipow2(level0);
+ double const dx = disc.get_pde2().domain().length(0) / num_cell;
+
+ double Ep = 0;
+ for (auto e : efield)
+ Ep += e * e;
+ Ep *= dx;
+
+ span2d
moments(num_moms * pdof, num_cell, mom_vec.data());
+
+ double Ek = 0;
+ for (int j : iindexof(num_cell))
+ Ek += moments[j][2 * pdof]; // integrating the third moment
+ Ek *= std::sqrt(disc.get_pde2().domain().length(0));
+
+ return Ep + Ek;
+ }
+
+ int const num_dims = disc.num_dims();
+
+ std::vector
const eref = disc.project_function(disc.get_pde2().ic_sep());
+
+ auto [space1d, timev] = [&]() -> std::array {
+ if constexpr (std::is_same_v) {
+ return {1.5 * PI, std::cos(disc.time_params().time())};
+ } else { // no analytic solution, code will be intercepted above
+ return {0, 0};
+ }
+ }();
+
+ double const enorm = asgard::fm::powi(space1d, num_dims) * timev * timev;
+
+ std::vector const &state = disc.current_state();
+ assert(eref.size() == state.size());
+
+ double nself = 0;
+ double ndiff = 0;
+ for (size_t i = 0; i < state.size(); i++)
+ {
+ double const e = eref[i] - state[i];
+ ndiff += e * e;
+ double const r = eref[i];
+ nself += r * r;
+ }
+
+ if (enorm < 1)
+ return std::sqrt(ndiff + enorm - nself);
+ else
+ return std::sqrt((ndiff + enorm - nself) / enorm);
+}
+
+/*!
+ * \internal
+ * \ingroup asgard_testing
+ * \brief Using the given PDE type and opts, integrate step-by-step and return max-L^2 error
+ *
+ * \endinternal
+ */
+template
+double get_time_error(int num_dims, std::string const &opts) {
+
+ auto options = make_opts(opts);
+
+ discretization_manager disc(make_testpde(num_dims, options));
+
+ double max_err = 0;
+
+ while (disc.time_params().num_remain() > 0)
+ {
+ advance_time(disc, 1);
+
+ max_err = std::max(max_err, get_qoi_indicator(disc));
+ }
+
+ return max_err;
+}
+
+} // namespace asgard
diff --git a/src/continuity.cpp b/src/continuity.cpp
index f28a31483..5f27aedd7 100644
--- a/src/continuity.cpp
+++ b/src/continuity.cpp
@@ -485,7 +485,7 @@ void self_test() {
dolongtest(0.02, 2, "-l 5 -t 10");
// adaptivity is tricky near the time-period when the solution vanishes
- dolongtest(0.03, 2, "-l 4 -m 8 -t 10 -a 1.E-2");
+ dolongtest(0.03, 2, "-l 4 -m 8 -t 10 -a 5.E-3");
// different explicit time-stepping
dotest(0.05, 2, "-s rk2 -l 5 -n 20");
diff --git a/src/pde/asgard_pde_base.hpp b/src/pde/asgard_pde_base.hpp
index b7eb3a3aa..535324488 100644
--- a/src/pde/asgard_pde_base.hpp
+++ b/src/pde/asgard_pde_base.hpp
@@ -2008,6 +2008,8 @@ class time_data
if (num_remain_ == 0)
num_remain_ = 1;
// readjust dt to minimize rounding error
+ if (dt_ * num_remain_ < stop_time_)
+ num_remain_ += 1;
dt_ = stop_time_ / static_cast(num_remain_);
}
//! specify number of steps and final time
diff --git a/testing/tests_general.hpp b/testing/tests_general.hpp
index e123b06f0..777e8be15 100644
--- a/testing/tests_general.hpp
+++ b/testing/tests_general.hpp
@@ -7,10 +7,12 @@
#pragma once
#include "asgard.hpp"
+#include "asgard_testpdes.hpp"
+
#include
/*!
- * \defgroup AsgardTesting Miscellaneous testing utilities
+ * \defgroup asgard_testing Miscellaneous testing utilities
*
* Helper functions to facilitate testing.
*/