From 042719cfbe01ac9fe647b70539f24b835231b321 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?K=C3=BChle=2C=20Laura=20Christine=20=28lakue103=29?=
<laura.kuehle@uni-duesseldorf.de>
Date: Thu, 2 Jun 2022 18:42:05 +0200
Subject: [PATCH] Contained interval boundaries in mesh.
---
DG_Approximation.py | 48 ++++++++++++++++---------------------
Initial_Condition.py | 30 +++++++++--------------
Plotting.py | 8 ++-----
Troubled_Cell_Detector.py | 17 ++++---------
projection_utils.py | 2 +-
workflows/ANN_data.smk | 3 +--
workflows/approximation.smk | 6 ++---
7 files changed, 42 insertions(+), 72 deletions(-)
diff --git a/DG_Approximation.py b/DG_Approximation.py
index 938de93..fd33f4c 100644
--- a/DG_Approximation.py
+++ b/DG_Approximation.py
@@ -14,9 +14,9 @@ TODO: Contemplate extracting boundary condition from InitialCondition
Urgent:
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 bounds in mesh -> Done
TODO: Contain interval length in mesh
+TODO: Contain number of grid cells in mesh
TODO: Create data dict for mesh separately
TODO: Create data dict for basis separately
TODO: Refactor eval_points in do_initial_projection()
@@ -167,8 +167,6 @@ class DGScheme:
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)
- self._left_bound = kwargs.pop('left_bound', -1)
- self._right_bound = kwargs.pop('right_bound', 1)
self._verbose = kwargs.pop('verbose', False)
self._history_threshold = kwargs.pop('history_threshold',
math.ceil(0.2/self._cfl_number))
@@ -185,6 +183,18 @@ class DGScheme:
polynomial_degree = kwargs.pop('polynomial_degree', 2)
self._basis = OrthonormalLegendre(polynomial_degree)
+ # Initialize mesh with two ghost cells on each side
+ left_bound = kwargs.pop('left_bound', -1)
+ right_bound = kwargs.pop('right_bound', 1)
+ self._mesh = Mesh(num_grid_cells=self._num_grid_cells,
+ num_ghost_cells=2, left_bound=left_bound,
+ right_bound=right_bound)
+ print(len(self._mesh.cells))
+ print(type(self._mesh.cells))
+
+ # Set additional necessary instance variables
+ self._interval_len = right_bound-left_bound
+
# Throw an error if there are extra keyword arguments
if len(kwargs) > 0:
extra = ', '.join('"%s"' % k for k in list(kwargs.keys()))
@@ -209,7 +219,6 @@ class DGScheme:
# Replace the string names with the actual class instances
# (and add the instance variables for the quadrature)
self._init_cond = getattr(Initial_Condition, self._init_cond)(
- left_bound=self._left_bound, right_bound=self._right_bound,
config=self._init_config)
self._limiter = getattr(Limiter, self._limiter)(
config=self._limiter_config)
@@ -218,8 +227,7 @@ 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,
- final_time=self._final_time, left_bound=self._left_bound,
- right_bound=self._right_bound, basis=self._basis,
+ final_time=self._final_time, basis=self._basis,
init_cond=self._init_cond, quadrature=self._quadrature)
self._update_scheme = getattr(Update_Scheme, self._update_scheme)(
polynomial_degree=polynomial_degree,
@@ -244,8 +252,7 @@ class DGScheme:
projection = do_initial_projection(
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)
+ num_grid_cells=self._num_grid_cells)
time_step = abs(self._cfl_number * self._mesh.cell_len /
self._wave_speed)
@@ -290,17 +297,6 @@ class DGScheme:
def _reset(self):
"""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
-
- # Initialize mesh with two ghost cells on each side
- self._mesh = Mesh(num_grid_cells=self._num_grid_cells,
- num_ghost_cells=2, left_bound=self._left_bound,
- right_bound=self._right_bound)
- print(len(self._mesh.cells))
- print(type(self._mesh.cells))
-
# Set additional necessary config parameters
self._limiter_config['cell_len'] = self._mesh.cell_len
@@ -335,8 +331,7 @@ class DGScheme:
projection = do_initial_projection(
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)
+ num_grid_cells=self._num_grid_cells, adjustment=adjustment)
return self._basis.calculate_cell_average(
projection=projection[:, 1:-1], stencil_length=stencil_length,
@@ -344,7 +339,7 @@ class DGScheme:
def do_initial_projection(initial_condition, mesh, basis, quadrature,
- num_grid_cells, left_bound, adjustment=0):
+ num_grid_cells, adjustment=0):
"""Calculates initial projection.
Calculates a projection at time step 0 and adds ghost cells on both
@@ -362,8 +357,6 @@ def do_initial_projection(initial_condition, mesh, basis, quadrature,
Quadrature used for evaluation.
num_grid_cells : int
Number of cells in the mesh. Usually exponential of 2.
- left_bound : float
- Left boundary of interval.
adjustment: float, optional
Extent of adjustment of each evaluation point in x-direction.
Default: 0.
@@ -379,14 +372,13 @@ def do_initial_projection(initial_condition, mesh, basis, quadrature,
output_matrix = [0]
basis_vector = basis.basis
- # 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)*mesh.cell_len
+ eval_point = mesh.bounds[0] + (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 + mesh.cell_len/2
+ mesh, 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/Initial_Condition.py b/Initial_Condition.py
index 1dd6420..a5a3f1a 100644
--- a/Initial_Condition.py
+++ b/Initial_Condition.py
@@ -10,11 +10,6 @@ import numpy as np
class InitialCondition(ABC):
"""Abstract class for initial condition function.
- Attributes
- ----------
- interval_len : float
- Length of the interval between left and right boundary.
-
Methods
-------
get_name()
@@ -29,22 +24,15 @@ class InitialCondition(ABC):
Evaluates function at given x-value.
"""
- def __init__(self, left_bound, right_bound, config):
+ def __init__(self, config):
"""Initializes InitialCondition.
Parameters
----------
- left_bound : float
- Left boundary of interval.
- right_bound : float
- Right boundary of interval.
config : dict
Additional parameters for initial condition.
"""
- self._left_bound = left_bound
- self._right_bound = right_bound
-
self._reset(config)
def _reset(self, config):
@@ -56,7 +44,7 @@ class InitialCondition(ABC):
Additional parameters for initial condition.
"""
- self._interval_len = self._right_bound-self._left_bound
+ pass
def get_name(self):
"""Returns string of class name."""
@@ -89,7 +77,7 @@ class InitialCondition(ABC):
"""
pass
- def calculate(self, x):
+ def calculate(self, mesh, x):
"""Evaluates function at given x-value.
Projects x-value into interval of the periodic function and evaluates
@@ -97,6 +85,8 @@ class InitialCondition(ABC):
Parameters
----------
+ mesh : Mesh
+ Mesh for calculation.
x : float
Evaluation point of function.
@@ -106,10 +96,12 @@ class InitialCondition(ABC):
Value of function evaluates at x-value.
"""
- while x < self._left_bound:
- x += self._interval_len
- while x > self._right_bound:
- x -= self._interval_len
+ left_bound, right_bound = mesh.bounds
+ interval_len = right_bound-left_bound
+ while x < left_bound:
+ x += interval_len
+ while x > right_bound:
+ x -= interval_len
return self._get_point(x)
@abstractmethod
diff --git a/Plotting.py b/Plotting.py
index f58d547..66b932e 100644
--- a/Plotting.py
+++ b/Plotting.py
@@ -370,8 +370,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: Mesh, num_grid_cells: int,
- wave_speed: float, final_time: float,
- left_bound: float, right_bound: float, basis: Basis,
+ wave_speed: float, final_time: float, basis: Basis,
quadrature: Quadrature, init_cond: InitialCondition,
colors: dict = None, coarse_projection: ndarray = None,
multiwavelet_coeffs: ndarray = None) -> None:
@@ -399,10 +398,6 @@ def plot_results(projection: ndarray, troubled_cell_history: list,
Speed of wave in rightward direction.
final_time : float
Final time for which approximation is calculated.
- left_bound : float
- Left boundary of interval.
- right_bound : float
- Right boundary of interval.
basis: Vector object
Basis used for calculation.
quadrature: Quadrature object
@@ -424,6 +419,7 @@ def plot_results(projection: ndarray, troubled_cell_history: list,
colors = _check_colors(colors)
# Calculate needed variables
+ left_bound, right_bound = mesh.bounds
interval_len = right_bound-left_bound
# Plot troubled cells
diff --git a/Troubled_Cell_Detector.py b/Troubled_Cell_Detector.py
index 209dee4..e1ac2a0 100644
--- a/Troubled_Cell_Detector.py
+++ b/Troubled_Cell_Detector.py
@@ -35,8 +35,7 @@ class TroubledCellDetector(ABC):
"""
def __init__(self, config, init_cond, quadrature, basis, mesh,
- wave_speed=1, num_grid_cells=64,
- final_time=1, left_bound=-1, right_bound=1):
+ wave_speed=1, num_grid_cells=64, final_time=1):
"""Initializes TroubledCellDetector.
Parameters
@@ -57,18 +56,13 @@ class TroubledCellDetector(ABC):
Number of cells in the mesh. Usually exponential of 2. Default: 64.
final_time : float, optional
Final time for which approximation is calculated. Default: 1.
- left_bound : float, optional
- Left boundary of interval. Default: -1.
- right_bound : float, optional
- Right boundary of interval. Default: 1.
"""
self._mesh = mesh
self._wave_speed = wave_speed
self._num_grid_cells = num_grid_cells
self._final_time = final_time
- self._left_bound = left_bound
- self._right_bound = right_bound
+ left_bound, right_bound = mesh.bounds
self._interval_len = right_bound - left_bound
self._basis = basis
self._init_cond = init_cond
@@ -104,15 +98,14 @@ class TroubledCellDetector(ABC):
pass
def create_data_dict(self, projection):
+ left_bound, right_bound = self._mesh.bounds
return {'projection': projection, 'wave_speed': self._wave_speed,
'num_grid_cells': self._num_grid_cells,
'final_time': self._final_time,
- 'left_bound': self._left_bound,
- 'right_bound': self._right_bound,
'basis': {'polynomial_degree': self._basis.polynomial_degree},
'mesh': {'num_grid_cells': self._num_grid_cells,
- 'left_bound': self._left_bound,
- 'right_bound': self._right_bound,
+ 'left_bound': left_bound,
+ 'right_bound': right_bound,
'num_ghost_cells': 2}
}
diff --git a/projection_utils.py b/projection_utils.py
index b4df91f..0b6f0b6 100644
--- a/projection_utils.py
+++ b/projection_utils.py
@@ -169,7 +169,7 @@ def calculate_exact_solution(
eval_values = []
for eval_point in eval_points:
- new_entry = init_cond.calculate(eval_point
+ new_entry = init_cond.calculate(mesh, eval_point
- wave_speed * final_time
+ num_periods * interval_len)
eval_values.append(new_entry)
diff --git a/workflows/ANN_data.smk b/workflows/ANN_data.smk
index c0af4d0..df966e1 100644
--- a/workflows/ANN_data.smk
+++ b/workflows/ANN_data.smk
@@ -30,8 +30,7 @@ rule generate_data:
initial_conditions = []
for function in params.functions:
initial_conditions.append({
- 'function': getattr(Initial_Condition, function)(
- params.left_bound, params.right_bound, {}),
+ 'function': getattr(Initial_Condition, function)({}),
'config': {} if function in config['functions']
else config['functions'][function]
})
diff --git a/workflows/approximation.smk b/workflows/approximation.smk
index 00b89a3..cee3af3 100644
--- a/workflows/approximation.smk
+++ b/workflows/approximation.smk
@@ -6,6 +6,7 @@ import Initial_Condition
import Quadrature
from DG_Approximation import DGScheme
from Plotting import plot_approximation_results
+from projection_utils import Mesh
configfile: 'config.yaml'
@@ -82,11 +83,8 @@ rule plot_approximation_results:
detector_dict = params.dg_params.copy()
- left_bound = detector_dict.pop('left_bound', -1)
- right_bound = detector_dict.pop('right_bound', 1)
init_cond = getattr(Initial_Condition, detector_dict.pop(
- 'init_cond', 'Sine'))(left_bound=left_bound,
- right_bound=right_bound,
+ 'init_cond', 'Sine'))(
config=detector_dict.pop('init_config', {}))
quadrature = getattr(Quadrature, detector_dict.pop(
'quadrature', 'Gauss'))(detector_dict.pop(
--
GitLab