diff --git a/src/standalone/benchmark.cc b/src/standalone/benchmark.cc
index 24f7f15..1a6a919 100644
--- a/src/standalone/benchmark.cc
+++ b/src/standalone/benchmark.cc
@@ -197,7 +197,6 @@ run_benchmark(void)
                double(results[int(NTH_PERCENTILE * NUM_ITERS)]), int(NTH_PERCENTILE * 100),
                mesh_info.int_params[AC_nx]);
 
-        acStore(mesh);
         acQuit();
         acmesh_destroy(mesh);
     }
@@ -269,7 +268,6 @@ run_benchmark(void)
     write_result_to_file(wallclock * 1e3f / steps);
     #endif
 
-    acStore(mesh);
     acQuit();
     acmesh_destroy(mesh);
 
diff --git a/src/standalone/renderer.cc b/src/standalone/renderer.cc
index 57292c1..2e8945f 100644
--- a/src/standalone/renderer.cc
+++ b/src/standalone/renderer.cc
@@ -32,13 +32,13 @@
 #include <string.h> // memcpy
 
 #include "config_loader.h"
-#include "src/core/errchk.h"
-#include "src/core/math_utils.h"
 #include "model/host_forcing.h"
 #include "model/host_memory.h"
 #include "model/host_timestep.h"
 #include "model/model_reduce.h"
 #include "model/model_rk3.h"
+#include "src/core/errchk.h"
+#include "src/core/math_utils.h"
 #include "timer_hires.h"
 
 // Window
@@ -410,9 +410,10 @@ run_renderer(void)
         /* Render */
         const float timer_diff_sec = timer_diff_nsec(frame_timer) / 1e9f;
         if (timer_diff_sec >= desired_frame_time) {
-            // acStore(mesh);
             const int num_vertices = mesh->info.int_params[AC_mxy];
             const int3 dst         = (int3){0, 0, k_slice};
+            acBoundcondStep();
+            // acStore(mesh);
             acStoreWithOffset(dst, num_vertices, mesh);
             acSynchronizeStream(STREAM_ALL);
             renderer_draw(*mesh); // Bottleneck is here
@@ -422,7 +423,6 @@ run_renderer(void)
     }
     printf("Wallclock time %f s\n", double(timer_diff_nsec(wallclock) / 1e9f));
 
-    acStore(mesh);
     acQuit();
     acmesh_destroy(mesh);
 
diff --git a/src/standalone/simulation.cc b/src/standalone/simulation.cc
index ebd527e..51a9a76 100644
--- a/src/standalone/simulation.cc
+++ b/src/standalone/simulation.cc
@@ -27,13 +27,13 @@
 #include "run.h"
 
 #include "config_loader.h"
-#include "src/core/errchk.h"
-#include "src/core/math_utils.h"
 #include "model/host_forcing.h"
 #include "model/host_memory.h"
 #include "model/host_timestep.h"
 #include "model/model_reduce.h"
 #include "model/model_rk3.h"
+#include "src/core/errchk.h"
+#include "src/core/math_utils.h"
 #include "timer_hires.h"
 
 #include <string.h>
@@ -205,6 +205,7 @@ run_simulation(void)
     write_mesh_info(&mesh_info);
     print_diagnostics(0, AcReal(.0), t_step, diag_file);
 
+    acBoundcondStep();
     acStore(mesh);
     save_mesh(*mesh, 0, t_step);
 
@@ -218,13 +219,13 @@ run_simulation(void)
     /* initialize random seed: */
     srand(312256655);
 
-    //TODO_SINK. init_sink_particle()
+    // TODO_SINK. init_sink_particle()
     //  Initialize the basic variables of the sink particle to a suitable initial value.
     //  1. Location of the particle
     //  2. Mass of the particle
     //  (3. Velocity of the particle)
-    //  This at the level of Host in this case. 
-    //  acUpdate_sink_particle() will do the similar trick to the device. 
+    //  This at the level of Host in this case.
+    //  acUpdate_sink_particle() will do the similar trick to the device.
 
     /* Step the simulation */
     for (int i = 1; i < max_steps; ++i) {
@@ -236,26 +237,26 @@ run_simulation(void)
         loadForcingParamsToDevice(forcing_params);
 #endif
 
-       //TODO_SINK acUpdate_sink_particle()
-       //  Update properties of the sing particle for acIntegrate(). Essentially:
-       //  1. Location of the particle
-       //  2. Mass of the particle
-       //  (3. Velocity of the particle)
-       //  These can be used for calculating he gravitational field. 
+        // TODO_SINK acUpdate_sink_particle()
+        //  Update properties of the sing particle for acIntegrate(). Essentially:
+        //  1. Location of the particle
+        //  2. Mass of the particle
+        //  (3. Velocity of the particle)
+        //  These can be used for calculating he gravitational field.
 
         acIntegrate(dt);
-       //TODO_SINK acAdvect_sink_particle()
-       //  THIS IS OPTIONAL. We will start from unmoving particle. 
-       //  1. Calculate the equation of motion for the sink particle. 
-       //  NOTE: Might require embedding with acIntegrate(dt). 
+        // TODO_SINK acAdvect_sink_particle()
+        //  THIS IS OPTIONAL. We will start from unmoving particle.
+        //  1. Calculate the equation of motion for the sink particle.
+        //  NOTE: Might require embedding with acIntegrate(dt).
 
-       //TODO_SINK acAccrete_sink_particle()
-       //  Calculate accretion of the sink particle from the surrounding medium 
-       //  1. Transfer density into sink particle mass
-       //  2. Transfer momentum into sink particle 
-       //  (OPTIONAL: Affection the motion of the particle)
-       //  NOTE: Might require embedding with acIntegrate(dt).
-       //  This is the hardest part. Please see Lee et al. ApJ 783 (2014) for reference. 
+        // TODO_SINK acAccrete_sink_particle()
+        //  Calculate accretion of the sink particle from the surrounding medium
+        //  1. Transfer density into sink particle mass
+        //  2. Transfer momentum into sink particle
+        //  (OPTIONAL: Affection the motion of the particle)
+        //  NOTE: Might require embedding with acIntegrate(dt).
+        //  This is the hardest part. Please see Lee et al. ApJ 783 (2014) for reference.
 
         t_step += dt;
 
@@ -298,9 +299,10 @@ run_simulation(void)
                 assumed to be asynchronous, so the meshes must be also synchronized
                 before transferring the data to the CPU. Like so:
 
-                acSynchronize();
+                acBoundcondStep();
                 acStore(mesh);
             */
+            acBoundcondStep();
             acStore(mesh);
 
             save_mesh(*mesh, i, t_step);