From 6bf3a1cae5151f97aa508ad2de7bf382738c3608 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?K=C3=BChle=2C=20Laura=20Christine=20=28lakue103=29?=
 <laura.kuehle@uni-duesseldorf.de>
Date: Fri, 10 Mar 2023 09:32:11 +0100
Subject: [PATCH] Enforced boundary condition for initial condition in exact
 solution only.

---
 Snakefile                        |  2 +-
 scripts/tcd/Equation.py          | 17 +++++++++++++----
 scripts/tcd/Initial_Condition.py |  9 ---------
 scripts/tcd/projection_utils.py  | 16 +++++++++++++---
 4 files changed, 27 insertions(+), 17 deletions(-)

diff --git a/Snakefile b/Snakefile
index d74be81..94421e2 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 62c9ba3..be8de94 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 fa69c85..e52c68d 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 202f50e..b7a2bed 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))
-- 
GitLab