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

Set all instance variables to protected.

parent c0bb51db
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
DG_Approximation.py
+ 79
− 78
View file @ f897dc57
...@@ -12,20 +12,21 @@ TODO: Combine initial projection and approx solution somehow ...@@ -12,20 +12,21 @@ TODO: Combine initial projection and approx solution somehow
TODO: Investigate why there are no weights in approx calc -> Ask TODO: Investigate why there are no weights in approx calc -> Ask
TODO: Change order of methods TODO: Change order of methods
TODO: Change code to not include self. but global variables; NO, do private \ 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: Find better names for A, B, anything pertaining Gauss
TODO: Write documentation for all methods 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: Contemplate moving A and B to Vectors_of_Polynomials
TODO: Improve access to names of classes TODO: Improve access to names of classes
TODO: Check plots for NoDetection class -> Done (Adjusted scale of x-axis for shock tube) 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: 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: Improve saving of plots
TODO: Extend color options TODO: Extend color options
TODO: Implement type check for all kwargs and configs TODO: Implement type check for all kwargs and configs
TODO: Make sure that self.cfl_number is not changed -> Done 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 import numpy as np
...@@ -50,24 +51,24 @@ xi = Symbol('z') ...@@ -50,24 +51,24 @@ xi = Symbol('z')
class DGScheme(object): class DGScheme(object):
def __init__(self, detector, **kwargs): def __init__(self, detector, **kwargs):
# Unpack keyword arguments # Unpack keyword arguments
self.wave_speed = kwargs.pop('wave_speed', 1) self._wave_speed = kwargs.pop('wave_speed', 1)
self.polynom_degree = kwargs.pop('polynom_degree', 2) self._polynom_degree = kwargs.pop('polynom_degree', 2)
self.cfl_number = kwargs.pop('cfl_number', 0.2) self._cfl_number = kwargs.pop('cfl_number', 0.2)
self.num_grid_cells = kwargs.pop('num_grid_cells', 64) self._num_grid_cells = kwargs.pop('num_grid_cells', 64)
self.final_time = kwargs.pop('final_time', 1) self._final_time = kwargs.pop('final_time', 1)
self.left_bound = kwargs.pop('left_bound', -1) self._left_bound = kwargs.pop('left_bound', -1)
self.right_bound = kwargs.pop('right_bound', 1) self._right_bound = kwargs.pop('right_bound', 1)
self.verbose = kwargs.pop('verbose', False) self._verbose = kwargs.pop('verbose', False)
self.history_threshold = kwargs.pop('history_threshold', math.ceil(0.2/self.cfl_number)) self._history_threshold = kwargs.pop('history_threshold', math.ceil(0.2/self._cfl_number))
self.detector = detector self._detector = detector
self.detector_config = kwargs.pop('detector_config', {}) self._detector_config = kwargs.pop('detector_config', {})
self.init_cond = kwargs.pop('init_cond', 'Sine') self._init_cond = kwargs.pop('init_cond', 'Sine')
self.init_config = kwargs.pop('init_config', {}) self._init_config = kwargs.pop('init_config', {})
self.limiter = kwargs.pop('limiter', 'ModifiedMinMod') self._limiter = kwargs.pop('limiter', 'ModifiedMinMod')
self.limiter_config = kwargs.pop('limiter_config', {}) self._limiter_config = kwargs.pop('limiter_config', {})
self.quadrature = kwargs.pop('quadrature', 'Gauss') self._quadrature = kwargs.pop('quadrature', 'Gauss')
self.quadrature_config = kwargs.pop('quadrature_config', {}) self._quadrature_config = kwargs.pop('quadrature_config', {})
self.update_scheme = kwargs.pop('update_scheme', 'SSPRK3') self._update_scheme = kwargs.pop('update_scheme', 'SSPRK3')
# Throw an error if there are extra keyword arguments # Throw an error if there are extra keyword arguments
if len(kwargs) > 0: if len(kwargs) > 0:
...@@ -75,35 +76,35 @@ class DGScheme(object): ...@@ -75,35 +76,35 @@ class DGScheme(object):
raise ValueError('Unrecognized arguments: %s' % extra) raise ValueError('Unrecognized arguments: %s' % extra)
# Make sure all classes actually exist # Make sure all classes actually exist
if not hasattr(Troubled_Cell_Detector, self.detector): 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): 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): 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): 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): 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)
self._reset() self._reset()
# Replace the string names with the actual class instances # Replace the string names with the actual class instances
# (and add the instance variables for the quadrature) # (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._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._limiter = getattr(Limiter, self._limiter)(self._limiter_config)
self.quadrature = getattr(Quadrature, self.quadrature)(self.quadrature_config) self._quadrature = getattr(Quadrature, self._quadrature)(self._quadrature_config)
self.detector = getattr(Troubled_Cell_Detector, self.detector)( self._detector = getattr(Troubled_Cell_Detector, self._detector)(
self.detector_config, self.mesh, self.wave_speed, self.polynom_degree, self.num_grid_cells, 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._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._update_scheme = getattr(Update_Scheme, self._update_scheme)(
self.detector, self.limiter, self.init_cond, self.mesh, self.wave_speed, self.polynom_degree, 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) self._num_grid_cells, self._final_time, self._history_threshold, self._left_bound, self._right_bound)
# print(self.detector.get_name()) # print(self._detector.get_name())
# print(self.detector.__class__.__name__) # print(self._detector.__class__.__name__)
# print(Troubled_Cell_Detector.Boxplot.__bases__) # 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__) # print(Troubled_Cell_Detector.WaveletDetector in Troubled_Cell_Detector.Boxplot.__mro__)
def approximate(self): def approximate(self):
...@@ -115,41 +116,41 @@ class DGScheme(object): ...@@ -115,41 +116,41 @@ class DGScheme(object):
projection = self._do_initial_projection() 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 current_time = 0
iteration = 0 iteration = 0
troubled_cell_history = [] troubled_cell_history = []
time_history = [] time_history = []
while current_time < self.final_time: while current_time < self._final_time:
# Adjust for last cell # Adjust for last cell
cfl_number = self.cfl_number cfl_number = self._cfl_number
if current_time+time_step > self.final_time: if current_time+time_step > self._final_time:
time_step = self.final_time-current_time time_step = self._final_time-current_time
cfl_number = self.wave_speed * time_step / self.num_grid_cells cfl_number = self._wave_speed * time_step / self._num_grid_cells
# Update projection # 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 iteration += 1
if (iteration % self.history_threshold) == 0: if (iteration % self._history_threshold) == 0:
troubled_cell_history.append(troubled_cells) troubled_cell_history.append(troubled_cells)
time_history.append(current_time) time_history.append(current_time)
current_time += time_step current_time += time_step
# Plot exact/approximate results, errors, shock tubes and any detector-dependant plots # 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() plt.show()
def save_plots(self): def save_plots(self):
name = self.init_cond.get_name() + '__' + self.detector.get_name() + '__' + self.limiter.get_name() \ 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_' \ + '__' + 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._final_time) + '__wave_speed_' + str(self._wave_speed) + '__number_of_cells_' \
+ str(self.num_grid_cells) + '__polynomial_degree_' + str(self.polynom_degree) + str(self._num_grid_cells) + '__polynomial_degree_' + str(self._polynom_degree)
save_fig_1 = True save_fig_1 = True
if save_fig_1: if save_fig_1:
...@@ -198,50 +199,50 @@ class DGScheme(object): ...@@ -198,50 +199,50 @@ class DGScheme(object):
os.makedirs('testfig/shock_tube') os.makedirs('testfig/shock_tube')
# Set additional necessary instance variables # Set additional necessary instance variables
self.interval_len = self.right_bound-self.left_bound self._interval_len = self._right_bound-self._left_bound
self.cell_len = self.interval_len / self.num_grid_cells 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 # 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 # 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 self._cell_len) # +3/2
# Set inverse mass matrix # Set inverse mass matrix
mass_matrix = [] mass_matrix = []
for i in range(self.polynom_degree+1): for i in range(self._polynom_degree+1):
new_row = [] new_row = []
for j in range(self.polynom_degree+1): for j in range(self._polynom_degree+1):
new_entry = 0.0 new_entry = 0.0
if i == j: if i == j:
new_entry = 1.0 new_entry = 1.0
new_row.append(new_entry) new_row.append(new_entry)
mass_matrix.append(new_row) mass_matrix.append(new_row)
self.inv_mass = np.array(mass_matrix) self._inv_mass = np.array(mass_matrix)
def _do_initial_projection(self): def _do_initial_projection(self):
# Initialize matrix and set first entry to accommodate for ghost cell # Initialize matrix and set first entry to accommodate for ghost cell
output_matrix = [0] output_matrix = [0]
for cell in range(self.num_grid_cells): for cell in range(self._num_grid_cells):
new_row = [] new_row = []
eval_point = self.left_bound + (cell+0.5)*self.cell_len eval_point = self._left_bound + (cell+0.5)*self._cell_len
for degree in range(self.polynom_degree + 1): for degree in range(self._polynom_degree + 1):
new_entry = sum(self.init_cond.calculate(eval_point new_entry = sum(self._init_cond.calculate(eval_point
+ self.cell_len/2 * self.quadrature.get_eval_points()[point]) + self._cell_len/2 * self._quadrature.get_eval_points()[point])
* self.basis[degree].subs(x, self.quadrature.get_eval_points()[point]) * self._basis[degree].subs(x, self._quadrature.get_eval_points()[point])
* self.quadrature.get_weights()[point] * self._quadrature.get_weights()[point]
for point in range(self.quadrature.get_num_points())) for point in range(self._quadrature.get_num_points()))
new_row.append(np.float64(new_entry)) new_row.append(np.float64(new_entry))
new_row = np.array(new_row) 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 # 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]) output_matrix.append(output_matrix[1])
print(np.array(output_matrix).shape) print(np.array(output_matrix).shape)
......
Initial_Condition.py
+ 26
− 26
View file @ f897dc57
...@@ -8,9 +8,9 @@ import numpy as np ...@@ -8,9 +8,9 @@ import numpy as np
class InitialCondition(object): class InitialCondition(object):
def __init__(self, left_bound, right_bound, config): def __init__(self, left_bound, right_bound, config):
self.left_bound = left_bound self._left_bound = left_bound
self.right_bound = right_bound self._right_bound = right_bound
self.interval_len = self.right_bound-self.left_bound self._interval_len = self._right_bound-self._left_bound
self.function_name = 'None' self.function_name = 'None'
...@@ -18,10 +18,10 @@ class InitialCondition(object): ...@@ -18,10 +18,10 @@ class InitialCondition(object):
return self.function_name return self.function_name
def calculate(self, x): def calculate(self, x):
while x < self.left_bound: while x < self._left_bound:
x = x + self.interval_len x = x + self._interval_len
while x > self.right_bound: while x > self._right_bound:
x = x - self.interval_len x = x - self._interval_len
return self._get_point(x) return self._get_point(x)
def _get_point(self, x): def _get_point(self, x):
...@@ -36,10 +36,10 @@ class Sine(InitialCondition): ...@@ -36,10 +36,10 @@ class Sine(InitialCondition):
self.function_name = 'Sine' self.function_name = 'Sine'
# Unpack necessary configurations # Unpack necessary configurations
self.factor = config.pop('factor', 2) self._factor = config.pop('factor', 2)
def _get_point(self, x): def _get_point(self, x):
return np.sin(self.factor * np.pi * x) return np.sin(self._factor * np.pi * x)
class Box(InitialCondition): class Box(InitialCondition):
...@@ -68,28 +68,28 @@ class FourPeakWave(InitialCondition): ...@@ -68,28 +68,28 @@ class FourPeakWave(InitialCondition):
self.function_name = 'FourPeakWave' self.function_name = 'FourPeakWave'
# Set additional necessary parameter # Set additional necessary parameter
self.alpha = 10 self._alpha = 10
self.delta = 0.005 self._delta = 0.005
self.beta = np.log(2) / (36 * self.delta**2) self._beta = np.log(2) / (36 * self._delta**2)
self.a = 0.5 self._a = 0.5
self.z = -0.7 self._z = -0.7
def _get_point(self, x): def _get_point(self, x):
if (x >= -0.8) & (x <= -0.6): 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): if (x >= -0.4) & (x <= -0.2):
return 1 return 1
if (x >= 0) & (x <= 0.2): if (x >= 0) & (x <= 0.2):
return 1 - abs(10 * (x-0.1)) return 1 - abs(10 * (x-0.1))
if (x >= 0.4) & (x <= 0.6): 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 return 0
def _G(self, x, z): 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): 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): class Linear(InitialCondition):
...@@ -100,10 +100,10 @@ class Linear(InitialCondition): ...@@ -100,10 +100,10 @@ class Linear(InitialCondition):
self.function_name = 'Linear' self.function_name = 'Linear'
# Unpack necessary configurations # Unpack necessary configurations
self.factor = config.pop('factor', 1) self._factor = config.pop('factor', 1)
def _get_point(self, x): def _get_point(self, x):
return self.factor * x return self._factor * x
class LinearAbsolut(InitialCondition): class LinearAbsolut(InitialCondition):
...@@ -114,10 +114,10 @@ class LinearAbsolut(InitialCondition): ...@@ -114,10 +114,10 @@ class LinearAbsolut(InitialCondition):
self.function_name = 'LinearAbsolut' self.function_name = 'LinearAbsolut'
# Unpack necessary configurations # Unpack necessary configurations
self.factor = config.pop('factor', 1) self._factor = config.pop('factor', 1)
def _get_point(self, x): def _get_point(self, x):
return self.factor * abs(x) return self._factor * abs(x)
class DiscontinuousConstant(InitialCondition): class DiscontinuousConstant(InitialCondition):
...@@ -128,9 +128,9 @@ class DiscontinuousConstant(InitialCondition): ...@@ -128,9 +128,9 @@ class DiscontinuousConstant(InitialCondition):
self.function_name = 'DiscontinuousConstant' self.function_name = 'DiscontinuousConstant'
# Unpack necessary configurations # Unpack necessary configurations
self.x0 = config.pop('x0', 0) self._x0 = config.pop('x0', 0)
self.left_factor = config.pop('left_factor', 1) self._left_factor = config.pop('left_factor', 1)
self.right_factor = config.pop('right_factor', 0.5) self._right_factor = config.pop('right_factor', 0.5)
def _get_point(self, x): 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)
Limiter.py
+ 8
− 8
View file @ f897dc57
...@@ -32,10 +32,10 @@ class MinMod(Limiter): ...@@ -32,10 +32,10 @@ class MinMod(Limiter):
super().__init__(config) super().__init__(config)
# Unpack necessary configurations # Unpack necessary configurations
self.erase_degree = config.pop('erase_degree', 0) self._erase_degree = config.pop('erase_degree', 0)
# Set name of function # 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): def apply(self, projection, cell):
cell_slope = self._set_cell_slope(projection, cell) cell_slope = self._set_cell_slope(projection, cell)
...@@ -46,7 +46,7 @@ class MinMod(Limiter): ...@@ -46,7 +46,7 @@ class MinMod(Limiter):
adapted_projection = projection[:, cell].copy() adapted_projection = projection[:, cell].copy()
for i in range(len(adapted_projection)): for i in range(len(adapted_projection)):
if i > self.erase_degree: if i > self._erase_degree:
adapted_projection[i] = 0 adapted_projection[i] = 0
return adapted_projection return adapted_projection
...@@ -72,16 +72,16 @@ class ModifiedMinMod(MinMod): ...@@ -72,16 +72,16 @@ class ModifiedMinMod(MinMod):
super().__init__(config) super().__init__(config)
# Unpack necessary configurations # Unpack necessary configurations
self.cell_len = config.pop('cell_len') self._cell_len = config.pop('cell_len')
self.mod_factor = config.pop('mod_factor', 0) self._mod_factor = config.pop('mod_factor', 0)
# Set name of function # 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): 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 True
return super()._determine_modification(projection, cell, cell_slope) return super()._determine_modification(projection, cell, cell_slope)
Quadrature.py
+ 9
− 9
View file @ f897dc57
...@@ -9,21 +9,21 @@ import numpy.polynomial.legendre as leg ...@@ -9,21 +9,21 @@ import numpy.polynomial.legendre as leg
class Quadrature(object): class Quadrature(object):
def __init__(self, config): def __init__(self, config):
self.function_name = 'None' self.function_name = 'None'
self.eval_points = None self._eval_points = None
self.weights = None self._weights = None
self.num_eval_points = None self._num_eval_points = None
def get_name(self): def get_name(self):
return self.function_name return self.function_name
def get_eval_points(self): def get_eval_points(self):
return self.eval_points return self._eval_points
def get_weights(self): def get_weights(self):
return self.weights return self._weights
def get_num_points(self): def get_num_points(self):
return self.num_eval_points return self._num_eval_points
class Gauss(Quadrature): class Gauss(Quadrature):
...@@ -31,9 +31,9 @@ class Gauss(Quadrature): ...@@ -31,9 +31,9 @@ class Gauss(Quadrature):
super().__init__(config) super().__init__(config)
# Unpack necessary configurations # 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 # Set name of function
self.function_name = 'Gauss' + str(self.num_eval_points) self.function_name = 'Gauss' + str(self._num_eval_points)
Troubled_Cell_Detector.py
+ 69
− 69
View file @ f897dc57
...@@ -19,18 +19,18 @@ xi = Symbol('z') ...@@ -19,18 +19,18 @@ xi = Symbol('z')
class TroubledCellDetector(object): class TroubledCellDetector(object):
def __init__(self, config, mesh, wave_speed, polynom_degree, num_grid_cells, final_time, left_bound, right_bound, def __init__(self, config, mesh, wave_speed, polynom_degree, num_grid_cells, final_time, left_bound, right_bound,
basis, init_cond, quadrature): basis, init_cond, quadrature):
self.mesh = mesh self._mesh = mesh
self.wave_speed = wave_speed self._wave_speed = wave_speed
self.polynom_degree = polynom_degree self._polynom_degree = polynom_degree
self.num_grid_cells = num_grid_cells self._num_grid_cells = num_grid_cells
self.final_time = final_time self._final_time = final_time
self.left_bound = left_bound self._left_bound = left_bound
self.right_bound = right_bound self._right_bound = right_bound
self.interval_len = right_bound - left_bound self._interval_len = right_bound - left_bound
self.cell_len = self.interval_len / num_grid_cells self._cell_len = self._interval_len / num_grid_cells
self.basis = basis self._basis = basis
self.init_cond = init_cond self._init_cond = init_cond
self.quadrature = quadrature self._quadrature = quadrature
self.function_name = 'None' self.function_name = 'None'
def get_name(self): def get_name(self):
...@@ -43,12 +43,12 @@ class TroubledCellDetector(object): ...@@ -43,12 +43,12 @@ class TroubledCellDetector(object):
self._plot_shock_tube(troubled_cell_history, time_history) self._plot_shock_tube(troubled_cell_history, time_history)
max_error = self._plot_mesh(projection, color_exact, color_approx) max_error = self._plot_mesh(projection, color_exact, color_approx)
print("p =", self.polynom_degree) print("p =", self._polynom_degree)
print("N =", self.num_grid_cells) print("N =", self._num_grid_cells)
print("maximum error =", max_error) print("maximum error =", max_error)
def _plot_mesh(self, projection, 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]) approx = self._calculate_approximate_solution(projection[:, 1:-1])
pointwise_error = np.abs(exact-approx) pointwise_error = np.abs(exact-approx)
...@@ -90,7 +90,7 @@ class TroubledCellDetector(object): ...@@ -90,7 +90,7 @@ class TroubledCellDetector(object):
current_cells = troubled_cell_history[pos] current_cells = troubled_cell_history[pos]
for cell in current_cells: for cell in current_cells:
plt.plot(cell, time_history[pos], 'k.') 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.xlabel('Cell')
plt.ylabel('Time') plt.ylabel('Time')
plt.title('Shock Tubes') plt.title('Shock Tubes')
...@@ -98,15 +98,15 @@ class TroubledCellDetector(object): ...@@ -98,15 +98,15 @@ class TroubledCellDetector(object):
def _calculate_exact_solution(self, mesh, cell_len): def _calculate_exact_solution(self, mesh, cell_len):
grid = [] grid = []
exact = [] 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)): 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 = [] eval_values = []
for point in range(len(eval_points)): 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) + num_periods*self._interval_len)
eval_values.append(new_entry) eval_values.append(new_entry)
grid.append(eval_points) grid.append(eval_points)
...@@ -118,15 +118,15 @@ class TroubledCellDetector(object): ...@@ -118,15 +118,15 @@ class TroubledCellDetector(object):
return grid, exact return grid, exact
def _calculate_approximate_solution(self, projection): def _calculate_approximate_solution(self, projection):
points = self.quadrature.get_eval_points() points = self._quadrature.get_eval_points()
num_points = self.quadrature.get_num_points() num_points = self._quadrature.get_num_points()
basis = self.basis 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)] for degree in range(self._polynom_degree+1)]
approx = [[sum(projection[degree][cell] * basis_matrix[degree][point] 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 point in range(num_points)]
for cell in range(len(projection[0]))] for cell in range(len(projection[0]))]
...@@ -153,23 +153,23 @@ class WaveletDetector(TroubledCellDetector): ...@@ -153,23 +153,23 @@ class WaveletDetector(TroubledCellDetector):
basis, init_cond, quadrature) basis, init_cond, quadrature)
# Set fixed basis and wavelet vectors # 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) self._wavelet = AlpertsWavelet(self._polynom_degree).get_vector(x)
# Set additional necessary parameter # Set additional necessary parameter
self.num_coarse_grid_cells = self.num_grid_cells//2 self._num_coarse_grid_cells = self._num_grid_cells//2
self.M1 = self._set_multiwavelet_matrix(xi, -0.5*(xi-1), True) 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._M2 = self._set_multiwavelet_matrix(xi, 0.5*(xi+1), False)
def _set_multiwavelet_matrix(self, first_param, second_param, is_M1): def _set_multiwavelet_matrix(self, first_param, second_param, is_M1):
matrix = [] matrix = []
for i in range(self.polynom_degree+1): for i in range(self._polynom_degree+1):
row = [] row = []
for j in range(self.polynom_degree+1): 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)) (xi, -1, 1))
if is_M1: if is_M1:
entry = entry * (-1)**(j + self.polynom_degree + 1) entry = entry * (-1)**(j + self._polynom_degree + 1)
row.append(np.float64(entry)) row.append(np.float64(entry))
matrix.append(row) matrix.append(row)
return matrix return matrix
...@@ -183,8 +183,8 @@ class WaveletDetector(TroubledCellDetector): ...@@ -183,8 +183,8 @@ class WaveletDetector(TroubledCellDetector):
def _calculate_wavelet_coeffs(self, projection): def _calculate_wavelet_coeffs(self, projection):
output_matrix = [] 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] @ self.M1 + projection[:, 2*i+1] @ self.M2) new_entry = 0.5*(projection[:, 2*i] @ self._M1 + projection[:, 2*i+1] @ self._M2)
output_matrix.append(new_entry) output_matrix.append(new_entry)
return np.transpose(np.array(output_matrix)) return np.transpose(np.array(output_matrix))
...@@ -193,7 +193,7 @@ class WaveletDetector(TroubledCellDetector): ...@@ -193,7 +193,7 @@ class WaveletDetector(TroubledCellDetector):
super().plot_results(projection, troubled_cell_history, time_history, color_exact, color_approx) super().plot_results(projection, troubled_cell_history, time_history, color_exact, color_approx)
def _plot_mesh(self, projection, 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]) approx = self._calculate_approximate_solution(projection[:, 1:-1])
pointwise_error = np.abs(exact-approx) pointwise_error = np.abs(exact-approx)
...@@ -208,8 +208,8 @@ class WaveletDetector(TroubledCellDetector): ...@@ -208,8 +208,8 @@ class WaveletDetector(TroubledCellDetector):
return max_error return max_error
def _plot_coarse_mesh(self, projection, color_exact, color_approx): def _plot_coarse_mesh(self, projection, color_exact, color_approx):
coarse_cell_len = 2*self.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_mesh = np.arange(self._left_bound - (0.5*coarse_cell_len), self._right_bound + (1.5*coarse_cell_len),
coarse_cell_len) coarse_cell_len)
coarse_projection = self._calculate_coarse_projection(projection) coarse_projection = self._calculate_coarse_projection(projection)
...@@ -220,12 +220,12 @@ class WaveletDetector(TroubledCellDetector): ...@@ -220,12 +220,12 @@ class WaveletDetector(TroubledCellDetector):
self._plot_solution_and_approx(grid, exact, approx, color_exact, color_approx) self._plot_solution_and_approx(grid, exact, approx, color_exact, color_approx)
def _plot_details(self, projection): def _plot_details(self, projection):
fine_mesh = self.mesh[2:-2] fine_mesh = self._mesh[2:-2]
fine_projection = projection[:, 1:-1] fine_projection = projection[:, 1:-1]
coarse_projection = self._calculate_coarse_projection(projection) coarse_projection = self._calculate_coarse_projection(projection)
multiwavelet_coeffs = self._calculate_wavelet_coeffs(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 # For later consideration
...@@ -233,11 +233,11 @@ class WaveletDetector(TroubledCellDetector): ...@@ -233,11 +233,11 @@ class WaveletDetector(TroubledCellDetector):
# tic = timeit.default_timer() # tic = timeit.default_timer()
# averaged_projection1 = [] # averaged_projection1 = []
# wavelet_projection1 = [] # wavelet_projection1 = []
# for degree in range(self.polynom_degree + 1): # for degree in range(self._polynom_degree + 1):
# leftMesh = coarse_projection[degree] * self.basis[degree].subs(x, -1 / 2) # leftMesh = coarse_projection[degree] * self._basis[degree].subs(x, -1 / 2)
# rightMesh = 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) \ # 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) # rightTest = multiwavelet_coeffs[degree] * wavelet[degree].subs(xi, 1 / 2)
# newRowMesh = [] # newRowMesh = []
# newRowTest = [] # newRowTest = []
...@@ -253,15 +253,15 @@ class WaveletDetector(TroubledCellDetector): ...@@ -253,15 +253,15 @@ class WaveletDetector(TroubledCellDetector):
######################################################################################################################## ########################################################################################################################
# tic = timeit.default_timer() # tic = timeit.default_timer()
averaged_projection = [[coarse_projection[degree][cell] * self.basis[degree].subs(x, value) averaged_projection = [[coarse_projection[degree][cell] * self._basis[degree].subs(x, value)
for cell in range(self.num_coarse_grid_cells) for cell in range(self._num_coarse_grid_cells)
for value in [-0.5, 0.5]] 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 wavelet_projection = [[multiwavelet_coeffs[degree][cell] * wavelet[degree].subs(xi, 0.5) * value
for cell in range(self.num_coarse_grid_cells) for cell in range(self._num_coarse_grid_cells)
for value in [(-1) ** (self.polynom_degree + degree + 1), 1]] for value in [(-1) ** (self._polynom_degree + degree + 1), 1]]
for degree in range(self.polynom_degree + 1)] for degree in range(self._polynom_degree + 1)]
# toc = timeit.default_timer() # toc = timeit.default_timer()
# print("List:", toc-tic) # print("List:", toc-tic)
...@@ -269,8 +269,8 @@ class WaveletDetector(TroubledCellDetector): ...@@ -269,8 +269,8 @@ class WaveletDetector(TroubledCellDetector):
# print(wavelet_projection1 == wavelet_projection) # print(wavelet_projection1 == wavelet_projection)
projected_coarse = np.sum(averaged_projection, axis=0) projected_coarse = np.sum(averaged_projection, axis=0)
projected_fine = np.sum([fine_projection[degree] * self.basis[degree].subs(x, 0) projected_fine = np.sum([fine_projection[degree] * self._basis[degree].subs(x, 0)
for degree in range(self.polynom_degree + 1)], axis=0) for degree in range(self._polynom_degree + 1)], axis=0)
projected_wavelet_coeffs = np.sum(wavelet_projection, axis=0) projected_wavelet_coeffs = np.sum(wavelet_projection, axis=0)
plt.figure(4) plt.figure(4)
...@@ -290,7 +290,7 @@ class WaveletDetector(TroubledCellDetector): ...@@ -290,7 +290,7 @@ class WaveletDetector(TroubledCellDetector):
# Calculate projection on coarse mesh # Calculate projection on coarse mesh
output_matrix = [] 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) new_entry = 0.5*(projection[:, 2*i] @ basis_matrix_left + projection[:, 2*i+1] @ basis_matrix_right)
output_matrix.append(new_entry) output_matrix.append(new_entry)
coarse_projection = np.transpose(np.array(output_matrix)) coarse_projection = np.transpose(np.array(output_matrix))
...@@ -299,10 +299,10 @@ class WaveletDetector(TroubledCellDetector): ...@@ -299,10 +299,10 @@ class WaveletDetector(TroubledCellDetector):
def _set_basis_matrix(self, first_param, second_param): def _set_basis_matrix(self, first_param, second_param):
matrix = [] matrix = []
for i in range(self.polynom_degree+1): for i in range(self._polynom_degree+1):
row = [] row = []
for j in range(self.polynom_degree+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)) entry = integrate(self._basis[i].subs(x, first_param) * self._basis[j].subs(x, second_param), (xi, -1, 1))
row.append(np.float64(entry)) row.append(np.float64(entry))
matrix.append(row) matrix.append(row)
return matrix return matrix
...@@ -319,15 +319,15 @@ class Boxplot(WaveletDetector): ...@@ -319,15 +319,15 @@ class Boxplot(WaveletDetector):
# Unpack necessary configurations # Unpack necessary configurations
self.fold_len = config.pop('fold_len', 16) 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): 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: if self._num_coarse_grid_cells < self.fold_len:
self.fold_len = self.num_coarse_grid_cells 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 = [] troubled_cells = []
for fold in range(num_folds): for fold in range(num_folds):
...@@ -341,8 +341,8 @@ class Boxplot(WaveletDetector): ...@@ -341,8 +341,8 @@ class Boxplot(WaveletDetector):
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] + balance_factor * sorted_fold[3*boundary_index][0]
lower_bound = first_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) upper_bound = third_quartile + self._whisker_len * (third_quartile-first_quartile)
# Check for lower extreme outliers and add respective cells # Check for lower extreme outliers and add respective cells
for cell in sorted_fold: for cell in sorted_fold:
...@@ -371,14 +371,14 @@ class Theoretical(WaveletDetector): ...@@ -371,14 +371,14 @@ class Theoretical(WaveletDetector):
self.function_name = 'Theoretical' self.function_name = 'Theoretical'
# Unpack necessary configurations # 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 # comment to line above: or 2 or 3
def _get_cells(self, multiwavelet_coeffs, projection): def _get_cells(self, multiwavelet_coeffs, projection):
troubled_cells = [] 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): if self._is_troubled_cell(multiwavelet_coeffs, cell, max_avg):
troubled_cells.append(cell) troubled_cells.append(cell)
...@@ -386,7 +386,7 @@ class Theoretical(WaveletDetector): ...@@ -386,7 +386,7 @@ class Theoretical(WaveletDetector):
def _is_troubled_cell(self, multiwavelet_coeffs, cell, max_avg): def _is_troubled_cell(self, multiwavelet_coeffs, cell, max_avg):
max_value = max(abs(multiwavelet_coeffs[degree][cell]) max_value = max(abs(multiwavelet_coeffs[degree][cell])
for degree in range(self.polynom_degree+1))/max_avg for degree in range(self._polynom_degree+1))/max_avg
eps = self.cutoff_factor / (self.cell_len*self.num_coarse_grid_cells*2) eps = self._cutoff_factor / (self._cell_len*self._num_coarse_grid_cells*2)
return max_value > eps return max_value > eps
Update_Scheme.py
+ 27
− 27
View file @ f897dc57
...@@ -20,17 +20,17 @@ class UpdateScheme(object): ...@@ -20,17 +20,17 @@ class UpdateScheme(object):
def __init__(self, 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): history_threshold, left_bound, right_bound):
# Unpack positional arguments # Unpack positional arguments
self.detector = detector self._detector = detector
self.limiter = limiter self._limiter = limiter
self.init_cond = init_cond self._init_cond = init_cond
self.mesh = mesh self._mesh = mesh
self.wave_speed = wave_speed self._wave_speed = wave_speed
self.polynom_degree = polynom_degree self._polynom_degree = polynom_degree
self.num_grid_cells = num_grid_cells self._num_grid_cells = num_grid_cells
self.final_time = final_time self._final_time = final_time
self.history_threshold = history_threshold self._history_threshold = history_threshold
self.left_bound = left_bound self._left_bound = left_bound
self.right_bound = right_bound self._right_bound = right_bound
self._reset() self._reset()
...@@ -48,43 +48,43 @@ class UpdateScheme(object): ...@@ -48,43 +48,43 @@ class UpdateScheme(object):
def _reset(self): def _reset(self):
# Set additional necessary fixed instance variables # Set additional necessary fixed instance variables
self.name = 'None' self.name = 'None'
self.interval_len = self.right_bound-self.left_bound self._interval_len = self._right_bound-self._left_bound
self.cell_len = self.interval_len / self.num_grid_cells self._cell_len = self._interval_len / self._num_grid_cells
# Set matrix A # Set matrix A
matrix = [] matrix = []
for i in range(self.polynom_degree+1): for i in range(self._polynom_degree+1):
new_row = [] new_row = []
for j in range(self.polynom_degree+1): for j in range(self._polynom_degree+1):
new_entry = -1.0 new_entry = -1.0
if (j < i) & ((i+j) % 2 == 1): if (j < i) & ((i+j) % 2 == 1):
new_entry = 1.0 new_entry = 1.0
new_row.append(new_entry*np.sqrt((i+0.5) * (j+0.5))) new_row.append(new_entry*np.sqrt((i+0.5) * (j+0.5)))
matrix.append(new_row) 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 # Set matrix B
matrix = [] matrix = []
for i in range(self.polynom_degree+1): for i in range(self._polynom_degree+1):
new_row = [] 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_entry = np.sqrt((i+0.5) * (j+0.5)) * (-1.0)**i
new_row.append(new_entry) new_row.append(new_entry)
matrix.append(new_row) 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): 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() new_projection = current_projection.copy()
for cell in troubled_cells: 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 return new_projection, troubled_cells
def _enforce_boundary_condition(self, current_projection): def _enforce_boundary_condition(self, current_projection):
current_projection[:, 0] = current_projection[:, self.num_grid_cells] current_projection[:, 0] = current_projection[:, self._num_grid_cells]
current_projection[:, self.num_grid_cells+1] = current_projection[:, 1] current_projection[:, self._num_grid_cells+1] = current_projection[:, 1]
return current_projection return current_projection
...@@ -119,12 +119,12 @@ class SSPRK3(UpdateScheme): ...@@ -119,12 +119,12 @@ class SSPRK3(UpdateScheme):
# Initialize vector and set first entry to accommodate for ghost cell # Initialize vector and set first entry to accommodate for ghost cell
right_hand_side = [0] right_hand_side = [0]
for j in range(self.num_grid_cells): for j in range(self._num_grid_cells):
right_hand_side.append(2*(self.A @ current_projection[:, j+1] right_hand_side.append(2*(self._A @ current_projection[:, j+1]
+ self.B @ current_projection[:, j])) + self._B @ current_projection[:, j]))
# Set ghost cells to respective value # 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]) right_hand_side.append(right_hand_side[1])
return np.transpose(right_hand_side) return np.transpose(right_hand_side)
......
Vectors_of_Polynomials.py
+ 4
− 4
View file @ f897dc57
...@@ -8,7 +8,7 @@ import numpy as np ...@@ -8,7 +8,7 @@ import numpy as np
class Vector(object): class Vector(object):
def __init__(self, polynom_degree): def __init__(self, polynom_degree):
self.polynom_degree = polynom_degree self._polynom_degree = polynom_degree
def get_vector(self, eval_points): def get_vector(self, eval_points):
pass pass
...@@ -20,7 +20,7 @@ class Legendre(Vector): ...@@ -20,7 +20,7 @@ class Legendre(Vector):
def _calculate_legendre_vector(self, eval_point): def _calculate_legendre_vector(self, eval_point):
vector = [] vector = []
for degree in range(self.polynom_degree+1): for degree in range(self._polynom_degree+1):
if degree == 0: if degree == 0:
vector.append(1.0 + 0*eval_point) vector.append(1.0 + 0*eval_point)
else: else:
...@@ -35,12 +35,12 @@ class Legendre(Vector): ...@@ -35,12 +35,12 @@ class Legendre(Vector):
class OrthonormalLegendre(Legendre): class OrthonormalLegendre(Legendre):
def get_vector(self, eval_point): def get_vector(self, eval_point):
leg_vector = self._calculate_legendre_vector(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): class AlpertsWavelet(Vector):
def get_vector(self, eval_point): def get_vector(self, eval_point):
degree = self.polynom_degree degree = self._polynom_degree
if degree == 0: if degree == 0:
return [np.sqrt(0.5) + eval_point*0] return [np.sqrt(0.5) + eval_point*0]
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment