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

Adjusted line breaks. Removed redundant passes and variables. Shortened True/False returns.

parent c3efcb99
No related branches found
No related tags found
No related merge requests found
......@@ -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.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
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,12 +244,9 @@ 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])
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))
......@@ -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,17 +280,14 @@ 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
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)
......@@ -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)):
......
......@@ -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)
......@@ -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()
......@@ -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()
......
......@@ -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
......
......@@ -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
......@@ -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 = 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))
......@@ -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,8 +27,7 @@ class Legendre(Vector):
if degree == 1:
vector.append(eval_point)
else:
poly = (2.0*degree - 1) / degree\
* eval_point * vector[len(vector)-1]\
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))]
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')
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment