transport_dg.hh

Go to the documentation of this file.

/*!
 *
 * Copyright (C) 2015 Technical University of Liberec.  All rights reserved.
 * 
 * This program is free software; you can redistribute it and/or modify it under
 * the terms of the GNU General Public License version 3 as published by the
 * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
 * 
 * This program is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 * FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
 *
 * 
 * @file    transport_dg.hh
 * @brief   Discontinuous Galerkin method for equation of transport with dispersion.
 * @author  Jan Stebel
 */
 
#ifndef TRANSPORT_DG_HH_
#define TRANSPORT_DG_HH_
 
#include <math.h>                              // for fabs
#include <string.h>                            // for memcpy
#include <algorithm>                           // for max
 
#include <string>                              // for operator<<
#include <vector>                              // for vector
#include <armadillo>
#include "fem/update_flags.hh"                 // for operator|
#include "fields/field_values.hh"              // for FieldValue<>::Scalar
#include "fields/field.hh"
#include "fields/field_fe.hh"
#include "fields/multi_field.hh"
#include "la/vector_mpi.hh"
#include "fields/equation_output.hh"
#include "la/linsys.hh"
#include "input/accessors.hh"                  // for ExcAccessorForNullStorage
#include "input/accessors_impl.hh"             // for Record::val
#include "input/storage.hh"                    // for ExcStorageTypeMismatch
#include "input/type_base.hh"                  // for Array
#include "input/type_generic.hh"               // for Instance
#include "input/type_record.hh"                // for Record::ExcRecordKeyNo...
#include "system/index_types.hh"               // for LongIdx
#include "mesh/accessors.hh"                   // for ElementAccessor, SIdeIter
#include "mesh/elements.h"                     // for Element::dim, Element:...
#include "mesh/neighbours.h"                   // for Neighbour::element
#include "mpi.h"                               // for MPI_Comm_rank
#include "petscmat.h"                          // for Mat, MatDestroy
#include "petscvec.h"                          // for Vec, VecDestroy, VecSc...
#include "transport/concentration_model.hh"    // for ConcentrationTransport...
#include "transport/heat_model.hh"             // for HeatTransferModel, Hea...
#include "tools/mixed.hh"
class DiscreteSpace;
class Distribution;
class OutputTime;
class DOFHandlerMultiDim;
class GenericAssemblyBase;
template<unsigned int dim, class Model> class AssemblyDG;
template<unsigned int dim, class Model> class MassAssemblyDG;
template<unsigned int dim, class Model> class StiffnessAssemblyDG;
template<unsigned int dim, class Model> class SourcesAssemblyDG;
template<unsigned int dim, class Model> class BdrConditionAssemblyDG;
template<unsigned int dim, class Model> class InitConditionAssemblyDG;
template<unsigned int dim, class Model> class InitProjectionAssemblyDG;
template< template<IntDim...> class DimAssembly> class GenericAssembly;
template<unsigned int dim, unsigned int spacedim> class FEValuesBase;
template<unsigned int dim> class FiniteElement;
template<unsigned int dim, unsigned int spacedim> class Mapping;
class Quadrature;
namespace Input { namespace Type { class Selection; } }
class ElementCacheMap;
 
/*class FEObjects {
public:
 
    inline FEObjects(Mesh *mesh_, unsigned int fe_order)
    {
        fe = MixedPtr<FE_P_disc>(fe_order);
        fe_rt = MixedPtr<FE_RT0>();
        q = MixedPtr<QGauss>(2*fe_order);
 
        auto ds = std::make_shared<EqualOrderDiscreteSpace>(mesh_, fe);
        dh = std::make_shared<DOFHandlerMultiDim>(*mesh_);
        dh->distribute_dofs(ds_);
    }
 
    ~FEObjects();
 
    MixedPtr<FiniteElement> fe;
    MixedPtr<FiniteElement> fe_rt;
    MixedPtr<Quadrature> q;
 
    /// Object for distribution of dofs.
    std::shared_ptr<DOFHandlerMultiDim> dh;
};*/
 
 
 
/**
 * @brief Transport with dispersion implemented using discontinuous Galerkin method.
 *
 * TransportDG implements the discontinuous Galerkin method for the transport and diffusion of substances.
 * The concentration @f$ c_i ~[kg/m^3]@f$ of the i-th substance is governed by the advection-diffusion equation
 * @f[
 *      \partial_t c_i + \mathbf v\cdot\nabla c_i - \mathrm{div}(D\nabla c_i) = F \mbox{ in }\Omega^d,
 * @f]
 * where @f$\mathbf v@f$ is the fluid velocity and @f$\Omega^d@f$ the @f$d@f$-dimensional domain, respectively.
 * The hydrodynamic dispersivity tensor @f$\mathbf D ~[m^2/s]@f$ is given by:
 * @f[
 *      \mathbf D = D_m\mathbf I + |\mathbf v|\left(\alpha_T\mathbf I + (\alpha_L-\alpha_T)\frac{\mathbf v\otimes\mathbf v}{|\mathbf v|^2}\right).
 * @f]
 * The molecular dispersivity @f$D_m~[m^2/s]@f$, as well as the longitudal and transversal dispersivity @f$\alpha_L,~\alpha_T~[m]@f$ are input parameters of the model.
 *
 * For lower dimensions @f$d=1,2@f$ the advection-diffusion equation is multiplied by the fracture cross-cut @f$\delta^d~[m^{3-d}]@f$.
 *
 * The boundary @f$\partial\Omega^d@f$ is divided into three disjoint parts @f$\Gamma^d_D\cup\Gamma^d_N\cup\Gamma^d_F@f$.
 * We prescribe the following boundary conditions:
 * @f{eqnarray*}{
 *      c_i^d &= c_{iD}^d &\mbox{ on }\Gamma^d_D \mbox{ (Dirichlet)},\\
 *      \mathbf D^d\nabla c_i^d\cdot\mathbf n &= 0 &\mbox{ on }\Gamma^d_N \mbox{ (Neumann)},
 * @f}
 * The transfer of mass through fractures is described by the transmission conditions on @f$\Gamma^d_F@f$:
 * @f[
 *      -\mathbf D^d\nabla c_i^d\cdot\mathbf n = \sigma(c_i^d-c_i^{d-1}) + \left\{\begin{array}{cl}0 &\mbox{ if }\mathbf v^d\cdot\mathbf n\ge 0\\\mathbf v^d\cdot\mathbf n(c_i^{d-1}-c_i^d) & \mbox{ if }\mathbf v^d\cdot\mathbf n<0\end{array}\right.,\qquad
 *      F^{d-1} = (\sigma + |\mathbf v^d\cdot\mathbf n|)(c_i^d-c_i^{d-1}).
 * @f]
 * Here @f$\mathbf n@f$ stands for the unit outward normal vector to @f$\partial\Omega^d@f$.
 * The coefficient @f$\sigma@f$ determines the transfer of mass through fractures due to diffusion.
 *
 * @ingroup transport_mod
 *
 */
template<class Model>
class TransportDG : public Model
{
public:
 
    template<unsigned int dim> using MassAssemblyDim = MassAssemblyDG<dim, Model>;
    template<unsigned int dim> using StiffnessAssemblyDim = StiffnessAssemblyDG<dim, Model>;
    template<unsigned int dim> using SourcesAssemblyDim = SourcesAssemblyDG<dim, Model>;
    template<unsigned int dim> using BdrConditionAssemblyDim = BdrConditionAssemblyDG<dim, Model>;
    template<unsigned int dim> using InitConditionAssemblyDim = InitConditionAssemblyDG<dim, Model>;
    template<unsigned int dim> using InitProjectionAssemblyDim = InitProjectionAssemblyDG<dim, Model>;
 
    typedef std::vector<std::shared_ptr<FieldFE< 3, FieldValue<3>::Scalar>>> FieldFEScalarVec;
 
    class EqFields : public Model::ModelEqFields {
    public:
 
        EqFields();
 
        MultiField<3, FieldValue<3>::Scalar> fracture_sigma;    ///< Transition parameter for diffusive transfer on fractures (for each substance).
        MultiField<3, FieldValue<3>::Scalar> dg_penalty;        ///< Penalty enforcing inter-element continuity of solution (for each substance).
        Field<3, FieldValue<3>::Scalar> region_id;
        Field<3, FieldValue<3>::Scalar> subdomain;
 
        EquationOutput output_fields;
    };
 
 
    class EqData : public Model::ModelEqData {
    public:
 
        EqData() {}
 
        /**
         * @brief Sets up parameters of the DG method on a given boundary edge.
         *
         * Assumption is that the edge consists of only 1 side.
         * @param side              The boundary side.
         * @param K_size            Size of vector of tensors K.
         * @param K                 Dispersivity tensor.
         * @param ad_vector         Advection vector.
         * @param normal_vector     Normal vector (assumed constant along the edge).
         * @param alpha             Penalty parameter that influences the continuity
         *                          of the solution (large value=more continuity).
         * @param gamma             Computed penalty parameters.
         */
        void set_DG_parameters_boundary(Side side,
                    const int K_size,
                    const std::vector<arma::mat33> &K,
                    const double flux,
                    const arma::vec3 &normal_vector,
                    const double alpha,
                    double &gamma);
 
 
        inline void set_time_governor(TimeGovernor *time) {
            ASSERT_PTR(time);
            this->time_ = time;
        }
 
 
        /// Compute and return anisotropy of given element
        double elem_anisotropy(ElementAccessor<3> e) const;
 
 
        /// @name Parameters of the numerical method
        // @{
 
        /// Penalty parameters.
        std::vector<std::vector<double> > gamma;
 
        /// DG variant ((non-)symmetric/incomplete
        int dg_variant;
 
        /// Polynomial order of finite elements.
        unsigned int dg_order;
 
        // @}
 
 
        /// Auxiliary vectors for calculation of sources in balance due to retardation (e.g. sorption).
        std::vector<Vec> ret_vec;
 
        /// Linear algebra system for the transport equation.
        LinSys **ls;
 
        /// Linear algebra system for the time derivative (actually it is used only for handling the matrix structures).
        LinSys **ls_dt;
 
        /// @name Auxiliary fields used during assembly
        // @{
 
        /// Diffusion coefficients.
        vector<vector<arma::mat33> > dif_coef;
        /// Maximal number of edge sides (evaluate from dim 1,2,3)
        unsigned int max_edg_sides;
 
        // @}
 
        /// Object for distribution of dofs.
        std::shared_ptr<DOFHandlerMultiDim> dh_;
 
        /// Vector of solution data.
        std::vector<VectorMPI> output_vec;
 
        FieldFEScalarVec conc_fe;
        std::shared_ptr<DOFHandlerMultiDim> dh_p0;
        TimeGovernor *time_;
        std::shared_ptr<Balance> balance_;
    };
 
    enum DGVariant {
        // Non-symmetric weighted interior penalty DG
        non_symmetric = -1,
 
        // Incomplete weighted interior penalty DG
        incomplete = 0,
 
        // Symmetric weighted interior penalty DG
        symmetric = 1
    };
 
    /**
     * @brief Constructor.
     * @param init_mesh         computational mesh
     * @param in_rec            input record
     */
    TransportDG(Mesh &init_mesh, const Input::Record in_rec);
    /**
 
     * @brief Declare input record type for the equation TransportDG.
     */
    static const Input::Type::Record & get_input_type();
 
    /**
     * @brief Input type for the DG variant selection.
     */
    static const Input::Type::Selection & get_dg_variant_selection_input_type();
 
    /**
     * @brief Initialize solution in the zero time.
     */
    void zero_time_step() override;
    
    bool evaluate_time_constraint(double &)
    { return false; }
 
    /**
     * @brief Computes the solution in one time instant.
     */
    void update_solution() override;
 
    /**
     * @brief Postprocesses the solution and writes to output file.
     */
    void output_data();
 
    /**
     * @brief Destructor.
     */
    ~TransportDG() override;
 
    void initialize() override;
 
    void calculate_cumulative_balance();
 
    /// Return PETSc vector with solution for sbi-th component.
    Vec get_component_vec(unsigned int sbi)
    { return eq_data_->ls[sbi]->get_solution(); }
 
    /// Getter for P0 interpolation by FieldFE.
    FieldFEScalarVec& get_p0_interpolation()
    { return eq_data_->conc_fe;}
 
    /// Compute P0 interpolation of the solution (used in reaction term).
    void compute_p0_interpolation();
 
    void update_after_reactions(bool solution_changed);
 
    /// Access to balance object of Model
    inline std::shared_ptr<Balance> balance() const {
        return Model::balance_;
    }
 
    inline typename Model::ModelEqFields &eq_fields() { return *eq_fields_; }
 
    inline typename Model::ModelEqData &eq_data() { return *eq_data_; }
 
private:
    /// Registrar of class to factory
    static const int registrar;
 
    void preallocate();
 
    /**
     * @brief Calculates the dispersivity (diffusivity) tensor from the velocity field.
     *
     * @param K        The computed dispersivity tensor.
     * @param velocity The velocity field (at quadrature points).
     * @param Dm       Molecular diffusivities.
     * @param alphaL   Longitudal dispersivities.
     * @param alphaT   Transversal dispersivities.
     * @param porosity  Porosities.
     * @param cross_cut Cross-cuts of higher dimension.
     */
//  void calculate_dispersivity_tensor(arma::mat33 &K, const arma::vec3 &velocity,
//          double Dm, double alphaL, double alphaT, double porosity,
//          double cross_cut);
 
    /**
     * @brief Sets the initial condition.
     *
     * This method just calls AssemblyDG::prepare_initial_condition() for each elements.
     */
    void set_initial_condition();
 
    
    
    void output_region_statistics();
 
 
 
    /// @name Physical parameters
    // @{
 
    /// Fields for model parameters.
    std::shared_ptr<EqFields> eq_fields_;
 
    /// Data for model parameters.
    std::shared_ptr<EqData> eq_data_;
 
    // @}
 
 
 
    /// @name Solution of algebraic system
    // @{
 
    /// Vector of right hand side.
    std::vector<Vec> rhs;
 
    /// The stiffness matrix.
    std::vector<Mat> stiffness_matrix;
 
    /// The mass matrix.
    std::vector<Mat> mass_matrix;
    
    /// Mass from previous time instant (necessary when coefficients of mass matrix change in time).
    std::vector<Vec> mass_vec;
    // @}
 
 
    /// @name Output to file
    // @{
 
    /// Array for storing the output solution data.
    //vector<double*> output_solution;
 
    /// Record with input specification.
    Input::Record input_rec;
    
    ofstream reg_stat_stream;
 
 
    // @}
 
 
    /// @name Auxiliary fields used during assembly
    // @{
 
    /// Temporary values of increments due to retardation (e.g. sorption)
    vector<double> ret_sources, ret_sources_prev;
 
    // @}
 
 
 
 
    /// @name Other
    // @{
 
    /// Indicates whether matrices have been preallocated.
    bool allocation_done;
 
    bool init_projection;
    // @}
 
    /// general assembly objects, hold assembly objects of appropriate dimension
    GenericAssembly< MassAssemblyDim > * mass_assembly_;
    GenericAssembly< StiffnessAssemblyDim > * stiffness_assembly_;
    GenericAssembly< SourcesAssemblyDim > * sources_assembly_;
    GenericAssembly< BdrConditionAssemblyDim > * bdr_cond_assembly_;
    GenericAssemblyBase * init_assembly_;
};
 
 
 
 
 
 
#endif /* TRANSPORT_DG_HH_ */