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/core/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 <dune/grid/CpGrid.hpp>
28 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
42 #include <opm/elasticity/mortar_schur_precond.hpp>
45 #include <opm/parser/eclipse/Parser/Parser.hpp>
46 #include <opm/parser/eclipse/Deck/Deck.hpp>
49 namespace Elasticity {
66 enum MultiplierPreconditioner {
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);
261 void addMPC(Direction dir,
int slavenode,
267 void periodicBCs(
const double* min,
const double* max);
277 const double* max,
int n1,
int n2,
283 void fixCorners(
const double* min,
const double* max);
288 void assemble(
int loadcase,
bool matrix);
296 const Vector&
u,
int loadcase);
300 void solve(
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);
360 BoundaryGrid extractMasterFace(Direction dir,
ctype coord,
361 SIDE side=LEFT,
bool dc=
false);
366 void determineSideFaces(
const double* min,
const double* max);
373 void periodicPlane(Direction planenormal, Direction dir,
374 const std::vector<BoundaryGrid::Vertex>& slavepointgrid,
375 const BoundaryGrid& mastergrid);
393 typedef std::map<int,BoundaryGrid::Vertex> SlaveGrid;
403 int addBBlockMortar(
const BoundaryGrid& b1,
404 const BoundaryGrid& interface,
int dir,
405 int n1,
int n2,
int colofs);
415 void assembleBBlockMortar(
const BoundaryGrid& b1,
416 const BoundaryGrid& interface,
int dir,
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;
436 AdjacencyPattern Bpatt;
442 typedef std::shared_ptr<Dune::InverseOperator<Vector, Vector> > SolverPtr;
443 std::vector<SolverPtr> tsolver;
446 std::shared_ptr<Operator> op;
449 typedef MortarBlockEvaluator<Dune::Preconditioner<Vector, Vector> > SchurPreconditioner;
452 typedef MortarBlockEvaluator<Dune::InverseOperator<Vector, Vector> > SchurEvaluator;
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;
475 Elasticity<GridType> E;
bool report
Give a report at end of solution phase.
Definition: elasticity_upscale.hpp:111
void averageStress(Dune::FieldVector< ctype, comp > &sigma, const Vector &u, int loadcase)
Calculate the average stress vector for the given loadcase.
int zcells
Number of cells in z to collapse in each cell.
Definition: elasticity_upscale.hpp:108
bool uzawa
Use a Uzawa approach.
Definition: elasticity_upscale.hpp:99
static const int dim
Dimension of our grid.
Definition: elasticity_upscale.hpp:194
Helper class with some matrix operations.
void periodicBCs(const double *min, const double *max)
Establish periodic boundaries using the MPC approach.
Preconditioners for elasticity upscaling.
Vector b[6]
The load vectors.
Definition: elasticity_upscale.hpp:223
Solver type
The linear solver to employ.
Definition: elasticity_upscale.hpp:84
Uzawa scheme helper class.
Representation of multi-point constraint (MPC) equations.
void solve(int loadcase)
Solve Au = b for u.
PC::type PCType
Our preconditioner type.
Definition: elasticity_upscale.hpp:212
void assemble(int loadcase, bool matrix)
Assemble (optionally) stiffness matrix A and load vector.
Class describing 2D quadrilateral grids.
ASMHandler< GridType > A
The linear operator.
Definition: elasticity_upscale.hpp:218
The main driver class.
Definition: elasticity_upscale.hpp:190
Generate a coloring of a mesh suitable for threaded assembly.
Definition: meshcolorizer.hpp:26
bool bySat
Are volume fractions grouped by SATNUM?
Definition: elasticity_upscale.hpp:228
GridType::LeafGridView::template Codim< 0 >::Iterator LeafIterator
An iterator over grid cells.
Definition: elasticity_upscale.hpp:209
Class handling finite element assembly.
int restart
Number of iterations in GMRES before restart.
Definition: elasticity_upscale.hpp:87
Preconditioner pre
Preconditioner for elasticity block.
Definition: elasticity_upscale.hpp:117
Linear operator for a Mortar block.
void periodicBCsMortar(const double *min, const double *max, int n1, int n2, int p1, int p2)
Establish periodic boundaries using the mortar approach.
A class describing a 2D vertex.
Definition: boundarygrid.hh:64
void findBoundaries(double *min, double *max)
Find boundary coordinates.
int maxit
Max number of iterations.
Definition: elasticity_upscale.hpp:90
void fixCorners(const double *min, const double *max)
Fix corner nodes.
Definition: elasticity_upscale.hpp:82
Class handling finite element assembly.
Definition: asmhandler.hpp:35
GridType::LeafGridView::ctype ctype
A basic number.
Definition: elasticity_upscale.hpp:197
Smoother smoother
Smoother type used in the AMG.
Definition: elasticity_upscale.hpp:114
int coarsen_target
Coarsening target in the AMG.
Definition: elasticity_upscale.hpp:105
void addMPC(Direction dir, int slavenode, const BoundaryGrid::Vertex &m)
Add a MPC equation.
bool symmetric
Use MINRES instead of GMRES (and thus symmetric preconditioning)
Definition: elasticity_upscale.hpp:96
std::vector< double > volumeFractions
Vector holding the volume fractions for materials (grouped by SATNUM)
Definition: elasticity_upscale.hpp:226
GridType::LeafGridView::IndexSet LeafIndexSet
A set of indices.
Definition: elasticity_upscale.hpp:206
Schur based linear operator for a Mortar block.
MultiplierPreconditioner mortarpre
Preconditioner for mortar block.
Definition: elasticity_upscale.hpp:120
Classes for shape functions. Loosely based on code in dune-grid-howto.
GridType::LeafGridView::template Codim< 1 >::Geometry::GlobalCoordinate GlobalCoordinate
A global coordinate.
Definition: elasticity_upscale.hpp:203
void setupSolvers(const LinSolParams ¶ms)
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
std::shared_ptr< typename PC::type > PCPtr
A pointer to our preconditioner.
Definition: elasticity_upscale.hpp:215
void parse(Opm::ParameterGroup ¶m)
Parse command line parameters.
Definition: elasticity_upscale.hpp:124
Vector u[6]
The solution vectors.
Definition: elasticity_upscale.hpp:221
Logging helper utilities.
Dune::FieldVector< double, dim > NodeValue
A vectorial node value.
Definition: elasticity_upscale.hpp:200
double tol
The tolerance for the iterative linear solver.
Definition: elasticity_upscale.hpp:93
Elasticity upscale class - template implementations.
int steps[2]
The number of pre/post steps in the AMG.
Definition: elasticity_upscale.hpp:102
double upscaledRho
Upscaled density.
Definition: elasticity_upscale.hpp:231