diff --git a/DG_Approximation.py b/DG_Approximation.py index 9e6d1a6ef9fdd0aa15a635147af9878268f1a5c0..ec4a529b1dac09b5f583a98814bda550a0ec13f4 100644 --- a/DG_Approximation.py +++ b/DG_Approximation.py @@ -10,15 +10,9 @@ 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 -> Done -TODO: Restructure Vectors_of_Polynomials +TODO: Restructure Vectors_of_Polynomials -> Done 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 @@ -32,7 +26,6 @@ import Limiter import Quadrature import Update_Scheme from Vectors_of_Polynomials import OrthonormalLegendre -from Vectors_of_Polynomials import AlpertsWavelet x = Symbol('x') xi = Symbol('z') @@ -166,7 +159,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) + basis_vector = self._basis.get_vector() for cell in range(self._num_grid_cells): new_row = [] diff --git a/Troubled_Cell_Detector.py b/Troubled_Cell_Detector.py index 9bf9751f30366d4d8e2d316d864bf8430d42810e..6d5636fd5af0d2054b1f14aa669969649b2b8089 100644 --- a/Troubled_Cell_Detector.py +++ b/Troubled_Cell_Detector.py @@ -11,8 +11,6 @@ import numpy as np import matplotlib.pyplot as plt from sympy import Symbol, integrate -from Vectors_of_Polynomials import OrthonormalLegendre, AlpertsWavelet - x = Symbol('x') xi = Symbol('z') @@ -137,7 +135,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.get_vector(x) + basis = self._basis.get_vector() basis_matrix = [[basis[degree].subs(x, points[point]) for point in range(num_points)] for degree in range(self._polynom_degree+1)] @@ -221,8 +219,8 @@ class WaveletDetector(TroubledCellDetector): fine_projection = projection[:, 1:-1] 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) + basis = self._basis.get_vector() + wavelet = self._basis.get_wavelet_vector() ######################################################################################################################## # For later consideration diff --git a/Vectors_of_Polynomials.py b/Vectors_of_Polynomials.py index 4c5d9a8023a13cdf9524c5a2f589cee76dbf33ec..8663c759f9cb8b71a6f80ef667ca59dedd0ea6f3 100644 --- a/Vectors_of_Polynomials.py +++ b/Vectors_of_Polynomials.py @@ -14,14 +14,17 @@ 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) + self._wavelet = self._set_wavelet_vector(xi) - def get_vector(self, eval_points): - pass + def get_vector(self): + return self._basis def _set_basis_vector(self, eval_point): return [] + def get_wavelet_vector(self): + return self._wavelet + def _set_wavelet_vector(self, eval_point): return [] @@ -33,9 +36,6 @@ class Vector(object): 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) @@ -59,8 +59,35 @@ class OrthonormalLegendre(Legendre): 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) + degree = self._polynom_degree + + if degree == 0: + return [np.sqrt(0.5) + eval_point*0] + if degree == 1: + return [np.sqrt(1.5) * (-1 + 2*eval_point), np.sqrt(0.5) * (-2 + 3*eval_point)] + if degree == 2: + return [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.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/34) * (-16 + 105*eval_point - 192*(eval_point**2) + 105*(eval_point**3))] + if degree == 4: + return [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 get_basis_matrices(self): matrix_left = self._set_basis_matrix(xi, 0.5 * (xi - 1)) @@ -89,7 +116,7 @@ class OrthonormalLegendre(Legendre): 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), + entry = integrate(self._basis[i].subs(x, first_param) * self._wavelet[j].subs(xi, second_param), (xi, -1, 1)) if is_M1: entry = entry * (-1)**(j + self._polynom_degree + 1) @@ -97,35 +124,3 @@ class OrthonormalLegendre(Legendre): matrix.append(row) return matrix - -class AlpertsWavelet(Vector): - def get_vector(self, eval_point): - degree = self._polynom_degree - - if degree == 0: - return [np.sqrt(0.5) + eval_point*0] - if degree == 1: - return [np.sqrt(1.5) * (-1 + 2*eval_point), np.sqrt(0.5) * (-2 + 3*eval_point)] - if degree == 2: - return [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.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/34) * (-16 + 105*eval_point - 192*(eval_point**2) + 105*(eval_point**3))] - if degree == 4: - return [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')