examples: boris: tweaking logging and data output

torbjoernk · torbjoernk · commit 10aaf71ba2ab · 2014-12-19T08:35:30.000+01:00
diff --git a/examples/boris/CMakeLists.txt b/examples/boris/CMakeLists.txt
@@ -17,7 +17,10 @@ foreach(example ${boris_examples})
         ${3rdparty_DEPENDEND_LIBS}
         simple_physics_solver
     )
-    add_dependencies(${example} ${pfasst_DEPENDEND_TARGETS} simple_physics_solver)
+    if(${pfasst_NUM_DEPENDEND_TARGETS} GREATER 0)
+        add_dependencies(${example} ${pfasst_DEPENDEND_TARGETS})
+    endif()
+    add_dependencies(${example} simple_physics_solver)
     if(pfasst_INSTALL_EXAMPLES)
         install(TARGETS ${example} RUNTIME DESTINATION bin)
     endif()
diff --git a/examples/boris/bindings/CMakeLists.txt b/examples/boris/bindings/CMakeLists.txt
@@ -1,4 +1,3 @@
-
 add_library(simple_physics_solver
     simple_physics_solver.cpp
 )
diff --git a/examples/boris/bindings/simple_physics_solver.cpp b/examples/boris/bindings/simple_physics_solver.cpp
@@ -7,6 +7,7 @@ using namespace std;
 
 #define UNUSED(expr) (void)(expr)
 
