diff --git a/DG_Approximation.py b/DG_Approximation.py
index 8571c5d208b2a6a1d8f2bdb303053793a02c6377..9e6d1a6ef9fdd0aa15a635147af9878268f1a5c0 100644
--- a/DG_Approximation.py
+++ b/DG_Approximation.py
@@ -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))
diff --git a/Troubled_Cell_Detector.py b/Troubled_Cell_Detector.py
index e88600f672f76d764b61047f47fa5a9febe6c66a..9bf9751f30366d4d8e2d316d864bf8430d42810e 100644
--- a/Troubled_Cell_Detector.py
+++ b/Troubled_Cell_Detector.py
@@ -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])
diff --git a/Vectors_of_Polynomials.py b/Vectors_of_Polynomials.py
index 60bc2a7758129a5f4528a5314c9714e63a41af8a..4c5d9a8023a13cdf9524c5a2f589cee76dbf33ec 100644
--- a/Vectors_of_Polynomials.py
+++ b/Vectors_of_Polynomials.py
@@ -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):