Select Git revision
projection_utils.py
Laura Christine Kühle authored
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
projection_utils.py 6.61 KiB
# -*- coding: utf-8 -*-
"""
@author: Laura C. Kühle
"""
from __future__ import annotations
from functools import cache
from typing import Tuple
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/grid.
Each cell is characterized by its center.
Attributes
----------
num_grid_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_grid_cells: int, num_ghost_cells: int,
left_bound: float, right_bound: float) -> None:
"""Initialize Mesh.
Parameters
----------
num_grid_cells : int
Number of cells in the mesh (ghost cells notwithstanding). Usually
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.
"""
self._num_grid_cells = num_grid_cells
self._num_ghost_cells = num_ghost_cells
self._left_bound = left_bound
self._right_bound = right_bound
@property
def num_grid_cells(self) -> int:
"""Return number of grid cells."""
return self._num_grid_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_grid_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_grid_cells': self._num_grid_cells,
'left_bound': self._left_bound,
'right_bound': self._right_bound,
'num_ghost_cells': self._num_ghost_cells}
def random_stencil(self, stencil_length: 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 grid spacing to be within interval if necessary
grid_spacing = self.cell_len
max_spacing = 2/stencil_length*min(point-self._left_bound,
self._right_bound-point)
while grid_spacing > max_spacing:
grid_spacing /= 2
# Return new mesh instance
return Mesh(left_bound=point - stencil_length/2 * grid_spacing,
right_bound=point + stencil_length/2 * grid_spacing,
num_grid_cells=stencil_length, num_ghost_cells=2)
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