Added 'update_time_step()' for Burgers class.

Snakefile

@@ -21,9 +21,18 @@ TODO: Discuss name for quadrature mesh (now: grid)
 TODO: Contemplate using lambdify for basis
 TODO: Contemplate allowing vector input for ICs
 TODO: Discuss how wavelet details should be plotted
+TODO: Ask whether we are dealing with inviscid Burgers only
+TODO: Ask whether it is helpful to enforce boundary at every SSPRK3 step
+TODO: Ask why we limit the initial time step
+TODO: ASk whether cfl number should be absolute value
 
 Urgent:
-TODO: Add Burger class
+TODO: Add Burgers class completely
+TODO: Add 'update_time_step()' to Burgers -> Done
+TODO: Add 'solve_exactly()' to Burgers
+TODO: Add 'update_right_hand_side()' to Burgers
+TODO: Add initialization to Burgers
+TODO: Add Dirichlet boundary
 TODO: Use cfl_number for updating, not just time (equation-related?)
 
 Critical, but not urgent:

scripts/tcd/Equation.py

@@ -131,12 +131,14 @@ class Equation(ABC):
                 }
 
     @abstractmethod
-    def update_time_step(self, current_time: float,
+    def update_time_step(self, projection: ndarray, current_time: float,
                          time_step: float) -> Tuple[float, float]:
         """Update time step.
 
         Parameters
         ----------
+        projection : ndarray
+            Current projection during approximation.
         current_time : float
             Current time during approximation.
         time_step : float
@@ -242,12 +244,14 @@ class LinearAdvection(Equation):
                                      mesh=self._mesh, basis=self._basis,
                                      quadrature=self._quadrature)
 
-    def update_time_step(self, current_time: float,
+    def update_time_step(self, projection: ndarray, current_time: float,
                          time_step: float) -> Tuple[float, float]:
         """Update time step.
 
         Parameters
         ----------
+        projection : ndarray
+            Current projection during approximation.
         current_time : float
             Current time during approximation.
         time_step : float
@@ -261,9 +265,12 @@ class LinearAdvection(Equation):
 
         """
         cfl_number = self._cfl_number
+
+        # Adjust for final time-step
         if current_time+time_step > self.final_time:
             time_step = self.final_time-current_time
             cfl_number = self.wave_speed * time_step / self._mesh.cell_len
+
         return cfl_number, time_step
 
     def solve_exactly(self, mesh: Mesh) -> Tuple[ndarray, ndarray]:
@@ -341,3 +348,171 @@ class LinearAdvection(Equation):
                                 -self._mesh.num_ghost_cells])
 
         return right_hand_side
+
+
+class Burgers(Equation):
+    """Class for Burgers' equation.
+
+    Attributes
+    ----------
+    volume_integral_matrix : ndarray
+        Volume integral matrix.
+    flux_matrix : ndarray
+        Flux matrix.
+
+    Methods
+    -------
+    update_right_hand_side(projection)
+        Update right-hand side.
+
+    Notes
+    -----
+    .. math:: u_t + u*u_x = 0
+
+    """
+
+    def _reset(self) -> None:
+        """Reset instance variables."""
+        pass
+        # matrix_shape = (self._basis.polynomial_degree+1,
+        #                 self._basis.polynomial_degree+1)
+        # root_vector = np.array([np.sqrt(degree+0.5) for degree in range(
+        #     self._basis.polynomial_degree+1)])
+        # degree_matrix = np.matmul(root_vector[:, np.newaxis],
+        #                           root_vector[:, np.newaxis].T)
+        #
+        # # Set volume integral matrix
+        # matrix = np.ones(matrix_shape)
+        # if self._wave_speed > 0:
+        #     matrix[np.fromfunction(
+        #         lambda i, j: (j >= i) | ((i+j) % 2 == 0), matrix_shape)] = -1.0
+        # else:
+        #     matrix[np.fromfunction(
+        #         lambda i, j: (j > i) & ((i+j) % 2 == 1), matrix_shape)] = -1.0
+        # self._volume_integral_matrix = matrix * degree_matrix
+        #
+        # # Set flux matrix
+        # matrix = np.fromfunction(lambda i, j: (-1.0)**i, matrix_shape) \
+        #     if self._wave_speed > 0 \
+        #     else np.fromfunction(lambda i, j: (-1.0)**j, matrix_shape)
+        # self._flux_matrix = matrix * degree_matrix
+
+    @enforce_boundary()
+    def _initialize_projection(self) -> ndarray:
+        """Initialize projection."""
+        return do_initial_projection(init_cond=self._init_cond,
+                                     mesh=self._mesh, basis=self._basis,
+                                     quadrature=self._quadrature)
+
+    def update_time_step(self, projection: ndarray, current_time: float,
+                         time_step: float) -> Tuple[float, float]:
+        """Update time step.
+
+        Adapt CFL number to ensure right-hand side is multiplied with dt/dx.
+
+        Parameters
+        ----------
+        projection : ndarray
+            Current projection during approximation.
+        current_time : float
+            Current time during approximation.
+        time_step : float
+            Length of time step during approximation.
+        Returns
+        -------
+        cfl_number : float
+            Updated CFL number to ensure stability.
+        time_step : float
+            Updated time step for approximation.
+
+        """
+        cfl_number = self._cfl_number
+
+        max_velocity = max(abs(projection[0, :] * np.sqrt(0.5)))
+        time_step = cfl_number * self._mesh.cell_len / max_velocity
+        cfl_number = time_step / self._mesh.cell_len
+
+        # Adjust for final time-step
+        if current_time+time_step > self._final_time:
+            time_step = self._final_time-current_time
+            cfl_number = time_step / self._mesh.cell_len
+
+        return cfl_number, time_step
+
+    def solve_exactly(self, mesh: Mesh) -> Tuple[ndarray, ndarray]:
+        """Calculate exact solution.
+
+        Parameters
+        ----------
+        mesh : Mesh
+            Mesh for evaluation.
+
+        Returns
+        -------
+        grid : ndarray
+            Array containing evaluation grid for a function.
+        exact : ndarray
+            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)
+        #
+        # # 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 = 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([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))
+        #
+        # return grid, exact
+
+    @enforce_boundary()
+    def update_right_hand_side(self, projection: ndarray) -> ndarray:
+        """Update right-hand side.
+
+        Parameter
+        ---------
+        projection : ndarray
+            Matrix of projection for each polynomial degree.
+
+        Returns
+        -------
+        ndarray
+            Matrix of right-hand side.
+
+        """
+        pass
+        # right_hand_side = np.zeros_like(projection)
+        # if self._wave_speed > 0:
+        #     right_hand_side[:, self._mesh.num_ghost_cells:
+        #                     -self._mesh.num_ghost_cells] = \
+        #         2 * (self._flux_matrix @
+        #              projection[:, self._mesh.num_ghost_cells-1:
+        #                         -self._mesh.num_ghost_cells-1] +
+        #              self._volume_integral_matrix @
+        #              projection[:, self._mesh.num_ghost_cells:
+        #                         -self._mesh.num_ghost_cells])
+        # else:
+        #     right_hand_side[:, self._mesh.num_ghost_cells:
+        #                     -self._mesh.num_ghost_cells] = \
+        #         2 * (self._flux_matrix @
+        #              projection[:, self._mesh.num_ghost_cells+1:] +
+        #              self._volume_integral_matrix @
+        #              projection[:, self._mesh.num_ghost_cells:
+        #                         -self._mesh.num_ghost_cells])
+        #
+        # return right_hand_side