Select Git revision
BAstCreator.java
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
Basis_Function.py 17.51 KiB
# -*- coding: utf-8 -*-
"""Module for polynomial basis.
@author: Laura C. Kühle
"""
from functools import cache
from typing import Tuple
import numpy as np
from numpy import ndarray
from sympy import Symbol, integrate
from projection_utils import calculate_approximate_solution
x = Symbol('x')
z = Symbol('z')
class Basis:
"""Class for polynomial basis.
Attributes
----------
polynomial degree : int
Polynomial degree.
basis : ndarray
Array of basis.
wavelet : ndarray
Array of wavelet.
inv_mass : ndarray
Inverse mass matrix.
basis_projection : Tuple[ndarray, ndarray]
Two arrays containing integrals of all basis vector
combinations evaluated on the left and right cell boundary,
respectively.
wavelet_projection : Tuple[ndarray, ndarray]
Two arrays containing integrals of all basis vector/
wavelet vector combinations evaluated on the left and right cell
boundary, respectively.
Methods
-------
calculate_cell_average(projection, stencil_length, add_reconstructions)
Calculate cell averages for a given projection.
"""
def __init__(self, polynomial_degree: int) -> None:
"""Initialize Basis.
Parameters
----------
polynomial_degree : int
Polynomial degree.
"""
self._polynomial_degree = polynomial_degree
@property
def polynomial_degree(self) -> int:
"""Return polynomial degree."""
return self._polynomial_degree
@property
@cache
def basis(self) -> ndarray:
"""Return basis vector."""
return self._build_basis_vector(x)
def _build_basis_vector(self, eval_point: float) -> ndarray:
"""Construct basis vector.
Parameters
----------
eval_point : float
Evaluation point.
Returns
-------
ndarray
Vector containing basis evaluated at evaluation point.
"""
pass
@property
@cache
def wavelet(self) -> ndarray:
"""Return wavelet vector."""
return self._build_wavelet_vector(z)
def _build_wavelet_vector(self, eval_point: float) -> ndarray:
"""Construct wavelet vector.
Parameters
----------
eval_point : float
Evaluation point.
Returns
-------
ndarray
Vector containing wavelet evaluated at evaluation point.
"""
pass
@property
@cache
def inverse_mass_matrix(self) -> ndarray:
"""Return inverse mass matrix."""
return self._build_inverse_mass_matrix()
def _build_inverse_mass_matrix(self) -> ndarray:
"""Construct inverse mass matrix.
Returns
-------
ndarray
Inverse mass matrix.
"""
pass
@property
@cache
def basis_projection(self) -> Tuple[ndarray, ndarray]:
"""Return basis projection."""
return np.array([]), np.array([])
@property
@cache
def wavelet_projection(self) -> Tuple[ndarray, ndarray]:
"""Return wavelet projection."""
return np.array([]), np.array([])
def calculate_cell_average(self, projection: ndarray, stencil_length: int,
add_reconstructions: bool = True) -> ndarray:
"""Calculate cell averages for a given projection.
Calculate the cell averages of all cells in a projection.
If desired, reconstructions are calculated for the middle cell
and added left and right to it, respectively.
Parameters
----------
projection : ndarray
Matrix of projection for each polynomial degree.
stencil_length : int
Size of data array.
add_reconstructions : bool, optional
Flag whether reconstructions of the middle cell are included.
Default: True.
Returns
-------
ndarray
Matrix containing cell averages (and reconstructions) for given
projection.
"""
cell_averages = calculate_approximate_solution(
projection, np.array([0]), 0, self.basis)
if add_reconstructions:
middle_idx = stencil_length // 2
left_reconstructions, right_reconstructions = \
self._calculate_reconstructions(
projection[:, middle_idx:middle_idx+1])
return np.array(list(map(
np.float64, zip(cell_averages[:, :middle_idx],
left_reconstructions,
cell_averages[:, middle_idx],
right_reconstructions,
cell_averages[:, middle_idx+1:]))))
return np.array(list(map(np.float64, cell_averages)))
def _calculate_reconstructions(self, projection: ndarray) \
-> Tuple[ndarray, ndarray]:
"""Calculate left and right reconstructions for a given projection.
Parameters
----------
projection : ndarray
Matrix of projection for each polynomial degree.
Returns
-------
left_reconstruction : ndarray
List containing left reconstructions for given projection.
right_reconstruction : ndarray
List containing right reconstructions for given projection.
"""
left_reconstructions = calculate_approximate_solution(
projection, np.array([-1]), self._polynomial_degree,
self.basis)
right_reconstructions = calculate_approximate_solution(
projection, np.array([1]), self._polynomial_degree,
self.basis)
return left_reconstructions, right_reconstructions
class Legendre(Basis):
"""Class for Legendre basis."""
def _build_basis_vector(self, eval_point: float) -> ndarray:
"""Construct basis vector.
Parameters
----------
eval_point : float
Evaluation point.
Returns
-------
ndarray
Vector containing basis evaluated at evaluation point.
"""
return self._calculate_legendre_vector(eval_point)
def _calculate_legendre_vector(self, eval_point: float) -> ndarray:
"""Construct Legendre vector.
Parameters
----------
eval_point : float
Evaluation point.
Returns
-------
ndarray
Vector containing Legendre polynomial evaluated at evaluation
point.
"""
vector = []
for degree in range(self._polynomial_degree+1):
if degree == 0:
vector.append(1.0 + 0*eval_point)
else:
if degree == 1:
vector.append(eval_point)
else:
poly = (2.0*degree - 1)/degree * eval_point * vector[-1] \
- (degree-1)/degree * vector[-2]
vector.append(poly)
return np.array(vector)
class OrthonormalLegendre(Legendre):
"""Class for orthonormal Legendre basis."""
def _build_basis_vector(self, eval_point: float) -> ndarray:
"""Construct basis vector.
Parameters
----------
eval_point : float
Evaluation point.
Returns
-------
ndarray
Vector containing basis evaluated at evaluation point.
"""
leg_vector = self._calculate_legendre_vector(eval_point)
return np.array([leg_vector[degree] * np.sqrt(degree+0.5)
for degree in range(self._polynomial_degree+1)])
def _build_wavelet_vector(self, eval_point: float) -> ndarray:
"""Construct wavelet vector.
Parameters
----------
eval_point : float
Evaluation point.
Returns
-------
ndarray
Vector containing wavelet evaluated at evaluation point.
Notes
-----
Hardcoded version only for now.
"""
degree = self._polynomial_degree
if degree == 0:
return np.array([np.sqrt(0.5) + eval_point*0])
if degree == 1:
return np.array([np.sqrt(1.5) * (-1 + 2*eval_point),
np.sqrt(0.5) * (-2 + 3*eval_point)])
if degree == 2:
return np.array([1/3 * np.sqrt(0.5) *
(1 - 24*eval_point + 30*(eval_point**2)),
1/2 * np.sqrt(1.5) *
(3 - 16*eval_point + 15*(eval_point**2)),
1/3 * np.sqrt(2.5) *
(4 - 15*eval_point + 12*(eval_point**2))])
if degree == 3:
return np.array([np.sqrt(15/34) *
(1 + 4*eval_point - 30*(eval_point**2)
+ 28*(eval_point**3)),
np.sqrt(1/42) *
(-4 + 105*eval_point - 300*(eval_point**2)
+ 210*(eval_point**3)),
1/2 * np.sqrt(35/34) *
(-5 + 48*eval_point - 105*(eval_point**2)
+ 64*(eval_point**3)),
1/2 * np.sqrt(5/42) *
(-16 + 105*eval_point - 192*(eval_point**2)
+ 105*(eval_point**3))])
if degree == 4:
return np.array([np.sqrt(1/186) *
(1 + 30*eval_point + 210*(eval_point**2)
- 840*(eval_point**3) + 630*(eval_point**4)),
0.5 * np.sqrt(1/38) *
(-5 - 144*eval_point + 1155*(eval_point**2)
- 2240*(eval_point**3) + 1260*(eval_point**4)),
np.sqrt(35/14694) *
(22 - 735*eval_point + 3504*(eval_point**2)
- 5460*(eval_point**3) + 2700*(eval_point**4)),
1/8 * np.sqrt(21/38) *
(35 - 512*eval_point + 1890*(eval_point**2)
- 2560*(eval_point**3) + 1155*(eval_point**4)),
0.5 * np.sqrt(7/158) *
(32 - 315*eval_point + 960*(eval_point**2)
- 1155*(eval_point**3) + 480*(eval_point**4))])
raise ValueError('Invalid value: Alpert\'s wavelet is only available \
up to degree 4 for this application')
def _build_inverse_mass_matrix(self) -> ndarray:
"""Construct inverse mass matrix.
Returns
-------
ndarray
Inverse mass matrix.
"""
mass_matrix = []
for i in range(self._polynomial_degree+1):
new_row = []
for j in range(self._polynomial_degree+1):
new_entry = 0.0
if i == j:
new_entry = 1.0
new_row.append(new_entry)
mass_matrix.append(new_row)
return np.array(mass_matrix)
@property
@cache
def basis_projections(self) -> Tuple[ndarray, ndarray]:
"""Return basis projection.
Construct matrices containing the integrals of the
product of two basis vectors for every degree combination evaluated
at the left and right cell boundary.
Returns
-------
left_basis_projection : ndarray
Array containing the left basis projection.
right_basis_projection : ndarray
Array containing the right basis projection.
"""
left_basis_projection = self._build_basis_matrix(z, 0.5 * (z - 1))
right_basis_projection = self._build_basis_matrix(z, 0.5 * (z + 1))
return left_basis_projection, right_basis_projection
def _build_basis_matrix(self, first_param: float, second_param: float) \
-> ndarray:
"""Construct a basis matrix.
Parameters
----------
first_param : float
First parameter.
second_param : float
Second parameter.
Returns
-------
ndarray
Matrix containing the integral of basis products.
"""
matrix = []
for i in range(self._polynomial_degree + 1):
row = []
for j in range(self._polynomial_degree + 1):
entry = integrate(self.basis[i].subs(x, first_param)
* self.basis[j].subs(x, second_param),
(z, -1, 1))
row.append(np.float64(entry))
matrix.append(row)
return np.array(matrix)
@property
@cache
def multiwavelet_projections(self) -> Tuple[ndarray, ndarray]:
"""Return wavelet projection.
Construct matrices containing the integrals of the
product of a basis vector and a wavelet vector for every degree
combination evaluated at the left and right cell boundary.
Returns
-------
left_wavelet_projection : ndarray
Array containing the left multiwavelet projection.
right_wavelet_projection : ndarray
Array containing the right multiwavelet projection.
"""
left_wavelet_projection = self._build_multiwavelet_matrix(
z, -0.5*(z-1), True)
right_wavelet_projection = self._build_multiwavelet_matrix(
z, 0.5*(z+1), False)
return left_wavelet_projection, right_wavelet_projection
def _build_multiwavelet_matrix(self, first_param: float,
second_param: float, is_left_matrix: bool) \
-> ndarray:
"""Construct a multiwavelet matrix.
Parameters
----------
first_param : float
First parameter.
second_param : float
Second parameter.
is_left_matrix : bool
Flag whether the left matrix is calculated.
Returns
-------
ndarray
Matrix containing the integral of products of a basis and a wavelet
vector.
"""
matrix = []
for i in range(self._polynomial_degree+1):
row = []
for j in range(self._polynomial_degree+1):
entry = integrate(self.basis[i].subs(x, first_param)
* self.wavelet[j].subs(z, second_param),
(z, -1, 1))
if is_left_matrix:
entry = entry * (-1)**(j + self._polynomial_degree + 1)
row.append(np.float64(entry))
matrix.append(row)
return np.array(matrix)
def calculate_cell_average(self, projection: ndarray, stencil_length: int,
add_reconstructions: bool = True) -> ndarray:
"""Calculate cell averages for a given projection.
Calculate the cell averages of all cells in a projection.
If desired, reconstructions are calculated for the middle cell
and added left and right to it, respectively.
Notes
-----
To increase speed, this function uses a simplified calculation
specific to the orthonormal Legendre polynomial basis.
Parameters
----------
projection : ndarray
Matrix of projection for each polynomial degree.
stencil_length : int
Size of data array.
add_reconstructions : bool, optional
Flag whether reconstructions of the middle cell are included.
Default: True.
Returns
-------
ndarray
Matrix containing cell averages (and reconstructions) for given
projection.
"""
cell_averages = np.array([projection[0] / np.sqrt(2)])
if add_reconstructions:
middle_idx = stencil_length // 2
left_reconstructions, right_reconstructions = \
self._calculate_reconstructions(
projection[:, middle_idx:middle_idx+1])
return np.array(list(map(
np.float64, zip(cell_averages[:, :middle_idx],
left_reconstructions,
cell_averages[:, middle_idx],
right_reconstructions,
cell_averages[:, middle_idx+1:]))))
return np.array(list(map(np.float64, cell_averages)))
def _calculate_reconstructions(self, projection: ndarray) \
-> Tuple[list, list]:
"""Calculate left and right reconstructions for a given projection.
Parameters
----------
projection : ndarray
Matrix of projection for each polynomial degree.
Returns
-------
left_reconstruction : list
List containing left reconstructions for given projection.
right_reconstruction : list
List containing right reconstructions for given projection.
Notes
-----
To increase speed. this function uses a simplified calculation
specific to the orthonormal Legendre polynomial basis.
"""
left_reconstructions = [
sum(projection[degree][cell] * (-1)**degree
* np.sqrt(degree + 0.5)
for degree in range(self._polynomial_degree+1))
for cell in range(len(projection[0]))]
right_reconstructions = [
sum(projection[degree][cell] * np.sqrt(degree + 0.5)
for degree in range(self._polynomial_degree+1))
for cell in range(len(projection[0]))]
return left_reconstructions, right_reconstructions