12#ifndef ELASTICITY_UPSCALE_HPP_
13#define ELASTICITY_UPSCALE_HPP_
16#include <opm/common/utility/platform_dependent/disable_warnings.h>
18#include <dune/common/fmatrix.hh>
19#include <opm/common/utility/parameters/ParameterGroup.hpp>
20#include <dune/grid/common/mcmgmapper.hh>
21#include <dune/geometry/quadraturerules.hh>
22#include <dune/istl/ilu.hh>
23#include <dune/istl/solvers.hh>
24#include <dune/istl/preconditioners.hh>
25#include <opm/grid/CpGrid.hpp>
28#include <opm/common/utility/platform_dependent/reenable_warnings.h>
42#include <opm/elasticity/mortar_schur_precond.hpp>
45#include <opm/input/eclipse/Parser/Parser.hpp>
46#include <opm/input/eclipse/Deck/Deck.hpp>
124 void parse(Opm::ParameterGroup& param)
126 std::string solver = param.getDefault<std::string>(
"linsolver_type",
"iterative");
127 if (solver ==
"iterative")
131 restart = param.getDefault<
int>(
"linsolver_restart", 1000);
132 tol = param.getDefault<
double>(
"ltol",1.e-8);
133 maxit = param.getDefault<
int>(
"linsolver_maxit", 10000);
134 steps[0] = param.getDefault<
int>(
"linsolver_presteps", 2);
135 steps[1] = param.getDefault<
int>(
"linsolver_poststeps", 2);
136 coarsen_target = param.getDefault<
int>(
"linsolver_coarsen", 5000);
137 symmetric = param.getDefault<
bool>(
"linsolver_symmetric",
true);
138 report = param.getDefault<
bool>(
"linsolver_report",
false);
139 solver = param.getDefault<std::string>(
"linsolver_pre",
"");
144 if (solver ==
"schwarz")
146 else if (solver ==
"fastamg")
148 else if (solver ==
"twolevel")
150 else if (solver ==
"heuristic")
155 solver = param.getDefault<std::string>(
"linsolver_mortarpre",
"schur");
156 if (solver ==
"simple")
158 else if (solver ==
"amg")
160 else if (solver ==
"schwarz")
162 else if (solver ==
"twolevel")
167 uzawa = param.getDefault<
bool>(
"linsolver_uzawa",
false);
168 zcells = param.getDefault<
int>(
"linsolver_zcells", 2);
170 solver = param.getDefault<std::string>(
"linsolver_smoother",
"ssor");
171 if (solver ==
"schwarz")
173 else if (solver ==
"ilu")
175 else if (solver ==
"jacobi")
178 if (solver !=
"ssor")
179 std::cerr <<
"WARNING: Invalid smoother specified, falling back to SSOR" << std::endl;
189 template<
class Gr
idType,
class PC>
194 static const int dim = GridType::dimension;
197 typedef typename GridType::LeafGridView::ctype
ctype;
203 typedef typename GridType::LeafGridView::template Codim<1>::Geometry::GlobalCoordinate
GlobalCoordinate;
209 typedef typename GridType::LeafGridView::template Codim<0>::Iterator
LeafIterator;
215 typedef std::shared_ptr<typename PC::type>
PCPtr;
241 const std::string& file,
const std::string& rocklist,
243 :
A(gv_), gv(gv_), tol(tol_), Escale(Escale_), E(gv_), verbose(verbose_),
246 if (rocklist.empty())
247 loadMaterialsFromGrid(file);
249 loadMaterialsFromRocklist(file,rocklist);
277 const double* max,
int n1,
int n2,
296 const Vector&
u,
int loadcase);
307 typedef typename GridType::LeafGridView::template Codim<dim>::Iterator LeafVertexIterator;
340 std::vector< std::shared_ptr<Material> > materials;
346 std::vector<BoundaryGrid::Vertex> extractFace(
Direction dir,
ctype coord);
361 SIDE side=LEFT,
bool dc=
false);
366 void determineSideFaces(
const double* min,
const double* max);
374 const std::vector<BoundaryGrid::Vertex>& slavepointgrid,
393 typedef std::map<int,BoundaryGrid::Vertex> SlaveGrid;
405 int n1,
int n2,
int colofs);
417 int n1,
int n2,
int colofs,
double alpha=1.0);
421 void loadMaterialsFromGrid(
const std::string& file);
426 void loadMaterialsFromRocklist(
const std::string& file,
427 const std::string& rocklist);
430 std::vector<BoundaryGrid> master;
433 std::vector< std::vector<BoundaryGrid::Vertex> > slave;
442 typedef std::shared_ptr<Dune::InverseOperator<Vector, Vector> > SolverPtr;
443 std::vector<SolverPtr> tsolver;
446 std::shared_ptr<Operator> op;
455 std::shared_ptr<SchurEvaluator> op2;
458 std::vector<PCPtr> upre;
461 typedef Dune::InverseOperator2Preconditioner<LUSolver,
462 Dune::SolverCategory::sequential> SeqLU;
464 std::shared_ptr<SeqLU> lpre;
465 std::shared_ptr<LUSolver> lprep;
468 typedef std::shared_ptr< MortarSchurPre<PCType> > MortarAmgPtr;
469 std::vector<MortarAmgPtr> tmpre;
472 std::shared_ptr<MortarEvaluator> meval;
Class handling finite element assembly.
Class describing 2D quadrilateral grids.
Generate a coloring of a mesh suitable for threaded assembly.
Definition meshcolorizer.hpp:26
Class handling finite element assembly.
Definition asmhandler.hpp:35
A class describing a 2D vertex.
Definition boundarygrid.hh:62
A class describing a quad grid.
Definition boundarygrid.hh:33
The main driver class.
Definition elasticity_upscale.hpp:191
GridType::LeafGridView::IndexSet LeafIndexSet
A set of indices.
Definition elasticity_upscale.hpp:206
bool bySat
Are volume fractions grouped by SATNUM?
Definition elasticity_upscale.hpp:228
GridType::LeafGridView::template Codim< 1 >::Geometry::GlobalCoordinate GlobalCoordinate
A global coordinate.
Definition elasticity_upscale.hpp:203
double upscaledRho
Upscaled density.
Definition elasticity_upscale.hpp:231
void solve(int loadcase)
Solve Au = b for u.
std::shared_ptr< typename PC::type > PCPtr
A pointer to our preconditioner.
Definition elasticity_upscale.hpp:215
Vector u[6]
The solution vectors.
Definition elasticity_upscale.hpp:221
GridType::LeafGridView::template Codim< 0 >::Iterator LeafIterator
An iterator over grid cells.
Definition elasticity_upscale.hpp:209
ASMHandler< GridType > A
The linear operator.
Definition elasticity_upscale.hpp:218
void periodicBCsMortar(const double *min, const double *max, int n1, int n2, int p1, int p2)
Establish periodic boundaries using the mortar approach.
void findBoundaries(double *min, double *max)
Find boundary coordinates.
void setupSolvers(const LinSolParams ¶ms)
void addMPC(Direction dir, int slavenode, const BoundaryGrid::Vertex &m)
Add a MPC equation.
void averageStress(Dune::FieldVector< ctype, comp > &sigma, const Vector &u, int loadcase)
Calculate the average stress vector for the given loadcase.
static const int dim
Dimension of our grid.
Definition elasticity_upscale.hpp:194
GridType::LeafGridView::ctype ctype
A basic number.
Definition elasticity_upscale.hpp:197
void periodicBCs(const double *min, const double *max)
Establish periodic boundaries using the MPC approach.
Dune::FieldVector< double, dim > NodeValue
A vectorial node value.
Definition elasticity_upscale.hpp:200
void fixCorners(const double *min, const double *max)
Fix corner nodes.
Vector b[6]
The load vectors.
Definition elasticity_upscale.hpp:223
void assemble(int loadcase, bool matrix)
Assemble (optionally) stiffness matrix A and load vector.
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
PC::type PCType
Our preconditioner type.
Definition elasticity_upscale.hpp:212
std::vector< double > volumeFractions
Vector holding the volume fractions for materials (grouped by SATNUM)
Definition elasticity_upscale.hpp:226
Elasticity helper class.
Definition elasticity.hpp:22
This implements a operator evaluation for the schur mortar-block S = B^T*A^-1*B !
Definition mortar_schur.hpp:27
Preconditioners for elasticity upscaling.
Solver
An enumeration of available linear solver classes.
Definition elasticity_upscale.hpp:52
MultiplierPreconditioner
An enumeration of the available preconditioners for multiplier block.
Definition elasticity_upscale.hpp:66
@ SCHURAMG
schur + amg
Definition elasticity_upscale.hpp:69
@ SCHUR
schur + primary preconditioner
Definition elasticity_upscale.hpp:68
@ SCHURTWOLEVEL
schur + twolevel
Definition elasticity_upscale.hpp:71
@ SIMPLE
diagonal approximation of A
Definition elasticity_upscale.hpp:67
@ SCHURSCHWARZ
schur + schwarz+lu
Definition elasticity_upscale.hpp:70
Smoother
Smoother used in the AMG.
Definition elasticity_upscale.hpp:75
Elasticity upscale class - template implementations.
Logging helper utilities.
Helper class with some matrix operations.
Dune::BCRSMatrix< Dune::FieldMatrix< double, 1, 1 > > Matrix
A sparse matrix holding our operator.
Definition matrixops.hpp:27
Dune::BlockVector< Dune::FieldVector< double, 1 > > Vector
A vector holding our RHS.
Definition matrixops.hpp:33
std::vector< std::set< int > > AdjacencyPattern
For storing matrix adjacency/sparsity patterns.
Definition matrixops.hpp:30
Linear operator for a Mortar block.
Schur based linear operator for a Mortar block.
Representation of multi-point constraint (MPC) equations.
Direction
An enum for specification of global coordinate directions.
Definition mpc.hh:24
Inverting small matrices.
Definition ImplicitAssembly.hpp:43
Classes for shape functions. Loosely based on code in dune-grid-howto.
Definition elasticity_upscale.hpp:82
Preconditioner pre
Preconditioner for elasticity block.
Definition elasticity_upscale.hpp:117
bool report
Give a report at end of solution phase.
Definition elasticity_upscale.hpp:111
int maxit
Max number of iterations.
Definition elasticity_upscale.hpp:90
MultiplierPreconditioner mortarpre
Preconditioner for mortar block.
Definition elasticity_upscale.hpp:120
int steps[2]
The number of pre/post steps in the AMG.
Definition elasticity_upscale.hpp:102
Solver type
The linear solver to employ.
Definition elasticity_upscale.hpp:84
int coarsen_target
Coarsening target in the AMG.
Definition elasticity_upscale.hpp:105
void parse(Opm::ParameterGroup ¶m)
Parse command line parameters.
Definition elasticity_upscale.hpp:124
int restart
Number of iterations in GMRES before restart.
Definition elasticity_upscale.hpp:87
bool uzawa
Use a Uzawa approach.
Definition elasticity_upscale.hpp:99
double tol
The tolerance for the iterative linear solver.
Definition elasticity_upscale.hpp:93
int zcells
Number of cells in z to collapse in each cell.
Definition elasticity_upscale.hpp:108
Smoother smoother
Smoother type used in the AMG.
Definition elasticity_upscale.hpp:114
bool symmetric
Use MINRES instead of GMRES (and thus symmetric preconditioning)
Definition elasticity_upscale.hpp:96
Uzawa scheme helper class.