Skip to content
Snippets Groups Projects
Commit dc0d74f0 authored by Laura Christine Kühle's avatar Laura Christine Kühle
Browse files

Moved basis/wavelet matrices to Vectors_of_Polynomials.

parent 849dae30
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
DG_Approximation.py
+ 6
− 3
View file @ dc0d74f0
......@@ -10,13 +10,15 @@ TODO: Replace loops with list comprehension if feasible
TODO: Find better names for A, B, anything pertaining Gauss
TODO: Write documentation for all methods
TODO: Contemplate moving basis/wavelet matrices to Vectors_of_Polynomials
TODO: Contemplate moving basis/wavelet matrices to Vectors_of_Polynomials -> Done
TODO: Restructure Vectors_of_Polynomials
TODO: Contemplate how to make shock tubes comparable
TODO: Improve saving of plots -> Done (moved to Troubled_Cell_Detector)
TODO: Added option to set plot directory -> Done
TODO: Extend color options -> Done
TODO: Implement type check for all kwargs and configs
TODO: Add option to not save plots -> Done (not call function in Main)
TODO: Set methods to static if possible -> Done
"""
import numpy as np
......@@ -140,7 +142,7 @@ class DGScheme(object):
# Set additional necessary instance variables
self._interval_len = self._right_bound-self._left_bound
self._cell_len = self._interval_len / self._num_grid_cells
self._basis = OrthonormalLegendre(self._polynom_degree).get_vector(x)
self._basis = OrthonormalLegendre(self._polynom_degree)
# Set additional necessary config parameters
self._limiter_config['cell_len'] = self._cell_len
......@@ -164,6 +166,7 @@ class DGScheme(object):
def _do_initial_projection(self):
# Initialize matrix and set first entry to accommodate for ghost cell
output_matrix = [0]
basis_vector = self._basis.get_vector(x)
for cell in range(self._num_grid_cells):
new_row = []
......@@ -172,7 +175,7 @@ class DGScheme(object):
for degree in range(self._polynom_degree + 1):
new_entry = sum(self._init_cond.calculate(
eval_point + self._cell_len/2 * self._quadrature.get_eval_points()[point])
* self._basis[degree].subs(x, self._quadrature.get_eval_points()[point])
* basis_vector[degree].subs(x, self._quadrature.get_eval_points()[point])
* self._quadrature.get_weights()[point]
for point in range(self._quadrature.get_num_points()))
new_row.append(np.float64(new_entry))
......
Troubled_Cell_Detector.py
+ 8
− 37
View file @ dc0d74f0
......@@ -137,7 +137,7 @@ class TroubledCellDetector(object):
def _calculate_approximate_solution(self, projection):
points = self._quadrature.get_eval_points()
num_points = self._quadrature.get_num_points()
basis = self._basis
basis = self._basis.get_vector(x)
basis_matrix = [[basis[degree].subs(x, points[point]) for point in range(num_points)]
for degree in range(self._polynom_degree+1)]
......@@ -193,27 +193,9 @@ class WaveletDetector(TroubledCellDetector):
self._colors['coarse_approx'] = self._colors.get('coarse_approx', 'y')
def _reset(self, config):
# Set fixed basis and wavelet vectors
self._basis = OrthonormalLegendre(self._polynom_degree).get_vector(x)
self._wavelet = AlpertsWavelet(self._polynom_degree).get_vector(x)
# Set additional necessary parameter
self._num_coarse_grid_cells = self._num_grid_cells//2
self._M1 = self._set_multiwavelet_matrix(xi, -0.5*(xi-1), True)
self._M2 = self._set_multiwavelet_matrix(xi, 0.5*(xi+1), False)
def _set_multiwavelet_matrix(self, first_param, second_param, is_M1):
matrix = []
for i in range(self._polynom_degree+1):
row = []
for j in range(self._polynom_degree+1):
entry = integrate(self._basis[i].subs(x, first_param) * self._wavelet[j].subs(x, second_param),
(xi, -1, 1))
if is_M1:
entry = entry * (-1)**(j + self._polynom_degree + 1)
row.append(np.float64(entry))
matrix.append(row)
return matrix
self._M1, self._M2 = self._basis.get_multiwavelet_matrices()
def get_cells(self, projection):
multiwavelet_coeffs = self._calculate_wavelet_coeffs(projection[:, 1: -1])
......@@ -240,6 +222,7 @@ class WaveletDetector(TroubledCellDetector):
coarse_projection = self._calculate_coarse_projection(projection)
multiwavelet_coeffs = self._calculate_wavelet_coeffs(projection)
wavelet = AlpertsWavelet(self._polynom_degree).get_vector(xi)
basis = self._basis.get_vector(x)
########################################################################################################################
# For later consideration
......@@ -248,8 +231,8 @@ class WaveletDetector(TroubledCellDetector):
# averaged_projection1 = []
# wavelet_projection1 = []
# for degree in range(self._polynom_degree + 1):
# leftMesh = coarse_projection[degree] * self._basis[degree].subs(x, -1 / 2)
# rightMesh = coarse_projection[degree] * self._basis[degree].subs(x, 1 / 2)
# leftMesh = coarse_projection[degree] * basis[degree].subs(x, -1 / 2)
# rightMesh = coarse_projection[degree] * basis[degree].subs(x, 1 / 2)
# leftTest = multiwavelet_coeffs[degree] * wavelet[degree].subs(xi, 1 / 2) \
# * (-1)**(self._polynom_degree + 1 + degree)
# rightTest = multiwavelet_coeffs[degree] * wavelet[degree].subs(xi, 1 / 2)
......@@ -267,7 +250,7 @@ class WaveletDetector(TroubledCellDetector):
########################################################################################################################
# tic = timeit.default_timer()
averaged_projection = [[coarse_projection[degree][cell] * self._basis[degree].subs(x, value)
averaged_projection = [[coarse_projection[degree][cell] * basis[degree].subs(x, value)
for cell in range(self._num_coarse_grid_cells)
for value in [-0.5, 0.5]]
for degree in range(self._polynom_degree + 1)]
......@@ -283,7 +266,7 @@ class WaveletDetector(TroubledCellDetector):
# print(wavelet_projection1 == wavelet_projection)
projected_coarse = np.sum(averaged_projection, axis=0)
projected_fine = np.sum([fine_projection[degree] * self._basis[degree].subs(x, 0)
projected_fine = np.sum([fine_projection[degree] * basis[degree].subs(x, 0)
for degree in range(self._polynom_degree + 1)], axis=0)
projected_wavelet_coeffs = np.sum(wavelet_projection, axis=0)
......@@ -296,8 +279,7 @@ class WaveletDetector(TroubledCellDetector):
plt.title('Wavelet Coefficients')
def _calculate_coarse_projection(self, projection):
basis_matrix_left = self._set_basis_matrix(xi, 0.5 * (xi - 1))
basis_matrix_right = self._set_basis_matrix(xi, 0.5 * (xi + 1))
basis_matrix_left, basis_matrix_right = self._basis.get_basis_matrices()
# Remove ghost cells
projection = projection[:, 1:-1]
......@@ -312,17 +294,6 @@ class WaveletDetector(TroubledCellDetector):
return coarse_projection
def _set_basis_matrix(self, first_param, second_param):
matrix = []
for i in range(self._polynom_degree + 1):
row = []
for j in range(self._polynom_degree + 1):
entry = integrate(self._basis[i].subs(x, first_param) * self._basis[j].subs(x, second_param),
(xi, -1, 1))
row.append(np.float64(entry))
matrix.append(row)
return matrix
def _plot_mesh(self, projection):
grid, exact = self._calculate_exact_solution(self._mesh[2:-2], self._cell_len)
approx = self._calculate_approximate_solution(projection[:, 1:-1])
......
Vectors_of_Polynomials.py
+ 61
− 1
View file @ dc0d74f0
......@@ -4,18 +4,39 @@
"""
import numpy as np
from sympy import Symbol, integrate
x = Symbol('x')
xi = Symbol('z')
class Vector(object):
def __init__(self, polynom_degree):
self._polynom_degree = polynom_degree
self._basis = self._set_basis_vector(x)
self._wavelet = self._set_wavelet_vector(x)
def get_vector(self, eval_points):
pass
def _set_basis_vector(self, eval_point):
return []
def _set_wavelet_vector(self, eval_point):
return []
def get_basis_matrices(self):
pass
def get_multiwavelet_matrices(self):
pass
class Legendre(Vector):
def get_vector(self, eval_point):
return self._basis
def _set_basis_vector(self, eval_point):
return self._calculate_legendre_vector(eval_point)
def _calculate_legendre_vector(self, eval_point):
......@@ -33,10 +54,49 @@ class Legendre(Vector):
class OrthonormalLegendre(Legendre):
def get_vector(self, eval_point):
def _set_basis_vector(self, eval_point):
leg_vector = self._calculate_legendre_vector(eval_point)
return [leg_vector[degree] * np.sqrt(degree+0.5) for degree in range(self._polynom_degree+1)]
def _set_wavelet_vector(self, eval_point):
alpert = AlpertsWavelet(self._polynom_degree)
return alpert.get_vector(eval_point)
def get_basis_matrices(self):
matrix_left = self._set_basis_matrix(xi, 0.5 * (xi - 1))
matrix_right = self._set_basis_matrix(xi, 0.5 * (xi + 1))
return matrix_left, matrix_right
def _set_basis_matrix(self, first_param, second_param):
matrix = []
for i in range(self._polynom_degree + 1):
row = []
for j in range(self._polynom_degree + 1):
entry = integrate(self._basis[i].subs(x, first_param) *
self._basis[j].subs(x, second_param),
(xi, -1, 1))
row.append(np.float64(entry))
matrix.append(row)
return matrix
def get_multiwavelet_matrices(self):
M1 = self._set_multiwavelet_matrix(xi, -0.5*(xi-1), True)
M2 = self._set_multiwavelet_matrix(xi, 0.5*(xi+1), False)
return M1, M2
def _set_multiwavelet_matrix(self, first_param, second_param, is_M1):
matrix = []
for i in range(self._polynom_degree+1):
row = []
for j in range(self._polynom_degree+1):
entry = integrate(self._basis[i].subs(x, first_param) * self._wavelet[j].subs(x, second_param),
(xi, -1, 1))
if is_M1:
entry = entry * (-1)**(j + self._polynom_degree + 1)
row.append(np.float64(entry))
matrix.append(row)
return matrix
class AlpertsWavelet(Vector):
def get_vector(self, eval_point):
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment