diff --git a/Snakefile b/Snakefile
index e5ddb4e8195e4e2be82eb7ff57640517697ad9ea..fc36903fa8434af112128b923bd116addcf9689f 100644
--- a/Snakefile
+++ b/Snakefile
@@ -23,6 +23,7 @@ TODO: Contemplate allowing vector input for ICs
TODO: Discuss how wavelet details should be plotted
TODO: Ask why implicit solver is always over the interval [0,1] with 200 cells
+TODO: Ask why the zero-entry of the boundary matrix is chosen for Burgers
Urgent:
TODO: Mark Burgers' equation as inviscid -> Done
@@ -31,6 +32,8 @@ TODO: Make sure only ghost cells are limited for Dirichlet boundary -> Done
TODO: Move height adjustment and stretch factor into implicit solver
function -> Done
TODO: Enable choice of equation -> Done
+TODO: Ensure termination of Burgers with DiscontinuousConstant function -> Done
+TODO: Ensure termination of Burgers with Sine function
TODO: Rework burgers_volume_integral()
TODO: Rework/refactor implicit solver
TODO: Add subclasses of Burgers (rarefaction and shock case)
diff --git a/config.yaml b/config.yaml
index 83eceaa28d71e5c331fe604669f082f08a8fb23a..31af1a68ce03ea81f30736f38beafe7aa29c2d16 100644
--- a/config.yaml
+++ b/config.yaml
@@ -44,9 +44,10 @@ Approximation:
detector_config:
fold_len: 16
init_cond: 'FourPeakWave'
- Burgers:
+ Burgers_Test:
equation: 'Burgers'
detector: 'Theoretical'
+ init_cond: 'DiscontinuousConstant'
# Parameter for Training Data Generation
ANN_Data:
diff --git a/scripts/approximate_solution.py b/scripts/approximate_solution.py
index b7a06df2e25713050dfd659b8f1e390bd75c43ad..9c6c0081002c0541aba2969dd9cc0826b4115437 100644
--- a/scripts/approximate_solution.py
+++ b/scripts/approximate_solution.py
@@ -15,7 +15,6 @@ from tcd import Update_Scheme
from tcd.Basis_Function import OrthonormalLegendre
from tcd.DG_Approximation import DGScheme
from tcd.Mesh import Mesh
-# from tcd.Equation import LinearAdvection
def main() -> None:
@@ -35,39 +34,19 @@ def main() -> None:
# Adapt number of ghost cells
if 'stencil_len' in params['detector_config']:
- if num_ghost_cells < params['detector_config']['stencil_len']//2:
- num_ghost_cells = params['detector_config']['stencil_len']//2
+ if num_ghost_cells < \
+ params['detector_config']['stencil_len'] // 2:
+ num_ghost_cells = \
+ params['detector_config']['stencil_len'] // 2
print(f'Number of ghost cells was increased to '
- f'{num_ghost_cells} to accommodate the stencil length.')
+ f'{num_ghost_cells} to accommodate the stencil '
+ f'length.')
print(params)
- # Initialize mesh with two ghost cells on each side
- mesh = Mesh(num_cells=params.pop('num_mesh_cells', 64),
- left_bound=params.pop('left_bound', -1),
- right_bound=params.pop('right_bound', 1),
- num_ghost_cells=num_ghost_cells,
- boundary_config={'type': 'periodic'})
-
# Initialize basis
- basis = OrthonormalLegendre(params.pop('polynomial_degree', 2))
-
- # Initialize limiter after checking its legitimacy
- limiter_name = params.pop('limiter', 'ModifiedMinMod')
- limiter_config = params.pop('limiter_config', {})
- if not hasattr(Limiter, limiter_name):
- raise ValueError('Invalid limiter: "%s"' % limiter_name)
- limiter = getattr(Limiter, limiter_name)(mesh=mesh,
- config=limiter_config)
-
- # Initialize troubled cell detector after checking its legitimacy
- detector_name = params.pop('detector')
- detector_config = params.pop('detector_config', {})
- if not hasattr(Troubled_Cell_Detector, detector_name):
- raise ValueError('Invalid detector: "%s"' % detector_name)
- detector = getattr(Troubled_Cell_Detector, detector_name)(
- config=detector_config,
- mesh=mesh, basis=basis)
+ polynomial_degree = params.pop('polynomial_degree', 2)
+ basis = OrthonormalLegendre(polynomial_degree)
# Initialize quadrature after checking its legitimacy
quadrature_name = params.pop('quadrature', 'Gauss')
@@ -85,12 +64,26 @@ def main() -> None:
% init_name)
init_cond = getattr(Initial_Condition, init_name)(config=init_config)
- # Initialize equation after checking its legitimacy
wave_speed = params.pop('wave_speed', 1)
final_time = params.pop('final_time', 1)
cfl_number = params.pop('cfl_number', 0.2)
-
equation_name = params.pop('equation', 'LinearAdvection')
+
+ # Initialize mesh with two ghost cells on each side
+ boundary_config = {'type': 'periodic'} \
+ if not (equation_name == 'Burgers' and init_name ==
+ 'DiscontinuousConstant') \
+ else {'type': 'dirichlet', 'polynomial_degree': polynomial_degree,
+ 'final_time': final_time, 'num_ghost_cells': num_ghost_cells,
+ 'left_factor': init_cond._left_factor,
+ 'right_factor': init_cond._right_factor}
+ mesh = Mesh(num_cells=params.pop('num_mesh_cells', 64),
+ left_bound=params.pop('left_bound', -1),
+ right_bound=params.pop('right_bound', 1),
+ num_ghost_cells=num_ghost_cells,
+ boundary_config=boundary_config)
+
+ # Initialize equation after checking its legitimacy
if not hasattr(Equation, equation_name):
raise ValueError('Invalid equation: "%s"' % equation_name)
equation = getattr(Equation, equation_name)(
@@ -98,6 +91,23 @@ def main() -> None:
init_cond=init_cond, quadrature=quadrature, mesh=mesh,
wave_speed=wave_speed)
+ # Initialize limiter after checking its legitimacy
+ limiter_name = params.pop('limiter', 'ModifiedMinMod')
+ limiter_config = params.pop('limiter_config', {})
+ if not hasattr(Limiter, limiter_name):
+ raise ValueError('Invalid limiter: "%s"' % limiter_name)
+ limiter = getattr(Limiter, limiter_name)(mesh=mesh,
+ config=limiter_config)
+
+ # Initialize troubled cell detector after checking its legitimacy
+ detector_name = params.pop('detector')
+ detector_config = params.pop('detector_config', {})
+ if not hasattr(Troubled_Cell_Detector, detector_name):
+ raise ValueError('Invalid detector: "%s"' % detector_name)
+ detector = getattr(Troubled_Cell_Detector, detector_name)(
+ config=detector_config,
+ mesh=mesh, basis=basis)
+
# Initialize update scheme after checking its legitimacy
scheme_name = params.pop('update_scheme', 'SSPRK3')
if not hasattr(Update_Scheme, scheme_name):
diff --git a/scripts/tcd/Boundary_Condition.py b/scripts/tcd/Boundary_Condition.py
index 72f3a8cc0c4d7077c75ea78b25029f0aea4649fb..54d85d481105e9937e4043d811a243baeebb14a4 100644
--- a/scripts/tcd/Boundary_Condition.py
+++ b/scripts/tcd/Boundary_Condition.py
@@ -67,8 +67,10 @@ def dirichlet_boundary(projection: ndarray, config: dict) -> ndarray:
projection_values_left[0] = np.sqrt(2) * left_factor
projection_values_right[0] = np.sqrt(2) * right_factor
- projection[:, :num_ghost_cells] = np.repeat(projection_values_left)
- projection[:, -num_ghost_cells:] = np.repeat(projection_values_right)
+ projection[:, :num_ghost_cells] = np.repeat(
+ projection_values_left, num_ghost_cells).reshape(-1, num_ghost_cells)
+ projection[:, -num_ghost_cells:] = np.repeat(
+ projection_values_right, num_ghost_cells).reshape(-1, num_ghost_cells)
# projection[:, 0] = projection_values_left
# projection[:, 1] = projection_values_left
# projection[:, -2] = projection_values_right
diff --git a/scripts/tcd/Equation.py b/scripts/tcd/Equation.py
index a00670e05836fbd3e3937ea77978684e86ebdb18..201a03dd519a92f824ce1fd62294e075a7c72138 100644
--- a/scripts/tcd/Equation.py
+++ b/scripts/tcd/Equation.py
@@ -607,8 +607,8 @@ class Burgers(Equation):
# Set ghost cells to respective value
# (Periodic, Updated to Dirichlet in enforce_boundary_condition
# for DiscontinuousConstant problem)
- # right_hand_side[0] = right_hand_side[self._mesh.num_cells]
- # right_hand_side.append(right_hand_side[1])
+ right_hand_side[0] = right_hand_side[self._mesh.num_cells]
+ right_hand_side.append(right_hand_side[1])
return np.transpose(right_hand_side)
# right_hand_side = np.zeros_like(projection)
@@ -642,13 +642,12 @@ class Burgers(Equation):
for cell in range(self._mesh.num_cells):
new_row = []
# Approximation u at Gauss-Legendre points
- approx_eval = [sum(projection[degree][cell + 1]
- * basis_vector[degree].subs(
- x, self._quadrature.nodes[point])
- for degree in range(
- self._basis.polynomial_degree + 1))
- for point in range(self._quadrature.nodes)]
- approx_eval = np.array(approx_eval)
+ approx_eval = np.array(
+ [sum(projection[degree][cell + 1] *
+ basis_vector[degree].subs(x,
+ self._quadrature.nodes[point])
+ for degree in range(self._basis.polynomial_degree + 1))
+ for point in range(self._quadrature.num_nodes)])
# print("u", approx_eval)
# Quadrature Evaluation
@@ -658,7 +657,7 @@ class Burgers(Equation):
x, self._quadrature.nodes[point])
* self._quadrature.weights[point]
for point in range(self._quadrature.num_nodes)]
- new_entry = sum(list(approx_eval**2 * weighted_derivatives)) / 2
+ new_entry = sum(list(approx_eval**2 * weighted_derivatives))/2
new_row.append(np.float64(new_entry))
@@ -668,9 +667,9 @@ class Burgers(Equation):
# Set ghost cells to respective value
# (Periodic, Updated to Dirichlet in enforce_boundary_condition
# for DiscontinuousConstant problem)
- # volume_integral_matrix[0] = volume_integral_matrix[
- # self._mesh.num_cells]
- # volume_integral_matrix.append(volume_integral_matrix[1])
+ volume_integral_matrix[0] = volume_integral_matrix[
+ self._mesh.num_cells]
+ volume_integral_matrix.append(volume_integral_matrix[1])
return np.transpose(np.array(volume_integral_matrix))
@@ -726,8 +725,8 @@ class Burgers(Equation):
# Set Ghost cells to respective value
# (Periodic, Updated to Dirichlet in enforce_boundary_condition
# for DiscontinuousConstant problem)
- # boundary_matrix[0] = boundary_matrix[self._mesh.num_cells]
- # boundary_matrix.append(boundary_matrix[1])
+ boundary_matrix[0] = boundary_matrix[self._mesh.num_cells]
+ boundary_matrix.append(boundary_matrix[1])
boundary_matrix = np.transpose(np.array(boundary_matrix))