diff --git a/DG_Approximation.py b/DG_Approximation.py
index d21d194a918b4c43fbf8df9b94329efd5b6cfbf7..e8487ee58702b264c596d1156dfc848cd7878f9f 100644
--- a/DG_Approximation.py
+++ b/DG_Approximation.py
@@ -3,7 +3,8 @@
@author: Laura C. Kühle
Urgent:
-TODO: Extract do_initial_projection() from DGScheme
+TODO: Extract do_initial_projection() from DGScheme -> Done
+TODO: Move inverse mass matrix to basis
TODO: Extract calculate_cell_average() from TCD
TODO: Extract calculate_[...]_solution() from Plotting
TODO: Extract plotting from TCD completely
@@ -19,7 +20,6 @@ TODO: Add an environment file for Snakemake
Critical, but not urgent:
TODO: Force input_size for each ANN model to be stencil length
-TODO: Move inverse mass matrix to basis
TODO: Use full path for ANN model state
TODO: Enforce abstract classes/methods (abc.ABC, abc.abstractmethod)
TODO: Extract object initialization from DGScheme
@@ -232,7 +232,11 @@ class DGScheme:
Path to file in which data will be saved.
"""
- projection = self._do_initial_projection(self._init_cond)
+ projection = do_initial_projection(
+ initial_condition=self._init_cond, basis=self._basis,
+ quadrature=self._quadrature, num_grid_cells=self._num_grid_cells,
+ left_bound=self._left_bound, right_bound=self._right_bound,
+ polynomial_degree=self._polynomial_degree)
time_step = abs(self._cfl_number * self._cell_len / self._wave_speed)
@@ -286,69 +290,6 @@ class DGScheme:
self._right_bound + (5/2*self._cell_len),
self._cell_len) # +3/2
- # Set inverse mass matrix
- mass_matrix = []
- for i in range(self._polynomial_degree+1):
- new_row = []
- for j in range(self._polynomial_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)
-
- def _do_initial_projection(self, initial_condition, adjustment=0):
- """Calculates initial projection.
-
- Calculates a projection at time step 0 and adds ghost cells on both
- sides of the array.
-
- Parameters
- ----------
- initial_condition : InitialCondition object
- Initial condition used for calculation. May differ from instance
- variable.
- adjustment: float
- Extent of adjustment of each evaluation point in x-direction.
-
- Returns
- -------
- ndarray
- Matrix containing projection of size (N+2, p+1) with N being the
- number of grid cells and p being the polynomial degree.
-
- """
- # Initialize matrix and set first entry to accommodate for ghost cell
- output_matrix = [0]
- basis_vector = self._basis.get_basis_vector()
-
- 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._polynomial_degree + 1):
- new_entry = sum(
- initial_condition.calculate(
- eval_point + self._cell_len/2
- * self._quadrature.get_eval_points()[point]
- - adjustment)
- * basis_vector[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)
-
- # Set ghost cells to respective value
- output_matrix[0] = output_matrix[self._num_grid_cells]
- output_matrix.append(output_matrix[1])
-
- # print(np.array(output_matrix).shape)
- return np.transpose(np.array(output_matrix))
-
def build_training_data(self, adjustment, stencil_length,
add_reconstructions, initial_condition=None):
"""Builds training data set.
@@ -377,13 +318,96 @@ class DGScheme:
"""
if initial_condition is None:
initial_condition = self._init_cond
- projection = self._do_initial_projection(initial_condition, adjustment)
+ projection = do_initial_projection(
+ initial_condition=initial_condition, basis=self._basis,
+ quadrature=self._quadrature, num_grid_cells=self._num_grid_cells,
+ left_bound=self._left_bound, right_bound=self._right_bound,
+ polynomial_degree=self._polynomial_degree, adjustment=adjustment)
return self._detector.calculate_cell_average(projection[:, 1:-1],
stencil_length,
add_reconstructions)
+def do_initial_projection(initial_condition, basis, quadrature,
+ num_grid_cells, left_bound, right_bound,
+ polynomial_degree, adjustment=0):
+ """Calculates initial projection.
+
+ Calculates a projection at time step 0 and adds ghost cells on both
+ sides of the array.
+
+ Parameters
+ ----------
+ initial_condition : InitialCondition object
+ Initial condition used for calculation.
+ basis: Vector object
+ Basis used for calculation.
+ quadrature: Quadrature object
+ Quadrature fused for evaluation.
+ num_grid_cells : int
+ Number of cells in the mesh. Usually exponential of 2.
+ left_bound : float
+ Left boundary of interval.
+ right_bound : float
+ Right boundary of interval.
+ polynomial_degree : int
+ Polynomial degree.
+ adjustment: float, optional
+ Extent of adjustment of each evaluation point in x-direction.
+ Default: 0.
+
+ Returns
+ -------
+ ndarray
+ Matrix containing projection of size (N+2, p+1) with N being the
+ number of grid cells and p being the polynomial degree.
+
+ """
+ # Set inverse mass matrix
+ mass_matrix = []
+ for i in range(polynomial_degree+1):
+ new_row = []
+ for j in range(polynomial_degree+1):
+ new_entry = 0.0
+ if i == j:
+ new_entry = 1.0
+ new_row.append(new_entry)
+ mass_matrix.append(new_row)
+ inv_mass = np.array(mass_matrix)
+
+ # Initialize matrix and set first entry to accommodate for ghost cell
+ output_matrix = [0]
+ basis_vector = basis.get_basis_vector()
+
+ cell_len = (right_bound-left_bound)/num_grid_cells
+ for cell in range(num_grid_cells):
+ new_row = []
+ eval_point = left_bound + (cell+0.5)*cell_len
+
+ for degree in range(polynomial_degree + 1):
+ new_entry = sum(
+ initial_condition.calculate(
+ eval_point + cell_len/2
+ * quadrature.get_eval_points()[point]
+ - adjustment)
+ * basis_vector[degree].subs(
+ x, quadrature.get_eval_points()[point])
+ * quadrature.get_weights()[point]
+ for point in range(quadrature.get_num_points()))
+ new_row.append(np.float64(new_entry))
+
+ new_row = np.array(new_row)
+ output_matrix.append(inv_mass @ new_row)
+
+ # Set ghost cells to respective value
+ output_matrix[0] = output_matrix[num_grid_cells]
+ output_matrix.append(output_matrix[1])
+
+ # print(np.array(output_matrix).shape)
+ return np.transpose(np.array(output_matrix))
+
+
def plot_approximation_results(detector, data_file, directory, plot_name):
"""Plots given approximation results.