Skip to content

Commit cd7a1b5

Browse files
memmettmemmett
committed
Merge pull request #61 from memmett/feature/scalar2
Add DRs scalar example.
2 parents d5494a0 + ed178fb commit cd7a1b5

File tree

4 files changed

+180
-0
lines changed

4 files changed

+180
-0
lines changed

‎examples/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,7 @@
11
set(all_example_runs)
22

33
add_subdirectory(advection_diffusion)
4+
add_subdirectory(scalar)
45

56
add_custom_target(run_example_all
67
DEPENDS ${all_example_runs}

‎examples/scalar/CMakeLists.txt‎

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,16 @@
1+
message(STATUS " scalar")
2+
include_directories(
3+
${pfasst_INCLUDES}
4+
${3rdparty_INCLUDES}
5+
)
6+
7+
set(scalar_examples
8+
scalar_sdc
9+
)
10+
11+
foreach(example ${scalar_examples})
12+
add_executable(${example} ${CMAKE_CURRENT_SOURCE_DIR}/${example}.cpp)
13+
set_target_properties(${example}
14+
PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/examples/scalar
15+
)
16+
endforeach(example)

‎examples/scalar/scalar_sdc.cpp‎

Lines changed: 37 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,37 @@
1+
#include<complex>
2+
3+
#include<pfasst.hpp>
4+
#include<pfasst/sdc.hpp>
5+
#include<pfasst/encap/vector.hpp>
6+
7+
#include "scalar_sweeper.hpp"
8+
9+
int main(int, char**)
10+
{
11+
pfasst::SDC<> sdc;
12+
13+
const size_t nsteps = 2;
14+
const double dt = 1.0;
15+
const size_t nnodes = 4;
16+
const size_t niters = 6;
17+
const complex<double> lambda = complex<double>(-1.0, 1.0);
18+
const complex<double> y0 = complex<double>(1.0, 0.0);
19+
20+
// so far, only lobatto nodes appear to be working
21+
auto nodes = pfasst::compute_nodes(nnodes, "gauss-lobatto");
22+
auto factory = make_shared<pfasst::encap::VectorFactory<complex<double>>>(1);
23+
auto sweeper = make_shared<ScalarSweeper<>>(lambda, y0);
24+
25+
sweeper->set_nodes(nodes);
26+
sweeper->set_factory(factory);
27+
28+
sdc.add_level(sweeper);
29+
sdc.set_duration(0.0, dt * nsteps, dt, niters);
30+
sdc.setup();
31+
32+
auto q0 = sweeper->get_state(0);
33+
sweeper->exact(q0, 0.0);
34+
35+
sdc.run();
36+
37+
}

‎examples/scalar/scalar_sweeper.hpp‎

Lines changed: 126 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,126 @@
1+
/*
2+
* Sweeper for scalar test equation
3+
*/
4+
5+
#ifndef _SCALAR_SWEEPER_HPP_
6+
#define _SCALAR_SWEEPER_HPP_
7+
8+
#include <complex>
9+
#include <vector>
10+
#include <pfasst/encap/imex_sweeper.hpp>
11+
#include <pfasst/encap/vector.hpp>
12+
13+
using namespace std;
14+
15+
template<typename time = pfasst::time_precision>
16+
class ScalarSweeper : public pfasst::encap::IMEXSweeper<time>
17+
{
18+
19+
typedef pfasst::encap::Encapsulation<time> Encapsulation;
20+
typedef pfasst::encap::VectorEncapsulation<complex<double>> CVectorT;
21+
22+
public:
23+
24+
ScalarSweeper(complex<double> lambda, complex<double> y0)
25+
{
26+
this->_lambda = lambda;
27+
this->_y0 = y0;
28+
this->_nf1eval = 0;
29+
this->_nf2eval = 0;
30+
this->_nf2comp = 0;
31+
}
32+
33+
virtual ~ScalarSweeper()
34+
{
35+
cout << "Number of calls to f1eval: " << this->_nf1eval << endl;
36+
cout << "Number of calls to f2eval: " << this->_nf2eval << endl;
37+
cout << "Number of calls to f2comp: " << this->_nf2comp << endl;
38+
}
39+
40+
void echo_error(time t)
41+
{
42+
auto& qend = pfasst::encap::as_vector<complex<double>,time>(this->get_end_state());
43+
44+
CVectorT qex(qend.size());
45+
46+
this->exact(qex, t);
47+
double max_err = abs(qend[0] - qex[0]) / abs(qex[0]);
48+
cout << "err: " << scientific << max_err << endl;
49+
}
50+
51+
void predict(bool initial) override
52+
{
53+
pfasst::encap::IMEXSweeper<time>::predict(initial);
54+
time t = this->get_controller()->get_time();
55+
time dt = this->get_controller()->get_time_step();
56+
this->echo_error(t + dt);
57+
}
58+
59+
void sweep() override
60+
{
61+
pfasst::encap::IMEXSweeper<time>::sweep();
62+
time t = this->get_controller()->get_time();
63+
time dt = this->get_controller()->get_time_step();
64+
this->echo_error(t + dt);
65+
}
66+
67+
void exact(CVectorT& q, time t)
68+
{
69+
q[0] = _y0 * exp(_lambda * t);
70+
}
71+
72+
void exact(shared_ptr<Encapsulation> q_encap, time t)
73+
{
74+
auto& q = pfasst::encap::as_vector<complex<double>,time>(q_encap);
75+
exact(q, t);
76+
}
77+
78+
void f_expl_eval(shared_ptr<Encapsulation> f_encap,
79+
shared_ptr<Encapsulation> q_encap, time t) override
80+
{
81+
UNUSED(t);
82+
auto& f = pfasst::encap::as_vector<complex<double>,time>(f_encap);
83+
auto& q = pfasst::encap::as_vector<complex<double>,time>(q_encap);
84+
85+
// f1 = multiply with imaginary part of lambda
86+
f[0] = i_complex * imag(this->_lambda) * q[0];
87+
this->_nf1eval++;
88+
}
89+
90+
void f_impl_eval(shared_ptr<Encapsulation> f_encap,
91+
shared_ptr<Encapsulation> q_encap, time t) override
92+
{
93+
UNUSED(t);
94+
auto& f = pfasst::encap::as_vector<complex<double>,time>(f_encap);
95+
auto& q = pfasst::encap::as_vector<complex<double>,time>(q_encap);
96+
97+
// f2 = multiply with real part of lambda
98+
f[0] = real(this->_lambda) * q[0];
99+
this->_nf2eval++;
100+
}
101+
102+
void impl_solve(shared_ptr<Encapsulation> f_encap,
103+
shared_ptr<Encapsulation> q_encap, time t, time dt,
104+
shared_ptr<Encapsulation> rhs_encap) override
105+
{
106+
UNUSED(t);
107+
auto& f = pfasst::encap::as_vector<complex<double>,time>(f_encap);
108+
auto& q = pfasst::encap::as_vector<complex<double>,time>(q_encap);
109+
auto& rhs = pfasst::encap::as_vector<complex<double>,time>(rhs_encap);
110+
111+
// invert f2=multiply with inverse of real part of lambda
112+
double inv = 1 / (1 - double(dt) * real(this->_lambda));
113+
q[0] = inv * rhs[0];
114+
f[0] = real(this->_lambda) * q[0];
115+
this->_nf2comp++;
116+
}
117+
118+
119+
private:
120+
121+
complex<double> _lambda, _y0;
122+
int _nf1eval, _nf2eval, _nf2comp;
123+
const complex<double> i_complex = complex<double>(0, 1);
124+
125+
};
126+
#endif

0 commit comments

Comments
 (0)