Skip to content

Commit ef2ad09

Browse files
memmettmemmett
committed
MLSDC: Get FAS correction working.
Signed-off-by: Matthew Emmett <[email protected]>
1 parent 10950e1 commit ef2ad09

File tree

7 files changed

+109
-41
lines changed

7 files changed

+109
-41
lines changed

‎examples/advection_diffusion/ex2.cpp‎

Lines changed: 23 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
/*
22
* Advection/diffusion example using an encapsulated IMEX sweeper.
33
*
4-
* This example uses a multi-level SDC sweeper.
4+
* This example uses a (serial) multi-level SDC sweeper.
55
*/
66

77
#include <pfasst.hpp>
@@ -21,9 +21,15 @@ int main(int argc, char **argv)
2121
const int xrat = 2;
2222
const int trat = 2;
2323

24-
int nnodes = 5;
25-
int ndofs = 128;
24+
int nnodes = 5;
25+
int ndofs = 128;
2626

27+
/*
28+
* build space/time discretisation levels and add them to mlsdc
29+
* controller. this loop adds the finest level first, and
30+
* subsequently refines in time (accoring to 'trat') and space
31+
* (according to 'xrat').
32+
*/
2733
for (int l=0; l<nlevs; l++) {
2834
auto nodes = pfasst::compute_nodes<double>(nnodes, "gauss-lobatto");
2935
auto* factory = new pfasst::encap::VectorFactory<double,double>(ndofs);
@@ -33,21 +39,32 @@ int main(int argc, char **argv)
3339
sweeper->set_nodes(nodes);
3440
sweeper->set_factory(factory);
3541

42+
mlsdc.add_level(sweeper, transfer);
43+
3644
ndofs = ndofs / xrat;
3745
nnodes = (nnodes-1) / trat + 1;
38-
39-
mlsdc.add_level(sweeper, transfer);
4046
}
4147

42-
mlsdc.set_duration(dt, nsteps, niters);
48+
/*
49+
* set up the mlsdc controller (which in turn calls 'setup' on the
50+
* sweepers that were added above). this stage typically
51+
* preallocates various buffers that the sweepers need.
52+
*/
4353
mlsdc.setup();
4454

55+
/*
56+
* set initial conditions on each level
57+
*/
4558
for (int l=0; l<nlevs; l++) {
4659
auto* sweeper = mlsdc.get_level<AdvectionDiffusionSweeper<double,double>>(l);
4760
auto* q0 = sweeper->get_state(0);
4861
sweeper->exact(q0, 0.0);
4962
}
5063

64+
/*
65+
* run mlsdc!
66+
*/
67+
mlsdc.set_duration(dt, nsteps, niters);
5168
mlsdc.run();
5269

5370
return 0;

‎include/pfasst/encap/encapsulation.hpp‎

Lines changed: 47 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -120,14 +120,16 @@ namespace pfasst {
120120
this->set_state(this->get_end_state(), 0);
121121
}
122122

123-
virtual void integrate(Encapsulation<scalar,time>* dst, time dt) {
123+
virtual void integrate(time dt, vector<Encapsulation<scalar,time>*> dst) const
124+
{
124125
throw NotImplementedYet("sweeper");
125126
}
126127
};
127128

