Merge pull request #172 from torbjoernk/feature/fix-boris-again

Sweet!

Author: memmett

examples/boris/bindings/simple_physics_solver.cpp

@@ -101,7 +101,7 @@ namespace simple_physics_solver
 //         cout << "      dist = " << dist << endl;
         r = sqrt(dist2 + config->sigma2);
 //         cout << "      r = " << r << " (= sqrt(dist^2+" << config->sigma2 << ")" << endl;
-//         phis[i] += charges[j] / r;
+        phis[i] += charges[j] / r;
 
         r3 = r * r * r;
         for (size_t d = 0; d < DIM; ++d) {

examples/boris/boris_mlsdc.cpp

@@ -22,7 +22,8 @@ namespace pfasst
     {
       template<typename scalar>
       error_map<scalar> run_boris_sdc(const size_t nsteps, const scalar dt, const size_t nnodes,
-                                      const size_t nparticles, const size_t niters)
+                                      const size_t nparticles, const size_t niters,
+                                      const double abs_res_tol, const double rel_res_tol)
       {
         MLSDC<> controller;
 
@@ -33,20 +34,29 @@ namespace pfasst
           make_shared<bindings::WrapperSimplePhysicsSolver<double, double>>();
         bindings::setup(dynamic_pointer_cast<bindings::WrapperSimplePhysicsSolver<double, double>>(impl_solver));
 
-        auto quad1    = quadrature::quadrature_factory<double>(nnodes, quadrature::QuadratureType::GaussLobatto);
-        auto factory1 = make_shared<ParticleCloudFactory<double>>(nparticles, 3, mass, charge);
-        string data_file1 = "s" + to_string(nsteps) + "_i" + to_string(niters) + "_dt" + to_string(dt) + "_m" + to_string(nnodes) + "_p" + to_string(nparticles) + "_level1.csv";
-        auto sweeper1 = make_shared<BorisSweeper<double, double>>(impl_solver, data_file1);
-        auto transfer1 = make_shared<InjectiveTransfer<double, double>>();
+        // fine level
+        auto quad1        = quadrature::quadrature_factory<double>(nnodes,
+                                                                   quadrature::QuadratureType::GaussLobatto);
+        auto factory1     = make_shared<ParticleCloudFactory<double>>(nparticles, 3, mass, charge);
+        string data_file1 = "s" + to_string(nsteps) + "_i" + to_string(niters)
+                            + "_dt" + to_string(dt) + "_m" + to_string(nnodes)
+                            + "_p" + to_string(nparticles) + "_level1.csv";
+        auto sweeper1     = make_shared<BorisSweeper<double, double>>(impl_solver, data_file1);
+        auto transfer1    = make_shared<InjectiveTransfer<double, double>>();
         sweeper1->set_quadrature(quad1);
         sweeper1->set_factory(factory1);
+        sweeper1->set_residual_tolerances(abs_res_tol, rel_res_tol);
         controller.add_level(sweeper1, transfer1);
 
-        auto quad2    = quadrature::quadrature_factory<double>(nnodes, quadrature::QuadratureType::GaussLobatto);
-        auto factory2 = make_shared<ParticleCloudFactory<double>>(nparticles, 3, mass, charge);
-        string data_file2 = "s" + to_string(nsteps) + "_i" + to_string(niters) + "_dt" + to_string(dt) + "_m" + to_string(nnodes) + "_p" + to_string(nparticles) + "_level2.csv";
-        auto sweeper2 = make_shared<BorisSweeper<double, double>>(impl_solver, data_file2);
-        auto transfer2 = make_shared<InjectiveTransfer<double, double>>();
+        // coarse level
+        auto quad2        = quadrature::quadrature_factory<double>(nnodes,
+                                                                   quadrature::QuadratureType::GaussLobatto);
+        auto factory2     = make_shared<ParticleCloudFactory<double>>(nparticles, 3, mass, charge);
+        string data_file2 = "s" + to_string(nsteps) + "_i" + to_string(niters)
+                            + "_dt" + to_string(dt) + "_m" + to_string(nnodes)
+                            + "_p" + to_string(nparticles) + "_level2.csv";
+        auto sweeper2     = make_shared<BorisSweeper<double, double>>(impl_solver, data_file2);
+        auto transfer2    = make_shared<InjectiveTransfer<double, double>>();
         sweeper2->set_quadrature(quad2);
         sweeper2->set_factory(factory2);
         controller.add_level(sweeper2, transfer2);
@@ -60,17 +70,13 @@ namespace pfasst
         center->vel()[2] = 100;
 
         auto fine_sweeper = controller.get_finest<BorisSweeper<double, double>>();
-        shared_ptr<ParticleCloud<double>> q0 = dynamic_pointer_cast<ParticleCloud<double>>(fine_sweeper->get_start_state());
+        shared_ptr<ParticleCloud<double>> q0 = \
+          dynamic_pointer_cast<ParticleCloud<double>>(fine_sweeper->get_start_state());
         q0->distribute_around_center(center);
-        CLOG(INFO, "Boris") << "Initial Particle (fine) : " << *(dynamic_pointer_cast<ParticleCloud<double>>(fine_sweeper->get_start_state()));
+        CLOG(INFO, "Boris") << "Initial Particle (fine) : "
+                            << *(dynamic_pointer_cast<ParticleCloud<double>>(fine_sweeper->get_start_state()));
         fine_sweeper->set_initial_energy();
 
-//         auto coarse_sweeper = controller.get_coarsest<BorisSweeper<double, double>>();
-//         shared_ptr<ParticleCloud<double>> q0_coarse = dynamic_pointer_cast<ParticleCloud<double>>(coarse_sweeper->get_start_state());
-//         q0_coarse->copy(q0);
-//         LOG(INFO) << OUT::green << "Initial Particle (coarse): " << *(dynamic_pointer_cast<ParticleCloud<double>>(coarse_sweeper->get_start_state()));
-//         coarse_sweeper->set_initial_energy();
-
         controller.run();
 
         return fine_sweeper->get_errors();
@@ -93,7 +99,17 @@ int main(int argc, char** argv)
   const size_t nnodes     = pfasst::config::get_value<size_t>("num_nodes", 5);
   const size_t nparticles = pfasst::config::get_value<size_t>("num_particles", 1);
   const size_t niters     = pfasst::config::get_value<size_t>("num_iter", 2);
-
-  pfasst::examples::boris::run_boris_sdc<double>(nsteps, dt, nnodes, nparticles, niters);
+  const double abs_res_tol = pfasst::config::get_value<double>("abs_res_tol", 0.0);
+  const double rel_res_tol = pfasst::config::get_value<double>("rel_res_tol", 0.0);
+
+  LOG(INFO) << "nsteps=" << nsteps << ", "
+            << "dt=" << dt << ", "
+            << "nnodes=" << nnodes << ", "
+            << "nparticles=" << nparticles << ", "
+            << "niter=" << niters << ", "
+            << "abs res=" << abs_res_tol << ", "
+            << "rel res=" << rel_res_tol;
+  
+  pfasst::examples::boris::run_boris_sdc<double>(nsteps, dt, nnodes, nparticles, niters, abs_res_tol, rel_res_tol);
 }
 #endif

examples/boris/boris_sdc.cpp

@@ -79,6 +79,14 @@ int main(int argc, char** argv)
   const double abs_res_tol = pfasst::config::get_value<double>("abs_res_tol", 0.0);
   const double rel_res_tol = pfasst::config::get_value<double>("rel_res_tol", 0.0);
 
+  LOG(INFO) << "nsteps=" << nsteps << ", "
+            << "dt=" << dt << ", "
+            << "nnodes=" << nnodes << ", "
+            << "nparticles=" << nparticles << ", "
+            << "niter=" << niters << ", "
+            << "abs res=" << abs_res_tol << ", "
+            << "rel res=" << rel_res_tol;
+
   pfasst::examples::boris::run_boris_sdc<double>(nsteps, dt, nnodes, nparticles, niters, abs_res_tol, rel_res_tol);
 }
 #endif

examples/boris/boris_sweeper_impl.hpp

@@ -546,11 +546,14 @@ namespace pfasst
           zero(dst_q[m]);
           zero(dst_qq[m]);
           for (size_t n = 0; n < nnodes; ++n) {
+            // for positions
+            *(dst_qq[m].get()) += rhs[n] * dt * dt * this->qq_mat(m, n)
+                                  + this->start_particles->velocities() * dt * this->q_mat(m, n);
+            // for velocities
             *(dst_q[m].get()) += rhs[n] * dt * this->q_mat(m, n);
-            *(dst_qq[m].get()) += rhs[n] * dt * dt * this->qq_mat(m, n);
           }
-          BCVLOG(6) << "integral(Q)[" << m << "]:  " << *(dst_q[m].get());
           BCVLOG(6) << "integral(QQ)[" << m << "]: " << *(dst_qq[m].get());
+          BCVLOG(6) << "integral(Q)[" << m << "]:  " << *(dst_q[m].get());
         }
         this->log_indent->decrement(6);
       }
@@ -564,21 +567,35 @@ namespace pfasst
         assert(dst.size() == nnodes);
 
         vector<shared_ptr<encap_type>> dst_cast(nnodes);
+        vector<shared_ptr<acceleration_type>> q_int(nnodes), qq_int(nnodes);
         for (size_t m = 0; m < nnodes; ++m) {
           dst_cast[m] = dynamic_pointer_cast<encap_type>(dst[m]);
           assert(dst_cast[m]);
+
+          qq_int[m] = make_shared<acceleration_type>(this->start_particles->size() * this->start_particles->dim());
+          q_int[m] = make_shared<acceleration_type>(this->start_particles->size() * this->start_particles->dim());
         }
 
+        // cf. pySDC: QF(u)
+        this->integrate(dt, q_int, qq_int);
+
         for (size_t m = 1; m < nnodes; ++m) {
           BCVLOG(8) << "for node " << m;
           zero(dst_cast[m]->positions());
           zero(dst_cast[m]->velocities());
-          for (size_t j = 0; j < nnodes; ++j) {
-            dst_cast[m]->positions() += this->particles[j]->velocities() * this->q_mat(m, j) * dt;
-            dst_cast[m]->velocities() += this->build_rhs(j) * this->q_mat(m, j) * dt;
-          }
+
+          dst_cast[m]->positions() = *(qq_int[m].get());
+          dst_cast[m]->velocities() = *(q_int[m].get());
+
+          // cf. pySDC: L.u[0] - L.u[m+1] (in boris_2nd_order#compute_residual())
           dst_cast[m]->positions() += this->start_particles->positions() - this->particles[m]->positions();
           dst_cast[m]->velocities() += this->start_particles->velocities() - this->particles[m]->velocities();
+
+          // add tau correction (if available)
+          if (this->tau_q_corrections.size() > 0 && this->tau_qq_corrections.size()) {
+            dst_cast[m]->positions() += *(this->tau_qq_corrections[m].get());
+            dst_cast[m]->velocities() += *(this->tau_q_corrections[m].get());
+          }
         }
 
         this->log_indent->decrement(8);
@@ -695,6 +712,12 @@ namespace pfasst
               BCVLOG(2) << "+= tau_qq[" << m << "] (<" << this->tau_qq_corrections[m] << ">" << *(this->tau_qq_corrections[m].get()) << ")";
               this->s_integrals[m] += *(this->tau_q_corrections[m].get());
               this->ss_integrals[m] += *(this->tau_qq_corrections[m].get());
+              if (m > 0) {
+                BCVLOG(2) << "-= tau_q[" << m-1 << "]  (<" << this->tau_q_corrections[m-1]  << ">" << *(this->tau_q_corrections[m-1].get()) << ")";
+                BCVLOG(2) << "-= tau_qq[" << m-1 << "] (<" << this->tau_qq_corrections[m-1] << ">" << *(this->tau_qq_corrections[m-1].get()) << ")";
+                this->s_integrals[m] -= *(this->tau_q_corrections[m-1].get());
+                this->ss_integrals[m] -= *(this->tau_qq_corrections[m-1].get());
+              }
             }
             this->log_indent->decrement(2);
           }

include/pfasst/config.hpp

@@ -238,6 +238,7 @@ namespace pfasst
 
       options::add_option<double>("Duration", "dt", "time step size");
       options::add_option<double>("Duration", "tend", "final time of simulation");
+      options::add_option<size_t>("Duration", "num_steps", "number time steps");
       options::add_option<size_t>("Duration", "num_iter", "number of iterations");
 
       options::add_option<size_t>("Quadrature", "num_nodes", "number of quadrature nodes");