OilfieldToolsNet/opm-upscaling
4
0
Fork 0
mirror of https://github.com/OPM/opm-upscaling.git synced 2025-02-25 18:45:23 -06:00
Code Issues Packages Projects Releases Wiki Activity

added: improve the linear algebra

Browse Source 
- change linsolver_type=cg to linsolver_type=iterative
- rename linsolver_type=slu to linsolver_type=direct
- adds a schur based preconditioner for the mortar system
- can use a symmetric iterative scheme (minres) or an asymmetric (gmres)
- allows employing a uzawa approach (mostly for debugging)

the elasticity block and (through schur-decomposition) multiplier
blocks are both preconditioned using an AMG. a diagonal preconditioner
for the multipliers are also available.
This commit is contained in:
Arne Morten Kvarving 2013-01-02 14:02:09 +01:00
parent a7dc718f2c
commit eaca1c97ad
15 changed files with 1146 additions and 536 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

47
dune/elasticity/applier.hpp Normal file
View File

@ -0,0 +1,47 @@
//==============================================================================
//!
//! \file applier.hpp
//!
//! \date Dec 22 2012
//!
//! \author Arne Morten Kvarving / SINTEF
//!
//! \brief Helper class for abstracting differences between inverse operators
//! and preconditioners in DUNE
//!
//==============================================================================
#pragma once
namespace Opm {
namespace Elasticity {
template<class T>
struct OperatorApplier
{
OperatorApplier(T& t) : A(t)
{
}
void apply(Vector& v, Vector& d);
T& A;
};
typedef OperatorApplier<Dune::InverseOperator<Vector, Vector> > InverseApplier;
typedef OperatorApplier<Dune::Preconditioner<Vector, Vector> > PreApplier;
template<>
void InverseApplier::apply(Vector& v, Vector& d)
{
Dune::InverseOperatorResult r;
A.apply(v, d, r);
}
template<>
void PreApplier::apply(Vector& v, Vector& d)
{
A.apply(v, d);
}
}
}

8
dune/elasticity/asmhandler.hpp
View File

@ -15,6 +15,7 @@
#include <dune/common/fmatrix.hh>
#include <dune/istl/bvector.hh>
#include <dune/common/fvector.hh>
#include <dune/elasticity/logutils.hpp>
#include <dune/elasticity/mpc.hh>
#include <dune/elasticity/matrixops.hpp>
@ -139,6 +140,13 @@ class ASMHandler {
//! \brief Print the current load vector
void printLoadVector() const;
//! \brief Access current adjacency pattern
//! \details Can be used to add extra entries, such as other blocks
AdjacencyPattern& getAdjacencyPattern()
{
return adjacencyPattern;
}
protected:
//! \brief Resolve chained MPCs
void resolveMPCChains()

22
dune/elasticity/asmhandler_impl.hpp
View File

@ -21,14 +21,16 @@ void ASMHandler<GridType>::initForAssembly()
A.setBuildMode(Matrix::random);
A.endrowsizes();
MatrixOps::fromAdjacency(A,adjacencyPattern,maxeqn,maxeqn);
b.resize(maxeqn);
MatrixOps::fromAdjacency(A,adjacencyPattern,
adjacencyPattern.size(),adjacencyPattern.size());
b.resize(adjacencyPattern.size());
b = 0;
adjacencyPattern.clear();
// print some information
std::cout << "Number of nodes: " << gv.size(dim) << std::endl;
std::cout << "Number of elements: " << gv.size(0) << std::endl;
std::cout << "Number of constraints: " << mpcs.size() << std::endl;
std::cout << "\tNumber of nodes: " << gv.size(dim) << std::endl;
std::cout << "\tNumber of elements: " << gv.size(0) << std::endl;
std::cout << "\tNumber of constraints: " << mpcs.size() << std::endl;
int fixedDofs=0;
for (fixIt it = fixedNodes.begin(); it != fixedNodes.end(); ++it) {
if (it->second.first & X)
@ -38,7 +40,7 @@ void ASMHandler<GridType>::initForAssembly()
if (it->second.first & Z)
fixedDofs++;
}
std::cout << "Number of fixed dofs: " << fixedDofs << std::endl;
std::cout << "\tNumber of fixed dofs: " << fixedDofs << std::endl;
}
template<class GridType>
@ -339,7 +341,7 @@ void ASMHandler<GridType>::preprocess()
}
}
}
std::cout << "number of equations: " << maxeqn << std::endl;
std::cout << "\tnumber of equations: " << maxeqn << std::endl;
}
template<class GridType>
@ -370,12 +372,15 @@ void ASMHandler<GridType>::nodeAdjacency(const LeafIterator& it,
void ASMHandler<GridType>::determineAdjacencyPattern()
{
adjacencyPattern.resize(maxeqn);
std::cout << "\tsetting up sparsity pattern..." << std::endl;
LoggerHelper help(gv.size(0), 5, 50000);
const LeafIndexSet& set = gv.leafView().indexSet();
LeafIterator itend = gv.leafView().template end<0>();
// iterate over cells
for (LeafIterator it = gv.leafView().template begin<0>(); it != itend; ++it) {
int cell=0;
for (LeafIterator it = gv.leafView().template begin<0>(); it != itend; ++it, ++cell) {
Dune::GeometryType gt = it->type();
const Dune::template GenericReferenceElement<double,dim>& ref =
Dune::GenericReferenceElements<double,dim>::general(gt);
@ -395,5 +400,6 @@ void ASMHandler<GridType>::determineAdjacencyPattern()
nodeAdjacency(it,vertexsize,meqn[indexi*dim+k]);
}
}
help.log(cell, "\t\t... still processing ... cell ");
}
}

250
dune/elasticity/elasticity_upscale.hpp
View File

