Skip to content

Commit 00462bd

Browse files
torbjoernktorbjoernk
committed
Merge remote-tracking branch 'memmett/feature/fix-quadrature' into development
* memmett/feature/fix-quadrature: Style. quadrature: Initialize matrix with zeros. imex: Support virtual nodes (0 and 1). quadrature: Add `pfasst::augment_nodes` and handle `is_proper` flags appropriately. fixes PR #68 Signed-off-by: Torbjörn Klatt <[email protected]>
2 parents 6279fdf + c12134d commit 00462bd

File tree

4 files changed

+104
-28
lines changed

4 files changed

+104
-28
lines changed

‎include/pfasst/encap/encap_sweeper.hpp‎

Lines changed: 9 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -25,6 +25,7 @@ namespace pfasst
2525
{
2626
//! @{
2727
vector<time> nodes;
28+
vector<bool> is_proper;
2829
shared_ptr<EncapFactory<time>> factory;
2930
//! @}
3031

@@ -69,14 +70,21 @@ namespace pfasst
6970
//! @{
7071
void set_nodes(vector<time> nodes)
7172
{
72-
this->nodes = nodes;
73+
auto augmented = pfasst::augment_nodes(nodes);
74+
this->nodes = get<0>(augmented);
75+
this->is_proper = get<1>(augmented);
7376
}
7477

7578
const vector<time> get_nodes() const
7679
{
7780
return nodes;
7881
}
7982

83+
const vector<bool> get_is_proper() const
84+
{
85+
return is_proper;
86+
}
87+
8088
void set_factory(shared_ptr<EncapFactory<time>> factory)
8189
{
8290
this->factory = shared_ptr<EncapFactory<time>>(factory);

‎include/pfasst/encap/imex_sweeper.hpp‎

Lines changed: 57 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -30,9 +30,9 @@ namespace pfasst
3030
* explicitly and the stiff part implicitly.
3131
* Therefore, we define the splitting \\( F(t,u) = F_{expl}(t,u) + F_{impl}(t,u) \\).
3232
*
33-
* This sweeper requires three interfaces to be implement for such ODEs: two routines to
34-
* evaluate the explicit \\( F_{\\rm expl} \\) and implicit \\( F_{\\rm impl} \\) pieces for a
35-
* given state, and one that solves (perhaps with an external solver) the backward-Euler
33+
* This sweeper requires three interfaces to be implement for such ODEs: two routines to
34+
* evaluate the explicit \\( F_{\\rm expl} \\) and implicit \\( F_{\\rm impl} \\) pieces for a
35+
* given state, and one that solves (perhaps with an external solver) the backward-Euler
3636
* equation \\( U^{n+1} - \\Delta t F_{\\rm impl}(U^{n+1}) = RHS \\) for \\( U^{n+1} \\).
3737
*
3838
* @tparam time precision type of the time dimension
@@ -86,6 +86,11 @@ namespace pfasst
8686
*/
8787
matrix<time> s_mat;
8888

89+
/**
90+
* quadrature matrix containing weights for 0-to-last node integration
91+
*/
92+
matrix<time> b_mat;
93+
8994
/**
9095
* quadrature matrix containing weights for node-to-node integration of explicit part
9196
*
@@ -101,6 +106,24 @@ namespace pfasst
101106
matrix<time> s_mat_impl;
102107
//! @}
103108

109+
110+
/**
111+
* Set end state to \\( U_0 + \\int F_{expl} + F_{expl} \\).
112+
*/
113+
void integrate_end_state(time dt)
114+
{
115+
vector<shared_ptr<Encapsulation<time>>> dst = { this->u_state.back() };
116+
dst[0]->copy(this->u_state.front());
117+
dst[0]->mat_apply(dst, dt, this->b_mat, this->fs_expl, false);
118+
dst[0]->mat_apply(dst, dt, this->b_mat, this->fs_impl, false);
119+
}
120+
121+
inline bool last_node_is_virtual()
122+
{
123+
return !this->get_is_proper().back();
124+
}
125+
126+
104127
public:
105128
//! @{
106129
virtual ~IMEXSweeper()
@@ -173,9 +196,16 @@ namespace pfasst
173196
virtual void setup(bool coarse) override
174197
{
175198
auto nodes = this->get_nodes();
199+
auto is_proper = this->get_is_proper();
176200
assert(nodes.size() >= 1);
177201

178-
this->s_mat = compute_quadrature(nodes, nodes, 's');
202+
this->s_mat = compute_quadrature(nodes, nodes, is_proper, 's');
203+
204+
auto q_mat = compute_quadrature(nodes, nodes, is_proper, 'q');
205+
this->b_mat = matrix<time>(1, nodes.size());
206+
for (size_t m = 0; m < nodes.size(); m++) {
207+
this->b_mat(0, m) = q_mat(nodes.size() - 2, m);
208+
}
179209

180210
this->s_mat_expl = this->s_mat;
181211
this->s_mat_impl = this->s_mat;
@@ -190,6 +220,7 @@ namespace pfasst
190220
if (coarse) {
191221
this->previous_u_state.push_back(this->get_factory()->create(pfasst::encap::solution));
192222
}
223+
// XXX: can prolly save some f's here for virtual nodes...
193224
this->fs_expl.push_back(this->get_factory()->create(pfasst::encap::function));
194225
this->fs_impl.push_back(this->get_factory()->create(pfasst::encap::function));
195226
}
@@ -211,14 +242,11 @@ namespace pfasst
211242
time dt = this->get_controller()->get_time_step();
212243
time t = this->get_controller()->get_time();
213244

214-
if (initial) {
215-
this->f_expl_eval(this->fs_expl[0], this->u_state[0], t);
216-
this->f_impl_eval(this->fs_impl[0], this->u_state[0], t);
217-
}
245+
if (initial) { this->evaluate(0); }
218246

219247
shared_ptr<Encapsulation<time>> rhs = this->get_factory()->create(pfasst::encap::solution);
220248

221-
for (size_t m = 0; m < nnodes - 1; m++) {
249+
for (size_t m = 0; m < nnodes - (this->last_node_is_virtual() ? 2 : 1); m++) {
222250
time ds = dt * (nodes[m + 1] - nodes[m]);
223251
rhs->copy(this->u_state[m]);
224252
rhs->saxpy(ds, this->fs_expl[m]);
@@ -227,6 +255,8 @@ namespace pfasst
227255

228256
t += ds;
229257
}
258+
259+
if (this->last_node_is_virtual()) { this->integrate_end_state(dt); }
230260
}
231261

232262
virtual void sweep()
@@ -250,7 +280,7 @@ namespace pfasst
250280
// sweep
251281
shared_ptr<Encapsulation<time>> rhs = this->get_factory()->create(pfasst::encap::solution);
252282

253-
for (size_t m = 0; m < nnodes - 1; m++) {
283+
for (size_t m = 0; m < nnodes - (this->last_node_is_virtual() ? 2 : 1); m++) {
254284
time ds = dt * (nodes[m + 1] - nodes[m]);
255285

256286
rhs->copy(this->u_state[m]);
@@ -261,6 +291,8 @@ namespace pfasst
261291

262292
t += ds;
263293
}
294+
295+
if (this->last_node_is_virtual()) { this->integrate_end_state(dt); }
264296
}
265297

266298
virtual void advance() override
@@ -281,11 +313,23 @@ namespace pfasst
281313
}
282314
}
283315

316+
/**
317+
* @copybrief EncapSweeper::evaluate()
318+
*
319+
* If the node `m` is virtual (not proper), we can save some
320+
* implicit/explicit evaluations.
321+
*/
284322
virtual void evaluate(size_t m) override
285323
{
286-
time t = this->get_nodes()[m]; // XXX
287-
this->f_expl_eval(this->fs_expl[m], this->u_state[m], t);
288-
this->f_impl_eval(this->fs_impl[m], this->u_state[m], t);
324+
time t0 = this->get_controller()->get_time();
325+
time dt = this->get_controller()->get_time_step();
326+
time t = t0 + dt * this->get_nodes()[m];
327+
if ((m == 0) || this->get_is_proper()[m]) {
328+
this->f_expl_eval(this->fs_expl[m], this->u_state[m], t);
329+
}
330+
if (this->get_is_proper()[m]) {
331+
this->f_impl_eval(this->fs_impl[m], this->u_state[m], t);
332+
}
289333
}
290334

291335
/**

‎include/pfasst/quadrature.hpp‎

Lines changed: 26 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -195,12 +195,12 @@ namespace pfasst
195195

196196
} else if (qtype == "clenshaw-curtis") {
197197
for (size_t j = 0; j < nnodes; j++) {
198-
nodes[j] = 0.5 * (1.0 - cos(j * PI / (nnodes-1)));
198+
nodes[j] = 0.5 * (1.0 - cos(j * PI / (nnodes - 1)));
199199
}
200200

201201
} else if (qtype == "uniform") {
202202
for (size_t j = 0; j < nnodes; j++) {
203-
nodes[j] = node(j) / (nnodes-1);
203+
nodes[j] = node(j) / (nnodes - 1);
204204
}
205205

206206
} else {
@@ -210,23 +210,44 @@ namespace pfasst
210210
return nodes;
211211
}
212212

213+
template<typename node>
214+
auto augment_nodes(vector<node> const orig) -> pair<vector<node>, vector<bool>> {
215+
vector<node> nodes = orig;
216+
217+
bool left = nodes.front() == node(0.0);
218+
bool right = nodes.back() == node(1.0);
219+
220+
if (!left) { nodes.insert(nodes.begin(), node(0.0)); }
221+
if (!right) { nodes.insert(nodes.end(), node(1.0)); }
222+
223+
vector<bool> is_proper(nodes.size(), true);
224+
is_proper.front() = left;
225+
is_proper.back() = right;
226+
227+
return pair<vector<node>, vector<bool>>(nodes, is_proper);
228+
}
229+
230+
213231
template<typename node = time_precision>
214-
matrix<node> compute_quadrature(vector<node> dst, vector<node> src, char type)
232+
matrix<node> compute_quadrature(vector<node> dst, vector<node> src, vector<bool> is_proper,
233+
char type)
215234
{
216235
const size_t ndst = dst.size();
217236
const size_t nsrc = src.size();
218237

219238
assert(ndst >= 1);
220-
matrix<node> mat(ndst - 1, nsrc);
239+
matrix<node> mat(ndst - 1, nsrc, node(0.0));
221240

222241
Polynomial<node> p(nsrc + 1), p1(nsrc + 1);
223242

224243
for (size_t i = 0; i < nsrc; i++) {
244+
if (!is_proper[i]) { continue; }
245+
225246
// construct interpolating polynomial coefficients
226247
p[0] = 1.0;
227248
for (size_t j = 1; j < nsrc + 1; j++) { p[j] = 0.0; }
228249
for (size_t m = 0; m < nsrc; m++) {
229-
if (m == i) { continue; }
250+
if ((!is_proper[m]) || (m == i)) { continue; }
230251

231252
// p_{m+1}(x) = (x - x_j) * p_m(x)
232253
p1[0] = 0.0;

‎tests/test_quadrature.cpp‎

Lines changed: 12 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -195,19 +195,19 @@ TEST(NodesTest, UniformNodes)
195195
{
196196
const long double u2e[2] = { 0.0,
197197
1.0
198-
};
198+
};
199199

200200
const long double u3e[3] = { 0.0,
201201
0.5,
202202
1.0
203-
};
203+
};
204204

205205
const long double u5e[5] = { 0.0,
206206
0.25,
207207
0.5,
208208
0.75,
209209
1.0
210-
};
210+
};
211211

212212
auto u2 = pfasst::compute_nodes<long double>(2, "uniform");
213213
EXPECT_THAT(u2, testing::Pointwise(DoubleNear(), u2e));
@@ -222,7 +222,8 @@ TEST(NodesTest, UniformNodes)
222222
TEST(QuadratureTest, GaussLobattoNodes)
223223
{
224224
auto l3 = pfasst::compute_nodes<long double>(3, "gauss-lobatto");
225-
auto s3 = pfasst::compute_quadrature(l3, l3, 's');
225+
auto a3 = pfasst::augment_nodes(l3);
226+
auto s3 = pfasst::compute_quadrature(get<0>(a3), get<0>(a3), get<1>(a3), 's');
226227
const long double s3e[6] = { 0.20833333333333333,
227228
0.33333333333333333,
228229
-0.04166666666666666,
@@ -234,7 +235,8 @@ TEST(QuadratureTest, GaussLobattoNodes)
234235
EXPECT_THAT(s3.data(), testing::Pointwise(DoubleNear(), s3e));
235236

236237
auto l5 = pfasst::compute_nodes<long double>(5, "gauss-lobatto");
237-
auto s5 = pfasst::compute_quadrature(l5, l5, 's');
238+
auto a5 = pfasst::augment_nodes(l5);
239+
auto s5 = pfasst::compute_quadrature(get<0>(a5), get<0>(a5), get<1>(a5), 's');
238240
const long double s5e[] = { 0.067728432186156897969267419174073482,
239241
0.11974476934341168251615379970493965,
240242
-0.021735721866558113665511351745074292,
@@ -255,16 +257,17 @@ TEST(QuadratureTest, GaussLobattoNodes)
255257
-0.021735721866558113665511351745074289,
256258
0.11974476934341168251615379970493965,
257259
0.067728432186156897969267419174073482
258-
};
260+
};
259261
EXPECT_THAT(s5.data(), testing::Pointwise(DoubleNear(), s5e));
260262
}
261263

262264
TEST(QuadratureTest, ClenshawCurtisNodes)
263265
{
264266
auto c4 = pfasst::compute_nodes<long double>(4, "clenshaw-curtis");
265-
auto s4 = pfasst::compute_quadrature(c4, c4, 's');
267+
auto a4 = pfasst::augment_nodes(c4);
268+
auto s4 = pfasst::compute_quadrature(get<0>(a4), get<0>(a4), get<1>(a4), 's');
266269
const long double s4e[] = { 0.10243055555555555555555555555555556,
267-
0.16319444444444444444444444444444444,
270+
0.16319444444444444444444444444444444,
268271
-0.024305555555555555555555555555555556,
269272
0.0086805555555555555555555555555555557,
270273
-0.055555555555555555555555555555555556,
@@ -275,7 +278,7 @@ TEST(QuadratureTest, ClenshawCurtisNodes)
275278
-0.024305555555555555555555555555555554,
276279
0.16319444444444444444444444444444444,
277280
0.10243055555555555555555555555555556
278-
};
281+
};
279282
EXPECT_THAT(s4.data(), testing::Pointwise(DoubleNear(), s4e));
280283
}
281284

0 commit comments

Comments
 (0)