Skip to content

Commit 9bd90a8

Browse files
memmettmemmett
committed
mlsdc: Begin add 'initial' flags. Tidy various.
Signed-off-by: Matthew Emmett <[email protected]>
1 parent 2e8f2df commit 9bd90a8

File tree

10 files changed

+126
-105
lines changed

10 files changed

+126
-105
lines changed

‎examples/advection_diffusion/CMakeLists.txt‎

Lines changed: 15 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -9,10 +9,6 @@ set(ex1_SOURCES
99
${CMAKE_CURRENT_SOURCE_DIR}/ex1.cpp
1010
)
1111

12-
set(ex2_SOURCES
13-
${CMAKE_CURRENT_SOURCE_DIR}/ex2.cpp
14-
)
15-
1612
add_executable(ex1 ${ex1_SOURCES})
1713
add_dependencies(ex1 fftw3)
1814
target_link_libraries(ex1
@@ -22,6 +18,9 @@ set_target_properties(ex1
2218
PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${pfasst_BINARY_DIR}/examples
2319
)
2420

21+
set(ex2_SOURCES
22+
${CMAKE_CURRENT_SOURCE_DIR}/ex2.cpp
23+
)
2524
add_executable(ex2 ${ex2_SOURCES})
2625
add_dependencies(ex2 fftw3)
2726
target_link_libraries(ex2
@@ -30,3 +29,15 @@ target_link_libraries(ex2
3029
set_target_properties(ex2
3130
PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${pfasst_BINARY_DIR}/examples
3231
)
32+
33+
set(ex3_SOURCES
34+
${CMAKE_CURRENT_SOURCE_DIR}/ex3.cpp
35+
)
36+
add_executable(ex3 ${ex3_SOURCES})
37+
add_dependencies(ex3 fftw3)
38+
target_link_libraries(ex3
39+
${fftw3_LIBS}
40+
)
41+
set_target_properties(ex3
42+
PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${pfasst_BINARY_DIR}/examples
43+
)

‎examples/advection_diffusion/advection_diffusion_sweeper.hpp‎

Lines changed: 8 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,7 @@ class AdvectionDiffusionSweeper : public pfasst::imex::IMEX<scalar,time> {
3030
scalar v = 1.0;
3131
time t0 = 1.0;
3232
scalar nu = 0.02;
33+
int nf1evals = 0;
3334

3435
public:
3536

@@ -42,7 +43,6 @@ class AdvectionDiffusionSweeper : public pfasst::imex::IMEX<scalar,time> {
4243
ddx[i] = complex<scalar>(0.0, 1.0) * kx;
4344
lap[i] = (kx*kx < 1e-13) ? 0.0 : -kx*kx;
4445
}
45-
4646
}
4747

4848
void exact(Encapsulation<scalar,time>* q, scalar t)
@@ -66,7 +66,7 @@ class AdvectionDiffusionSweeper : public pfasst::imex::IMEX<scalar,time> {
6666
}
6767
}
6868

69-
void echo_error(time t)
69+
void echo_error(time t, bool predict=false)
7070
{
7171
auto& qend = *dynamic_cast<dvector*>(this->get_state(this->get_nodes().size()-1));
7272
auto qex = dvector(qend.size());
@@ -79,13 +79,13 @@ class AdvectionDiffusionSweeper : public pfasst::imex::IMEX<scalar,time> {
7979
if (d > max)
8080
max = d;
8181
}
82-
cout << "err: " << scientific << max << " (" << qend.size() << ")" << endl;
82+
cout << "err: " << scientific << max << " (" << qend.size() << ", " << predict << ")" << endl;
8383
}
8484

85-
void predict(time t, time dt)
85+
void predict(time t, time dt, bool initial)
8686
{
87-
pfasst::imex::IMEX<scalar,time>::predict(t, dt);
88-
echo_error(t+dt);
87+
pfasst::imex::IMEX<scalar,time>::predict(t, dt, initial);
88+
echo_error(t+dt, true);
8989
}
9090

9191
void sweep(time t, time dt)
@@ -105,6 +105,8 @@ class AdvectionDiffusionSweeper : public pfasst::imex::IMEX<scalar,time> {
105105
for (int i=0; i<q.size(); i++)
106106
z[i] *= c * ddx[i];
107107
fft.backward(f);
108+
109+
nf1evals++;
108110
}
109111

110112
void f2eval(Encapsulation<scalar,time> *F, Encapsulation<scalar,time> *Q, time t)

‎examples/advection_diffusion/ex1.cpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@
55
*/
66

77
#include <pfasst.hpp>
8+
#include <pfasst/sdc.hpp>
89
#include <pfasst/encap/vector.hpp>
910

1011
#include "advection_diffusion_sweeper.hpp"

‎examples/advection_diffusion/ex2.cpp‎

Lines changed: 4 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@
55
*/
66

77
#include <pfasst.hpp>
8+
#include <pfasst/mlsdc.hpp>
89
#include <pfasst/encap/vector.hpp>
910

1011
#include "advection_diffusion_sweeper.hpp"
@@ -58,11 +59,9 @@ int main(int argc, char **argv)
5859
/*
5960
* set initial conditions on each level
6061
*/
61-
for (int l=0; l<nlevs; l++) {
62-
auto* sweeper = mlsdc.get_level<AdvectionDiffusionSweeper<double,double>>(l);
63-
auto* q0 = sweeper->get_state(0);
64-
sweeper->exact(q0, 0.0);
65-
}
62+
auto* sweeper = mlsdc.get_level<AdvectionDiffusionSweeper<double,double>>(mlsdc.nlevels()-1);
63+
auto* q0 = sweeper->get_state(0);
64+
sweeper->exact(q0, 0.0);
6665

6766
/*
6867
* run mlsdc!

‎examples/advection_diffusion/ex3.cpp‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -53,7 +53,8 @@ int main(int argc, char **argv)
5353
* the 'initial' function is called once for each level to set the
5454
* intial conditions.
5555
*/
56-
auto initial = [] (EncapSweeper<double,double> *sweeper, Encapsulation<double,double> *q0) {
56+
auto initial = [] (EncapSweeper<double,double> *sweeper,
57+
Encapsulation<double,double> *q0) {
5758
auto* ad = dynamic_cast<AdvectionDiffusionSweeper<double,double>*>(sweeper);
5859
ad->exact(q0, 0.0);
5960
};

‎include/pfasst.hpp‎

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -2,8 +2,5 @@
22
#define _PFASST_HPP_
33

44
#include "pfasst/interfaces.hpp"
5-
#include "pfasst/sdc.hpp"
6-
#include "pfasst/mlsdc.hpp"
7-
#include "pfasst/pfasst.hpp"
85

96
#endif

‎include/pfasst/encap/encapsulation.hpp‎

Lines changed: 36 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -117,7 +117,7 @@ namespace pfasst {
117117

118118
virtual void advance()
119119
{
120-
this->set_state(this->get_end_state(), 0);
120+
throw NotImplementedYet("sweeper");
121121
}
122122

123123
virtual void integrate(time dt, vector<Encapsulation<scalar,time>*> dst) const
@@ -134,66 +134,67 @@ namespace pfasst {
134134

135135
virtual ~PolyInterpMixin() { }
136136

137-
virtual void interpolate(ISweeper *DST, const ISweeper *SRC, bool initial)
137+
virtual void interpolate(ISweeper *dst, const ISweeper *src, bool initial)
138138
{
139-
auto* dst = dynamic_cast<EncapSweeper<scalar,time>*>(DST);
140-
auto* src = dynamic_cast<const EncapSweeper<scalar,time>*>(SRC);
139+
auto* fine = dynamic_cast<EncapSweeper<scalar,time>*>(dst);
140+
auto* crse = dynamic_cast<const EncapSweeper<scalar,time>*>(src);
141141

142142
if (tmat.size1() == 0)
143-
tmat = pfasst::compute_interp<time>(dst->get_nodes(), src->get_nodes());
143+
tmat = pfasst::compute_interp<time>(fine->get_nodes(), crse->get_nodes());
144144

145-
int ndst = dst->get_nodes().size();
146-
int nsrc = src->get_nodes().size();
145+
int nfine = fine->get_nodes().size();
146+
int ncrse = crse->get_nodes().size();
147147

148-
auto* crse_factory = src->get_factory();
149-
auto* fine_factory = dst->get_factory();
148+
auto* crse_factory = crse->get_factory();
149+
auto* fine_factory = fine->get_factory();
150150

151-
vector<Encapsulation<scalar,time>*> fine_q(ndst), fine_tmp(nsrc);
151+
vector<Encapsulation<scalar,time>*> fine_state(nfine), fine_delta(ncrse);
152152

153-
for (int m=0; m<ndst; m++) fine_q[m] = dst->get_state(m);
154-
for (int m=0; m<nsrc; m++) fine_tmp[m] = fine_factory->create(solution);
153+
for (int m=0; m<nfine; m++) fine_state[m] = fine->get_state(m);
154+
for (int m=0; m<ncrse; m++) fine_delta[m] = fine_factory->create(solution);
155155

156156
if (initial)
157-
for (int m=1; m<ndst; m++)
158-
fine_q[m]->copy(fine_q[0]);
157+
for (int m=1; m<nfine; m++)
158+
fine_state[m]->copy(fine_state[0]);
159159

160-
auto* crse_tmp = crse_factory->create(solution);
161-
for (int m=0; m<nsrc; m++) {
162-
crse_tmp->copy(src->get_state(m));
160+
auto* crse_delta = crse_factory->create(solution);
161+
for (int m=0; m<ncrse; m++) {
162+
crse_delta->copy(crse->get_state(m));
163163
if (initial)
164-
crse_tmp->saxpy(-1.0, src->get_state(0));
164+
crse_delta->saxpy(-1.0, crse->get_state(0));
165165
else
166-
crse_tmp->saxpy(-1.0, src->get_saved_state(m));
167-
interpolate(fine_tmp[m], crse_tmp);
166+
crse_delta->saxpy(-1.0, crse->get_saved_state(m));
167+
interpolate(fine_delta[m], crse_delta);
168168
}
169-
delete crse_tmp;
169+
delete crse_delta;
170170

171-
dst->get_state(0)->mat_apply(fine_q, 1.0, tmat, fine_tmp, false);
171+
fine->get_state(0)->mat_apply(fine_state, 1.0, tmat, fine_delta, false);
172172

173-
for (int m=0; m<nsrc; m++) delete fine_tmp[m];
174-
for (int m=0; m<ndst; m++) dst->evaluate(m);
173+
for (int m=0; m<ncrse; m++) delete fine_delta[m];
174+
for (int m=0; m<nfine; m++) fine->evaluate(m);
175175
}
176176

177-
virtual void restrict(ISweeper *DST, const ISweeper *SRC)
177+
virtual void restrict(ISweeper *dst, const ISweeper *src, bool restrict_initial)
178178
{
179-
auto* dst = dynamic_cast<EncapSweeper<scalar,time>*>(DST);
180-
auto* src = dynamic_cast<const EncapSweeper<scalar,time>*>(SRC);
179+
auto* crse = dynamic_cast<EncapSweeper<scalar,time>*>(dst);
180+
auto* fine = dynamic_cast<const EncapSweeper<scalar,time>*>(src);
181181

182-
auto dnodes = dst->get_nodes();
183-
auto snodes = src->get_nodes();
182+
auto dnodes = crse->get_nodes();
183+
auto snodes = fine->get_nodes();
184184

185-
int ndst = dst->get_nodes().size();
186-
int nsrc = src->get_nodes().size();
185+
int ncrse = crse->get_nodes().size();
186+
int nfine = fine->get_nodes().size();
187187

188-
int trat = (nsrc - 1) / (ndst - 1);
188+
int trat = (nfine - 1) / (ncrse - 1);
189189

190-
for (int m=0; m<ndst; m++) {
190+
int m0 = restrict_initial ? 0 : 1;
191+
for (int m=m0; m<ncrse; m++) {
191192
if (dnodes[m] != snodes[m*trat])
192193
throw NotImplementedYet("coarse nodes must be nested");
193-
this->restrict(dst->get_state(m), src->get_state(m*trat));
194+
this->restrict(crse->get_state(m), fine->get_state(m*trat));
194195
}
195196

196-
for (int m=0; m<ndst; m++) dst->evaluate(m);
197+
for (int m=m0; m<ncrse; m++) crse->evaluate(m);
197198
}
198199

199200
virtual void fas(time dt, ISweeper *dst, const ISweeper *src)

‎include/pfasst/encap/imex.hpp‎

Lines changed: 16 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -50,6 +50,14 @@ namespace pfasst {
5050
return pQ[m];
5151
}
5252

53+
virtual void advance()
54+
{
55+
Q[0]->copy(Q[Q.size()-1]);
56+
Fe[0]->copy(Fe[Fe.size()-1]);
57+
Fi[0]->copy(Fi[Fi.size()-1]);
58+
}
59+
60+
5361
virtual void integrate(time dt, vector<Encapsulation<scalar,time>*> dst) const
5462
{
5563
auto* encap = dst[0];
@@ -118,13 +126,15 @@ namespace pfasst {
118126
delete rhs;
119127
}
120128

121-
virtual void predict(time t0, time dt)
129+
virtual void predict(time t0, time dt, bool initial)
122130
{
123131
const auto nodes = this->get_nodes();
124132
const int nnodes = nodes.size();
125133

126-
f1eval(Fe[0], Q[0], t0);
127-
f2eval(Fi[0], Q[0], t0);
134+
if (initial) {
135+
f1eval(Fe[0], Q[0], t0);
136+
f2eval(Fi[0], Q[0], t0);
137+
}
128138

129139
Encapsulation<scalar,time> *rhs = this->get_factory()->create(pfasst::encap::solution);
130140

@@ -150,9 +160,9 @@ namespace pfasst {
150160

151161
virtual void evaluate(int m)
152162
{
153-
// XXX: time
154-
f1eval(Fe[m], Q[m], 0.0);
155-
f2eval(Fi[m], Q[m], 0.0);
163+
time t = this->get_nodes()[m]; // XXX
164+
f1eval(Fe[m], Q[m], t);
165+
f2eval(Fi[m], Q[m], t);
156166
}
157167

158168

‎include/pfasst/interfaces.hpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -31,16 +31,16 @@ namespace pfasst {
3131
virtual ~ISweeper() { }
3232

3333
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
34+
virtual void predict(double t, double dt, bool initial) = 0; // XXX: this needs to be templated
3535
virtual void advance() = 0;
3636
virtual void save() { NotImplementedYet("mlsdc/pfasst"); }
3737
};
3838

3939
class ITransfer {
4040
public:
4141
// XXX: pass level iterator to these routines as well
42-
virtual void interpolate(ISweeper *dst, const ISweeper *src, bool initial) = 0;
43-
virtual void restrict(ISweeper *dst, const ISweeper *src) = 0;
42+
virtual void interpolate(ISweeper *dst, const ISweeper *src, bool initial=false) = 0;
43+
virtual void restrict(ISweeper *dst, const ISweeper *src, bool initial=false) = 0;
4444
virtual void fas(double dt, ISweeper *dst, const ISweeper *src) = 0; // XXX: this needs to be templated
4545
};
4646

0 commit comments

Comments
 (0)