# -*- coding: utf-8 -*-
"""
@author: Laura C. Kühle
TODO: Discuss descriptions (matrices, cfl number, right-hand side,
limiting slope, basis, wavelet, etc.)
TODO: Discuss referencing info on SSPRK3
"""
from abc import ABC, abstractmethod
import numpy as np
import time
class UpdateScheme(ABC):
"""Abstract class for updating projections at a time step.
Attributes
----------
stiffness_matrix : ndarray
Matrix
boundary_matrix : ndarray
Matrix
Methods
-------
get_name()
Returns string of class name.
step(projection, cfl_number)
Performs time step.
"""
def __init__(self, polynomial_degree, num_grid_cells, detector, limiter):
"""Initializes UpdateScheme.
Parameters
----------
polynomial_degree : int
Polynomial degree.
num_grid_cells : int
Number of cells in the mesh. Usually exponential of 2.
detector : TroubledCellDetector object
Troubled cell detector for evaluation.
limiter : Limiter object
Limiter for evaluation.
"""
# Unpack positional arguments
self._polynomial_degree = polynomial_degree
self._num_grid_cells = num_grid_cells
self._detector = detector
self._limiter = limiter
self._reset()
def _reset(self):
"""Resets instance variables."""
# Set stiffness matrix
matrix = []
for i in range(self._polynomial_degree+1):
new_row = []
for j in range(self._polynomial_degree+1):
new_entry = -1.0
if (j < i) & ((i+j) % 2 == 1):
new_entry = 1.0
new_row.append(new_entry*np.sqrt((i+0.5) * (j+0.5)))
matrix.append(new_row)
self._stiffness_matrix = np.array(matrix)
# former: inv_mass @ np.array(matrix)
# Set boundary matrix
matrix = []
for i in range(self._polynomial_degree+1):
new_row = []
for j in range(self._polynomial_degree+1):
new_entry = np.sqrt((i+0.5) * (j+0.5)) * (-1.0)**i
new_row.append(new_entry)
matrix.append(new_row)
self._boundary_matrix = np.array(matrix)
# former: inv_mass @ np.array(matrix)
def get_name(self):
"""Returns string of class name."""
return self.__class__.__name__
def step(self, projection, cfl_number):
"""Performs time step.
Parameters
----------
projection : ndarray
Matrix of projection for each polynomial degree.
cfl_number : float
CFL number to ensure stability.
Returns
-------
current_projection : ndarray
Matrix of projection of current update step for each polynomial
degree.
troubled_cells : list
List of indices for all detected troubled cells.
"""
current_projection, troubled_cells = self._apply_stability_method(
projection, cfl_number)
return current_projection, troubled_cells
@abstractmethod
def _apply_stability_method(self, projection, cfl_number):
"""Applies stability method.
Parameters
----------
projection : ndarray
Matrix of projection for each polynomial degree.
cfl_number : float
CFL number to ensure stability.
Returns
-------
current_projection : ndarray
Matrix of projection of current update step for each polynomial
degree.
troubled_cells : list
List of indices for all detected troubled cells.
"""
pass
def _apply_limiter(self, current_projection):
"""Applies limiter on troubled cells.
Parameters
----------
current_projection : ndarray
Matrix of projection of current update step for each polynomial
degree.
Returns
-------
new_projection : ndarray
Matrix of updated projection for each polynomial degree.
troubled_cells : list
List of indices for all detected troubled cells.
"""
troubled_cells = self._detector.get_cells(current_projection)
new_projection = current_projection.copy()
for cell in troubled_cells:
new_projection[:, cell] = self._limiter.apply(current_projection,
cell)
return new_projection, troubled_cells
def _enforce_boundary_condition(self, current_projection):
"""Enforces boundary condition.
Adjusts ghost cells to ensure periodic boundary condition.
Parameters
----------
current_projection : ndarray
Matrix of projection of current update step for each polynomial
degree.
Returns
-------
current_projection : ndarray
Matrix of projection of current update step for each polynomial
degree.
"""
current_projection[:, 0] = current_projection[:, self._num_grid_cells]
current_projection[:, self._num_grid_cells+1] \
= current_projection[:, 1]
return current_projection
class SSPRK3(UpdateScheme):
"""Class for strong stability-preserving Runge Kutta of order 3.
Notes
-----
Reference (?)
"""
# Override method of superclass
def _apply_stability_method(self, projection, cfl_number):
"""Applies stability method.
Parameters
----------
projection : ndarray
Matrix of projection for each polynomial degree.
cfl_number : float
CFL number to ensure stability.
Returns
-------
current_projection : ndarray
Matrix of projection of current update step for each polynomial
degree.
troubled_cells : list
List of indices for all detected troubled cells.
"""
original_projection = projection
current_projection = self._apply_first_step(original_projection,
cfl_number)
current_projection, __ = self._apply_limiter(current_projection)
current_projection = self._enforce_boundary_condition(
current_projection)
current_projection = self._apply_second_step(original_projection,
current_projection,
cfl_number)
current_projection, __ = self._apply_limiter(current_projection)
current_projection = self._enforce_boundary_condition(
current_projection)
current_projection = self._apply_third_step(original_projection,
current_projection,
cfl_number)
current_projection, troubled_cells = self._apply_limiter(
current_projection)
current_projection = self._enforce_boundary_condition(
current_projection)
return current_projection, troubled_cells
def _apply_first_step(self, original_projection, cfl_number):
"""Applies first step of SSPRK3.
Parameters
----------
original_projection : ndarray
Matrix of original projection for each polynomial degree.
cfl_number : float
CFL number to ensure stability.
Returns
-------
ndarray
Matrix of updated projection for each polynomial degree.
"""
right_hand_side = self._update_right_hand_side(original_projection)
return original_projection + (cfl_number*right_hand_side)
def _apply_second_step(self, original_projection, current_projection,
cfl_number):
"""Applies second step of SSPRK3.
Parameters
----------
original_projection : ndarray
Matrix of original projection for each polynomial degree.
current_projection : ndarray
Matrix of projection of current update step for each polynomial
degree.
cfl_number : float
CFL number to ensure stability.
Returns
-------
ndarray
Matrix of updated projection for each polynomial degree.
"""
right_hand_side = self._update_right_hand_side(current_projection)
return 1/4 * (3*original_projection
+ (current_projection + cfl_number*right_hand_side))
def _apply_third_step(self, original_projection, current_projection,
cfl_number):
"""Applies third step of SSPRK3.
Parameter
---------
original_projection : ndarray
Matrix of original projection for each polynomial degree.
current_projection : ndarray
Matrix of projection of current update step for each polynomial
degree.
cfl_number : float
CFL number to ensure stability.
Returns
-------
ndarray
Matrix of updated projection for each polynomial degree.
"""
right_hand_side = self._update_right_hand_side(current_projection)
return 1/3 * (original_projection
+ 2*(current_projection + cfl_number*right_hand_side))
def _update_right_hand_side(self, current_projection):
"""Updates right-hand side.
Parameter
---------
current_projection : ndarray
Matrix of projection of current update step for each polynomial
degree.
Returns
-------
ndarray
Matrix of right-hand side.
"""
# Initialize vector and set first entry to accommodate for ghost cell
right_hand_side = [0]
for j in range(self._num_grid_cells):
right_hand_side.append(2*(self._stiffness_matrix
@ current_projection[:, j+1]
+ self._boundary_matrix
@ current_projection[:, j]))
# Set ghost cells to respective value
right_hand_side[0] = right_hand_side[self._num_grid_cells]
right_hand_side.append(right_hand_side[1])
return np.transpose(right_hand_side)