128129
template<typename scalar, typename time>
129130
class PolyInterpMixin : public pfasst::ITransfer {
130-
matrix<time> tmat;
131+
matrix<time> tmat, fmat;
132+
131133
public:
132134

133135
virtual void interpolate(ISweeper *DST, const ISweeper *SRC, bool initial)
@@ -192,22 +194,55 @@ namespace pfasst {
192194
for (int m=0; m<ndst; m++) dst->evaluate(m);
193195
}
194196

195-
virtual void fas(ISweeper *DST, const ISweeper *SRC)
197+
virtual void fas(time dt, ISweeper *dst, const ISweeper *src)
196198
{
197-
auto* dst = dynamic_cast<EncapsulatedSweeperMixin<scalar,time>*>(DST);
198-
auto* src = dynamic_cast<const EncapsulatedSweeperMixin<scalar,time>*>(SRC);
199+
auto* crse = dynamic_cast<EncapsulatedSweeperMixin<scalar,time>*>(dst);
200+
auto* fine = dynamic_cast<const EncapsulatedSweeperMixin<scalar,time>*>(src);
199201

200-
int ndst = dst->get_nodes().size();
201-
int nsrc = src->get_nodes().size();
202+
int ncrse = crse->get_nodes().size();
203+
int nfine = fine->get_nodes().size();
202204

203-
auto* crse_factory = src->get_factory();
204-
auto* fine_factory = dst->get_factory();
205+
auto* crse_factory = crse->get_factory();
206+
auto* fine_factory = fine->get_factory();
205207

206-
vector<Encapsulation<scalar,time>*> fine_q(ndst), fine_tmp(nsrc);
208+
vector<Encapsulation<scalar,time>*> crse_z2n(ncrse-1), fine_z2n(nfine-1), rstr_z2n(ncrse-1);
209+
for (int m=0; m<ncrse-1; m++) crse_z2n[m] = crse_factory->create(solution);
210+
for (int m=0; m<ncrse-1; m++) rstr_z2n[m] = crse_factory->create(solution);
211+
for (int m=0; m<nfine-1; m++) fine_z2n[m] = fine_factory->create(solution);
207212

208-
for (int m=0; m<ndst; m++) fine_q[m] = dst->get_state(m);
209-
for (int m=0; m<nsrc; m++) fine_tmp[m] = fine_factory->create(solution);
213+
// compute '0 to node' integral on the coarse level
214+
crse->integrate(dt, crse_z2n);
215+
for (int m=1; m<ncrse-1; m++)
216+
crse_z2n[m]->saxpy(1.0, crse_z2n[m-1]);
217+
218+
// compute '0 to node' integral on the fine level
219+
fine->integrate(dt, fine_z2n);
220+
for (int m=1; m<nfine-1; m++)
221+
fine_z2n[m]->saxpy(1.0, fine_z2n[m-1]);
222+
223+
// restrict '0 to node' fine integral
224+
int trat = (nfine - 1) / (ncrse - 1);
225+
for (int m=1; m<ncrse; m++)
226+
this->restrict(rstr_z2n[m-1], fine_z2n[m*trat-1]);
227+
228+
// compute 'node to node' tau correction
229+
vector<Encapsulation<scalar,time>*> tau(ncrse-1);
230+
for (int m=0; m<ncrse-1; m++) tau[m] = crse->get_tau(m);
231+
232+
tau[0]->copy(rstr_z2n[0]);
233+
tau[0]->saxpy(-1.0, crse_z2n[0]);
234+
235+
for (int m=1; m<ncrse-1; m++) {
236+
tau[m]->copy(rstr_z2n[m]);
237+
tau[m]->saxpy(-1.0, rstr_z2n[m-1]);
238+
239+
tau[m]->saxpy(-1.0, crse_z2n[m]);
240+
tau[m]->saxpy(1.0, crse_z2n[m-1]);
241+
}
210242

243+
for (int m=0; m<ncrse-1; m++) delete crse_z2n[m];
244+
for (int m=0; m<ncrse-1; m++) delete rstr_z2n[m];
245+
for (int m=0; m<nfine-1; m++) delete fine_z2n[m];
211246
}
212247

213248
// required for interp/restrict helpers

‎include/pfasst/encap/imex.hpp‎

Lines changed: 15 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -49,6 +49,13 @@ namespace pfasst {
4949
return pQ[m];
5050
}
5151

52+
virtual void integrate(time dt, vector<Encapsulation<scalar,time>*> dst) const
53+
{
54+
auto* encap = dst[0];
55+
encap->mat_apply(dst, dt, Smat, Fe, true);
56+
encap->mat_apply(dst, dt, Smat, Fi, false);
57+
}
58+
5259
void setup(bool coarse)
5360
{
5461
auto nodes = this->get_nodes();
@@ -84,15 +91,12 @@ namespace pfasst {
8491
const int nnodes = nodes.size();
8592

8693
// integrate
87-
for (int n=0; n<nnodes-1; n++) {
88-
S[n]->setval(0.0);
89-
for (int m=0; m<nnodes; m++) {
90-
S[n]->saxpy(dt * SEmat(n,m), Fe[m]);
91-
S[n]->saxpy(dt * SImat(n,m), Fi[m]);
92-
}
93-
if (T.size() > 0)
94-
S[n]->saxpy(1.0, T[n]);
95-
}
94+
S[0]->mat_apply(S, dt, SEmat, Fe, true);
95+
S[0]->mat_apply(S, dt, SImat, Fi, false);
96+
if (T.size() > 0)
97+
for (int m=0; m<nnodes-1; m++)
98+
S[m]->saxpy(1.0, T[m]);
99+
96100

97101
// sweep
98102
Encapsulation<scalar,time> *rhs = this->get_factory()->create(pfasst::encap::solution);
@@ -143,7 +147,8 @@ namespace pfasst {
143147
pQ[m]->copy(Q[m]);
144148
}
145149

146-
virtual void evaluate(int m) {
150+
virtual void evaluate(int m)
151+
{
147152
// XXX: time
148153
f1eval(Fe[m], Q[m], 0.0);
149154
f2eval(Fi[m], Q[m], 0.0);

‎include/pfasst/encap/vector.hpp‎

Lines changed: 19 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -26,18 +26,24 @@ namespace pfasst {
2626
class VectorEncapsulation : public vector<scalar>, public Encapsulation<scalar,time> {
2727

2828
public:
29-
VectorEncapsulation(int size) : vector<scalar>(size) { }
29+
VectorEncapsulation(int size) : vector<scalar>(size)
30+
{
31+
setval(0.0);
32+
}
3033

31-
void setval(scalar v) {
34+
void setval(scalar v)
35+
{
3236
std::fill(this->begin(), this->end(), v);
3337
}
3438

35-
void copy(const Encapsulation<scalar,time>* X) {
39+
void copy(const Encapsulation<scalar,time>* X)
40+
{
3641
const auto* x = dynamic_cast<const VectorEncapsulation*>(X);
3742
std::copy(x->begin(), x->end(), this->begin());
3843
}
3944

40-
void saxpy(time a, const Encapsulation<scalar,time> *X) {
45+
void saxpy(time a, const Encapsulation<scalar,time> *X)
46+
{
4147
const auto& x = *dynamic_cast<const VectorEncapsulation*>(X);
4248
auto& y = *this;
4349

@@ -52,20 +58,24 @@ namespace pfasst {
5258
int nsrc = SRC.size();
5359

5460
vector<VectorEncapsulation<scalar,time>*> dst(ndst), src(nsrc);
55-
for (int n=0; n<ndst; n++)
56-
dst[n] = dynamic_cast<VectorEncapsulation*>(DST[n]);
57-
for (int m=0; m<nsrc; m++)
58-
src[m] = dynamic_cast<VectorEncapsulation*>(SRC[m]);
61+
for (int n=0; n<ndst; n++) dst[n] = dynamic_cast<VectorEncapsulation<scalar,time>*>(DST[n]);
62+
for (int m=0; m<nsrc; m++) src[m] = dynamic_cast<VectorEncapsulation<scalar,time>*>(SRC[m]);
5963

6064
if (zero)
6165
for (int n=0; n<ndst; n++)
6266
dst[n]->setval(0.0);
6367

68+
// for (int m=0; m<ndst; m++) cout << DST[m] << endl;
69+
// for (int m=0; m<nsrc; m++) cout << SRC[m] << endl;
70+
71+
// for (int m=0; m<ndst; m++) cout << dst[m]->size() << endl;
72+
// for (int m=0; m<nsrc; m++) cout << src[m]->size() << endl;
73+
6474
int ndofs = (*dst[0]).size();
6575
for (int i=0; i<ndofs; i++)
6676
for (int n=0; n<ndst; n++)
6777
for (int m=0; m<nsrc; m++)
68-
(*dst[n])[i] += a * mat(n,m) * (*src[m])[i];
78+
dst[n]->data()[i] += a * mat(n,m) * src[m]->data()[i];
6979
}
7080

7181
scalar norm0() const

‎include/pfasst/interfaces.hpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -30,8 +30,8 @@ namespace pfasst {
3030
virtual void setup(bool coarse=false) { }
3131
virtual ~ISweeper() { }
3232

33-
virtual void sweep(double t, double dt) = 0;
34-
virtual void predict(double t, double dt) = 0;
33+
virtual void sweep(double t, double dt) = 0; // XXX: this needs to be a templated
34+
virtual void predict(double t, double dt) = 0; // XXX: this needs to be templated
3535
virtual void advance() = 0;
3636
virtual void save() { NotImplementedYet("mlsdc/pfasst"); }
3737
};
@@ -41,7 +41,7 @@ namespace pfasst {
4141
// XXX: pass level iterator to these routines as well
4242
virtual void interpolate(ISweeper *dst, const ISweeper *src, bool initial) = 0;
4343
virtual void restrict(ISweeper *dst, const ISweeper *src) = 0;
44-
virtual void fas(ISweeper *dst, const ISweeper *src) = 0;
44+
virtual void fas(double dt, ISweeper *dst, const ISweeper *src) = 0; // XXX: this needs to be templated
4545
};
4646

4747
class ICommunicator {

‎include/pfasst/mlsdc.hpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -79,7 +79,7 @@ namespace pfasst {
7979
fine->sweep(t, dt);
8080

8181
trns->restrict(crse, fine);
82-
trns->fas(crse, fine);
82+
trns->fas(dt, crse, fine);
8383
crse->save();
8484

8585
return leviter - 1;

‎include/pfasst/quadrature.hpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,7 @@
99
#include <vector>
1010

1111
#include <boost/numeric/ublas/matrix.hpp>
12+
#include <boost/numeric/ublas/io.hpp>
1213

1314
using std::complex;
1415
using std::string;

0 commit comments

Comments
 (0)