diff --git a/Snakefile b/Snakefile
index d74be814f4ef7cf1585c5fa7edd6ef6346a8c55f..94421e233f25667e0d4de4849b3dbde795cc3717 100644
--- a/Snakefile
+++ b/Snakefile
@@ -27,7 +27,7 @@ TODO: Refactor Boxplot class -> Done
 TODO: Enforce periodic boundary condition for projection with decorator -> Done
 TODO: Enforce Boxplot bounds with decorator -> Done
 TODO: Enforce Boxplot folds with decorator -> Done
-TODO: Enforce boundary for initial condition in exact solution only
+TODO: Enforce boundary for initial condition in exact solution only -> Done
 TODO: Adapt number of ghost cells based on ANN stencil
 TODO: Ensure exact solution is calculated in Equation class
 TODO: Extract objects from UpdateScheme
diff --git a/scripts/tcd/Equation.py b/scripts/tcd/Equation.py
index 62c9ba32463770813780f1db504264501b2e7b3e..be8de94e8fe269a92dd847b94d8b2b60da4e8ca7 100644
--- a/scripts/tcd/Equation.py
+++ b/scripts/tcd/Equation.py
@@ -254,10 +254,19 @@ class LinearAdvection(Equation):
 
         grid = np.repeat(mesh.non_ghost_cells, self._quadrature.num_nodes) + \
             mesh.cell_len / 2 * np.tile(self._quadrature.nodes, mesh.num_cells)
-        exact = np.array([self._init_cond.calculate(
-            mesh=mesh, x=point-self.wave_speed * self.final_time+num_periods *
-                         mesh.interval_len)
-            for point in grid])
+
+        # Project points into correct periodic interval
+        points = np.array([point-self.wave_speed *
+                           self.final_time+num_periods * mesh.interval_len
+                           for point in grid])
+        left_bound, right_bound = self._mesh.bounds
+        while np.any(points < left_bound):
+            points[points < left_bound] += self._mesh.interval_len
+        while np.any(points) > right_bound:
+            points[points > right_bound] -= self._mesh.interval_len
+
+        exact = np.array([self._init_cond.calculate(mesh=mesh, x=point) for
+                          point in points])
 
         grid = np.reshape(grid, (1, grid.size))
         exact = np.reshape(exact, (1, exact.size))
diff --git a/scripts/tcd/Initial_Condition.py b/scripts/tcd/Initial_Condition.py
index fa69c85d43767edb8828a24c6a0d23bb2babd266..e52c68da676e737c9d6817171832998741a983a8 100644
--- a/scripts/tcd/Initial_Condition.py
+++ b/scripts/tcd/Initial_Condition.py
@@ -88,9 +88,6 @@ class InitialCondition(ABC):
     def calculate(self, x, mesh):
         """Evaluates function at given x-value.
 
-        In evaluation (not training) mode, the x-value is projected
-        into its periodic interval before evaluating the function.
-
         Parameters
         ----------
         x : float
@@ -104,12 +101,6 @@ class InitialCondition(ABC):
             Value of function evaluates at x-value.
 
         """
-        if mesh.mode == 'evaluation':
-            left_bound, right_bound = mesh.bounds
-            while x < left_bound:
-                x += mesh.interval_len
-            while x > right_bound:
-                x -= mesh.interval_len
         return self._get_point(x, mesh)
 
     @abstractmethod
diff --git a/scripts/tcd/projection_utils.py b/scripts/tcd/projection_utils.py
index 202f50e199d34a6252568272d44d07670c196999..b7a2bed4462a36bb494110abdefcf10345cdfaf1 100644
--- a/scripts/tcd/projection_utils.py
+++ b/scripts/tcd/projection_utils.py
@@ -83,9 +83,19 @@ def calculate_exact_solution(
 
     grid = np.repeat(mesh.non_ghost_cells, quadrature.num_nodes) + \
         mesh.cell_len/2 * np.tile(quadrature.nodes, mesh.num_cells)
-    exact = np.array([init_cond.calculate(
-        mesh=mesh, x=point-wave_speed*final_time+num_periods*mesh.interval_len)
-        for point in grid])
+
+    # Project points into correct periodic interval
+    points = np.array([point-wave_speed *
+                       final_time+num_periods * mesh.interval_len
+                       for point in grid])
+    left_bound, right_bound = mesh.bounds
+    while np.any(points < left_bound):
+        points[points < left_bound] += mesh.interval_len
+    while np.any(points) > right_bound:
+        points[points > right_bound] -= mesh.interval_len
+
+    exact = np.array([init_cond.calculate(mesh=mesh, x=point) for
+                      point in points])
 
     grid = np.reshape(grid, (1, grid.size))
     exact = np.reshape(exact, (1, exact.size))