From 16ce084490c876a1b4ce4c8a3d28a90f91b2b447 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, 31 Mar 2023 12:49:25 +0200
Subject: [PATCH] Cleaned up exact solution for Burgers class.

---
 Snakefile               |  3 +-
 scripts/tcd/Equation.py | 71 ++++++++++++++++++-----------------------
 2 files changed, 33 insertions(+), 41 deletions(-)

diff --git a/Snakefile b/Snakefile
index b25379d..6af5e2c 100644
--- a/Snakefile
+++ b/Snakefile
@@ -34,7 +34,8 @@ TODO: Ask whether Burgers has compatible ICs other than Sine and
 
 Urgent:
 TODO: Add Burgers class completely
-TODO: Clean up exact solution for Burgers
+TODO: Clean up exact solution for Burgers -> Done
+TODO: Enabale choice of equation
 TODO: Adapt code to accommodate boundary decorator -> Done
 TODO: Add subclasses of Burgers
 TODO: Ignore ghost domains by flagging all cells with it for extreme outlier
diff --git a/scripts/tcd/Equation.py b/scripts/tcd/Equation.py
index 6b170e2..c91f798 100644
--- a/scripts/tcd/Equation.py
+++ b/scripts/tcd/Equation.py
@@ -459,12 +459,8 @@ class Burgers(Equation):
             Array containing exact evaluation of a function.
 
         """
-        pass
-        num_periods = np.floor(self.wave_speed * self.final_time /
-                               mesh.interval_len)
-        #
-        # grid = np.repeat(mesh.non_ghost_cells, self._quadrature.num_nodes) + \
-        #     mesh.cell_len / 2 * np.tile(self._quadrature.nodes, mesh.num_cells)
+        grid = np.repeat(mesh.non_ghost_cells, self._quadrature.num_nodes) + \
+            mesh.cell_len / 2 * np.tile(self._quadrature.nodes, mesh.num_cells)
         #
         # # Project points into correct periodic interval
         # points = np.array([point-self.wave_speed *
@@ -479,58 +475,53 @@ class Burgers(Equation):
         # exact = np.array([self._init_cond.calculate(mesh=mesh, x=point) for
         #                   point in points])
         #
-        # grid = np.reshape(grid, (1, grid.size))
+        grid = np.reshape(grid, (1, grid.size))
         # exact = np.reshape(exact, (1, exact.size))
         #
         # return grid, exact
-        grid = []
         exact = []
 
-        for cell in range(len(mesh)):
-            eval_points = mesh[cell] + cell_len / 2 * self._quadrature.get_eval_points()
-            grid.append(eval_points)
-            eval_values = []
-            if scheme_label == 'Linear':
-                for point in range(len(eval_points)):
-                    new_entry = self._init_cond.calculate(eval_points[point] - self._wave_speed * self._final_time
-                                                        + num_periods * self._interval_len)
-                    eval_values.append(new_entry)
-                exact.append(eval_values)
-
-        grid = np.reshape(np.array(grid), (1, len(grid) * len(grid[0])))
-
-        if  scheme_label == 'Linear':
-            # u(x,t) = u(x - wavespeed*time, 0)
-            exact = np.reshape(np.array(exact), (1, len(exact) * len(exact[0])))
-
-        if scheme_label == 'Burgers' and initial_condition == 'Sine':
+        if self._init_cond.__name__ == 'Sine':
             # u(x,t) = u(x - u0*time, 0)
             exact = self._init_cond._height_adjustment + \
                     self._init_cond._stretch_factor * \
-                    self.implicit_burgers_solver(
-                grid, self._final_time)
-        elif scheme_label == 'Burgers' and initial_condition == 'DiscontinuousConstant':
+                    self.implicit_burgers_solver(grid, mesh)
+        elif self._init_cond.__name__ == 'DiscontinuousConstant':
             # u(x,t) = u(x - u0*time, 0)
-            exact = self.rarefaction_wave(grid, self._final_time)
+            exact = self.rarefaction_wave(grid)
 
-    def rarefaction_wave(self, grid_values, final_time):
+        return grid, exact
+
+    def rarefaction_wave(self, grid_values):
         uexact = 0 * grid_values
         N = np.size(grid_values)
 
         for i in range(N):
-            if grid_values[0, i] <= - final_time:
+            if grid_values[0, i] <= - self._final_time:
                 uexact[0, i] = -1
-            elif - final_time < grid_values[0, i] < final_time:
-                uexact[0, i] = grid_values[0, i] / final_time
+            elif - self._final_time < grid_values[0, i] < self._final_time:
+                uexact[0, i] = grid_values[0, i] / self._final_time
             else:
                 uexact[0, i] = 1
 
         return uexact
 
-    def implicit_burgers_solver(self, grid_values, final_time):
+    def implicit_burgers_solver(self, grid_values, mesh):
+        """
+
+        Parameters
+        ----------
+        grid_values
+        mesh
+
+        Returns
+        -------
 
-        ## Code Adapted from Hesthaven DG Code
+        Notes
+        -----
+        Code adapted from DG code by Hesthaven.
 
+        """
         sspeed = self._init_cond._height_adjustment
         uexact = np.zeros(np.size(grid_values))
 
@@ -538,9 +529,9 @@ class Burgers(Equation):
         xx = np.linspace(0, 1, 200)
         uu = np.sin(self._init_cond._factor * np.pi * xx)
 
-        te = final_time * (2 / self._interval_len)
-        xt0 = (2 / self._interval_len) * (
-                    grid_values-sspeed * final_time).reshape(
+        te = self._final_time * (2 / mesh.interval_len)
+        xt0 = (2 / mesh.interval_len) * (
+                    grid_values-sspeed * self._final_time).reshape(
             np.size(grid_values))
         for ix in range(len(xt0)):
             if xt0[ix] > 1:
@@ -574,7 +565,7 @@ class Burgers(Equation):
             y = np.sign(xt0[ix]) * un
             if abs(ay-1) < 1e-15:
                 y = 0
-            if abs(final_time) < 1e-15:
+            if abs(self._final_time) < 1e-15:
                 y = np.sin(self._init_cond._factor * np.pi * xt0[ix])
             uexact[ix] = y
 
-- 
GitLab