All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
elasticity_upscale.hpp
Go to the documentation of this file.
1 //==============================================================================
11 //==============================================================================
12 #ifndef ELASTICITY_UPSCALE_HPP_
13 #define ELASTICITY_UPSCALE_HPP_
14 
15 
16 #include <opm/common/utility/platform_dependent/disable_warnings.h>
17 
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>
27 
28 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
29 
38 #include <opm/elasticity/mpc.hh>
42 #include <opm/elasticity/mortar_schur_precond.hpp>
44 
45 #include <opm/parser/eclipse/Parser/Parser.hpp>
46 #include <opm/parser/eclipse/Deck/Deck.hpp>
47 
48 namespace Opm {
49 namespace Elasticity {
50 
52 enum Solver {
53  DIRECT,
54  ITERATIVE
55 };
56 
57 enum Preconditioner {
58  AMG,
59  FASTAMG,
60  SCHWARZ,
61  TWOLEVEL,
62  UNDETERMINED
63 };
64 
66 enum MultiplierPreconditioner {
67  SIMPLE,
68  SCHUR,
69  SCHURAMG,
70  SCHURSCHWARZ,
71  SCHURTWOLEVEL
72 };
73 
75 enum Smoother {
76  SMOOTH_SSOR = 0,
77  SMOOTH_SCHWARZ = 1,
78  SMOOTH_JACOBI = 2,
79  SMOOTH_ILU = 4
80 };
81 
82 struct LinSolParams {
84  Solver type;
85 
87  int restart;
88 
90  int maxit;
91 
93  double tol;
94 
96  bool symmetric;
97 
99  bool uzawa;
100 
102  int steps[2];
103 
106 
108  int zcells;
109 
111  bool report;
112 
114  Smoother smoother;
115 
117  Preconditioner pre;
118 
120  MultiplierPreconditioner mortarpre;
121 
124  void parse(Opm::ParameterGroup& param)
125  {
126  std::string solver = param.getDefault<std::string>("linsolver_type","iterative");
127  if (solver == "iterative")
128  type = ITERATIVE;
129  else
130  type = DIRECT;
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","");
140 
141  if (solver == "") // set default based on aspect ratio heuristic
142  solver="heuristic";
143 
144  if (solver == "schwarz")
145  pre = SCHWARZ;
146  else if (solver == "fastamg")
147  pre = FASTAMG;
148  else if (solver == "twolevel")
149  pre = TWOLEVEL;
150  else if (solver == "heuristic")
151  pre = UNDETERMINED;
152  else
153  pre = AMG;
154 
155  solver = param.getDefault<std::string>("linsolver_mortarpre","schur");
156  if (solver == "simple")
157  mortarpre = SIMPLE;
158  else if (solver == "amg")
159  mortarpre = SCHURAMG;
160  else if (solver == "schwarz")
161  mortarpre = SCHURSCHWARZ;
162  else if (solver == "twolevel")
163  mortarpre = SCHURTWOLEVEL;
164  else
165  mortarpre = SCHUR;
166 
167  uzawa = param.getDefault<bool>("linsolver_uzawa", false);
168  zcells = param.getDefault<int>("linsolver_zcells", 2);
169 
170  solver = param.getDefault<std::string>("linsolver_smoother","ssor");
171  if (solver == "schwarz")
172  smoother = SMOOTH_SCHWARZ;
173  else if (solver == "ilu")
174  smoother = SMOOTH_ILU;
175  else if (solver == "jacobi")
176  smoother = SMOOTH_JACOBI;
177  else {
178  if (solver != "ssor")
179  std::cerr << "WARNING: Invalid smoother specified, falling back to SSOR" << std::endl;
180  smoother = SMOOTH_SSOR;
181  }
182 
183  if (symmetric)
184  steps[1] = steps[0];
185  }
186 };
187 
189  template<class GridType, class PC>
191 {
192  public:
194  static const int dim = GridType::dimension;
195 
197  typedef typename GridType::LeafGridView::ctype ctype;
198 
200  typedef Dune::FieldVector<double,dim> NodeValue;
201 
203  typedef typename GridType::LeafGridView::template Codim<1>::Geometry::GlobalCoordinate GlobalCoordinate;
204 
206  typedef typename GridType::LeafGridView::IndexSet LeafIndexSet;
207 
209  typedef typename GridType::LeafGridView::template Codim<0>::Iterator LeafIterator;
210 
212  typedef typename PC::type PCType;
213 
215  typedef std::shared_ptr<typename PC::type> PCPtr;
216 
219 
221  Vector u[6];
223  Vector b[6];
224 
226  std::vector<double> volumeFractions;
228  bool bySat;
229 
231  double upscaledRho;
232 
240  ElasticityUpscale(const GridType& gv_, ctype tol_, ctype Escale_,
241  const std::string& file, const std::string& rocklist,
242  bool verbose_)
243  : A(gv_), gv(gv_), tol(tol_), Escale(Escale_), E(gv_), verbose(verbose_),
244  color(gv_)
245  {
246  if (rocklist.empty())
247  loadMaterialsFromGrid(file);
248  else
249  loadMaterialsFromRocklist(file,rocklist);
250  }
251 
255  void findBoundaries(double* min, double* max);
256 
261  void addMPC(Direction dir, int slavenode,
262  const BoundaryGrid::Vertex& m);
263 
267  void periodicBCs(const double* min, const double* max);
268 
276  void periodicBCsMortar(const double* min,
277  const double* max, int n1, int n2,
278  int p1, int p2);
279 
283  void fixCorners(const double* min, const double* max);
284 
288  void assemble(int loadcase, bool matrix);
289 
294  template <int comp>
295  void averageStress(Dune::FieldVector<ctype,comp>& sigma,
296  const Vector& u, int loadcase);
297 
300  void solve(int loadcase);
301 
303  void setupSolvers(const LinSolParams& params);
304 
305  private:
307  typedef typename GridType::LeafGridView::template Codim<dim>::Iterator LeafVertexIterator;
308 
310  const GridType& gv;
311 
313  ctype tol;
314 
316  ctype Escale;
317 
319  ctype Emin;
320 
325  bool isOnPlane(Direction plane, GlobalCoordinate coord, ctype value);
326 
332  bool isOnLine(Direction dir, GlobalCoordinate coord, ctype x, ctype y);
333 
337  bool isOnPoint(GlobalCoordinate coord, GlobalCoordinate point);
338 
340  std::vector< std::shared_ptr<Material> > materials;
341 
346  std::vector<BoundaryGrid::Vertex> extractFace(Direction dir, ctype coord);
347 
349  enum SIDE {
350  LEFT,
351  RIGHT
352  };
353 
360  BoundaryGrid extractMasterFace(Direction dir, ctype coord,
361  SIDE side=LEFT, bool dc=false);
362 
366  void determineSideFaces(const double* min, const double* max);
367 
373  void periodicPlane(Direction planenormal, Direction dir,
374  const std::vector<BoundaryGrid::Vertex>& slavepointgrid,
375  const BoundaryGrid& mastergrid);
376 
381  void fixPoint(Direction dir, GlobalCoordinate coord,
382  const NodeValue& value = NodeValue(0));
383 
389  void fixLine(Direction dir, ctype x, ctype y,
390  const NodeValue& value = NodeValue(0));
391 
393  typedef std::map<int,BoundaryGrid::Vertex> SlaveGrid;
394 
403  int addBBlockMortar(const BoundaryGrid& b1,
404  const BoundaryGrid& interface, int dir,
405  int n1, int n2, int colofs);
406 
415  void assembleBBlockMortar(const BoundaryGrid& b1,
416  const BoundaryGrid& interface, int dir,
417  int n1, int n2, int colofs, double alpha=1.0);
418 
421  void loadMaterialsFromGrid(const std::string& file);
422 
426  void loadMaterialsFromRocklist(const std::string& file,
427  const std::string& rocklist);
428 
430  std::vector<BoundaryGrid> master;
431 
433  std::vector< std::vector<BoundaryGrid::Vertex> > slave;
434 
436  AdjacencyPattern Bpatt;
437  Matrix B;
438 
439  Matrix P;
440 
442  typedef std::shared_ptr<Dune::InverseOperator<Vector, Vector> > SolverPtr;
443  std::vector<SolverPtr> tsolver;
444 
446  std::shared_ptr<Operator> op;
447 
449  typedef MortarBlockEvaluator<Dune::Preconditioner<Vector, Vector> > SchurPreconditioner;
450 
452  typedef MortarBlockEvaluator<Dune::InverseOperator<Vector, Vector> > SchurEvaluator;
453 
455  std::shared_ptr<SchurEvaluator> op2;
456 
458  std::vector<PCPtr> upre;
459 
461  typedef Dune::InverseOperator2Preconditioner<LUSolver,
462  Dune::SolverCategory::sequential> SeqLU;
464  std::shared_ptr<SeqLU> lpre;
465  std::shared_ptr<LUSolver> lprep;
466 
468  typedef std::shared_ptr< MortarSchurPre<PCType> > MortarAmgPtr;
469  std::vector<MortarAmgPtr> tmpre;
470 
472  std::shared_ptr<MortarEvaluator> meval;
473 
475  Elasticity<GridType> E;
476 
478  bool verbose;
479 
482 };
483 
484 }} // namespace Opm, Elasticity
485 
487 
488 #endif
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.
Material properties.
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.
Mortar 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
Mesh colorizer class.
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.
Elasticity helper class.
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 &params)
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 &param)
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