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

Contained polynomial degree in basis.

parent ed0a0574
No related branches found
No related tags found
No related merge requests found
......@@ -21,6 +21,8 @@ class Basis:
Attributes
----------
polynomial degree : int
Polynomial degree.
basis : ndarray
Array of basis.
wavelet : ndarray
......@@ -47,6 +49,11 @@ class Basis:
"""
self._polynomial_degree = polynomial_degree
@property
def polynomial_degree(self):
"""Return polynomial degree."""
return self._polynomial_degree
@property
@cache
def basis(self):
......
......@@ -7,13 +7,15 @@ TODO: Contemplate saving 5-CV split and evaluating models separately
TODO: Contemplate separating cell average and reconstruction calculations
completely
TODO: Contemplate removing Methods section from class docstring
TODO: Contemplate moving stiffness and boundary matrix from update scheme
to basis
Urgent:
TODO: Investigate protected vs. public variables -> Done
TODO: Investigate decorators and properties -> Done
TODO: Investigate how to make variables read-only -> Done
TODO: Replace getter with property attributes for basis -> Done
TODO: Contain polynomial degree in basis
TODO: Contain polynomial degree in basis -> Done
TODO: Improve docstring for Basis class
TODO: Enforce abstract classes/methods (abc.ABC, abc.abstractmethod)
TODO: Introduce Mesh class
......@@ -154,7 +156,6 @@ class DGScheme:
"""
# Unpack keyword arguments
self._wave_speed = kwargs.pop('wave_speed', 1)
self._polynomial_degree = kwargs.pop('polynomial_degree', 2)
self._cfl_number = kwargs.pop('cfl_number', 0.2)
self._num_grid_cells = kwargs.pop('num_grid_cells', 64)
self._final_time = kwargs.pop('final_time', 1)
......@@ -173,6 +174,9 @@ class DGScheme:
self._quadrature_config = kwargs.pop('quadrature_config', {})
self._update_scheme = kwargs.pop('update_scheme', 'SSPRK3')
polynomial_degree = kwargs.pop('polynomial_degree', 2)
self._basis = OrthonormalLegendre(polynomial_degree)
# Throw an error if there are extra keyword arguments
if len(kwargs) > 0:
extra = ', '.join('"%s"' % k for k in list(kwargs.keys()))
......@@ -206,12 +210,11 @@ class DGScheme:
self._detector = getattr(Troubled_Cell_Detector, self._detector)(
config=self._detector_config, mesh=self._mesh,
wave_speed=self._wave_speed, num_grid_cells=self._num_grid_cells,
polynomial_degree=self._polynomial_degree,
final_time=self._final_time, left_bound=self._left_bound,
right_bound=self._right_bound, basis=self._basis,
init_cond=self._init_cond, quadrature=self._quadrature)
self._update_scheme = getattr(Update_Scheme, self._update_scheme)(
polynomial_degree=self._polynomial_degree,
polynomial_degree=polynomial_degree,
num_grid_cells=self._num_grid_cells, detector=self._detector,
limiter=self._limiter)
......@@ -233,8 +236,7 @@ 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,
polynomial_degree=self._polynomial_degree)
left_bound=self._left_bound, right_bound=self._right_bound)
time_step = abs(self._cfl_number * self._cell_len / self._wave_speed)
......@@ -281,7 +283,6 @@ class DGScheme:
# Set additional necessary instance variables
self._interval_len = self._right_bound-self._left_bound
self._cell_len = self._interval_len / self._num_grid_cells
self._basis = OrthonormalLegendre(self._polynomial_degree)
# Set additional necessary config parameters
self._limiter_config['cell_len'] = self._cell_len
......@@ -323,7 +324,7 @@ class DGScheme:
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)
adjustment=adjustment)
return self._basis.calculate_cell_average(
projection=projection[:, 1:-1], stencil_length=stencil_length,
......@@ -332,7 +333,7 @@ class DGScheme:
def do_initial_projection(initial_condition, basis, quadrature,
num_grid_cells, left_bound, right_bound,
polynomial_degree, adjustment=0):
adjustment=0):
"""Calculates initial projection.
Calculates a projection at time step 0 and adds ghost cells on both
......@@ -352,8 +353,6 @@ def do_initial_projection(initial_condition, basis, quadrature,
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.
......@@ -362,7 +361,7 @@ def do_initial_projection(initial_condition, basis, quadrature,
-------
ndarray
Matrix containing projection of size (N+2, p+1) with N being the
number of grid cells and p being the polynomial degree.
number of grid cells and p being the polynomial degree of the basis.
"""
# Initialize matrix and set first entry to accommodate for ghost cell
......@@ -374,7 +373,7 @@ def do_initial_projection(initial_condition, basis, quadrature,
new_row = []
eval_point = left_bound + (cell+0.5)*cell_len
for degree in range(polynomial_degree + 1):
for degree in range(basis.polynomial_degree + 1):
new_row.append(np.float64(sum(initial_condition.calculate(
eval_point + cell_len/2
* quadrature.get_eval_points()[point] - adjustment)
......
......@@ -124,10 +124,9 @@ def plot_shock_tube(num_grid_cells: int, troubled_cell_history: list,
plt.title('Shock Tubes')
def plot_details(fine_projection: ndarray, fine_mesh: ndarray,
coarse_projection: ndarray, basis: ndarray, wavelet: ndarray,
multiwavelet_coeffs: ndarray, num_coarse_grid_cells: int,
polynomial_degree: int) -> None:
def plot_details(fine_projection: ndarray, fine_mesh: ndarray, basis: Basis,
coarse_projection: ndarray, multiwavelet_coeffs: ndarray,
num_coarse_grid_cells: int) -> None:
"""Plots details of projection to coarser mesh.
Parameters
......@@ -136,35 +135,32 @@ def plot_details(fine_projection: ndarray, fine_mesh: ndarray,
Matrix of projection for each polynomial degree.
fine_mesh : ndarray
List of evaluation points for fine mesh.
basis : ndarray
Basis vector for calculation.
wavelet : ndarray
Wavelet vector for calculation.
basis: Basis object
Basis used for calculation.
multiwavelet_coeffs : ndarray
Matrix of multiwavelet coefficients.
num_coarse_grid_cells : int
Number of cells in the coarse mesh (half the cells of the fine mesh).
Usually exponential of 2.
polynomial_degree : int
Polynomial degree.
"""
averaged_projection = [[coarse_projection[degree][cell]
* basis[degree].subs(x, value)
* basis.basis[degree].subs(x, value)
for cell in range(num_coarse_grid_cells)
for value in [-0.5, 0.5]]
for degree in range(polynomial_degree + 1)]
for degree in range(basis.polynomial_degree + 1)]
wavelet_projection = [[multiwavelet_coeffs[degree][cell]
* wavelet[degree].subs(z, 0.5) * value
* basis.wavelet[degree].subs(z, 0.5) * value
for cell in range(num_coarse_grid_cells)
for value in [(-1) ** (polynomial_degree
for value in [(-1) ** (basis.polynomial_degree
+ degree + 1), 1]]
for degree in range(polynomial_degree + 1)]
for degree in range(basis.polynomial_degree + 1)]
projected_coarse = np.sum(averaged_projection, axis=0)
projected_fine = np.sum([fine_projection[degree] * basis[degree].subs(x, 0)
for degree in range(polynomial_degree + 1)],
projected_fine = np.sum([fine_projection[degree]
* basis.basis[degree].subs(x, 0)
for degree in range(basis.polynomial_degree + 1)],
axis=0)
projected_wavelet_coeffs = np.sum(wavelet_projection, axis=0)
......@@ -352,6 +348,7 @@ def plot_approximation_results(data_file: str, directory: str, plot_name: str,
# Decode all ndarrays by converting lists
approx_stats = {key: decode_ndarray(approx_stats[key])
for key in approx_stats.keys()}
approx_stats.pop('polynomial_degree')
# Plot exact/approximate results, errors, shock tubes,
# and any detector-dependant plots
......@@ -375,7 +372,7 @@ def plot_approximation_results(data_file: str, directory: str, plot_name: str,
def plot_results(projection: ndarray, troubled_cell_history: list,
time_history: list, mesh: ndarray, num_grid_cells: int,
polynomial_degree: int, wave_speed: float, final_time: float,
wave_speed: float, final_time: float,
left_bound: float, right_bound: float, basis: Basis,
quadrature: Quadrature, init_cond: InitialCondition,
colors: dict = None, coarse_projection: ndarray = None,
......@@ -400,8 +397,6 @@ def plot_results(projection: ndarray, troubled_cell_history: list,
List of mesh valuation points.
num_grid_cells : int
Number of cells in the mesh. Usually exponential of 2.
polynomial_degree : int
Polynomial degree.
wave_speed : float
Speed of wave in rightward direction.
final_time : float
......@@ -444,7 +439,7 @@ def plot_results(projection: ndarray, troubled_cell_history: list,
init_cond)
approx = calculate_approximate_solution(
projection[:, 1:-1], quadrature.get_eval_points(),
polynomial_degree, basis.basis)
basis.polynomial_degree, basis.basis)
# Plot multiwavelet solution (fine and coarse grid)
if coarse_projection is not None:
......@@ -460,17 +455,15 @@ def plot_results(projection: ndarray, troubled_cell_history: list,
init_cond)
coarse_approx = calculate_approximate_solution(
coarse_projection, quadrature.get_eval_points(),
polynomial_degree, basis.basis)
basis.polynomial_degree, basis.basis)
plot_solution_and_approx(
coarse_grid, coarse_exact, coarse_approx, colors['coarse_exact'],
colors['coarse_approx'])
# Plot multiwavelet details
num_coarse_grid_cells = num_grid_cells//2
plot_details(projection[:, 1:-1], mesh[2:-2], coarse_projection,
basis.basis, basis.wavelet, multiwavelet_coeffs,
num_coarse_grid_cells,
polynomial_degree)
plot_details(projection[:, 1:-1], mesh[2:-2], basis, coarse_projection,
multiwavelet_coeffs, num_coarse_grid_cells)
plot_solution_and_approx(grid, exact, approx,
colors['fine_exact'],
......@@ -491,7 +484,7 @@ def plot_results(projection: ndarray, troubled_cell_history: list,
plot_semilog_error(grid, pointwise_error)
plot_error(grid, exact, approx)
print('p =', polynomial_degree)
print('p =', basis.polynomial_degree)
print('N =', num_grid_cells)
print('maximum error =', max_error)
......
......@@ -37,7 +37,7 @@ class TroubledCellDetector:
"""
def __init__(self, config, init_cond, quadrature, basis, mesh,
wave_speed=1, polynomial_degree=2, num_grid_cells=64,
wave_speed=1, num_grid_cells=64,
final_time=1, left_bound=-1, right_bound=1):
"""Initializes TroubledCellDetector.
......@@ -55,8 +55,6 @@ class TroubledCellDetector:
List of mesh valuation points.
wave_speed : float, optional
Speed of wave in rightward direction. Default: 1.
polynomial_degree : int, optional
Polynomial degree. Default: 2.
num_grid_cells : int, optional
Number of cells in the mesh. Usually exponential of 2. Default: 64.
final_time : float, optional
......@@ -69,7 +67,6 @@ class TroubledCellDetector:
"""
self._mesh = mesh
self._wave_speed = wave_speed
self._polynomial_degree = polynomial_degree
self._num_grid_cells = num_grid_cells
self._final_time = final_time
self._left_bound = left_bound
......@@ -113,7 +110,7 @@ class TroubledCellDetector:
'num_grid_cells': self._num_grid_cells, 'mesh': self._mesh,
'final_time': self._final_time, 'left_bound':
self._left_bound, 'right_bound': self._right_bound,
'polynomial_degree': self._polynomial_degree
'polynomial_degree': self._basis.polynomial_degree
}
......@@ -513,7 +510,8 @@ class Theoretical(WaveletDetector):
"""
max_value = max(abs(multiwavelet_coeffs[degree][cell])
for degree in range(self._polynomial_degree+1))/max_avg
for degree in range(
self._basis.polynomial_degree+1))/max_avg
eps = self._cutoff_factor\
/ (self._cell_len*self._num_coarse_grid_cells*2)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment