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

Merge pull request #19 from akva2/high_order_multipliers

Browse Source 
Elasticity upscaling - autumn update
This commit is contained in:
Bård Skaflestad 2013-01-07 06:32:54 -08:00
commit ffadae0f34
28 changed files with 2116 additions and 615 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
dune/elasticity/Makefile.am
View File

@ -2,6 +2,7 @@ lib_LTLIBRARIES = libelasticityupscale.la
libelasticityupscale_la_SOURCES = \
boundarygrid.cpp \
dynmatrixev.cpp \
material.cpp \
materials.cpp \
mpc.cpp

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

@ -0,0 +1,50 @@
//==============================================================================
//!
//! \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
//!
//==============================================================================
#ifndef APPLIER_H_
#define APPLIER_H_
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);
}
}
}
#endif

18
dune/elasticity/asmhandler.hpp
View File

@ -9,15 +9,14 @@
//! \brief Class handling finite element assembly
//!
//==============================================================================
#pragma once
//! \brief A vector holding our RHS
typedef Dune::BlockVector<Dune::FieldVector<double,1> > Vector;
#ifndef ASMHANDLER_HPP_
#define ASMHANDLER_HPP_
#include <dune/geometry/referenceelements.hh>
#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>
@ -58,7 +57,7 @@ class ASMHandler {
//! \brief Get the number of equations in the system
//! \returns The number of equations in the system
int getEqns() const
size_t getEqns() const
{
return maxeqn;
}
@ -142,6 +141,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()
@ -228,3 +234,5 @@ class ASMHandler {
}
}
#endif

44
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>
@ -59,7 +61,7 @@ void ASMHandler<GridType>::addDOF(int row, int erow,
for (int l=0;l<dim;++l) {
MPC* mpc = getMPC(index2,l);
if (mpc) {
for (int n=0;n<mpc->getNoMaster();++n) {
for (size_t n=0;n<mpc->getNoMaster();++n) {
int idx = meqn[mpc->getMaster(n).node*dim+mpc->getMaster(n).dof-1];
if (idx != -1)
A[row][idx] += scale*mpc->getMaster(n).coeff*(*K)[erow][j*dim+l];
@ -70,7 +72,7 @@ void ASMHandler<GridType>::addDOF(int row, int erow,
}
}
}
if (S)
if (S && b)
(*b)[row] += scale*(*S)[erow];
}
@ -93,7 +95,7 @@ void ASMHandler<GridType>::addElement(
for (int k=0;k<dim;++k) {
MPC* mpc = getMPC(index1,k);
if (mpc) {
for (int n=0;n<mpc->getNoMaster();++n) {
for (size_t n=0;n<mpc->getNoMaster();++n) {
int idx = meqn[mpc->getMaster(n).node*dim+mpc->getMaster(n).dof-1];
addDOF(idx,i*dim+k,K,S,set,cell,b2,mpc->getMaster(n).coeff);
}
@ -124,7 +126,7 @@ void ASMHandler<GridType>::extractValues(Dune::FieldVector<double,comp>& v,
if (it2 != fixedNodes.end() && it2->second.first & (1 << n))
v[l++] = it2->second.second[n];
else if (mpc) {
for (int m=0;m<mpc->getNoMaster();++m) {
for (size_t m=0;m<mpc->getNoMaster();++m) {
int idx = meqn[mpc->getMaster(m).node*dim+mpc->getMaster(m).dof-1];
if (idx != -1)
v[l] += u[idx]*mpc->getMaster(m).coeff;
@ -143,7 +145,7 @@ void ASMHandler<GridType>::expandSolution(Vector& result, const Vector& u)
result.resize(nodes*dim);
result = 0;
int l=0;
for (size_t i=0;i<nodes;++i) {
for (int i=0;i<nodes;++i) {
fixIt it = fixedNodes.find(i);
Direction dir;
if (it == fixedNodes.end())
@ -163,11 +165,11 @@ void ASMHandler<GridType>::expandSolution(Vector& result, const Vector& u)
}
// second loop - handle MPC couplings
l = 0;
for (size_t i=0;i<nodes;++i) {
for (int i=0;i<nodes;++i) {
for (int j=0;j<dim;++j) {
MPC* mpc = getMPC(i,j);
if (mpc) {
for (int n=0;n<mpc->getNoMaster();++n) {
for (size_t n=0;n<mpc->getNoMaster();++n) {
int idx = mpc->getMaster(n).node*dim+mpc->getMaster(n).dof-1;
if (meqn[idx] != -1)
result[l] += u[meqn[idx]]*mpc->getMaster(n).coeff;
@ -188,7 +190,7 @@ void ASMHandler<GridType>::addMPC(MPC* mpc)
fixIt it = fixedNodes.find(mpc->getSlave().node);
int flag = 1 << (mpc->getSlave().dof-1);
if (it == fixedNodes.end() ||
!(it->second.first & (1 << (mpc->getSlave().dof-1)))) {
!(it->second.first & flag)) {
mpcs.insert(std::make_pair(slaveNode,mpc));
return;
}
@ -257,7 +259,7 @@ void ASMHandler<GridType>::resolveMPCChain(MPC* mpc)
coeff[master.dof-1] = master.coeff;
int removeOld = 0;
for (size_t d = 1; d <= dim; d++)
for (int d = 1; d <= dim; d++)
if (fabs(coeff[d-1]) > 1.0e-8)
{
MPC* mpc2 = getMPC(mpc->getMaster(i).node,d-1);
@ -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>
@ -354,7 +356,7 @@ void ASMHandler<GridType>::nodeAdjacency(const LeafIterator& it,
for (int l=0;l<dim;++l) {
MPC* mpc = getMPC(indexj,l);
if (mpc) {
for (int i=0;i<mpc->getNoMaster();++i) {
for (size_t i=0;i<mpc->getNoMaster();++i) {
int idx = meqn[mpc->getMaster(i).node*dim+
mpc->getMaster(i).dof-1];
if (idx != -1)
@ -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);
@ -386,7 +391,7 @@ void ASMHandler<GridType>::determineAdjacencyPattern()
for (int k=0;k<dim;++k) {
MPC* mpc = getMPC(indexi,k);
if (mpc) {
for (int l=0;l<mpc->getNoMaster();++l) {
for (size_t l=0;l<mpc->getNoMaster();++l) {
nodeAdjacency(it,vertexsize,
meqn[mpc->getMaster(l).node*dim+
mpc->getMaster(l).dof-1]);
@ -395,5 +400,6 @@ void ASMHandler<GridType>::determineAdjacencyPattern()
nodeAdjacency(it,vertexsize,meqn[indexi*dim+k]);
}
}
help.log(cell, "\t\t... still processing ... cell ");
}
}

8
dune/elasticity/boundarygrid.cpp
View File

@ -175,7 +175,7 @@ int BoundaryGrid::Q4inv(FaceCoord& res, const Quad& q,
// check that obtained solutions are inside element
double tol = 1+epsOut;
int nInside=0;
for (int i=0;i<xi.size();++i) {
for (size_t i=0;i<xi.size();++i) {
if (xi[i] < tol && eta[i] < tol) {
if (++nInside > 1) {
std::cout << "multiple solutions" << std::endl;
@ -232,7 +232,7 @@ bool BoundaryGrid::bilinearSolve(double epsilon, double order,
double Q2 = A[0]*B[1]-B[0]*A[1];
std::vector<double> Z;
cubicSolve(epsilon,0,Q2,Q1,Q0,Z);
for (int i=0;i<Z.size();++i) {
for (size_t i=0;i<Z.size();++i) {
Q0 = A[0]*Z[i]+A[2];
if (fabs(Q0) > tol) {
X.push_back(Z[i]);
@ -244,10 +244,10 @@ bool BoundaryGrid::bilinearSolve(double epsilon, double order,
Q2 = A[0]*B[2]-B[0]*A[2];
Z.clear();
cubicSolve(epsilon,0,Q2,Q1,Q0,Z);
for (int i=0;i<Z.size();++i) {
for (size_t i=0;i<Z.size();++i) {
Q0 = A[0]*Z[i]+A[1];
if (fabs(Q0) > tol) {
int j=0;
size_t j=0;
for (j=0;j<Y.size();++j)
if (fabs(Y[j]-Z[i]) <= epsilon*order) break;
if (j == Y.size()) {

28
dune/elasticity/boundarygrid.hh
View File

@ -9,7 +9,8 @@
//! \brief Class describing 2D quadrilateral grids
//!
//==============================================================================
#pragma once
#ifndef BOUNDARYGRID_HH_
#define BOUNDARYGRID_HH_
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
@ -113,6 +114,13 @@ class BoundaryGrid {
//! \param[in] quad The quad to add
void add(const Quad& quad);
void addToColumn(size_t col, const Quad& quad)
{
if (col >= colGrids.size())
colGrids.resize(col+1);
colGrids[col].push_back(quad);
}
//! \brief Obtain a reference to a quad
//! \param[in] index The index of the requested quad
//! \returns A reference to the requested quad
@ -129,12 +137,22 @@ class BoundaryGrid {
return grid[index];
}
const Quad& getQuad(int col, int index) const
{
return colGrids[col][index];
}
//! \brief Obtain the number of quads in the grid
size_t size() const
{
return grid.size();
}
size_t colSize(int i) const
{
return colGrids[i].size();
}
//! \brief Return the total number of nodes on the grid when known
//! \sa uniform
size_t totalNodes() const
@ -209,6 +227,7 @@ class BoundaryGrid {
protected:
//! \brief Our quadrilateral elements
std::vector<Quad> grid;
std::vector<std::vector<Quad> > colGrids;
//! \brief Whether or not a given node is marked as fixed
std::vector<bool> fixNodes;
@ -219,7 +238,7 @@ class BoundaryGrid {
//! \brief Print to a stream
friend std::ostream& operator <<(std::ostream& os, const BoundaryGrid& g)
{
for (int i=0;i<g.size();++i)
for (size_t i=0;i<g.size();++i)
os << g[i] << std::endl;
return os;
}
@ -308,7 +327,7 @@ class HexGeometry<2, cdim, GridImp>
const typename GridImp::LeafGridView::template Codim<3>::Iterator itend = gv.leafView().template end<3>();
for (int i=0;i<4;++i) {
typename GridImp::LeafGridView::template Codim<3>::Iterator it=start;
for (it ; it != itend; ++it) {
for (; it != itend; ++it) {
if (mapper.map(*it) == q.v[i].i)
break;
}
@ -396,7 +415,6 @@ class HexGeometry<2, cdim, GridImp>
Dune::GenericReferenceElements< ctype, 2 >::general(type());
LocalCoordinate x = refElement.position(0,0);
LocalCoordinate dx;
int i=0;
do {
using namespace Dune::GenericGeometry;
// DF^n dx^n = F^n, x^{n+1} -= dx^n
@ -483,3 +501,5 @@ BoundaryGrid::Vertex minXmaxY(std::vector<BoundaryGrid::Vertex>& in);
}
}
#endif

152
dune/elasticity/dynmatrixev.cpp Normal file
View File

@ -0,0 +1,152 @@
#ifndef DUNE_FMATRIXEIGENVALUES_EXT_CC
#define DUNE_FMATRIXEIGENVALUES_EXT_CC
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <iostream>
#include <cmath>
#include <cassert>
#include <dune/common/exceptions.hh>
#if HAVE_LAPACK
// nonsymetric matrices
#define DGEEV_FORTRAN FC_FUNC (dgeev, DGEEV)
// dsyev declaration (in liblapack)
extern "C" {
/*
*
** purpose
** =======
**
** xgeev computes for an N-by-N BASE DATA TYPE nonsymmetric matrix A, the
** eigenvalues and, optionally, the left and/or right eigenvectors.
**
** The right eigenvector v(j) of A satisfies
** A * v(j) = lambda(j) * v(j)
** where lambda(j) is its eigenvalue.
** The left eigenvector u(j) of A satisfies
** u(j)**T * A = lambda(j) * u(j)**T
** where u(j)**T denotes the transpose of u(j).
**
** The computed eigenvectors are normalized to have Euclidean norm
** equal to 1 and largest component real.
**
** arguments
** =========
**
** jobvl (input) char
** = 'n': left eigenvectors of a are not computed;
** = 'v': left eigenvectors of a are computed.
**
** jobvr (input) char
** = 'n': right eigenvectors of a are not computed;
** = 'v': right eigenvectors of a are computed.
**
** n (input) long int
** the order of the matrix v. v >= 0.
**
** a (input/output) BASE DATA TYPE array, dimension (lda,n)
** on entry, the n-by-n matrix a.
** on exit, a has been overwritten.
**
** lda (input) long int
** the leading dimension of the array a. lda >= max(1,n).
**
** wr (output) BASE DATA TYPE array, dimension (n)
** wi (output) BASE DATA TYPE array, dimension (n)
** wr and wi contain the real and imaginary parts,
** respectively, of the computed eigenvalues. complex
** conjugate pairs of eigenvalues appear consecutively
** with the eigenvalue having the positive imaginary part
** first.
**
** vl (output) COMPLEX DATA TYPE array, dimension (ldvl,n)
** if jobvl = 'v', the left eigenvectors u(j) are stored one
** after another in the columns of vl, in the same order
** as their eigenvalues.
** if jobvl = 'n', vl is not referenced.
** if the j-th eigenvalue is real, then u(j) = vl(:,j),
** the j-th column of vl.
** if the j-th and (j+1)-st eigenvalues form a complex
** conjugate pair, then u(j) = vl(:,j) + i*vl(:,j+1) and
** u(j+1) = vl(:,j) - i*vl(:,j+1).
**
** ldvl (input) long int
** the leading dimension of the array vl. ldvl >= 1; if
** jobvl = 'v', ldvl >= n.
**
** vr (output) COMPLEX DATA TYPE array, dimension (ldvr,n)
** if jobvr = 'v', the right eigenvectors v(j) are stored one
** after another in the columns of vr, in the same order
** as their eigenvalues.
** if jobvr = 'n', vr is not referenced.
** if the j-th eigenvalue is real, then v(j) = vr(:,j),
** the j-th column of vr.
** if the j-th and (j+1)-st eigenvalues form a complex
** conjugate pair, then v(j) = vr(:,j) + i*vr(:,j+1) and
** v(j+1) = vr(:,j) - i*vr(:,j+1).
**
** ldvr (input) long int
** the leading dimension of the array vr. ldvr >= 1; if
** jobvr = 'v', ldvr >= n.
**
** work (workspace/output) BASE DATA TYPE array, dimension (max(1,lwork))
** on exit, if info = 0, work(1) returns the optimal lwork.
**
** lwork (input) long int
** the dimension of the array work. lwork >= max(1,3*n), and
** if jobvl = 'v' or jobvr = 'v', lwork >= 4*n. for good
** performance, lwork must generally be larger.
**
** if lwork = -1, then a workspace query is assumed; the routine
** only calculates the optimal size of the work array, returns
** this value as the first entry of the work array, and no error
** message related to lwork is issued by xerbla.
**
** info (output) long int
** = 0: successful exit
** < 0: if info = -i, the i-th argument had an illegal value.
** > 0: if info = i, the qr algorithm failed to compute all the
** eigenvalues, and no eigenvectors have been computed;
** elements i+1:n of wr and wi contain eigenvalues which
** have converged.
**
**/
extern void DGEEV_FORTRAN(const char* jobvl, const char* jobvr, const long
int* n, double* a, const long int* lda, double* wr, double* wi, double* vl,
const long int* ldvl, double* vr, const long int* ldvr, double* work,
const long int* lwork, const long int* info);
} // end extern C
#endif
namespace Dune {
namespace DynamicMatrixHelp {
void eigenValuesNonsymLapackCall(
const char* jobvl, const char* jobvr, const long
int* n, double* a, const long int* lda, double* wr, double* wi, double* vl,
const long int* ldvl, double* vr, const long int* ldvr, double* work,
const long int* lwork, const long int* info)
{
#if HAVE_LAPACK
// call LAPACK dgeev
DGEEV_FORTRAN(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr,
work, lwork, info);
#else
DUNE_THROW(NotImplemented,"eigenValuesNonsymLapackCall: LAPACK not found!");
#endif
}
} // end namespace FMatrixHelp
} // end namespace Dune
#endif

67
dune/elasticity/dynmatrixev.hh Normal file
View File

@ -0,0 +1,67 @@
#ifndef DUNE_DYNMATRIXEIGENVALUES_HH
#define DUNE_DYNMATRIXEIGENVALUES_HH
#include <dune/common/dynmatrix.hh>
namespace Dune {
namespace DynamicMatrixHelp {
// defined in fmatrixev_ext.cpp
extern void eigenValuesNonsymLapackCall(
const char* jobvl, const char* jobvr, const long
int* n, double* a, const long int* lda, double* wr, double* wi, double* vl,
const long int* ldvl, double* vr, const long int* ldvr, double* work,
const long int* lwork, const long int* info);
template <typename K, class C>
static void eigenValuesNonSym(const DynamicMatrix<K>& matrix,
DynamicVector<C>& eigenValues)
{
{
const long int N = matrix.rows();
const char jobvl = 'n';
const char jobvr = 'n';
// matrix to put into dgeev
double matrixVector[N * N];
// copy matrix
int row = 0;
for(int i=0; i<N; ++i)
{
for(int j=0; j<N; ++j, ++row)
{
matrixVector[ row ] = matrix[ i ][ j ];
}
}
// working memory
double eigenR[N];
double eigenI[N];
double work[3*N];
// return value information
long int info = 0;
long int lwork = 3*N;
// call LAPACK routine (see fmatrixev_ext.cc)
eigenValuesNonsymLapackCall(&jobvl, &jobvr, &N, &matrixVector[0], &N,
&eigenR[0], &eigenI[0], 0, &N, 0, &N, &work[0],
&lwork, &info);
if( info != 0 )
{
std::cerr << "For matrix " << matrix << " eigenvalue calculation failed! " << std::endl;
DUNE_THROW(InvalidStateException,"eigenValues: Eigenvalue calculation failed!");
}
for (int i=0;i<N;++i)
eigenValues[i] = std::complex<double>(eigenR[i], eigenI[i]);
}
}
}
}
#endif

5
dune/elasticity/elasticity.hpp
View File

@ -9,7 +9,8 @@
//! \brief Elasticity helper class
//!
//==============================================================================
#pragma once
#ifndef ELASTICITY_HPP_
#define ELASTICITY_HPP_
namespace Opm {
namespace Elasticity {
@ -69,3 +70,5 @@ class Elasticity {
}
}
#endif

6
dune/elasticity/elasticity_impl.hpp
View File

@ -30,7 +30,7 @@ void Elasticity<GridType>::getBmatrix(Dune::FieldMatrix<ctype,components,funcdim
if (dim == 3) {
#define INDEX(i,j) i+6*j
for (size_t i=0; i < funcs; ++i) {
for (int i=0; i < funcs; ++i) {
// normal strain part
temp[INDEX(0,0)][i] = dNdX[i][0];
temp[INDEX(1,1)][i] = dNdX[i][1];
@ -48,7 +48,7 @@ void Elasticity<GridType>::getBmatrix(Dune::FieldMatrix<ctype,components,funcdim
} else if (dim == 2) {
#undef INDEX
#define INDEX(i,j) i+3*j
for (size_t i=0; i < funcs; ++i) {
for (int i=0; i < funcs; ++i) {
// normal strain part
temp[INDEX(0,0)][i] = dNdX[i][0];
temp[INDEX(1,1)][i] = dNdX[i][1];
@ -59,7 +59,7 @@ void Elasticity<GridType>::getBmatrix(Dune::FieldMatrix<ctype,components,funcdim
}
}
size_t k=0;
for (size_t j=0;j<funcs*dim;++j)
for (int j=0;j<funcs*dim;++j)
for (size_t i=0;i<components;++i,++k)
B[i][j] = temp[k%rows][k/rows];
}

271
dune/elasticity/elasticity_upscale.hpp
View File

@ -9,7 +9,8 @@
//! \brief Elasticity upscale class
//!
//==============================================================================
#pragma once
#ifndef ELASTICITY_UPSCALE_HPP_
#define ELASTICITY_UPSCALE_HPP_
#include <dune/common/fmatrix.hh>
#include <opm/core/eclipse/EclipseGridParser.hpp>
@ -17,30 +18,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 +153,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 +185,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,22 +211,21 @@ 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
//! \param[in] n1 The number of elements on the lambda grid in the X direction
//! \param[in] n2 The number of elements on the lambda grid in the Y direction
//! \param[in] p1 The order of multipliers in the X direction
//! \param[in] p2 The order of multipliers in the Y direction
void periodicBCsMortar(const double* min,
const double* max, int n1, int n2);
const double* max, int n1, int n2,
int p1, int p2);
//! \brief Fix corner nodes
//! \param[in] min The minimum coordinates on the grid
//! \param[in] max The maximum coordinates on the grid
void fixCorners(const double* min, const double* max);
//! \brief Assemble (optionally) stiffness matrix A and load vector
//! \param[in] loadcase The strain load case. Set to -1 to skip
@ -168,15 +241,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;
@ -210,11 +284,6 @@ class ElasticityUpscale
//! \brief Vector holding material parameters for each active grid cell
std::vector<Material*> materials;
//! \brief Fix corner nodes
//! \param[in] min The minimum coordinates on the grid
//! \param[in] max The maximum coordinates on the grid
void fixCorners(const double* min, const double* max);
//! \brief Extract the vertices on a given face
//! \param[in] dir The direction of the face normal
//! \param[in] coord The coordinate of the face plane
@ -268,23 +337,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)
Matrix findLMatrixMortar(const BoundaryGrid& b1,
const BoundaryGrid& interface, int dir);
//! \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 (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
@ -296,9 +371,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;
@ -306,31 +384,66 @@ 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"
}
}
#endif

727
dune/elasticity/elasticity_upscale_impl.hpp
View File

File diff suppressed because it is too large Load Diff

152
dune/elasticity/fmatrixev_ext.cc Normal file
View File

@ -0,0 +1,152 @@
#ifndef DUNE_FMATRIXEIGENVALUES_EXT_CC
#define DUNE_FMATRIXEIGENVALUES_EXT_CC
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <iostream>
#include <cmath>
#include <cassert>
#include <dune/common/exceptions.hh>
#if HAVE_LAPACK
// nonsymetric matrices
#define DGEEV_FORTRAN FC_FUNC (dgeev, DGEEV)
// dsyev declaration (in liblapack)
extern "C" {
/*
*
** purpose
** =======
**
** xgeev computes for an N-by-N BASE DATA TYPE nonsymmetric matrix A, the
** eigenvalues and, optionally, the left and/or right eigenvectors.
**
** The right eigenvector v(j) of A satisfies
** A * v(j) = lambda(j) * v(j)
** where lambda(j) is its eigenvalue.
** The left eigenvector u(j) of A satisfies
** u(j)**T * A = lambda(j) * u(j)**T
** where u(j)**T denotes the transpose of u(j).
**
** The computed eigenvectors are normalized to have Euclidean norm
** equal to 1 and largest component real.
**
** arguments
** =========
**
** jobvl (input) char
** = 'n': left eigenvectors of a are not computed;
** = 'v': left eigenvectors of a are computed.
**
** jobvr (input) char
** = 'n': right eigenvectors of a are not computed;
** = 'v': right eigenvectors of a are computed.
**
** n (input) long int
** the order of the matrix v. v >= 0.
**
** a (input/output) BASE DATA TYPE array, dimension (lda,n)
** on entry, the n-by-n matrix a.
** on exit, a has been overwritten.
**
** lda (input) long int
** the leading dimension of the array a. lda >= max(1,n).
**
** wr (output) BASE DATA TYPE array, dimension (n)
** wi (output) BASE DATA TYPE array, dimension (n)
** wr and wi contain the real and imaginary parts,
** respectively, of the computed eigenvalues. complex
** conjugate pairs of eigenvalues appear consecutively
** with the eigenvalue having the positive imaginary part
** first.
**
** vl (output) COMPLEX DATA TYPE array, dimension (ldvl,n)
** if jobvl = 'v', the left eigenvectors u(j) are stored one
** after another in the columns of vl, in the same order
** as their eigenvalues.
** if jobvl = 'n', vl is not referenced.
** if the j-th eigenvalue is real, then u(j) = vl(:,j),
** the j-th column of vl.
** if the j-th and (j+1)-st eigenvalues form a complex
** conjugate pair, then u(j) = vl(:,j) + i*vl(:,j+1) and
** u(j+1) = vl(:,j) - i*vl(:,j+1).
**
** ldvl (input) long int
** the leading dimension of the array vl. ldvl >= 1; if
** jobvl = 'v', ldvl >= n.
**
** vr (output) COMPLEX DATA TYPE array, dimension (ldvr,n)
** if jobvr = 'v', the right eigenvectors v(j) are stored one
** after another in the columns of vr, in the same order
** as their eigenvalues.
** if jobvr = 'n', vr is not referenced.
** if the j-th eigenvalue is real, then v(j) = vr(:,j),
** the j-th column of vr.
** if the j-th and (j+1)-st eigenvalues form a complex
** conjugate pair, then v(j) = vr(:,j) + i*vr(:,j+1) and
** v(j+1) = vr(:,j) - i*vr(:,j+1).
**
** ldvr (input) long int
** the leading dimension of the array vr. ldvr >= 1; if
** jobvr = 'v', ldvr >= n.
**
** work (workspace/output) BASE DATA TYPE array, dimension (max(1,lwork))
** on exit, if info = 0, work(1) returns the optimal lwork.
**
** lwork (input) long int
** the dimension of the array work. lwork >= max(1,3*n), and
** if jobvl = 'v' or jobvr = 'v', lwork >= 4*n. for good
** performance, lwork must generally be larger.
**
** if lwork = -1, then a workspace query is assumed; the routine
** only calculates the optimal size of the work array, returns
** this value as the first entry of the work array, and no error
** message related to lwork is issued by xerbla.
**
** info (output) long int
** = 0: successful exit
** < 0: if info = -i, the i-th argument had an illegal value.
** > 0: if info = i, the qr algorithm failed to compute all the
** eigenvalues, and no eigenvectors have been computed;
** elements i+1:n of wr and wi contain eigenvalues which
** have converged.
**
**/
extern void DGEEV_FORTRAN(const char* jobvl, const char* jobvr, const long
int* n, double* a, const long int* lda, double* wr, double* wi, double* vl,
const long int* ldvl, double* vr, const long int* ldvr, double* work,
const long int* lwork, const long int* info);
} // end extern C
#endif
namespace Dune {
namespace FMatrixHelp {
void eigenValuesNonsymLapackCall(
const char* jobvl, const char* jobvr, const long
int* n, double* a, const long int* lda, double* wr, double* wi, double* vl,
const long int* ldvl, double* vr, const long int* ldvr, double* work,
const long int* lwork, const long int* info)
{
#if HAVE_LAPACK
// call LAPACK dgeev
DGEEV_FORTRAN(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr,
work, lwork, info);
#else
DUNE_THROW(NotImplemented,"eigenValuesNonsymLapackCall: LAPACK not found!");
#endif
}
} // end namespace FMatrixHelp
} // end namespace Dune
#endif

70
dune/elasticity/fmatrixev_ext.hh Normal file
View File

@ -0,0 +1,70 @@
#ifndef DUNE_FMATRIXEIGENVALUES_EXT_HH
#define DUNE_FMATRIXEIGENVALUES_EXT_HH
namespace Dune {
namespace FMatrixHelp {
// defined in fmatrixev_ext.cpp
extern void eigenValuesNonsymLapackCall(
const char* jobvl, const char* jobvr, const long
int* n, double* a, const long int* lda, double* wr, double* wi, double* vl,
const long int* ldvl, double* vr, const long int* ldvr, double* work,
const long int* lwork, const long int* info);
template <int dim, typename K, class C>
static void eigenValuesNonSym(const FieldMatrix<K, dim, dim>& matrix,
FieldVector<C, dim>& eigenValues)
{
{
const long int N = dim ;
const char jobvl = 'n';
const char jobvr = 'n';
// length of matrix vector
const long int w = N * N ;
// matrix to put into dgeev
double matrixVector[dim * dim];
// copy matrix
int row = 0;
for(int i=0; i<dim; ++i)
{
for(int j=0; j<dim; ++j, ++row)
{
matrixVector[ row ] = matrix[ i ][ j ];
}
}
// working memory
double eigenR[dim];
double eigenI[dim];
double work[3*dim];
// return value information
long int info = 0;
long int lwork = 3*dim;
// call LAPACK routine (see fmatrixev_ext.cc)
eigenValuesNonsymLapackCall(&jobvl, &jobvr, &N, &matrixVector[0], &N,
&eigenR[0], &eigenI[0], 0, &N, 0, &N, &work[0],
&lwork, &info);
if( info != 0 )
{
std::cerr << "For matrix " << matrix << " eigenvalue calculation failed! " << std::endl;
DUNE_THROW(InvalidStateException,"eigenValues: Eigenvalue calculation failed!");
}
for (int i=0;i<N;++i) {
eigenValues[i].real = eigenR[i];
eigenValues[i].imag = eigenI[i];
}
}
}
}
}
#endif

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

@ -0,0 +1,34 @@
//==============================================================================
//!
//! \file logutils.hpp
//!
//! \date Nov 9 2011
//!
//! \author Arne Morten Kvarving / SINTEF
//!
//! \brief Logging helper utilities
//!
//==============================================================================
#ifndef LOGUTILS_HPP_
#define LOGUTILS_HPP_
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
};
#endif

5
dune/elasticity/material.hh
View File

@ -9,7 +9,8 @@
//! \brief Material interface.
//!
//==============================================================================
#pragma once
#ifndef MATERIAL_HH_
#define MATERIAL_HH_
#include <dune/common/fmatrix.hh>
#include <dune/common/dynvector.hh>
@ -97,3 +98,5 @@ private:
}
}
#endif

5
dune/elasticity/materials.hh
View File

@ -9,7 +9,8 @@
//! \brief Material properties.
//!
//==============================================================================
#pragma once
#ifndef MATERIALS_HH_
#define MATERIALS_HH_
#include "material.hh"
@ -173,3 +174,5 @@ private:
}
}
#endif

34
dune/elasticity/matrixops.hpp
View File

@ -9,7 +9,8 @@
//! \brief Helper class with some matrix operations
//!
//==============================================================================
#pragma once
#ifndef MATRIXOPS_HPP_
#define MATRIXOPS_HPP_
#include <dune/istl/bcrsmatrix.hh>
@ -22,6 +23,9 @@ typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> > Matrix;
//! \brief For storing matrix adjacency/sparsity patterns
typedef std::vector< std::set<int> > AdjacencyPattern;
//! \brief A vector holding our RHS
typedef Dune::BlockVector<Dune::FieldVector<double,1> > Vector;
//! \brief Helper class with some matrix operations
class MatrixOps {
public:
@ -33,6 +37,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);
@ -52,9 +61,32 @@ class MatrixOps {
//! \param[in] symmetric If true, augment symmetrically
static Matrix augment(const Matrix& A, const Matrix& B,
size_t r0, size_t c0, bool symmetric);
//! \brief Extract the diagonal of a matrix into a new matrix
//! \param[in] The matrix to extract the diagonal from
//! \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
//! \details This is only useful for debugging as the files grow very big
static void saveAsc(const Matrix& A, const std::string& file);
};
#include "matrixops_impl.hpp"
}
}
#endif

125
dune/elasticity/matrixops_impl.hpp
View File

@ -12,14 +12,14 @@
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();
A.setSize(rows, cols, sum);
A.setBuildMode(Matrix::random);
for (int i = 0; i < adj.size(); i++)
for (size_t i = 0; i < adj.size(); i++)
A.setrowsize(i,adj[i].size());
A.endrowsizes();
@ -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();
@ -53,7 +77,7 @@ Matrix MatrixOps::Axpy(const Matrix& A, const Matrix& B, double alpha)
assert(A.M() == B.M() && A.N() == B.N());
// establish union adjacency pattern
std::vector<std::set<int> > adj;
AdjacencyPattern adj;
adj.resize(A.N());
for (Matrix::ConstIterator it = A.begin();
it != A.end(); ++it) {
@ -120,7 +144,7 @@ Matrix MatrixOps::augment(const Matrix& A, const Matrix& B,
}
if (symmetric) {
// always establish diagonal elements or superLU crashes
for (int i=0;i<nrow;++i)
for (size_t i=0;i<nrow;++i)
adj[i].insert(i);
}
Matrix result;
@ -144,3 +168,96 @@ Matrix MatrixOps::augment(const Matrix& A, const Matrix& B,
return result;
}
Matrix MatrixOps::extractDiagonal(const Matrix& A)
{
AdjacencyPattern adj;
adj.resize(A.M());
for (size_t i=0;i<A.M();++i)
adj[i].insert(i);
Matrix result;
fromAdjacency(result,adj,A.M(),A.M());
for (size_t i=0;i<A.M();++i)
result[i][i] = A[i][i];
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;
f.open(file.c_str());
f << "% " << A.N() << " " << A.M() << std::endl;
int prevrow=-1;
for (Matrix::ConstIterator it = A.begin();
it != A.end(); ++it) {
for (int i=0;i<int(it.index())-prevrow-1;++i) {
for (size_t j=0;j<A.M();++j)
f << "0 ";
f << std::endl;
}
int prevcol=-1;
for (Matrix::ConstColIterator it2 = it->begin();
it2 != it->end();++it2) {
for (int j=0;j<int(it2.index())-prevcol-1;++j)
f << "0 ";
double val = *it2;
f << val << " ";
prevcol = it2.index();
}
for (int j=0;j<int(A.M())-prevcol-1;++j)
f << "0 ";
prevrow = it.index();
f << std::endl;
}
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;
}

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

@ -0,0 +1,82 @@
//==============================================================================
//!
//! \file mortar_evaluator.hpp
//!
//! \date Nov 15 2012
//!
//! \author Arne Morten Kvarving / SINTEF
//!
//! \brief Linear operator for a Mortar block
//!
//==============================================================================
#ifndef MORTAR_EVALUATOR_HPP_
#define MORTAR_EVALUATOR_HPP_
#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
};
}
}
#endif

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

@ -0,0 +1,76 @@
//==============================================================================
//!
//! \file mortar_schur.hpp
//!
//! \date Nov 15 2012
//!
//! \author Arne Morten Kvarving / SINTEF
//!
//! \brief Schur based linear operator for a Mortar block
//!
//==============================================================================
#ifndef MORTAR_SCHUR_HPP_
#define MORTAR_SCHUR_HPP_
#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
};
}
}
#endif

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

@ -0,0 +1,126 @@
//==============================================================================
//!
//! \file schur_precond.hpp
//!
//! \date Nov 15 2012
//!
//! \author Arne Morten Kvarving / SINTEF
//!
//! \brief Schur based preconditioner for saddle point systems
//!
//==============================================================================
#ifndef SCHUR_PRECOND_HPP_
#define SCHUR_PRECOND_HPP_
#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), 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;
};
}
}
#endif

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

@ -0,0 +1,45 @@
//==============================================================================
//!
//! \file mortar_utils.hpp
//!
//! \date Nov 9 2012
//!
//! \author Arne Morten Kvarving / SINTEF
//!
//! \brief Mortar helper class
//!
//==============================================================================
#ifndef MORTAR_UTILS_HPP_
#define MORTAR_UTILS_HPP_
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);
}
};
}
}
#endif

5
dune/elasticity/mpc.hh
View File

@ -9,7 +9,8 @@
//! \brief Representation of multi-point constraint (MPC) equations.
//!
//==============================================================================
#pragma once
#ifndef MPC_HH_
#define MPC_HH_
#include <map>
#include <vector>
@ -159,3 +160,5 @@ typedef std::map<int,MPC*> MPCMap;
}
}
#endif

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

@ -0,0 +1,79 @@
//==============================================================================
//!
//! \file seqlu.hpp
//!
//! \date Dec 22 2012
//!
//! \author Arne Morten Kvarving / SINTEF
//!
//! \brief SuperLU preconditioner interface
//!
//==============================================================================
#ifndef SEQLU_HPP_
#define SEQLU_HPP_
#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
#endif

250
dune/elasticity/shapefunctions.hpp
View File

@ -6,12 +6,16 @@
//!
//! \author Arne Morten Kvarving / SINTEF
//!
//! \brief Class for linear shape functions. Loosely based on code in dune-grid-howto
//! \brief Classes for shape functions. Loosely based on code in dune-grid-howto
//!
//==============================================================================
#pragma once
#ifndef SHAPEFUNCTIONS_HPP_
#define SHAPEFUNCTIONS_HPP_
#include <dune/common/fvector.hh>
#include "dynmatrixev.hh"
#include <complex>
namespace Opm {
namespace Elasticity {
@ -80,6 +84,96 @@ class LinearShapeFunction
Dune::FieldVector<rtype,dim> coeff1;
};
//! \brief Represents a cardinal function on a line
template<class ctype, class rtype>
class LagrangeCardinalFunction
{
public:
//! \brief Empty default constructor
LagrangeCardinalFunction() {}
//! \brief Construct a cardinal function with the given nodes
//! \param[in] nodes_ The nodes
//! \param[in] i The node this function is associated with
LagrangeCardinalFunction(const std::vector<rtype>& nodes_,
size_t i)
: nodes(nodes_), node(i) {}
//! \brief Evaluate the shape function
//! \param[in] local The local coordinates
rtype evaluateFunction(const ctype& local) const
{
rtype result = 1;
for (size_t i=0; i < nodes.size(); ++i) {
if (i != node)
result *= (local-nodes[i])/(nodes[node]-nodes[i]);
}
return result;
}
//! \brief Evaluate the derivative of the cardinal function
//! \param[in] local The local coordinates
rtype evaluateGradient(const ctype& local) const
{
rtype result = 0;
for (size_t i=0; i < nodes.size(); ++i) {
rtype f = 1;
for (int j=0; j < nodes.size(); ++j) {
if (i != j && j != node)
f *= (local-nodes[j])/(nodes[node]-nodes[j]);
}
result += f/(nodes[node]-nodes[i]);
}
return result;
}
private:
//! \brief The nodes
std::vector<rtype> nodes;
size_t node;
};
//! \brief Represents a tensor-product of 1D functions
template<class rtype, class ctype, class ftype, int dim>
class TensorProductFunction
{
public:
//! \brief The dimension of the function
enum { dimension = dim };
//! \brief Empty default constructor
TensorProductFunction() {}
//! \brief Construct a tensor-product function
//! \param[in] funcs_ The functions
TensorProductFunction(const Dune::FieldVector<ftype, dim>& funcs_)
: funcs(funcs_) {}
//! \brief Evaluate the function
//! \param[in] local The local coordinates
rtype evaluateFunction(const Dune::FieldVector<ctype,dim>& local) const
{
rtype result = 1;
for (int i=0; i < dim; ++i)
result *= funcs[i].evaluateFunction(local[i]);
return result;
}
Dune::FieldVector<rtype, dim>
evaluateGradient(const Dune::FieldVector<ctype,dim>& local) const
{
Dune::FieldVector<rtype, dim> result;
for (int i=0; i < dim; ++i)
result[i] = funcs[i].evaluateGradient(local[i]);
}
private:
Dune::FieldVector<ftype, dim> funcs;
};
//! \brief Singleton handler for the set of LinearShapeFunction
template<class ctype, class rtype, int dim>
class P1ShapeFunctionSet
@ -121,8 +215,8 @@ private:
0, 1,
1, 0,
0, 0};
static rtype coeffs22[] = {-1,-1,
1,-1,
static rtype coeffs22[] = {-1,-1,
1,-1,
-1, 1,
1, 1};
@ -134,8 +228,8 @@ private:
0, 1, 0,
1, 0, 0,
0, 0, 0};
static rtype coeffs32[] = {-1,-1,-1,
1,-1,-1,
static rtype coeffs32[] = {-1,-1,-1,
1,-1,-1,
-1, 1,-1,
1, 1,-1,
-1,-1, 1,
@ -173,5 +267,149 @@ private:
ShapeFunction f[n];
};
template<int dim>
class PNShapeFunctionSet
{
public:
typedef LagrangeCardinalFunction<double, double> CardinalFunction;
typedef TensorProductFunction<double, double, CardinalFunction, dim>
ShapeFunction;
PNShapeFunctionSet(int n1, int n2, int n3=0)
{
int dims[3] = {n1, n2, n3};
cfuncs.resize(dim);
for (int i=0; i < dim; ++i) {
std::vector<double> grid;
grid = gaussLobattoLegendreGrid(dims[i]);
for (int j=0;j<dims[i];++j)
cfuncs[i].push_back(CardinalFunction(grid,j));
}
int l=0;
Dune::FieldVector<CardinalFunction,dim> fs;
if (dim == 3) {
f.resize(n1*n2*n3);
for (int k=0; k < n3; ++k) {
for (int j=0; j < n2; ++j)
for (int i=0; i< n1; ++i) {
fs[0] = cfuncs[0][i];
fs[1] = cfuncs[1][j];
fs[2] = cfuncs[2][k];
f[l++] = ShapeFunction(fs);
}
}
} else {
f.resize(n1*n2);
for (int j=0; j < n2; ++j) {
for (int i=0; i< n1; ++i) {
fs[0] = cfuncs[0][i];
fs[1] = cfuncs[1][j];
f[l++] = ShapeFunction(fs);
}
}
}
}
//! \brief Obtain a given shape function
//! \param[in] i The requested shape function
const ShapeFunction& operator[](int i) const
{
return f[i];
}
int size()
{
return f.size();
}
protected:
std::vector< std::vector<CardinalFunction> > cfuncs;
std::vector<ShapeFunction> f;
double legendre(double x, int n)
{
std::vector<double> Ln;
Ln.resize(n+1);
Ln[0] = 1.f;
Ln[1] = x;
if( n > 1 ) {
for( int i=1;i<n;i++ )
Ln[i+1] = (2*i+1.0)/(i+1.0)*x*Ln[i]-i/(i+1.0)*Ln[i-1];
}
return Ln[n];
}
double legendreDerivative(double x, int n)
{
std::vector<double> Ln;
Ln.resize(n+1);
Ln[0] = 1.0; Ln[1] = x;
if( (x == 1.0) || (x == -1.0) )
return( pow(x,n-1)*n*(n+1.0)/2.0 );
else {
for( int i=1;i<n;i++ )
Ln[i+1] = (2.0*i+1.0)/(i+1.0)*x*Ln[i]-(double)i/(i+1.0)*Ln[i-1];
return( (double)n/(1.0-x*x)*Ln[n-1]-n*x/(1-x*x)*Ln[n] );
}
}
std::vector<double> gaussLegendreGrid(int n)
{
Dune::DynamicMatrix<double> A(n,n,0.0);
A[0][1] = 1.f;
for (int i=1;i<n-1;++i) {
A[i][i-1] = (double)i/(2.0*(i+1.0)-1.0);
A[i][i+1] = (double)(i+1.0)/(2*(i+1.0)-1.0);
}
A[n-1][n-2] = (n-1.0)/(2.0*n-1.0);
Dune::DynamicVector<std::complex<double> > eigenValues(n);
Dune::DynamicMatrixHelp::eigenValuesNonSym(A, eigenValues);
std::vector<double> result(n);
for (int i=0; i < n; ++i)
result[i] = std::real(eigenValues[i]);
std::sort(result.begin(),result.begin()+n);
return result;
}
std::vector<double> gaussLobattoLegendreGrid(int n)
{
assert(n > 1);
const double tolerance = 1.e-15;
std::vector<double> result(n);
result[0] = 0.0;
result[n-1] = 1.0;
if (n == 3)
result[1] = 0.5;
if (n < 4)
return result;
std::vector<double> glgrid = gaussLegendreGrid(n-1);
for (int i=1;i<n-1;++i) {
result[i] = (glgrid[i-1]+glgrid[i])/2.0;
double old = 0.0;
while (std::abs(old-result[i]) > tolerance) {
old = result[i];
double L = legendre(old, n-1);
double Ld = legendreDerivative(old, n-1);
result[i] += (1.0-old*old)*Ld/((n-1.0)*n*L);
}
result[i] = (result[i]+1.0)/2.0;
}
return result;
}
};
}
}
#endif

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

@ -0,0 +1,70 @@
#ifndef UZAWA_SOLVER_HPP_
#define UZAWA_SOLVER_HPP_
//==============================================================================
//!
//! \file uzawa_solver.hpp
//!
//! \date Nov 9 2012
//!
//! \author Arne Morten Kvarving / SINTEF
//!
//! \brief Uzawa scheme helper class
//!
//==============================================================================
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;
};
}
}
#endif

196
examples/upscale_elasticity.cpp
View File

@ -17,6 +17,7 @@
#include <iostream>
#include <unistd.h>
#include <cstring>
#include <dune/common/exceptions.hh> // We use exceptions
#include <opm/core/utility/StopWatch.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
@ -27,6 +28,7 @@
#endif
#include <dune/elasticity/elasticity_upscale.hpp>
#include <dune/elasticity/matrixops.hpp>
//! \brief Display the available command line parameters
@ -34,23 +36,38 @@ 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_MORTAR = 2
};
//! \brief Structure holding parameters configurable from command line
struct Params {
//! \brief The eclipse grid file
@ -61,74 +78,72 @@ struct Params {
std::string vtufile;
//! \brief Text output file
std::string output;
//! \brief Whether or not to use LLM couplings
bool LLM;
//! \brief Whether or not to use Mortar couplings
bool mortar;
//! \brief Whether or not to use MPC couplings
bool mpc;
//! \brief Method
UpscaleMethod method;
//! \brief A scaling parameter for the E-modulus to avoid numerical issues
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 verbose output
bool verbose;
};
//! \brief Parse the command line arguments
void parseCommandLine(int argc, char** argv, Params& p)
{
Opm::parameter::ParameterGroup param(argc, argv);
p.max[0] = param.getDefault("xmax",-1);
p.max[1] = param.getDefault("ymax",-1);
p.max[2] = param.getDefault("zmax",-1);
p.min[0] = param.getDefault("xmin",-1);
p.min[1] = param.getDefault("ymin",-1);
p.min[2] = param.getDefault("zmin",-1);
//std::string method = param.getDefault<std::string>("method","mpc");
p.max[0] = param.getDefault("xmax",-1);
p.max[1] = param.getDefault("ymax",-1);
p.max[2] = param.getDefault("zmax",-1);
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", 2);
p.lambda[1] = param.getDefault("lambday", 2);
std::string method = param.getDefault<std::string>("method","mortar");
p.LLM = (strcasecmp(method.c_str(),"llm") == 0);
//p.mortar = (strcasecmp(method.c_str(),"mortar") == 0);
p.mpc = (strcasecmp(method.c_str(),"mpc") == 0);
p.Emin = param.getDefault<double>("Emin",0.f);
p.ctol = param.getDefault<double>("ctol",1.e-8);
p.ltol = param.getDefault<double>("ltol",1.e-10);
if (method == "mpc")
p.method = UPSCALE_MPC;
if (method == "mortar")
p.method = UPSCALE_MORTAR;
if (method == "none")
p.method = UPSCALE_NONE;
p.Emin = param.getDefault<double>("Emin",0.0);
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;
}
@ -149,10 +164,10 @@ void writeOutput(const Params& p, Opm::time::StopWatch& watch, int cells,
gethostname(hostname,1024);
std::string method = "mortar";
if (p.LLM)
method = "llm";
if (p.mpc)
if (p.method == UPSCALE_MPC)
method = "mpc";
if (p.method == UPSCALE_NONE)
method = "none";
// write log
std::ofstream f;
@ -180,10 +195,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
@ -191,7 +206,7 @@ void writeOutput(const Params& p, Opm::time::StopWatch& watch, int cells,
}
f << "#" << std::endl
<<"# Materials: " << volume.size() << std::endl;
for (int i=0;i<volume.size();++i)
for (size_t i=0;i<volume.size();++i)
f << "#\t Material" << i+1 << ": " << volume[i]*100 << "%" << std::endl;
f << "#" << std::endl
<< "######################################################################" << std::endl
@ -203,13 +218,12 @@ int main(int argc, char** argv)
{
try {
static const int dim = 3;
static const int bfunc = 4+(dim-2)*4;
typedef Dune::CpGrid GridType;
if (argc < 2 || strcasecmp(argv[1],"-h") == 0
|| strcasecmp(argv[1],"--help") == 0
|| strcasecmp(argv[1],"-?") == 0) {
if (argc < 2 || strcmp(argv[1],"-h") == 0
|| strcmp(argv[1],"--help") == 0
|| strcmp(argv[1],"-?") == 0) {
syntax(argv);
exit(1);
}
@ -232,7 +246,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);
@ -243,33 +260,55 @@ 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.LLM) {
std::cout << "using LLM couplings..." << std::endl;
upscale.periodicBCsLLM(p.min,p.max,p.n1,p.n2);
} else if (p.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 {
} else if (p.method == UPSCALE_MORTAR) {
std::cout << "using Mortar couplings.." << std::endl;
upscale.periodicBCsMortar(p.min,p.max,p.n1,p.n2);
}
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());
Vector field[6];
std::cout << "assembling..." << "\n";
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 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 << "..." << "\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)
for (int j=0;j<field[i].size();++j) {
#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);
}
@ -281,10 +320,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]);