+
 namespace simple_physics_solver
 {
   SimplePhysicsSolverConfig::SimplePhysicsSolverConfig(const double omega_e, const double omega_b,
@@ -48,22 +49,21 @@ namespace simple_physics_solver
                                  double* forces)
   {
     UNUSED(t);
-    cout << "SimplePhysicsSolver::evaluate_external_e_field(t=" << t << ")" << endl;
+//     cout << "SimplePhysicsSolver::evaluate_external_e_field(t=" << t << ")" << endl;
     double pre_factor = (- config->epsilon) * (config->omega_e * config->omega_e);
     double factor = double(0.0);
     for (size_t i = 0; i < num_particles; ++i) {
-      cout << "  particle " << i << endl;
+//       cout << "  particle " << i << endl;
       factor = pre_factor / (charges[i] / masses[i]);
-      // BUG!!!
-      cout << "    position: ";
-      internal::print_vec(positions + (i * DIM));
-      cout << endl << "    factor: " << factor << endl;
+//       cout << "    position: ";
+//       internal::print_vec(positions + (i * DIM));
+//       cout << endl << "    factor: " << factor << endl;
       internal::scale_mat_mul_vec(config->external_e_field_matrix, positions + (i * DIM), factor, forces + (i * DIM));
-      cout << "    force: ";
-      internal::print_vec(forces + (i * DIM));
-      cout << endl;
+//       cout << "    force: ";
+//       internal::print_vec(forces + (i * DIM));
+//       cout << endl;
     }
-    cout << "DONE SimplePhysicsSolver::evaluate_external_e_field(t=" << t << ")" << endl;
+//     cout << "DONE SimplePhysicsSolver::evaluate_external_e_field(t=" << t << ")" << endl;
   }
 
   void evaluate_internal_e_field(const double* positions, const double* charges, const double* masses,
@@ -72,48 +72,48 @@ namespace simple_physics_solver
                                  double* exyz, double* phis)
   {
     UNUSED(masses); UNUSED(t);
-    cout << "SimplePhysicsSolver::evaluate_internal_e_field(t=" << t << ")" << endl;
+//     cout << "SimplePhysicsSolver::evaluate_internal_e_field(t=" << t << ")" << endl;
     double r = double(0.0),
            r3 = double(0.0),
            dist = double(0.0);
 
     // computing forces on particle i
     for (size_t i = 0; i < num_particles; ++i) {
-      cout << "  particle " << i << endl;
+//       cout << "  particle " << i << endl;
 
       // null result values
       phis[i] = double(0.0);
       for (size_t d = 0; d < DIM; ++d) { exyz[i * DIM + d] = double(0.0); }
 
       // effects of particle j on particle i
       for (size_t j = 0; j < num_particles; ++j) {
-        cout << "    particle " << j << endl;
+//         cout << "    particle " << j << endl;
         if (j == i) {
-          cout << "      skipping" << endl;
+//           cout << "      skipping" << endl;
           continue;
         }
         dist = double(0.0);
         for (size_t d = 0; d < DIM; ++d) {
           dist += (positions[i * DIM + d] - positions[j * DIM + d]) * (positions[i * DIM + d] - positions[j * DIM + d]);
         }
-        cout << "      dist = " << dist << endl;
+//         cout << "      dist = " << dist << endl;
         r = sqrt(dist * dist + config->sigma2);
-        cout << "      r = " << r << " (= sqrt(dist^2+" << config->sigma2 << ")" << endl;
+//         cout << "      r = " << r << " (= sqrt(dist^2+" << config->sigma2 << ")" << endl;
         phis[i] += charges[j] / r;
 
         r3 = r * r * r;
         for (size_t d = 0; d < DIM; ++d) {
           exyz[i * DIM + d] += positions[j * DIM + d] / r3 * charges[j];
-          cout << "      exyz[" << i * DIM + d << "] += " << positions[j * DIM + d] << " / " << r3 << " * " << charges[j] << "(q) => " << exyz[i * DIM + d] << endl;
+//           cout << "      exyz[" << i * DIM + d << "] += " << positions[j * DIM + d] << " / " << r3 << " * " << charges[j] << "(q) => " << exyz[i * DIM + d] << endl;
         }
       }
 
-      cout << "    exyz = ";
-      internal::print_vec(exyz+(i*DIM));
-      cout << endl
-           << "    phi_i = " << phis[i] << endl;
+//       cout << "    exyz = ";
+//       internal::print_vec(exyz+(i*DIM));
+//       cout << endl
+//            << "    phi_i = " << phis[i] << endl;
     }
-    cout << "DONE SimplePhysicsSolver::evaluate_internal_e_field(t=" << t << ")" << endl;
+//     cout << "DONE SimplePhysicsSolver::evaluate_internal_e_field(t=" << t << ")" << endl;
   }
 
 
@@ -122,22 +122,22 @@ namespace simple_physics_solver
                         const SimplePhysicsSolverConfig* config,
                         double* forces)
   {
-    cout << "SimplePhysicsSolver::evaluate_e_field(t=" << t << ")" << endl;
+//     cout << "SimplePhysicsSolver::evaluate_e_field(t=" << t << ")" << endl;
     double* external = new double[num_particles * DIM];
     double* internal = new double[num_particles * DIM];
     double* phis = new double[num_particles];
     evaluate_external_e_field(positions, charges, masses, num_particles, t, config, external);
     evaluate_internal_e_field(positions, charges, masses, num_particles, t, config, internal, phis);
-    cout << "  -> e_forces = [ " << endl;
+//     cout << "  -> e_forces = [ " << endl;
     for (size_t i = 0; i < num_particles; ++i) {
-      cout << "                  ";
+//       cout << "                  ";
       for (size_t d = 0; d < DIM; ++d) {
         forces[i * DIM + d] = external[i * DIM + d] + internal[i * DIM + d];
       }
-      internal::print_vec(forces + (i * DIM));
-      cout << endl;
+//       internal::print_vec(forces + (i * DIM));
+//       cout << endl;
     }
-    cout << "                ]" << endl;
+//     cout << "                ]" << endl;
     delete[] external;
     delete[] internal;
     delete[] phis;
@@ -178,7 +178,7 @@ namespace simple_physics_solver
                         const size_t num_particles, const double t,
                         const SimplePhysicsSolverConfig* config)
   {
-    cout << "SimplePhysicsSolver::compute_energy(t=" << t << ")" << endl;
+//     cout << "SimplePhysicsSolver::compute_energy(t=" << t << ")" << endl;
     double e_kin = double(0.0);
     double e_pot = double(0.0);
 
@@ -190,33 +190,33 @@ namespace simple_physics_solver
     evaluate_internal_e_field(positions, charges, masses, num_particles, t, config, exyz, phis);
 
     for (size_t i = 0; i < num_particles; ++i) {
-      cout << "  particle " << i << endl;
+//       cout << "  particle " << i << endl;
       // potential energy (induced by external electric field, position and internal electric field (i.e. phis))
       internal::scale_mat_mul_vec(config->external_e_field_matrix, positions + (i * DIM),
                                   (- config->epsilon * config->omega_e * config->omega_e / double(2.0) * (charges[i] / masses[i])),
                                   temp);
       e_pot += (charges[i] * phis[i]) - internal::scalar_prod(temp, positions + (i * DIM));
-      cout << "    e_pot += " << charges[i] << " * " << phis[i] << " - <";
-      internal::print_vec(temp);
-      cout << ", ";
-      internal::print_vec(positions + (i * DIM));
-      cout << "> (=" << internal::scalar_prod(temp, positions + (i * DIM)) << ")" << endl;
+//       cout << "    e_pot += " << charges[i] << " * " << phis[i] << " - <";
+//       internal::print_vec(temp);
+//       cout << ", ";
+//       internal::print_vec(positions + (i * DIM));
+//       cout << "> (=" << internal::scalar_prod(temp, positions + (i * DIM)) << ")" << endl;
 
       // kinetic energy (induced by velocity)
       v2 = internal::scalar_prod(velocities + (i * DIM), velocities + (i * DIM));
       e_kin += masses[i] / double(2.0) * v2;
-      cout << "    e_kin += " << masses[i] << "(m) / 2.0" << " * <";
-      internal::print_vec(velocities + (i * DIM));
-      cout << " , ";
-      internal::print_vec(velocities + (i * DIM));
-      cout << "> (=" << v2 << ")" << endl;
+//       cout << "    e_kin += " << masses[i] << "(m) / 2.0" << " * <";
+//       internal::print_vec(velocities + (i * DIM));
+//       cout << " , ";
+//       internal::print_vec(velocities + (i * DIM));
+//       cout << "> (=" << v2 << ")" << endl;
     }
 
     delete[] exyz;
     delete[] phis;
     delete[] temp;
-    cout << "  -> energy = " << e_kin << " + " << e_pot << endl;
-    cout << "DONE SimplePhysicsSolver::compute_energy(t=" << t << ")" << endl;
+//     cout << "  -> energy = " << e_kin << " + " << e_pot << endl;
+//     cout << "DONE SimplePhysicsSolver::compute_energy(t=" << t << ")" << endl;
     return e_kin + e_pot;
   }
 
diff --git a/examples/boris/bindings/simple_physics_solver.hpp b/examples/boris/bindings/simple_physics_solver.hpp
@@ -8,6 +8,7 @@ using namespace std;
   #define DIM 3
 #endif
 
+
 namespace simple_physics_solver
 {
   class SimplePhysicsSolverConfig
diff --git a/examples/boris/boris_sdc.cpp b/examples/boris/boris_sdc.cpp
@@ -26,15 +26,15 @@ namespace pfasst
         const double mass = 1.0;
         const double charge = 1.0;
 
-        const double epsilon = config::get_value<scalar>("epsilon", -1.0);
-
         auto quad    = quadrature::quadrature_factory<double>(nnodes, quadrature::QuadratureType::GaussLobatto);
         auto factory = make_shared<ParticleCloudFactory<double>>(nparticles, 3, mass, charge);
 
         shared_ptr<bindings::WrapperInterface<double, double>> impl_solver = \
           make_shared<bindings::WrapperSimplePhysicsSolver<double, double>>();
         bindings::setup(dynamic_pointer_cast<bindings::WrapperSimplePhysicsSolver<double, double>>(impl_solver));
-        auto sweeper = make_shared<BorisSweeper<double, double>>(impl_solver);
+
+        string data_file = "s" + to_string(nsteps) + "_i" + to_string(niters) + "_dt" + to_string(dt) + "_m" + to_string(nnodes) + "_p" + to_string(nparticles) + ".csv";
+        auto sweeper = make_shared<BorisSweeper<double, double>>(impl_solver, data_file);
 
         sweeper->set_quadrature(quad);
         sweeper->set_factory(factory);
@@ -50,9 +50,9 @@ namespace pfasst
 
         shared_ptr<ParticleCloud<double>> q0 = dynamic_pointer_cast<ParticleCloud<double>>(sweeper->get_state(0));
         q0->distribute_around_center(center);
+        LOG(INFO) << OUT::green << "Initial Particle: " << *(dynamic_pointer_cast<ParticleCloud<double>>(sweeper->get_state(0)));
 
         sweeper->set_initial_energy();
-        LOG(INFO) << OUT::green << "Initial Particle: " << *(dynamic_pointer_cast<ParticleCloud<double>>(sweeper->get_state(0)));
         sdc.run();
 
         return sweeper->get_errors();
diff --git a/examples/boris/boris_sweeper.hpp b/examples/boris/boris_sweeper.hpp
@@ -5,8 +5,9 @@
 #include <cmath>
 #include <complex>
 #include <cstdlib>
-#include <vector>
+#include <iostream>
 #include <map>
+#include <vector>
 using namespace std;
 
 #include <Eigen/Core>
@@ -137,6 +138,8 @@ namespace pfasst
           Matrix<time> qx_mat;
           Matrix<time> qt_mat;
 
+          ofstream data_stream;
+
           acceleration_type build_rhs(const size_t m, const bool previous = false) const
           {
 //             VLOG_FUNC_START("BorisSweeper") << "m=" << m << ", previous=" << boolalpha << previous;
@@ -188,13 +191,16 @@ namespace pfasst
 
         public:
           //! @{
-          BorisSweeper(shared_ptr<bindings::WrapperInterface<scalar, time>>& impl_solver)
+          BorisSweeper(shared_ptr<bindings::WrapperInterface<scalar, time>>& impl_solver,
+                       const string& data_file)
             :   impl_solver(impl_solver)
               , errors()
               , exact_updated(false)
               , f_evals(0)
+              , data_stream(data_file, ios::out | ios::trunc)
           {
             VLOG_FUNC_START("BorisSweeper");
+            assert(data_stream.is_open() && data_stream.good());
           }
 
           BorisSweeper(const BorisSweeper<scalar, time>& other) = delete;
@@ -240,11 +246,17 @@ namespace pfasst
             VLOG_FUNC_START("BorisSweeper");
             auto p0 = this->particles.front();
             this->initial_energy = this->impl_solver->energy(p0, this->get_controller()->get_time());
+            LOG(INFO) << OUT::green << "Initial Energy:" << this->initial_energy;
             for (size_t p = 0; p < p0->size(); ++p) {
-              cout << "data,0,0," << p << "," << fixed << setprecision(LOG_PRECISION)
-                   << p0->positions()[p][0] << "," << p0->positions()[p][1] << "," << p0->positions()[p][2] << ","
-                   << p0->velocities()[p][0] << "," << p0->velocities()[p][1] << "," << p0->velocities()[p][2] << endl;
+              this->data_stream << "0,0," << p << ","
+                                << fixed << setprecision(LOG_PRECISION)
+                                << p0->positions()[p][0] << "," << p0->positions()[p][1] << "," << p0->positions()[p][2] << ","
+                                << p0->velocities()[p][0] << "," << p0->velocities()[p][1] << "," << p0->velocities()[p][2] << endl;
             }
+            auto center = p0->center_of_mass();
+            this->data_stream << "0,0,-1,"
+                              << fixed << setprecision(LOG_PRECISION)
+                              << center[0] << "," << center[1] << "," << center[2] << ",0,0,0" << endl;
             VLOG_FUNC_END("BorisSweeper");
           }
           //! @}
@@ -356,13 +368,15 @@ namespace pfasst
             this->errors.insert(pair<error_index, ErrorTuple<scalar>>(nk, e_tuple));
 
             for (size_t p = 0; p < end->size(); ++p) {
-              cout << "data," << n+1 << "," << k << "," << p << "," << fixed << setprecision(LOG_PRECISION)
-                   << end->positions()[p][0] << "," << end->positions()[p][1] << "," << end->positions()[p][2] << ","
-                   << end->velocities()[p][0] << "," << end->velocities()[p][1] << "," << end->velocities()[p][2] << endl;
+              this->data_stream << n+1 << "," << k << "," << p << ","
+                                << fixed << setprecision(LOG_PRECISION)
+                                << end->positions()[p][0] << "," << end->positions()[p][1] << "," << end->positions()[p][2] << ","
+                                << end->velocities()[p][0] << "," << end->velocities()[p][1] << "," << end->velocities()[p][2] << endl;
             }
             auto center = end->center_of_mass();
-              cout << "data," << n+1 << "," << k << ",-1," << fixed << setprecision(LOG_PRECISION)
-                   << center[0] << "," << center[1] << "," << center[2] << ",0,0,0" << endl;
+            this->data_stream << n+1 << "," << k << ",-1,"
+                              << fixed << setprecision(LOG_PRECISION)
+                              << center[0] << "," << center[1] << "," << center[2] << ",0,0,0" << endl;
           }
 
           virtual error_map<scalar> get_errors() const

Original file line number	Diff line number	Diff line change
`@@ -1,4 +1,3 @@`
`1`		`-`
`2`	`1`	`add_library(simple_physics_solver`
`3`	`2`	`simple_physics_solver.cpp`
`4`	`3`	`)`
Original file line number	Diff line number	Diff line change
`@@ -8,6 +8,7 @@ using namespace std;`
`8`	`8`	`#define DIM 3`
`9`	`9`	`#endif`
`10`	`10`
	`11`	`+`
`11`	`12`	`namespace simple_physics_solver`
`12`	`13`	`{`
`13`	`14`	`class SimplePhysicsSolverConfig`