diff --git a/DG_Approximation.py b/DG_Approximation.py
index 0e58a6965c14ba5178478636aa39f7bd4ef96950..938de9305433bea6d1a8efd7c69575b693ef971f 100644
--- a/DG_Approximation.py
+++ b/DG_Approximation.py
@@ -9,17 +9,19 @@ TODO: Contemplate separating cell average and reconstruction calculations
TODO: Contemplate removing Methods section from class docstring
TODO: Contemplate containing the quadrature application for plots in Mesh
TODO: Contemplate containing coarse mesh generation in Mesh
+TODO: Contemplate extracting boundary condition from InitialCondition
Urgent:
-TODO: Put basis initialization for plots in function
-TODO: Contain cell length in mesh
+TODO: Put basis initialization for plots in function -> Done
+TODO: Contain cell length in mesh -> Done
TODO: Contain bounds in mesh
TODO: Contain number of grid cells in mesh
TODO: Contain interval length in mesh
TODO: Create data dict for mesh separately
+TODO: Create data dict for basis separately
+TODO: Refactor eval_points in do_initial_projection()
TODO: Check whether ghost cells are handled/set correctly
TODO: Ensure uniform use of mesh and grid
-TODO: Check whether eval_point in initial projection is set correctly -> Done
TODO: Replace getter with property attributes for quadrature
TODO: Remove use of DGScheme from ANN_Data_Generator
TODO: Find error in centering for ANN training
@@ -27,12 +29,15 @@ TODO: Adapt TCD from Soraya
(Dropbox->...->TEST_troubled-cell-detector->Troubled_Cell_Detector)
TODO: Add TC condition to only flag cell if left-adjacent one is flagged as
well
+TODO: Move plot_approximation_results() into plotting script
TODO: Add verbose output
TODO: Improve file naming (e.g. use '.' instead of '__')
TODO: Combine ANN workflows
TODO: Add an environment file for Snakemake
Critical, but not urgent:
+TODO: Investigate profiling for speed up
+TODO: Rework Theoretical TC for efficiency
TODO: Check sign change in stiffness matrix to accommodate negative wave speed
TODO: Investigate g-mesh(?)
TODO: Create g-mesh with Mesh class
@@ -93,8 +98,6 @@ class DGScheme:
----------
interval_len : float
Length of the interval between left and right boundary.
- cell_len : float
- Length of a cell in mesh.
basis : Basis object
Basis for calculation.
mesh : Mesh
@@ -239,11 +242,13 @@ class DGScheme:
"""
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)
+ initial_condition=self._init_cond, mesh=self._mesh,
+ basis=self._basis, quadrature=self._quadrature,
+ num_grid_cells=self._num_grid_cells,
+ left_bound=self._left_bound)
- time_step = abs(self._cfl_number * self._cell_len / self._wave_speed)
+ time_step = abs(self._cfl_number * self._mesh.cell_len /
+ self._wave_speed)
current_time = 0
iteration = 0
@@ -254,7 +259,7 @@ class DGScheme:
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._cell_len
+ cfl_number = self._wave_speed * time_step / self._mesh.cell_len
# Update projection
projection, troubled_cells = self._update_scheme.step(projection,
@@ -287,10 +292,7 @@ class DGScheme:
"""Resets instance variables."""
# Set additional necessary instance variables
self._interval_len = self._right_bound-self._left_bound
- self._cell_len = self._interval_len / self._num_grid_cells
-
- # Set additional necessary config parameters
- self._limiter_config['cell_len'] = self._cell_len
+ # self._cell_len = self._interval_len / self._num_grid_cells
# Initialize mesh with two ghost cells on each side
self._mesh = Mesh(num_grid_cells=self._num_grid_cells,
@@ -299,6 +301,9 @@ class DGScheme:
print(len(self._mesh.cells))
print(type(self._mesh.cells))
+ # Set additional necessary config parameters
+ self._limiter_config['cell_len'] = self._mesh.cell_len
+
def build_training_data(self, adjustment, stencil_length,
add_reconstructions, initial_condition=None):
"""Builds training data set.
@@ -328,19 +333,18 @@ class DGScheme:
if initial_condition is None:
initial_condition = self._init_cond
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,
- adjustment=adjustment)
+ initial_condition=initial_condition, mesh=self._mesh,
+ basis=self._basis, quadrature=self._quadrature,
+ num_grid_cells=self._num_grid_cells,
+ left_bound=self._left_bound, adjustment=adjustment)
return self._basis.calculate_cell_average(
projection=projection[:, 1:-1], stencil_length=stencil_length,
add_reconstructions=add_reconstructions)
-def do_initial_projection(initial_condition, basis, quadrature,
- num_grid_cells, left_bound, right_bound,
- adjustment=0):
+def do_initial_projection(initial_condition, mesh, basis, quadrature,
+ num_grid_cells, left_bound, adjustment=0):
"""Calculates initial projection.
Calculates a projection at time step 0 and adds ghost cells on both
@@ -350,6 +354,8 @@ def do_initial_projection(initial_condition, basis, quadrature,
----------
initial_condition : InitialCondition object
Initial condition used for calculation.
+ mesh : Mesh
+ Mesh for calculation.
basis: Basis object
Basis used for calculation.
quadrature: Quadrature object
@@ -358,8 +364,6 @@ def do_initial_projection(initial_condition, basis, quadrature,
Number of cells in the mesh. Usually exponential of 2.
left_bound : float
Left boundary of interval.
- right_bound : float
- Right boundary of interval.
adjustment: float, optional
Extent of adjustment of each evaluation point in x-direction.
Default: 0.
@@ -375,14 +379,14 @@ def do_initial_projection(initial_condition, basis, quadrature,
output_matrix = [0]
basis_vector = basis.basis
- cell_len = (right_bound-left_bound)/num_grid_cells
+ # 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
+ eval_point = left_bound + (cell+0.5)*mesh.cell_len
for degree in range(basis.polynomial_degree + 1):
new_row.append(np.float64(sum(initial_condition.calculate(
- eval_point + cell_len/2
+ eval_point + mesh.cell_len/2
* quadrature.get_eval_points()[point] - adjustment)
* basis_vector[degree].subs(
x, quadrature.get_eval_points()[point])
diff --git a/Plotting.py b/Plotting.py
index cfed1efae84b5312a87bc3b9bd42bab8301ddef5..f58d547e4d360675b8fed64b54aee6e91310994b 100644
--- a/Plotting.py
+++ b/Plotting.py
@@ -425,15 +425,13 @@ def plot_results(projection: ndarray, troubled_cell_history: list,
# Calculate needed variables
interval_len = right_bound-left_bound
- cell_len = interval_len/num_grid_cells
# Plot troubled cells
plot_shock_tube(num_grid_cells, troubled_cell_history, time_history)
# Determine exact and approximate solution
grid, exact = calculate_exact_solution(
- mesh, cell_len, wave_speed,
- final_time, interval_len, quadrature,
+ mesh, wave_speed, final_time, interval_len, quadrature,
init_cond)
approx = calculate_approximate_solution(
projection[:, 1:-1], quadrature.get_eval_points(),
@@ -441,18 +439,13 @@ def plot_results(projection: ndarray, troubled_cell_history: list,
# Plot multiwavelet solution (fine and coarse grid)
if coarse_projection is not None:
- coarse_cell_len = 2*cell_len
coarse_mesh = Mesh(num_grid_cells=num_grid_cells//2,
num_ghost_cells=1, left_bound=left_bound,
right_bound=right_bound)
- # coarse_mesh = np.arange(left_bound - (0.5*coarse_cell_len),
- # right_bound + (1.5*coarse_cell_len),
- # coarse_cell_len)
# Plot exact and approximate solutions for coarse mesh
coarse_grid, coarse_exact = calculate_exact_solution(
- coarse_mesh, coarse_cell_len, wave_speed,
- final_time, interval_len, quadrature,
+ coarse_mesh, wave_speed, final_time, interval_len, quadrature,
init_cond)
coarse_approx = calculate_approximate_solution(
coarse_projection, quadrature.get_eval_points(),
diff --git a/Troubled_Cell_Detector.py b/Troubled_Cell_Detector.py
index ed420b966214d58067ed4948d85d485d22aa3d83..209dee40d3b2599ab3c0715e520b6b98566deb8f 100644
--- a/Troubled_Cell_Detector.py
+++ b/Troubled_Cell_Detector.py
@@ -25,8 +25,6 @@ class TroubledCellDetector(ABC):
----------
interval_len : float
Length of the interval between left and right boundary.
- cell_len : float
- Length of a cell in mesh.
Methods
-------
@@ -72,7 +70,6 @@ class TroubledCellDetector(ABC):
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
@@ -468,7 +465,7 @@ class Theoretical(WaveletDetector):
# Unpack necessary configurations
self._cutoff_factor = config.pop('cutoff_factor',
- np.sqrt(2) * self._cell_len)
+ np.sqrt(2) * self._mesh.cell_len)
# comment to line above: or 2 or 3
def _get_cells(self, multiwavelet_coeffs, projection):
@@ -520,6 +517,6 @@ class Theoretical(WaveletDetector):
for degree in range(
self._basis.polynomial_degree+1))/max_avg
eps = self._cutoff_factor\
- / (self._cell_len*self._num_coarse_grid_cells*2)
+ / (self._mesh.cell_len*self._num_coarse_grid_cells*2)
return max_value > eps
diff --git a/projection_utils.py b/projection_utils.py
index 41d66544fa8498620a86238f86f77948fb11f052..b4df91f359ca12446c945a2a72e326e4b37ecc27 100644
--- a/projection_utils.py
+++ b/projection_utils.py
@@ -131,7 +131,7 @@ def calculate_approximate_solution(
def calculate_exact_solution(
- mesh: Mesh, cell_len: float, wave_speed: float, final_time:
+ mesh: Mesh, wave_speed: float, final_time:
float, interval_len: float, quadrature: Quadrature, init_cond:
InitialCondition) -> Tuple[ndarray, ndarray]:
"""Calculate exact solution.
@@ -140,8 +140,6 @@ def calculate_exact_solution(
----------
mesh : Mesh
Mesh for evaluation.
- cell_len : float
- Length of a cell in mesh.
wave_speed : float
Speed of wave in rightward direction.
final_time : float
@@ -166,7 +164,8 @@ def calculate_exact_solution(
num_periods = np.floor(wave_speed * final_time / interval_len)
for cell_center in mesh.non_ghost_cells:
- eval_points = cell_center+cell_len / 2 * quadrature.get_eval_points()
+ eval_points = cell_center+mesh.cell_len / 2 * \
+ quadrature.get_eval_points()
eval_values = []
for eval_point in eval_points: