elasticity_upscale.hpp

Go to the documentation of this file.

//==============================================================================
//==============================================================================
#ifndef ELASTICITY_UPSCALE_HPP_
#define ELASTICITY_UPSCALE_HPP_

#include <dune/common/fmatrix.hh>
#include <opm/core/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 <dune/grid/CpGrid.hpp>
#include <opm/elasticity/shapefunctions.hpp>

#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/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/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::parameter::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 slave, 
                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 plane, Direction dir,
                       const std::vector<BoundaryGrid::Vertex>& slave,
                       const BoundaryGrid& master);

    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