diff --git a/DG_Approximation.py b/DG_Approximation.py
index b785c93801208ebf6997283dc5708a53391c3baa..ee7ed1e14b9eb80134631d0e303949e741b6866b 100644
--- a/DG_Approximation.py
+++ b/DG_Approximation.py
@@ -17,8 +17,8 @@ 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: Make sure all instance variables are actually necessary
-TODO: Make sure instance variables are only set in __init__()
+TODO: Make sure all instance variables are actually necessary -> Done
+TODO: Make sure instance variables are only set in __init__() -> Done
TODO: Contemplate moving plots to pertaining files
TODO: Contemplate moving A and B to Vectors_of_Polynomials
TODO: Combine plot for coarse and fine approximation for wavelet detectors
diff --git a/Limiter.py b/Limiter.py
index 8cf1501cfc6585bf1708fe9819c31d76c61f08e9..a5033e70fb44d3cb9a519fb8f90ed8a1db3b0e9a 100644
--- a/Limiter.py
+++ b/Limiter.py
@@ -37,35 +37,34 @@ class MinMod(Limiter):
# Set name of function
self.function_name = 'MinMod' + str(self.erase_degree)
- self.cell = 0
-
def apply(self, projection, cell):
- self.cell = cell
- self._set_cell_slope(projection)
- is_good_cell = self._determine_modification(projection)
+ cell_slope = self._set_cell_slope(projection, cell)
+ is_good_cell = self._determine_modification(projection, cell, cell_slope)
if is_good_cell:
- return projection[:, self.cell]
+ return projection[:, cell]
- adapted_projection = projection[:, self.cell].copy()
+ adapted_projection = projection[:, cell].copy()
for i in range(len(adapted_projection)):
if i > self.erase_degree:
adapted_projection[i] = 0
return adapted_projection
- def _set_cell_slope(self, projection):
+ @staticmethod
+ def _set_cell_slope(projection, cell):
slope = []
- for cell in range(len(projection[0])):
- new_entry = sum(projection[degree][cell] * (degree+0.5)**0.5 for degree in range(1, len(projection)))
+ for current_cell in range(len(projection[0])):
+ new_entry = sum(projection[degree][current_cell] * (degree+0.5)**0.5
+ for degree in range(1, len(projection)))
slope.append(new_entry)
- self.cell_slope = slope[self.cell]
+ return slope[cell]
- def _determine_modification(self, projection):
- forward_slope = (projection[0][self.cell+1] - projection[0][self.cell]) * (0.5**0.5)
- backward_slope = (projection[0][self.cell] - projection[0][self.cell-1]) * (0.5**0.5)
+ def _determine_modification(self, projection, cell, cell_slope):
+ forward_slope = (projection[0][cell+1] - projection[0][cell]) * (0.5**0.5)
+ backward_slope = (projection[0][cell] - projection[0][cell-1]) * (0.5**0.5)
- return (forward_slope <= 0) & (backward_slope <= 0) & (self.cell_slope <= 0) \
- | (forward_slope >= 0) & (backward_slope >= 0) & (self.cell_slope >= 0)
+ return (forward_slope <= 0) & (backward_slope <= 0) & (cell_slope <= 0) \
+ | (forward_slope >= 0) & (backward_slope >= 0) & (cell_slope >= 0)
class ModifiedMinMod(MinMod):
@@ -81,8 +80,8 @@ class ModifiedMinMod(MinMod):
self.cell = 0
- def _determine_modification(self, projection):
- if abs(self.cell_slope) <= self.mod_factor*self.cell_len**2:
+ def _determine_modification(self, projection, cell, cell_slope):
+ if abs(cell_slope) <= self.mod_factor*self.cell_len**2:
return True
- return super()._determine_modification(projection)
+ return super()._determine_modification(projection, cell, cell_slope)
diff --git a/Troubled_Cell_Detector.py b/Troubled_Cell_Detector.py
index 9953b7f908231057d9990b431803a41345f6716d..f4d1d7bcbc5ad1a9e4a1075c46bfcc61827cfe51 100644
--- a/Troubled_Cell_Detector.py
+++ b/Troubled_Cell_Detector.py
@@ -144,23 +144,19 @@ class Theoretical(WaveletDetector):
self.cutoff_factor = config.pop('cutoff_factor', np.sqrt(2) * self.cell_len)
# comment to line above: or 2 or 3
- self.multiwavelet_coeffs = []
- self.max_avg = None
-
def _get_cells(self, multiwavelet_coeffs, projection):
- self.multiwavelet_coeffs = multiwavelet_coeffs
troubled_cells = []
- self.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):
- if self._is_troubled_cell(cell):
+ if self._is_troubled_cell(multiwavelet_coeffs, cell, max_avg):
troubled_cells.append(cell)
return troubled_cells
- 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
+ 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)
return max_value > eps
diff --git a/Update_Scheme.py b/Update_Scheme.py
index e0f56544fd8bc144713f7c3aa7f7a33a333f6efc..a3eea802a526321cee5e8589c1ef1aafc165f428 100644
--- a/Update_Scheme.py
+++ b/Update_Scheme.py
@@ -44,29 +44,27 @@ class UpdateScheme(object):
return self.time_history
def step(self, projection, cfl_number, current_time):
- self.original_projection = projection
- self.current_projection = projection
- self.cfl_number = cfl_number
- self.time = current_time
-
- self._apply_stability_method()
+ current_projection, troubled_cells = self._apply_stability_method(projection, cfl_number)
self.iteration += 1
if (self.iteration % self.history_threshold) == 0:
- self.troubled_cell_history.append(self.troubled_cells)
- self.time_history.append(self.time)
+ self.troubled_cell_history.append(troubled_cells)
+ self.time_history.append(current_time)
- return self.current_projection, self.troubled_cells
+ return current_projection, troubled_cells
- def _apply_stability_method(self):
- pass
+ def _apply_stability_method(self, projection, cfl_number):
+ return projection, []
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.troubled_cell_history = []
+ self.time_history = []
+ self.iteration = 0
# Set matrix A
matrix = []
@@ -90,29 +88,19 @@ class UpdateScheme(object):
matrix.append(new_row)
self.B = np.array(matrix) # former: inv_mass @ np.array(matrix)
- # Initialize temporary instance variables
- self.original_projection = []
- self.current_projection = []
- self.right_hand_side = []
- self.troubled_cells = []
- self.troubled_cell_history = []
- self.time_history = []
- self.cfl_number = 0
- self.time = 0
- self.iteration = 0
-
- def _apply_limiter(self):
- self.troubled_cells = self.detector.get_cells(self.current_projection)
+ def _apply_limiter(self, current_projection):
+ troubled_cells = self.detector.get_cells(current_projection)
- new_projection = self.current_projection.copy()
- for cell in self.troubled_cells:
- new_projection[:, cell] = self.limiter.apply(self.current_projection, cell)
+ new_projection = current_projection.copy()
+ for cell in troubled_cells:
+ new_projection[:, cell] = self.limiter.apply(current_projection, cell)
- self.current_projection = new_projection
+ return new_projection, troubled_cells
- def _enforce_boundary_condition(self):
- self.current_projection[:, 0] = self.current_projection[:, self.num_grid_cells]
- self.current_projection[:, self.num_grid_cells+1] = self.current_projection[:, 1]
+ 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]
+ return current_projection
class SSPRK3(UpdateScheme):
@@ -125,43 +113,45 @@ class SSPRK3(UpdateScheme):
self.name = 'SSPRK3'
# Override method of superclass
- def _apply_stability_method(self):
- self._apply_first_step()
- self._apply_limiter()
- self._enforce_boundary_condition()
+ def _apply_stability_method(self, projection, cfl_number):
+ original_projection = projection
+
+ current_projection = self._apply_first_step(original_projection, cfl_number)
+ current_projection, __ = self._apply_limiter(current_projection)
+ current_projection = self._enforce_boundary_condition(current_projection)
+
+ current_projection = self._apply_second_step(original_projection, current_projection, cfl_number)
+ current_projection, __ = self._apply_limiter(current_projection)
+ current_projection = self._enforce_boundary_condition(current_projection)
- self._apply_second_step()
- self._apply_limiter()
- self._enforce_boundary_condition()
+ current_projection = self._apply_third_step(original_projection, current_projection, cfl_number)
+ current_projection, troubled_cells = self._apply_limiter(current_projection)
+ current_projection = self._enforce_boundary_condition(current_projection)
- self._apply_third_step()
- self._apply_limiter()
- self._enforce_boundary_condition()
+ return current_projection, troubled_cells
- def _update_right_hand_side(self):
+ def _update_right_hand_side(self, current_projection):
# 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 @ self.current_projection[:, j+1]
- + self.B @ self.current_projection[:, j]))
+ 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.append(right_hand_side[1])
- self.right_hand_side = np.transpose(right_hand_side)
+ return np.transpose(right_hand_side)
- def _apply_first_step(self):
- self._update_right_hand_side()
- self.current_projection = self.original_projection + (self.cfl_number*self.right_hand_side)
+ def _apply_first_step(self, original_projection, cfl_number):
+ right_hand_side = self._update_right_hand_side(original_projection)
+ return original_projection + (cfl_number*right_hand_side)
- def _apply_second_step(self):
- self._update_right_hand_side()
- self.current_projection = 1/4 * (3 * self.original_projection
- + (self.current_projection + self.cfl_number*self.right_hand_side))
+ def _apply_second_step(self, original_projection, current_projection, cfl_number):
+ right_hand_side = self._update_right_hand_side(current_projection)
+ return 1/4 * (3*original_projection + (current_projection + cfl_number*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))
+ def _apply_third_step(self, original_projection, current_projection, cfl_number):
+ right_hand_side = self._update_right_hand_side(current_projection)
+ return 1/3 * (original_projection + 2*(current_projection + cfl_number*right_hand_side))