elasticity_upscale.hpp

Go to the documentation of this file.

//==============================================================================
//==============================================================================
#ifndef ELASTICITY_UPSCALE_HPP_
#define ELASTICITY_UPSCALE_HPP_
 
 
#include <opm/common/utility/platform_dependent/disable_warnings.h>
 
#include <dune/common/fmatrix.hh>
#include <opm/common/utility/parameters/ParameterGroup.hpp>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/istl/ilu.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/preconditioners.hh>
#include <opm/grid/CpGrid.hpp>
#include <opm/elasticity/shapefunctions.hpp>
 
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
 
#include <opm/elasticity/asmhandler.hpp>
#include <opm/elasticity/boundarygrid.hh>
#include <opm/elasticity/elasticity.hpp>
#include <opm/elasticity/elasticity_preconditioners.hpp>
#include <opm/elasticity/logutils.hpp>
#include <opm/elasticity/materials.hh>
#include <opm/elasticity/matrixops.hpp>
#include <opm/elasticity/meshcolorizer.hpp>
#include <opm/elasticity/mpc.hh>
#include <opm/elasticity/mortar_schur.hpp>
#include <opm/elasticity/mortar_utils.hpp>
#include <opm/elasticity/mortar_evaluator.hpp>
#include <opm/elasticity/mortar_schur_precond.hpp>
#include <opm/elasticity/uzawa_solver.hpp>
 
#include <opm/input/eclipse/Parser/Parser.hpp>
#include <opm/input/eclipse/Deck/Deck.hpp>
 
namespace Opm {
namespace Elasticity {
 
enum Solver {
  DIRECT,
  ITERATIVE 
};
 
enum Preconditioner {
  AMG,
  FASTAMG,
  SCHWARZ,
  TWOLEVEL,
  UNDETERMINED
};
 
enum MultiplierPreconditioner {
  SIMPLE,        
  SCHUR,         
  SCHURAMG,      
  SCHURSCHWARZ,  
  SCHURTWOLEVEL  
};
 
enum Smoother {
  SMOOTH_SSOR    = 0,
  SMOOTH_SCHWARZ = 1,
  SMOOTH_JACOBI  = 2,
  SMOOTH_ILU     = 4
};
 
struct LinSolParams {
  Solver type;
 
  int restart;
 
  int maxit;
 
  double tol;
 
  bool symmetric;
 
  bool uzawa;
 
  int steps[2];
 
  int coarsen_target;
 
  int zcells;
 
  bool report;
 
  Smoother smoother;
 
  Preconditioner pre;
 
  MultiplierPreconditioner mortarpre;
 
  void parse(Opm::ParameterGroup& param)
  {
    std::string solver = param.getDefault<std::string>("linsolver_type","iterative");
    if (solver == "iterative")
      type = ITERATIVE;
    else
      type = DIRECT;
    restart = param.getDefault<int>("linsolver_restart", 1000);
    tol    = param.getDefault<double>("ltol",1.e-8);
    maxit   = param.getDefault<int>("linsolver_maxit", 10000);
    steps[0] = param.getDefault<int>("linsolver_presteps", 2);
    steps[1] = param.getDefault<int>("linsolver_poststeps", 2);
    coarsen_target = param.getDefault<int>("linsolver_coarsen", 5000);
    symmetric = param.getDefault<bool>("linsolver_symmetric", true);
    report = param.getDefault<bool>("linsolver_report", false);
    solver = param.getDefault<std::string>("linsolver_pre","");
 
    if (solver == "") // set default based on aspect ratio heuristic
      solver="heuristic";
 
    if (solver == "schwarz")
      pre = SCHWARZ;
    else if (solver == "fastamg")
      pre = FASTAMG;
    else if (solver == "twolevel")
      pre = TWOLEVEL;
    else if (solver == "heuristic")
      pre = UNDETERMINED;
    else
      pre = AMG;
 
    solver = param.getDefault<std::string>("linsolver_mortarpre","schur");
    if (solver == "simple")
      mortarpre = SIMPLE;
    else if (solver == "amg")
      mortarpre = SCHURAMG;
    else if (solver == "schwarz")
      mortarpre = SCHURSCHWARZ;
    else if (solver == "twolevel")
      mortarpre = SCHURTWOLEVEL;
    else
      mortarpre = SCHUR;
 
    uzawa = param.getDefault<bool>("linsolver_uzawa", false);
    zcells = param.getDefault<int>("linsolver_zcells", 2);
 
    solver = param.getDefault<std::string>("linsolver_smoother","ssor");
    if (solver == "schwarz")
      smoother = SMOOTH_SCHWARZ;
    else if (solver == "ilu")
      smoother = SMOOTH_ILU;
    else if (solver == "jacobi")
      smoother = SMOOTH_JACOBI;
    else {
      if (solver != "ssor")
        std::cerr << "WARNING: Invalid smoother specified, falling back to SSOR" << std::endl;
      smoother = SMOOTH_SSOR;
    }
 
    if (symmetric)
      steps[1] = steps[0];
  }
};
 
  template<class GridType, class PC>
class ElasticityUpscale
{
  public:
    static const int dim = GridType::dimension;
 
    typedef typename GridType::LeafGridView::ctype ctype;
 
    typedef Dune::FieldVector<double,dim> NodeValue;
 
    typedef typename GridType::LeafGridView::template Codim<1>::Geometry::GlobalCoordinate GlobalCoordinate;
 
    typedef typename GridType::LeafGridView::IndexSet LeafIndexSet;
 
    typedef typename GridType::LeafGridView::template Codim<0>::Iterator LeafIterator;
 
    typedef typename PC::type PCType;
 
    typedef std::shared_ptr<typename PC::type> PCPtr;
 
    ASMHandler<GridType> A;
 
    Vector u[6];
    Vector b[6];
 
    std::vector<double> volumeFractions;
    bool bySat; 
 
    double upscaledRho;
 
    ElasticityUpscale(const GridType& gv_, ctype tol_, ctype Escale_, 
                      const std::string& file, const std::string& rocklist,
                      bool verbose_)
      :  A(gv_), gv(gv_), tol(tol_), Escale(Escale_), E(gv_), verbose(verbose_),
         color(gv_)
    {
      if (rocklist.empty())
        loadMaterialsFromGrid(file);
      else
        loadMaterialsFromRocklist(file,rocklist);
    }
 
    void findBoundaries(double* min, double* max);
 
    void addMPC(Direction dir, int slavenode,
                const BoundaryGrid::Vertex& m);
 
    void periodicBCs(const double* min, const double* max);
 
    void periodicBCsMortar(const double* min,
                           const double* max, int n1, int n2,
                           int p1, int p2);
 
    void fixCorners(const double* min, const double* max);
 
    void assemble(int loadcase, bool matrix);
 
      template <int comp>
    void averageStress(Dune::FieldVector<ctype,comp>& sigma,
                       const Vector& u, int loadcase);
 
    void solve(int loadcase);
 
    void setupSolvers(const LinSolParams& params);
 
  private:
    typedef typename GridType::LeafGridView::template Codim<dim>::Iterator LeafVertexIterator;
 
    const GridType& gv;
 
    ctype tol;
 
    ctype Escale;
 
    ctype Emin;
 
    bool isOnPlane(Direction plane, GlobalCoordinate coord, ctype value);
 
    bool isOnLine(Direction dir, GlobalCoordinate coord, ctype x, ctype y);
 
    bool isOnPoint(GlobalCoordinate coord, GlobalCoordinate point);
 
    std::vector< std::shared_ptr<Material> > materials;
 
    std::vector<BoundaryGrid::Vertex> extractFace(Direction dir, ctype coord);
 
    enum SIDE {
      LEFT,
      RIGHT
    };
 
    BoundaryGrid extractMasterFace(Direction dir, ctype coord,
                                   SIDE side=LEFT, bool dc=false);
 
    void determineSideFaces(const double* min, const double* max);
 
    void periodicPlane(Direction planenormal, Direction dir,
                       const std::vector<BoundaryGrid::Vertex>& slavepointgrid,
                       const BoundaryGrid& mastergrid);
 
    void fixPoint(Direction dir, GlobalCoordinate coord,
                  const NodeValue& value = NodeValue(0));
 
    void fixLine(Direction dir, ctype x, ctype y,
                 const NodeValue& value = NodeValue(0));
 
    typedef std::map<int,BoundaryGrid::Vertex> SlaveGrid;
 
    int addBBlockMortar(const BoundaryGrid& b1,
                        const BoundaryGrid& interface, int dir,
                        int n1, int n2, int colofs);
 
    void assembleBBlockMortar(const BoundaryGrid& b1,
                              const BoundaryGrid& interface, int dir,
                              int n1, int n2, int colofs, double alpha=1.0);
 
    void loadMaterialsFromGrid(const std::string& file);
 
    void loadMaterialsFromRocklist(const std::string& file,
                                   const std::string& rocklist);
 
    std::vector<BoundaryGrid> master;
 
    std::vector< std::vector<BoundaryGrid::Vertex> > slave;
 
    AdjacencyPattern Bpatt;
    Matrix B;
 
    Matrix P; 
 
    typedef std::shared_ptr<Dune::InverseOperator<Vector, Vector> > SolverPtr;
    std::vector<SolverPtr> tsolver;
 
    std::shared_ptr<Operator> op;
 
    typedef MortarBlockEvaluator<Dune::Preconditioner<Vector, Vector> > SchurPreconditioner;
 
    typedef MortarBlockEvaluator<Dune::InverseOperator<Vector, Vector> > SchurEvaluator;
 
    std::shared_ptr<SchurEvaluator> op2;
 
    std::vector<PCPtr> upre;
 
    typedef Dune::InverseOperator2Preconditioner<LUSolver,
                                        Dune::SolverCategory::sequential> SeqLU;
    std::shared_ptr<SeqLU> lpre;
    std::shared_ptr<LUSolver> lprep;
 
    typedef std::shared_ptr< MortarSchurPre<PCType> > MortarAmgPtr;
    std::vector<MortarAmgPtr> tmpre;
 
    std::shared_ptr<MortarEvaluator> meval;
 
    Elasticity<GridType> E;
 
    bool verbose;
 
    MeshColorizer<GridType> color;
};
 
}} // namespace Opm, Elasticity
 
#include "elasticity_upscale_impl.hpp"
 
#endif