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

Contained cell length in mesh.

parent cbf98d27
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
DG_Approximation.py
+ 30
− 26
View file @ ec4d6b35
......@@ -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])
......
Plotting.py
+ 2
− 9
View file @ ec4d6b35
......@@ -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(),
......
Troubled_Cell_Detector.py
+ 2
− 5
View file @ ec4d6b35
......@@ -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
projection_utils.py
+ 3
− 4
View file @ ec4d6b35
......@@ -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:
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment