diff --git a/.gitignore b/.gitignore
index fcbb18e667efb2eb74151572917015f302a33eff..4205259e44f53e38f85702960d2499ffcf11103b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,7 +2,7 @@
 functions.dot
 functions.pdf
 **/*.vtk
-build
+build*
 dist
 pairs.egg-info
 output*
diff --git a/runtime/vtk.cpp b/runtime/vtk.cpp
index 801b0e6ca5dcca7a30d02c24ede63f91193b4689..f7d1e235e61643796123ae48accc0c5425ac9c8c 100644
--- a/runtime/vtk.cpp
+++ b/runtime/vtk.cpp
@@ -81,6 +81,10 @@ void vtk_write_data(
         out_file << "\n\n";
         out_file.close();
     }
+    else {
+        std::cerr << "vtk_write_data: Failed to open " << filename_oss.str() << std::endl;
+        exit(-1);
+    }
 }
 
 }
diff --git a/src/pairs/analysis/devices.py b/src/pairs/analysis/devices.py
index 5a15a13ea57ba34ce0f019ec0f202e040e5bd396..9fdddc2e53a8a50a6ad50338fff9857d84b30814 100644
--- a/src/pairs/analysis/devices.py
+++ b/src/pairs/analysis/devices.py
@@ -12,28 +12,18 @@ from pairs.ir.vectors import VectorOp
 class MarkCandidateLoops(Visitor):
     def __init__(self, ast=None):
         super().__init__(ast)
+        self.device_module = False
 
-    def visit_Module(self, ast_node):
-        possible_candidates = []
-        for stmt in ast_node._block.stmts:
-            if stmt is not None:
-                if isinstance(stmt, Branch):
-                    for branch_stmt in stmt.block_if.stmts:
-                        if isinstance(branch_stmt, For):
-                            possible_candidates.append(branch_stmt)
-
-                    if stmt.block_else is not None:
-                        for branch_stmt in stmt.block_else.stmts:
-                            if isinstance(branch_stmt, For):
-                                possible_candidates.append(branch_stmt)
-
-                if isinstance(stmt, For) and not stmt.not_kernel:
-                    possible_candidates.append(stmt)
-
-        for stmt in possible_candidates:
-            if not isinstance(stmt.min, Lit) or not isinstance(stmt.max, Lit):
-                stmt.mark_as_kernel_candidate()
+    def visit_For(self, ast_node):
+        if self.device_module:
+            if ast_node.not_kernel:
+                self.visit(ast_node.block)
+            else:
+                if not isinstance(ast_node.min, Lit) or not isinstance(ast_node.max, Lit):
+                    ast_node.mark_as_kernel_candidate()
 
+    def visit_Module(self, ast_node):
+        self.device_module = ast_node.run_on_device
         self.visit_children(ast_node)
 
 
diff --git a/src/pairs/code_gen/cgen.py b/src/pairs/code_gen/cgen.py
index 18d4771e891f05b69b8d1814e2402e0f92d20229..a7166be8709c99e30455318baad083d7d65e7876 100644
--- a/src/pairs/code_gen/cgen.py
+++ b/src/pairs/code_gen/cgen.py
@@ -370,7 +370,7 @@ class CGen:
 
         self.print.end()
 
-    def generate_library(self, initialize_module, create_domain_module, setup_sim_module,  do_timestep_module, reverse_comm_module, communicate_module):
+    def generate_library(self, initialize_module, create_domain_module, setup_sim_module,  do_timestep_module, reverse_comm_module, communicate_module, reset_volatiles_module):
         self.generate_interfaces()
         # Generate CUDA/CPP file with modules
         ext = ".cu" if self.target.is_gpu() else ".cpp"
@@ -398,7 +398,7 @@ class CGen:
             self.generate_kernel(kernel)
 
         for module in self.sim.modules():
-            if module.name not in ['initialize', 'create_domain', 'setup_sim', 'do_timestep', 'reverse_comm', 'communicate']:
+            if module.name not in ['initialize', 'create_domain', 'setup_sim', 'do_timestep', 'reverse_comm', 'communicate', 'reset_volatiles']:
                 self.generate_module(module)
 
         self.print.end()
@@ -495,6 +495,13 @@ class CGen:
         self.print("}")
         self.print("")
 
+        self.print("void reset_volatiles() {")
+        self.print.add_indent(4)
+        self.generate_statement(reset_volatiles_module.block)
+        self.print.add_indent(-4)
+        self.print("}")
+        self.print("")
+
         self.print("void end() {")
         self.print("    pairs::print_timers(pairs_runtime);")
         self.print("    pairs::print_stats(pairs_runtime, pobj->nlocal, pobj->nghost);")
@@ -507,6 +514,7 @@ class CGen:
 
         if self.sim.partitioner()==DomainPartitioners.BlockForest:
             self.generate_host_pairs_accessor_class()
+        # self.generate_host_pairs_accessor_class()
         
         self.print.end()
         self.generate_full_object_names = False
diff --git a/src/pairs/sim/domain_partitioning.py b/src/pairs/sim/domain_partitioning.py
index cf30351a6ca22b621b4f91cbc745d9e7a003cd13..e3d95c06421ea9fab2ac09a3f111b564fe907f17 100644
--- a/src/pairs/sim/domain_partitioning.py
+++ b/src/pairs/sim/domain_partitioning.py
@@ -17,6 +17,7 @@ class DimensionRanges:
         self.neighbor_ranks     = sim.add_static_array('neighbor_ranks', [sim.ndims() * 2], Types.Int32)
         self.pbc                = sim.add_static_array('pbc', [sim.ndims() * 2], Types.Int32)
         self.subdom             = sim.add_static_array('subdom', [sim.ndims() * 2], Types.Real)
+        self.rank               = sim.add_var('rank', Types.Int32)
 
     def min(self, dim):
         return self.subdom[dim * 2 + 0]
@@ -49,6 +50,8 @@ class DimensionRanges:
 
     def update(self):
         Call_Void(self.sim, "pairs_runtime->updateDomain", [])
+        Assign(self.sim, self.rank, Call_Int(self.sim, "pairs_runtime->getDomainPartitioner()->getRank", []))
+
 
     def ghost_particles(self, step, position, offset=0.0):
         # Particles with one of the following flags are ignored
@@ -162,7 +165,7 @@ class BlockForest:
         # Particles with one of the following flags are ignored
         flags_to_exclude = (Flags.Infinite | Flags.Global)
 
-        for r in For(self.sim, 0, self.nranks):     # for every neighbor rank
+        for r in self.step_indexes(0):     # for every neighbor rank
             for i in For(self.sim, 0, self.sim.nlocal):     # for every local particle in this rank
                 particle_flags = self.sim.particle_flags
 
diff --git a/src/pairs/sim/simulation.py b/src/pairs/sim/simulation.py
index 870231bd9f83e09c911f49c330e71b35bf53a817..044a1357dd6a4bd739cfa5639ee559e70bf83da9 100644
--- a/src/pairs/sim/simulation.py
+++ b/src/pairs/sim/simulation.py
@@ -482,12 +482,13 @@ class Simulation:
 
             timestep_procedures.append(ResetContactHistoryUsageStatus(self, self._contact_history))
 
+        # Reset volatile properties
+        if self._generate_whole_program:
+            timestep_procedures += [ResetVolatileProperties(self)]
+
         # add computational kernels
         timestep_procedures += self.functions
 
-        # Reset volatile properties
-        timestep_procedures += [ResetVolatileProperties(self)]
-
         # For whole-program-generation, add reverse_comm wherever needed in the timestep loop (eg: after computational kernels) like this:
         if self._generate_whole_program:
             timestep_procedures += [reverse_comm_module]
@@ -560,13 +561,14 @@ class Simulation:
             setup_sim_module = Module(self, name='setup_sim', block=setup_sim)
             do_timestep_module = Module(self, name='do_timestep', block=timestep.as_block())
             communicate_module = Module(self, name='communicate', block=Timestep(self, 0, comm_routine).as_block())
+            reset_volatiles_module = Module(self, name='reset_volatiles', block=Block(self, ResetVolatileProperties(self)))
 
-            modules_list = [initialize_module, create_domain_module, setup_sim_module, do_timestep_module, reverse_comm_module, communicate_module]
+            modules_list = [initialize_module, create_domain_module, setup_sim_module, do_timestep_module, reverse_comm_module, communicate_module, reset_volatiles_module]
 
             transformations = Transformations(modules_list, self._target)
             transformations.apply_all()
 
             # Generate library
-            self.code_gen.generate_library(initialize_module, create_domain_module, setup_sim_module, do_timestep_module, reverse_comm_module, communicate_module)
+            self.code_gen.generate_library(initialize_module, create_domain_module, setup_sim_module, do_timestep_module, reverse_comm_module, communicate_module, reset_volatiles_module)
 
         self.code_gen.generate_interfaces()
diff --git a/src/pairs/transformations/devices.py b/src/pairs/transformations/devices.py
index e050f9ed52c82b253970ec35c17097b0f710fc93..dceb1e3a6757f00a1c3b60de583625bd8d5938cf 100644
--- a/src/pairs/transformations/devices.py
+++ b/src/pairs/transformations/devices.py
@@ -99,62 +99,25 @@ class AddDeviceKernels(Mutator):
             self._kernel_id += 1
 
         return kernel
+    
+    def mutate_For(self, ast_node):
+        if ast_node.is_kernel_candidate():
+            kernel = self.create_kernel(ast_node.sim, ast_node.iterator, ast_node.max, ast_node.block)
+            ast_node = KernelLaunch(ast_node.sim, kernel, ast_node.iterator, ast_node.min, ast_node.max)
+
+        else:
+            ast_node.block = self.mutate(ast_node.block)
+        
+        return ast_node
 
     def mutate_Module(self, ast_node):
         if ast_node.run_on_device:
             self._module_name = ast_node.name
             self._kernel_id = 0
 
-            new_stmts = []
-            for stmt in ast_node._block.stmts:
-                if stmt is not None:
-                    if isinstance(stmt, For) and stmt.is_kernel_candidate():
-                        kernel = self.create_kernel(ast_node.sim, stmt.iterator, stmt.max, stmt.block)
-                        new_stmts.append(
-                            KernelLaunch(ast_node.sim, kernel, stmt.iterator, stmt.min, stmt.max))
-
-                    else:
-                        if isinstance(stmt, Branch):
-                            stmt = self.check_and_mutate_branch(stmt)
-
-                        new_stmts.append(stmt)
-
-            ast_node._block.stmts = new_stmts
-
         ast_node._block = self.mutate(ast_node._block)
         return ast_node
 
-    def check_and_mutate_branch(self, ast_node):
-        new_stmts = []
-        for stmt in ast_node.block_if.stmts:
-            if stmt is not None:
-                if isinstance(stmt, For) and stmt.is_kernel_candidate():
-                    kernel = self.create_kernel(ast_node.sim, stmt.iterator, stmt.max, stmt.block)
-                    new_stmts.append(
-                        KernelLaunch(ast_node.sim, kernel, stmt.iterator, stmt.min, stmt.max))
-
-                else:
-                    new_stmts.append(stmt)
-
-        ast_node.block_if.stmts = new_stmts
-
-        if ast_node.block_else is not None:
-            new_stmts = []
-            for stmt in ast_node.block_else.stmts:
-                if stmt is not None:
-                    if isinstance(stmt, For) and stmt.is_kernel_candidate():
-                        kernel = self.create_kernel(ast_node.sim, stmt.iterator, stmt.max, stmt.block)
-                        new_stmts.append(
-                            KernelLaunch(ast_node.sim, kernel, stmt.iterator, stmt.min, stmt.max))
-
-                    else:
-                        new_stmts.append(stmt)
-
-            ast_node.block_else.stmts = new_stmts
-
-        return ast_node
-
-
 class AddHostReferencesToModules(Mutator):
     def __init__(self, ast=None):
         super().__init__(ast)