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

Moved everything pertaining wavelets to Troubled_Cell_Detector

parent 2a22685d
No related branches found
No related tags found
No related merge requests found
......@@ -25,6 +25,9 @@ TODO: Write documentation for all methods
TODO: Add a verbose option
TODO: Check whether consistency is given/possible for each class instance
TODO: Make projection local variable
TODO: Check whether current and original projection in Update_Scheme are identical
"""
import numpy as np
from sympy import Symbol, integrate
......@@ -90,6 +93,10 @@ class DGScheme(object):
raise ValueError('Invalid update scheme: "%s"'
% self.update_scheme)
# Adjust config for detectors using wavelets
if self.detector in ['Boxplot', 'Theoretical']:
self.detector_config['polynom_degree'] = self.polynom_degree
self._reset()
# Replace the string names with the actual class instances
......@@ -135,16 +142,19 @@ class DGScheme(object):
self._plot_shock_tube()
self._plot_coarse_mesh('k-', 'y')
approx = self._plot_fine_mesh('k-.', 'b')
troubled_cells = self.update_scheme.get_troubled_cell_history()
wavelet = self.update_scheme.get_wavelet_coeffs(
self.projection[:, 1: -1])
approx = self._plot_fine_mesh('k-.', 'b')
return approx
# What is that??? What is it for?
return self.projection[:, 1:-1], self.mesh[2:-2], troubled_cells, \
wavelet
# troubled_cells = self.update_scheme.get_troubled_cell_history()
# wavelet = self.update_scheme.get_wavelet_coeffs(
# self.projection[:, 1: -1])
#
# return approx
#
# return self.projection[:, 1:-1], self.mesh[2:-2], troubled_cells, \
# wavelet
def save_plots(self):
name = self.init_cond.get_name() + '__' \
......@@ -245,7 +255,7 @@ class DGScheme(object):
return matrix
def _do_initial_projection(self):
# Initialze matrix and set first entry to accommodate for ghost cell
# Initialize matrix and set first entry to accommodate for ghost cell
output_matrix = [0]
for cell in range(self.num_grid_cells):
......@@ -274,7 +284,7 @@ class DGScheme(object):
print(np.array(output_matrix).shape)
self.projection = np.transpose(np.array(output_matrix))
def _calculate_approximative_solution(self, projection):
def _calculate_approximate_solution(self, projection):
points = self.quadrature.get_eval_points()
num_points = self.quadrature.get_num_points()
basis = self.basis
......@@ -357,10 +367,10 @@ class DGScheme(object):
print(coarse_projection == coarse_projection1)
# Plot exact and approximative solutions for coarse mesh
# Plot exact and approximate solutions for coarse mesh
grid, exact = self._calculate_exact_solution(coarse_mesh[1:-1],
coarse_cell_len)
approx = self._calculate_approximative_solution(coarse_projection)
approx = self._calculate_approximate_solution(coarse_projection)
self._plot_solution_and_approx(grid, exact, approx,
color_exact, color_approx)
# plt.legend(['Exact-Coarse', 'Approx-Coarse'])
......@@ -368,7 +378,7 @@ class DGScheme(object):
def _plot_fine_mesh(self, color_exact, color_approx):
grid, exact = self._calculate_exact_solution(self.mesh[2:-2],
self.cell_len)
approx = self._calculate_approximative_solution(
approx = self._calculate_approximate_solution(
self.projection[:, 1:-1])
pointwise_error = np.abs(exact-approx)
......
......@@ -64,3 +64,6 @@ dg_scheme.save_plots()
# print(dg_scheme.limiter)
# # print(dg_scheme.projection)
# =============================================================================
# if __name__ == '__main__':
# pass
......@@ -7,16 +7,23 @@ TODO: Check cutoff_factor/whether _is_troubled_cell works correctly \
"""
import numpy as np
from sympy import Symbol, integrate
from Vectors_of_Polynomials import OrthonormalLegendre, AlpertsWavelet
x = Symbol('x')
xi = Symbol('z')
class TroubledCellDetector(object):
def __init__(self, config):
self.function_name = 'None'
pass
def get_name(self):
return self.function_name
def get_cells(self, multiwavelet_coeffs, projection):
def get_cells(self, projection):
pass
......@@ -26,12 +33,65 @@ class NoDetection(TroubledCellDetector):
self.function_name = 'NoDetection'
pass
def get_cells(self, multiwavelet_coeffs, projection):
def get_cells(self, projection):
return []
class WaveletDetector(TroubledCellDetector):
def __init__(self, config):
# Unpack necessary configurations
self.polynom_degree = config.pop('polynom_degree')
# 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.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
def _calculate_wavelet_coeffs(self, projection):
transposed_vector = np.transpose(projection)
output_matrix = []
for i in range(int(len(projection[0])/2)):
new_entry = 0.5*(transposed_vector[2*i] @ self.M1
+ transposed_vector[2*i+1] @ self.M2)
output_matrix.append(new_entry)
return np.transpose(np.array(output_matrix))
# def get_wavelet_coeffs(self, projection=None):
# if projection is None:
# return self.multiwavelet_coeffs
# else:
# return self._calculate_wavelet_coeffs(projection)
def get_cells(self, projection):
multiwavelet_coeffs = self._calculate_wavelet_coeffs(projection[:, 1: -1])
return self._get_cells(multiwavelet_coeffs, projection)
def _get_cells(self, multiwavelet_coeffs, projection):
return []
class Boxplot(TroubledCellDetector):
class Boxplot(WaveletDetector):
def __init__(self, config):
super().__init__(config)
# Set name of function
self.function_name = 'Boxplot'
......@@ -44,7 +104,7 @@ class Boxplot(TroubledCellDetector):
config.pop('polynom_degree')
config.pop('cell_len')
def get_cells(self, multiwavelet_coeffs, projection):
def _get_cells(self, multiwavelet_coeffs, projection):
self.num_grid_cells = len(multiwavelet_coeffs[0])
indexed_coeffs = [[multiwavelet_coeffs[0, i], i]
for i in range(self.num_grid_cells)]
......@@ -92,13 +152,14 @@ class Boxplot(TroubledCellDetector):
return sorted(troubled_cells)
class Theoretical(TroubledCellDetector):
class Theoretical(WaveletDetector):
def __init__(self, config):
super().__init__(config)
# Set name of function
self.function_name = 'Theoretical'
# Unpack necessary configurations
self.polynom_degree = config.pop('polynom_degree')
self.num_grid_cells = config.pop('num_grid_cells')
self.cell_len = config.pop('cell_len')
self.cutoff_factor = config.pop('cutoff_factor',
......@@ -109,7 +170,7 @@ class Theoretical(TroubledCellDetector):
self.projection = []
self.max_avg = None
def get_cells(self, multiwavelet_coeffs, projection):
def _get_cells(self, multiwavelet_coeffs, projection):
self.multiwavelet_coeffs = multiwavelet_coeffs
self.projection = projection
troubled_cells = []
......
......@@ -13,12 +13,6 @@ TODO: Find better names for A, B, M1, and M2
"""
import numpy as np
from sympy import Symbol, integrate
from Vectors_of_Polynomials import OrthonormalLegendre, AlpertsWavelet
x = Symbol('x')
xi = Symbol('z')
class UpdateScheme(object):
......@@ -50,12 +44,6 @@ class UpdateScheme(object):
def get_time_history(self):
return self.time_history
def get_wavelet_coeffs(self, projection=None):
if projection is None:
return self.multiwavelet_coeffs
else:
return self._calculate_wavelet_coeffs(projection)
def step(self, projection, cfl_number, time):
self.original_projection = projection
self.current_projection = projection
......@@ -76,17 +64,10 @@ class UpdateScheme(object):
pass
def _reset(self):
# 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 fixed instance variables
self.name = 'None'
self.interval_len = self.right_bound-self.left_bound
self.cell_len = self.interval_len / self.num_grid_cells
self.M1 = self._set_multiwavelet_matrix(xi, -0.5*(xi-1), True)
self.M2 = self._set_multiwavelet_matrix(xi, 0.5*(xi+1), False)
# Set matrix A
matrix = []
......@@ -114,7 +95,6 @@ class UpdateScheme(object):
self.original_projection = []
self.current_projection = []
self.right_hand_side = []
self.multiwavelet_coeffs = []
self.troubled_cells = []
self.troubled_cell_history = []
self.time_history = []
......@@ -122,37 +102,8 @@ class UpdateScheme(object):
self.time = 0
self.iteration = 0
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
def _update_wavelet_coeffs(self):
projection = self.current_projection[:, 1: -1]
self.multiwavelet_coeffs = self._calculate_wavelet_coeffs(projection)
def _calculate_wavelet_coeffs(self, projection):
transposed_vector = np.transpose(projection)
output_matrix = []
for i in range(int(len(projection[0])/2)):
new_entry = 0.5*(transposed_vector[2*i] @ self.M1
+ transposed_vector[2*i+1] @ self.M2)
output_matrix.append(new_entry)
return np.transpose(np.array(output_matrix))
def _apply_limiter(self):
self._update_wavelet_coeffs()
self.troubled_cells = self.detector.get_cells(self.multiwavelet_coeffs,
self.current_projection)
self.troubled_cells = self.detector.get_cells(self.current_projection)
new_projection = self.current_projection.copy()
for cell in self.troubled_cells:
......@@ -174,13 +125,13 @@ class SSPRK3(UpdateScheme):
def __init__(self, detector, limiter, init_cond, mesh, wave_speed,
polynom_degree, num_grid_cells, final_time, history_threshold,
left_bound, right_bound):
# Set name of update scheme
self.name = 'SSPRK3'
super().__init__(detector, limiter, init_cond, mesh, wave_speed,
polynom_degree, num_grid_cells, final_time,
history_threshold, left_bound, right_bound)
# Set name of update scheme
self.name = 'SSPRK3'
# Override method of superclass
def _apply_stability_method(self):
self._apply_first_step()
......@@ -199,7 +150,7 @@ class SSPRK3(UpdateScheme):
# Transpose projection for easier calculation
transposed_projection = np.transpose(self.current_projection)
# Initialze vector and set first entry to accommodate for ghost cell
# Initialize vector and set first entry to accommodate for ghost cell
right_hand_side = [0]
for j in range(self.num_grid_cells):
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment