diff --git a/Snakefile b/Snakefile
index 9a4e7911895b3eadf9025bdc537510b1230771e0..8252eea6e3d394955c568a8f0ba5653d0630af49 100644
--- a/Snakefile
+++ b/Snakefile
@@ -9,8 +9,6 @@ TODO: Contemplate separating cell average and reconstruction calculations
TODO: Contemplate removing Methods section from class docstring
TODO: Contemplate containing the quadrature application for plots in Mesh
TODO: Contemplate containing coarse mesh generation in Mesh
-TODO: Contemplate extracting boundary condition from InitialCondition
-TODO: Contemplate containing boundary condition in Mesh
TODO: Ask whether all quadratures depend on freely chosen num_nodes
TODO: Contemplate saving training data for each IC separately
TODO: Contemplate removing TrainingDataGenerator class
@@ -24,16 +22,34 @@ TODO: Contemplate using lambdify for basis
TODO: Contemplate allowing vector input for ICs
TODO: Discuss how wavelet details should be plotted
+TODO: Ask after reasoning behind plotting approx for this specific Burgers case
+TODO: Ask what derivative basis represents
+TODO: Ask after Dirichlet boundary setup
+TODO: Ask after exact solution for different equations
+TODO: Ask where Euler 1D exact solution comes from
+TODO: Discuss first calculating plotting data and then plotting in separate
+ step
+TODO: Ask whether the option of multiple ghost cells has any potential use
+ cases
+TODO: Discuss necessary object dependencies
+TODO: Contemplate extracting boundary condition from InitialCondition
+TODO: Contemplate containing boundary condition in Mesh
+TODO: Discuss extracting all objects from UpdateScheme and giving equation as
+ input instead
+
Urgent:
-TODO: Check whether 'projection' is always a ndarray
-TODO: Check whether ghost cells are handled/set correctly
-TODO: Make sure that the cell indices are the same over all TCDs
-TODO: Make sure TCs are reported as ndarray
-TODO: Extract equation-related functionality into LinearAdvection class
-TODO: Use cfl_number for updating, not just time (equation-related?)
+TODO: Check whether 'projection' is always a ndarray -> Done
+TODO: Add property attribute for num_ghost_cells -> Done
+TODO: Ensure compatibility with arbitrary number of ghost cells -> Done
+TODO: Move 'do_initial_projection()' to 'projection_utils.py' -> Done
+TODO: Extract equation-related functionality into LinearAdvection class -> Done
+TODO: Enforce num_ghost_cells to be positive integer for DG (not training)
TODO: Add Burger class
+TODO: Use cfl_number for updating, not just time (equation-related?)
Critical, but not urgent:
+TODO: Make sure TCs are reported as ndarray
+TODO: Make sure that the cell indices are the same over all TCDs
TODO: Combine ANN workflows if feasible
TODO: Check whether all instance variables are sensible
TODO: Check whether all ValueError are set (correctly)
@@ -49,6 +65,7 @@ TODO: Enforce even number of ghost cells on each side on fine mesh (?)
TODO: Create g-mesh with Mesh class
TODO: Add images to report
TODO: Add verbose output
+TODO: Add tests for all functions
Not feasible yet or doc-related:
TODO: Move plot_approximation_results() into plotting script
diff --git a/scripts/approximate_solution.py b/scripts/approximate_solution.py
index bbda5a4d0bebe7ed6c75bcc7b99047349f4c1499..d1c6c28d48561a569c097a9eac01cd80101b0867 100644
--- a/scripts/approximate_solution.py
+++ b/scripts/approximate_solution.py
@@ -14,6 +14,7 @@ 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:
@@ -57,15 +58,6 @@ def main() -> None:
config=detector_config,
mesh=mesh, basis=basis)
- # Initialize update scheme after checking its legitimacy
- scheme_name = params.pop('update_scheme', 'SSPRK3')
- wave_speed = params.pop('wave_speed', 1)
- if not hasattr(Update_Scheme, scheme_name):
- raise ValueError('Invalid update scheme: "%s"' % scheme_name)
- update_scheme = getattr(Update_Scheme, scheme_name)(
- polynomial_degree=basis.polynomial_degree, mesh=mesh,
- detector=detector, limiter=limiter, wave_speed=wave_speed)
-
# Initialize quadrature after checking its legitimacy
quadrature_name = params.pop('quadrature', 'Gauss')
quadrature_config = params.pop('quadrature_config', {})
@@ -82,10 +74,25 @@ def main() -> None:
% init_name)
init_cond = getattr(Initial_Condition, init_name)(config=init_config)
- dg_scheme = DGScheme(detector=detector, quadrature=quadrature,
- init_cond=init_cond, update_scheme=update_scheme,
- mesh=mesh, basis=basis, wave_speed=wave_speed,
- **params)
+ # Initialize linear advection equation
+ wave_speed = params.pop('wave_speed', 1)
+ final_time = params.pop('final_time', 1)
+ cfl_number = params.pop('cfl_number', 0.2)
+ equation = LinearAdvection(final_time=final_time,
+ cfl_number=cfl_number,
+ basis=basis, init_cond=init_cond,
+ quadrature=quadrature,
+ mesh=mesh, wave_speed=wave_speed)
+
+ # Initialize update scheme after checking its legitimacy
+ scheme_name = params.pop('update_scheme', 'SSPRK3')
+ if not hasattr(Update_Scheme, scheme_name):
+ raise ValueError('Invalid update scheme: "%s"' % scheme_name)
+ update_scheme = getattr(Update_Scheme, scheme_name)(
+ detector=detector, limiter=limiter, equation=equation)
+
+ dg_scheme = DGScheme(detector=detector, update_scheme=update_scheme,
+ equation=equation, **params)
dg_scheme.approximate(
data_file=snakemake.params['plot_dir'] + '/' +
diff --git a/scripts/tcd/DG_Approximation.py b/scripts/tcd/DG_Approximation.py
index a4cf1df3f853e6af8aa3ca26f05bb7d8eb994ff2..5f38c2649e03f8a3cb3f920c848ccdc6eadee4ba 100644
--- a/scripts/tcd/DG_Approximation.py
+++ b/scripts/tcd/DG_Approximation.py
@@ -6,12 +6,8 @@
import json
import math
-from .Basis_Function import Basis
from .encoding_utils import encode_ndarray
-from .Initial_Condition import InitialCondition
-from .Mesh import Mesh
-from .projection_utils import do_initial_projection
-from .Quadrature import Quadrature
+from .Equation import Equation
from .Troubled_Cell_Detector import TroubledCellDetector
from .Update_Scheme import UpdateScheme
@@ -33,34 +29,20 @@ class DGScheme:
"""
def __init__(self, detector: TroubledCellDetector,
- quadrature: Quadrature, init_cond: InitialCondition,
- update_scheme: UpdateScheme, mesh: Mesh, basis: Basis,
- wave_speed, **kwargs):
+ update_scheme: UpdateScheme, equation: Equation, **kwargs):
"""Initializes DGScheme.
Parameters
----------
detector : TroubledCellDetector object
Troubled cell detector.
- quadrature : Quadrature object
- Quadrature for evaluation.
- init_cond : InitialCondition object
- Initial condition for evaluation.
update_scheme : UpdateScheme object
Update scheme for evaluation.
- mesh : Mesh object
- Mesh for calculation.
- basis : Basis object
- Basis for calculation.
- wave_speed : float, optional
- Speed of wave in rightward direction.
+ equation : Equation object
+ Equation to approximate.
Other Parameters
----------------
- cfl_number : float, optional
- CFL number to ensure stability. Default: 0.2.
- final_time : float, optional
- Final time for which approximation is calculated. Default: 1.
verbose : bool, optional
Flag whether commentary in console is wanted. Default: False.
history_threshold : float, optional
@@ -69,19 +51,13 @@ class DGScheme:
"""
self._detector = detector
- self._quadrature = quadrature
- self._init_cond = init_cond
self._update_scheme = update_scheme
- self._mesh = mesh
- self._basis = basis
- self._wave_speed = wave_speed
+ self._equation = equation
# Unpack keyword arguments
- self._cfl_number = kwargs.pop('cfl_number', 0.2)
- self._final_time = kwargs.pop('final_time', 1)
self._verbose = kwargs.pop('verbose', False)
- self._history_threshold = kwargs.pop('history_threshold',
- math.ceil(0.2/self._cfl_number))
+ self._history_threshold = kwargs.pop(
+ 'history_threshold', math.ceil(0.2/self._equation.cfl_number))
# Throw an error if there are extra keyword arguments
if len(kwargs) > 0:
@@ -103,23 +79,14 @@ class DGScheme:
Path to file in which data will be saved.
"""
- projection = do_initial_projection(
- init_cond=self._init_cond, mesh=self._mesh,
- basis=self._basis, quadrature=self._quadrature)
-
- time_step = abs(self._cfl_number * self._mesh.cell_len /
- self._wave_speed)
-
+ time_step, projection = self._equation.initialize_approximation()
current_time = 0
iteration = 0
troubled_cell_history = []
time_history = []
- while current_time < self._final_time:
- # Adjust for last cell
- cfl_number = self._cfl_number
- 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
+ while current_time < self._equation.final_time:
+ cfl_number, time_step = self._equation.update_time_step(
+ current_time=current_time, time_step=time_step)
# Update projection
projection, troubled_cells = self._update_scheme.step(projection,
@@ -136,8 +103,8 @@ class DGScheme:
approx_stats = self._detector.create_data_dict(projection)
# Save approximation results in dictionary
- approx_stats['wave_speed'] = self._wave_speed
- approx_stats['final_time'] = self._final_time
+ approx_stats['wave_speed'] = self._equation.wave_speed
+ approx_stats['final_time'] = self._equation.final_time
approx_stats['time_history'] = time_history
approx_stats['troubled_cell_history'] = troubled_cell_history
diff --git a/scripts/tcd/Equation.py b/scripts/tcd/Equation.py
new file mode 100644
index 0000000000000000000000000000000000000000..e0e7c21cfcc19befeb0ee50dff89d14b70908485
--- /dev/null
+++ b/scripts/tcd/Equation.py
@@ -0,0 +1,347 @@
+# -*- coding: utf-8 -*-
+"""Module for equation class.
+
+@author: Laura C. Kühle
+"""
+from typing import Tuple
+from abc import ABC, abstractmethod
+import numpy as np
+from numpy import ndarray
+
+from .Basis_Function import Basis
+from .Initial_Condition import InitialCondition
+from .Mesh import Mesh
+from .projection_utils import do_initial_projection
+from .Quadrature import Quadrature
+
+
+class Equation(ABC):
+ """Abstract class for equation.
+
+ Methods
+ -------
+ initialize_approximation()
+ Initialize projection and time step for approximation.
+ update_time_step(current_time, time_step)
+ Update time step.
+ solve_exactly(mesh)
+ Calculate exact solution.
+ enforce_boundary_condition(projection)
+ Enforce boundary condition.
+ update_right_hand_side(projection)
+ Update right-hand side.
+ """
+
+ def __init__(self, quadrature: Quadrature, init_cond: InitialCondition,
+ basis: Basis, mesh: Mesh, final_time: float,
+ wave_speed: float, cfl_number: float) -> None:
+ """Initialize equation class.
+
+ Parameters
+ ----------
+ quadrature : Quadrature object
+ Quadrature for evaluation.
+ init_cond : InitialCondition object
+ Initial condition for evaluation.
+ basis : Basis object
+ Basis for calculation.
+ mesh : Mesh object
+ Mesh for calculation.
+ final_time : float
+ Final time for which approximation is calculated.
+ wave_speed : float
+ Speed of wave in rightward direction.
+ cfl_number : float
+ CFL number to ensure stability.
+
+ """
+ self._quadrature = quadrature
+ self._init_cond = init_cond
+ self._basis = basis
+ self._mesh = mesh
+ self._final_time = final_time
+ self._wave_speed = wave_speed
+ self._cfl_number = cfl_number
+
+ self._reset()
+
+ @property
+ def final_time(self) -> float:
+ """Return final time."""
+ return self._final_time
+
+ @property
+ def wave_speed(self) -> float:
+ """Return wave speed."""
+ return self._wave_speed
+
+ @property
+ def cfl_number(self) -> float:
+ """Return CFL number."""
+ return self._cfl_number
+
+ def _reset(self) -> None:
+ """Reset instance variables."""
+ pass
+
+ def initialize_approximation(self) -> Tuple[float, ndarray]:
+ """Initialize projection and time step for approximation."""
+ return self._initialize_time_step(), self._initialize_projection()
+
+ def _initialize_time_step(self):
+ """Initialize time step."""
+ return abs(self._cfl_number * self._mesh.cell_len / self._wave_speed)
+
+ @abstractmethod
+ def _initialize_projection(self) -> ndarray:
+ """Initialize projection."""
+ pass
+
+ @abstractmethod
+ def update_time_step(self, current_time: float,
+ time_step: float) -> Tuple[float, float]:
+ """Update time step.
+
+ Parameters
+ ----------
+ 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.
+
+ """
+ pass
+
+ @abstractmethod
+ 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
+
+ @abstractmethod
+ def enforce_boundary_condition(self, projection: ndarray) -> ndarray:
+ """Enforce boundary condition.
+
+ Adjust ghost cells to ensure periodic boundary condition.
+
+ Parameters
+ ----------
+ projection : ndarray
+ Matrix of projection for each polynomial degree.
+
+ Returns
+ -------
+ projection : ndarray
+ Matrix of projection for each polynomial degree.
+
+ """
+ pass
+
+ @abstractmethod
+ 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
+
+
+class LinearAdvection(Equation):
+ """Class for linear advection equation.
+
+ Attributes
+ ----------
+ stiffness_matrix : ndarray
+ Stiffness matrix.
+ boundary_matrix : ndarray
+ Boundary matrix.
+
+ Methods
+ -------
+ update_right_hand_side(projection)
+ Update right-hand side.
+ enforce_boundary_condition(projection)
+ Enforce boundary condition.
+
+ Notes
+ -----
+ .. math:: u_t + u_x = 0
+
+ """
+
+ def _reset(self) -> None:
+ """Reset instance variables."""
+ 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 stiffness 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._stiffness_matrix = matrix * degree_matrix
+
+ # Set boundary 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._boundary_matrix = matrix * degree_matrix
+
+ 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, current_time: float,
+ time_step: float) -> Tuple[float, float]:
+ """Update time step.
+
+ Parameters
+ ----------
+ 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
+ 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]:
+ """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.
+
+ """
+ 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)
+ exact = np.array([self._init_cond.calculate(
+ mesh=mesh, x=point-self.wave_speed * self.final_time+num_periods *
+ mesh.interval_len)
+ for point in grid])
+
+ grid = np.reshape(grid, (1, grid.size))
+ exact = np.reshape(exact, (1, exact.size))
+
+ return grid, exact
+
+ def enforce_boundary_condition(self, projection: ndarray) -> ndarray:
+ """Enforce boundary condition.
+
+ Adjust ghost cells to ensure periodic boundary condition.
+
+ Parameters
+ ----------
+ projection : ndarray
+ Matrix of projection for each polynomial degree.
+
+ Returns
+ -------
+ projection : ndarray
+ Matrix of projection for each polynomial degree.
+
+ """
+ projection[:, :self._mesh.num_ghost_cells] = projection[
+ :, -2*self._mesh.num_ghost_cells:-self._mesh.num_ghost_cells]
+ projection[:, -self._mesh.num_ghost_cells:] = projection[
+ :, self._mesh.num_ghost_cells:2*self._mesh.num_ghost_cells]
+ return projection
+
+ 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.
+
+ """
+ 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._boundary_matrix @
+ projection[:, self._mesh.num_ghost_cells-1:
+ -self._mesh.num_ghost_cells-1] +
+ self._stiffness_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._boundary_matrix @
+ projection[:, self._mesh.num_ghost_cells+1:] +
+ self._stiffness_matrix @
+ projection[:, self._mesh.num_ghost_cells:
+ -self._mesh.num_ghost_cells])
+
+ # Set ghost cells to respective value
+ right_hand_side[:, :self._mesh.num_ghost_cells] = right_hand_side[
+ :, -2*self._mesh.num_ghost_cells:-self._mesh.num_ghost_cells]
+ right_hand_side[:, -self._mesh.num_ghost_cells:] = right_hand_side[
+ :, self._mesh.num_ghost_cells:2*self._mesh.num_ghost_cells]
+
+ return right_hand_side
diff --git a/scripts/tcd/Update_Scheme.py b/scripts/tcd/Update_Scheme.py
index 19cd20165f8e365ee44350e2b901a3d6aecab55b..e7c9122b13803a0d41fecd91bfac7177ebd21abe 100644
--- a/scripts/tcd/Update_Scheme.py
+++ b/scripts/tcd/Update_Scheme.py
@@ -4,20 +4,14 @@
"""
from abc import ABC, abstractmethod
-import numpy as np
import time
+from .Equation import Equation
+
class UpdateScheme(ABC):
"""Abstract class for updating projections at a time step.
- Attributes
- ----------
- stiffness_matrix : ndarray
- Matrix
- boundary_matrix : ndarray
- Matrix
-
Methods
-------
get_name()
@@ -26,56 +20,28 @@ class UpdateScheme(ABC):
Performs time step.
"""
- def __init__(self, polynomial_degree, mesh, detector, limiter,
- wave_speed):
+ def __init__(self, detector, limiter, equation):
"""Initializes UpdateScheme.
Parameters
----------
- polynomial_degree : int
- Polynomial degree.
- mesh : Mesh
- Mesh for calculation.
detector : TroubledCellDetector object
Troubled cell detector for evaluation.
limiter : Limiter object
Limiter for evaluation.
- wave_speed : float
- Speed of wave in rightward direction.
+ equation: Equation
+ Equation.
"""
# Unpack positional arguments
- self._polynomial_degree = polynomial_degree
- self._mesh = mesh
self._detector = detector
self._limiter = limiter
- self._wave_speed = wave_speed
+ self._equation = equation
self._reset()
def _reset(self):
"""Resets instance variables."""
- matrix_shape = (self._polynomial_degree+1, self._polynomial_degree+1)
- root_vector = np.array([np.sqrt(degree+0.5) for degree in range(
- self._polynomial_degree+1)])
- degree_matrix = np.matmul(root_vector[:, np.newaxis],
- root_vector[:, np.newaxis].T)
-
- # Set stiffness 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._stiffness_matrix = matrix * degree_matrix
-
- # Set boundary 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._boundary_matrix = matrix * degree_matrix
def get_name(self):
"""Returns string of class name."""
@@ -150,32 +116,6 @@ class UpdateScheme(ABC):
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[:, :self._mesh.num_ghost_cells] = \
- current_projection[
- :, -2*self._mesh.num_ghost_cells:-self._mesh.num_ghost_cells]
- current_projection[:, -self._mesh.num_ghost_cells:] = \
- current_projection[
- :, self._mesh.num_ghost_cells:2*self._mesh.num_ghost_cells]
- return current_projection
-
class SSPRK3(UpdateScheme):
"""Class for strong stability-preserving Runge Kutta of order 3.
@@ -210,14 +150,14 @@ class SSPRK3(UpdateScheme):
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 = self._equation.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 = self._equation.enforce_boundary_condition(
current_projection)
current_projection = self._apply_third_step(original_projection,
@@ -225,7 +165,7 @@ class SSPRK3(UpdateScheme):
cfl_number)
current_projection, troubled_cells = self._apply_limiter(
current_projection)
- current_projection = self._enforce_boundary_condition(
+ current_projection = self._equation.enforce_boundary_condition(
current_projection)
return current_projection, troubled_cells
@@ -246,7 +186,8 @@ class SSPRK3(UpdateScheme):
Matrix of updated projection for each polynomial degree.
"""
- right_hand_side = self._update_right_hand_side(original_projection)
+ right_hand_side = self._equation.update_right_hand_side(
+ original_projection)
return original_projection + (cfl_number*right_hand_side)
def _apply_second_step(self, original_projection, current_projection,
@@ -269,7 +210,8 @@ class SSPRK3(UpdateScheme):
Matrix of updated projection for each polynomial degree.
"""
- right_hand_side = self._update_right_hand_side(current_projection)
+ right_hand_side = self._equation.update_right_hand_side(
+ current_projection)
return 1/4 * (3*original_projection
+ (current_projection + cfl_number*right_hand_side))
@@ -293,48 +235,7 @@ class SSPRK3(UpdateScheme):
Matrix of updated projection for each polynomial degree.
"""
- right_hand_side = self._update_right_hand_side(current_projection)
+ right_hand_side = self._equation.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.
-
- """
- right_hand_side = np.zeros_like(current_projection)
- if self._wave_speed > 0:
- right_hand_side[:, self._mesh.num_ghost_cells:
- -self._mesh.num_ghost_cells] = \
- 2 * (self._boundary_matrix @
- current_projection[:, self._mesh.num_ghost_cells-1:
- -self._mesh.num_ghost_cells-1] +
- self._stiffness_matrix @
- current_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._boundary_matrix @
- current_projection[:, self._mesh.num_ghost_cells+1:] +
- self._stiffness_matrix @
- current_projection[:, self._mesh.num_ghost_cells:
- -self._mesh.num_ghost_cells])
-
- # Set ghost cells to respective value
- right_hand_side[:, :self._mesh.num_ghost_cells] = right_hand_side[
- :, -2*self._mesh.num_ghost_cells:-self._mesh.num_ghost_cells]
- right_hand_side[:, -self._mesh.num_ghost_cells:] = right_hand_side[
- :, self._mesh.num_ghost_cells:2*self._mesh.num_ghost_cells]
-
- return right_hand_side