darcy_flow_lmh.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    darcy_flow_lmh.hh
 * @brief   Lumped mixed-hybrid model of linear Darcy flow, possibly unsteady.
 * @author  Jan Brezina
 *
 * Main object for mixed-hybrid discretization of the linear elliptic PDE (Laplace)
 * on  a multidimensional domain. Discretization of saturated Darcy flow using
 * RT0 approximation for the velocity
 */
 
/*
 * list of files dependent on this one:
 *
 * posprocess.cc
 * problem.cc
 * main.hh
 * transport.cc
 */
 
 
#ifndef DARCY_FLOW_LMH_HH
#define DARCY_FLOW_LMH_HH
 
#include <petscmat.h>                           // for Mat
#include <string.h>                             // for memcpy
#include <algorithm>                            // for swap
 
#include <memory>                               // for shared_ptr, allocator...
#include <new>                                  // for operator new[]
#include <string>                               // for string, operator<<
#include <vector>                               // for vector, vector<>::con...
#include <armadillo>
#include "fields/bc_field.hh"                   // for BCField
#include "fields/field.hh"                      // for Field
#include "fields/field_set.hh"                  // for FieldSet
#include "fields/field_values.hh"               // for FieldValue<>::Scalar
#include "flow/darcy_flow_interface.hh"         // for DarcyFlowInterface
#include "input/input_exception.hh"             // for DECLARE_INPUT_EXCEPTION
#include "input/type_base.hh"                   // for Array
#include "input/type_generic.hh"                // for Instance
#include "mesh/mesh.h"                          // for Mesh
#include "petscvec.h"                           // for Vec, _p_Vec, VecScatter
#include "system/exceptions.hh"                 // for ExcStream, operator<<
#include "tools/time_governor.hh"               // for TimeGovernor
#include "la/vector_mpi.hh"                     // for VectorMPI
 
 
class Balance;
class DarcyFlowMHOutput;
class Element;
class Intersection;
class LinSys;
// class LinSys_BDDC;
class LocalSystem;
class LocalConstraint;
namespace Input {
    class AbstractRecord;
    class Record;
    namespace Type {
        class Record;
        class Selection;
    }
}
template<unsigned int dim, class TEqData> class ReadInitCondAssemblyLMH;
template<unsigned int dim, class TEqData> class MHMatrixAssemblyLMH;
template<unsigned int dim, class TEqData> class ReconstructSchurAssemblyLMH;
class GenericAssemblyBase;
template< template<IntDim...> class DimAssembly> class GenericAssembly;
 
 
/**
 * @brief Mixed-hybrid model of linear Darcy flow, possibly unsteady.
 *
 * Abstract class for various implementations of Darcy flow. In future there should be
 * one further level of abstraction for general time dependent problem.
 *
 * maybe TODO:
 * split compute_one_step to :
 * 1) prepare_next_timestep
 * 2) actualize_solution - this is for iterative nonlinear solvers
 *
 */
 
 
/**
 * @brief Mixed-hybrid of steady Darcy flow with sources and variable density.
 *
 * solve equations:
 * @f[
 *      q= -{\mathbf{K}} \nabla h -{\mathbf{K}} R \nabla z
 * @f]
 * @f[
 *      \mathrm{div} q = f
 * @f]
 *
 * where
 * - @f$ q @f$ is flux @f$[ms^{-1}]@f$ for 3d, @f$[m^2s^{-1}]@f$ for 2d and @f$[m^3s^{-1}]@f$ for 1d.
 * - @f$ \mathbf{K} @f$ is hydraulic tensor ( its orientation for 2d, 1d case is questionable )
 * - @f$ h = \frac{\pi}{\rho_0 g}+z @f$ is pressure head, @f$ \pi, \rho_0, g @f$ are the pressure, water density, and acceleration of gravity , respectively.
 *   Assumes gravity force acts counter to the direction of the @f$ z @f$ axis.
 * - @f$ R @f$ is destity or gravity variability coefficient. For density driven flow it should be
 * @f[
 *    R = \frac{\rho}{\rho_0} -1 = \rho_0^{-1}\sum_{i=1}^s c_i
 * @f]
 *   where @f$ c_i @f$ is concentration in @f$ kg m^{-3} @f$.
 *
 * The time key is optional, when not specified the equation is forced to steady regime. Using Steady TimeGovernor which have no dt constraints.
 *
 *
 * TODO:
 * Make solution regular field (need FeSeystem and parallel DofHandler for edge pressures), then remove get_solution_vector from
 * Equation interface.
 */
/**
 * Model for transition coefficients due to Martin, Jaffre, Roberts (see manual for full reference)
 *
 * TODO:
 * - how we can reuse field values computed during assembly
 *
 */
 
class DarcyLMH : public DarcyFlowInterface
{
public:
    TYPEDEF_ERR_INFO( EI_Reason, string);
    DECLARE_EXCEPTION(ExcSolverDiverge,
            << "Diverged nonlinear solver. Reason: " << EI_Reason::val
             );
    DECLARE_INPUT_EXCEPTION(ExcMissingTimeGovernor,
            << "Missing the key 'time', obligatory for the transient problems.");
 
    /** Class with all fields used in the equation DarcyFlow.
    * This is common to all implementations since this provides interface
    * to this equation for possible coupling.
    * 
    * This class uses the common output class DarcyFlowMHOutput.
    * It is base class of RichardsLMH::EqFields.
    * */
    class EqFields : public FieldSet {
    public:
        /**
         * For compatibility with old BCD file we have to assign integer codes starting from 1.
         */
        enum BC_Type {
            none=0,
            dirichlet=1,
            total_flux=4,
            seepage=5,
            river=6
        };
 
        /// Return a Selection corresponding to enum BC_Type.
        static const Input::Type::Selection & get_bc_type_selection();
 
        /// Creation of all fields.
        EqFields();
 
        /// Return coords field
        FieldCoords &X() {
            return this->X_;
        }
 
 
        Field<3, FieldValue<3>::TensorFixed > anisotropy;
        Field<3, FieldValue<3>::Scalar > conductivity;
        Field<3, FieldValue<3>::Scalar > cross_section;
        Field<3, FieldValue<3>::Scalar > water_source_density;
        Field<3, FieldValue<3>::Scalar > water_source_sigma;
        Field<3, FieldValue<3>::Scalar > water_source_ref_pressure;
        Field<3, FieldValue<3>::Scalar > sigma;
 
        BCField<3, FieldValue<3>::Enum > bc_type; // Discrete need Selection for initialization
        BCField<3, FieldValue<3>::Scalar > bc_pressure;
        BCField<3, FieldValue<3>::Scalar > bc_flux;
        BCField<3, FieldValue<3>::Scalar > bc_robin_sigma;
        BCField<3, FieldValue<3>::Scalar > bc_switch_pressure;
 
        Field<3, FieldValue<3>::Scalar > init_pressure;
        Field<3, FieldValue<3>::Scalar > storativity;
        Field<3, FieldValue<3>::Scalar > extra_storativity; /// Externally added storativity.
        Field<3, FieldValue<3>::Scalar > extra_source; /// Externally added water source.
 
        Field<3, FieldValue<3>::Scalar> field_ele_pressure;
        Field<3, FieldValue<3>::Scalar> field_ele_piezo_head;
        Field<3, FieldValue<3>::VectorFixed > field_ele_velocity;
        Field<3, FieldValue<3>::VectorFixed > flux;
        Field<3, FieldValue<3>::Scalar> field_edge_pressure;
 
        Field<3, FieldValue<3>::VectorFixed > gravity_field; /// Holds gravity vector acceptable in FieldModel
        BCField<3, FieldValue<3>::VectorFixed > bc_gravity; /// Same as previous but used in boundary fields
        Field<3, FieldValue<3>::Scalar> init_piezo_head;
        BCField<3, FieldValue<3>::Scalar> bc_piezo_head;
        BCField<3, FieldValue<3>::Scalar> bc_switch_piezo_head;
 
        Field<3, FieldValue<3>::Scalar> ref_pressure; /// Precompute l2 difference outputs
        Field<3, FieldValue<3>::VectorFixed> ref_velocity;
        Field<3, FieldValue<3>::Scalar> ref_divergence;
    };
 
    class EqData {
    public:
        typedef DarcyLMH::EqFields EqFields;
 
        EqData(shared_ptr<EqFields> eq_fields);
        
        void init();     ///< Initialize vectors, ...
        void reset();    ///< Reset data members
 
        /// Shared pointer of EqFields
        std::shared_ptr<EqFields> eq_fields_;
 
        /**
         * Gravity vector and constant shift of pressure potential. Used to convert piezometric head
         * to pressure head and vice versa.
         */
        arma::vec4 gravity_;
        arma::vec3 gravity_vec_;
 
        // Mirroring the following members of DarcyLMH:
        Mesh *mesh;
        std::shared_ptr<DOFHandlerMultiDim> dh_;         ///< full DOF handler represents DOFs of sides, elements and edges
        std::shared_ptr<SubDOFHandlerMultiDim> dh_cr_;   ///< DOF handler represents DOFs of edges
        std::shared_ptr<DOFHandlerMultiDim> dh_cr_disc_; ///< DOF handler represents DOFs of sides
        std::shared_ptr<SubDOFHandlerMultiDim> dh_p_;    ///< DOF handler represents DOFs of element pressure
 
 
        uint water_balance_idx;
 
        // TODO: check usage of MortarMethod in LMH
        MortarMethod mortar_method_;
 
        int is_linear;              ///< Hack fo BDDC solver.
        bool force_no_neumann_bc;       ///< auxiliary flag for switchting Dirichlet like BC
 
        /// Idicator of dirichlet or neumann type of switch boundary conditions.
        std::vector<char> bc_switch_dirichlet;
 
        VectorMPI full_solution;     //< full solution [vel,press,lambda] from 2. Schur complement
        
        // Propagate test for the time term to the assembly.
        // This flag is necessary for switching BC to avoid setting zero neumann on the whole boundary in the steady case.
        bool use_steady_assembly_;
        
        // for time term assembly
        double time_step_;
        
        std::shared_ptr<LinSys> lin_sys_schur;  //< Linear system of the 2. Schur complement.
        VectorMPI p_edge_solution;               //< 2. Schur complement solution
        VectorMPI p_edge_solution_previous;      //< 2. Schur complement previous solution (iterative)
        VectorMPI p_edge_solution_previous_time; //< 2. Schur complement previous solution (time)
 
        std::map<LongIdx, LocalSystem> seepage_bc_systems;
 
        /// Shared Balance object
        std::shared_ptr<Balance> balance_;
 
        unsigned int nonlinear_iteration_; //< Actual number of completed nonlinear iterations, need to pass this information into assembly.
 
        /// Following data members are stored in vectors, one item for every cell.
//        /** TODO: Investigate why the hell do we need this flag.
//        *  If removed, it does not break any of the integration tests,
//        * however it must influence the Dirichlet rows in matrix.
//        */
//        std::vector<unsigned int> dirichlet_edge;
 
        std::vector<LocalSystem> loc_system_;
        std::vector<LocalConstraint> loc_constraint_;
        std::vector<arma::vec> postprocess_solution_;
        std::array<std::vector<unsigned int>, 3> loc_side_dofs;
        std::array<std::vector<unsigned int>, 3> loc_edge_dofs;
        std::array<unsigned int, 3> loc_ele_dof;
 
//        // std::shared_ptr<MortarAssemblyBase> mortar_assembly;
 
        std::vector<bool> save_local_system_;       ///< Flag for saving the local system. Currently used only in case of seepage BC.
        std::vector<bool> bc_fluxes_reconstruted;   ///< Flag indicating whether the fluxes for seepage BC has been reconstructed already.
        std::array<unsigned int, 3> schur_offset_;  ///< Index offset in the local system for the Schur complement (of dim = 1,2,3).
    };
 
    template<unsigned int dim> using ReadInitCondAssemblyLMHDim = ReadInitCondAssemblyLMH<dim, EqData>;
    template<unsigned int dim> using MHMatrixAssemblyLMHDim = MHMatrixAssemblyLMH<dim, EqData>;
    template<unsigned int dim> using ReconstructSchurAssemblyLMHDim = ReconstructSchurAssemblyLMH<dim, EqData>;
 
    /// Selection for enum MortarMethod.
    static const Input::Type::Selection & get_mh_mortar_selection();
 
 
 
 
 
    DarcyLMH(Mesh &mesh, const Input::Record in_rec, TimeGovernor *tm = nullptr);
 
    static const Input::Type::Record & type_field_descriptor();
    static const Input::Type::Record & get_input_type();
 
    void init_eq_data();
    void initialize() override;
    virtual void initialize_specific();
    void zero_time_step() override;
    void update_solution() override;
 
    /// Solve the problem without moving to next time and without output.
    void solve_time_step(bool output = true);
    
    /// postprocess velocity field (add sources)
    virtual void accept_time_step();
    virtual void postprocess();
    virtual void output_data() override;
 
    virtual double solved_time() override;
 
    inline EqFields &eq_fields() { return *eq_fields_; }
    inline EqData &eq_data() { return *eq_data_; }
 
    /// Sets external storarivity field (coupling with other equation).
    void set_extra_storativity(const Field<3, FieldValue<3>::Scalar> &extra_stor)
    { eq_fields_->extra_storativity = extra_stor; }
 
    /// Sets external source field (coupling with other equation).
    void set_extra_source(const Field<3, FieldValue<3>::Scalar> &extra_src)
    { eq_fields_->extra_source = extra_src; }
 
    virtual ~DarcyLMH() override;
 
 
protected:
    /**
     * Returns true is the fields involved in the time term have values that makes the time term zero.
     * For time_global==true, it returns true if there are no field descriptors in the input list, so the
     * fields )of the time ter) have their default values for whole simulation.
     * If time_global==false (default), only the actual values are considered.
     */
    virtual bool zero_time_term(bool time_global=false);
 
    /// Solve method common to zero_time_step and update solution.
    void solve_nonlinear();
 
    /**
     * Create and preallocate MH linear system (including matrix, rhs and solution vectors)
     */
    void create_linear_system(Input::AbstractRecord rec);
 
    /**
     * Read initial condition into solution vector.
     * Must be called after create_linear_system.
     *
     * For the LMH scheme we have to be able to save edge pressures in order to
     * restart simulation or use results of one simulation as initial condition for other one.
     */
//    void read_initial_condition();
    
    /**
     * In some circumstances, the intial condition must be processed.
     * It is called at the end of @p read_initial_condition().
     * This is used in Richards equation due the update of water content.
     */
//    virtual void initial_condition_postprocess();
    
    /**
     * Allocates linear system matrix for MH.
     * TODO:
     * - use general preallocation methods in DofHandler
     */
    void allocate_mh_matrix();
 
    /**
     * Assembles linear system matrix for MH.
     * Element by element assembly is done using dim-template assembly class.
     * Assembles only steady part of the equation.
     * TODO:
     * - include time term - DONE
     * - add support for Robin type sources
     * - support for nonlinear solvers - assembly either residual vector, matrix, or both (using FADBAD++)
     */
//    void assembly_mh_matrix(MultidimAssembly& assembler);
    
//    void reconstruct_solution_from_schur(MultidimAssembly& assembler);
 
    /**
     * Assembly or update whole linear system.
     */
    virtual void assembly_linear_system();
 
//     void set_mesh_data_for_bddc(LinSys_BDDC * bddc_ls);
    /**
     * Return a norm of residual vector.
     * TODO: Introduce Equation::compute_residual() updating
     * residual field, standard part of EqData.
     */
//    virtual double solution_precision() const;
    
    /// Print darcy flow matrix in matlab format into a file.
    void print_matlab_matrix(string matlab_file);
 
    /// Get vector of all DOF indices of given component (0..side, 1..element, 2..edge)
    std::vector<int> get_component_indices_vec(unsigned int component) const;
 
    /// Getter for the linear system of the 2. Schur complement.
    LinSys& lin_sys_schur()
    { return *(eq_data_->lin_sys_schur); }
 
    /// Create and initialize assembly objects
    virtual void initialize_asm();
 
    /// Call assemble of read_init_cond_assembly_
    virtual void read_init_cond_asm();
 
    std::shared_ptr<Balance> balance_;
 
    DarcyFlowMHOutput *output_object;
 
    int size;                   // global size of MH matrix
 
    bool data_changed_;
 
    // Setting of the nonlinear solver. TODO: Move to the solver class later on.
    double tolerance_;
    unsigned int min_n_it_;
    unsigned int max_n_it_;
 
    std::shared_ptr<EqFields> eq_fields_;
    std::shared_ptr<EqData> eq_data_;
 
 
    friend class DarcyFlowMHOutput;
    //friend class P0_CouplingAssembler;
    //friend class P1_CouplingAssembler;
 
    /// general assembly objects, hold assembly objects of appropriate dimension
    GenericAssembly< ReadInitCondAssemblyLMHDim > * read_init_cond_assembly_;
    GenericAssemblyBase * mh_matrix_assembly_;
    GenericAssemblyBase * reconstruct_schur_assembly_;
private:
    /// Registrar of class to factory
    static const int registrar;
 
    static std::string equation_name()
    { return "Flow_Darcy_LMH";}
};
 
#endif  //DARCY_FLOW_LMH_HH
//-----------------------------------------------------------------------------
// vim: set cindent: