Skip to content
Snippets Groups Projects
Select Git revision
  • 141632da14ea2a2fc0b388927b8424f141abb04d
  • master default protected
2 results

DG_Approximation.py

Blame
  • Code owners
    Assign users and groups as approvers for specific file changes. Learn more.
    DG_Approximation.py 15.14 KiB
    # -*- coding: utf-8 -*-
    """
    @author: Laura C. Kühle
    
    Urgent:
    TODO: Move plotting into separate function
    TODO: Adapt TCD from Soraya
        (Dropbox->...->TEST_troubled-cell-detector->Troubled_Cell_Detector)
    TODO: Add verbose output
    TODO: Improve file naming (e.g. use '.' instead of '__')
    
    Critical, but not urgent:
    TODO: Use cfl_number for updating, not just time
    TODO: Adjust code to allow classes for all equations
        (Burger, linear advection, 1D Euler)
    
    Currently not critical:
    TODO: Replace loops with list comprehension if feasible
    TODO: Check whether 'projection' is always a np.array()
    TODO: Check whether all instance variables are sensible
    TODO: Rename files according to standard
    TODO: Outsource scripts into separate directory
    TODO: Allow comparison between ANN training datasets
    TODO: Add a default model state
    TODO: Add an environment file for Snakemake
    
    Not feasible yet or doc-related:
    TODO: Double-check everything! (also with pylint, pytype, pydoc, pycodestyle)
    TODO: Check whether documentation style is correct
    TODO: Check whether all types in doc are correct
    TODO: Discuss adding kwargs to attributes in documentation
    TODO: Add type annotations to function heads
    
    """
    import os
    import json
    import numpy as np
    from sympy import Symbol
    import math
    import matplotlib
    from matplotlib import pyplot as plt
    
    import Troubled_Cell_Detector
    import Initial_Condition
    import Limiter
    import Quadrature
    import Update_Scheme
    from Basis_Function import OrthonormalLegendre
    
    matplotlib.use('Agg')
    x = Symbol('x')
    
    
    def encode_ndarray(obj):
        if isinstance(obj, np.ndarray):
            return obj.tolist()
        return obj
    
    
    def decode_ndarray(obj):
        if isinstance(obj, list):
            return np.asarray(obj)
        return obj
    
    
    class DGScheme:
        """Class for Discontinuous Galerkin Method.
    
        Approximates linear advection equation using Discontinuous Galerkin Method
        with troubled-cell-based limiting.
    
        Attributes
        ----------
        interval_len : float
            Length of the interval between left and right boundary.
        cell_len : float
            Length of a cell in mesh.
        basis : Basis object
            Basis for calculation.
        mesh : ndarray
            List of mesh valuation points.
        inv_mass : ndarray
            Inverse mass matrix.
    
        Methods
        -------
        approximate()
            Approximates projection.
        save_plots()
            Saves plots generated during approximation process.
        build_training_data(adjustment, stencil_length, initial_condition=None)
            Builds training data set.
    
        """
        def __init__(self, detector, **kwargs):
            """Initializes DGScheme.
    
            Parameters
            ----------
            detector : str
                Name of troubled cell detector class.
    
            Other Parameters
            ----------------
            wave_speed : float, optional
                Speed of wave in rightward direction. Default: 1.
            polynomial_degree : int, optional
                Polynomial degree. Default: 2.
            cfl_number : float, optional
                CFL number to ensure stability. Default: 0.2.
            num_grid_cells : int, optional
                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.
            verbose : bool, optional
                Flag whether commentary in console is wanted. Default: False.
            plot_dir : str, optional
                Path to directory in which plots are saved. Default: 'test'.
            history_threshold : float, optional
                Threshold when history will be recorded.
                Default: math.ceil(0.2/cfl_number).
            detector_config : dict, optional
                Additional parameters for detector object. Default: {}.
            init_cond : str, optional
                Name of initial condition for evaluation. Default: 'Sine'
            init_config : dict, optional
                Additional parameters for initial condition object. Default: {}.
            limiter : str, optional
                Name of limiter for evaluation. Default: 'ModifiedMinMod'.
            limiter_config : dict, optional
                Additional parameters for limiter. object. Default: {}:
            quadrature : str, optional
                Name of quadrature for evaluation. Default: 'Gauss'.
            quadrature_config : dict, optional
                Additional parameters for quadrature object. Default: {}.
            update_scheme : str, optional
                Name of update scheme for evaluation. Default: 'SSPRK3'.
    
            """
            # Unpack keyword arguments
            self._wave_speed = kwargs.pop('wave_speed', 1)
            self._polynomial_degree = kwargs.pop('polynomial_degree', 2)
            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._plot_dir = kwargs.pop('plot_dir', 'testing')
            self._history_threshold = kwargs.pop('history_threshold',
                                                 math.ceil(0.2/self._cfl_number))
            self._detector = detector
            self._detector_config = kwargs.pop('detector_config', {})
            self._init_cond = kwargs.pop('init_cond', 'Sine')
            self._init_config = kwargs.pop('init_config', {})
            self._limiter = kwargs.pop('limiter', 'ModifiedMinMod')
            self._limiter_config = kwargs.pop('limiter_config', {})
            self._quadrature = kwargs.pop('quadrature', 'Gauss')
            self._quadrature_config = kwargs.pop('quadrature_config', {})
            self._update_scheme = kwargs.pop('update_scheme', 'SSPRK3')
    
            # Throw an error if there are extra keyword arguments
            if len(kwargs) > 0:
                extra = ', '.join('"%s"' % k for k in list(kwargs.keys()))
                raise ValueError('Unrecognized arguments: %s' % extra)
    
            # Make sure all classes actually exist
            if not hasattr(Troubled_Cell_Detector, self._detector):
                raise ValueError('Invalid detector: "%s"' % self._detector)
            if not hasattr(Initial_Condition, self._init_cond):
                raise ValueError('Invalid initial condition: "%s"'
                                 % self._init_cond)
            if not hasattr(Limiter, self._limiter):
                raise ValueError('Invalid limiter: "%s"' % self._limiter)
            if not hasattr(Quadrature, self._quadrature):
                raise ValueError('Invalid quadrature: "%s"' % self._quadrature)
            if not hasattr(Update_Scheme, self._update_scheme):
                raise ValueError('Invalid update scheme: "%s"'
                                 % self._update_scheme)
    
            self._reset()
    
            # 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)(
                self._left_bound, self._right_bound, self._init_config)
            self._limiter = getattr(Limiter, self._limiter)(self._limiter_config)
            self._quadrature = getattr(Quadrature, self._quadrature)(
                self._quadrature_config)
            self._detector = getattr(Troubled_Cell_Detector, self._detector)(
                self._detector_config, self._mesh, self._wave_speed,
                self._polynomial_degree, self._num_grid_cells, self._final_time,
                self._left_bound, self._right_bound, self._basis,
                self._init_cond, self._quadrature)
            self._update_scheme = getattr(Update_Scheme, self._update_scheme)(
                self._polynomial_degree, self._num_grid_cells, self._detector,
                self._limiter)
    
        def approximate(self, data_name):
            """Approximates projection.
    
            Initializes projection and evolves it in time. Each time step consists
            of three parts: A projection update, a troubled-cell detection,
            and limiting based on the detected cells.
    
            At final time, result and error plots are
            generated and, if verbose flag is set, also displayed.
    
            Attributes
            ----------
            data_name : str
                Name of data.
    
            """
            projection = self._do_initial_projection(self._init_cond)
    
            time_step = abs(self._cfl_number * self._cell_len / self._wave_speed)
    
            current_time = 0
            iteration = 0
            troubled_cell_history = []
            time_history = []
            while current_time < self._final_time:
                # Adjust for last cell
                cfl_number = self._cfl_number
                if current_time+time_step > self._final_time:
                    time_step = self._final_time-current_time
                    cfl_number = self._wave_speed * time_step / self._cell_len
    
                # Update projection
                projection, troubled_cells = self._update_scheme.step(projection,
                                                                      cfl_number)
                iteration += 1
    
                if (iteration % self._history_threshold) == 0:
                    troubled_cell_history.append(troubled_cells)
                    time_history.append(current_time)
    
                current_time += time_step
    
            # Save approximation results in dictionary
            approx_stats = {'projection': projection, 'time_history': time_history,
                            'troubled_cell_history': troubled_cell_history}
    
            # Encode all ndarrays to fit JSON format
            approx_stats = {key: encode_ndarray(approx_stats[key])
                            for key in approx_stats.keys()}
    
            # Save approximation results in JSON format
            with open(self._plot_dir+'/' + data_name + '.json', 'w') \
                    as json_file:
                json_file.write(json.dumps(approx_stats))
    
            # Read approximation results
            with open(self._plot_dir+'/' + data_name + '.json') as json_file:
                approx_stats = json.load(json_file)
    
            # Decode all ndarrays by converting lists
            approx_stats = {key: decode_ndarray(approx_stats[key])
                            for key in approx_stats.keys()}
    
            # Plot exact/approximate results, errors, shock tubes,
            # and any detector-dependant plots
            self._detector.plot_results(**approx_stats)
    
        def save_plots(self, plot_name):
            """Saves plotted results.
    
            Sets plot directory, if not already existing, and saves plots
            generated during the last approximation.
    
            Parameters
            ----------
            plot_name : str
                Name of plot.
    
            """
            # Set paths for plot files if not existing already
            if not os.path.exists(self._plot_dir):
                os.makedirs(self._plot_dir)
    
            # Save plots
            for identifier in plt.get_figlabels():
                # Set path for figure directory if not existing already
                if not os.path.exists(self._plot_dir + '/' + identifier):
                    os.makedirs(self._plot_dir + '/' + identifier)
    
                plt.figure(identifier)
                plt.savefig(self._plot_dir + '/' + identifier + '/' +
                            plot_name + '.pdf')
    
        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
            self._basis = OrthonormalLegendre(self._polynomial_degree)
    
            # Set additional necessary config parameters
            self._limiter_config['cell_len'] = self._cell_len
    
            # Set mesh with one ghost point on each side
            self._mesh = np.arange(self._left_bound - (3/2*self._cell_len),
                                   self._right_bound + (5/2*self._cell_len),
                                   self._cell_len)  # +3/2
    
            # Set inverse mass matrix
            mass_matrix = []
            for i in range(self._polynomial_degree+1):
                new_row = []
                for j in range(self._polynomial_degree+1):
                    new_entry = 0.0
                    if i == j:
                        new_entry = 1.0
                    new_row.append(new_entry)
                mass_matrix.append(new_row)
            self._inv_mass = np.array(mass_matrix)
    
        def _do_initial_projection(self, initial_condition, adjustment=0):
            """Calculates initial projection.
    
            Calculates a projection at time step 0 and adds ghost cells on both
            sides of the array.
    
            Parameters
            ----------
            initial_condition : InitialCondition object
                Initial condition used for calculation. May differ from instance
                variable.
            adjustment: float
                Extent of adjustment of each evaluation point in x-direction.
    
            Returns
            -------
            ndarray
                Matrix containing projection of size (N+2, p+1) with N being the
                number of grid cells and p being the polynomial degree.
    
            """
            # Initialize matrix and set first entry to accommodate for ghost cell
            output_matrix = [0]
            basis_vector = self._basis.get_basis_vector()
    
            for cell in range(self._num_grid_cells):
                new_row = []
                eval_point = self._left_bound + (cell+0.5)*self._cell_len
    
                for degree in range(self._polynomial_degree + 1):
                    new_entry = sum(
                        initial_condition.calculate(
                            eval_point + self._cell_len/2
                            * self._quadrature.get_eval_points()[point]
                            - adjustment)
                        * basis_vector[degree].subs(
                            x, self._quadrature.get_eval_points()[point])
                        * self._quadrature.get_weights()[point]
                        for point in range(self._quadrature.get_num_points()))
                    new_row.append(np.float64(new_entry))
    
                new_row = np.array(new_row)
                output_matrix.append(self._inv_mass @ new_row)
    
            # Set ghost cells to respective value
            output_matrix[0] = output_matrix[self._num_grid_cells]
            output_matrix.append(output_matrix[1])
    
            # print(np.array(output_matrix).shape)
            return np.transpose(np.array(output_matrix))
    
        def build_training_data(self, adjustment, stencil_length,
                                add_reconstructions, initial_condition=None):
            """Builds training data set.
    
            Initializes projection and calculates cell averages and
            reconstructions for it.
    
            Parameters
            ----------
            adjustment : float
                Extent of adjustment of each evaluation point in x-direction.
            stencil_length : int
                Size of training data array.
            initial_condition : InitialCondition object, optional
                Initial condition used for calculation.
                Default: None (i.e. instance variable).
    
            Returns
            -------
            ndarray
                Matrix containing cell averages and reconstructions for initial
                projection.
    
            """
            if initial_condition is None:
                initial_condition = self._init_cond
            projection = self._do_initial_projection(initial_condition, adjustment)
    
            return self._detector.calculate_cell_average_and_reconstructions(
                projection[:, 1:-1], stencil_length)