Commit 1a0fc322 authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

port geomechanics models to stable. This completes FS#191. Big thanks to...

port geomechanics models to stable. This completes FS#191. Big thanks to Melanie for developing these models.

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@10903 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 0245e8c7
......@@ -38,6 +38,10 @@ AC_CONFIG_FILES([dumux.pc
dumux/freeflow/stokes/Makefile
dumux/freeflow/stokes2c/Makefile
dumux/freeflow/stokes2cni/Makefile
dumux/geomechanics/Makefile
dumux/geomechanics/elastic/Makefile
dumux/geomechanics/el1p2c/Makefile
dumux/geomechanics/el2p/Makefile
dumux/implicit/Makefile
dumux/implicit/1p/Makefile
dumux/implicit/1p2c/Makefile
......@@ -105,6 +109,10 @@ AC_CONFIG_FILES([dumux.pc
test/freeflow/stokes/Makefile
test/freeflow/stokes2c/Makefile
test/freeflow/stokes2cni/Makefile
test/geomechanics/Makefile
test/geomechanics/elastic/Makefile
test/geomechanics/el1p2c/Makefile
test/geomechanics/el2p/Makefile
test/material/Makefile
test/material/fluidsystems/Makefile
test/material/immiscibleflash/Makefile
......
......@@ -2,6 +2,7 @@ SUBDIRS = \
common \
decoupled \
freeflow \
geomechanics \
implicit \
io \
linear \
......
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \brief Defines a function to calculate eigenvalues of n x n matrices. For n > 2 a cyclic jacobi method is used.
* This implementation is not efficient for larger matrices!
*/
#ifndef DUMUX_EIGENVALUES_HH
#define DUMUX_EIGENVALUES_HH
#include <math.h>
namespace Dumux
{
template<class ValueType>
int sign(const ValueType& value)
{
return (value < 0 ? -1 : 1);
}
template<int dim, class Matrix>
void identityMatrix(Matrix& matrix)
{
for (int i = 0; i < dim; i++)
for (int j = 0; j < dim; j++)
matrix[i][j] = (i == j) ? 1.0 : 0.0;
}
template<int dim, class Matrix>
double calcOffDiagonalNorm(Matrix& matrix)
{
double norm = 0;
for (int i = 0; i < dim; i++)
for (int j = 0; j < dim; j++)
if (i != j)
{
norm += matrix[i][j] * matrix[i][j];
}
return std::sqrt(norm);
}
template<int dim, class EVVectorType, class MatrixType>
bool calculateEigenValues(EVVectorType &eigVel, MatrixType& matrix)
{
if (dim == 2)
{
eigVel = 0;
double b = -(matrix[0][0] + matrix[1][1]);
double c = matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0];
eigVel[0] = (-b + sqrt(b * b - 4.0 * c)) / 2.0;
eigVel[1] = (-b - sqrt(b * b - 4.0 * c)) / 2.0;
if (std::isnan(eigVel[0]) || std::isinf(eigVel[0]))
return false;
if (std::isnan(eigVel[1]) || std::isinf(eigVel[1]))
return false;
return true;
}
else if (dim > 2)
{
int maxIter = 100;
int iter = 0;
double matrixNorm = 0;
double offDiagonalNorm = 1;
MatrixType rotationMatrix(0.0);
MatrixType evMatrix(matrix);
while (iter < maxIter && offDiagonalNorm > 0.01 * matrixNorm)
{
for (int i = 0; i < dim - 1; i++)
{
for (int j = i + 1; j < dim; j++)
{
identityMatrix<dim>(rotationMatrix);
double theta = (evMatrix[i][i] - evMatrix[j][j])
/ (2 * evMatrix[i][j]);
double t = sign(theta)
/ (std::abs(theta) + std::sqrt(1 + theta * theta));
double c = 1 / std::sqrt(1 + t * t);
double s = c * t;
rotationMatrix[i][i] = c;
rotationMatrix[j][j] = c;
rotationMatrix[i][j] = s;
rotationMatrix[j][i] = -s;
evMatrix.leftmultiply(rotationMatrix);
rotationMatrix[i][j] = -s;
rotationMatrix[j][i] = s;
evMatrix.rightmultiply(rotationMatrix);
}
}
matrixNorm = evMatrix.frobenius_norm();
offDiagonalNorm = calcOffDiagonalNorm<dim>(evMatrix);
iter++;
}
for (int i = 0; i < dim; i++)
{
eigVel[i] = evMatrix[i][i];
if (std::isnan(eigVel[i]) || std::isinf(eigVel[i]))
return false;
}
return true;
}
else
{
return false;
}
}
;
}
#endif
SUBDIRS = elastic el1p2c el2p
geomechanicsdir = $(includedir)/dumux/geomechanics
include $(top_srcdir)/am/global-rules
el1p2cdir = $(includedir)/dumux/geomechanics/el1p2c
el1p2c_HEADERS = *.hh
include $(top_srcdir)/am/global-rules
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \brief Volume variables gathered on an element
*/
#ifndef DUMUX_BOX_EL1P2C_ELEMENT_VOLUME_VARIABLES_HH
#define DUMUX_BOX_EL1P2C_ELEMENT_VOLUME_VARIABLES_HH
#include <dumux/implicit/box/boxproperties.hh>
#include <dumux/implicit/box/boxelementvolumevariables.hh>
namespace Dumux
{
/*!
* \ingroup ElOnePTwoCBoxModel
*
* \brief This class stores an array of VolumeVariables objects, one
* volume variables object for each of the element's vertices
*/
template<class TypeTag>
class ElOnePTwoCElementVolumeVariables : public BoxElementVolumeVariables<TypeTag>
{
typedef BoxElementVolumeVariables<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, VertexMapper) VertexMapper;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GridView::template Codim<0>::Entity Element;
enum { dim = GridView::dimension };
typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
public:
/*!
* \brief The constructor.
*/
ElOnePTwoCElementVolumeVariables()
{ }
/*!
* \brief Construct the volume variables for all vertices of an element.
*
* \param problem The problem which needs to be simulated.
* \param element The DUNE Codim<0> entity for which the volume variables ought to be calculated
* \param fvGeometry The finite volume geometry of the element
* \param oldSol Tells whether the model's previous or current solution should be used.
*
* This class is required for the update of the effective porosity values at the
* vertices since it is a function of the divergence of the solid displacement
* at the integration points
*/
void update(const Problem &problem,
const Element &element,
const FVElementGeometry &fvGeometry,
bool oldSol)
{
ParentType::update(problem, element, fvGeometry, oldSol);
this->updateEffPorosity(problem, element, fvGeometry);
};
/*!
* \brief Update the effective porosities and the volumetric strain divU for all vertices of an element.
*
* \param problem The problem which needs to be simulated.
* \param element The DUNE Codim<0> entity for which the volume variables ought to be calculated
* \param fvGeometry The finite volume geometry of the element
*
* This function is required for the update of the effective porosity / divU values at the
* vertices.
*
* During the partial derivative calculation, changes of the solid displacement
* at vertex i can affect effective porosities / divU of all element vertices.
* To correctly update the effective porosities / divU of all element vertices
* an iteration over all scv faces is required.
* The remaining volvars are only updated for the vertex whose primary variable
* is changed for the derivative calculation.
*/
void updateEffPorosity(const Problem &problem,
const Element &element,
const FVElementGeometry &fvGeometry)
{
// we assert that the i-th shape function is
// associated to the i-th vert of the element.
int numScv = element.template count<dim>();
// number of faces which contribute to the porosity value in the sub-control volume
std::vector<double> numContributingFaces;
numContributingFaces.resize(numScv);
for (int scvIdx = 0; scvIdx < numScv; scvIdx++) {
(*this)[scvIdx].effPorosity = 0.0;
(*this)[scvIdx].divU = 0.0;
numContributingFaces[scvIdx] = 0.0;
}
for (int faceIdx = 0; faceIdx < fvGeometry.numScvf; faceIdx++)
{
// evaluate the gradients at the IPs for each subcontrol volume face
FluxVariables fluxVars(problem,
element,
fvGeometry,
faceIdx,
*this);
numContributingFaces[fluxVars.face().i] += 1;
numContributingFaces[fluxVars.face().j] += 1;
// average value for the effective porosity
(*this)[fluxVars.face().i].effPorosity += fluxVars.effPorosity();
(*this)[fluxVars.face().j].effPorosity += fluxVars.effPorosity();
// average value for the volumetric strain
(*this)[fluxVars.face().i].divU += fluxVars.divU();
(*this)[fluxVars.face().j].divU += fluxVars.divU();
}
for (int scvIdx = 0; scvIdx < numScv; scvIdx++) {
(*this)[scvIdx].effPorosity /= numContributingFaces[scvIdx];
(*this)[scvIdx].divU /= numContributingFaces[scvIdx];
}
};
};
} // namespace Dumux
#endif
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
*
* \brief This file contains the calculation of all the fluxes over the surface of the
* finite volume that make up the volume, the mass and the momentum balance
* for the one-phase two-component linear-elastic model.
*
* This means pressure, concentration and solid-displacement gradients, phase densities at
* the integration point, etc.
*
* This class inherits from the one-phase two-component model FluxVariables and from the
* linear elasticity model FluxVariables
*/
#ifndef DUMUX_ELASTIC1P2C_FLUX_VARIABLES_HH
#define DUMUX_ELASTIC1P2C_FLUX_VARIABLES_HH
#include <dumux/geomechanics/elastic/elasticfluxvariables.hh>
#include <dumux/implicit/1p2c/1p2cfluxvariables.hh>
namespace Dumux
{
/*!
* \ingroup ElOnePTwoCBoxModel
* \ingroup ImplicitFluxVariables
* \brief This template class contains the data which is required to
* calculate the fluxes over the surface of the
* finite volume that make up the volume, the mass and the momentum balance
* for the one-phase two-component linear-elastic model.
*
* This means pressure, concentration and solid-displacement gradients, phase densities at
* the integration point, etc.
*
*/
template<class TypeTag>
class ElOnePTwoCFluxVariables: public ElasticFluxVariablesBase<TypeTag> ,
public OnePTwoCFluxVariables<TypeTag>
{
typedef ElasticFluxVariablesBase<TypeTag> ElasticBase;
typedef OnePTwoCFluxVariables<TypeTag> OnePTwoCBase;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GridView::template Codim<0>::Entity Element;
enum
{
dim = GridView::dimension,
dimWorld = GridView::dimensionworld
};
typedef typename GridView::ctype CoordScalar;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef Dune::FieldVector<CoordScalar, dim> DimVector;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
public:
/*
* \brief The constructor
*
* \param problem The problem
* \param element The finite element
* \param fvGeometry The finite-volume geometry in the fully implicit scheme
* \param faceIdx The local index of the SCV (sub-control-volume) face
* \param elemVolVars The volume variables of the current element
* \param onBoundary A boolean variable to specify whether the flux variables
* are calculated for interior SCV faces or boundary faces, default=false
*/
ElOnePTwoCFluxVariables(const Problem &problem,
const Element &element,
const FVElementGeometry &fvGeometry,
int faceIdx,
const ElementVolumeVariables &elemVolVars,
const bool onBoundary = false)
: ElasticBase(problem, element, fvGeometry, faceIdx, elemVolVars),
OnePTwoCBase(problem, element, fvGeometry, faceIdx, elemVolVars),
fvGeometry_(fvGeometry), faceIdx_(faceIdx), onBoundary_(onBoundary)
{
dU_ = Scalar(0);
dGradP_ = Scalar(0);
porosity_ = 0;
effPorosity_ = 0;
pressure_ = 0;
timeDerivUNormal_ = 0;
elOnePTwoCGradients_(problem, element, elemVolVars);
calculateEffectiveValues_(problem, element, elemVolVars);
calculateDiffCoeffPM_(problem, element, elemVolVars);
calculateDDt_(problem, element, elemVolVars);
}
;
public:
/*!
* \brief Return porosity [-] at the integration point.
*/
Scalar porosity() const
{
return porosity_;
}
/*!
* \brief Return effective porosity [-] at the integration point.
*/
Scalar effPorosity() const
{
return effPorosity_;
}
/*!
* \brief Return pressure [Pa] at the integration
* point.
*/
Scalar pressure() const
{
return pressure_;
}
/*!
* \brief Return change of pressure gradient with time [Pa/m] at
* integration point.
*/
Scalar dGradP(int dimIdx) const
{
return dGradP_[dimIdx];
}
/*!
* \brief Return gradient of time derivative of pressure [Pa].
*/
Scalar timeDerivGradPNormal() const
{
return timeDerivGradPNormal_;
}
/*!
* \brief Return change of u [m] with time at integration point
* point.
*/
Scalar dU(int dimIdx) const
{
return dU_[dimIdx];
}
/*!
* \brief Return time derivative of u [m/s] in normal direction at integration point
*/
Scalar timeDerivUNormal() const
{
return timeDerivUNormal_;
}
/*!
* \brief Return porous medium diffusion coefficient [m^2]
*/
Scalar diffCoeffPM() const
{
return diffCoeffPM_;
}
const SCVFace &face() const
{
return fvGeometry_.subContVolFace[faceIdx_];
}
protected:
/*!
* \brief Calculation of the solid displacement and pressure gradients.
*
* \param problem The considered problem file
* \param element The considered element of the grid
* \param elemVolVars The parameters stored in the considered element
*/
void elOnePTwoCGradients_(const Problem &problem,
const Element &element,
const ElementVolumeVariables &elemVolVars)
{
// calculate gradients
GlobalPosition tmp(0.0);
for (int idx = 0; idx < fvGeometry_.numScv; idx++) // loop over adjacent vertices
{
// FE gradient at vertex idx
const DimVector &feGrad = face().grad[idx];
// the gradient of the temporal pressure change (needed for stabilization term)
tmp = feGrad;
tmp *= elemVolVars[idx].dPressure();
dGradP_ += tmp;
// average the pressure at integration point
pressure_ += elemVolVars[idx].pressure()
* face().shapeValue[idx];
// average temporal displacement change at integration point (for calculation of solid displacement velocity)
for (int i = 0; i < dim; ++i)
dU_[i] += elemVolVars[idx].dU(i)
* face().shapeValue[idx];
// average porosity at integration point
porosity_ += elemVolVars[idx].porosity()
* face().shapeValue[idx];
}
}
/*!
* \brief Calculation of the effective porosity.
*
* \param problem The considered problem file
* \param element The considered element of the grid
* \param elemVolVars The parameters stored in the considered element
*/
void calculateEffectiveValues_(const Problem &problem,
const Element &element,
const ElementVolumeVariables &elemVolVars)
{
// the effective porosity is calculated as a function of solid displacement and initial porosity
// according to Han & Dusseault (2003)
// calculate effective porosity as a function of solid displacement and initial porosity
effPorosity_ = (porosity_ + this->divU())
/ (1 + this->divU());
}
/*!
* \brief Calculation of the effective porous media diffusion coefficient.
*
* \param problem The considered problem file
* \param element The considered element of the grid
* \param elemVolVars The parameters stored in the considered element
*/
void calculateDiffCoeffPM_(const Problem &problem,
const Element &element,
const ElementVolumeVariables &elemVolVars)
{
const VolumeVariables &volVarsI = elemVolVars[face().i];
const VolumeVariables &volVarsJ = elemVolVars[face().j];
// Effective diffusion coefficient in the porous medium
diffCoeffPM_ = effPorosity_ * harmonicMean(volVarsI.tortuosity() * volVarsI.diffCoeff(),
volVarsJ.tortuosity() * volVarsJ.diffCoeff());
}
/*!