Select Git revision
configuration.py
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
projection_utils.py 7.29 KiB
# -*- coding: utf-8 -*-
"""
@author: Laura C. Kühle
"""
from __future__ import annotations
from functools import cache
from typing import Tuple
import math
import numpy as np
from numpy import ndarray
from sympy import Symbol
from Quadrature import Quadrature
from Initial_Condition import InitialCondition
x = Symbol('x')
class Mesh:
"""Class for mesh.
Each cell is characterized by its center.
Attributes
----------
mode : str
Mode for mesh use. Either 'training' or 'evaluation'.
num_cells : int
Number of cells in the mesh (ghost cells notwithstanding). Usually
exponential of 2.
bounds : Tuple[float, float]
Left and right boundary of the mesh interval.
interval_len : float
Length of the mesh interval.
cell_len : float
Length of a mesh cell.
cells : ndarray
Array of cell centers in mesh.
"""
def __init__(self, num_cells: int, num_ghost_cells: int,
left_bound: float, right_bound: float,
training_data_mode: bool = False) -> None:
"""Initialize Mesh.
Parameters
----------
num_cells : int
Number of cells in the mesh (ghost cells notwithstanding). Has
to be an exponential of 2.
num_ghost_cells : int
Number of ghost cells on each side of the mesh.
left_bound : float
Left boundary of the mesh interval.
right_bound : float
Right boundary of the mesh interval.
training_data_mode : bool, optional
Flag indicating whether the mesh is used for training data
generation. Default: False.
"""
self._num_cells = num_cells
self._mode = 'training' if training_data_mode else 'evaluation'
if not training_data_mode:
if not math.log(self._num_cells, 2).is_integer():
raise ValueError('The number of cells in the mesh has to be '
'an exponential of 2')
self._num_ghost_cells = num_ghost_cells
self._left_bound = left_bound
self._right_bound = right_bound
@property
def mode(self) -> str:
"""Return mode ('training' or 'evaluation')."""
return self._mode
@property
def num_cells(self) -> int:
"""Return number of mesh cells."""
return self._num_cells
@property
def bounds(self) -> Tuple[float, float]:
"""Return left and right boundary of the mesh interval."""
return self._left_bound, self._right_bound
@property
def interval_len(self) -> float:
"""Return the length of the mesh interval."""
return self._right_bound - self._left_bound
@property
def cell_len(self) -> float:
"""Return the length of a mesh cell."""
return self.interval_len/self.num_cells
@property
@cache
def cells(self) -> ndarray:
"""Return the cell centers of the mesh (including ghost cells)."""
return np.arange(
self._left_bound - (self._num_ghost_cells*2-1)/2*self.cell_len,
self._right_bound + (self._num_ghost_cells*2+1)/2*self.cell_len,
self.cell_len)
@property
def non_ghost_cells(self) -> ndarray:
"""Return the cell centers of the mesh (excluding ghost cells)."""
return self.cells[self._num_ghost_cells:-self._num_ghost_cells]
def create_data_dict(self) -> dict:
"""Return dictionary with data necessary to construct mesh."""
return {'num_cells': self._num_cells,
'left_bound': self._left_bound,
'right_bound': self._right_bound,
'num_ghost_cells': self._num_ghost_cells}
def random_stencil(self, stencil_len: int) -> Mesh:
"""Return random stencil.
Build mesh with given number of cell centers around a random point
in the underlying interval.
Returns
-------
Mesh object
Mesh of given size around random point.
"""
# Pick random point between left and right bound
point = np.random.uniform(self._left_bound, self._right_bound)
# Adjust mesh spacing to be within interval if necessary
mesh_spacing = self.cell_len
max_spacing = 2/stencil_len*min(point-self._left_bound,
self._right_bound-point)
while mesh_spacing > max_spacing:
mesh_spacing /= 2
# Return new mesh instance
return Mesh(left_bound=point - stencil_len/2 * mesh_spacing,
right_bound=point + stencil_len/2 * mesh_spacing,
num_cells=stencil_len, num_ghost_cells=2,
training_data_mode=True)
def calculate_approximate_solution(
projection: ndarray, points: ndarray, polynomial_degree: int,
basis: ndarray) -> ndarray:
"""Calculate approximate solution.
Parameters
----------
projection : ndarray
Matrix of projection for each polynomial degree.
points : ndarray
List of evaluation points for mesh.
polynomial_degree : int
Polynomial degree.
basis : ndarray
Basis vector for calculation.
Returns
-------
ndarray
Array containing approximate evaluation of a function.
"""
num_points = len(points)
basis_matrix = [[basis[degree].subs(x, points[point])
for point in range(num_points)]
for degree in range(polynomial_degree+1)]
approx = [[sum(projection[degree][cell] * basis_matrix[degree][point]
for degree in range(polynomial_degree+1))
for point in range(num_points)]
for cell in range(len(projection[0]))]
return np.reshape(np.array(approx), (1, len(approx) * num_points))
def calculate_exact_solution(
mesh: Mesh, wave_speed: float, final_time:
float, quadrature: Quadrature,
init_cond: InitialCondition) -> Tuple[ndarray, ndarray]:
"""Calculate exact solution.
Parameters
----------
mesh : Mesh
Mesh for evaluation.
wave_speed : float
Speed of wave in rightward direction.
final_time : float
Final time for which approximation is calculated.
quadrature : Quadrature object
Quadrature for evaluation.
init_cond : InitialCondition object
Initial condition for evaluation.
Returns
-------
grid : ndarray
Array containing evaluation grid for a function.
exact : ndarray
Array containing exact evaluation of a function.
"""
grid = []
exact = []
num_periods = np.floor(wave_speed * final_time / mesh.interval_len)
for cell_center in mesh.non_ghost_cells:
eval_points = cell_center+mesh.cell_len / 2 * \
quadrature.nodes
eval_values = []
for eval_point in eval_points:
new_entry = init_cond.calculate(mesh=mesh, x=eval_point
- wave_speed * final_time
+ num_periods * mesh.interval_len)
eval_values.append(new_entry)
grid.append(eval_points)
exact.append(eval_values)
exact = np.reshape(np.array(exact), (1, len(exact) * len(exact[0])))
grid = np.reshape(np.array(grid), (1, len(grid) * len(grid[0])))
return grid, exact