|
| 1 | +/* |
| 2 | + * Tests for the scalar example solving the test equation |
| 3 | + */ |
| 4 | +#include <cmath> |
| 5 | + |
| 6 | +#include <gtest/gtest.h> |
| 7 | +#include <gmock/gmock.h> |
| 8 | + |
| 9 | +using namespace ::testing; |
| 10 | + |
| 11 | +#include <pfasst/quadrature.hpp> |
| 12 | + |
| 13 | +#define PFASST_UNIT_TESTING |
| 14 | +#include "../examples/scalar/scalar_sdc.cpp" |
| 15 | +#undef PFASST_UNIT_TESTING |
| 16 | + |
| 17 | +/* |
| 18 | + * parameterized test fixture with number of nodes as parameter |
| 19 | + */ |
| 20 | +class ConvergenceTest |
| 21 | + : public TestWithParam<tuple<size_t, pfasst::QuadratureType>> |
| 22 | +{ |
| 23 | + protected: |
| 24 | + size_t nnodes; // parameter |
| 25 | + |
| 26 | + complex<double> lambda; |
| 27 | + double Tend; |
| 28 | + vector<size_t> nsteps; |
| 29 | + size_t nsteps_l; |
| 30 | + vector<double> err; |
| 31 | + vector<double> convrate; |
| 32 | + double dt; |
| 33 | + size_t niters; |
| 34 | + pfasst::QuadratureType nodetype; |
| 35 | + size_t nnodes_in_call; |
| 36 | + |
| 37 | + // set parameters base on node type |
| 38 | + void set_parameters() |
| 39 | + { |
| 40 | + switch (this->nodetype) |
| 41 | + { |
| 42 | + case pfasst::QuadratureType::GaussLobatto: |
| 43 | + this->niters = 2 * this->nnodes - 2; |
| 44 | + this->Tend = 4.0; |
| 45 | + this->lambda = complex<double>(-1.0, 1.0); |
| 46 | + this->nsteps = { 2, 5, 10, 15, 20 }; |
| 47 | + this->nnodes_in_call = this->nnodes; |
| 48 | + break; |
| 49 | + |
| 50 | + case pfasst::QuadratureType::GaussLegendre: |
| 51 | + this->niters = 2 * this->nnodes; |
| 52 | + this->Tend = 6.0; |
| 53 | + this->lambda = complex<double>(-1.0, 2.0); |
| 54 | + this->nsteps = { 2, 4, 6, 8, 10 }; |
| 55 | + this->nnodes_in_call = this->nnodes + 2; |
| 56 | + break; |
| 57 | + |
| 58 | + case pfasst::QuadratureType::GaussRadau: |
| 59 | + this->niters = 2 * this->nnodes - 1; |
| 60 | + this->Tend = 5.0; |
| 61 | + this->lambda = complex<double>(-1.0, 2.0); |
| 62 | + this->nsteps = { 4, 6, 8, 10, 12 }; |
| 63 | + this->nnodes_in_call = this->nnodes + 1; |
| 64 | + break; |
| 65 | + |
| 66 | + // NOTE: At the moment, both Clenshaw Curtis and equidistant nodes do not |
| 67 | + // reproduce the expected convergence rate... something is wrong, either with |
| 68 | + // the test or with the nodes |
| 69 | + |
| 70 | + // Also: What is the ACTUAL number of quadrature nodes in both cases? |
| 71 | + case pfasst::QuadratureType::ClenshawCurtis: |
| 72 | + this->niters = this->nnodes; |
| 73 | + this->Tend = 1.0; |
| 74 | + this->lambda = complex<double>(-1.0, 1.0); |
| 75 | + this->nsteps = {3, 5, 7, 9, 11}; |
| 76 | + this->nnodes_in_call = this->nnodes + 1; |
| 77 | + break; |
| 78 | + |
| 79 | + case pfasst::QuadratureType::Uniform: |
| 80 | + this->niters = this->nnodes; |
| 81 | + this->Tend = 5.0; |
| 82 | + this->lambda = complex<double>(-1.0, 1.0); |
| 83 | + this->nsteps = {3, 5, 7, 9, 11}; |
| 84 | + this->nnodes_in_call = this->nnodes; |
| 85 | + break; |
| 86 | + |
| 87 | + |
| 88 | + default: |
| 89 | + break; |
| 90 | + } |
| 91 | + } |
| 92 | + |
| 93 | + public: |
| 94 | + virtual void SetUp() |
| 95 | + { |
| 96 | + this->nnodes = get<0>(GetParam()); |
| 97 | + this->nodetype = get<1>(GetParam()); |
| 98 | + this->set_parameters(); |
| 99 | + this->nsteps_l = this->nsteps.size(); |
| 100 | + this->err.resize(this->nsteps.size()); |
| 101 | + this->convrate.resize(this->nsteps.size()); |
| 102 | + |
| 103 | + // run to compute errors |
| 104 | + for (size_t i = 0; i <= this->nsteps_l - 1; ++i) { |
| 105 | + this->dt = this->Tend / double(this->nsteps[i]); |
| 106 | + this->err[i] = run_scalar_sdc(this->nsteps[i], this->dt, this->nnodes_in_call, |
| 107 | + this->niters, this->lambda, this->nodetype); |
| 108 | + } |
| 109 | + |
| 110 | + // compute convergence rates |
| 111 | + for (size_t i = 0; i <= this->nsteps_l - 2; ++i) { |
| 112 | + this->convrate[i] = log10(this->err[i+1] / this->err[i]) / |
| 113 | + log10(double(this->nsteps[i]) / double(this->nsteps[i + 1])); |
| 114 | + } |
| 115 | + } |
| 116 | + |
| 117 | + virtual void TearDown() |
| 118 | + {} |
| 119 | +}; |
| 120 | + |
| 121 | +/* |
| 122 | + * The test below verifies that the code approximately (up to a safety factor) reproduces |
| 123 | + * the theoretically expected rate of convergence |
| 124 | + */ |
| 125 | +TEST_P(ConvergenceTest, AllNodes) |
| 126 | +{ |
| 127 | + for (size_t i = 0; i <= nsteps_l - 2; ++i) { |
| 128 | + switch (this->nodetype) |
| 129 | + { |
| 130 | + case pfasst::QuadratureType::GaussLobatto: |
| 131 | + // Expect convergence rate of 2*nodes-2 from collocation formula, doing an identical number |
| 132 | + // of iteration should suffice to reach this as each iteration should increase order by one |
| 133 | + EXPECT_THAT(convrate[i], |
| 134 | + DoubleNear(double(2 * nnodes - 2), 0.99)) << "Convergence rate for " |
| 135 | + << this->nnodes |
| 136 | + << " Gauss-Lobatto nodes " |
| 137 | + << " at node " << i |
| 138 | + << " not within expected range."; |
| 139 | + break; |
| 140 | + |
| 141 | + case pfasst::QuadratureType::GaussLegendre: |
| 142 | + // convergence rates for Legendre nodes should be 2*nodes but are actually better, so |
| 143 | + // use Ge here |
| 144 | + EXPECT_THAT(convrate[i], Ge<double>(2 * this->nnodes)) << "Convergence rate for " |
| 145 | + << this->nnodes |
| 146 | + << " Gauss-Legendre nodes " |
| 147 | + << " at node " << i |
| 148 | + << " not within expected range."; |
| 149 | + break; |
| 150 | + |
| 151 | + |
| 152 | + case pfasst::QuadratureType::GaussRadau: |
| 153 | + // convergence rate for Radau nodes should be 2*nodes-1 |
| 154 | + // For some case, the convergence rate is ALMOST that value, hence put in the 0.99 |
| 155 | + EXPECT_THAT(convrate[i], Ge<double>(0.99* 2 * this->nnodes - 1)) << "Convergence rate for " |
| 156 | + << this->nnodes |
| 157 | + << " Gauss-Radu nodes " |
| 158 | + << " at node " << i |
| 159 | + << " not within expected range."; |
| 160 | + break; |
| 161 | + |
| 162 | + case pfasst::QuadratureType::ClenshawCurtis: |
| 163 | + // Clenshaw Curtis should be of order nnodes |
| 164 | + EXPECT_THAT(convrate[i], Ge<double>(this->nnodes)) << "Convergence rate for " |
| 165 | + << this->nnodes |
| 166 | + << " Clenshaw-Curtis nodes " |
| 167 | + << " at node " << i |
| 168 | + << " not within expected range."; |
| 169 | + break; |
| 170 | + |
| 171 | + case pfasst::QuadratureType::Uniform: |
| 172 | + // Equidistant nodes should be of order nnodes |
| 173 | + EXPECT_THAT(convrate[i], Ge<double>(this->nnodes)) << "Convergence rate for " |
| 174 | + << this->nnodes |
| 175 | + << " equidistant nodes " |
| 176 | + << " at node " << i |
| 177 | + << " not within expected range."; |
| 178 | + break; |
| 179 | + |
| 180 | + |
| 181 | + default: |
| 182 | + EXPECT_TRUE(false); |
| 183 | + break; |
| 184 | + } |
| 185 | + } |
| 186 | +} |
| 187 | + |
| 188 | +INSTANTIATE_TEST_CASE_P(ScalarSDC, ConvergenceTest, |
| 189 | + Combine(Range<size_t>(2, 7), |
| 190 | + Values(pfasst::QuadratureType::GaussLobatto, |
| 191 | + pfasst::QuadratureType::GaussLegendre, |
| 192 | + pfasst::QuadratureType::GaussRadau/*, |
| 193 | + pfasst::QuadratureType::ClenshawCurtis, |
| 194 | + pfasst::QuadratureType::Uniform*/)) |
| 195 | +); |
| 196 | + |
| 197 | +int main(int argc, char** argv) |
| 198 | +{ |
| 199 | + testing::InitGoogleTest(&argc, argv); |
| 200 | + return RUN_ALL_TESTS(); |
| 201 | +} |
| 202 | + |
| 203 | + |
0 commit comments