Skip to content

Commit 8792d7b

Browse files
torbjoernktorbjoernk
committed
Merge pull request #141 from memmett/feature/add_fas_residual
IMEX: Add FAS corrections to residual.
2 parents b7ffc09 + b8b5c76 commit 8792d7b

File tree

6 files changed

+97
-18
lines changed

6 files changed

+97
-18
lines changed

‎examples/advection_diffusion/advection_diffusion_sweeper.hpp‎

Lines changed: 38 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -35,11 +35,13 @@ namespace pfasst
3535
namespace advection_diffusion
3636
{
3737
/**
38-
* errors at different iterations and time nodes
39-
*
40-
* Mapping a pair of step/iteration indices onto the error of the solution.
38+
* Containers for errors/residuals etc.
4139
*/
42-
typedef map<pair<size_t, size_t>, double> error_map;
40+
typedef map<tuple<size_t, size_t>, double> error_map; // step, iteration -> error
41+
typedef map<size_t, error_map> residual_map; // level, (step, iteration) -> residual
42+
43+
typedef error_map::value_type vtype;
44+
typedef error_map::key_type ktype;
4345

4446
template<typename time = pfasst::time_precision>
4547
class AdvectionDiffusionSweeper
@@ -70,6 +72,7 @@ namespace pfasst
7072

7173
//! @{
7274
error_map errors;
75+
error_map residuals;
7376
//! @}
7477

7578
//! @{
@@ -138,18 +141,42 @@ namespace pfasst
138141

139142
auto n = this->get_controller()->get_step();
140143
auto k = this->get_controller()->get_iteration();
141-
LOG(INFO) << "err:" << n << k << max << "(" << qend.size() << "," << predict << ")";
144+
LOG(INFO) << "err: " << n << " " << k << " " << max << " (" << qend.size() << "," << predict << ")";
142145

143-
this->errors.insert(pair<pair<size_t, size_t>, double>(pair<size_t, size_t>(n, k), max));
146+
this->errors.insert(vtype(ktype(n, k), max));
147+
}
148+
149+
void echo_residual()
150+
{
151+
vector<shared_ptr<Encapsulation<time>>> residuals;
152+
153+
for (size_t m = 0; m < this->get_nodes().size(); m++) {
154+
residuals.push_back(this->get_factory()->create(pfasst::encap::solution));
155+
}
156+
this->residual(this->get_controller()->get_time_step(), residuals);
157+
158+
vector<time> rnorms;
159+
for (auto r: residuals) {
160+
rnorms.push_back(r->norm0());
161+
}
162+
auto rmax = *std::max_element(rnorms.begin(), rnorms.end());
163+
164+
auto n = this->get_controller()->get_step();
165+
auto k = this->get_controller()->get_iteration();
166+
LOG(INFO) << "res: " << n << " " << k << " " << rmax << " (" << residuals.size() << ")";
167+
168+
this->residuals[ktype(n, k)] = rmax;
144169
}
145170

146-
/**
147-
* retrieve errors at iterations and time nodes
148-
*/
149171
error_map get_errors()
150172
{
151173
return this->errors;
152174
}
175+
176+
error_map get_residuals()
177+
{
178+
return this->residuals;
179+
}
153180
//! @}
154181

155182
//! @{
@@ -161,6 +188,7 @@ namespace pfasst
161188
time t = this->get_controller()->get_time();
162189
time dt = this->get_controller()->get_time_step();
163190
this->echo_error(t + dt, true);
191+
this->echo_residual();
164192
}
165193

