Merge pull request #152 from torbjoernk/feature/boris-linearized-cloud

Boris-(ML)SDC and MPI preparations

Author: memmett

.gitignore

@@ -6,6 +6,8 @@
 /include/pfasst/site_config.hpp
 
 .python-version
+__pycache__
+.ipynb_checkpoints
 /run
 
 # Created by http://www.gitignore.io

examples/advection_diffusion/serial_mlsdc.cpp

@@ -39,6 +39,9 @@ namespace pfasst
         size_t nnodes = config::get_value<size_t>("num_nodes", 5);
         size_t ndofs  = config::get_value<size_t>("spatial_dofs", 128);
 
+        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);
+
         /*
          * build space/time discretisation levels and add them to mlsdc
          * controller.  this loop adds the finest level first, and
@@ -53,6 +56,7 @@ namespace pfasst
 
           sweeper->set_quadrature(quad);
           sweeper->set_factory(factory);
+          sweeper->set_residual_tolerances(abs_res_tol, rel_res_tol);
 
           mlsdc.add_level(sweeper, transfer);
 

examples/boris/Boris_SDC_Plots.ipynb

@@ -1,6 +1,9 @@
 set(examples_to_install ${examples_to_install})
 
 message(STATUS "  boris")
+
+option(pfasst_BORIS_SAME_LEVELS "Use identical levels with MLSDC" ON)
+
 include_directories(
     ${3rdparty_INCLUDES}
     ${pfasst_INCLUDES}
@@ -9,6 +12,7 @@ add_subdirectory(bindings)
 
 set(boris_examples
     boris_sdc
+    boris_mlsdc
 )
 
 foreach(example ${boris_examples})
@@ -26,4 +30,19 @@ foreach(example ${boris_examples})
     endif()
 endforeach(example)
 
+if(pfasst_BORIS_SAME_LEVELS)
+    get_target_property(boris_mlsdc_definitions boris_mlsdc COMPILE_DEFINITIONS)
+    if(${boris_mlsdc_definitions})
+        list(APPEND boris_mlsdc_definitions "BORIS_SAME_LEVELS")
+    else()
+        set(boris_mlsdc_definitions "BORIS_SAME_LEVELS")
+    endif()
+    set_target_properties(boris_mlsdc PROPERTIES COMPILE_DEFINITIONS ${boris_mlsdc_definitions})
+endif()
+
+add_feature_info(Dummy-Boris-MLSDC
+    pfasst_BORIS_SAME_LEVELS
+    "use identical levels for Boris-MLSDC"
+)
+
 set(examples_to_install ${examples_to_install} ${boris_examples} PARENT_SCOPE)

examples/boris/CMakeLists.txt

@@ -36,7 +36,7 @@ namespace simple_physics_solver
                                  const size_t num_particles, const double t,
                                  const SimplePhysicsSolverConfig* config,
                                  double* exyz, double* phis);
-  
+
   void evaluate_e_field(const double* positions, const double* charges, const double* masses,
                         const size_t num_particles, const double t,
                         const SimplePhysicsSolverConfig* config,

examples/boris/bindings/simple_physics_solver.hpp

@@ -43,6 +43,9 @@ namespace pfasst
           public:
             virtual ~WrapperInterface() {}
 
+            virtual ParticleCloudComponent<scalar>
+            external_e_field_evaluate(const particle_cloud_type& particles, const time t) = 0;
+
             virtual ParticleCloudComponent<scalar>
             e_field_evaluate(const particle_cloud_type& particles, const time t) = 0;
 

examples/boris/bindings/wrapper_interface.hpp

@@ -36,7 +36,7 @@ namespace pfasst
 
           protected:
             virtual size_t vector_to_array(const vector<scalar>& vec, scalar* arr);
-            virtual size_t vector2d_to_array(const vector<vector<scalar>>& vec, scalar* arr);
+            virtual size_t vector2d_to_array(const vector<scalar>& vec, scalar* arr);
 
             virtual size_t pack_positions(const particle_cloud_type& particles, scalar* packed) override;
 
@@ -58,6 +58,7 @@ namespace pfasst
             WrapperSimplePhysicsSolver();
             virtual ~WrapperSimplePhysicsSolver();
 
+            virtual ParticleCloudComponent<scalar> external_e_field_evaluate(const particle_cloud_type& particles, const time t) override;
             virtual ParticleCloudComponent<scalar> e_field_evaluate(const particle_cloud_type& particles, const time t) override;
             virtual ParticleCloudComponent<scalar> b_field_evaluate(const particle_cloud_type& particles, const time t) override;
             virtual ParticleCloudComponent<scalar> force_evaluate(const particle_cloud_type& particles, const time t) override;

examples/boris/bindings/wrapper_simple_physics_solver.hpp

@@ -27,12 +27,12 @@ namespace pfasst
 
         template<typename scalar, typename time>
         size_t
-        WrapperSimplePhysicsSolver<scalar, time>::vector2d_to_array(const vector<vector<scalar>>& vec,
-                                                                    scalar* arr)
+        WrapperSimplePhysicsSolver<scalar, time>::vector2d_to_array(const vector<scalar>& vec, scalar* arr)
         {
-          for (size_t p = 0; p < vec.size(); ++p) {
+          const size_t npart = vec.size() / DIM;
+          for (size_t p = 0; p < npart; ++p) {
             for (size_t d = 0; d < DIM; ++d) {
-              arr[p * DIM + d] = vec[p][d];
+              arr[p * DIM + d] = vec[p * DIM + d];
             }
           }
           return vec.size() * DIM;
@@ -95,11 +95,7 @@ namespace pfasst
         WrapperSimplePhysicsSolver<scalar, time>::unpack_2d(const scalar* packed,
                                                             const size_t num_particles)
         {
-          ParticleCloudComponent<scalar> out(num_particles);
-          for (size_t p = 0; p < num_particles; ++p) {
-            out[p] = this->unpack_1d(packed + (p * DIM), DIM);
-          }
-          return out;
+          return this->unpack_1d(packed, num_particles * DIM);
         }
 
 
@@ -111,11 +107,38 @@ namespace pfasst
         WrapperSimplePhysicsSolver<scalar, time>::~WrapperSimplePhysicsSolver()
         {}
 
+        template<typename scalar, typename time>
+        ParticleCloudComponent<scalar>
+        WrapperSimplePhysicsSolver<scalar, time>::external_e_field_evaluate(const particle_cloud_type& particles, const time t)
+        {
+          CVLOG(6, "Solver") << "evaluating external E-Field at t=" << t;
+          size_t num_particles = particles->size();
+          assert(DIM == particles->dim());
+
+          scalar* packed_positions = new scalar[num_particles * DIM];
+          scalar* packed_masses = new scalar[num_particles];
+          scalar* packed_charges = new scalar[num_particles];
+          this->pack_positions(particles, packed_positions);
+          this->pack_masses(particles, packed_masses);
+          this->pack_charges(particles, packed_charges);
+
+          scalar* packed_forces = new scalar[num_particles * DIM];
+          solver::evaluate_external_e_field(packed_positions, packed_charges, packed_masses, num_particles, t,
+                                            this->config.get(), packed_forces);
+
+          delete[] packed_positions;
+          delete[] packed_masses;
+          delete[] packed_charges;
+          auto forces = this->unpack_2d(packed_forces, num_particles);
+          delete[] packed_forces;
+          return forces;
+        }
+
         template<typename scalar, typename time>
         ParticleCloudComponent<scalar>
         WrapperSimplePhysicsSolver<scalar, time>::e_field_evaluate(const particle_cloud_type& particles, const time t)
         {
-          VLOG_FUNC_START("WrapperSimplePhysicsSolver") << " time=" << t;
+          CVLOG(6, "Solver") << "evaluating complete E-Field at t=" << t;
           size_t num_particles = particles->size();
           assert(DIM == particles->dim());
 
@@ -135,15 +158,14 @@ namespace pfasst
           delete[] packed_charges;
           auto forces = this->unpack_2d(packed_forces, num_particles);
           delete[] packed_forces;
-          VLOG_FUNC_END("WrapperSimplePhysicsSolver");
           return forces;
         }
 
         template<typename scalar, typename time>
         ParticleCloudComponent<scalar>
         WrapperSimplePhysicsSolver<scalar, time>::b_field_evaluate(const particle_cloud_type& particles, const time t)
         {
-          VLOG_FUNC_START("WrapperSimplePhysicsSolver") << " time=" << t;
+          CVLOG(6, "Solver") << "evaluating B-Field at t=" << t;
           size_t num_particles = particles->size();
           assert(DIM == particles->dim());
 
@@ -163,39 +185,36 @@ namespace pfasst
           delete[] packed_charges;
           auto forces = this->unpack_2d(packed_forces, num_particles);
           delete[] packed_forces;
-          VLOG_FUNC_END("WrapperSimplePhysicsSolver");
           return forces;
         }
 
         template<typename scalar, typename time>
         ParticleCloudComponent<scalar>
         WrapperSimplePhysicsSolver<scalar, time>::b_field_vecs(const particle_cloud_type& particles, const time t)
         {
-          VLOG_FUNC_START("WrapperSimplePhysicsSolver") << " time=" << t;
           auto b_vecs = cloud_component_factory<scalar>(particles->size(), particles->dim());
           for (size_t p = 0; p < particles->size(); ++p) {
-            b_vecs[p] = this->get_b_field_vector() / particles->charges()[p] / particles->masses()[p];
+            auto bvec = this->get_b_field_vector() / particles->charges()[p] / particles->masses()[p];
+            std::copy(bvec.cbegin(), bvec.cend(), b_vecs.begin() + (p * DIM));
           }
-          VLOG_FUNC_END("WrapperSimplePhysicsSolver");
           return b_vecs;
         }
 
         template<typename scalar, typename time>
         ParticleCloudComponent<scalar>
         WrapperSimplePhysicsSolver<scalar, time>::force_evaluate(const particle_cloud_type& particles, const time t)
         {
-          VLOG_FUNC_START("WrapperSimplePhysicsSolver") << " time=" << t;
+          CVLOG(6, "Solver") << "compute total force at t=" << t;
           auto e_force = this->e_field_evaluate(particles, t);
           auto b_force = this->b_field_evaluate(particles, t);
-          VLOG_FUNC_END("WrapperSimplePhysicsSolver");
           return e_force + b_force;
         }
 
         template<typename scalar, typename time>
         scalar
         WrapperSimplePhysicsSolver<scalar, time>::energy(const particle_cloud_type& particles, const time t)
         {
-          VLOG_FUNC_START("WrapperSimplePhysicsSolver") << " time=" << t;
+          CVLOG(6, "Solver") << "computing system's total energy at t=" << t;
           size_t num_particles = particles->size();
           assert(DIM == particles->dim());
 
@@ -214,7 +233,6 @@ namespace pfasst
           delete[] packed_positions;
           delete[] packed_masses;
           delete[] packed_charges;
-          VLOG_FUNC_END("WrapperSimplePhysicsSolver")
           return energy;
         }