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

Restructured 'Vectors_of_Polynomials'.

parent dc0d74f0
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
DG_Approximation.py
+ 2
− 9
View file @ f0ff4436
......@@ -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 = []
......
Troubled_Cell_Detector.py
+ 3
− 5
View file @ f0ff4436
......@@ -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
......
Vectors_of_Polynomials.py
+ 36
− 41
View file @ f0ff4436
......@@ -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')
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment