Skip to content
Snippets Groups Projects
Select Git revision
  • 990f670ca9230e8e731c33470afcc0883d8bd506
  • master default protected
2 results

projection_utils.py

Blame
  • 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