166194
/**
@@ -171,6 +199,7 @@ namespace pfasst
171199
time t = this->get_controller()->get_time();
172200
time dt = this->get_controller()->get_time_step();
173201
this->echo_error(t + dt);
202+
this->echo_residual();
174203
}
175204
//! @}
176205

‎examples/advection_diffusion/serial_mlsdc.cpp‎

Lines changed: 9 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -24,14 +24,13 @@ namespace pfasst
2424
{
2525
namespace advection_diffusion
2626
{
27-
error_map run_serial_mlsdc()
27+
tuple<error_map, residual_map> run_serial_mlsdc(size_t nlevs)
2828
{
2929
MLSDC<> mlsdc;
3030

31-
const size_t nlevs = 2;
3231
const size_t nsteps = 4;
3332
const double dt = 0.01;
34-
const size_t niters = 4;
33+
const size_t niters = 8;
3534
const int xrat = 2;
3635
const int trat = 2;
3736

@@ -82,7 +81,12 @@ namespace pfasst
8281

8382
fftw_cleanup();
8483

85-
return sweeper->get_errors();
84+
tuple<error_map, residual_map> rinfo;
85+
get<0>(rinfo) = mlsdc.get_finest<AdvectionDiffusionSweeper<>>()->get_errors();
86+
for (auto l = mlsdc.coarsest(); l <= mlsdc.finest(); ++l) {
87+
get<1>(rinfo).insert(pair<size_t, error_map>(l.level, l.current<AdvectionDiffusionSweeper<>>()->get_residuals()));
88+
}
89+
return rinfo;
8690
}
8791
} // ::pfasst::examples::advection_diffusion
8892
} // ::pfasst::examples
@@ -91,6 +95,6 @@ namespace pfasst
9195
#ifndef PFASST_UNIT_TESTING
9296
int main(int /*argc*/, char** /*argv*/)
9397
{
94-
pfasst::examples::advection_diffusion::run_serial_mlsdc();
98+
pfasst::examples::advection_diffusion::run_serial_mlsdc(3);
9599
}
96100
#endif

‎include/pfasst/encap/imex_sweeper.hpp‎

Lines changed: 8 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -206,7 +206,14 @@ namespace pfasst
206206
for (size_t m = 0; m < this->quadrature->get_num_nodes(); m++) {
207207
dst[m]->copy(this->start_state);
208208
dst[m]->saxpy(-1.0, this->state[m]);
209-
// XXX: add tau corrections
209+
}
210+
if (this->fas_corrections.size() > 0) {
211+
// XXX: build a matrix and use mat_apply to do this
212+
for (size_t m = 0; m < this->quadrature->get_num_nodes(); m++) {
213+
for (size_t n = 0; n <= m; n++) {
214+
dst[m]->saxpy(1.0, this->fas_corrections[n]);
215+
}
216+
}
210217
}
211218
dst[0]->mat_apply(dst, dt, this->quadrature->get_q_mat(), this->fs_expl, false);
212219
dst[0]->mat_apply(dst, dt, this->quadrature->get_q_mat(), this->fs_impl, false);

‎include/pfasst/mlsdc.hpp‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -70,7 +70,9 @@ namespace pfasst
7070

7171
perform_sweeps(this->finest().level);
7272

73-
this->get_finest()->post_step();
73+
for (auto l = this->finest(); l >= this->coarsest(); --l) {
74+
l.current()->post_step();
75+
}
7476

7577
if (this->get_time() + this->get_time_step() < this->get_end_time()) {
7678
this->get_finest()->advance();

‎include/pfasst/pfasst.hpp‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -78,7 +78,9 @@ namespace pfasst
7878
cycle_v(this->finest());
7979
}
8080

81-
this->get_finest()->post_step();
81+
for (auto l = this->finest(); l >= this->coarsest(); --l) {
82+
l.current()->post_step();
83+
}
8284

8385
if (nblock < nblocks - 1) {
8486
broadcast();

‎tests/examples/advection_diffusion/test_advection_diffusion.cpp‎

Lines changed: 36 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -74,7 +74,8 @@ TEST(ErrorTest, SerialMLSDC)
7474
{
7575
typedef error_map::value_type vtype;
7676

77-
auto errors = run_serial_mlsdc();
77+
auto errors_and_residuals = run_serial_mlsdc(2);
78+
auto errors = get<0>(errors_and_residuals);
7879
auto get_iter = [](const vtype x) { return get<1>(get<0>(x)); };
7980
auto get_error = [](const vtype x) { return get<1>(x); };
8081

@@ -92,6 +93,40 @@ TEST(ErrorTest, SerialMLSDC)
9293
EXPECT_THAT(err, testing::Pointwise(DoubleLess(), tol));
9394
}
9495

96+
TEST(FASTest, SerialMLSDC)
97+
{
98+
typedef error_map::key_type ktype;
99+
100+
auto errors_and_residuals = run_serial_mlsdc(3);
101+
auto residuals = get<1>(errors_and_residuals);
102+
103+
ASSERT_NEAR(residuals[2][ktype(3, 0)], 0.000667207, 1.e-8);
104+
ASSERT_NEAR(residuals[0][ktype(3, 0)], 6.23966e-07, 1.e-12);
105+
ASSERT_NEAR(residuals[1][ktype(3, 0)], 1.27783e-08, 1.e-12);
106+
ASSERT_NEAR(residuals[2][ktype(3, 1)], 6.60607e-07, 1.e-12);
107+
ASSERT_NEAR(residuals[0][ktype(3, 1)], 5.19702e-10, 1.e-14);
108+
ASSERT_NEAR(residuals[1][ktype(3, 1)], 2.59963e-10, 1.e-12);
109+
ASSERT_NEAR(residuals[2][ktype(3, 2)], 8.89424e-09, 1.e-12);
110+
ASSERT_NEAR(residuals[0][ktype(3, 2)], 8.28716e-11, 1.e-14);
111+
ASSERT_NEAR(residuals[1][ktype(3, 2)], 4.54949e-11, 1.e-14);
112+
ASSERT_NEAR(residuals[2][ktype(3, 3)], 1.04101e-10, 1.e-12);
113+
ASSERT_NEAR(residuals[0][ktype(3, 3)], 8.35953e-11, 1.e-15);
114+
ASSERT_NEAR(residuals[1][ktype(3, 3)], 4.20877e-11, 1.e-15);
115+
ASSERT_NEAR(residuals[2][ktype(3, 4)], 2.18056e-12, 1.e-15);
116+
ASSERT_NEAR(residuals[0][ktype(3, 4)], 8.34365e-11, 1.e-15);
117+
ASSERT_NEAR(residuals[1][ktype(3, 4)], 4.19699e-11, 1.e-15);
118+
ASSERT_NEAR(residuals[2][ktype(3, 5)], 7.18701e-13, 1.e-15);
119+
ASSERT_NEAR(residuals[0][ktype(3, 5)], 8.34336e-11, 1.e-15);
120+
ASSERT_NEAR(residuals[1][ktype(3, 5)], 4.19691e-11, 1.e-15);
121+
ASSERT_NEAR(residuals[2][ktype(3, 6)], 7.07797e-13, 1.e-15);
122+
ASSERT_NEAR(residuals[0][ktype(3, 6)], 8.34340e-11, 1.e-15);
123+
ASSERT_NEAR(residuals[1][ktype(3, 6)], 4.19693e-11, 1.e-15);
124+
ASSERT_NEAR(residuals[2][ktype(3, 7)], 7.07356e-13, 1.e-15);
125+
ASSERT_NEAR(residuals[0][ktype(3, 7)], 8.34338e-11, 1.e-15);
126+
ASSERT_NEAR(residuals[1][ktype(3, 7)], 4.19698e-11, 1.e-15);
127+
ASSERT_NEAR(residuals[2][ktype(3, 8)], 7.07458e-13, 1.e-15);
128+
}
129+
95130
int main(int argc, char** argv)
96131
{
97132
testing::InitGoogleTest(&argc, argv);

0 commit comments

Comments
 (0)