Select Git revision
gradle-wrapper.properties
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
Update_Scheme.py 10.78 KiB
# -*- coding: utf-8 -*-
"""
@author: Laura C. Kühle
"""
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_cells, detector, limiter,
wave_speed):
"""Initializes UpdateScheme.
Parameters
----------
polynomial_degree : int
Polynomial degree.
num_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.
wave_speed : float
Speed of wave in rightward direction.
"""
# Unpack positional arguments
self._polynomial_degree = polynomial_degree
self._num_cells = num_cells
self._detector = detector
self._limiter = limiter
self._wave_speed = wave_speed
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):
if self._wave_speed > 0:
new_entry = -1.0
if (j < i) & ((i+j) % 2 == 1):
new_entry = 1.0
else:
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 \
if self._wave_speed > 0 \
else np.sqrt((i+0.5) * (j+0.5)) * (-1.0)**j
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 = self._limiter.apply(current_projection,
troubled_cells)
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_cells]
current_projection[:, self._num_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_cells):
if self._wave_speed > 0:
right_hand_side.append(2*(self._stiffness_matrix
@ current_projection[:, j+1]
+ self._boundary_matrix
@ current_projection[:, j]))
else:
right_hand_side.append(2*(self._stiffness_matrix
@ current_projection[:, j+1]
+ self._boundary_matrix
@ current_projection[:, j+2]))
# Set ghost cells to respective value
right_hand_side[0] = right_hand_side[self._num_cells]
right_hand_side.append(right_hand_side[1])
return np.transpose(right_hand_side)