diff --git a/Update_Scheme.py b/Update_Scheme.py
new file mode 100644
index 0000000000000000000000000000000000000000..dfb518400b808255c0805bbfcab3a55e75ef7f04
--- /dev/null
+++ b/Update_Scheme.py
@@ -0,0 +1,235 @@
+# -*- coding: utf-8 -*-
+"""
+@author: Laura C. Kühle
+
+d = detail coefficient (rename?)
+other A (from M) = ? (Is it the same???)
+A = basis_projection_left
+M1 = wavelet_projection_left
+phi = DG basis vector
+psi = wavelet vector
+
+TODO: Find better names for A, B, M1, and M2
+TODO: Contemplate how to handle closing brackets \
+ (end of line?, indented in new line?) -> Done (in line seems to be fine)
+TODO: Contemplate giving option to change number of iterations for history \
+ -> Done (history_threshold)
+TODO: Rename counter -> Done (iteration)
+TODO: Investigate whether list comprehension improves runtime \
+ (compared to for-loops) -> Done (not really but more compact)
+"""
+import numpy as np
+from sympy import Symbol, integrate
+
+from Vectors_of_Polynomials import OrthonormalLegendre, AlpertsWavelet
+
+x = Symbol('x')
+xi = Symbol('z')
+
+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._reset()
+ pass
+
+ def get_name(self):
+ return self.name
+
+ def get_troubled_cell_history(self):
+ return self.troubled_cell_history
+
+ def get_time_history(self):
+ return self.time_history
+
+ def get_wavelet_coeffs(self, projection=None):
+ if projection is None:
+ return self.multiwavelet_coeffs
+ else:
+ return self._calculate_wavelet_coeffs(projection)
+
+ def step(self, projection, cfl_number, time):
+ self.original_projection = projection
+ self.current_projection = projection
+ self.cfl_number = cfl_number
+ self.time = time
+
+ self._apply_stability_method()
+
+ self.iteration += 1
+
+ if (self.iteration % self.history_threshold) == 0:
+ self.troubled_cell_history.append(self.troubled_cells)
+ self.time_history.append(self.time)
+
+ return self.current_projection, self.troubled_cells
+
+ def _reset(self):
+ # Set fixed basis and wavelet vectors
+ self.basis = OrthonormalLegendre(
+ self.polynom_degree).get_vector(x)
+ self.wavelet = AlpertsWavelet(self.polynom_degree).get_vector(x)
+
+ # Set additional necessary fixed instance variables
+ self.interval_len = self.right_bound-self.left_bound
+ self.cell_len = self.interval_len / self.num_grid_cells
+ self.M1 = self._set_multiwavelet_matrix(xi, -0.5*(xi-1), True)
+ self.M2 = self._set_multiwavelet_matrix(xi, 0.5*(xi+1), False)
+
+ # Set matrix A
+ matrix = []
+ for i in range(self.polynom_degree+1):
+ new_row = []
+ 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)
+
+ # Set matrix B
+ matrix = []
+ for i in range(self.polynom_degree+1):
+ new_row = []
+ 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)
+
+ # Initialize temporary instance variables
+ self.original_projection = []
+ self.current_projection = []
+ self.right_hand_side = []
+ self.multiwavelet_coeffs = []
+ self.troubled_cells = []
+ self.troubled_cell_history = []
+ self.time_history = []
+ self.cfl_number = 0
+ self.time = 0
+ self.iteration = 0
+
+ def _set_multiwavelet_matrix(self, first_param, second_param, is_M1):
+ matrix = []
+ 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),
+ (xi, -1, 1))
+ if (is_M1):
+ entry = entry * (-1)**(j + self.polynom_degree + 1)
+ row.append(np.float64(entry))
+ matrix.append(row)
+ return matrix
+
+ def _update_wavelet_coeffs(self):
+ projection = self.current_projection[:, 1: -1]
+ self.multiwavelet_coeffs = self._calculate_wavelet_coeffs(projection)
+
+ def _calculate_wavelet_coeffs(self, projection):
+ transposed_vector = np.transpose(projection)
+ output_matrix = []
+ for i in range(int(len(projection[0])/2)):
+ new_entry = 0.5*(transposed_vector[2*i] @ self.M1
+ + transposed_vector[2*i+1] @ self.M2)
+ output_matrix.append(new_entry)
+ return np.transpose(np.array(output_matrix))
+
+ def _apply_limiter(self):
+ self._update_wavelet_coeffs()
+ self.troubled_cells = self.detector.get_cells(self.multiwavelet_coeffs,
+ self.current_projection)
+
+ new_projection = self.current_projection.copy()
+ for cell in self.troubled_cells:
+ np.transpose(new_projection)[cell] = self.limiter.apply(
+ self.current_projection, cell)
+
+ self.current_projection = new_projection
+
+ def _enforce_boundary_condition(self):
+ transposed_projection = np.transpose(self.current_projection)
+
+ transposed_projection[0] = transposed_projection[self.num_grid_cells]
+ transposed_projection[self.num_grid_cells+1] = transposed_projection[1]
+
+ self.current_projection = np.transpose(transposed_projection)
+
+
+class SSPRK3(UpdateScheme):
+
+ def __init__(self, detector, limiter, init_cond, mesh, wave_speed,
+ polynom_degree, num_grid_cells, final_time, history_threshold,
+ left_bound, right_bound):
+ # Set name of update scheme
+ self.name = 'SSPRK3'
+
+ super().__init__(detector, limiter, init_cond, mesh, wave_speed,
+ polynom_degree, num_grid_cells, final_time,
+ history_threshold, left_bound, right_bound)
+
+ # Override method of superclass
+ def _apply_stability_method(self):
+ self._apply_first_step()
+ self._apply_limiter()
+ self._enforce_boundary_condition()
+
+ self._apply_second_step()
+ self._apply_limiter()
+ self._enforce_boundary_condition()
+
+ self._apply_third_step()
+ self._apply_limiter()
+ self._enforce_boundary_condition()
+
+ def _update_right_hand_side(self):
+ # Transpose projection for easier calculation
+ transposed_projection = np.transpose(self.current_projection)
+
+ # Initialze 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 @ transposed_projection[j+1]
+ + self.B @ transposed_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)
+
+ 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_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_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))