Skip to content

Commit 4bb9e5b

Browse files
torbjoernktorbjoernk
committed
examples: boris: use and override void residual(dt, dst) const 
1 parent d6e12b0 commit 4bb9e5b

File tree

5 files changed

+69
-35
lines changed

5 files changed

+69
-35
lines changed

‎examples/boris/boris_sdc.cpp‎

Lines changed: 11 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,8 @@ namespace pfasst
2020
{
2121
template<typename scalar>
2222
error_map<scalar> run_boris_sdc(const size_t nsteps, const scalar dt, const size_t nnodes,
23-
const size_t nparticles, const size_t niters)
23+
const size_t nparticles, const size_t niters,
24+
const double abs_res_tol, const double rel_res_tol)
2425
{
2526
SDC<> sdc;
2627

@@ -39,6 +40,7 @@ namespace pfasst
3940

4041
sweeper->set_quadrature(quad);
4142
sweeper->set_factory(factory);
43+
sweeper->set_residual_tolerances(abs_res_tol, rel_res_tol);
4244

4345
sdc.add_level(sweeper);
4446
sdc.set_duration(0.0, nsteps*dt, dt, niters);
@@ -70,12 +72,14 @@ int main(int argc, char** argv)
7072
pfasst::log::add_custom_logger("Boris");
7173
pfasst::log::add_custom_logger("Solver");
7274

73-
const size_t nsteps = pfasst::config::get_value<size_t>("num_steps", 1);
74-
const double dt = pfasst::config::get_value<double>("delta_step", 0.015625);
75-
const size_t nnodes = pfasst::config::get_value<size_t>("num_nodes", 5);
76-
const size_t nparticles = pfasst::config::get_value<size_t>("num_particles", 1);
77-
const size_t niters = pfasst::config::get_value<size_t>("num_iter", 2);
75+
const size_t nsteps = pfasst::config::get_value<size_t>("num_steps", 1);
76+
const double dt = pfasst::config::get_value<double>("delta_step", 0.015625);
77+
const size_t nnodes = pfasst::config::get_value<size_t>("num_nodes", 5);
78+
const size_t nparticles = pfasst::config::get_value<size_t>("num_particles", 1);
79+
const size_t niters = pfasst::config::get_value<size_t>("num_iter", 2);
80+
const double abs_res_tol = pfasst::config::get_value<double>("abs_res_tol", 0.0);
81+
const double rel_res_tol = pfasst::config::get_value<double>("rel_res_tol", 0.0);
7882

79-
pfasst::examples::boris::run_boris_sdc<double>(nsteps, dt, nnodes, nparticles, niters);
83+
pfasst::examples::boris::run_boris_sdc<double>(nsteps, dt, nnodes, nparticles, niters, abs_res_tol, rel_res_tol);
8084
}
8185
#endif

‎examples/boris/boris_sweeper.hpp‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -158,7 +158,7 @@ namespace pfasst
158158
boost::format data_stream_fmt;
159159

160160
acceleration_type build_rhs(const size_t m, const bool previous = false) const;
161-
scalar compute_residual();
161+
scalar compute_residual_max();
162162
void write_center_to_file(const size_t iter, const size_t sweep, const ParticleComponent<scalar>& center,
163163
const scalar energy, const scalar drift, const scalar residual);
164164
void write_particle_cloud_to_file(const size_t iter, const size_t sweep, const shared_ptr<encap_type>& cloud,
@@ -199,7 +199,9 @@ namespace pfasst
199199
//! @{
200200
virtual void setup(bool coarse = false) override;
201201
virtual void integrate(time dt, vector<shared_ptr<Encapsulation<time>>> dst) const override;
202-
virtual void integrate(time dt, vector<shared_ptr<acceleration_type>> dst_q, vector<shared_ptr<acceleration_type>> dst_qq) const;
202+
virtual void integrate(time dt, vector<shared_ptr<acceleration_type>> dst_q,
203+
vector<shared_ptr<acceleration_type>> dst_qq) const;
204+
virtual void residual(time dt, vector<shared_ptr<Encapsulation<time>>> dst) const override;
203205
virtual void advance() override;
204206
virtual void evaluate(size_t m);
205207
virtual void predict(bool initial) override;

‎examples/boris/boris_sweeper_impl.hpp‎

Lines changed: 39 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -123,37 +123,20 @@ namespace pfasst
123123
}
124124

125125
template<typename scalar, typename time>
126-
scalar BorisSweeper<scalar, time>::compute_residual()
126+
scalar BorisSweeper<scalar, time>::compute_residual_max()
127127
{
128-
BCVLOG(8) << "computing residual";
128+
BCVLOG(8) << "computing max residual";
129129
this->log_indent->increment(8);
130-
const auto nodes = this->get_nodes();
131-
const size_t nnodes = nodes.size();
132-
time dt = this->get_controller()->get_time_step();
133130

134-
scalar max_residual = 0.0;
131+
this->residual(this->get_controller()->get_time_step(), this->residuals);
135132

136-
for (size_t m = 1; m < nnodes; ++m) {
137-
BCVLOG(8) << "for node " << m;
138-
position_type pos = cloud_component_factory<scalar>(this->particles[0]->size(), this->particles[0]->dim());
139-
velocity_type vel = cloud_component_factory<scalar>(this->particles[0]->size(), this->particles[0]->dim());
140-
for (size_t j = 0; j < nnodes; ++j) {
141-
pos += this->q_mat(m, j) * dt * this->particles[j]->velocities();
142-
vel += this->q_mat(m, j) * dt * this->build_rhs(j);
143-
}
144-
pos += this->start_particles->positions() - this->particles[m]->positions();
145-
vel += this->start_particles->velocities() - this->particles[m]->velocities();
133+
scalar max_residual = scalar(0.0);
146134

147-
for (size_t j = 0; j < pos.size(); ++j) {
148-
for (size_t d = 0; d < pos[j].size(); ++d) {
149-
max_residual = std::max(max_residual, abs(pos[j][d]));
150-
max_residual = std::max(max_residual, abs(vel[j][d]));
151-
}
152-
}
135+
for (size_t m = 1; m < this->residuals.size(); ++m) {
136+
max_residual = std::max(max_residual, this->residuals[m]->norm0());
153137
}
154138

155-
BCVLOG(8) << "=> residual: " << max_residual;
156-
this->log_indent->decrement(8);
139+
BCVLOG(8) << "=> max residual: " << max_residual;
157140
return max_residual;
158141
}
159142

@@ -412,7 +395,7 @@ namespace pfasst
412395
ErrorTuple<scalar> e_tuple;
413396
scalar e_end = this->impl_solver->energy(end, t);
414397
e_tuple.e_drift = abs(this->initial_energy - e_end);
415-
e_tuple.res = this->compute_residual();
398+
e_tuple.res = this->compute_residual_max();
416399

417400
size_t n = this->get_controller()->get_step();
418401
size_t k = this->get_controller()->get_iteration();
@@ -471,6 +454,7 @@ namespace pfasst
471454
this->energy_evals.resize(nnodes);
472455
for (size_t m = 0; m < nnodes; ++m) {
473456
this->particles.push_back(dynamic_pointer_cast<encap_type>(this->get_factory()->create(pfasst::encap::solution)));
457+
this->residuals.push_back(dynamic_pointer_cast<encap_type>(this->get_factory()->create(pfasst::encap::solution)));
474458
this->saved_particles.push_back(dynamic_pointer_cast<encap_type>(this->get_factory()->create(pfasst::encap::solution)));
475459
this->forces.push_back(cloud_component_factory<scalar>(this->particles[m]->size(), this->particles[m]->dim()));
476460
this->saved_forces.push_back(cloud_component_factory<scalar>(this->particles[m]->size(), this->particles[m]->dim()));
@@ -580,6 +564,35 @@ namespace pfasst
580564
this->log_indent->decrement(6);
581565
}
582566

567+
template<typename scalar, typename time>
568+
void BorisSweeper<scalar, time>::residual(time dt, vector<shared_ptr<Encapsulation<time>>> dst) const
569+
{
570+
BCVLOG(8) << "computing residual";
571+
const auto nodes = this->get_nodes();
572+
const size_t nnodes = nodes.size();
573+
assert(dst.size() == nnodes);
574+
575+
vector<shared_ptr<encap_type>> dst_cast(nnodes);
576+
for (size_t m = 0; m < nnodes; ++m) {
577+
dst_cast[m] = dynamic_pointer_cast<encap_type>(dst[m]);
578+
assert(dst_cast[m]);
579+
}
580+
581+
for (size_t m = 1; m < nnodes; ++m) {
582+
BCVLOG(8) << "for node " << m;
583+
zero(dst_cast[m]->positions());
584+
zero(dst_cast[m]->velocities());
585+
for (size_t j = 0; j < nnodes; ++j) {
586+
dst_cast[m]->positions() += this->q_mat(m, j) * dt * this->particles[j]->velocities();
587+
dst_cast[m]->velocities() += this->q_mat(m, j) * dt * this->build_rhs(j);
588+
}
589+
dst_cast[m]->positions() += this->start_particles->positions() - this->particles[m]->positions();
590+
dst_cast[m]->velocities() += this->start_particles->velocities() - this->particles[m]->velocities();
591+
}
592+
593+
this->log_indent->decrement(8);
594+
}
595+
583596
template<typename scalar, typename time>
584597
void BorisSweeper<scalar, time>::advance()
585598
{
@@ -805,7 +818,7 @@ namespace pfasst
805818
{
806819
UNUSED(comm); UNUSED(tag);
807820
// TODO: implement BorisSweeper::post
808-
};
821+
}
809822

810823
template<typename scalar, typename time>
811824
void BorisSweeper<scalar, time>::send(ICommunicator* comm, int tag, bool blocking)

‎examples/boris/particle_cloud.hpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -76,6 +76,8 @@ namespace pfasst
7676

7777
void distribute_around_center(const shared_ptr<Particle<precision>>& center);
7878

79+
virtual precision norm0() const;
80+
7981
virtual void log(el::base::type::ostream_t& os) const;
8082
};
8183

‎examples/boris/particle_cloud_impl.hpp‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -284,6 +284,19 @@ namespace pfasst
284284
VLOG(3) << LOG_INDENT << "center after distribute: " << this->center_of_mass();
285285
}
286286

287+
template<typename precision>
288+
precision ParticleCloud<precision>::norm0() const
289+
{
290+
precision norm = precision(0.0);
291+
for (size_t p = 0; p < this->size(); ++p) {
292+
for (size_t d = 0; d < this->dim(); ++d) {
293+
norm = std::max(norm, abs(this->positions()[p][d]));
294+
norm = std::max(norm, abs(this->velocities()[p][d]));
295+
}
296+
}
297+
return norm;
298+
}
299+
287300
template<typename precision>
288301
void ParticleCloud<precision>::log(el::base::type::ostream_t& os) const
289302
{

0 commit comments

Comments
 (0)