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

Renamed variables and functions.

parent 75ce2e3d
No related branches found
No related tags found
No related merge requests found
...@@ -11,37 +11,37 @@ xi = Symbol('z') ...@@ -11,37 +11,37 @@ xi = Symbol('z')
class Vector(object): class Vector(object):
def __init__(self, polynom_degree): def __init__(self, polynomial_degree):
self._polynom_degree = polynom_degree self._polynomial_degree = polynomial_degree
self._basis = self._set_basis_vector(x) self._basis = self._build_basis_vector(x)
self._wavelet = self._set_wavelet_vector(xi) self._wavelet = self._build_wavelet_vector(xi)
def get_vector(self): def get_basis_vector(self):
return self._basis return self._basis
def _set_basis_vector(self, eval_point): def _build_basis_vector(self, eval_point):
return [] return []
def get_wavelet_vector(self): def get_wavelet_vector(self):
return self._wavelet return self._wavelet
def _set_wavelet_vector(self, eval_point): def _build_wavelet_vector(self, eval_point):
return [] return []
def get_basis_matrices(self): def get_basis_projections(self):
pass pass
def get_multiwavelet_matrices(self): def get_multiwavelet_projections(self):
pass pass
class Legendre(Vector): class Legendre(Vector):
def _set_basis_vector(self, eval_point): def _build_basis_vector(self, eval_point):
return self._calculate_legendre_vector(eval_point) return self._calculate_legendre_vector(eval_point)
def _calculate_legendre_vector(self, eval_point): def _calculate_legendre_vector(self, eval_point):
vector = [] vector = []
for degree in range(self._polynom_degree+1): for degree in range(self._polynomial_degree+1):
if degree == 0: if degree == 0:
vector.append(1.0 + 0*eval_point) vector.append(1.0 + 0*eval_point)
else: else:
...@@ -54,12 +54,12 @@ class Legendre(Vector): ...@@ -54,12 +54,12 @@ class Legendre(Vector):
class OrthonormalLegendre(Legendre): class OrthonormalLegendre(Legendre):
def _set_basis_vector(self, eval_point): def _build_basis_vector(self, eval_point):
leg_vector = self._calculate_legendre_vector(eval_point) leg_vector = self._calculate_legendre_vector(eval_point)
return [leg_vector[degree] * np.sqrt(degree+0.5) for degree in range(self._polynom_degree+1)] return [leg_vector[degree] * np.sqrt(degree+0.5) for degree in range(self._polynomial_degree+1)]
def _set_wavelet_vector(self, eval_point): def _build_wavelet_vector(self, eval_point):
degree = self._polynom_degree degree = self._polynomial_degree
if degree == 0: if degree == 0:
return [np.sqrt(0.5) + eval_point*0] return [np.sqrt(0.5) + eval_point*0]
...@@ -89,16 +89,16 @@ class OrthonormalLegendre(Legendre): ...@@ -89,16 +89,16 @@ class OrthonormalLegendre(Legendre):
raise ValueError('Invalid value: Alpert\'s wavelet is only available \ raise ValueError('Invalid value: Alpert\'s wavelet is only available \
up to degree 4 for this application') up to degree 4 for this application')
def get_basis_matrices(self): def get_basis_projections(self):
matrix_left = self._set_basis_matrix(xi, 0.5 * (xi - 1)) basis_projection_left = self._build_basis_matrix(xi, 0.5 * (xi - 1))
matrix_right = self._set_basis_matrix(xi, 0.5 * (xi + 1)) basis_projection_right = self._build_basis_matrix(xi, 0.5 * (xi + 1))
return matrix_left, matrix_right return basis_projection_left, basis_projection_right
def _set_basis_matrix(self, first_param, second_param): def _build_basis_matrix(self, first_param, second_param):
matrix = [] matrix = []
for i in range(self._polynom_degree + 1): for i in range(self._polynomial_degree + 1):
row = [] row = []
for j in range(self._polynom_degree + 1): for j in range(self._polynomial_degree + 1):
entry = integrate(self._basis[i].subs(x, first_param) * entry = integrate(self._basis[i].subs(x, first_param) *
self._basis[j].subs(x, second_param), self._basis[j].subs(x, second_param),
(xi, -1, 1)) (xi, -1, 1))
...@@ -106,21 +106,20 @@ class OrthonormalLegendre(Legendre): ...@@ -106,21 +106,20 @@ class OrthonormalLegendre(Legendre):
matrix.append(row) matrix.append(row)
return matrix return matrix
def get_multiwavelet_matrices(self): def get_multiwavelet_projections(self):
M1 = self._set_multiwavelet_matrix(xi, -0.5*(xi-1), True) wavelet_projection_left = self._build_multiwavelet_matrix(xi, -0.5*(xi-1), True)
M2 = self._set_multiwavelet_matrix(xi, 0.5*(xi+1), False) wavelet_projection_right = self._build_multiwavelet_matrix(xi, 0.5*(xi+1), False)
return M1, M2 return wavelet_projection_left, wavelet_projection_right
def _set_multiwavelet_matrix(self, first_param, second_param, is_M1): def _build_multiwavelet_matrix(self, first_param, second_param, is_left_matrix):
matrix = [] matrix = []
for i in range(self._polynom_degree+1): for i in range(self._polynomial_degree+1):
row = [] row = []
for j in range(self._polynom_degree+1): for j in range(self._polynomial_degree+1):
entry = integrate(self._basis[i].subs(x, first_param) * self._wavelet[j].subs(xi, second_param), entry = integrate(self._basis[i].subs(x, first_param) * self._wavelet[j].subs(xi, second_param),
(xi, -1, 1)) (xi, -1, 1))
if is_M1: if is_left_matrix:
entry = entry * (-1)**(j + self._polynom_degree + 1) entry = entry * (-1)**(j + self._polynomial_degree + 1)
row.append(np.float64(entry)) row.append(np.float64(entry))
matrix.append(row) matrix.append(row)
return matrix return matrix
...@@ -7,10 +7,14 @@ TODO: Contemplate using Seaborn instead of matplotlib ...@@ -7,10 +7,14 @@ TODO: Contemplate using Seaborn instead of matplotlib
TODO: Double-check everything! TODO: Double-check everything!
TODO: Replace loops with list comprehension if feasible TODO: Replace loops with list comprehension if feasible
TODO: Find better names for A, B, anything pertaining Gauss TODO: Find better names for A, B, anything pertaining Gauss -> Done
TODO: Write documentation for all methods TODO: Write documentation for all methods
TODO: Restructure Vectors_of_Polynomials -> Done TODO: Restructure Vectors_of_Polynomials -> Done
TODO: Rename Vectors_of_Polynomials -> Done (Basis_Function)
TODO: Ask about Legendre as scaling vector
TODO: Ask whether basis and wavelet should both depend on x (or x and z)
TODO: Ask about renaming ModifiedMinMod to Cockburn-Shu (like Gerhard et al.)
TODO: Contemplate how to make shock tubes comparable TODO: Contemplate how to make shock tubes comparable
TODO: Implement type check for all kwargs and configs TODO: Implement type check for all kwargs and configs
...@@ -25,7 +29,7 @@ import Initial_Condition ...@@ -25,7 +29,7 @@ import Initial_Condition
import Limiter import Limiter
import Quadrature import Quadrature
import Update_Scheme import Update_Scheme
from Vectors_of_Polynomials import OrthonormalLegendre from Basis_Function import OrthonormalLegendre
x = Symbol('x') x = Symbol('x')
xi = Symbol('z') xi = Symbol('z')
...@@ -35,7 +39,7 @@ class DGScheme(object): ...@@ -35,7 +39,7 @@ class DGScheme(object):
def __init__(self, detector, **kwargs): def __init__(self, detector, **kwargs):
# Unpack keyword arguments # Unpack keyword arguments
self._wave_speed = kwargs.pop('wave_speed', 1) self._wave_speed = kwargs.pop('wave_speed', 1)
self._polynom_degree = kwargs.pop('polynom_degree', 2) self._polynomial_degree = kwargs.pop('polynomial_degree', 2)
self._cfl_number = kwargs.pop('cfl_number', 0.2) self._cfl_number = kwargs.pop('cfl_number', 0.2)
self._num_grid_cells = kwargs.pop('num_grid_cells', 64) self._num_grid_cells = kwargs.pop('num_grid_cells', 64)
self._final_time = kwargs.pop('final_time', 1) self._final_time = kwargs.pop('final_time', 1)
...@@ -79,9 +83,9 @@ class DGScheme(object): ...@@ -79,9 +83,9 @@ class DGScheme(object):
self._limiter = getattr(Limiter, self._limiter)(self._limiter_config) self._limiter = getattr(Limiter, self._limiter)(self._limiter_config)
self._quadrature = getattr(Quadrature, self._quadrature)(self._quadrature_config) self._quadrature = getattr(Quadrature, self._quadrature)(self._quadrature_config)
self._detector = getattr(Troubled_Cell_Detector, self._detector)( self._detector = getattr(Troubled_Cell_Detector, self._detector)(
self._detector_config, self._mesh, self._wave_speed, self._polynom_degree, self._num_grid_cells, self._detector_config, self._mesh, self._wave_speed, self._polynomial_degree, self._num_grid_cells,
self._final_time, self._left_bound, self._right_bound, self._basis, self._init_cond, self._quadrature) self._final_time, self._left_bound, self._right_bound, self._basis, self._init_cond, self._quadrature)
self._update_scheme = getattr(Update_Scheme, self._update_scheme)(self._polynom_degree, self._num_grid_cells, self._update_scheme = getattr(Update_Scheme, self._update_scheme)(self._polynomial_degree, self._num_grid_cells,
self._detector, self._limiter) self._detector, self._limiter)
def approximate(self): def approximate(self):
...@@ -127,7 +131,7 @@ class DGScheme(object): ...@@ -127,7 +131,7 @@ class DGScheme(object):
name = self._init_cond.get_name() + '__' + self._detector.get_name() + '__' + self._limiter.get_name() \ name = self._init_cond.get_name() + '__' + self._detector.get_name() + '__' + self._limiter.get_name() \
+ '__' + self._update_scheme.get_name() + '__' + self._quadrature.get_name() + '__final_time_' \ + '__' + self._update_scheme.get_name() + '__' + self._quadrature.get_name() + '__final_time_' \
+ str(self._final_time) + '__wave_speed_' + str(self._wave_speed) + '__number_of_cells_' \ + str(self._final_time) + '__wave_speed_' + str(self._wave_speed) + '__number_of_cells_' \
+ str(self._num_grid_cells) + '__polynomial_degree_' + str(self._polynom_degree) + str(self._num_grid_cells) + '__polynomial_degree_' + str(self._polynomial_degree)
self._detector.save_plots(name) self._detector.save_plots(name)
...@@ -135,7 +139,7 @@ class DGScheme(object): ...@@ -135,7 +139,7 @@ class DGScheme(object):
# Set additional necessary instance variables # Set additional necessary instance variables
self._interval_len = self._right_bound-self._left_bound self._interval_len = self._right_bound-self._left_bound
self._cell_len = self._interval_len / self._num_grid_cells self._cell_len = self._interval_len / self._num_grid_cells
self._basis = OrthonormalLegendre(self._polynom_degree) self._basis = OrthonormalLegendre(self._polynomial_degree)
# Set additional necessary config parameters # Set additional necessary config parameters
self._limiter_config['cell_len'] = self._cell_len self._limiter_config['cell_len'] = self._cell_len
...@@ -146,9 +150,9 @@ class DGScheme(object): ...@@ -146,9 +150,9 @@ class DGScheme(object):
# Set inverse mass matrix # Set inverse mass matrix
mass_matrix = [] mass_matrix = []
for i in range(self._polynom_degree+1): for i in range(self._polynomial_degree+1):
new_row = [] new_row = []
for j in range(self._polynom_degree+1): for j in range(self._polynomial_degree+1):
new_entry = 0.0 new_entry = 0.0
if i == j: if i == j:
new_entry = 1.0 new_entry = 1.0
...@@ -159,15 +163,15 @@ class DGScheme(object): ...@@ -159,15 +163,15 @@ class DGScheme(object):
def _do_initial_projection(self): def _do_initial_projection(self):
# Initialize matrix and set first entry to accommodate for ghost cell # Initialize matrix and set first entry to accommodate for ghost cell
output_matrix = [0] output_matrix = [0]
basis_vector = self._basis.get_vector() basis_vector = self._basis.get_basis_vector()
for cell in range(self._num_grid_cells): for cell in range(self._num_grid_cells):
new_row = [] new_row = []
eval_point = self._left_bound + (cell+0.5)*self._cell_len eval_point = self._left_bound + (cell+0.5)*self._cell_len
for degree in range(self._polynom_degree + 1): for degree in range(self._polynomial_degree + 1):
new_entry = sum(self._init_cond.calculate( new_entry = sum(
eval_point + self._cell_len/2 * self._quadrature.get_eval_points()[point]) self._init_cond.calculate(eval_point + self._cell_len/2 * self._quadrature.get_eval_points()[point])
* basis_vector[degree].subs(x, self._quadrature.get_eval_points()[point]) * basis_vector[degree].subs(x, self._quadrature.get_eval_points()[point])
* self._quadrature.get_weights()[point] * self._quadrature.get_weights()[point]
for point in range(self._quadrature.get_num_points())) for point in range(self._quadrature.get_num_points()))
......
...@@ -66,19 +66,23 @@ class FourPeakWave(InitialCondition): ...@@ -66,19 +66,23 @@ class FourPeakWave(InitialCondition):
def _get_point(self, x): def _get_point(self, x):
if (x >= -0.8) & (x <= -0.6): if (x >= -0.8) & (x <= -0.6):
return 1/6 * (self._G(x, self._z-self._delta) + self._G(x, self._z+self._delta) + 4 * self._G(x, self._z)) return 1/6 * (self._gaussian_function(x, self._z-self._delta) +
self._gaussian_function(x, self._z+self._delta) +
4 * self._gaussian_function(x, self._z))
if (x >= -0.4) & (x <= -0.2): if (x >= -0.4) & (x <= -0.2):
return 1 return 1
if (x >= 0) & (x <= 0.2): if (x >= 0) & (x <= 0.2):
return 1 - abs(10 * (x-0.1)) return 1 - abs(10 * (x-0.1))
if (x >= 0.4) & (x <= 0.6): if (x >= 0.4) & (x <= 0.6):
return 1/6 * (self._F(x, self._a-self._delta) + self._F(x, self._a+self._delta) + 4 * self._F(x, self._a)) return 1/6 * (self._elliptic_function(x, self._a-self._delta) +
self._elliptic_function(x, self._a+self._delta) +
4 * self._elliptic_function(x, self._a))
return 0 return 0
def _G(self, x, z): def _gaussian_function(self, x, z):
return np.exp(-self._beta * (x-z)**2) return np.exp(-self._beta * (x-z)**2)
def _F(self, x, a): def _elliptic_function(self, x, a):
return np.sqrt(max(1 - self._alpha**2 * (x-a)**2, 0)) return np.sqrt(max(1 - self._alpha**2 * (x-a)**2, 0))
......
...@@ -45,7 +45,7 @@ quadrature = 'Gauss' ...@@ -45,7 +45,7 @@ quadrature = 'Gauss'
dg_scheme = DGScheme(detector, detector_config=detector_config, init_cond=init_cond, init_config=init_config, dg_scheme = DGScheme(detector, detector_config=detector_config, init_cond=init_cond, init_config=init_config,
limiter=limiter, limiter_config=limiter_config, quadrature=quadrature, limiter=limiter, limiter_config=limiter_config, quadrature=quadrature,
quadrature_config=quadrature_config, update_scheme=update_scheme, wave_speed=alpha, quadrature_config=quadrature_config, update_scheme=update_scheme, wave_speed=alpha,
polynom_degree=p, cfl_number=cfl, num_grid_cells=N, final_time=finalTime, polynomial_degree=p, cfl_number=cfl, num_grid_cells=N, final_time=finalTime,
verbose=False, left_bound=xL, right_bound=xR) verbose=False, left_bound=xL, right_bound=xR)
# __, __, troubled_cells, __ = # __, __, troubled_cells, __ =
......
...@@ -2,9 +2,6 @@ ...@@ -2,9 +2,6 @@
""" """
@author: Laura C. Kühle @author: Laura C. Kühle
TODO: Check cutoff_factor/whether _is_troubled_cell works correctly \
-> Discuss! -> Seems alright; check paper (Siegfried Müller, Nils Gerhard?)
""" """
import os import os
import numpy as np import numpy as np
...@@ -16,11 +13,11 @@ xi = Symbol('z') ...@@ -16,11 +13,11 @@ xi = Symbol('z')
class TroubledCellDetector(object): class TroubledCellDetector(object):
def __init__(self, config, mesh, wave_speed, polynom_degree, num_grid_cells, final_time, left_bound, right_bound, def __init__(self, config, mesh, wave_speed, polynomial_degree, num_grid_cells, final_time, left_bound, right_bound,
basis, init_cond, quadrature): basis, init_cond, quadrature):
self._mesh = mesh self._mesh = mesh
self._wave_speed = wave_speed self._wave_speed = wave_speed
self._polynom_degree = polynom_degree self._polynomial_degree = polynomial_degree
self._num_grid_cells = num_grid_cells self._num_grid_cells = num_grid_cells
self._final_time = final_time self._final_time = final_time
self._left_bound = left_bound self._left_bound = left_bound
...@@ -55,7 +52,7 @@ class TroubledCellDetector(object): ...@@ -55,7 +52,7 @@ class TroubledCellDetector(object):
self._plot_shock_tube(troubled_cell_history, time_history) self._plot_shock_tube(troubled_cell_history, time_history)
max_error = self._plot_mesh(projection) max_error = self._plot_mesh(projection)
print("p =", self._polynom_degree) print("p =", self._polynomial_degree)
print("N =", self._num_grid_cells) print("N =", self._num_grid_cells)
print("maximum error =", max_error) print("maximum error =", max_error)
...@@ -120,8 +117,8 @@ class TroubledCellDetector(object): ...@@ -120,8 +117,8 @@ class TroubledCellDetector(object):
eval_values = [] eval_values = []
for point in range(len(eval_points)): for point in range(len(eval_points)):
new_entry = self._init_cond.calculate(eval_points[point] - self._wave_speed*self._final_time new_entry = self._init_cond.calculate(eval_points[point] - self._wave_speed*self._final_time +
+ num_periods*self._interval_len) num_periods*self._interval_len)
eval_values.append(new_entry) eval_values.append(new_entry)
grid.append(eval_points) grid.append(eval_points)
...@@ -135,13 +132,13 @@ class TroubledCellDetector(object): ...@@ -135,13 +132,13 @@ class TroubledCellDetector(object):
def _calculate_approximate_solution(self, projection): def _calculate_approximate_solution(self, projection):
points = self._quadrature.get_eval_points() points = self._quadrature.get_eval_points()
num_points = self._quadrature.get_num_points() num_points = self._quadrature.get_num_points()
basis = self._basis.get_vector() basis = self._basis.get_basis_vector()
basis_matrix = [[basis[degree].subs(x, points[point]) for point in range(num_points)] basis_matrix = [[basis[degree].subs(x, points[point]) for point in range(num_points)]
for degree in range(self._polynom_degree+1)] for degree in range(self._polynomial_degree+1)]
approx = [[sum(projection[degree][cell] * basis_matrix[degree][point] approx = [[sum(projection[degree][cell] * basis_matrix[degree][point]
for degree in range(self._polynom_degree+1)) for degree in range(self._polynomial_degree+1))
for point in range(num_points)] for point in range(num_points)]
for cell in range(len(projection[0]))] for cell in range(len(projection[0]))]
...@@ -193,7 +190,7 @@ class WaveletDetector(TroubledCellDetector): ...@@ -193,7 +190,7 @@ class WaveletDetector(TroubledCellDetector):
def _reset(self, config): def _reset(self, config):
# Set additional necessary parameter # Set additional necessary parameter
self._num_coarse_grid_cells = self._num_grid_cells//2 self._num_coarse_grid_cells = self._num_grid_cells//2
self._M1, self._M2 = self._basis.get_multiwavelet_matrices() self._wavelet_projection_left, self._wavelet_projection_right = self._basis.get_multiwavelet_projections()
def get_cells(self, projection): def get_cells(self, projection):
multiwavelet_coeffs = self._calculate_wavelet_coeffs(projection[:, 1: -1]) multiwavelet_coeffs = self._calculate_wavelet_coeffs(projection[:, 1: -1])
...@@ -202,7 +199,8 @@ class WaveletDetector(TroubledCellDetector): ...@@ -202,7 +199,8 @@ class WaveletDetector(TroubledCellDetector):
def _calculate_wavelet_coeffs(self, projection): def _calculate_wavelet_coeffs(self, projection):
output_matrix = [] output_matrix = []
for i in range(self._num_coarse_grid_cells): for i in range(self._num_coarse_grid_cells):
new_entry = 0.5*(projection[:, 2*i] @ self._M1 + projection[:, 2*i+1] @ self._M2) new_entry = 0.5*(projection[:, 2*i] @ self._wavelet_projection_left +
projection[:, 2*i+1] @ self._wavelet_projection_right)
output_matrix.append(new_entry) output_matrix.append(new_entry)
return np.transpose(np.array(output_matrix)) return np.transpose(np.array(output_matrix))
...@@ -219,7 +217,7 @@ class WaveletDetector(TroubledCellDetector): ...@@ -219,7 +217,7 @@ class WaveletDetector(TroubledCellDetector):
fine_projection = projection[:, 1:-1] fine_projection = projection[:, 1:-1]
coarse_projection = self._calculate_coarse_projection(projection) coarse_projection = self._calculate_coarse_projection(projection)
multiwavelet_coeffs = self._calculate_wavelet_coeffs(projection) multiwavelet_coeffs = self._calculate_wavelet_coeffs(projection)
basis = self._basis.get_vector() basis = self._basis.get_basis_vector()
wavelet = self._basis.get_wavelet_vector() wavelet = self._basis.get_wavelet_vector()
######################################################################################################################## ########################################################################################################################
...@@ -228,11 +226,11 @@ class WaveletDetector(TroubledCellDetector): ...@@ -228,11 +226,11 @@ class WaveletDetector(TroubledCellDetector):
# tic = timeit.default_timer() # tic = timeit.default_timer()
# averaged_projection1 = [] # averaged_projection1 = []
# wavelet_projection1 = [] # wavelet_projection1 = []
# for degree in range(self._polynom_degree + 1): # for degree in range(self._polynomial_degree + 1):
# leftMesh = coarse_projection[degree] * basis[degree].subs(x, -1 / 2) # leftMesh = coarse_projection[degree] * basis[degree].subs(x, -1 / 2)
# rightMesh = coarse_projection[degree] * basis[degree].subs(x, 1 / 2) # rightMesh = coarse_projection[degree] * basis[degree].subs(x, 1 / 2)
# leftTest = multiwavelet_coeffs[degree] * wavelet[degree].subs(xi, 1 / 2) \ # leftTest = multiwavelet_coeffs[degree] * wavelet[degree].subs(xi, 1 / 2) \
# * (-1)**(self._polynom_degree + 1 + degree) # * (-1)**(self._polynomial_degree + 1 + degree)
# rightTest = multiwavelet_coeffs[degree] * wavelet[degree].subs(xi, 1 / 2) # rightTest = multiwavelet_coeffs[degree] * wavelet[degree].subs(xi, 1 / 2)
# newRowMesh = [] # newRowMesh = []
# newRowTest = [] # newRowTest = []
...@@ -251,12 +249,12 @@ class WaveletDetector(TroubledCellDetector): ...@@ -251,12 +249,12 @@ class WaveletDetector(TroubledCellDetector):
averaged_projection = [[coarse_projection[degree][cell] * basis[degree].subs(x, value) averaged_projection = [[coarse_projection[degree][cell] * basis[degree].subs(x, value)
for cell in range(self._num_coarse_grid_cells) for cell in range(self._num_coarse_grid_cells)
for value in [-0.5, 0.5]] for value in [-0.5, 0.5]]
for degree in range(self._polynom_degree + 1)] for degree in range(self._polynomial_degree + 1)]
wavelet_projection = [[multiwavelet_coeffs[degree][cell] * wavelet[degree].subs(xi, 0.5) * value wavelet_projection = [[multiwavelet_coeffs[degree][cell] * wavelet[degree].subs(xi, 0.5) * value
for cell in range(self._num_coarse_grid_cells) for cell in range(self._num_coarse_grid_cells)
for value in [(-1) ** (self._polynom_degree + degree + 1), 1]] for value in [(-1) ** (self._polynomial_degree + degree + 1), 1]]
for degree in range(self._polynom_degree + 1)] for degree in range(self._polynomial_degree + 1)]
# toc = timeit.default_timer() # toc = timeit.default_timer()
# print("List:", toc-tic) # print("List:", toc-tic)
...@@ -265,7 +263,7 @@ class WaveletDetector(TroubledCellDetector): ...@@ -265,7 +263,7 @@ class WaveletDetector(TroubledCellDetector):
projected_coarse = np.sum(averaged_projection, axis=0) projected_coarse = np.sum(averaged_projection, axis=0)
projected_fine = np.sum([fine_projection[degree] * basis[degree].subs(x, 0) projected_fine = np.sum([fine_projection[degree] * basis[degree].subs(x, 0)
for degree in range(self._polynom_degree + 1)], axis=0) for degree in range(self._polynomial_degree + 1)], axis=0)
projected_wavelet_coeffs = np.sum(wavelet_projection, axis=0) projected_wavelet_coeffs = np.sum(wavelet_projection, axis=0)
plt.figure(4) plt.figure(4)
...@@ -277,7 +275,7 @@ class WaveletDetector(TroubledCellDetector): ...@@ -277,7 +275,7 @@ class WaveletDetector(TroubledCellDetector):
plt.title('Wavelet Coefficients') plt.title('Wavelet Coefficients')
def _calculate_coarse_projection(self, projection): def _calculate_coarse_projection(self, projection):
basis_matrix_left, basis_matrix_right = self._basis.get_basis_matrices() basis_projection_left, basis_projection_right = self._basis.get_basis_projections()
# Remove ghost cells # Remove ghost cells
projection = projection[:, 1:-1] projection = projection[:, 1:-1]
...@@ -285,8 +283,8 @@ class WaveletDetector(TroubledCellDetector): ...@@ -285,8 +283,8 @@ class WaveletDetector(TroubledCellDetector):
# Calculate projection on coarse mesh # Calculate projection on coarse mesh
output_matrix = [] output_matrix = []
for i in range(self._num_coarse_grid_cells): for i in range(self._num_coarse_grid_cells):
new_entry = 0.5 * ( new_entry = 0.5 * (projection[:, 2 * i] @ basis_projection_left +
projection[:, 2 * i] @ basis_matrix_left + projection[:, 2 * i + 1] @ basis_matrix_right) projection[:, 2 * i + 1] @ basis_projection_right)
output_matrix.append(new_entry) output_matrix.append(new_entry)
coarse_projection = np.transpose(np.array(output_matrix)) coarse_projection = np.transpose(np.array(output_matrix))
...@@ -388,7 +386,7 @@ class Theoretical(WaveletDetector): ...@@ -388,7 +386,7 @@ class Theoretical(WaveletDetector):
def _get_cells(self, multiwavelet_coeffs, projection): def _get_cells(self, multiwavelet_coeffs, projection):
troubled_cells = [] troubled_cells = []
max_avg = max(1, max(abs(projection[0][degree]) for degree in range(self._polynom_degree+1))) max_avg = max(1, max(abs(projection[0][degree]) for degree in range(self._polynomial_degree+1)))
for cell in range(self._num_coarse_grid_cells): for cell in range(self._num_coarse_grid_cells):
if self._is_troubled_cell(multiwavelet_coeffs, cell, max_avg): if self._is_troubled_cell(multiwavelet_coeffs, cell, max_avg):
...@@ -398,7 +396,7 @@ class Theoretical(WaveletDetector): ...@@ -398,7 +396,7 @@ class Theoretical(WaveletDetector):
def _is_troubled_cell(self, multiwavelet_coeffs, cell, max_avg): def _is_troubled_cell(self, multiwavelet_coeffs, cell, max_avg):
max_value = max(abs(multiwavelet_coeffs[degree][cell]) max_value = max(abs(multiwavelet_coeffs[degree][cell])
for degree in range(self._polynom_degree+1))/max_avg for degree in range(self._polynomial_degree+1))/max_avg
eps = self._cutoff_factor / (self._cell_len*self._num_coarse_grid_cells*2) eps = self._cutoff_factor / (self._cell_len*self._num_coarse_grid_cells*2)
return max_value > eps return max_value > eps
...@@ -2,14 +2,7 @@ ...@@ -2,14 +2,7 @@
""" """
@author: Laura C. Kühle @author: Laura C. Kühle
d = detail coefficient (rename?) TODO: Find better names for A, B, M1, and M2 -> Done
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
""" """
import numpy as np import numpy as np
...@@ -17,9 +10,9 @@ import timeit ...@@ -17,9 +10,9 @@ import timeit
class UpdateScheme(object): class UpdateScheme(object):
def __init__(self, polynom_degree, num_grid_cells, detector, limiter): def __init__(self, polynomial_degree, num_grid_cells, detector, limiter):
# Unpack positional arguments # Unpack positional arguments
self._polynom_degree = polynom_degree self._polynomial_degree = polynomial_degree
self._num_grid_cells = num_grid_cells self._num_grid_cells = num_grid_cells
self._detector = detector self._detector = detector
self._limiter = limiter self._limiter = limiter
...@@ -27,27 +20,27 @@ class UpdateScheme(object): ...@@ -27,27 +20,27 @@ class UpdateScheme(object):
self._reset() self._reset()
def _reset(self): def _reset(self):
# Set matrix A # Set stiffness matrix
matrix = [] matrix = []
for i in range(self._polynom_degree+1): for i in range(self._polynomial_degree+1):
new_row = [] new_row = []
for j in range(self._polynom_degree+1): for j in range(self._polynomial_degree+1):
new_entry = -1.0 new_entry = -1.0
if (j < i) & ((i+j) % 2 == 1): if (j < i) & ((i+j) % 2 == 1):
new_entry = 1.0 new_entry = 1.0
new_row.append(new_entry*np.sqrt((i+0.5) * (j+0.5))) new_row.append(new_entry*np.sqrt((i+0.5) * (j+0.5)))
matrix.append(new_row) matrix.append(new_row)
self._A = np.array(matrix) # former: inv_mass @ np.array(matrix) self._stiffness_matrix = np.array(matrix) # former: inv_mass @ np.array(matrix)
# Set matrix B # Set boundary matrix
matrix = [] matrix = []
for i in range(self._polynom_degree+1): for i in range(self._polynomial_degree+1):
new_row = [] new_row = []
for j in range(self._polynom_degree+1): for j in range(self._polynomial_degree+1):
new_entry = np.sqrt((i+0.5) * (j+0.5)) * (-1.0)**i new_entry = np.sqrt((i+0.5) * (j+0.5)) * (-1.0)**i
new_row.append(new_entry) new_row.append(new_entry)
matrix.append(new_row) matrix.append(new_row)
self._B = np.array(matrix) # former: inv_mass @ np.array(matrix) self._boundary_matrix = np.array(matrix) # former: inv_mass @ np.array(matrix)
def get_name(self): def get_name(self):
return self.__class__.__name__ return self.__class__.__name__
...@@ -111,8 +104,8 @@ class SSPRK3(UpdateScheme): ...@@ -111,8 +104,8 @@ class SSPRK3(UpdateScheme):
right_hand_side = [0] right_hand_side = [0]
for j in range(self._num_grid_cells): for j in range(self._num_grid_cells):
right_hand_side.append(2*(self._A @ current_projection[:, j+1] right_hand_side.append(2*(self._stiffness_matrix @ current_projection[:, j+1]
+ self._B @ current_projection[:, j])) + self._boundary_matrix @ current_projection[:, j]))
# Set ghost cells to respective value # Set ghost cells to respective value
right_hand_side[0] = right_hand_side[self._num_grid_cells] right_hand_side[0] = right_hand_side[self._num_grid_cells]
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment