diff --git a/Snakefile b/Snakefile
index 2fbf55308db60ce38218907736a7713f869c34c6..14b18dd9705eb9aa714297ddbf60813702d79eee 100644
--- a/Snakefile
+++ b/Snakefile
@@ -21,9 +21,16 @@ TODO: Discuss descriptions (matrices, cfl number, right-hand side,
TODO: Discuss referencing info on SSPRK3
TODO: Discuss name for quadrature mesh (now: grid)
TODO: Contemplate using lambdify for basis
+TODO: Ask why MinMod slope is only calculated from degree 1, not 0
Urgent:
-TODO: Vectorize '_calculate_reconstructions()' in OrthonormalLegendre
+TODO: Vectorize '_calculate_reconstructions()' in OrthonormalLegendre -> Done
+TODO: Vectorize 'do_initial_projection()' -> Done
+TODO: Rework UpdateScheme to limit over array instead of cell -> Done
+TODO: Vectorize '_determine_modification()' in ModifiedMinMod
+TODO: Vectorize '_set_cell_slope()' in MinMod
+TODO: Vectorize '_determine_modification()' in MinMod
+TODO: Vectorize 'apply()' in MinMod
TODO: Replace loops/list comprehension with vectorization if feasible
TODO: Replace loops with list comprehension if feasible
TODO: Rework ICs to allow vector input
diff --git a/scripts/tcd/Limiter.py b/scripts/tcd/Limiter.py
index 8b797896e05c205e199c09c3dac0f0d51ad40fc1..bc3b6c04e208ebfe9582c05c1b72713f7b142609 100644
--- a/scripts/tcd/Limiter.py
+++ b/scripts/tcd/Limiter.py
@@ -4,6 +4,7 @@
"""
from abc import ABC, abstractmethod
+import numpy as np
class Limiter(ABC):
@@ -13,8 +14,8 @@ class Limiter(ABC):
-------
get_name()
Returns string of class name.
- apply(projection, cell)
- Applies limiting to cell.
+ apply(projection, cells)
+ Applies limiting to cells.
"""
def __init__(self, config):
@@ -45,15 +46,15 @@ class Limiter(ABC):
return self.__class__.__name__
@abstractmethod
- def apply(self, projection, cell):
- """Applies limiting to cell.
+ def apply(self, projection, cells):
+ """Applies limiting to cells.
Parameters
----------
projection : ndarray
Matrix of projection for each polynomial degree.
- cell: int
- Index of cell.
+ cells : list
+ Index of cells to limit.
Returns
-------
@@ -71,13 +72,13 @@ class NoLimiter(Limiter):
-------
get_name()
Returns string of class name.
- apply(projection, cell)
- Applies no limiting to cell.
+ apply(projection, cells)
+ Applies no limiting to cells.
"""
- def apply(self, projection, cell):
- """Returns projection of cell without limiting."""
- return projection[:, cell]
+ def apply(self, projection, cells):
+ """Returns projection without limiting."""
+ return projection
class MinMod(Limiter):
@@ -96,8 +97,8 @@ class MinMod(Limiter):
-------
get_name()
Returns string of class name.
- apply(projection, cell)
- Applies limiting to cell.
+ apply(projection, cells)
+ Applies limiting to cells.
"""
def _reset(self, config):
@@ -116,34 +117,35 @@ class MinMod(Limiter):
"""Returns string of class name concatenated with the erase-degree."""
return self.__class__.__name__ + str(self._erase_degree)
- def apply(self, projection, cell):
- """Applies limiting to cell.
+ def apply(self, projection, cells):
+ """Applies limiting to cells.
Parameters
----------
projection : ndarray
Matrix of projection for each polynomial degree.
- cell : int
- Index of cell.
+ cells : list
+ Index of cells to limit.
Returns
-------
- adapted_projection : ndarray
+ new_projection : ndarray
Matrix of updated projection for each polynomial degree.
"""
- cell_slope = self._set_cell_slope(projection, cell)
- is_good_cell = self._determine_modification(projection, cell,
- cell_slope)
-
- if is_good_cell:
- return projection[:, cell]
-
- adapted_projection = projection[:, cell].copy()
- for i in range(len(adapted_projection)):
- if i > self._erase_degree:
- adapted_projection[i] = 0
- return adapted_projection
+ new_projection = projection.copy()
+ for cell in cells:
+ cell_slope = self._set_cell_slope(projection, cell)
+ is_good_cell = self._determine_modification(projection, cell,
+ cell_slope)
+
+ if not is_good_cell:
+ adapted_projection = projection[:, cell].copy()
+ for i in range(len(adapted_projection)):
+ if i > self._erase_degree:
+ adapted_projection[i] = 0
+ new_projection[:, cell] = adapted_projection
+ return new_projection
def _determine_modification(self, projection, cell, cell_slope):
"""Determines whether limiting is applied.
@@ -197,6 +199,16 @@ class MinMod(Limiter):
slope.append(new_entry)
return slope[cell]
+ # # print(np.array(slope).shape)
+ # # print(slope)
+ # root_vector = np.array([np.sqrt(degree+0.5)
+ # for degree in range(len(projection))])
+ # test = root_vector[1:] @ projection[1:]
+ # # print(test.shape)
+ # # print(np.isclose(test, slope, rtol=1e-16))
+ #
+ # return test[cell]
+
class ModifiedMinMod(MinMod):
"""Class for modified minmod limiting function.
@@ -235,6 +247,9 @@ class ModifiedMinMod(MinMod):
super()._reset(config)
# Unpack necessary configurations
+ # cell_len = config.pop('cell_len')
+ # mod_factor = config.pop('mod_factor', 0)
+ # self._threshold = mod_factor * cell_len**2
self._cell_len = config.pop('cell_len')
self._mod_factor = config.pop('mod_factor', 0)
@@ -260,6 +275,7 @@ class ModifiedMinMod(MinMod):
Flag whether cell should be adjusted.
"""
+ # if abs(cell_slope) <= self._threshold:
if abs(cell_slope) <= self._mod_factor*self._cell_len**2:
return True
diff --git a/scripts/tcd/Update_Scheme.py b/scripts/tcd/Update_Scheme.py
index d789576e53c8e7463f0a3e2dfe92557ed83b8677..bc740d9bd3c62a456787278014284b817d674484 100644
--- a/scripts/tcd/Update_Scheme.py
+++ b/scripts/tcd/Update_Scheme.py
@@ -154,11 +154,8 @@ class UpdateScheme(ABC):
"""
troubled_cells = self._detector.get_cells(current_projection)
-
- new_projection = current_projection.copy()
- for cell in troubled_cells:
- new_projection[:, cell] = self._limiter.apply(current_projection,
- cell)
+ new_projection = self._limiter.apply(current_projection,
+ troubled_cells)
return new_projection, troubled_cells