Coverage for src/pymor/discretizers/advection : 24%

# This file is part of the pyMOR project (http://www.pymor.org). 

# Copyright Holders: Rene Milk, Stephan Rave, Felix Schindler 

# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause) 

from __future__ import absolute_import, division, print_function 

import numpy as np 

from pymor.algorithms.timestepping import ExplicitEulerTimeStepper 

from pymor.analyticalproblems.advection import InstationaryAdvectionProblem 

from pymor.core import inject_sid 

from pymor.discretizations import InstationaryDiscretization 

from pymor.domaindiscretizers import discretize_domain_default 

from pymor.gui.qt import PatchVisualizer, Matplotlib1DVisualizer 

from pymor.la import NumpyVectorArray 

from pymor.operators import NumpyGenericOperator 

from pymor.operators.fv import (nonlinear_advection_lax_friedrichs_operator, 

                                nonlinear_advection_engquist_osher_operator, 

                                nonlinear_advection_simplified_engquist_osher_operator, 

                                L2Product, L2ProductFunctional) 

def discretize_nonlinear_instationary_advection_fv(analytical_problem, diameter=None, nt=100, num_flux='lax_friedrichs', 

                                                   lxf_lambda=1., eo_gausspoints=5, eo_intervals=1, num_values=None, 

                                                   domain_discretizer=None, grid=None, boundary_info=None): 

    '''Discretizes an |InstationaryAdvectionProblem| using the finite volume method. 

    Simple explicit Euler time-stepping is used for time-discretization. 

    Parameters 

    ---------- 

    analytical_problem 

        The |InstationaryAdvectionProblem| to discretize. 

    diameter 

        If not None, is passed to the domain_discretizer. 

    nt 

        The number of time-steps. 

    num_flux 

        The numerical flux to use in the finite volume formulation. Allowed 

        values are `'lax_friedrichs'`, `'engquist_osher'`, `'simplified_engquist_osher'`. 

        (See :mod:`pymor.operators.fv`.) 

    lxf_lambda 

        The stabilization parameter for the Lax-Friedrichs numerical flux. 

        (Ignored, if different flux is chosen.) 

    eo_gausspoints 

        Number of Gauss points for the Engquist-Osher numerical flux. 

        (Ignored, if different flux is chosen.) 

    eo_intervals 

        Number of sub-intervals to use for integration when using Engquist-Osher 

        numerical flux. (Ignored, if different flux is chosen.) 

    num_values 

        The number of returned vectors of the solution trajectory. If `None`, each 

        intermediate vector that is calculated is returned. 

    domain_discretizer 

        Discretizer to be used for discretizing the analytical domain. This has 

        to be a function `domain_discretizer(domain_description, diameter, ...)`. 

        If further arguments should be passed to the discretizer, use 

        :func:`functools.partial`. If `None`, |discretize_domain_default| is used. 

    grid 

        Instead of using a domain discretizer, the |Grid| can also be passed directly 

        using this parameter. 

    boundary_info 

        A |BoundaryInfo| specifying the boundary types of the grid boundary entities. 

        Must be provided if `grid` is provided. 

    Returns 

    ------- 

    discretization 

        The discretization that has been generated. 

    data 

        Dictionary with the following entries: 

            :grid:           The generated |Grid|. 

            :boundary_info:  The generated |BoundaryInfo|. 

    ''' 

    assert isinstance(analytical_problem, InstationaryAdvectionProblem) 

    assert grid is None or boundary_info is not None 

    assert boundary_info is None or grid is not None 

    assert grid is None or domain_discretizer is None 

    assert num_flux in ('lax_friedrichs', 'engquist_osher', 'simplified_engquist_osher') 

    if grid is None: 

        domain_discretizer = domain_discretizer or discretize_domain_default 

        if diameter is None: 

            grid, boundary_info = domain_discretizer(analytical_problem.domain) 

        else: 

            grid, boundary_info = domain_discretizer(analytical_problem.domain, diameter=diameter) 

    p = analytical_problem 

    if num_flux == 'lax_friedrichs': 

        L = nonlinear_advection_lax_friedrichs_operator(grid, boundary_info, p.flux_function, 

                                                        dirichlet_data=p.dirichlet_data, lxf_lambda=lxf_lambda) 

    elif num_flux == 'engquist_osher': 

        L = nonlinear_advection_engquist_osher_operator(grid, boundary_info, p.flux_function, 

                                                        p.flux_function_derivative, 

                                                        gausspoints=eo_gausspoints, intervals=eo_intervals, 

                                                        dirichlet_data=p.dirichlet_data) 

    else: 

        L = nonlinear_advection_simplified_engquist_osher_operator(grid, boundary_info, p.flux_function, 

                                                                   p.flux_function_derivative, 

                                                                   dirichlet_data=p.dirichlet_data) 

    F = None if p.rhs is None else L2ProductFunctional(grid, p.rhs) 

    if p.initial_data.parametric: 

        def initial_projection(U, mu): 

            I = p.initial_data.evaluate(grid.quadrature_points(0, order=2), mu).squeeze() 

            I = np.sum(I * grid.reference_element.quadrature(order=2)[1], axis=1) * (1. / grid.reference_element.volume) 

            I = NumpyVectorArray(I, copy=False) 

            return I.lincomb(U).data 

        I = NumpyGenericOperator(initial_projection, dim_range=grid.size(0), linear=True, 

                                 parameter_type=p.initial_data.parameter_type) 

    else: 

        I = p.initial_data.evaluate(grid.quadrature_points(0, order=2)).squeeze() 

        I = np.sum(I * grid.reference_element.quadrature(order=2)[1], axis=1) * (1. / grid.reference_element.volume) 

        I = NumpyVectorArray(I, copy=False) 

    inject_sid(I, __name__ + '.discretize_nonlinear_instationary_advection_fv.initial_data', p.initial_data, grid) 

    products = {'l2': L2Product(grid, boundary_info)} 

    if grid.dim == 2: 

        visualizer = PatchVisualizer(grid=grid, bounding_box=grid.domain, codim=0) 

    elif grid.dim == 1: 

        visualizer = Matplotlib1DVisualizer(grid, codim=0) 

    else: 

        visualizer = None 

    parameter_space = p.parameter_space if hasattr(p, 'parameter_space') else None 

    time_stepper = ExplicitEulerTimeStepper(nt=nt) 

    discretization = InstationaryDiscretization(operator=L, rhs=F, initial_data=I, T=p.T, products=products, 

                                                time_stepper=time_stepper, 

                                                parameter_space=parameter_space, visualizer=visualizer, 

                                                num_values=num_values, name='{}_FV'.format(p.name)) 

    return discretization, {'grid': grid, 'boundary_info': boundary_info}