diff --git a/DG_Approximation.py b/DG_Approximation.py
index 1b666e46acdec5474f362f0f3e3cc18f319ca9cd..b31aa5218ace8af33c000263ad08b5f3f06bc0f6 100644
--- a/DG_Approximation.py
+++ b/DG_Approximation.py
@@ -12,20 +12,21 @@ TODO: Combine initial projection and approx solution somehow
TODO: Investigate why there are no weights in approx calc -> Ask
TODO: Change order of methods
TODO: Change code to not include self. but global variables; NO, do private \
- or better protected instance variables instead
+ or better protected instance variables instead -> Done (all protected)
TODO: Find better names for A, B, anything pertaining Gauss
TODO: Write documentation for all methods
-TODO: Check whether consistency is given/possible for each class instance
+TODO: Check whether consistency is given/possible for each class instance -> Done (should be okay when protected)
TODO: Contemplate moving A and B to Vectors_of_Polynomials
TODO: Improve access to names of classes
TODO: Check plots for NoDetection class -> Done (Adjusted scale of x-axis for shock tube)
TODO: Ask about shock tubes (course mesh vs fine mesh)
-TODO: Check instance variables in DG_Approximation -> Done
+TODO: Check instance variables in DG_Approximation -> Done (do for all else)
TODO: Improve saving of plots
TODO: Extend color options
TODO: Implement type check for all kwargs and configs
TODO: Make sure that self.cfl_number is not changed -> Done
+TODO: Contemplate having a dict with a params and just getting the necessary ones
"""
import numpy as np
@@ -50,24 +51,24 @@ xi = Symbol('z')
class DGScheme(object):
def __init__(self, detector, **kwargs):
# Unpack keyword arguments
- self.wave_speed = kwargs.pop('wave_speed', 1)
- self.polynom_degree = kwargs.pop('polynom_degree', 2)
- self.cfl_number = kwargs.pop('cfl_number', 0.2)
- self.num_grid_cells = kwargs.pop('num_grid_cells', 64)
- self.final_time = kwargs.pop('final_time', 1)
- self.left_bound = kwargs.pop('left_bound', -1)
- self.right_bound = kwargs.pop('right_bound', 1)
- self.verbose = kwargs.pop('verbose', False)
- self.history_threshold = kwargs.pop('history_threshold', math.ceil(0.2/self.cfl_number))
- self.detector = detector
- self.detector_config = kwargs.pop('detector_config', {})
- self.init_cond = kwargs.pop('init_cond', 'Sine')
- self.init_config = kwargs.pop('init_config', {})
- self.limiter = kwargs.pop('limiter', 'ModifiedMinMod')
- self.limiter_config = kwargs.pop('limiter_config', {})
- self.quadrature = kwargs.pop('quadrature', 'Gauss')
- self.quadrature_config = kwargs.pop('quadrature_config', {})
- self.update_scheme = kwargs.pop('update_scheme', 'SSPRK3')
+ self._wave_speed = kwargs.pop('wave_speed', 1)
+ self._polynom_degree = kwargs.pop('polynom_degree', 2)
+ self._cfl_number = kwargs.pop('cfl_number', 0.2)
+ self._num_grid_cells = kwargs.pop('num_grid_cells', 64)
+ self._final_time = kwargs.pop('final_time', 1)
+ self._left_bound = kwargs.pop('left_bound', -1)
+ self._right_bound = kwargs.pop('right_bound', 1)
+ self._verbose = kwargs.pop('verbose', False)
+ self._history_threshold = kwargs.pop('history_threshold', math.ceil(0.2/self._cfl_number))
+ self._detector = detector
+ self._detector_config = kwargs.pop('detector_config', {})
+ self._init_cond = kwargs.pop('init_cond', 'Sine')
+ self._init_config = kwargs.pop('init_config', {})
+ self._limiter = kwargs.pop('limiter', 'ModifiedMinMod')
+ self._limiter_config = kwargs.pop('limiter_config', {})
+ self._quadrature = kwargs.pop('quadrature', 'Gauss')
+ self._quadrature_config = kwargs.pop('quadrature_config', {})
+ self._update_scheme = kwargs.pop('update_scheme', 'SSPRK3')
# Throw an error if there are extra keyword arguments
if len(kwargs) > 0:
@@ -75,35 +76,35 @@ class DGScheme(object):
raise ValueError('Unrecognized arguments: %s' % extra)
# Make sure all classes actually exist
- if not hasattr(Troubled_Cell_Detector, 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)
- if not hasattr(Limiter, self.limiter):
- raise ValueError('Invalid limiter: "%s"' % self.limiter)
- if not hasattr(Quadrature, 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)
+ if not hasattr(Troubled_Cell_Detector, 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)
+ if not hasattr(Limiter, self._limiter):
+ raise ValueError('Invalid limiter: "%s"' % self._limiter)
+ if not hasattr(Quadrature, 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)
self._reset()
# Replace the string names with the actual class instances
# (and add the instance variables for the quadrature)
- 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.detector = getattr(Troubled_Cell_Detector, self.detector)(
- self.detector_config, self.mesh, self.wave_speed, self.polynom_degree, self.num_grid_cells,
- self.final_time, self.left_bound, self.right_bound, self.basis, self.init_cond, self.quadrature)
- 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)
-
- # print(self.detector.get_name())
- # print(self.detector.__class__.__name__)
+ 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._detector = getattr(Troubled_Cell_Detector, self._detector)(
+ self._detector_config, self._mesh, self._wave_speed, self._polynom_degree, self._num_grid_cells,
+ self._final_time, self._left_bound, self._right_bound, self._basis, self._init_cond, self._quadrature)
+ 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)
+
+ # print(self._detector.get_name())
+ # print(self._detector.__class__.__name__)
# print(Troubled_Cell_Detector.Boxplot.__bases__)
- # print(self.detector.__class__)
+ # print(self._detector.__class__)
# print(Troubled_Cell_Detector.WaveletDetector in Troubled_Cell_Detector.Boxplot.__mro__)
def approximate(self):
@@ -115,41 +116,41 @@ class DGScheme(object):
projection = self._do_initial_projection()
- time_step = abs(self.cfl_number * self.cell_len / self.wave_speed)
+ time_step = abs(self._cfl_number * self._cell_len / self._wave_speed)
current_time = 0
iteration = 0
troubled_cell_history = []
time_history = []
- while current_time < self.final_time:
+ while current_time < self._final_time:
# Adjust for last cell
- cfl_number = self.cfl_number
- if current_time+time_step > self.final_time:
- time_step = self.final_time-current_time
- cfl_number = self.wave_speed * time_step / self.num_grid_cells
+ cfl_number = self._cfl_number
+ if current_time+time_step > self._final_time:
+ time_step = self._final_time-current_time
+ cfl_number = self._wave_speed * time_step / self._num_grid_cells
# Update projection
- projection, troubled_cells = self.update_scheme.step(projection, cfl_number, current_time)
+ projection, troubled_cells = self._update_scheme.step(projection, cfl_number, current_time)
iteration += 1
- if (iteration % self.history_threshold) == 0:
+ if (iteration % self._history_threshold) == 0:
troubled_cell_history.append(troubled_cells)
time_history.append(current_time)
current_time += time_step
# Plot exact/approximate results, errors, shock tubes and any detector-dependant plots
- self.detector.plot_results(projection, troubled_cell_history, time_history, 'k-', 'y')
+ self._detector.plot_results(projection, troubled_cell_history, time_history, 'k-', 'y')
- if self.verbose:
+ if self._verbose:
plt.show()
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)
+ 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:
@@ -198,50 +199,50 @@ class DGScheme(object):
os.makedirs('testfig/shock_tube')
# 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._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)
# Set additional necessary config parameters
- self.limiter_config['cell_len'] = self.cell_len
+ 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.cell_len) # +3/2
+ 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
mass_matrix = []
- for i in range(self.polynom_degree+1):
+ for i in range(self._polynom_degree+1):
new_row = []
- for j in range(self.polynom_degree+1):
+ for j in range(self._polynom_degree+1):
new_entry = 0.0
if i == j:
new_entry = 1.0
new_row.append(new_entry)
mass_matrix.append(new_row)
- self.inv_mass = np.array(mass_matrix)
+ self._inv_mass = np.array(mass_matrix)
def _do_initial_projection(self):
# Initialize matrix and set first entry to accommodate for ghost cell
output_matrix = [0]
- for cell in range(self.num_grid_cells):
+ for cell in range(self._num_grid_cells):
new_row = []
- eval_point = self.left_bound + (cell+0.5)*self.cell_len
-
- 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()))
+ eval_point = self._left_bound + (cell+0.5)*self._cell_len
+
+ 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_row.append(np.float64(new_entry))
new_row = np.array(new_row)
- output_matrix.append(self.inv_mass @ new_row)
+ output_matrix.append(self._inv_mass @ new_row)
# Set ghost cells to respective value
- output_matrix[0] = output_matrix[self.num_grid_cells]
+ output_matrix[0] = output_matrix[self._num_grid_cells]
output_matrix.append(output_matrix[1])
print(np.array(output_matrix).shape)
diff --git a/Initial_Condition.py b/Initial_Condition.py
index 04c80c8fb7f1826c350046ba696f0b41678461fd..834f8fe5ac7cca0ccdef6c075a1671c1473bec45 100644
--- a/Initial_Condition.py
+++ b/Initial_Condition.py
@@ -8,9 +8,9 @@ import numpy as np
class InitialCondition(object):
def __init__(self, left_bound, right_bound, config):
- self.left_bound = left_bound
- self.right_bound = right_bound
- self.interval_len = self.right_bound-self.left_bound
+ self._left_bound = left_bound
+ self._right_bound = right_bound
+ self._interval_len = self._right_bound-self._left_bound
self.function_name = 'None'
@@ -18,10 +18,10 @@ class InitialCondition(object):
return self.function_name
def calculate(self, x):
- while x < self.left_bound:
- x = x + self.interval_len
- while x > self.right_bound:
- x = x - self.interval_len
+ while x < self._left_bound:
+ x = x + self._interval_len
+ while x > self._right_bound:
+ x = x - self._interval_len
return self._get_point(x)
def _get_point(self, x):
@@ -36,10 +36,10 @@ class Sine(InitialCondition):
self.function_name = 'Sine'
# Unpack necessary configurations
- self.factor = config.pop('factor', 2)
+ self._factor = config.pop('factor', 2)
def _get_point(self, x):
- return np.sin(self.factor * np.pi * x)
+ return np.sin(self._factor * np.pi * x)
class Box(InitialCondition):
@@ -68,28 +68,28 @@ class FourPeakWave(InitialCondition):
self.function_name = 'FourPeakWave'
# Set additional necessary parameter
- self.alpha = 10
- self.delta = 0.005
- self.beta = np.log(2) / (36 * self.delta**2)
- self.a = 0.5
- self.z = -0.7
+ self._alpha = 10
+ self._delta = 0.005
+ self._beta = np.log(2) / (36 * self._delta**2)
+ self._a = 0.5
+ self._z = -0.7
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):
- return np.exp(-self.beta * (x-z)**2)
+ return np.exp(-self._beta * (x-z)**2)
def _F(self, x, a):
- return np.sqrt(max(1 - self.alpha**2 * (x-a)**2, 0))
+ return np.sqrt(max(1 - self._alpha**2 * (x-a)**2, 0))
class Linear(InitialCondition):
@@ -100,10 +100,10 @@ class Linear(InitialCondition):
self.function_name = 'Linear'
# Unpack necessary configurations
- self.factor = config.pop('factor', 1)
+ self._factor = config.pop('factor', 1)
def _get_point(self, x):
- return self.factor * x
+ return self._factor * x
class LinearAbsolut(InitialCondition):
@@ -114,10 +114,10 @@ class LinearAbsolut(InitialCondition):
self.function_name = 'LinearAbsolut'
# Unpack necessary configurations
- self.factor = config.pop('factor', 1)
+ self._factor = config.pop('factor', 1)
def _get_point(self, x):
- return self.factor * abs(x)
+ return self._factor * abs(x)
class DiscontinuousConstant(InitialCondition):
@@ -128,9 +128,9 @@ class DiscontinuousConstant(InitialCondition):
self.function_name = 'DiscontinuousConstant'
# Unpack necessary configurations
- self.x0 = config.pop('x0', 0)
- self.left_factor = config.pop('left_factor', 1)
- self.right_factor = config.pop('right_factor', 0.5)
+ self._x0 = config.pop('x0', 0)
+ self._left_factor = config.pop('left_factor', 1)
+ self._right_factor = config.pop('right_factor', 0.5)
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 a5033e70fb44d3cb9a519fb8f90ed8a1db3b0e9a..bedb510b0bbc06b311de3e317e71e7214ed98ee9 100644
--- a/Limiter.py
+++ b/Limiter.py
@@ -32,10 +32,10 @@ class MinMod(Limiter):
super().__init__(config)
# Unpack necessary configurations
- self.erase_degree = config.pop('erase_degree', 0)
+ self._erase_degree = config.pop('erase_degree', 0)
# Set name of function
- self.function_name = 'MinMod' + str(self.erase_degree)
+ self.function_name = 'MinMod' + str(self._erase_degree)
def apply(self, projection, cell):
cell_slope = self._set_cell_slope(projection, cell)
@@ -46,7 +46,7 @@ class MinMod(Limiter):
adapted_projection = projection[:, cell].copy()
for i in range(len(adapted_projection)):
- if i > self.erase_degree:
+ if i > self._erase_degree:
adapted_projection[i] = 0
return adapted_projection
@@ -72,16 +72,16 @@ class ModifiedMinMod(MinMod):
super().__init__(config)
# Unpack necessary configurations
- self.cell_len = config.pop('cell_len')
- self.mod_factor = config.pop('mod_factor', 0)
+ self._cell_len = config.pop('cell_len')
+ self._mod_factor = config.pop('mod_factor', 0)
# Set name of function
- self.function_name = 'ModifiedMinMod' + str(self.erase_degree)
+ self.function_name = 'ModifiedMinMod' + str(self._erase_degree)
- self.cell = 0
+ self._cell = 0
def _determine_modification(self, projection, cell, cell_slope):
- if abs(cell_slope) <= self.mod_factor*self.cell_len**2:
+ if abs(cell_slope) <= self._mod_factor*self._cell_len**2:
return True
return super()._determine_modification(projection, cell, cell_slope)
diff --git a/Quadrature.py b/Quadrature.py
index f7c6b19e818d1beecb065b32e77ef47054a7b48a..9cf793eb396922a201243db614f993bb58055d04 100644
--- a/Quadrature.py
+++ b/Quadrature.py
@@ -9,21 +9,21 @@ import numpy.polynomial.legendre as leg
class Quadrature(object):
def __init__(self, config):
self.function_name = 'None'
- self.eval_points = None
- self.weights = None
- self.num_eval_points = None
+ self._eval_points = None
+ self._weights = None
+ self._num_eval_points = None
def get_name(self):
return self.function_name
def get_eval_points(self):
- return self.eval_points
+ return self._eval_points
def get_weights(self):
- return self.weights
+ return self._weights
def get_num_points(self):
- return self.num_eval_points
+ return self._num_eval_points
class Gauss(Quadrature):
@@ -31,9 +31,9 @@ class Gauss(Quadrature):
super().__init__(config)
# Unpack necessary configurations
- self.num_eval_points = config.pop('num_eval_points', 6)
+ self._num_eval_points = config.pop('num_eval_points', 6)
- self.eval_points, self.weights = leg.leggauss(self.num_eval_points)
+ self._eval_points, self._weights = leg.leggauss(self._num_eval_points)
# Set name of function
- self.function_name = 'Gauss' + str(self.num_eval_points)
+ self.function_name = 'Gauss' + str(self._num_eval_points)
diff --git a/Troubled_Cell_Detector.py b/Troubled_Cell_Detector.py
index 1c210b9d8b6be3fff66398902a0ef18e8b5509d6..327e24b97dfd0e17220e5415418d7c6135427612 100644
--- a/Troubled_Cell_Detector.py
+++ b/Troubled_Cell_Detector.py
@@ -19,18 +19,18 @@ xi = Symbol('z')
class TroubledCellDetector(object):
def __init__(self, config, mesh, wave_speed, polynom_degree, num_grid_cells, final_time, left_bound, right_bound,
basis, init_cond, quadrature):
- self.mesh = mesh
- self.wave_speed = wave_speed
- self.polynom_degree = polynom_degree
- self.num_grid_cells = num_grid_cells
- self.final_time = final_time
- self.left_bound = left_bound
- self.right_bound = right_bound
- self.interval_len = right_bound - left_bound
- self.cell_len = self.interval_len / num_grid_cells
- self.basis = basis
- self.init_cond = init_cond
- self.quadrature = quadrature
+ self._mesh = mesh
+ self._wave_speed = wave_speed
+ self._polynom_degree = polynom_degree
+ self._num_grid_cells = num_grid_cells
+ self._final_time = final_time
+ self._left_bound = left_bound
+ self._right_bound = right_bound
+ self._interval_len = right_bound - left_bound
+ self._cell_len = self._interval_len / num_grid_cells
+ self._basis = basis
+ self._init_cond = init_cond
+ self._quadrature = quadrature
self.function_name = 'None'
def get_name(self):
@@ -43,12 +43,12 @@ class TroubledCellDetector(object):
self._plot_shock_tube(troubled_cell_history, time_history)
max_error = self._plot_mesh(projection, color_exact, color_approx)
- print("p =", self.polynom_degree)
- print("N =", self.num_grid_cells)
+ print("p =", self._polynom_degree)
+ print("N =", self._num_grid_cells)
print("maximum error =", max_error)
def _plot_mesh(self, projection, color_exact, color_approx):
- grid, exact = self._calculate_exact_solution(self.mesh[2:-2], self.cell_len)
+ grid, exact = self._calculate_exact_solution(self._mesh[2:-2], self._cell_len)
approx = self._calculate_approximate_solution(projection[:, 1:-1])
pointwise_error = np.abs(exact-approx)
@@ -90,7 +90,7 @@ class TroubledCellDetector(object):
current_cells = troubled_cell_history[pos]
for cell in current_cells:
plt.plot(cell, time_history[pos], 'k.')
- plt.xlim((0, self.num_grid_cells//2))
+ plt.xlim((0, self._num_grid_cells//2))
plt.xlabel('Cell')
plt.ylabel('Time')
plt.title('Shock Tubes')
@@ -98,15 +98,15 @@ class TroubledCellDetector(object):
def _calculate_exact_solution(self, 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)
@@ -118,15 +118,15 @@ class TroubledCellDetector(object):
return grid, exact
def _calculate_approximate_solution(self, projection):
- points = self.quadrature.get_eval_points()
- num_points = self.quadrature.get_num_points()
- basis = self.basis
+ points = self._quadrature.get_eval_points()
+ num_points = self._quadrature.get_num_points()
+ basis = self._basis
basis_matrix = [[basis[degree].subs(x, points[point]) for point in range(num_points)]
- for degree in range(self.polynom_degree+1)]
+ for degree in range(self._polynom_degree+1)]
approx = [[sum(projection[degree][cell] * basis_matrix[degree][point]
- for degree in range(self.polynom_degree+1))
+ for degree in range(self._polynom_degree+1))
for point in range(num_points)]
for cell in range(len(projection[0]))]
@@ -153,23 +153,23 @@ class WaveletDetector(TroubledCellDetector):
basis, init_cond, quadrature)
# Set fixed basis and wavelet vectors
- self.basis = OrthonormalLegendre(self.polynom_degree).get_vector(x)
- self.wavelet = AlpertsWavelet(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
- self.num_coarse_grid_cells = self.num_grid_cells//2
- self.M1 = self._set_multiwavelet_matrix(xi, -0.5*(xi-1), True)
- self.M2 = self._set_multiwavelet_matrix(xi, 0.5*(xi+1), False)
+ self._num_coarse_grid_cells = self._num_grid_cells//2
+ 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):
+ 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),
+ 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)
+ entry = entry * (-1)**(j + self._polynom_degree + 1)
row.append(np.float64(entry))
matrix.append(row)
return matrix
@@ -183,8 +183,8 @@ class WaveletDetector(TroubledCellDetector):
def _calculate_wavelet_coeffs(self, projection):
output_matrix = []
- for i in range(self.num_coarse_grid_cells):
- new_entry = 0.5*(projection[:, 2*i] @ self.M1 + projection[:, 2*i+1] @ self.M2)
+ for i in range(self._num_coarse_grid_cells):
+ new_entry = 0.5*(projection[:, 2*i] @ self._M1 + projection[:, 2*i+1] @ self._M2)
output_matrix.append(new_entry)
return np.transpose(np.array(output_matrix))
@@ -193,7 +193,7 @@ class WaveletDetector(TroubledCellDetector):
super().plot_results(projection, troubled_cell_history, time_history, color_exact, color_approx)
def _plot_mesh(self, projection, color_exact, color_approx):
- grid, exact = self._calculate_exact_solution(self.mesh[2:-2], self.cell_len)
+ grid, exact = self._calculate_exact_solution(self._mesh[2:-2], self._cell_len)
approx = self._calculate_approximate_solution(projection[:, 1:-1])
pointwise_error = np.abs(exact-approx)
@@ -208,8 +208,8 @@ class WaveletDetector(TroubledCellDetector):
return max_error
def _plot_coarse_mesh(self, projection, color_exact, color_approx):
- 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_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_cell_len)
coarse_projection = self._calculate_coarse_projection(projection)
@@ -220,12 +220,12 @@ class WaveletDetector(TroubledCellDetector):
self._plot_solution_and_approx(grid, exact, approx, color_exact, color_approx)
def _plot_details(self, projection):
- fine_mesh = self.mesh[2:-2]
+ fine_mesh = self._mesh[2:-2]
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)
+ wavelet = AlpertsWavelet(self._polynom_degree).get_vector(xi)
########################################################################################################################
# For later consideration
@@ -233,11 +233,11 @@ class WaveletDetector(TroubledCellDetector):
# tic = timeit.default_timer()
# averaged_projection1 = []
# wavelet_projection1 = []
- # for degree in range(self.polynom_degree + 1):
- # leftMesh = coarse_projection[degree] * self.basis[degree].subs(x, -1 / 2)
- # rightMesh = coarse_projection[degree] * self.basis[degree].subs(x, 1 / 2)
+ # for degree in range(self._polynom_degree + 1):
+ # leftMesh = coarse_projection[degree] * self._basis[degree].subs(x, -1 / 2)
+ # rightMesh = coarse_projection[degree] * self._basis[degree].subs(x, 1 / 2)
# leftTest = multiwavelet_coeffs[degree] * wavelet[degree].subs(xi, 1 / 2) \
- # * (-1)**(self.polynom_degree + 1 + degree)
+ # * (-1)**(self._polynom_degree + 1 + degree)
# rightTest = multiwavelet_coeffs[degree] * wavelet[degree].subs(xi, 1 / 2)
# newRowMesh = []
# newRowTest = []
@@ -253,15 +253,15 @@ class WaveletDetector(TroubledCellDetector):
########################################################################################################################
# tic = timeit.default_timer()
- averaged_projection = [[coarse_projection[degree][cell] * self.basis[degree].subs(x, value)
- for cell in range(self.num_coarse_grid_cells)
+ averaged_projection = [[coarse_projection[degree][cell] * self._basis[degree].subs(x, value)
+ for cell in range(self._num_coarse_grid_cells)
for value in [-0.5, 0.5]]
- for degree in range(self.polynom_degree + 1)]
+ for degree in range(self._polynom_degree + 1)]
wavelet_projection = [[multiwavelet_coeffs[degree][cell] * wavelet[degree].subs(xi, 0.5) * value
- for cell in range(self.num_coarse_grid_cells)
- for value in [(-1) ** (self.polynom_degree + degree + 1), 1]]
- for degree in range(self.polynom_degree + 1)]
+ for cell in range(self._num_coarse_grid_cells)
+ for value in [(-1) ** (self._polynom_degree + degree + 1), 1]]
+ for degree in range(self._polynom_degree + 1)]
# toc = timeit.default_timer()
# print("List:", toc-tic)
@@ -269,8 +269,8 @@ class WaveletDetector(TroubledCellDetector):
# print(wavelet_projection1 == wavelet_projection)
projected_coarse = np.sum(averaged_projection, axis=0)
- projected_fine = np.sum([fine_projection[degree] * self.basis[degree].subs(x, 0)
- for degree in range(self.polynom_degree + 1)], axis=0)
+ projected_fine = np.sum([fine_projection[degree] * self._basis[degree].subs(x, 0)
+ for degree in range(self._polynom_degree + 1)], axis=0)
projected_wavelet_coeffs = np.sum(wavelet_projection, axis=0)
plt.figure(4)
@@ -290,7 +290,7 @@ class WaveletDetector(TroubledCellDetector):
# Calculate projection on coarse mesh
output_matrix = []
- for i in range(self.num_coarse_grid_cells):
+ for i in range(self._num_coarse_grid_cells):
new_entry = 0.5*(projection[:, 2*i] @ basis_matrix_left + projection[:, 2*i+1] @ basis_matrix_right)
output_matrix.append(new_entry)
coarse_projection = np.transpose(np.array(output_matrix))
@@ -299,10 +299,10 @@ class WaveletDetector(TroubledCellDetector):
def _set_basis_matrix(self, first_param, second_param):
matrix = []
- for i in range(self.polynom_degree+1):
+ 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))
+ 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))
row.append(np.float64(entry))
matrix.append(row)
return matrix
@@ -319,15 +319,15 @@ class Boxplot(WaveletDetector):
# Unpack necessary configurations
self.fold_len = config.pop('fold_len', 16)
- self.whisker_len = config.pop('whisker_len', 3)
+ self._whisker_len = config.pop('whisker_len', 3)
def _get_cells(self, multiwavelet_coeffs, projection):
- indexed_coeffs = [[multiwavelet_coeffs[0, i], i]for i in range(self.num_coarse_grid_cells)]
+ indexed_coeffs = [[multiwavelet_coeffs[0, i], i]for i in range(self._num_coarse_grid_cells)]
- if self.num_coarse_grid_cells < self.fold_len:
- self.fold_len = self.num_coarse_grid_cells
+ if self._num_coarse_grid_cells < self.fold_len:
+ self.fold_len = self._num_coarse_grid_cells
- num_folds = self.num_coarse_grid_cells//self.fold_len
+ num_folds = self._num_coarse_grid_cells//self.fold_len
troubled_cells = []
for fold in range(num_folds):
@@ -341,8 +341,8 @@ class Boxplot(WaveletDetector):
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:
@@ -371,14 +371,14 @@ class Theoretical(WaveletDetector):
self.function_name = 'Theoretical'
# Unpack necessary configurations
- 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
def _get_cells(self, multiwavelet_coeffs, projection):
troubled_cells = []
- max_avg = max(1, max(abs(projection[0][degree]) for degree in range(self.polynom_degree+1)))
+ max_avg = max(1, max(abs(projection[0][degree]) for degree in range(self._polynom_degree+1)))
- for cell in range(self.num_coarse_grid_cells):
+ for cell in range(self._num_coarse_grid_cells):
if self._is_troubled_cell(multiwavelet_coeffs, cell, max_avg):
troubled_cells.append(cell)
@@ -386,7 +386,7 @@ class Theoretical(WaveletDetector):
def _is_troubled_cell(self, multiwavelet_coeffs, cell, max_avg):
max_value = max(abs(multiwavelet_coeffs[degree][cell])
- for degree in range(self.polynom_degree+1))/max_avg
- eps = self.cutoff_factor / (self.cell_len*self.num_coarse_grid_cells*2)
+ for degree in range(self._polynom_degree+1))/max_avg
+ eps = self._cutoff_factor / (self._cell_len*self._num_coarse_grid_cells*2)
return max_value > eps
diff --git a/Update_Scheme.py b/Update_Scheme.py
index 5510ed114bb2a1239532f07f6c13bc7beb7cee8d..599550e7e2e45dc384b768a91a2412f8fd2342e4 100644
--- a/Update_Scheme.py
+++ b/Update_Scheme.py
@@ -20,17 +20,17 @@ 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):
# Unpack positional arguments
- self.detector = detector
- self.limiter = limiter
- self.init_cond = init_cond
- self.mesh = mesh
- self.wave_speed = wave_speed
- self.polynom_degree = polynom_degree
- self.num_grid_cells = num_grid_cells
- self.final_time = final_time
- self.history_threshold = history_threshold
- self.left_bound = left_bound
- self.right_bound = right_bound
+ self._detector = detector
+ self._limiter = limiter
+ self._init_cond = init_cond
+ self._mesh = mesh
+ self._wave_speed = wave_speed
+ self._polynom_degree = polynom_degree
+ self._num_grid_cells = num_grid_cells
+ self._final_time = final_time
+ self._history_threshold = history_threshold
+ self._left_bound = left_bound
+ self._right_bound = right_bound
self._reset()
@@ -48,43 +48,43 @@ class UpdateScheme(object):
def _reset(self):
# 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._interval_len = self._right_bound-self._left_bound
+ self._cell_len = self._interval_len / self._num_grid_cells
# Set matrix A
matrix = []
- for i in range(self.polynom_degree+1):
+ for i in range(self._polynom_degree+1):
new_row = []
- for j in range(self.polynom_degree+1):
+ for j in range(self._polynom_degree+1):
new_entry = -1.0
if (j < i) & ((i+j) % 2 == 1):
new_entry = 1.0
new_row.append(new_entry*np.sqrt((i+0.5) * (j+0.5)))
matrix.append(new_row)
- self.A = np.array(matrix) # former: inv_mass @ np.array(matrix)
+ self._A = np.array(matrix) # former: inv_mass @ np.array(matrix)
# Set matrix B
matrix = []
- for i in range(self.polynom_degree+1):
+ for i in range(self._polynom_degree+1):
new_row = []
- for j in range(self.polynom_degree+1):
+ for j in range(self._polynom_degree+1):
new_entry = np.sqrt((i+0.5) * (j+0.5)) * (-1.0)**i
new_row.append(new_entry)
matrix.append(new_row)
- self.B = np.array(matrix) # former: inv_mass @ np.array(matrix)
+ self._B = np.array(matrix) # former: inv_mass @ np.array(matrix)
def _apply_limiter(self, current_projection):
- troubled_cells = self.detector.get_cells(current_projection)
+ troubled_cells = self._detector.get_cells(current_projection)
new_projection = current_projection.copy()
for cell in troubled_cells:
- new_projection[:, cell] = self.limiter.apply(current_projection, cell)
+ new_projection[:, cell] = self._limiter.apply(current_projection, cell)
return new_projection, troubled_cells
def _enforce_boundary_condition(self, current_projection):
- current_projection[:, 0] = current_projection[:, self.num_grid_cells]
- current_projection[:, self.num_grid_cells+1] = current_projection[:, 1]
+ current_projection[:, 0] = current_projection[:, self._num_grid_cells]
+ current_projection[:, self._num_grid_cells+1] = current_projection[:, 1]
return current_projection
@@ -119,12 +119,12 @@ class SSPRK3(UpdateScheme):
# Initialize vector and set first entry to accommodate for ghost cell
right_hand_side = [0]
- for j in range(self.num_grid_cells):
- right_hand_side.append(2*(self.A @ current_projection[:, j+1]
- + self.B @ current_projection[:, j]))
+ for j in range(self._num_grid_cells):
+ right_hand_side.append(2*(self._A @ current_projection[:, j+1]
+ + self._B @ current_projection[:, j]))
# Set ghost cells to respective value
- right_hand_side[0] = right_hand_side[self.num_grid_cells]
+ right_hand_side[0] = right_hand_side[self._num_grid_cells]
right_hand_side.append(right_hand_side[1])
return np.transpose(right_hand_side)
diff --git a/Vectors_of_Polynomials.py b/Vectors_of_Polynomials.py
index 07b597fbe3d53066a360fa9c4ef7037407313185..60bc2a7758129a5f4528a5314c9714e63a41af8a 100644
--- a/Vectors_of_Polynomials.py
+++ b/Vectors_of_Polynomials.py
@@ -8,7 +8,7 @@ import numpy as np
class Vector(object):
def __init__(self, polynom_degree):
- self.polynom_degree = polynom_degree
+ self._polynom_degree = polynom_degree
def get_vector(self, eval_points):
pass
@@ -20,7 +20,7 @@ class Legendre(Vector):
def _calculate_legendre_vector(self, eval_point):
vector = []
- for degree in range(self.polynom_degree+1):
+ for degree in range(self._polynom_degree+1):
if degree == 0:
vector.append(1.0 + 0*eval_point)
else:
@@ -35,12 +35,12 @@ class Legendre(Vector):
class OrthonormalLegendre(Legendre):
def get_vector(self, eval_point):
leg_vector = self._calculate_legendre_vector(eval_point)
- return [leg_vector[degree] * np.sqrt(degree+0.5) for degree in range(self.polynom_degree+1)]
+ return [leg_vector[degree] * np.sqrt(degree+0.5) for degree in range(self._polynom_degree+1)]
class AlpertsWavelet(Vector):
def get_vector(self, eval_point):
- degree = self.polynom_degree
+ degree = self._polynom_degree
if degree == 0:
return [np.sqrt(0.5) + eval_point*0]