Reworked code to ensure termination of Burgers' equation with DiscontinuousConstant function.

Snakefile

@@ -23,6 +23,7 @@ TODO: Contemplate allowing vector input for ICs
 TODO: Discuss how wavelet details should be plotted
 
 TODO: Ask why implicit solver is always over the interval [0,1] with 200 cells
+TODO: Ask why the zero-entry of the boundary matrix is chosen for Burgers
 
 Urgent:
 TODO: Mark Burgers' equation as inviscid -> Done
@@ -31,6 +32,8 @@ TODO: Make sure only ghost cells are limited for Dirichlet boundary -> Done
 TODO: Move height adjustment and stretch factor into implicit solver
     function -> Done
 TODO: Enable choice of equation -> Done
+TODO: Ensure termination of Burgers with DiscontinuousConstant function -> Done
+TODO: Ensure termination of Burgers with Sine function
 TODO: Rework burgers_volume_integral()
 TODO: Rework/refactor implicit solver
 TODO: Add subclasses of Burgers (rarefaction and shock case)

config.yaml

@@ -44,9 +44,10 @@ Approximation:
       detector_config:
           fold_len: 16
       init_cond: 'FourPeakWave'
-    Burgers:
+    Burgers_Test:
       equation: 'Burgers'
       detector: 'Theoretical'
+      init_cond: 'DiscontinuousConstant'
 
 # Parameter for Training Data Generation
 ANN_Data:

scripts/approximate_solution.py

@@ -15,7 +15,6 @@ from tcd import Update_Scheme
 from tcd.Basis_Function import OrthonormalLegendre
 from tcd.DG_Approximation import DGScheme
 from tcd.Mesh import Mesh
-# from tcd.Equation import LinearAdvection
 
 
 def main() -> None:
@@ -35,39 +34,19 @@ def main() -> None:
 
             # Adapt number of ghost cells
             if 'stencil_len' in params['detector_config']:
-                if num_ghost_cells < params['detector_config']['stencil_len']//2:
-                    num_ghost_cells = params['detector_config']['stencil_len']//2
+                if num_ghost_cells < \
+                        params['detector_config']['stencil_len'] // 2:
+                    num_ghost_cells = \
+                        params['detector_config']['stencil_len'] // 2
                     print(f'Number of ghost cells was increased to '
-                          f'{num_ghost_cells} to accommodate the stencil length.')
+                          f'{num_ghost_cells} to accommodate the stencil '
+                          f'length.')
 
         print(params)
 
-        # Initialize mesh with two ghost cells on each side
-        mesh = Mesh(num_cells=params.pop('num_mesh_cells', 64),
-                    left_bound=params.pop('left_bound', -1),
-                    right_bound=params.pop('right_bound', 1),
-                    num_ghost_cells=num_ghost_cells,
-                    boundary_config={'type': 'periodic'})
-
         # Initialize basis
-        basis = OrthonormalLegendre(params.pop('polynomial_degree', 2))
-
-        # Initialize limiter after checking its legitimacy
-        limiter_name = params.pop('limiter', 'ModifiedMinMod')
-        limiter_config = params.pop('limiter_config', {})
-        if not hasattr(Limiter, limiter_name):
-            raise ValueError('Invalid limiter: "%s"' % limiter_name)
-        limiter = getattr(Limiter, limiter_name)(mesh=mesh,
-                                                 config=limiter_config)
-
-        # Initialize troubled cell detector after checking its legitimacy
-        detector_name = params.pop('detector')
-        detector_config = params.pop('detector_config', {})
-        if not hasattr(Troubled_Cell_Detector, detector_name):
-            raise ValueError('Invalid detector: "%s"' % detector_name)
-        detector = getattr(Troubled_Cell_Detector, detector_name)(
-            config=detector_config,
-            mesh=mesh, basis=basis)
+        polynomial_degree = params.pop('polynomial_degree', 2)
+        basis = OrthonormalLegendre(polynomial_degree)
 
         # Initialize quadrature after checking its legitimacy
         quadrature_name = params.pop('quadrature', 'Gauss')
@@ -85,12 +64,26 @@ def main() -> None:
                              % init_name)
         init_cond = getattr(Initial_Condition, init_name)(config=init_config)
 
-        # Initialize equation after checking its legitimacy
         wave_speed = params.pop('wave_speed', 1)
         final_time = params.pop('final_time', 1)
         cfl_number = params.pop('cfl_number', 0.2)
-
         equation_name = params.pop('equation', 'LinearAdvection')
+
+        # Initialize mesh with two ghost cells on each side
+        boundary_config = {'type': 'periodic'} \
+            if not (equation_name == 'Burgers' and init_name ==
+                    'DiscontinuousConstant') \
+            else {'type': 'dirichlet', 'polynomial_degree': polynomial_degree,
+                  'final_time': final_time, 'num_ghost_cells': num_ghost_cells,
+                  'left_factor': init_cond._left_factor,
+                  'right_factor': init_cond._right_factor}
+        mesh = Mesh(num_cells=params.pop('num_mesh_cells', 64),
+                    left_bound=params.pop('left_bound', -1),
+                    right_bound=params.pop('right_bound', 1),
+                    num_ghost_cells=num_ghost_cells,
+                    boundary_config=boundary_config)
+
+        # Initialize equation after checking its legitimacy
         if not hasattr(Equation, equation_name):
             raise ValueError('Invalid equation: "%s"' % equation_name)
         equation = getattr(Equation, equation_name)(
@@ -98,6 +91,23 @@ def main() -> None:
             init_cond=init_cond, quadrature=quadrature, mesh=mesh,
             wave_speed=wave_speed)
 
+        # Initialize limiter after checking its legitimacy
+        limiter_name = params.pop('limiter', 'ModifiedMinMod')
+        limiter_config = params.pop('limiter_config', {})
+        if not hasattr(Limiter, limiter_name):
+            raise ValueError('Invalid limiter: "%s"' % limiter_name)
+        limiter = getattr(Limiter, limiter_name)(mesh=mesh,
+                                                 config=limiter_config)
+
+        # Initialize troubled cell detector after checking its legitimacy
+        detector_name = params.pop('detector')
+        detector_config = params.pop('detector_config', {})
+        if not hasattr(Troubled_Cell_Detector, detector_name):
+            raise ValueError('Invalid detector: "%s"' % detector_name)
+        detector = getattr(Troubled_Cell_Detector, detector_name)(
+            config=detector_config,
+            mesh=mesh, basis=basis)
+
         # Initialize update scheme after checking its legitimacy
         scheme_name = params.pop('update_scheme', 'SSPRK3')
         if not hasattr(Update_Scheme, scheme_name):

scripts/tcd/Boundary_Condition.py

@@ -67,8 +67,10 @@ def dirichlet_boundary(projection: ndarray, config: dict) -> ndarray:
         projection_values_left[0] = np.sqrt(2) * left_factor
         projection_values_right[0] = np.sqrt(2) * right_factor
 
-    projection[:, :num_ghost_cells] = np.repeat(projection_values_left)
-    projection[:, -num_ghost_cells:] = np.repeat(projection_values_right)
+    projection[:, :num_ghost_cells] = np.repeat(
+        projection_values_left, num_ghost_cells).reshape(-1, num_ghost_cells)
+    projection[:, -num_ghost_cells:] = np.repeat(
+        projection_values_right, num_ghost_cells).reshape(-1, num_ghost_cells)
     # projection[:, 0] = projection_values_left
     # projection[:, 1] = projection_values_left
     # projection[:, -2] = projection_values_right

scripts/tcd/Equation.py

@@ -607,8 +607,8 @@ class Burgers(Equation):
         # Set ghost cells to respective value
         # (Periodic, Updated to Dirichlet in enforce_boundary_condition
         # for DiscontinuousConstant problem)
-        # right_hand_side[0] = right_hand_side[self._mesh.num_cells]
-        # right_hand_side.append(right_hand_side[1])
+        right_hand_side[0] = right_hand_side[self._mesh.num_cells]
+        right_hand_side.append(right_hand_side[1])
 
         return np.transpose(right_hand_side)
         # right_hand_side = np.zeros_like(projection)
@@ -642,13 +642,12 @@ class Burgers(Equation):
         for cell in range(self._mesh.num_cells):
             new_row = []
             # Approximation u at Gauss-Legendre points
-            approx_eval = [sum(projection[degree][cell + 1]
-                               * basis_vector[degree].subs(
-                x, self._quadrature.nodes[point])
-                               for degree in range(
-                self._basis.polynomial_degree + 1))
-                           for point in range(self._quadrature.nodes)]
-            approx_eval = np.array(approx_eval)
+            approx_eval = np.array(
+                [sum(projection[degree][cell + 1] *
+                     basis_vector[degree].subs(x,
+                                               self._quadrature.nodes[point])
+                 for degree in range(self._basis.polynomial_degree + 1))
+                 for point in range(self._quadrature.num_nodes)])
 
             # print("u", approx_eval)
             # Quadrature Evaluation
@@ -668,9 +667,9 @@ class Burgers(Equation):
         # Set ghost cells to respective value
         # (Periodic, Updated to Dirichlet in enforce_boundary_condition
         # for DiscontinuousConstant problem)
-        # volume_integral_matrix[0] = volume_integral_matrix[
-        # self._mesh.num_cells]
-        # volume_integral_matrix.append(volume_integral_matrix[1])
+        volume_integral_matrix[0] = volume_integral_matrix[
+            self._mesh.num_cells]
+        volume_integral_matrix.append(volume_integral_matrix[1])
 
         return np.transpose(np.array(volume_integral_matrix))
 
@@ -726,8 +725,8 @@ class Burgers(Equation):
         # Set Ghost cells to respective value
         # (Periodic, Updated to Dirichlet in enforce_boundary_condition
         # for DiscontinuousConstant problem)
-        # boundary_matrix[0] = boundary_matrix[self._mesh.num_cells]
-        # boundary_matrix.append(boundary_matrix[1])
+        boundary_matrix[0] = boundary_matrix[self._mesh.num_cells]
+        boundary_matrix.append(boundary_matrix[1])
 
         boundary_matrix = np.transpose(np.array(boundary_matrix))