diff --git a/DG_Approximation.py b/DG_Approximation.py
index 219b88ef63e5045113fe68b299d25eda9e0cb6ee..ed405c52ffbdd53807d159dab7052cb802c3562b 100644
--- a/DG_Approximation.py
+++ b/DG_Approximation.py
@@ -26,7 +26,10 @@ TODO: Add a verbose option
TODO: Check whether consistency is given/possible for each class instance
TODO: Make projection local variable -> Done
-TODO: Check whether current and original projection in Update_Scheme are identical
+TODO: Check whether current and original projection in Update_Scheme are identical -> Done (No)
+TODO: Adjust line breaks to standard -> Done
+TODO: Remove redundant passes -> Done
+TODO: Shorten returns for True/False and redundant variables -> Done
"""
import numpy as np
@@ -76,20 +79,15 @@ class DGScheme(object):
# Make sure all classes actually exist
if not hasattr(Troubled_Cell_Detector, self.detector):
- raise ValueError('Invalid detector: "%s"'
- % self.detector)
+ raise ValueError('Invalid detector: "%s"' % self.detector)
if not hasattr(Initial_Condition, self.init_cond):
- raise ValueError('Invalid initial condition: "%s"'
- % self.init_cond)
+ raise ValueError('Invalid initial condition: "%s"' % self.init_cond)
if not hasattr(Limiter, self.limiter):
- raise ValueError('Invalid limiter: "%s"'
- % self.limiter)
+ raise ValueError('Invalid limiter: "%s"' % self.limiter)
if not hasattr(Quadrature, self.quadrature):
- raise ValueError('Invalid quadrature: "%s"'
- % self.quadrature)
+ raise ValueError('Invalid quadrature: "%s"' % self.quadrature)
if not hasattr(Update_Scheme, self.update_scheme):
- raise ValueError('Invalid update scheme: "%s"'
- % self.update_scheme)
+ raise ValueError('Invalid update scheme: "%s"' % self.update_scheme)
# Adjust config for detectors using wavelets
if self.detector in ['Boxplot', 'Theoretical']:
@@ -99,20 +97,15 @@ class DGScheme(object):
# Replace the string names with the actual class instances
# (and add the instance variables for the quadrature)
- self.detector = getattr(Troubled_Cell_Detector, self.detector)(
- self.detector_config)
- self.init_cond = getattr(Initial_Condition, self.init_cond)(
- self.left_bound, self.right_bound, self.init_config)
- self.limiter = getattr(Limiter, self.limiter)(
- self.limiter_config)
- self.quadrature = getattr(Quadrature, self.quadrature)(
- self.quadrature_config)
- self.update_scheme = getattr(Update_Scheme, self.update_scheme)(
- self.detector, self.limiter, self.init_cond, self.mesh,
- self.wave_speed, self.polynom_degree, self.num_grid_cells,
- self.final_time, self.history_threshold, self.left_bound,
- self.right_bound)
- pass
+ self.detector = getattr(Troubled_Cell_Detector, self.detector)(self.detector_config)
+ self.init_cond = getattr(Initial_Condition, self.init_cond)(self.left_bound, self.right_bound, self.init_config)
+ self.limiter = getattr(Limiter, self.limiter)(self.limiter_config)
+ self.quadrature = getattr(Quadrature, self.quadrature)(self.quadrature_config)
+ self.update_scheme = getattr(Update_Scheme, self.update_scheme)(self.detector, self.limiter, self.init_cond,
+ self.mesh, self.wave_speed, self.polynom_degree,
+ self.num_grid_cells, self.final_time,
+ self.history_threshold, self.left_bound,
+ self.right_bound)
def approximate(self):
"""
@@ -129,8 +122,7 @@ class DGScheme(object):
# Adjust for last cell
if current_time+time_step > self.final_time:
time_step = self.final_time-current_time
- self.cfl_number = self.wave_speed * time_step \
- / self.num_grid_cells
+ self.cfl_number = self.wave_speed * time_step / self.num_grid_cells
# Update projection
projection, troubled_cells = self.update_scheme.step(projection, self.cfl_number, current_time)
@@ -154,37 +146,34 @@ class DGScheme(object):
# wavelet
def save_plots(self):
- name = self.init_cond.get_name() + '__' \
- + self.detector.get_name() + '__' + self.limiter.get_name() \
- + '__' + self.update_scheme.get_name() + '__' \
- + self.quadrature.get_name() + '__final_time_' \
- + str(self.final_time) + '__wave_speed_' + str(self.wave_speed) \
- + '__number_of_cells_' + str(self.num_grid_cells) \
- + '__polynomial_degree_' + str(self.polynom_degree)
-
- saveFig1 = True
- if saveFig1:
+ name = self.init_cond.get_name() + '__' + self.detector.get_name() + '__' + self.limiter.get_name() \
+ + '__' + self.update_scheme.get_name() + '__' + self.quadrature.get_name() + '__final_time_' \
+ + str(self.final_time) + '__wave_speed_' + str(self.wave_speed) + '__number_of_cells_' \
+ + str(self.num_grid_cells) + '__polynomial_degree_' + str(self.polynom_degree)
+
+ save_fig_1 = True
+ if save_fig_1:
plt.figure(1)
plt.savefig('testfig/exact_and_approx/' + name + '.pdf')
# plt.clf()
- saveFig2 = True
- if saveFig2:
+ save_fig_2 = True
+ if save_fig_2:
plt.figure(2)
plt.savefig('testfig/semilog_error/' + name + '.pdf')
- saveFig3 = True
- if saveFig3:
+ save_fig_3 = True
+ if save_fig_3:
plt.figure(3)
plt.savefig('testfig/error/' + name + '.pdf')
- saveFig4 = False
- if saveFig4:
+ save_fig_4 = False
+ if save_fig_4:
plt.figure(4)
plt.savefig('testfig/coeff_details/' + name + '.pdf')
- saveFig6 = True
- if saveFig6:
+ save_fig_6 = True
+ if save_fig_6:
plt.figure(6)
plt.savefig('testfig/shock_tube/' + name + '.pdf')
@@ -211,8 +200,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).get_vector(x)
# Set additional necessary config parameters
self.detector_config['num_grid_cells'] = self.num_grid_cells
@@ -221,8 +209,7 @@ class DGScheme(object):
self.limiter_config['cell_len'] = self.cell_len
# Set mesh with one ghost point on each side
- self.mesh = np.arange(self.left_bound - (3/2*self.cell_len),
- self.right_bound + (5/2*self.cell_len),
+ self.mesh = np.arange(self.left_bound - (3/2*self.cell_len), self.right_bound + (5/2*self.cell_len),
self.cell_len) # +3/2
# Set inverse mass matrix
@@ -242,9 +229,7 @@ class DGScheme(object):
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))
+ 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
@@ -259,14 +244,11 @@ class DGScheme(object):
# former to line above: currentX
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])
- * self.quadrature.get_weights()[point]
- for point in range(self.quadrature.get_num_points()))
+ 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])
+ * self.quadrature.get_weights()[point]
+ for point in range(self.quadrature.get_num_points()))
new_row.append(np.float64(new_entry))
new_row = np.array(new_row)
@@ -284,8 +266,7 @@ class DGScheme(object):
num_points = self.quadrature.get_num_points()
basis = self.basis
- basis_matrix = [[basis[degree].subs(x, points[point])
- for point in range(num_points)]
+ basis_matrix = [[basis[degree].subs(x, points[point]) for point in range(num_points)]
for degree in range(self.polynom_degree+1)]
approx = [[sum(projection[degree][cell] * basis_matrix[degree][point]
@@ -299,18 +280,15 @@ class DGScheme(object):
# print('Exact: ', len(mesh), ' ', cell_len)
grid = []
exact = []
- num_periods = np.floor(self.wave_speed * self.final_time
- / self.interval_len)
+ num_periods = np.floor(self.wave_speed * self.final_time / self.interval_len)
for cell in range(len(mesh)):
- eval_points = mesh[cell] + cell_len/2 \
- * self.quadrature.get_eval_points()
+ eval_points = mesh[cell] + cell_len/2 * self.quadrature.get_eval_points()
eval_values = []
for point in range(len(eval_points)):
- new_entry = self.init_cond.calculate(
- eval_points[point] - self.wave_speed*self.final_time
- + num_periods*self.interval_len)
+ new_entry = self.init_cond.calculate(eval_points[point] - self.wave_speed*self.final_time
+ + num_periods*self.interval_len)
eval_values.append(new_entry)
grid.append(eval_points)
@@ -326,8 +304,7 @@ class DGScheme(object):
basis_matrix_left = self._set_basis_matrix(xi, 0.5*(xi-1))
basis_matrix_right = self._set_basis_matrix(xi, 0.5*(xi+1))
coarse_cell_len = 2*self.cell_len
- coarse_mesh = np.arange(self.left_bound - (0.5*coarse_cell_len),
- self.right_bound + (1.5*coarse_cell_len),
+ coarse_mesh = np.arange(self.left_bound - (0.5*coarse_cell_len), self.right_bound + (1.5*coarse_cell_len),
coarse_cell_len)
# Calculate projection on coarse mesh
@@ -336,8 +313,7 @@ class DGScheme(object):
output_matrix = []
for i in range(int(len(transposed_vector)/2)):
# print(transposed_vector[2*i].shape)
- new_entry = 0.5*(transposed_vector[2*i] @ basis_matrix_left
- + transposed_vector[2*i+1] @ basis_matrix_right)
+ new_entry = 0.5*(transposed_vector[2*i] @ basis_matrix_left + transposed_vector[2*i+1] @ basis_matrix_right)
output_matrix.append(new_entry)
# print(new_entry)
coarse_projection = np.transpose(np.array(output_matrix))
@@ -363,11 +339,9 @@ class DGScheme(object):
print(coarse_projection == coarse_projection1)
# Plot exact and approximate solutions for coarse mesh
- grid, exact = self._calculate_exact_solution(coarse_mesh[1:-1],
- coarse_cell_len)
+ grid, exact = self._calculate_exact_solution(coarse_mesh[1:-1], coarse_cell_len)
approx = self._calculate_approximate_solution(coarse_projection)
- self._plot_solution_and_approx(grid, exact, approx,
- color_exact, color_approx)
+ self._plot_solution_and_approx(grid, exact, approx, color_exact, color_approx)
# plt.legend(['Exact-Coarse', 'Approx-Coarse'])
def _plot_fine_mesh(self, projection, color_exact, color_approx):
@@ -377,10 +351,8 @@ class DGScheme(object):
pointwise_error = np.abs(exact-approx)
max_error = np.max(pointwise_error)
- self._plot_solution_and_approx(grid, exact, approx,
- color_exact, color_approx)
- plt.legend(['Exact (Coarse)', 'Approx (Coarse)',
- 'Exact (Fine)', 'Approx (Fine)'])
+ self._plot_solution_and_approx(grid, exact, approx, color_exact, color_approx)
+ plt.legend(['Exact (Coarse)', 'Approx (Coarse)', 'Exact (Fine)', 'Approx (Fine)'])
self._plot_semilog_error(grid, pointwise_error)
self._plot_error(grid, exact, approx)
@@ -390,8 +362,7 @@ class DGScheme(object):
return approx
- def _plot_solution_and_approx(self, grid, exact, approx,
- color_exact, color_approx):
+ def _plot_solution_and_approx(self, grid, exact, approx, color_exact, color_approx):
plt.figure(1)
plt.plot(grid[0], exact[0], color_exact)
plt.plot(grid[0], approx[0], color_approx)
@@ -436,8 +407,7 @@ class DGScheme(object):
def _plot_shock_tube(self):
time_history = self.update_scheme.get_time_history()
- troubled_cell_history = \
- self.update_scheme.get_troubled_cell_history()
+ troubled_cell_history = self.update_scheme.get_troubled_cell_history()
plt.figure(6)
for pos in range(len(time_history)):
diff --git a/Initial_Condition.py b/Initial_Condition.py
index d154f70938e4c06c641397f9bdaf9c5b590f0120..5c22b45c33b460398f66275ddff038e4771ec108 100644
--- a/Initial_Condition.py
+++ b/Initial_Condition.py
@@ -15,7 +15,6 @@ class InitialCondition(object):
self.interval_len = self.right_bound-self.left_bound
self.function_name = 'None'
- pass
def get_name(self):
return self.function_name
@@ -34,11 +33,11 @@ class InitialCondition(object):
class Sine(InitialCondition):
def __init__(self, left_bound, right_bound, config):
super().__init__(left_bound, right_bound, config)
+
# Set name of function
self.function_name = 'Sine'
self.factor = config.pop('factor', 2)
- pass
def _get_point(self, x):
return np.sin(self.factor * np.pi * x)
@@ -47,9 +46,9 @@ class Sine(InitialCondition):
class Box(InitialCondition):
def __init__(self, left_bound, right_bound, config):
super().__init__(left_bound, right_bound, config)
+
# Set name of function
self.function_name = 'Box'
- pass
def _get_point(self, x):
if x < -1:
@@ -63,9 +62,9 @@ class Box(InitialCondition):
class FourPeakWave(InitialCondition):
-
def __init__(self, left_bound, right_bound, config):
super().__init__(left_bound, right_bound, config)
+
# Set name of function
self.function_name = 'FourPeakWave'
@@ -74,21 +73,16 @@ class FourPeakWave(InitialCondition):
self.beta = np.log(2) / (36 * self.delta**2)
self.a = 0.5
self.z = -0.7
- pass
def _get_point(self, x):
if (x >= -0.8) & (x <= -0.6):
- return 1/6 * (self._G(x, self.z-self.delta)
- + self._G(x, self.z+self.delta)
- + 4 * self._G(x, self.z))
+ return 1/6 * (self._G(x, self.z-self.delta) + self._G(x, self.z+self.delta) + 4 * self._G(x, self.z))
if (x >= -0.4) & (x <= -0.2):
return 1
if (x >= 0) & (x <= 0.2):
return 1 - abs(10 * (x-0.1))
if (x >= 0.4) & (x <= 0.6):
- return 1/6 * (self._F(x, self.a-self.delta)
- + self._F(x, self.a+self.delta)
- + 4 * self._F(x, self.a))
+ return 1/6 * (self._F(x, self.a-self.delta) + self._F(x, self.a+self.delta) + 4 * self._F(x, self.a))
return 0
def _G(self, x, z):
@@ -99,45 +93,41 @@ class FourPeakWave(InitialCondition):
class Linear(InitialCondition):
-
def __init__(self, left_bound, right_bound, config):
super().__init__(left_bound, right_bound, config)
+
# Set name of function
self.function_name = 'Linear'
self.factor = config.pop('factor', 1)
- pass
def _get_point(self, x):
return self.factor * x
class LinearAbsolut(InitialCondition):
-
def __init__(self, left_bound, right_bound, config):
super().__init__(left_bound, right_bound, config)
+
# Set name of function
self.function_name = 'LinearAbsolut'
self.factor = config.pop('factor', 1)
- pass
def _get_point(self, x):
return self.factor * abs(x)
class DiscontinuousConstant(InitialCondition):
-
def __init__(self, left_bound, right_bound, config):
super().__init__(left_bound, right_bound, config)
+
# Set name of function
self.function_name = 'DiscontinuousConstant'
self.x0 = config.pop('x0', 0)
self.left_factor = config.pop('left_factor', 1)
self.right_factor = config.pop('right_factor', 0.5)
- pass
def _get_point(self, x):
- return self.left_factor * (x <= self.x0) \
- + self.right_factor * (x > self.x0)
+ return self.left_factor * (x <= self.x0) + self.right_factor * (x > self.x0)
diff --git a/Limiter.py b/Limiter.py
index 44900e83e630fb955d042e87d0f5358cebafada0..34adfb364743312d0a866bc07314a4d421da58d0 100644
--- a/Limiter.py
+++ b/Limiter.py
@@ -9,7 +9,6 @@ import numpy as np
class Limiter(object):
def __init__(self, config):
self.function_name = 'None'
- pass
def get_name(self):
return self.function_name
@@ -22,7 +21,6 @@ class NoLimiter(Limiter):
def __init__(self, config):
# Set name of function
self.function_name = 'NoLimiter'
- pass
def apply(self, projection, cell):
return np.transpose(projection)[cell]
@@ -41,7 +39,6 @@ class MinMod(Limiter):
self.projection = []
self.cell = 0
- pass
def apply(self, projection, cell):
self.projection = projection
@@ -62,25 +59,17 @@ class MinMod(Limiter):
slope = []
transposed_projection = np.transpose(self.projection)
for cell in range(len(transposed_projection)):
- new_entry = sum(transposed_projection[cell][degree]
- * (degree+0.5)**0.5
+ new_entry = sum(transposed_projection[cell][degree] * (degree+0.5)**0.5
for degree in range(1, len(self.projection)))
slope.append(new_entry)
self.cell_slope = slope[self.cell]
def _determine_modification(self):
- forward_slope = (self.projection[0][self.cell+1]
- - self.projection[0][self.cell]) * (0.5**0.5)
- backward_slope = (self.projection[0][self.cell]
- - self.projection[0][self.cell-1]) * (0.5**0.5)
+ forward_slope = (self.projection[0][self.cell+1] - self.projection[0][self.cell]) * (0.5**0.5)
+ backward_slope = (self.projection[0][self.cell] - self.projection[0][self.cell-1]) * (0.5**0.5)
- if ((forward_slope <= 0) & (backward_slope <= 0)
- & (self.cell_slope <= 0)):
- return True
- if ((forward_slope >= 0) & (backward_slope >= 0)
- & (self.cell_slope >= 0)):
- return True
- return False
+ return (forward_slope <= 0) & (backward_slope <= 0) & (self.cell_slope <= 0) \
+ | (forward_slope >= 0) & (backward_slope >= 0) & (self.cell_slope >= 0)
class ModifiedMinMod(MinMod):
@@ -95,10 +84,9 @@ class ModifiedMinMod(MinMod):
self.projection = []
self.cell = 0
- pass
def _determine_modification(self):
if abs(self.cell_slope) <= self.mod_factor*self.cell_len**2:
- return True # former: u_tilde_entry
+ return True
return super()._determine_modification()
diff --git a/Main.py b/Main.py
index 8796233f006c351f7306f3356a5a964c9ba42752..ffabcebb1d688cfcbbe69593dc6ceabe92010941 100644
--- a/Main.py
+++ b/Main.py
@@ -37,19 +37,14 @@ quadrature_config['num_eval_points'] = 12 # 12
detector = 'Theoretical'
update_scheme = 'SSPRK3'
-init_cond = 'InitialCondition2'
+init_cond = 'Box'
limiter = 'ModifiedMinMod'
quadrature = 'Gauss'
-dg_scheme = DGScheme(detector, detector_config=detector_config,
- init_cond=init_cond, init_config=init_config,
- limiter=limiter, limiter_config=limiter_config,
- quadrature=quadrature,
- quadrature_config=quadrature_config,
- update_scheme=update_scheme, wave_speed=alpha,
- polynom_degree=p, cfl_number=cfl, num_grid_cells=N,
- final_time=finalTime,
- show_plot=True, # apply_limiter=True,
- left_bound=xL, right_bound=xR)
+dg_scheme = DGScheme(detector, detector_config=detector_config, init_cond=init_cond, init_config=init_config,
+ limiter=limiter, limiter_config=limiter_config, quadrature=quadrature,
+ quadrature_config=quadrature_config, update_scheme=update_scheme, wave_speed=alpha,
+ polynom_degree=p, cfl_number=cfl, num_grid_cells=N, final_time=finalTime,
+ show_plot=True, left_bound=xL, right_bound=xR)
# __, __, troubled_cells, __ =
dg_scheme.approximate()
diff --git a/Quadrature.py b/Quadrature.py
index 73424302882337388533ec25129b4e843c7c990e..36f191181bbff2842a8fa005ef4d256f20de1ed7 100644
--- a/Quadrature.py
+++ b/Quadrature.py
@@ -12,7 +12,6 @@ class Quadrature(object):
self.eval_points = None
self.weights = None
self.num_eval_points = None
- pass
def get_name(self):
return self.function_name
diff --git a/Troubled_Cell_Detector.py b/Troubled_Cell_Detector.py
index 7b0a4d8d3dbfe1840f5b6486f9baf29084b9f8e1..44d2ccbbbe6f38521cfbb9ba27f2c8831eae94aa 100644
--- a/Troubled_Cell_Detector.py
+++ b/Troubled_Cell_Detector.py
@@ -18,7 +18,6 @@ xi = Symbol('z')
class TroubledCellDetector(object):
def __init__(self, config):
self.function_name = 'None'
- pass
def get_name(self):
return self.function_name
@@ -31,7 +30,6 @@ class NoDetection(TroubledCellDetector):
def __init__(self, config):
# Set name of function
self.function_name = 'NoDetection'
- pass
def get_cells(self, projection):
return []
@@ -43,8 +41,7 @@ class WaveletDetector(TroubledCellDetector):
self.polynom_degree = config.pop('polynom_degree')
# Set fixed basis and wavelet vectors
- self.basis = OrthonormalLegendre(
- self.polynom_degree).get_vector(x)
+ self.basis = OrthonormalLegendre(self.polynom_degree).get_vector(x)
self.wavelet = AlpertsWavelet(self.polynom_degree).get_vector(x)
# Set additional necessary parameter
@@ -56,8 +53,7 @@ class WaveletDetector(TroubledCellDetector):
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(x, second_param),
(xi, -1, 1))
if is_M1:
entry = entry * (-1)**(j + self.polynom_degree + 1)
@@ -69,8 +65,7 @@ class WaveletDetector(TroubledCellDetector):
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)
+ 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))
@@ -106,8 +101,7 @@ class Boxplot(WaveletDetector):
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)]
+ indexed_coeffs = [[multiwavelet_coeffs[0, i], i]for i in range(self.num_grid_cells)]
if self.num_grid_cells < self.fold_len:
self.fold_len = self.num_grid_cells
@@ -116,24 +110,19 @@ class Boxplot(WaveletDetector):
troubled_cells = []
for fold in range(num_folds):
- sorted_fold = sorted(indexed_coeffs[fold * self.fold_len:
- (fold+1) * self.fold_len])
+ sorted_fold = sorted(indexed_coeffs[fold * self.fold_len:(fold+1) * self.fold_len])
# fold_len/4.0 had comment: int((P+3)/2)/2.0
balance_factor = self.fold_len/4.0 % 1 # former: int(...)
boundary_index = int(self.fold_len/4.0 - balance_factor)
- first_quartile = (1-balance_factor)\
- * sorted_fold[boundary_index-1][0]\
+ first_quartile = (1-balance_factor) * sorted_fold[boundary_index-1][0] \
+ balance_factor * sorted_fold[boundary_index][0]
- third_quartile = (1-balance_factor)\
- * sorted_fold[3*boundary_index-1][0]\
+ third_quartile = (1-balance_factor) * sorted_fold[3*boundary_index-1][0]\
+ balance_factor * sorted_fold[3*boundary_index][0]
- lower_bound = first_quartile\
- - self.whisker_len * (third_quartile-first_quartile)
- upper_bound = third_quartile\
- + self.whisker_len * (third_quartile-first_quartile)
+ lower_bound = first_quartile - self.whisker_len * (third_quartile-first_quartile)
+ upper_bound = third_quartile + self.whisker_len * (third_quartile-first_quartile)
# Check for lower extreme outliers and add respective cells
for cell in sorted_fold:
@@ -162,8 +151,7 @@ class Theoretical(WaveletDetector):
# Unpack necessary configurations
self.num_grid_cells = config.pop('num_grid_cells')
self.cell_len = config.pop('cell_len')
- self.cutoff_factor = config.pop('cutoff_factor',
- np.sqrt(2) * self.cell_len)
+ self.cutoff_factor = config.pop('cutoff_factor', np.sqrt(2) * self.cell_len)
# comment to line above: or 2 or 3
self.multiwavelet_coeffs = []
@@ -174,8 +162,7 @@ class Theoretical(WaveletDetector):
self.multiwavelet_coeffs = multiwavelet_coeffs
self.projection = projection
troubled_cells = []
- self.max_avg = max(1, max(abs(self.projection[0][degree])
- for degree in range(self.polynom_degree+1)))
+ self.max_avg = max(1, max(abs(self.projection[0][degree]) for degree in range(self.polynom_degree+1)))
for cell in range(int(self.num_grid_cells/2)):
if self._is_troubled_cell(cell):
@@ -185,10 +172,7 @@ class Theoretical(WaveletDetector):
def _is_troubled_cell(self, cell):
max_value = max(abs(self.multiwavelet_coeffs[degree][cell])
- for degree in range(self.polynom_degree+1)) \
- / self.max_avg
+ for degree in range(self.polynom_degree+1))/self.max_avg
eps = self.cutoff_factor / (self.cell_len*self.num_grid_cells)
- if max_value > eps:
- return True
- return False
+ return max_value > eps
diff --git a/Update_Scheme.py b/Update_Scheme.py
index 135130733c03c73f0a5a098ef26129b791a309ef..f74a836d3d1b6b5b1eb724214de7c79e95194015 100644
--- a/Update_Scheme.py
+++ b/Update_Scheme.py
@@ -16,9 +16,8 @@ import numpy as np
class UpdateScheme(object):
- def __init__(self, detector, limiter, init_cond, mesh, wave_speed,
- polynom_degree, num_grid_cells, final_time, history_threshold,
- left_bound, right_bound):
+ def __init__(self, detector, limiter, init_cond, mesh, wave_speed, polynom_degree, num_grid_cells, final_time,
+ history_threshold, left_bound, right_bound):
# Unpack positional arguments
self.detector = detector
self.limiter = limiter
@@ -33,7 +32,6 @@ class UpdateScheme(object):
self.right_bound = right_bound
self._reset()
- pass
def get_name(self):
return self.name
@@ -107,8 +105,7 @@ class UpdateScheme(object):
new_projection = self.current_projection.copy()
for cell in self.troubled_cells:
- np.transpose(new_projection)[cell] = self.limiter.apply(
- self.current_projection, cell)
+ np.transpose(new_projection)[cell] = self.limiter.apply(self.current_projection, cell)
self.current_projection = new_projection
@@ -122,11 +119,9 @@ class UpdateScheme(object):
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):
- super().__init__(detector, limiter, init_cond, mesh, wave_speed,
- polynom_degree, num_grid_cells, final_time,
+ def __init__(self, detector, limiter, init_cond, mesh, wave_speed, polynom_degree, num_grid_cells, final_time,
+ history_threshold, left_bound, right_bound):
+ 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
@@ -166,18 +161,14 @@ class SSPRK3(UpdateScheme):
def _apply_first_step(self):
self._update_right_hand_side()
- self.current_projection = self.original_projection \
- + (self.cfl_number*self.right_hand_side)
+ self.current_projection = self.original_projection + (self.cfl_number*self.right_hand_side)
def _apply_second_step(self):
self._update_right_hand_side()
- self.current_projection = 1/4 * (
- 3*self.original_projection
- + (self.current_projection + self.cfl_number*self.right_hand_side))
+ self.current_projection = 1/4 * (3 * self.original_projection
+ + (self.current_projection + self.cfl_number*self.right_hand_side))
def _apply_third_step(self):
self._update_right_hand_side()
- self.current_projection = 1/3 * (
- self.original_projection
- + 2 * (self.current_projection
- + self.cfl_number*self.right_hand_side))
+ self.current_projection = 1/3 * (self.original_projection
+ + 2 * (self.current_projection + self.cfl_number*self.right_hand_side))
diff --git a/Vectors_of_Polynomials.py b/Vectors_of_Polynomials.py
index 43fa9608d13f6a1c4f7038f741bd426ba6202db6..767c036468f601fb5c89fec1f145c36b23a38512 100644
--- a/Vectors_of_Polynomials.py
+++ b/Vectors_of_Polynomials.py
@@ -9,7 +9,6 @@ import numpy as np
class Vector(object):
def __init__(self, polynom_degree):
self.polynom_degree = polynom_degree
- pass
def get_vector(self, eval_points):
pass
@@ -17,8 +16,7 @@ class Vector(object):
class Legendre(Vector):
def get_vector(self, eval_point):
- vector = self._calculate_legendre_vector(eval_point)
- return vector
+ return self._calculate_legendre_vector(eval_point)
def _calculate_legendre_vector(self, eval_point):
vector = []
@@ -29,9 +27,8 @@ class Legendre(Vector):
if degree == 1:
vector.append(eval_point)
else:
- poly = (2.0*degree - 1) / degree\
- * eval_point * vector[len(vector)-1]\
- - (degree-1) / degree * vector[len(vector)-2]
+ poly = (2.0*degree - 1)/degree * eval_point * vector[len(vector)-1] \
+ - (degree-1)/degree * vector[len(vector)-2]
vector.append(poly)
return vector
@@ -39,52 +36,37 @@ class Legendre(Vector):
class OrthonormalLegendre(Legendre):
def get_vector(self, eval_point):
leg_vector = self._calculate_legendre_vector(eval_point)
- vector = [leg_vector[degree] * np.sqrt(degree+0.5)
- for degree in range(len(leg_vector))]
- return vector
+ return [leg_vector[degree] * np.sqrt(degree+0.5) for degree in range(len(leg_vector))]
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)]
+ 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))]
+ 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))]
+ 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))]
+ - 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: Albert\'s wavelet is only available \
up to degree 4 for this application')