Loading...
Searching...
No Matches
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
asmhandler.hpp
Class handling finite element assembly.
boundarygrid.hh
Class describing 2D quadrilateral grids.
MeshColorizer
Generate a coloring of a mesh suitable for threaded assembly.
Definition meshcolorizer.hpp:26
Opm::Elasticity::ASMHandler
Class handling finite element assembly.
Definition asmhandler.hpp:35
Opm::Elasticity::BoundaryGrid::Vertex
A class describing a 2D vertex.
Definition boundarygrid.hh:62
Opm::Elasticity::BoundaryGrid
A class describing a quad grid.
Definition boundarygrid.hh:33
Opm::Elasticity::ElasticityUpscale::LeafIndexSet
GridType::LeafGridView::IndexSet LeafIndexSet
A set of indices.
Definition elasticity_upscale.hpp:206
Opm::Elasticity::ElasticityUpscale::bySat
bool bySat
Are volume fractions grouped by SATNUM?
Definition elasticity_upscale.hpp:228
Opm::Elasticity::ElasticityUpscale::GlobalCoordinate
GridType::LeafGridView::template Codim< 1 >::Geometry::GlobalCoordinate GlobalCoordinate
A global coordinate.
Definition elasticity_upscale.hpp:203
Opm::Elasticity::ElasticityUpscale::upscaledRho
double upscaledRho
Upscaled density.
Definition elasticity_upscale.hpp:231
Opm::Elasticity::ElasticityUpscale::solve
void solve(int loadcase)
Solve Au = b for u.
Opm::Elasticity::ElasticityUpscale::PCPtr
std::shared_ptr< typename PC::type > PCPtr
A pointer to our preconditioner.
Definition elasticity_upscale.hpp:215
Opm::Elasticity::ElasticityUpscale::u
Vector u[6]
The solution vectors.
Definition elasticity_upscale.hpp:221
Opm::Elasticity::ElasticityUpscale::LeafIterator
GridType::LeafGridView::template Codim< 0 >::Iterator LeafIterator
An iterator over grid cells.
Definition elasticity_upscale.hpp:209
Opm::Elasticity::ElasticityUpscale::A
ASMHandler< GridType > A
The linear operator.
Definition elasticity_upscale.hpp:218
Opm::Elasticity::ElasticityUpscale::periodicBCsMortar
void periodicBCsMortar(const double *min, const double *max, int n1, int n2, int p1, int p2)
Establish periodic boundaries using the mortar approach.
Opm::Elasticity::ElasticityUpscale::findBoundaries
void findBoundaries(double *min, double *max)
Find boundary coordinates.
Opm::Elasticity::ElasticityUpscale::setupSolvers
void setupSolvers(const LinSolParams &params)
Opm::Elasticity::ElasticityUpscale::addMPC
void addMPC(Direction dir, int slavenode, const BoundaryGrid::Vertex &m)
Add a MPC equation.
Opm::Elasticity::ElasticityUpscale::averageStress
void averageStress(Dune::FieldVector< ctype, comp > &sigma, const Vector &u, int loadcase)
Calculate the average stress vector for the given loadcase.
Opm::Elasticity::ElasticityUpscale::dim
static const int dim
Dimension of our grid.
Definition elasticity_upscale.hpp:194
Opm::Elasticity::ElasticityUpscale::ctype
GridType::LeafGridView::ctype ctype
A basic number.
Definition elasticity_upscale.hpp:197
Opm::Elasticity::ElasticityUpscale::periodicBCs
void periodicBCs(const double *min, const double *max)
Establish periodic boundaries using the MPC approach.
Opm::Elasticity::ElasticityUpscale::NodeValue
Dune::FieldVector< double, dim > NodeValue
A vectorial node value.
Definition elasticity_upscale.hpp:200
Opm::Elasticity::ElasticityUpscale::fixCorners
void fixCorners(const double *min, const double *max)
Fix corner nodes.
Opm::Elasticity::ElasticityUpscale::b
Vector b[6]
The load vectors.
Definition elasticity_upscale.hpp:223
Opm::Elasticity::ElasticityUpscale::assemble
void assemble(int loadcase, bool matrix)
Assemble (optionally) stiffness matrix A and load vector.
Opm::Elasticity::ElasticityUpscale::ElasticityUpscale
ElasticityUpscale(const GridType &gv_, ctype tol_, ctype Escale_, const std::string &file, const std::string &rocklist, bool verbose_)
Main constructor.
Definition elasticity_upscale.hpp:240
Opm::Elasticity::ElasticityUpscale::PCType
PC::type PCType
Our preconditioner type.
Definition elasticity_upscale.hpp:212
Opm::Elasticity::ElasticityUpscale::volumeFractions
std::vector< double > volumeFractions
Vector holding the volume fractions for materials (grouped by SATNUM).
Definition elasticity_upscale.hpp:226
Opm::Elasticity::Elasticity
Elasticity helper class.
Definition elasticity.hpp:22
Opm::Elasticity::MortarBlockEvaluator
This implements a operator evaluation for the schur mortar-block S = B^T*A^-1*B !
Definition mortar_schur.hpp:27
Opm::GridType
elasticity.hpp
Elasticity helper class.
elasticity_preconditioners.hpp
Preconditioners for elasticity upscaling.
Opm::Elasticity::Solver
Solver
An enumeration of available linear solver classes.
Definition elasticity_upscale.hpp:52
Opm::Elasticity::MultiplierPreconditioner
MultiplierPreconditioner
An enumeration of the available preconditioners for multiplier block.
Definition elasticity_upscale.hpp:66
Opm::Elasticity::SCHURAMG
@ SCHURAMG
schur + amg
Definition elasticity_upscale.hpp:69
Opm::Elasticity::SCHUR
@ SCHUR
schur + primary preconditioner
Definition elasticity_upscale.hpp:68
Opm::Elasticity::SCHURTWOLEVEL
@ SCHURTWOLEVEL
schur + twolevel
Definition elasticity_upscale.hpp:71
Opm::Elasticity::SIMPLE
@ SIMPLE
diagonal approximation of A
Definition elasticity_upscale.hpp:67
Opm::Elasticity::SCHURSCHWARZ
@ SCHURSCHWARZ
schur + schwarz+lu
Definition elasticity_upscale.hpp:70
Opm::Elasticity::Smoother
Smoother
Smoother used in the AMG.
Definition elasticity_upscale.hpp:75
elasticity_upscale_impl.hpp
Elasticity upscale class - template implementations.
logutils.hpp
Logging helper utilities.
materials.hh
Material properties.
matrixops.hpp
Helper class with some matrix operations.
Opm::Elasticity::Matrix
Dune::BCRSMatrix< Dune::FieldMatrix< double, 1, 1 > > Matrix
A sparse matrix holding our operator.
Definition matrixops.hpp:27
Opm::Elasticity::Vector
Dune::BlockVector< Dune::FieldVector< double, 1 > > Vector
A vector holding our RHS.
Definition matrixops.hpp:33
Opm::Elasticity::AdjacencyPattern
std::vector< std::set< int > > AdjacencyPattern
For storing matrix adjacency/sparsity patterns.
Definition matrixops.hpp:30
meshcolorizer.hpp
Mesh colorizer class.
mortar_evaluator.hpp
Linear operator for a Mortar block.
mortar_schur.hpp
Schur based linear operator for a Mortar block.
mortar_utils.hpp
Mortar helper class.
mpc.hh
Representation of multi-point constraint (MPC) equations.
Opm::Elasticity::Direction
Direction
An enum for specification of global coordinate directions.
Definition mpc.hh:24
Opm
Inverting small matrices.
Definition ImplicitAssembly.hpp:43
shapefunctions.hpp
Classes for shape functions. Loosely based on code in dune-grid-howto.
Opm::Elasticity::LinSolParams
Definition elasticity_upscale.hpp:82
Opm::Elasticity::LinSolParams::pre
Preconditioner pre
Preconditioner for elasticity block.
Definition elasticity_upscale.hpp:117
Opm::Elasticity::LinSolParams::report
bool report
Give a report at end of solution phase.
Definition elasticity_upscale.hpp:111
Opm::Elasticity::LinSolParams::maxit
int maxit
Max number of iterations.
Definition elasticity_upscale.hpp:90
Opm::Elasticity::LinSolParams::mortarpre
MultiplierPreconditioner mortarpre
Preconditioner for mortar block.
Definition elasticity_upscale.hpp:120
Opm::Elasticity::LinSolParams::steps
int steps[2]
The number of pre/post steps in the AMG.
Definition elasticity_upscale.hpp:102
Opm::Elasticity::LinSolParams::type
Solver type
The linear solver to employ.
Definition elasticity_upscale.hpp:84
Opm::Elasticity::LinSolParams::coarsen_target
int coarsen_target
Coarsening target in the AMG.
Definition elasticity_upscale.hpp:105
Opm::Elasticity::LinSolParams::parse
void parse(Opm::ParameterGroup &param)
Parse command line parameters.
Definition elasticity_upscale.hpp:124
Opm::Elasticity::LinSolParams::restart
int restart
Number of iterations in GMRES before restart.
Definition elasticity_upscale.hpp:87
Opm::Elasticity::LinSolParams::uzawa
bool uzawa
Use a Uzawa approach.
Definition elasticity_upscale.hpp:99
Opm::Elasticity::LinSolParams::tol
double tol
The tolerance for the iterative linear solver.
Definition elasticity_upscale.hpp:93
Opm::Elasticity::LinSolParams::zcells
int zcells
Number of cells in z to collapse in each cell.
Definition elasticity_upscale.hpp:108
Opm::Elasticity::LinSolParams::smoother
Smoother smoother
Smoother type used in the AMG.
Definition elasticity_upscale.hpp:114
Opm::Elasticity::LinSolParams::symmetric
bool symmetric
Use MINRES instead of GMRES (and thus symmetric preconditioning).
Definition elasticity_upscale.hpp:96
uzawa_solver.hpp
Uzawa scheme helper class.