diff --git a/bfps/cpp/full_code/NSVE.cpp b/bfps/cpp/full_code/NSVE.cpp
index c77220997ce9baa45bd831f775df82404024a490..0b82a781fb73f850fc39e4386d6dbb4a5e5ef5df 100644
--- a/bfps/cpp/full_code/NSVE.cpp
+++ b/bfps/cpp/full_code/NSVE.cpp
@@ -62,75 +62,19 @@ int NSVE<rnumber>::initialize(void)
 }
 
 template <typename rnumber>
-int NSVE<rnumber>::main_loop(void)
+int NSVE<rnumber>::step(void)
 {
-    clock_t time0, time1;
-    double time_difference, local_time_difference;
-    time0 = clock();
-    bool stop_code_now = false;
-    int max_iter = (this->iteration + this->niter_todo -
-                    (this->iteration % this->niter_todo));
-    for (; this->iteration < max_iter;)
-    {
-    #ifdef USE_TIMINGOUTPUT
-        const std::string loopLabel = ("code::main_start::loop-" +
-                                       std::to_string(this->iteration));
-        TIMEZONE(loopLabel.c_str());
-    #endif
-        if (this->iteration % this->niter_stat == 0)
-            this->do_stats();
-        this->fs->step(dt);
-        this->iteration = this->fs->iteration;
-        if (this->fs->iteration % this->niter_out == 0)
-        {
-            this->fs->io_checkpoint(false);
-            this->checkpoint = this->fs->checkpoint;
-            this->write_iteration();
-        }
-        if (stop_code_now)
-            break;
-        time1 = clock();
-        local_time_difference = ((
-                (unsigned int)(time1 - time0)) /
-                ((double)CLOCKS_PER_SEC));
-        time_difference = 0.0;
-        MPI_Allreduce(
-                &local_time_difference,
-                &time_difference,
-                1,
-                MPI_DOUBLE,
-                MPI_SUM,
-                MPI_COMM_WORLD);
-        if (this->myrank == 0)
-            std::cout << "iteration " << iteration << " took " << time_difference/nprocs << " seconds" << std::endl;
-        if (this->myrank == 0)
-            std::cerr << "iteration " << iteration << " took " << time_difference/nprocs << " seconds" << std::endl;
-        time0 = time1;
-    }
-    if (this->iteration % this->niter_stat == 0)
-        this->do_stats();
-    time1 = clock();
-    local_time_difference = ((
-            (unsigned int)(time1 - time0)) /
-            ((double)CLOCKS_PER_SEC));
-    time_difference = 0.0;
-    MPI_Allreduce(
-            &local_time_difference,
-            &time_difference,
-            1,
-            MPI_DOUBLE,
-            MPI_SUM,
-            MPI_COMM_WORLD);
-    if (this->myrank == 0)
-        std::cout << "iteration " << iteration << " took " << time_difference/nprocs << " seconds" << std::endl;
-    if (this->myrank == 0)
-        std::cerr << "iteration " << iteration << " took " << time_difference/nprocs << " seconds" << std::endl;
-    if (this->iteration % this->niter_out != 0)
-    {
-        this->fs->io_checkpoint(false);
-        this->checkpoint = this->fs->checkpoint;
-        this->write_iteration();
-    }
+    this->fs->step(this->dt);
+    this->iteration = this->fs->iteration;
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVE<rnumber>::write_checkpoint(void)
+{
+    this->fs->io_checkpoint(false);
+    this->checkpoint = this->fs->checkpoint;
+    this->write_iteration();
     return EXIT_SUCCESS;
 }
 
@@ -147,6 +91,8 @@ int NSVE<rnumber>::finalize(void)
 template <typename rnumber>
 int NSVE<rnumber>::do_stats()
 {
+    if (!(this->iteration % this->niter_stat == 0))
+        return EXIT_SUCCESS;
     hid_t stat_group;
     if (this->myrank == 0)
         stat_group = H5Gopen(
@@ -174,24 +120,6 @@ int NSVE<rnumber>::do_stats()
 
     if (this->myrank == 0)
         H5Gclose(stat_group);
-
-    if (myrank == 0)
-    {
-        std::string fname = (
-                std::string("stop_") +
-                std::string(simname));
-        {
-            struct stat file_buffer;
-            this->stop_code_now = (
-                    stat(fname.c_str(), &file_buffer) == 0);
-        }
-    }
-    MPI_Bcast(
-            &this->stop_code_now,
-            1,
-            MPI_C_BOOL,
-            0,
-            MPI_COMM_WORLD);
     return EXIT_SUCCESS;
 }
 
diff --git a/bfps/cpp/full_code/NSVE.hpp b/bfps/cpp/full_code/NSVE.hpp
index 27f9736585a327890c28048f73af70788ec349ea..f4f2406a7a09175e00524a38dcd29e8f917487b6 100644
--- a/bfps/cpp/full_code/NSVE.hpp
+++ b/bfps/cpp/full_code/NSVE.hpp
@@ -66,10 +66,11 @@ class NSVE: public direct_numerical_simulation
         ~NSVE(){}
 
         int initialize(void);
-        int main_loop(void);
+        int step(void);
         int finalize(void);
 
         int read_parameters(void);
+        int write_checkpoint(void);
         int do_stats(void);
 };
 
diff --git a/bfps/cpp/full_code/NSVEp.cpp b/bfps/cpp/full_code/NSVEp.cpp
index 8ce8478c635161690fc0dd2313e7e4b09267a16c..16afab62a5da1a86057bb7042737b9fd7d1ad5cc 100644
--- a/bfps/cpp/full_code/NSVEp.cpp
+++ b/bfps/cpp/full_code/NSVEp.cpp
@@ -80,96 +80,30 @@ int NSVEp<rnumber>::initialize(void)
 }
 
 template <typename rnumber>
-int NSVEp<rnumber>::main_loop(void)
+int NSVEp<rnumber>::step(void)
 {
-    clock_t time0, time1;
-    double time_difference, local_time_difference;
-    time0 = clock();
-    bool stop_code_now = false;
-    int max_iter = (this->iteration + this->niter_todo -
-                    (this->iteration % this->niter_todo));
-    for (; this->iteration < max_iter;)
-    {
-    #ifdef USE_TIMINGOUTPUT
-        const std::string loopLabel = ("code::main_start::loop-" +
-                                       std::to_string(this->iteration));
-        TIMEZONE(loopLabel.c_str());
-    #endif
-        if (this->iteration % this->niter_stat == 0)
-            this->do_stats();
-
-        this->fs->compute_velocity(fs->cvorticity);
-        this->fs->cvelocity->ift();
-        this->ps->completeLoop(dt);
+    this->fs->compute_velocity(fs->cvorticity);
+    this->fs->cvelocity->ift();
+    this->ps->completeLoop(this->dt);
+    this->fs->step(this->dt);
+    this->iteration = this->fs->iteration;
+    return EXIT_SUCCESS;
+}
 
-        this->fs->step(dt);
-        this->iteration = this->fs->iteration;
-        if (this->fs->iteration % this->niter_out == 0)
-        {
-            this->fs->io_checkpoint(false);
-            this->particles_output_writer_mpi->open_file(fs->get_current_fname());
-            this->particles_output_writer_mpi->save(
-                    this->ps->getParticlesPositions(),
-                    this->ps->getParticlesRhs(),
-                    this->ps->getParticlesIndexes(),
-                    this->ps->getLocalNbParticles(),
-                    this->fs->iteration);
-            this->particles_output_writer_mpi->close_file();
-            this->checkpoint = this->fs->checkpoint;
-            this->write_iteration();
-        }
-        if (stop_code_now)
-            break;
-        time1 = clock();
-        local_time_difference = ((
-                (unsigned int)(time1 - time0)) /
-                ((double)CLOCKS_PER_SEC));
-        time_difference = 0.0;
-        MPI_Allreduce(
-                &local_time_difference,
-                &time_difference,
-                1,
-                MPI_DOUBLE,
-                MPI_SUM,
-                MPI_COMM_WORLD);
-        if (this->myrank == 0)
-            std::cout << "iteration " << iteration << " took " << time_difference/nprocs << " seconds" << std::endl;
-        if (this->myrank == 0)
-            std::cerr << "iteration " << iteration << " took " << time_difference/nprocs << " seconds" << std::endl;
-        time0 = time1;
-    }
-    if (this->iteration % this->niter_stat == 0)
-        this->do_stats();
-    time1 = clock();
-    local_time_difference = ((
-            (unsigned int)(time1 - time0)) /
-            ((double)CLOCKS_PER_SEC));
-    time_difference = 0.0;
-    MPI_Allreduce(
-            &local_time_difference,
-            &time_difference,
-            1,
-            MPI_DOUBLE,
-            MPI_SUM,
-            MPI_COMM_WORLD);
-    if (this->myrank == 0)
-        std::cout << "iteration " << iteration << " took " << time_difference/nprocs << " seconds" << std::endl;
-    if (this->myrank == 0)
-        std::cerr << "iteration " << iteration << " took " << time_difference/nprocs << " seconds" << std::endl;
-    if (this->iteration % this->niter_out != 0)
-    {
-        this->fs->io_checkpoint(false);
-        this->particles_output_writer_mpi->open_file(fs->get_current_fname());
-        this->particles_output_writer_mpi->save(
-                this->ps->getParticlesPositions(),
-                this->ps->getParticlesRhs(),
-                this->ps->getParticlesIndexes(),
-                this->ps->getLocalNbParticles(),
-                this->fs->iteration);
-        this->particles_output_writer_mpi->close_file();
-        this->checkpoint = this->fs->checkpoint;
-        this->write_iteration();
-    }
+template <typename rnumber>
+int NSVEp<rnumber>::write_checkpoint(void)
+{
+    this->fs->io_checkpoint(false);
+    this->particles_output_writer_mpi->open_file(fs->get_current_fname());
+    this->particles_output_writer_mpi->save(
+            this->ps->getParticlesPositions(),
+            this->ps->getParticlesRhs(),
+            this->ps->getParticlesIndexes(),
+            this->ps->getLocalNbParticles(),
+            this->fs->iteration);
+    this->particles_output_writer_mpi->close_file();
+    this->checkpoint = this->fs->checkpoint;
+    this->write_iteration();
     return EXIT_SUCCESS;
 }
 
@@ -188,51 +122,41 @@ int NSVEp<rnumber>::finalize(void)
 template <typename rnumber>
 int NSVEp<rnumber>::do_stats()
 {
-    hid_t stat_group;
-    if (this->myrank == 0)
-        stat_group = H5Gopen(
-                this->stat_file,
-                "statistics",
-                H5P_DEFAULT);
-    else
-        stat_group = 0;
-    fs->compute_velocity(fs->cvorticity);
-    *tmp_vec_field = fs->cvelocity->get_cdata();
-    tmp_vec_field->compute_stats(
-            fs->kk,
-            stat_group,
-            "velocity",
-            fs->iteration / niter_stat,
-            max_velocity_estimate/sqrt(3));
-
-    *tmp_vec_field = fs->cvorticity->get_cdata();
-    tmp_vec_field->compute_stats(
-            fs->kk,
-            stat_group,
-            "vorticity",
-            fs->iteration / niter_stat,
-            max_vorticity_estimate/sqrt(3));
-
-    if (this->myrank == 0)
-        H5Gclose(stat_group);
-
-    if (myrank == 0)
+    // fluid stats go here
+    if (this->iteration % this->niter_stat == 0)
     {
-        std::string fname = (
-                std::string("stop_") +
-                std::string(simname));
-        {
-            struct stat file_buffer;
-            this->stop_code_now = (
-                    stat(fname.c_str(), &file_buffer) == 0);
-        }
+        hid_t stat_group;
+        if (this->myrank == 0)
+            stat_group = H5Gopen(
+                    this->stat_file,
+                    "statistics",
+                    H5P_DEFAULT);
+        else
+            stat_group = 0;
+        fs->compute_velocity(fs->cvorticity);
+        *tmp_vec_field = fs->cvelocity->get_cdata();
+        tmp_vec_field->compute_stats(
+                fs->kk,
+                stat_group,
+                "velocity",
+                fs->iteration / niter_stat,
+                max_velocity_estimate/sqrt(3));
+
+        *tmp_vec_field = fs->cvorticity->get_cdata();
+        tmp_vec_field->compute_stats(
+                fs->kk,
+                stat_group,
+                "vorticity",
+                fs->iteration / niter_stat,
+                max_vorticity_estimate/sqrt(3));
+
+        if (this->myrank == 0)
+            H5Gclose(stat_group);
     }
-    MPI_Bcast(
-            &this->stop_code_now,
-            1,
-            MPI_C_BOOL,
-            0,
-            MPI_COMM_WORLD);
+    // particle sampling should go here
+    //if (this->iteration % this->niter_part == 0)
+    //{
+    //}
     return EXIT_SUCCESS;
 }
 
diff --git a/bfps/cpp/full_code/NSVEp.hpp b/bfps/cpp/full_code/NSVEp.hpp
index 2b3d97ef505a59f857ce49b1ac16323007105ec3..e875be7992a4fb81fb40f564cb27b219ab3f3ba7 100644
--- a/bfps/cpp/full_code/NSVEp.hpp
+++ b/bfps/cpp/full_code/NSVEp.hpp
@@ -76,10 +76,11 @@ class NSVEp: public direct_numerical_simulation
         ~NSVEp(){}
 
         int initialize(void);
-        int main_loop(void);
+        int step(void);
         int finalize(void);
 
         int read_parameters(void);
+        int write_checkpoint(void);
         int do_stats(void);
 };
 
diff --git a/bfps/cpp/full_code/direct_numerical_simulation.cpp b/bfps/cpp/full_code/direct_numerical_simulation.cpp
index c55eaa58dff04b39f0e756d93d70f9252058021e..0a7ba43ae45fa85eb0f13eb6c848f5c147d0a876 100644
--- a/bfps/cpp/full_code/direct_numerical_simulation.cpp
+++ b/bfps/cpp/full_code/direct_numerical_simulation.cpp
@@ -1,3 +1,6 @@
+#include <cstdlib>
+#include <sys/types.h>
+#include <sys/stat.h>
 #include "direct_numerical_simulation.hpp"
 
 int grow_single_dataset(hid_t dset, int tincrement)
@@ -38,6 +41,7 @@ direct_numerical_simulation::direct_numerical_simulation(
 {
     MPI_Comm_rank(this->comm, &this->myrank);
     MPI_Comm_size(this->comm, &this->nprocs);
+    this->stop_code_now = false;
 }
 
 
@@ -130,3 +134,95 @@ int direct_numerical_simulation::write_iteration(void)
     return EXIT_SUCCESS;
 }
 
+int direct_numerical_simulation::main_loop(void)
+{
+    clock_t time0, time1;
+    double time_difference, local_time_difference;
+    time0 = clock();
+    int max_iter = (this->iteration + this->niter_todo -
+                    (this->iteration % this->niter_todo));
+    for (; this->iteration < max_iter;)
+    {
+    #ifdef USE_TIMINGOUTPUT
+        const std::string loopLabel = ("code::main_start::loop-" +
+                                       std::to_string(this->iteration));
+        TIMEZONE(loopLabel.c_str());
+    #endif
+        this->do_stats();
+
+        this->step();
+        if (this->iteration % this->niter_out == 0)
+            this->write_checkpoint();
+        this->check_stopping_condition();
+        if (this->stop_code_now)
+            break;
+        time1 = clock();
+        local_time_difference = ((
+                (unsigned int)(time1 - time0)) /
+                ((double)CLOCKS_PER_SEC));
+        time_difference = 0.0;
+        MPI_Allreduce(
+                &local_time_difference,
+                &time_difference,
+                1,
+                MPI_DOUBLE,
+                MPI_SUM,
+                MPI_COMM_WORLD);
+        if (this->myrank == 0)
+            std::cout << "iteration " << iteration <<
+                         " took " << time_difference/nprocs <<
+                         " seconds" << std::endl;
+        if (this->myrank == 0)
+            std::cerr << "iteration " << iteration <<
+                         " took " << time_difference/nprocs <<
+                         " seconds" << std::endl;
+        time0 = time1;
+    }
+    this->do_stats();
+    time1 = clock();
+    local_time_difference = ((
+            (unsigned int)(time1 - time0)) /
+            ((double)CLOCKS_PER_SEC));
+    time_difference = 0.0;
+    MPI_Allreduce(
+            &local_time_difference,
+            &time_difference,
+            1,
+            MPI_DOUBLE,
+            MPI_SUM,
+            MPI_COMM_WORLD);
+    if (this->myrank == 0)
+        std::cout << "iteration " << iteration <<
+                     " took " << time_difference/nprocs <<
+                     " seconds" << std::endl;
+    if (this->myrank == 0)
+        std::cerr << "iteration " << iteration <<
+                     " took " << time_difference/nprocs <<
+                     " seconds" << std::endl;
+    if (this->iteration % this->niter_out != 0)
+        this->write_checkpoint();
+    return EXIT_SUCCESS;
+}
+
+int direct_numerical_simulation::check_stopping_condition(void)
+{
+    if (myrank == 0)
+    {
+        std::string fname = (
+                std::string("stop_") +
+                std::string(this->simname));
+        {
+            struct stat file_buffer;
+            this->stop_code_now = (
+                    stat(fname.c_str(), &file_buffer) == 0);
+        }
+    }
+    MPI_Bcast(
+            &this->stop_code_now,
+            1,
+            MPI_C_BOOL,
+            0,
+            MPI_COMM_WORLD);
+    return EXIT_SUCCESS;
+}
+
diff --git a/bfps/cpp/full_code/direct_numerical_simulation.hpp b/bfps/cpp/full_code/direct_numerical_simulation.hpp
index 73b0168c781d36b86a9c2afc4257a03f1853eb94..0b435ea0155cdcf3a3ef4df39ce430ea1a7c1da5 100644
--- a/bfps/cpp/full_code/direct_numerical_simulation.hpp
+++ b/bfps/cpp/full_code/direct_numerical_simulation.hpp
@@ -68,13 +68,17 @@ class direct_numerical_simulation
                 const std::string &simulation_name);
         virtual ~direct_numerical_simulation(){}
 
+        virtual int write_checkpoint(void) = 0;
         virtual int initialize(void) = 0;
-        virtual int main_loop(void) = 0;
+        virtual int step(void) = 0;
+        virtual int do_stats(void) = 0;
         virtual int finalize(void) = 0;
 
+        int main_loop(void);
         int read_iteration(void);
         int write_iteration(void);
         int grow_file_datasets(void);
+        int check_stopping_condition(void);
 };
 
 #endif//DIRECT_NUMERICAL_SIMULATION_HPP