@ -17,30 +17,101 @@
#include <dune/geometry/quadraturerules.hh>
#include <dune/istl/ilu.hh>
#include <dune/istl/solvers.hh>
#if HAVE_SUPERLU
#include <dune/istl/superlu.hh>
#endif
#include <dune/istl/preconditioners.hh>
#include <dune/istl/overlappingschwarz.hh>
#include <dune/elasticity/seqlu.hpp>
#include <dune/grid/CpGrid.hpp>
#include <dune/elasticity/shapefunctions.hpp>
#include <dune/elasticity/asmhandler.hpp>
#include <dune/elasticity/boundarygrid.hh>
#include <dune/elasticity/elasticity.hpp>
#include <dune/elasticity/logutils.hpp>
#include <dune/elasticity/materials.hh>
#include <dune/elasticity/mpc.hh>
#include <dune/elasticity/mortar_schur.hpp>
#include <dune/elasticity/mortar_utils.hpp>
#include <dune/elasticity/mortar_evaluator.hpp>
#include <dune/elasticity/mortar_schur_precond.hpp>
#include <dune/elasticity/uzawa_solver.hpp>
#include <dune/istl/paamg/amg.hh>
namespace Opm {
namespace Elasticity {
//! \brief An enumeration of available linear solvers
//! \brief An enumeration of available linear solver classes
enum Solver {
SLU,
CG
DIRECT,
ITERATIVE
};
typedef Dune::SeqOverlappingSchwarz<Matrix,Vector,Dune::SymmetricMultiplicativeSchwarzMode> OverlappingSchwarz;
//! \brief An enumeration of the available preconditioners
enum Preconditioner {
SCHURAMG,
SCHURDIAG
};
struct LinSolParams {
//! \brief The linear solver to employ
Opm::Elasticity::Solver type;
//! \brief Number of iterations in GMRES before restart
int restart;
//! \brief Max number of iterations
int maxit;
//! \brief The tolerance for the iterative linear solver
double tol;
//! \brief Use MINRES instead of GMRES (and thus symmetric preconditioning)
bool symmetric;
//! \brief Use a Uzawa approach
bool uzawa;
//! \brief The number of pre/post steps in the AMG
int steps[2];
//! \brief Coarsening target in the AMG
int coarsen_target;
//! \brief Number of cells in z to collapse in each cell
int zcells;
//! \brief Preconditioner for mortar block
Opm::Elasticity::Preconditioner mortarpre;
//! \brief Parse command line parameters
//! \param[in] param The parameter group to parse
void parse(Opm::parameter::ParameterGroup& param)
{
std::string solver = param.getDefault<std::string>("linsolver_type","iterative");
if (solver == "iterative")
type = Opm::Elasticity::ITERATIVE;
else
type = Opm::Elasticity::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);
solver = param.getDefault<std::string>("linsolver_mortarpre","schuramg");
if (solver == "schuramg")
mortarpre = Opm::Elasticity::SCHURAMG;
else
mortarpre = Opm::Elasticity::SCHURDIAG;
uzawa = param.getDefault<bool>("linsolver_uzawa", false);
zcells = param.getDefault<int>("linsolver_zcells", 2);
if (symmetric)
steps[1] = steps[0];
}
};
//! \brief The main driver class
template<class GridType>
@ -81,22 +152,24 @@ class ElasticityUpscale
//! \param[in] tol_ The tolerance to use when deciding whether or not a coordinate falls on a plane/line/point. \sa tol
//! \param[in] Escale_ A scale value for E-moduluses to avoid numerical issues
//! \param[in] file The eclipse grid file
//! \param[in] rocklist If true, file is a rocklist
//! \param[in] rocklist If not blank, file is a rocklist
//! \param[in] verbose If true, give verbose output
ElasticityUpscale(const GridType& gv_, ctype tol_, ctype Escale_,
const std::string& file, const std::string& rocklist)
: gv(gv_), tol(tol_), A(gv_), E(gv_), Escale(Escale_)
const std::string& file, const std::string& rocklist,
bool verbose_)
: A(gv_), gv(gv_), tol(tol_), Escale(Escale_), E(gv_), verbose(verbose_)
{
if (rocklist.empty())
loadMaterialsFromGrid(file);
else
loadMaterialsFromRocklist(file,rocklist);
#if HAVE_SUPERLU
slu = 0;
#endif
cgsolver = 0;
solver = 0;
op = 0;
ilu = 0;
ovl = 0;
mpre = 0;
upre = 0;
op2 = 0;
meval = 0;
lpre = 0;
}
//! \brief The destructor
@ -111,13 +184,13 @@ class ElasticityUpscale
it != itend; ++it)
delete *it;
#if HAVE_SUPERLU
delete slu;
#endif
delete cgsolver;
delete ilu;
delete solver;
delete op;
delete ovl;
delete op2;
delete meval;
delete upre;
delete lpre;
delete mpre;
}
//! \brief Find boundary coordinates
@ -137,15 +210,6 @@ class ElasticityUpscale
//! \param[in] max The maximum coordinates of the grid
void periodicBCs(const double* min, const double* max);
//! \brief Establish periodic boundaries using the LLM approach
//! \param[in] min The minimum coordinates of the grid
//! \param[in] max The maximum coordinates of the grid
//! \param[in] n1 The number of elements on the interface grid in the X direction
//! \param[in] n2 The number of elements on the interface grid in the Y direction
void periodicBCsLLM(const double* min,
const double* max,
int n1, int n2);
//! \brief Establish periodic boundaries using the mortar approach
//! \param[in] min The minimum coordinates of the grid
//! \param[in] max The maximum coordinates of the grid
@ -176,15 +240,16 @@ class ElasticityUpscale
const Vector& u, int loadcase);
//! \brief Solve Au = b for u
//! \param[in] solver The linear equation solver to employ
//! \param[in] tol The tolerance for iterative solvers
void solve(Solver solver, double tol, int loadcase);
//! \param[in] loadcase The load case to solve
void solve(int loadcase);
//! \param[in] params The linear solver parameters
void setupSolvers(const LinSolParams& params);
void setupSolvers(Solver solver);
private:
//! \brief An iterator over grid vertices
typedef typename GridType::LeafGridView::template Codim<dim>::Iterator LeafVertexIterator;
//! \brief A reference to our grid
const GridType& gv;
@ -271,26 +336,29 @@ class ElasticityUpscale
//! \brief A point grid
typedef std::map<int,BoundaryGrid::Vertex> SlaveGrid;
//! \brief Find the B matrix associated with a particular slave grid (LLM)
//! \param[in] slave The slave grid to extract
//! \returns The B matrix associated with the given slave grid
Matrix findBMatrixLLM(const SlaveGrid& slave);
//! \brief Find the L matrix associated with a particular slave grid (LLM)
//! \param[in] slave The slave grid we want to couple
//! \param[in] master The interface grid we want to couple to
//! \returns The L matrix describing the coupling between slave and master
Matrix findLMatrixLLM(const SlaveGrid& slave, const BoundaryGrid& master);
//! \brief Find the L matrix associated with mortar couplings
//! \brief Add a block to the B matrix associated with mortar couplings
//! \param[in] b1 The primal boundary to match
//! \param[in] interface The interface/multiplier grid
//! \param[in] dir The normal direction on the boundary (0/1)
//! \param[in] n1 The multipler order in the first direction
//! \param[in] n2 The multipler order in the second direction
Matrix findLMatrixMortar(const BoundaryGrid& b1,
const BoundaryGrid& interface, int dir,
int n1, int n2);
//! \param[in] colofs The column offset (multiplier number)
//! \returns Number of multipliers DOFs added
int addBBlockMortar(const BoundaryGrid& b1,
const BoundaryGrid& interface, int dir,
int n1, int n2, int colofs);
//! \brief Assemble part of the B block associated with mortar couplings
//! \param[in] b1 The primal boundary to match
//! \param[in] interface The interface/multiplier grid
//! \param[in] dir The normal direction on the boundary (0/1)
//! \param[in] n1 The multipler order in the first direction
//! \param[in] n2 The multipler order in the second direction
//! \param[in] colofs The column offset (first multiplier unknown)
//! \param[in] alpha Scaling for matrix elements
void assembleBBlockMortar(const BoundaryGrid& b1,
const BoundaryGrid& interface, int dir,
int n1, int n2, int colofs, double alpha=1.0);
//! \brief This function loads and maps materials to active grid cells
//! \param[in] file The eclipse grid to read materials from
@ -302,9 +370,12 @@ class ElasticityUpscale
void loadMaterialsFromRocklist(const std::string& file,
const std::string& rocklist);
//! \brief Setup an overlapping schwarz preconditioner with one
//! subdomain per pillar.
void setupPreconditioner();
//! \brief Setup AMG preconditioner
//! \param[in] pre The number of pre-smoothing steps
//! \param[in] post The number of post-smoothing steps
//! \param[in] target The coarsening target
//! \param[in] zcells The wanted number of cells to collapse in z per level
void setupAMG(int pre, int post, int target, int zcells);
//! \brief Master grids
std::vector<BoundaryGrid> master;
@ -312,30 +383,63 @@ class ElasticityUpscale
//! \brief Slave point grids
std::vector< std::vector<BoundaryGrid::Vertex> > slave;
//! \brief Vector of matrices holding boolean coupling matrices (LLM)
std::vector<Matrix> B;
//! \brief Vector of matrices holding Lagrangian multipliers
std::vector<Matrix> L;
//! \brief Lagrangian multiplier block
AdjacencyPattern Bpatt;
Matrix B;
#if HAVE_SUPERLU
//! \brief SuperLU solver
Dune::SuperLU<Matrix>* slu;
#endif
//! \brief CG solver
Dune::CGSolver<Vector>* cgsolver;
//! \brief The linear operator
Dune::MatrixAdapter<Matrix,Vector,Vector>* op;
Matrix P; //!< Preconditioner for multiplier block
//! \brief ILU preconditioner
Dune::SeqILU0<Matrix,Vector,Vector>* ilu;
//! \brief Overlapping Schwarz preconditioner
OverlappingSchwarz* ovl;
//! \brief Linear solver
Dune::InverseOperator<Vector, Vector>* solver;
//! \brief The smoother used in the AMG
typedef Dune::SeqSSOR<Matrix, Vector, Vector> Smoother;
//! \brief The coupling metric used in the AMG
typedef Dune::Amg::RowSum CouplingMetric;
//! \brief The coupling criterion used in the AMG
typedef Dune::Amg::SymmetricCriterion<Matrix, CouplingMetric> CritBase;
//! \brief The coarsening criterion used in the AMG
typedef Dune::Amg::CoarsenCriterion<CritBase> Criterion;
//! \brief A linear operator
typedef Dune::MatrixAdapter<Matrix,Vector,Vector> Operator;
//! \brief A preconditioner for an elasticity operator
typedef Dune::Amg::AMG<Operator, Vector, Smoother> ElasticityAMG;
//! \brief Matrix adaptor for the elasticity block
Operator* op;
//! \brief Preconditioner for multiplier block
typedef MortarBlockEvaluator<Dune::Preconditioner<Vector, Vector> > SchurPreconditioner;
//! \brief Evaluator for multiplier block
typedef MortarBlockEvaluator<Dune::InverseOperator<Vector, Vector> > SchurEvaluator;
//! \brief Outer evaluator, used with uzawa
SchurEvaluator* op2;
//! \brief The preconditioner for the elasticity operator
ElasticityAMG* upre;
//! \brief The preconditioner for the multiplier block (used with uzawa)
SeqLU<Matrix, Vector, Vector>* lpre;
//! \brief Preconditioner for the Mortar system
MortarSchurPre<ElasticityAMG>* mpre;
//! \brief Evaluator for the Mortar system
MortarEvaluator* meval;
//! \brief Elasticity helper class
Elasticity<GridType> E;
};
//! \brief Verbose output
bool verbose;
};
#include "elasticity_upscale_impl.hpp"
}

655
dune/elasticity/elasticity_upscale_impl.hpp
View File

@ -45,13 +45,7 @@ BoundaryGrid ElasticityUpscale<GridType>::extractMasterFace(Direction dir,
{2,3,6,7},
{4,5,6,7}};
const LeafIndexSet& set = gv.leafView().indexSet();
const LeafVertexIterator itend = gv.leafView().template end<dim>();
// make a mapper for codim dim entities in the leaf grid
Dune::LeafMultipleCodimMultipleGeomTypeMapper<GridType,
Dune::MCMGVertexLayout> mapper(gv);
LeafVertexIterator start=gv.leafView().template begin<dim>();
LeafIterator cellend = gv.leafView().template end<0>();
int c = 0;
int i = log2(dir);
BoundaryGrid result;
@ -59,33 +53,25 @@ BoundaryGrid ElasticityUpscale<GridType>::extractMasterFace(Direction dir,
// vertex. we then split this up into pillars for easy processing later
std::map<double, std::vector<BoundaryGrid::Quad> > nodeMap;
for (LeafIterator cell = gv.leafView().template begin<0>();
cell != cellend; ++cell, ++c) {
cell != gv.leafView().template end<0>(); ++cell, ++c) {
std::vector<BoundaryGrid::Vertex> verts;
int idx;
int idx=0;
if (side == LEFT)
idx = set.subIndex(*cell,V1[i][0],dim);
else if (side == RIGHT)
idx = set.subIndex(*cell,V2[i][0],dim);
LeafVertexIterator it=start;
for (; it != itend; ++it) {
if (mapper.map(*it) == idx)
break;
}
if (isOnPlane(dir,it->geometry().corner(0),coord)) {
Dune::FieldVector<double, 3> pos = gv.vertexPosition(idx);
if (isOnPlane(dir,pos,coord)) {
for (int j=0;j<4;++j) {
if (side == LEFT)
idx = set.subIndex(*cell,V1[i][j],dim);
if (side == RIGHT)
idx = set.subIndex(*cell,V2[i][j],dim);
LeafVertexIterator it=start;
for (; it != itend; ++it) {
if (mapper.map(*it) == idx)
break;
}
if (!isOnPlane(dir,it->geometry().corner(0),coord))
pos = gv.vertexPosition(idx);
if (!isOnPlane(dir,pos,coord))
continue;
BoundaryGrid::Vertex v;
BoundaryGrid::extract(v,it->geometry().corner(0),i);
BoundaryGrid::extract(v,pos,i);
v.i = idx;
verts.push_back(v);
}
@ -173,9 +159,10 @@ void ElasticityUpscale<GridType>::addMPC(Direction dir, int slave,
}
template<class GridType>
void ElasticityUpscale<GridType>::periodicPlane(Direction plane, Direction dir,
const std::vector<BoundaryGrid::Vertex>& slave,
const BoundaryGrid& master)
void ElasticityUpscale<GridType>::periodicPlane(Direction plane,
Direction dir,
const std::vector<BoundaryGrid::Vertex>& slave,
const BoundaryGrid& master)
{
for (size_t i=0;i<slave.size();++i) {
BoundaryGrid::Vertex coord;
@ -187,152 +174,44 @@ void ElasticityUpscale<GridType>::periodicPlane(Direction plane, Direction dir,
}
}
template<class GridType>
Matrix ElasticityUpscale<GridType>::findBMatrixLLM(const SlaveGrid& slave)
{
std::vector< std::set<int> > adj;
adj.resize(A.getEqns());
int cols=0;
std::vector<int> colmap;
cols=0;
for (std::map<int,BoundaryGrid::Vertex>::const_iterator it = slave.begin();
it != slave.end();++it) {
for (int k=0;k<dim;++k) {
if (A.isFixed(it->first))
continue;
MPC* mpc = A.getMPC(it->first,k);
if (mpc) {
for (size_t n=0;n<mpc->getNoMaster();++n) {
int idx = A.getEquationForDof(mpc->getMaster(n).node,
mpc->getMaster(n).dof-1);
if (idx != -1)
adj[idx].insert(cols++);
}
} else {
int idx = A.getEquationForDof(it->first,k);
if (idx != -1)
adj[idx].insert(cols++);
}
}
}
Matrix B;
MatrixOps::fromAdjacency(B,adj,A.getEqns(),cols);
int col=0;
for (std::map<int,BoundaryGrid::Vertex>::const_iterator it = slave.begin();
it != slave.end();++it) {
for (int k=0;k<dim;++k) {
if (A.isFixed(it->first))
continue;
MPC* mpc = A.getMPC(it->first,k);
if (mpc) {
for (size_t n=0;n<mpc->getNoMaster();++n) {
int idx = A.getEquationForDof(mpc->getMaster(n).node,
mpc->getMaster(n).dof-1);
if (idx != -1)
B[idx][col++] += 1;
}
} else {
int idx = A.getEquationForDof(it->first,k);
if (idx != -1)
B[idx][col++] += 1;
}
}
}
return B;
}
template<class GridType>
Matrix ElasticityUpscale<GridType>::findLMatrixLLM(const SlaveGrid& slave,
const BoundaryGrid& master)
{
int nbeqn=0;
std::vector<int> dofmap(master.totalNodes()*dim,-1);
int col=0;
for (size_t i=0;i<master.totalNodes();++i) {
if (master.isFixed(i)) {
dofmap[i*dim ] = -1;
dofmap[i*dim+1] = -1;
dofmap[i*dim+2] = -1;
}
else {
dofmap[i*dim ] = col++;
dofmap[i*dim+1] = col++;
dofmap[i*dim+2] = col++;
}
}
for (SlaveGrid::const_iterator it = slave.begin();
it != slave.end(); ++it) {
if (!A.isFixed(it->first))
nbeqn += 3;
}
std::vector< std::set<int> > adj;
adj.resize(nbeqn);
int row=0;
for (SlaveGrid::const_iterator it = slave.begin();
it != slave.end();++it) {
if (A.isFixed(it->first))
continue;
for (int k=0;k<dim;++k) {
if (it->second.i == -1) {
for (int i=0;i<4;++i) {
int idx = dofmap[it->second.q->v[i].i*dim+k];
if (idx > -1)
adj[row].insert(idx);
}
}
else {
int idx = dofmap[it->second.i*dim+k];
if (idx > -1)
adj[row].insert(idx);
}
row++;
}
}
Matrix L;
MatrixOps::fromAdjacency(L,adj,nbeqn,col);
row = 0;
for (SlaveGrid::const_iterator it = slave.begin();
it != slave.end();++it) {
if (A.isFixed(it->first))
continue;
for (int k=0;k<dim;++k) {
if (it->second.i == -1) {
std::vector<double> N = it->second.q->evalBasis(it->second.c[0],
it->second.c[1]);
for (int i=0;i<4;++i) {
int idx = dofmap[it->second.q->v[i].i*dim+k];
if (idx > -1)
L[row][idx] -= N[i];
}
}
else {
int idx = dofmap[it->second.i*dim+k];
if (idx > -1)
L[row][idx] -= 1;
}
row++;
}
}
return L;
}
//! \brief Static helper to renumber a Q4 mesh to a P_1/P_N lag mesh.
//! \param[in] b The boundary grid describing the Q4 mesh
//! \param[in] n1 Number of DOFS in the first direction for each element
//! \param[in] n2 Number of DOFS in the first direction for each element
//! \param[out] totalDOFs The total number of free DOFs
static std::vector< std::vector<int> > renumber(const BoundaryGrid& b,
int n1, int n2)
int n1, int n2,
int& totalDOFs)
{
std::vector<std::vector<int> > nodes;
nodes.resize(b.size());
// loop over elements
int ofs = 0;
totalDOFs = 0;
// fix lower left multiplicator.
// will be "transfered" to all corners through periodic conditions
nodes[0].push_back(-1);
for (size_t e=0; e < b.size(); ++e) {
// first direction major ordered nodes within each element
for (int i2=0; i2 < n2; ++i2) {
if (e != 0)
nodes[e].push_back(nodes[e-1][i2*n1+n1-1]);
for (int i1=(e==0?0:1); i1 < n1; ++i1)
nodes[e].push_back(ofs++);
int start = (e==0 && i2 != 0)?0:1;
// slave the buggers
if (i2 == n2-1 && n2 > 2) {
for (int i1=(e==0?0:1); i1 < n1; ++i1) {
nodes[e].push_back(nodes[e][i1]);
}
} else {
for (int i1=start; i1 < n1; ++i1) {
if (e == b.size()-1)
nodes[e].push_back(nodes[0][i2*n1]);
else
nodes[e].push_back(totalDOFs++);
}
}
}
}
@ -340,44 +219,43 @@ static std::vector< std::vector<int> > renumber(const BoundaryGrid& b,
}
template<class GridType>
Matrix ElasticityUpscale<GridType>::findLMatrixMortar(const BoundaryGrid& b1,
const BoundaryGrid& interface,
int dir, int n1, int n2)
int ElasticityUpscale<GridType>::addBBlockMortar(const BoundaryGrid& b1,
const BoundaryGrid& interface,
int dir, int n1, int n2,
int colofs)
{
std::vector< std::set<int> > adj;
adj.resize(A.getEqns());
#define QINDEX(p,q,dir) (q)
// get a set of P1 shape functions for the displacements
P1ShapeFunctionSet<ctype,ctype,2> ubasis =
P1ShapeFunctionSet<ctype,ctype,2>::instance();
// get a set of PN shape functions for the multipliers
PNShapeFunctionSet<2> lbasis(n1+1, n2+1);
std::vector<std::vector<int> > lnodes = renumber(interface,n1,n2);
int totalNodes = lnodes.back().back()+1;
// renumber the linear grid to the real multiplier grid
int totalEqns;
std::vector<std::vector<int> > lnodes = renumber(interface, n1+1,
n2+1, totalEqns);
if (Bpatt.empty())
Bpatt.resize(A.getEqns());
// process pillar by pillar
for (size_t p=0;p<interface.size();++p) {
for (size_t q=0;q<b1.colSize(p);++q) {
for (size_t i=0;i<ubasis.n;++i) {
for (size_t i=0;i<4;++i) {
for (size_t d=0;d<3;++d) {
MPC* mpc = A.getMPC(b1.getQuad(p,QINDEX(p,q,dir)).v[i].i,d);
MPC* mpc = A.getMPC(b1.getQuad(p,q).v[i].i,d);
if (mpc) {
for (size_t n=0;n<mpc->getNoMaster();++n) {
int dof = A.getEquationForDof(mpc->getMaster(n).node,d);
if (dof > -1) {
for (size_t j=0;j<lnodes[p].size();++j)
adj[dof].insert(3*lnodes[p][j]+d);
for (size_t j=0;j<lnodes[p].size();++j) {
int indexj = 3*lnodes[p][j]+d;
if (indexj > -1)
Bpatt[dof].insert(indexj+colofs);
}
}
}
} else {
int dof = A.getEquationForDof(b1.getQuad(p,QINDEX(p,q,dir)).v[i].i,d);
int dof = A.getEquationForDof(b1.getQuad(p,q).v[i].i,d);
if (dof > -1) {
for (size_t j=0;j<lnodes[p].size();++j)
adj[dof].insert(3*lnodes[p][j]+d);
for (size_t j=0;j<lnodes[p].size();++j) {
int indexj = 3*lnodes[p][j]+d;
if (indexj > -1)
Bpatt[dof].insert(indexj+colofs);
}
}
}
}
@ -385,9 +263,27 @@ Matrix ElasticityUpscale<GridType>::findLMatrixMortar(const BoundaryGrid& b1,
}
}
Matrix B;
MatrixOps::fromAdjacency(B,adj,A.getEqns(),3*totalNodes);
return 3*totalEqns;
}
template<class GridType>
void ElasticityUpscale<GridType>::assembleBBlockMortar(const BoundaryGrid& b1,
const BoundaryGrid& interface,
int dir, int n1,
int n2, int colofs,
double alpha)
{
// get a set of P1 shape functions for the displacements
P1ShapeFunctionSet<ctype,ctype,2> ubasis =
P1ShapeFunctionSet<ctype,ctype,2>::instance();
// get a set of PN shape functions for the multipliers
PNShapeFunctionSet<2> lbasis(n1+1, n2+1);
// renumber the linear grid to the real multiplier grid
int totalEqns;
std::vector<std::vector<int> > lnodes = renumber(interface, n1+1,
n2+1, totalEqns);
// get a reference element
Dune::GeometryType gt;
gt.makeCube(2);
@ -400,6 +296,7 @@ Matrix ElasticityUpscale<GridType>::findLMatrixMortar(const BoundaryGrid& b1,
// do the assembly loop
typename Dune::QuadratureRule<ctype,2>::const_iterator r;
Dune::DynamicMatrix<ctype> E(ubasis.n,(n1+1)*(n2+1),0.0);
LoggerHelper help(interface.size(), 5, 1000);
for (size_t p=0;p<interface.size();++p) {
const BoundaryGrid::Quad& qi(interface[p]);
HexGeometry<2,2,GridType> lg(qi);
@ -409,8 +306,8 @@ Matrix ElasticityUpscale<GridType>::findLMatrixMortar(const BoundaryGrid& b1,
E = 0;
for (r = rule.begin(); r != rule.end();++r) {
ctype detJ = hex.integrationElement(r->position());
if (detJ < 1.e-4 && r == rule.begin())
std::cerr << "warning: interface cell (close to) degenerated, |J|=" << detJ << std::endl;
if (detJ < 0)
assert(0);
typename HexGeometry<2,2,GridType>::LocalCoordinate loc =
lg.local(hex.global(r->position()));
@ -422,6 +319,7 @@ Matrix ElasticityUpscale<GridType>::findLMatrixMortar(const BoundaryGrid& b1,
}
}
}
// and assemble element contributions
for (int d=0;d<3;++d) {
for (int i=0;i<4;++i) {
@ -432,7 +330,8 @@ Matrix ElasticityUpscale<GridType>::findLMatrixMortar(const BoundaryGrid& b1,
if (indexi > -1) {
for (size_t j=0;j<lnodes[p].size();++j) {
int indexj = lnodes[p][j]*3+d;
B[indexi][indexj] += E[i][j];
if (indexj > -1)
B[indexi][indexj+colofs] += alpha*E[i][j];
}
}
}
@ -441,16 +340,16 @@ Matrix ElasticityUpscale<GridType>::findLMatrixMortar(const BoundaryGrid& b1,
if (indexi > -1) {
for (size_t j=0;j<lnodes[p].size();++j) {
int indexj = lnodes[p][j]*3+d;
B[indexi][indexj] += E[i][j];
if (indexj > -1)
B[indexi][indexj+colofs] += alpha*E[i][j];
}
}
}
}
}
}
help.log(p, "\t\t\t... still processing ... pillar ");
}
return B;
}
template<class GridType>
@ -557,9 +456,11 @@ void ElasticityUpscale<GridType>::assemble(int loadcase, bool matrix)
int m=0;
Dune::FieldMatrix<ctype,dim*bfunc,dim*bfunc>* KP=0;
if (matrix) {
A.getOperator() = 0;
KP = &K;
A.getOperator() = 0;
}
LoggerHelper help(gv.size(0), 5, 50000);
for (LeafIterator it = gv.leafView().template begin<0>(); it != itend; ++it) {
materials[m++]->getConstitutiveMatrix(C);
// determine geometry type of the current element and get the matching reference element
@ -578,7 +479,7 @@ void ElasticityUpscale<GridType>::assemble(int loadcase, bool matrix)
it->geometry().jacobianInverseTransposed(r->position());
ctype detJ = it->geometry().integrationElement(r->position());
if (detJ <= 1.e-4) {
if (detJ <= 1.e-5 && verbose) {
std::cout << "cell " << m << " is (close to) degenerated, detJ " << detJ << std::endl;
double zdiff=0.0;
for (int i=0;i<4;++i)
@ -604,6 +505,7 @@ void ElasticityUpscale<GridType>::assemble(int loadcase, bool matrix)
}
A.addElement(KP,EP,it,(loadcase > -1)?&b[loadcase]:NULL);
help.log(m, "\t\t... still processing ... cell ");
}
}
@ -643,8 +545,6 @@ void ElasticityUpscale<GridType>::averageStress(Dune::FieldVector<ctype,comp>& s
it->geometry().jacobianInverseTransposed(r->position());
ctype detJ = it->geometry().integrationElement(r->position());
if (detJ <= 0) // wtf ?
continue;
volume += detJ*r->weight();
@ -720,12 +620,14 @@ void ElasticityUpscale<GridType>::loadMaterialsFromGrid(const std::string& file)
MaterialMap::iterator it;
if ((it = cache.find(std::make_pair(Emod[k],Poiss[k]))) != cache.end())
{
assert(gv.cellVolume(i) > 0);
volume[it->second] += gv.cellVolume(i);
materials.push_back(it->second);
}
else {
Material* mat = new Isotropic(j++,Emod[k],Poiss[k]);
cache.insert(std::make_pair(std::make_pair(Emod[k],Poiss[k]),mat));
assert(gv.cellVolume(i) > 0);
volume.insert(std::make_pair(mat,gv.cellVolume(i)));
materials.push_back(mat);
}
@ -782,7 +684,7 @@ void ElasticityUpscale<GridType>::loadMaterialsFromRocklist(const std::string& f
std::vector<int> cells = gv.globalCell();
for (size_t i=0;i<cells.size();++i) {
int k = cells[i];
if (satnum[k]-1 >= cache.size()) {
if (satnum[k]-1 >= (int)cache.size()) {
std::cerr << "Material " << satnum[k] << " referenced but not available. Check your rocklist." << std::endl;
exit(1);
}
@ -857,254 +759,215 @@ void ElasticityUpscale<GridType>::periodicBCsMortar(const double* min,
// 2. establishes strong couplings (MPC) on top and bottom
// 3. extracts and establishes a quad grid for the left/right/front/back sides
// 4. establishes grids for the dual dofs
// 5. calculates the coupling matrix L1 between the left/right sides
// 6. calculates the coupling matrix L2 between the front/back sides
// 5. calculates the coupling matrix between the left/right sides
// 6. calculates the coupling matrix between the front/back sides
//
// The Mortar linear system is of the form
// [A B] [u] [b]
// [B' 0] [l] = [0]
// step 1
fixCorners(min,max);
// step 2
std::cout << "\textracting nodes on top face..." << std::endl;
slave.push_back(extractFace(Z,max[2]));
std::cout << "\treconstructing bottom face..." << std::endl;
BoundaryGrid bottom = extractMasterFace(Z,min[2]);
std::cout << "\testablishing couplings on top/bottom..." << std::endl;
periodicPlane(Z,XYZ,slave[0],bottom);
std::cout << "\tinitializing matrix..." << std::endl;
A.initForAssembly();
// step 3
master.push_back(extractMasterFace(X,min[0],LEFT,true));
master.push_back(extractMasterFace(X,max[0],RIGHT,true));
master.push_back(extractMasterFace(Y,min[1],LEFT,true));
master.push_back(extractMasterFace(Y,max[1],RIGHT,true));
std::cout << "\treconstructing left face..." << std::endl;
master.push_back(extractMasterFace(X, min[0], LEFT, true));
std::cout << "\treconstructing right face..." << std::endl;
master.push_back(extractMasterFace(X, max[0], RIGHT, true));
std::cout << "\treconstructing front face..." << std::endl;
master.push_back(extractMasterFace(Y, min[1], LEFT, true));
std::cout << "\treconstructing back face..." << std::endl;
master.push_back(extractMasterFace(Y, max[1], RIGHT, true));
std::cout << "Xsides " << master[0].size() << " " << master[1].size() << std::endl
<< "Ysides " << master[2].size() << " " << master[3].size() << std::endl
<< "Zmaster " << bottom.size() << std::endl;
std::cout << "Establish YZ multiplier grid with " << n2 << "x1" << " elements" << std::endl;
std::cout << "\testablished YZ multiplier grid with " << n2 << "x1" << " elements" << std::endl;
// step 4
BoundaryGrid::FaceCoord fmin,fmax;
fmin[0] = min[1]; fmin[1] = min[2];
fmax[0] = max[1]; fmax[1] = max[2];
BoundaryGrid lambdax = BoundaryGrid::uniform(fmin,fmax,n2,1,true);
BoundaryGrid lambdax = BoundaryGrid::uniform(fmin, fmax, n2, 1, true);
fmin[0] = min[0]; fmin[1] = min[2];
fmax[0] = max[0]; fmax[1] = max[2];
std::cout << "Establish XZ multiplier grid with " << n1 << "x1" << " elements" << std::endl;
BoundaryGrid lambday = BoundaryGrid::uniform(fmin,fmax,n1,1,true);
std::cout << "\testablished XZ multiplier grid with " << n1 << "x1" << " elements" << std::endl;
BoundaryGrid lambday = BoundaryGrid::uniform(fmin, fmax, n1, 1, true);
addBBlockMortar(master[0], lambdax, 0, 1, p2, 0);
int eqns = addBBlockMortar(master[1], lambdax, 0, 1, p2, 0);
addBBlockMortar(master[2], lambday, 1, 1, p1, eqns);
int eqns2 = addBBlockMortar(master[3], lambday, 1, 1, p1, eqns);
MatrixOps::fromAdjacency(B, Bpatt, A.getEqns(), eqns+eqns2);
Bpatt.clear();
// step 5
Matrix L1 = findLMatrixMortar(master[0],lambdax,0, p2, 1);
Matrix L2 = findLMatrixMortar(master[1],lambdax,0, p2, 1);
L.push_back(MatrixOps::Axpy(L1,L2,-1));
std::cout << "\tassembling YZ mortar matrix..." << std::endl;
assembleBBlockMortar(master[0], lambdax, 0, 1, p2, 0);
assembleBBlockMortar(master[1], lambdax, 0, 1, p2, 0, -1.0);
// step 6
Matrix L3 = findLMatrixMortar(master[2],lambday,1, p1, 1);
Matrix L4 = findLMatrixMortar(master[3],lambday,1, p1, 1);
L.push_back(MatrixOps::Axpy(L3,L4,-1));
std::cout << "\tassembling XZ mortar matrix..." << std::endl;
assembleBBlockMortar(master[2], lambday, 1, 1, p1, eqns);
assembleBBlockMortar(master[3], lambday, 1, 1, p1, eqns, -1.0);
master.clear();
slave.clear();
}
template<class GridType>
void ElasticityUpscale<GridType>::periodicBCsLLM(const double* min,
const double* max,
int n1, int n2)
void ElasticityUpscale<GridType>::setupAMG(int pre, int post,
int target, int zcells)
{
// this method
// 1. fixes the primal corner dofs
// 2. fixes the primal dofs on the skeleton
// 3. establishes strong couplings (MPC) on top and bottom
// 4. extracts a point grid for the left/right/front/back sides,
// and establishes a uniform interface grid in each direction
// 5. calculates the coupling matrices B1-4 and L1-4
Criterion crit;
ElasticityAMG::SmootherArgs args;
args.relaxationFactor = 1.0;
crit.setCoarsenTarget(target);
crit.setGamma(1);
crit.setNoPreSmoothSteps(pre);
crit.setNoPostSmoothSteps(post);
crit.setDefaultValuesIsotropic(3, zcells);
// step 1
fixCorners(min,max);
std::cout << "\t collapsing 2x2x" << zcells << " cells per level" << std::endl;
op = new Operator(A.getOperator());
upre = new ElasticityAMG(*op, crit, args);
// step 2
fixLine(X,min[1],min[2]);
fixLine(X,max[1],min[2]);
fixLine(X,min[1],max[2]);
fixLine(X,max[1],max[2]);
fixLine(Y,min[0],min[2]);
fixLine(Y,max[0],min[2]);
fixLine(Y,min[0],max[2]);
fixLine(Y,max[0],max[2]);
fixLine(Z,min[0],min[1]);
fixLine(Z,max[0],min[1]);
fixLine(Z,min[0],max[1]);
fixLine(Z,max[0],max[1]);
// step 3
std::vector<BoundaryGrid::Vertex> Zs = extractFace(Z,max[2]);
BoundaryGrid Zm = extractMasterFace(Z,min[2]);
periodicPlane(Z,XYZ,Zs,Zm);
A.initForAssembly();
// step 4
slave.push_back(extractFace(X,min[0]));
slave.push_back(extractFace(X,max[0]));
slave.push_back(extractFace(Y,min[1]));
slave.push_back(extractFace(Y,max[1]));
Dune::LeafMultipleCodimMultipleGeomTypeMapper<GridType,
Dune::MCMGVertexLayout> mapper(gv);
BoundaryGrid::FaceCoord fmin,fmax;
// YZ plane
fmin[0] = min[1]; fmin[1] = min[2];
fmax[0] = max[1]; fmax[1] = max[2];
std::cout << "Establish YZ interface grid with " << n2 << "x1" << " elements" << std::endl;
master.push_back(BoundaryGrid::uniform(fmin,fmax,n2,1));
// XZ plane
fmin[0] = min[0]; fmin[1] = min[2];
fmax[0] = max[0]; fmax[1] = max[2];
std::cout << "Establish XZ interface grid with " << n1 << "x1" << " elements" << std::endl;
master.push_back(BoundaryGrid::uniform(fmin,fmax,n1,1));
// step 5
std::map<int,BoundaryGrid::Vertex> m;
// find matching coefficients
for (size_t i=0;i<master.size();++i) {
for (int j=0;j<2;++j) {
m.clear();
for (size_t k=0;k<slave[i*2+j].size();++k) {
BoundaryGrid::Vertex c;
if (master[i].find(c,slave[i*2+j][k])) {
m.insert(std::make_pair(slave[i*2+j][k].i,c));
} else
assert(0);
}
B.push_back(findBMatrixLLM(m));
L.push_back(findLMatrixLLM(m,master[i]));
}
}
/*
amg->addContext("Apre");
amg->setContext("Apre");
Vector x,y;
// this is done here to make sure we are in a single-threaded section
// will have to be redone when AMG is refactored upstream
amg->pre(x,y);
*/
}
template<class GridType>
void ElasticityUpscale<GridType>::setupPreconditioner()
template<class GridType>
void ElasticityUpscale<GridType>::setupSolvers(const LinSolParams& params)
{
// setup subdomain maps
typename OverlappingSchwarz::subdomain_vector rows;
const LeafIterator itend = gv.leafView().template end<0>();
const LeafIndexSet& set = gv.leafView().indexSet();
rows.resize(gv.logicalCartesianSize()[0]*gv.logicalCartesianSize()[1]);
int cell=0, currdomain=0;
for (LeafIterator it = gv.leafView().template begin<0>(); it != itend; ++it, ++cell) {
if (cell / gv.logicalCartesianSize()[2] > 0
&& cell % gv.logicalCartesianSize()[2] == 0)
currdomain++;
for (size_t i=0;i<8;++i) {
int idx=set.subIndex(*it,i,dim);
for (int d=0;d<3;++d) {
MPC* mpc = A.getMPC(idx,d);
if (mpc) {
for (size_t n=0;n<mpc->getNoMaster();++n) {
int row = A.getEquationForDof(mpc->getMaster(n).node,d);
if (row > -1)
rows[currdomain].insert(row);
int siz = A.getOperator().N(); // system size
if (params.type == ITERATIVE) {
setupAMG(params.steps[0], params.steps[1], params.coarsen_target,
params.zcells);
// Mortar in use
if (B.N()) {
siz += B.M();
// schur system: B'*diag(A)^-1*B
if (params.mortarpre == SCHURAMG) {
Vector v, v2, v3;
v.resize(B.N());
v2.resize(B.N());
v = 0;
v2 = 0;
Dune::DynamicMatrix<double> T(B.M(), B.M());
upre->pre(v, v);
std::cout << "\tBuilding preconditioner for multipliers..." << std::endl;
MortarBlockEvaluator<Dune::Preconditioner<Vector,Vector> > pre(*upre, B);
LoggerHelper help(B.M(), 10, 100);
for (size_t i=0; i < B.M(); ++i) {
v[i] = 1;
pre.apply(v, v2);
for (size_t j=0; j < B.M(); ++j)
T[j][i] = v2[j];
v[i] = 0;
help.log(i, "\t\t... still processing ... multiplier ");
}
upre->post(v);
P = MatrixOps::fromDense(T);
} else if (params.mortarpre == SCHURDIAG) {
Matrix D = MatrixOps::diagonal(A.getEqns());
// scale by row sums
size_t row=0;
for (Matrix::ConstRowIterator it = A.getOperator().begin();
it != A.getOperator().end(); ++it, ++row) {
double alpha=0;
for (Matrix::ConstColIterator it2 = it->begin();
it2 != it->end(); ++it2) {
if (it2.index() != row)
alpha += fabs(*it2);
}
D[row][row] = 1.0/(A.getOperator()[row][row]/alpha);
}
Matrix t1;
// t1 = Ad*B
Dune::matMultMat(t1, D, B);
// P = B'*t1 = B'*Ad*B
Dune::transposeMatMultMat(P, B, t1);
}
if (params.uzawa) {
Dune::CGSolver<Vector>* innersolver =
new Dune::CGSolver<Vector>(*op, *upre, params.tol,
params.maxit, verbose?2:0);
op2 = new SchurEvaluator(*innersolver, B);
lpre = new SeqLU<Matrix, Vector, Vector>(P);
Dune::CGSolver<Vector>* outersolver =
new Dune::CGSolver<Vector>(*op2, *lpre, params.tol*10,
params.maxit, verbose?2:0);
solver = new UzawaSolver<Vector, Vector>(innersolver, outersolver, B);
} else {
mpre = new MortarSchurPre<ElasticityAMG>(P, B, *upre, params.symmetric);
meval = new MortarEvaluator(A.getOperator(), B);
if (params.symmetric) {
solver = new Dune::MINRESSolver<Vector>(*meval, *mpre,
params.tol,
params.maxit,
verbose?2:0);
} else {
int row = A.getEquationForDof(idx,d);
if (row > -1)
rows[currdomain].insert(row);
solver = new Dune::RestartedGMResSolver<Vector>(*meval, *mpre,
params.tol,
params.restart,
params.maxit,
verbose?2:0, true);
}
}
} else {
solver = new Dune::CGSolver<Vector>(*op, *upre, params.tol,
params.maxit, verbose?2:0);
}
}
ovl = new OverlappingSchwarz(A.getOperator(),rows,1,false);
}
template<class GridType>
void ElasticityUpscale<GridType>::setupSolvers(Solver solver)
{
if (!B.empty() && A.getOperator().N() == A.getEqns()) { // LLM in use
// The LLM linear system is of the form
// [A B1 B2 B3 B4 0 0] [ u] [b]
// [B1' 0 0 0 0 L1 0] [ l_1] [0]
// [B2' 0 0 0 0 L2 0] [ l_2] [0]
// [B3' 0 0 0 0 0 L3] [ l_3] = [0]
// [B4' 0 0 0 0 0 L4] [ l_4] [0]
// [ 0 L1' L2' 0 0 0 0] [ub_1] [0]
// [ 0 0 0 L3' L4' 0 0] [ub_2] [0]
int r = A.getOperator().N();
int c = A.getOperator().M();
for (size_t i=0;i<B.size();++i) {
A.getOperator() = MatrixOps::augment(A.getOperator(),B[i],0,c,true);
c += B[i].M();
}
for (size_t i=0;i<L.size();++i) {
A.getOperator() = MatrixOps::augment(A.getOperator(),L[i],r,c,true);
r += L[i].N();
if (i % 2 == 1)
c += L[i].M();
}
int siz=A.getOperator().N();
for (int i=0;i<6;++i)
b[i].resize(A.getOperator().N());
}
// Mortar in use
if (B.empty() && !L.empty() && A.getOperator().N() == A.getEqns()) {
// The Mortar linear system is of the form
// [A L1 L2] [ u] [b]
// [L1' 0 0] [l_1] = [0]
// [L2' 0 0] [l_2] [0]
int c = A.getOperator().M();
for (size_t i=0;i<L.size();++i) {
A.getOperator() = MatrixOps::augment(A.getOperator(),L[i],0,c,true);
c += L[i].M();
}
int siz=A.getLoadVector().size();
for (int i=0;i<6;++i) {
b[i].resize(A.getOperator().N());
for (int j=siz;j<b[i].size();++j)
b[i][j] = 0;
}
if (L.empty() && B.empty()) { // MPC
// overlapping schwarz is much more memory efficient
// and usually more effective
setupPreconditioner();
}
}
if (solver == SLU) {
} else {
if (B.N())
A.getOperator() = MatrixOps::augment(A.getOperator(), B,
0, A.getOperator().M(), true);
#if HAVE_SUPERLU
if (!slu)
slu = new Dune::SuperLU<Matrix>(A.getOperator(),false);
solver = new Dune::SuperLU<Matrix>(A.getOperator(),verbose);
#else
std::cerr << "SuperLU solver not enabled" << std::endl;
exit(1);
#endif
} else if (solver == CG) {
if (!op)
op = new Dune::MatrixAdapter<Matrix,Vector,Vector>(A.getOperator());
if (!ovl && !ilu)
ilu = new Dune::SeqILU0<Matrix,Vector,Vector>(A.getOperator(),0.92);
if (!cgsolver) {
if (ovl)
cgsolver = new Dune::CGSolver<Vector>(*op, *ovl, tol, 5000, 1);
else
cgsolver = new Dune::CGSolver<Vector>(*op, *ilu, tol, 5000, 1);
}
siz = A.getOperator().N();
}
for (int i=0;i<6;++i)
b[i].resize(siz);
}
template<class GridType>
void ElasticityUpscale<GridType>::solve(Solver solver, double tol, int loadcase)
void ElasticityUpscale<GridType>::solve(int loadcase)
{
Dune::InverseOperatorResult r;
try {
Dune::InverseOperatorResult r;
u[loadcase].resize(b[loadcase].size(), false);
u[loadcase] = 0;
// initialize u to some arbitrary value
u[loadcase].resize(A.getOperator().N(), false);
u[loadcase] = 2.0;
solver->apply(u[loadcase], b[loadcase], r);
#if HAVE_SUPERLU
if (solver == SLU) {
slu->apply(u[loadcase], b[loadcase], r);
std::cout << "\tsolution norm: " << u[loadcase].two_norm() << std::endl;
} catch (Dune::ISTLError& e) {
std::cerr << "exception thrown " << e << std::endl;
}
else
#endif
if (solver == CG) {
cgsolver->apply(u[loadcase], b[loadcase], r);
}
std::cout << "\t solution norm: " << u[loadcase].two_norm() << std::endl;
}

31
dune/elasticity/logutils.hpp Normal file
View File

@ -0,0 +1,31 @@
//==============================================================================
//!
//! \file logutils.hpp
//!
//! \date Nov 9 2011
//!
//! \author Arne Morten Kvarving / SINTEF
//!
//! \brief Logging helper utilities
//!
//==============================================================================
#pragma once
class LoggerHelper {
public:
LoggerHelper(int max_, int chunks, int minsize)
: per_chunk(max_/chunks), max(max_)
{
if (max < minsize)
per_chunk=-1;
}
void log(int it, const std::string& prefix)
{
if (per_chunk > 0 && it > 0 && (it % per_chunk) == 0 && max-it > per_chunk)
std::cout << prefix << it << '/' << max << std::endl;
}
protected:
int per_chunk; //!< Will log for each per_chunk processed
int max; //!< Total number of its
};

15
dune/elasticity/matrixops.hpp
View File

@ -36,6 +36,11 @@ class MatrixOps {
static void fromAdjacency(Matrix& A, const AdjacencyPattern& adj,
int rows, int cols);
//! \brief Create a sparse matrix from a dense matrix
//! \param[in] T The dense matrix
//! \returns The sparse matrix
static Matrix fromDense(const Dune::DynamicMatrix<double>& T);
//! \brief Print a matrix to stdout
//! \param[in] A The matrix to print
static void print(const Matrix& A);
@ -61,6 +66,16 @@ class MatrixOps {
//! \returns M = diag(A)
static Matrix extractDiagonal(const Matrix& A);
//! \brief Returns a diagonal matrix
//! \param[in] N The dimension of the matrix
static Matrix diagonal(size_t N);
//! \brief Extract a subblock of a matrix into a new matrix
//! \param[in] The matrix to extract from
//! \returns The subblock
static Matrix extractBlock(const Matrix& A,
size_t r0, size_t N, size_t c0, size_t M);
//! \brief Save a matrix as a dense asc file
//! \param[in] A The matrix to save
//! \param[in] file File name

75
dune/elasticity/matrixops_impl.hpp
View File

@ -12,7 +12,7 @@
void MatrixOps::fromAdjacency(Matrix& A, const std::vector< std::set<int> >& adj,
int rows, int cols)
{
{
size_t sum=0;
for (size_t i=0;i<adj.size();++i)
sum += adj[i].size();
@ -34,6 +34,30 @@ void MatrixOps::fromAdjacency(Matrix& A, const std::vector< std::set<int> >& adj
A = 0;
}
Matrix MatrixOps::fromDense(const Dune::DynamicMatrix<double>& T)
{
AdjacencyPattern a;
a.resize(T.N());
for (size_t i=0; i < T.N(); ++i) {
for (size_t j=0; j < T.M(); ++j) {
if (fabs(T[i][j]) > 1.e-14)
a[i].insert(j);
}
}
Matrix result;
fromAdjacency(result, a, T.N(), T.M());
for (Matrix::ConstIterator it = result.begin();
it != result.end(); ++it) {
for (Matrix::ConstColIterator it2 = it->begin();
it2 != it->end();++it2)
result[it.index()][it2.index()] = T[it.index()][it2.index()];
}
return result;
}
void MatrixOps::print(const Matrix& A)
{
for (Matrix::ConstIterator it = A.begin();
@ -159,6 +183,18 @@ Matrix MatrixOps::extractDiagonal(const Matrix& A)
return result;
}
Matrix MatrixOps::diagonal(size_t N)
{
AdjacencyPattern adj;
adj.resize(N);
for (size_t i=0;i<N;++i)
adj[i].insert(i);
Matrix result;
fromAdjacency(result,adj,N,N);
return result;
}
void MatrixOps::saveAsc(const Matrix& A, const std::string& file)
{
std::ofstream f;
@ -188,3 +224,40 @@ void MatrixOps::saveAsc(const Matrix& A, const std::string& file)
}
f.close();
}
Matrix MatrixOps::extractBlock(const Matrix& A, size_t r0, size_t N,
size_t c0, size_t M)
{
// establish adjacency pattern
AdjacencyPattern adj;
adj.resize(N);
for (Matrix::ConstIterator it = A.begin();
it != A.end(); ++it) {
if (it.index() < r0 || it.index() >= N+r0)
continue;
for (Matrix::ConstColIterator it2 = it->begin();
it2 != it->end(); ++it2) {
if (it2.index() < c0 || it2.index() >= M+c0)
continue;
adj[it.index()-r0].insert(it2.index()-c0);
}
}
Matrix result;
fromAdjacency(result, adj, N, M);
// now insert elements from A
for (Matrix::ConstIterator it = A.begin();
it != A.end(); ++it) {
if (it.index() < r0 || it.index() >= N+r0)
continue;
for (Matrix::ConstColIterator it2 = it->begin();
it2 != it->end(); ++it2) {
if (it2.index() < c0 || it2.index() >= M+c0)
continue;
result[it.index()-r0][it2.index()-c0] = *it2;
}
}
return result;
}

79
dune/elasticity/mortar_evaluator.hpp Normal file
View File

@ -0,0 +1,79 @@
//==============================================================================
//!
//! \file mortar_evaluator.hpp
//!
//! \date Nov 15 2012
//!
//! \author Arne Morten Kvarving / SINTEF
//!
//! \brief Linear operator for a Mortar block
//!
//==============================================================================
#pragma once
#include <dune/istl/matrixmatrix.hh>
#include <dune/elasticity/matrixops.hpp>
#include <dune/elasticity/mortar_utils.hpp>
namespace Opm {
namespace Elasticity {
/*! This implements a operator evaluation for the
* schur mortar-block S = B^T*A^-1*B
!*/
class MortarEvaluator : public Dune::LinearOperator<Vector, Vector> {
public:
// define the category
enum {
//! \brief The category the preconditioner is part of.
category=Dune::SolverCategory::sequential
};
//! \brief Constructor
//! \param[in] Ai Evaluator for A^-1
//! \param[in] B The mortar coupling matrix
MortarEvaluator(const Matrix& A_,
const Matrix& B_) :
A(A_), B(B_)
{
}
//! \brief Apply the multiplier block
//! \param[in] x The vector to apply the operator to
//! \param[out] y The result of the operator evaluation
void apply(const Vector& x, Vector& y) const
{
Vector lambda, l2;
MortarUtils::extractBlock(lambda, x, B.M(), A.N());
A.mv(x, y);
B.umv(lambda, y);
l2.resize(lambda.size());
B.mtv(x, l2);
MortarUtils::injectBlock(y, l2, B.M(), A.N());
}
//! \brief Apply the multiplier block with an embedded axpy
//! \param[in] alpha The scalar to scale with
//! \param[in] x The vector to apply the operator to
//! \param[out] y The result of the operator evaluation
void applyscaleadd(field_type alpha, const Vector& x, Vector& y) const
{
Vector lambda, l2;
A.usmv(alpha, x, y);
MortarUtils::extractBlock(lambda, x, B.M(), A.N());
B.umv(lambda, y);
l2.resize(lambda.size());
B.mtv(x, l2);
for (size_t i=A.N(); i < y.size(); ++i)
y[i] += alpha*l2[i-A.N()];
}
protected:
const Matrix& A; //!< Reference to A matrix
const Matrix& B; //!< Reference to the mortar coupling matrix
};
}
}

73
dune/elasticity/mortar_schur.hpp Normal file
View File

@ -0,0 +1,73 @@
//==============================================================================
//!
//! \file mortar_schur.hpp
//!
//! \date Nov 15 2012
//!
//! \author Arne Morten Kvarving / SINTEF
//!
//! \brief Schur based linear operator for a Mortar block
//!
//==============================================================================
#pragma once
#include <dune/istl/matrixmatrix.hh>
#include <dune/elasticity/applier.hpp>
#include <dune/elasticity/matrixops.hpp>
namespace Opm {
namespace Elasticity {
/*! This implements a operator evaluation for the
* schur mortar-block S = B^T*A^-1*B
!*/
template<class T>
class MortarBlockEvaluator : public Dune::LinearOperator<Vector, Vector> {
public:
// define the category
enum {
//! \brief The category the preconditioner is part of.
category=Dune::SolverCategory::sequential
};
//! \brief Constructor
//! \param[in] Ai Solver or preconditioner for A^-1
//! \param[in] B The mortar coupling matrix
MortarBlockEvaluator(T& Ai_,
const Matrix& B_) :
Ai(Ai_), B(B_), op(Ai)
{
}
//! \brief Apply the multiplier block
//! \param[in] x The vector to apply the operator to
//! \param[out] y The result of the operator evaluation
void apply(const Vector& x, Vector& y) const
{
y = 0;
applyscaleadd(1.0, x, y);
}
//! \brief Apply the multiplier block with an embedded axpy
//! \param[in] alpha The scalar to scale with
//! \param[in] x The vector to apply the operator to
//! \param[out] y The result of the operator evaluation
void applyscaleadd(field_type alpha, const Vector& x, Vector& y) const
{
Vector temp1(B.N());
B.mv(x, temp1);
Vector temp(temp1.size());
Dune::InverseOperatorResult res;
temp = 0;
op.apply(temp, temp1);
B.usmtv(alpha, temp, y);
}
protected:
T& Ai; //!< Reference to solver or evaluator for inverse operator
const Matrix& B; //!< Reference to the mortar coupling matrix
mutable OperatorApplier<T> op; //!< Applier for the preconditioner / inverse solver
};
}
}

123
dune/elasticity/mortar_schur_precond.hpp Normal file
View File

@ -0,0 +1,123 @@
//==============================================================================
//!
//! \file schur_precond.hpp
//!
//! \date Nov 15 2012
//!
//! \author Arne Morten Kvarving / SINTEF
//!
//! \brief Schur based preconditioner for saddle point systems
//!
//==============================================================================
#pragma once
#include <dune/istl/preconditioners.hh>
#include <dune/istl/matrixmatrix.hh>
#include <dune/elasticity/matrixops.hpp>
#include <dune/elasticity/mortar_utils.hpp>
namespace Opm {
namespace Elasticity {
/*! This implements a Schur-decomposition based preconditioner for the
* mortar-elasticity system
* [A B]
* [B' ]
*
* The preconditioner is
* [Apre B]
* [ P]
* Here Apre is some preconditioner for A and P some preconditioner for
* S = B^TA^-1B
!*/
template<class PrecondElasticityBlock>
class MortarSchurPre : public Dune::Preconditioner<Vector,Vector> {
public:
// define the category
enum {
//! \brief The category the preconditioner is part of.
category=Dune::SolverCategory::sequential
};
//! \brief Constructor
//! \param[in] P The multiplier block with diagonal A approximation
//! \param[in] B The mortar coupling matrix
//! \param[in] Apre_ A preconfigured preconditioner for A
//! \param[in] symmetric If true, use symmetric preconditioning
MortarSchurPre(const Matrix& P_, const Matrix& B_,
PrecondElasticityBlock& Apre_, bool symmetric_=false) :
Apre(Apre_), B(B_), N(B.N()), M(B.M()),
Lpre(P_, false, false), symmetric(symmetric_)
{
}
//! \brief Destructor
virtual ~MortarSchurPre()
{
}
//! \brief Preprocess preconditioner
virtual void pre(Vector& x, Vector& b)
{
// extract displacement DOFs
Vector tempx, tempb;
MortarUtils::extractBlock(tempx, x, N);
MortarUtils::extractBlock(tempb, b, N);
Apre.pre(tempx, tempb);
MortarUtils::injectBlock(x, tempx, N);
MortarUtils::injectBlock(b, tempb, N);
}
//! \brief Applies the preconditioner
//! \param[out] v The resulting vector
//! \param[in] d The vector to apply the preconditioner to
virtual void apply(Vector& v, const Vector& d)
{
// multiplier block (second row)
Vector temp;
MortarUtils::extractBlock(temp, d, M, N);
Dune::InverseOperatorResult r;
Vector temp2;
temp2.resize(temp.size(), false);
Lpre.apply(temp2, temp, r);
MortarUtils::injectBlock(v, temp2, M, N);
// elasticity block (first row)
MortarUtils::extractBlock(temp, d, N);
if (!symmetric)
B.mmv(temp2, temp);
temp2.resize(N, false);
temp2 = 0;
Apre.apply(temp2, temp);
MortarUtils::injectBlock(v, temp2, N);
}
//! \brief Dummy post-process function
virtual void post(Vector& x)
{
Apre.post(x);
}
protected:
//! \brief The preconditioner for the elasticity operator
PrecondElasticityBlock& Apre;
//! \brief The mortar coupling matrix
const Matrix& B;
//! \brief Number of displacement DOFs
int N;
//! \brief Number of multiplier DOFs
int M;
//! \brief Linear solver for the multiplier block
Dune::SuperLU<Matrix> Lpre;
//! \brief Whether or not to use a symmetric preconditioner
bool symmetric;
};
}
}

31
dune/elasticity/mortar_utils.hpp Normal file
View File

@ -0,0 +1,31 @@
#pragma once
namespace Opm {
namespace Elasticity {
class MortarUtils {
public:
//! \brief Extract a range of indices from a vector
//! \param[out] x The vector with the extracted data
//! \param[in] y The vector
//! \param[in] len The number of indices to extract
//! \param[in] start The first index in the range
static void extractBlock(Vector& x, const Vector& y, int len, int start=0)
{
x.resize(len,false);
std::copy(y.begin()+start,y.begin()+len+start,x.begin());
}
//! \brief Inject a range of indices into a vector
//! \param[in/out] x The vector to inject into
//! \param[in] y The vector with the data to inject
//! \param[in] len The number of indices to inject
//! \param[in] start The first index in the range
static void injectBlock(Vector& x, const Vector& y, int len, int start=0)
{
std::copy(y.begin(),y.begin()+len,x.begin()+start);
}
};
}
}

77
dune/elasticity/seqlu.hpp Normal file
View File

@ -0,0 +1,77 @@
//==============================================================================
//!
//! \file seqlu.hpp
//!
//! \date Dec 22 2012
//!
//! \author Arne Morten Kvarving / SINTEF
//!
//! \brief SuperLU preconditioner interface
//!
//==============================================================================
#pragma once
#if HAVE_SUPERLU
#include <dune/istl/superlu.hh>
#endif
#if HAVE_SUPERLU
template<class M, class X, class Y, int l=1>
class SeqLU : public Dune::Preconditioner<X,Y> {
public:
//! \brief The matrix type the preconditioner is for.
typedef typename Dune::remove_const<M>::type matrix_type;
//! \brief The domain type of the preconditioner.
typedef X domain_type;
//! \brief The range type of the preconditioner.
typedef Y range_type;
//! \brief The field type of the preconditioner.
typedef typename X::field_type field_type;
// define the category
enum {
//! \brief The category the preconditioner is part of.
category=Dune::SolverCategory::sequential
};
/*! \brief Constructor.
Constructor gets all parameters to operate the prec.
\param A The matrix to operate on.
*/
SeqLU (const M& A) :
slu(A, false)
{
}
/*!
\brief Prepare the preconditioner.
\copydoc Preconditioner::pre(X&,Y&)
*/
virtual void pre (X& x, Y& b) {}
/*!
\brief Apply the precondioner.
\copydoc Preconditioner::apply(X&,const Y&)
*/
virtual void apply (X& v, const Y& d)
{
Dune::InverseOperatorResult res;
Y t(d);
slu.apply(v, t, res);
}
/*!
\brief Clean up.
\copydoc Preconditioner::post(X&)
*/
virtual void post (X& x) {}
private:
Dune::SuperLU<M> slu;
};
#endif

56
dune/elasticity/uzawa_solver.hpp Normal file
View File

@ -0,0 +1,56 @@
#pragma once
namespace Opm {
namespace Elasticity {
template<class X, class Y>
class UzawaSolver : public Dune::InverseOperator<X,Y>
{
public:
UzawaSolver(Dune::InverseOperator<X,Y>* innersolver_,
Dune::InverseOperator<X,Y>* outersolver_,
const Matrix& B_) :
innersolver(innersolver_), outersolver(outersolver_), B(B_)
{
}
virtual ~UzawaSolver()
{
delete innersolver;
delete outersolver;
}
void apply(X& x, Y& b, double reduction,
Dune::InverseOperatorResult& res)
{
apply(x, b, res);
}
void apply(X& x, Y& b, Dune::InverseOperatorResult& res)
{
Vector lambda, lambda2, u, u2;
MortarUtils::extractBlock(u, b, B.N());
Dune::InverseOperatorResult res2;
u2.resize(u.size());
u2 = 0;
innersolver->apply(u2, u, res2);
lambda.resize(B.M());
B.mtv(u2, lambda);
lambda2.resize(lambda.size());
lambda2 = 0;
outersolver->apply(lambda2, lambda, res2);
MortarUtils::extractBlock(u, b, B.N());
B.usmv(-1.0, lambda2, u);
u2 = 0;
innersolver->apply(u2, u, res);
MortarUtils::injectBlock(x, u2, B.N());
MortarUtils::injectBlock(x, lambda2, B.M(), B.N());
}
protected:
Dune::InverseOperator<X,Y>* innersolver;
Dune::InverseOperator<X,Y>* outersolver;
const Matrix& B;
};
}
}

140
examples/upscale_elasticity.cpp
View File

@ -27,6 +27,7 @@
#endif
#include <dune/elasticity/elasticity_upscale.hpp>
#include <dune/elasticity/matrixops.hpp>
//! \brief Display the available command line parameters
@ -34,29 +35,36 @@ void syntax(char** argv)
{
std::cerr << "Usage: " << argv[0] << " gridfilename=filename.grdecl [method=]" << std::endl
<< "\t[xmax=] [ymax=] [zmax=] [xmin=] [ymin=] [zmin=] [linsolver_type=]" << std::endl
<<" \t[Emin=] [ctol=] [ltol=] [rock_list=] [vtufilename=] [output=]" << std::endl << std::endl
<< "\t gridfilename - the grid file. can be 'uniform'" << std::endl
<< "\t vtufilename - save results to vtu file" << std::endl
<< "\t rock_list - use material information from given rocklist" << std::endl
<< "\t output - output results in text report format" << std::endl
<< "\t method - can be MPC, LLM or Mortar" << std::endl
<<" \t[Emin=] [ctol=] [ltol=] [rock_list=] [vtufilename=] [output=] [linsolver_verbose=]" << std::endl << std::endl
<< "\t gridfilename - the grid file. can be 'uniform'" << std::endl
<< "\t vtufilename - save results to vtu file" << std::endl
<< "\t rock_list - use material information from given rocklist" << std::endl
<< "\t output - output results in text report format" << std::endl
<< "\t method - can be MPC or Mortar" << std::endl
<< "\t\t if not specified, Mortar couplings are used" << std::endl
<< "\t (x|y|z)max/(x|y|z)min - maximum/minimum grid coordinates." << std::endl
<< "\t\t if not specified, coordinates are found by grid traversal" << std::endl
<< "\t cells(x|y|z) - number of cells if gridfilename=uniform." << std::endl
<< "\t Emin - minimum E modulus (to avoid numerical issues)" << std::endl
<< "\t ctol - collapse tolerance in grid parsing" << std::endl
<< "\t ltol - tolerance in iterative linear solvers" << std::endl
<< "\t linsolver_type=cg - use the conjugate gradient method" << std::endl
<< "\t linsolver_type=slu - use the SuperLU sparse direct solver" << std::endl;
<< "\t cells(x|y|z) - number of cells if gridfilename=uniform." << std::endl
<< "\t lambda(x|y) - order of lagrangian multipliers in first/second direction" << std::endl
<< "\t Emin - minimum E modulus (to avoid numerical issues)" << std::endl
<< "\t ctol - collapse tolerance in grid parsing" << std::endl
<< "\t ltol - tolerance in iterative linear solvers" << std::endl
<< "\t linsolver_type=iterative - use a suitable iterative method (cg or gmres)" << std::endl
<< "\t linsolver_type=direct - use the SuperLU sparse direct solver" << std::endl
<< "\t verbose - set to true to get verbose output" << std::endl
<< "\t linsolver_restart - number of iterations before gmres is restarted" << std::endl
<< "\t linsolver_presteps - number of pre-smooth steps in the AMG" << std::endl
<< "\t linsolver_poststeps - number of post-smooth steps in the AMG" << std::endl
<< "\t\t affects memory usage" << std::endl
<< "\t linsolver_symmetric - use symmetric linear solver. Defaults to true" << std::endl
<< "\t mortar_precond - preconditioner for mortar block. Defaults to schur-amg" << std::endl;
}
enum UpscaleMethod {
UPSCALE_NONE = 0,
UPSCALE_MPC = 1,
UPSCALE_LLM = 2,
UPSCALE_MORTAR = 3
UPSCALE_MORTAR = 2
};
//! \brief Structure holding parameters configurable from command line
@ -75,28 +83,28 @@ struct Params {
double Emin;
//! \brief The tolerance for collapsing nodes in the Z direction
double ctol;
//! \brief The tolerance for the iterative linear solver
double ltol;
//! \brief Minimum grid coordinates
double min[3];
//! \brief Maximum grid coordinates
double max[3];
//! \brief The linear solver to employ
Opm::Elasticity::Solver solver;
//! \brief Polynomial order of multipliers in first / second direction
int lambda[2];
//! \brief Number of elements on interface grid in the x direction
int n1;
//! \brief Number of elements on interface grid in the y direction
int n2;
Opm::Elasticity::LinSolParams linsolver;
// Debugging options
//! \brief Number of elements on interface grid in the x direction (LLM)
int n1;
//! \brief Number of elements on interface grid in the y direction (LLM)
int n2;
//! \brief cell in x
int cellsx;
//! \brief cell in y
int cellsy;
//! \brief cell in z
int cellsz;
//! \brief Polynomial order of multipliers in first / second direction
int lambda[2];
//! \brief verbose output
bool verbose;
};
//! \brief Parse the command line arguments
@ -109,38 +117,32 @@ void parseCommandLine(int argc, char** argv, Params& p)
p.min[0] = param.getDefault("xmin",-1);
p.min[1] = param.getDefault("ymin",-1);
p.min[2] = param.getDefault("zmin",-1);
p.lambda[0] = param.getDefault("lambdax", 1);
p.lambda[1] = param.getDefault("lambday", 1);
p.lambda[0] = param.getDefault("lambdax", 2);
p.lambda[1] = param.getDefault("lambday", 2);
std::string method = param.getDefault<std::string>("method","mortar");
if (!strcasecmp(method.c_str(),"mpc"))
p.method = UPSCALE_MPC;
if (!strcasecmp(method.c_str(),"llm"))
p.method = UPSCALE_LLM;
if (!strcasecmp(method.c_str(),"mortar"))
p.method = UPSCALE_MORTAR;
if (!strcasecmp(method.c_str(),"none"))
p.method = UPSCALE_NONE;
p.Emin = param.getDefault<double>("Emin",0.0);
p.ctol = param.getDefault<double>("ctol",1.e-8);
p.ltol = param.getDefault<double>("ltol",1.e-10);
p.ctol = param.getDefault<double>("ctol",1.e-6);
p.file = param.get<std::string>("gridfilename");
p.rocklist = param.getDefault<std::string>("rock_list","");
p.vtufile = param.getDefault<std::string>("vtufilename","");
p.output = param.getDefault<std::string>("output","");
p.verbose = param.getDefault<bool>("verbose",false);
size_t i;
if ((i=p.vtufile.find(".vtu")) != std::string::npos)
p.vtufile = p.vtufile.substr(0,i);
//std::string solver = param.getDefault<std::string>("linsolver_type","slu");
std::string solver = param.getDefault<std::string>("linsolver_type","cg");
if (solver == "cg")
p.solver = Opm::Elasticity::CG;
else
p.solver = Opm::Elasticity::SLU;
if (p.file == "uniform") {
p.cellsx = param.getDefault("cellsx",3);
p.cellsy = param.getDefault("cellsy",3);
p.cellsz = param.getDefault("cellsz",3);
}
p.linsolver.parse(param);
p.n1 = -1;
p.n2 = -1;
}
@ -161,8 +163,6 @@ void writeOutput(const Params& p, Opm::time::StopWatch& watch, int cells,
gethostname(hostname,1024);
std::string method = "mortar";
if (p.method == UPSCALE_LLM)
method = "llm";
if (p.method == UPSCALE_MPC)
method = "mpc";
if (p.method == UPSCALE_NONE)
@ -194,10 +194,10 @@ void writeOutput(const Params& p, Opm::time::StopWatch& watch, int cells,
}
f << "# Options used:" << std::endl
<< "#\t method: " << method << std::endl
<< "#\t linsolver_type: " << (p.solver==Opm::Elasticity::SLU?"slu":"cg")
<< "#\t linsolver_type: " << (p.linsolver.type==Opm::Elasticity::DIRECT?"direct":"iterative")
<< std::endl;
if (p.solver == Opm::Elasticity::CG)
f << "#\t ltol: " << p.ltol << std::endl;
if (p.linsolver.type == Opm::Elasticity::ITERATIVE)
f << "#\t ltol: " << p.linsolver.tol << std::endl;
if (p.file == "uniform") {
f << "#\t cellsx: " << p.cellsx << std::endl
<< "#\t cellsy: " << p.cellsy << std::endl
@ -245,7 +245,10 @@ int main(int argc, char** argv)
grid.readEclipseFormat(p.file,p.ctol,false);
typedef GridType::ctype ctype;
Opm::Elasticity::ElasticityUpscale<GridType> upscale(grid, p.ctol, p.Emin, p.file, p.rocklist);
Opm::Elasticity::ElasticityUpscale<GridType> upscale(grid, p.ctol,
p.Emin, p.file,
p.rocklist,
p.verbose);
if (p.max[0] < 0 || p.min[0] < 0) {
std::cout << "determine side coordinates..." << std::endl;
upscale.findBoundaries(p.min,p.max);
@ -256,36 +259,54 @@ int main(int argc, char** argv)
p.n1 = grid.logicalCartesianSize()[0];
p.n2 = grid.logicalCartesianSize()[1];
}
if (p.linsolver.zcells == -1) {
double lz = p.max[2]-p.min[2];
int nz = grid.logicalCartesianSize()[2];
double hz = lz/nz;
double lp = sqrt((double)(p.max[0]-p.min[0])*(p.max[1]-p.min[1]));
int np = std::max(grid.logicalCartesianSize()[0],
grid.logicalCartesianSize()[1]);
double hp = lp/np;
p.linsolver.zcells = (int)(2*hp/hz+0.5);
}
std::cout << "logical dimension: " << grid.logicalCartesianSize()[0]
<< "x" << grid.logicalCartesianSize()[1]
<< "x" << grid.logicalCartesianSize()[2]
<< std::endl;
if (p.method == UPSCALE_LLM) {
std::cout << "using LLM couplings..." << std::endl;
upscale.periodicBCsLLM(p.min,p.max,p.n1,p.n2);
} else if (p.method == UPSCALE_MPC) {
if (p.method == UPSCALE_MPC) {
std::cout << "using MPC couplings in all directions..." << std::endl;
upscale.periodicBCs(p.min,p.max);
upscale.periodicBCs(p.min, p.max);
std::cout << "preprocessing grid..." << std::endl;
upscale.A.initForAssembly();
} else if (p.method == UPSCALE_MORTAR) {
std::cout << "using Mortar couplings.." << std::endl;
upscale.periodicBCsMortar(p.min,p.max,p.n1,p.n2,p.lambda[0], p.lambda[1]);
upscale.periodicBCsMortar(p.min, p.max, p.n1, p.n2,
p.lambda[0], p.lambda[1]);
} else if (p.method == UPSCALE_NONE) {
std::cout << "no periodicity approach applied.." << std::endl;
upscale.fixCorners(p.min, p.max);
upscale.A.initForAssembly();
}
Dune::FieldMatrix<double,6,6> C;
Dune::VTKWriter<GridType::LeafGridView> vtkwriter(grid.leafView());
Dune::VTKWriter<GridType::LeafGridView>* vtkwriter=0;
if (!p.vtufile.empty())
vtkwriter = new Dune::VTKWriter<GridType::LeafGridView>(grid.leafView());
Opm::Elasticity::Vector field[6];
std::cout << "assembling..." << "\n";
std::cout << "assembling elasticity operator..." << "\n";
upscale.assemble(-1,true);
upscale.setupSolvers(p.solver);
#pragma omp parallel for schedule(static)
std::cout << "setting up linear solver..." << std::endl;
upscale.setupSolvers(p.linsolver);
//#pragma omp parallel for schedule(static)
for (int i=0;i<6;++i) {
std::cout << "processing case " << i+1 << "..." << std::endl;
std::cout << "\tassembling load vector..." << std::endl;
upscale.assemble(i,false);
std::cout << "solving case " << i+1 << "..." << "\n";
upscale.solve(p.solver,p.ltol,i);
std::cout << "\tsolving..." << std::endl;
upscale.solve(i);
upscale.A.expandSolution(field[i],upscale.u[i]);
#define CLAMP(x) (fabs(x)<1.e-5?0.f:x)
#define CLAMP(x) (fabs(x)<1.e-4?0.0:x)
for (size_t j=0;j<field[i].size();++j) {
double val = field[i][j];
field[i][j] = CLAMP(val);
@ -298,10 +319,13 @@ int main(int argc, char** argv)
for (int i=0;i<6;++i) {
std::stringstream str;
str << "sol " << i+1;
vtkwriter.addVertexData(field[i], str.str().c_str(), dim);
if (vtkwriter)
vtkwriter->addVertexData(field[i], str.str().c_str(), dim);
}
if (vtkwriter) {
vtkwriter->write(p.vtufile);
delete vtkwriter;
}
if (!p.vtufile.empty())
vtkwriter.write(p.vtufile);
// voigt notation
for (int j=0;j<6;++j)
std::swap(C[3][j],C[5][j]);