diff --git a/Snakefile b/Snakefile
index 2d1b95adfefafbc5cf442a7f75c2ef2145de5a04..d66ead605992637c30dad60575c7fac8e54508ac 100644
--- a/Snakefile
+++ b/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:
diff --git a/scripts/tcd/Equation.py b/scripts/tcd/Equation.py
index e2765d25e7f1831dd68588c53f4fa2c3978adb80..d0bd5ceb8ad76ebc9dc6b07a44ceb21398aca9ab 100644
--- a/scripts/tcd/Equation.py
+++ b/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