Skip to content
Snippets Groups Projects
Commit 418b9c4f authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

[implicit] generalize velocity output 

Before, velocity output only worked for cubes. This generalizes it to
simplices (box and cc) and prisms/pyramids (box only). Implements
FS#248.

Some reference solutions are adpated due to an improved boundary
condition handling for cell centered.

Reviewed by Alex.



git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@14085 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 70364c71
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dumux/implicit/common/implicitvelocityoutput.hh
+ 182
− 48
View file @ 418b9c4f
......@@ -19,12 +19,13 @@
/*!
* \file
*
* \brief Adaption of the BOX scheme to the two-phase two-component flow model.
* \brief velocity output for implicit (porous media) models
*/
#ifndef DUMUX_IMPLICIT_VELOCITYOUTPUT_HH
#define DUMUX_IMPLICIT_VELOCITYOUTPUT_HH
#include "implicitproperties.hh"
#include "implicitporousmediaproblem.hh"
#include <unordered_map>
#include <dune/common/fvector.hh>
......@@ -47,12 +48,14 @@ class ImplicitVelocityOutput
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GridView::ctype CoordScalar;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::IntersectionIterator IntersectionIterator;
typedef typename GridView::Intersection Intersection;
enum {
dim = GridView::dimension,
......@@ -63,6 +66,7 @@ class ImplicitVelocityOutput
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef typename Dune::ReferenceElements<CoordScalar, dim> ReferenceElements;
typedef typename Dune::ReferenceElement<CoordScalar, dim> ReferenceElement;
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
enum { dofCodim = isBox ? dim : 0 };
......@@ -78,41 +82,31 @@ public:
{
// check, if velocity output can be used (works only for cubes so far)
velocityOutput_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddVelocity);
if (velocityOutput_ && isBox)
{
cellNum_.assign(problem_.gridView().size(dofCodim), 0);
}
ElementIterator eIt = problem_.gridView().template begin<0>();
ElementIterator eEndIt = problem_.gridView().template end<0>();
for (; eIt != eEndIt; ++eIt)
if (velocityOutput_)
{
if (eIt->geometry().type().isCube() == false)
{
velocityOutput_ = false;
}
if (velocityOutput_ && isBox)
if (isBox)
{
FVElementGeometry fvGeometry;
fvGeometry.update(problem_.gridView(), *eIt);
cellNum_.assign(problem_.gridView().size(dofCodim), 0);
// transform vertex velocities from local to global coordinates
for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
ElementIterator eIt = problem_.gridView().template begin<0>();
ElementIterator eEndIt = problem_.gridView().template end<0>();
for (; eIt != eEndIt; ++eIt)
{
FVElementGeometry fvGeometry;
fvGeometry.update(problem_.gridView(), *eIt);
for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
{
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
int vIdxGlobal = problem_.vertexMapper().subIndex(*eIt, scvIdx, dofCodim);
int vIdxGlobal = problem_.vertexMapper().subIndex(*eIt, scvIdx, dofCodim);
#else
int vIdxGlobal = problem_.vertexMapper().map(*eIt, scvIdx, dofCodim);
int vIdxGlobal = problem_.vertexMapper().map(*eIt, scvIdx, dofCodim);
#endif
cellNum_[vIdxGlobal] += 1;
cellNum_[vIdxGlobal] += 1;
}
}
}
}
if (velocityOutput_ != GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddVelocity))
std::cout << "ATTENTION: Velocity output only works for cubes and is set to false for simplices\n";
}
bool enableOutput()
......@@ -120,20 +114,46 @@ public:
return velocityOutput_;
}
// The following SFINAE enable_if usage allows compilation, even if only a
//
// boundaryTypes(BoundaryTypes&, const Vertex&)
//
// is provided in the problem file. In that case, the compiler cannot detect
// (without additional measures like "using...") the signature
//
// boundaryTypes(BoundaryTypes&, const Intersection&)
//
// in the problem base class. Therefore, calls to this method trigger a
// compiler error. However, that call is needed for calculating velocities
// if the cell-centered discretization is used. By proceeding as in the
// following lines, that call will only be compiled if cell-centered
// actually is used.
template <class T = TypeTag>
void problemBoundaryTypes(BoundaryTypes& bcTypes,
const typename std::enable_if<!GET_PROP_VALUE(T, ImplicitIsBox), Intersection>::type& intersection) const
{
problem_.boundaryTypes(bcTypes, intersection);
}
template <class T = TypeTag>
void problemBoundaryTypes(BoundaryTypes& bcTypes,
const typename std::enable_if<GET_PROP_VALUE(T, ImplicitIsBox), Intersection>::type& intersection) const
{}
template<class VelocityVector>
void calculateVelocity(VelocityVector& velocity,
const ElementVolumeVariables& elemVolVars,
const FVElementGeometry& fvGeometry,
const Element& element,
const Element& element,
int phaseIdx)
{
if (velocityOutput_)
{
// calculate vertex velocities
GlobalPosition tmpVelocity(0.0);
Dune::GeometryType geomType = element.geometry().type();
const ReferenceElement &referenceElement
= ReferenceElements::general(geomType);
const Dune::FieldVector<Scalar, dim>& localPos =
ReferenceElements::general(element.geometry().type()).position(0, 0);
const Dune::FieldVector<Scalar, dim>& localPos
= referenceElement.position(0, 0);
// get the transposed Jacobian of the element mapping
const typename Element::Geometry::JacobianTransposed jacobianT2 =
......@@ -141,8 +161,8 @@ public:
if (isBox)
{
typedef Dune::BlockVector<Dune::FieldVector<Scalar, dim> > SCVVelocities;
SCVVelocities scvVelocities(pow(2,dim));
typedef Dune::BlockVector<Dune::FieldVector<Scalar, dim> > ScvVelocities;
ScvVelocities scvVelocities(fvGeometry.numScv);
scvVelocities = 0;
for (int fIdx = 0; fIdx < fvGeometry.numScvf; fIdx++)
......@@ -164,17 +184,17 @@ public:
GlobalPosition localNormal(0);
jacobianT1.mv(globalNormal, localNormal);
// note only works for cubes
const Scalar localArea = pow(2,-(dim-1));
localNormal /= localNormal.two_norm();
// Get the Darcy velocities. The Darcy velocities are divided by the area of the subcontrolvolume
// face in the reference element.
// area of the subcontrolvolume face in the reference element
Scalar localArea = scvfReferenceArea_(geomType, fIdx);
// get the volume flux divided by the area of the
// subcontrolvolume face in the reference element
Scalar flux = fluxVars.volumeFlux(phaseIdx) / localArea;
// transform the normal Darcy velocity into a vector
tmpVelocity = localNormal;
// transform the volume flux into a velocity vector
GlobalPosition tmpVelocity = localNormal;
tmpVelocity *= flux;
scvVelocities[fluxVars.face().i] += tmpVelocity;
......@@ -200,7 +220,11 @@ public:
}
else
{
Dune::FieldVector<Scalar, 2*dim> scvVelocities(0.0);
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
std::vector<Scalar> scvfFluxes(element.subEntities(1), 0);
#else
std::vector<Scalar> scvfFluxes(element.template count<1>(), 0);
#endif
int fIdxInner = 0;
IntersectionIterator isEndIt = problem_.gridView().iend(element);
......@@ -217,8 +241,7 @@ public:
fIdxInner,
elemVolVars);
Scalar flux = fluxVars.volumeFlux(phaseIdx);
scvVelocities[fIdx] += flux;
scvfFluxes[fIdx] += fluxVars.volumeFlux(phaseIdx);
fIdxInner++;
}
......@@ -230,14 +253,68 @@ public:
fIdx,
elemVolVars,true);
Scalar flux = fluxVars.volumeFlux(phaseIdx);
scvVelocities[fIdx] = flux;
scvfFluxes[fIdx] = fluxVars.volumeFlux(phaseIdx);
}
}
Dune::FieldVector < Scalar, dim > refVelocity(0);
for (int i = 0; i < dim; i++)
refVelocity[i] = 0.5 * (scvVelocities[2*i + 1] - scvVelocities[2*i]);
// Correct boundary fluxes in case of Neumann conditions.
// They are simply set to an average of the inner fluxes. In
// this general setting, it would be very difficult to
// calculate correct phase, i.e., volume, fluxes from arbitrary
// Neumann conditions.
if (element.hasBoundaryIntersections())
{
IntersectionIterator isEndIt = problem_.gridView().iend(element);
for (IntersectionIterator isIt = problem_.gridView().ibegin(element);
isIt != isEndIt; ++isIt)
{
if (isIt->boundary())
{
BoundaryTypes bcTypes;
problemBoundaryTypes(bcTypes, *isIt);
if (bcTypes.hasNeumann())
{
int fIdx = isIt->indexInInside();
// cubes
if (dim == 1 || geomType.isCube()){
int fIdxOpposite = fIdx%2 ? fIdx-1 : fIdx+1;
scvfFluxes[fIdx] = -scvfFluxes[fIdxOpposite];
}
// simplices
else if (geomType.isSimplex()) {
scvfFluxes[fIdx] = 0;
}
}
}
}
}
Dune::FieldVector < Scalar, dim > refVelocity;
// cubes
if (dim == 1 || geomType.isCube()){
for (int i = 0; i < dim; i++)
{
refVelocity[i] = 0.5 * (scvfFluxes[2*i + 1] - scvfFluxes[2*i]);
}
}
// simplices: Raviart-Thomas-0 interpolation evaluated at the cell center
else if (geomType.isSimplex()) {
for (int dimIdx = 0; dimIdx < dim; dimIdx++)
{
refVelocity[dimIdx] = -scvfFluxes[dim - 1 - dimIdx];
for (int fIdx = 0; fIdx < dim + 1; fIdx++)
{
refVelocity[dimIdx] += scvfFluxes[fIdx]/(dim + 1);
}
}
}
// 3D prism and pyramids
else {
DUNE_THROW(Dune::NotImplemented,
"velocity output for cell-centered and prism/pyramid");
}
Dune::FieldVector<Scalar, dimWorld> scvVelocity(0);
jacobianT2.mtv(refVelocity, scvVelocity);
......@@ -255,6 +332,63 @@ public:
} // velocity output
}
private:
// The area of a subcontrolvolume face in a reference element.
// The 3d non-cube values have been calculated with quadrilateralArea3D
// of boxfvelementgeometry.hh.
static Scalar scvfReferenceArea_(Dune::GeometryType geomType, int fIdx)
{
if (dim == 1 || geomType.isCube())
{
return 1.0/(1 << (dim-1));
}
else if (geomType.isTriangle())
{
static const Scalar faceToArea[] = {0.372677996249965,
0.372677996249965,
0.235702260395516};
return faceToArea[fIdx];
}
else if (geomType.isTetrahedron())
{
static const Scalar faceToArea[] = {0.102062072615966,
0.102062072615966,
0.058925565098879,
0.102062072615966,
0.058925565098879,
0.058925565098879};
return faceToArea[fIdx];
}
else if (geomType.isPyramid())
{
static const Scalar faceToArea[] = {0.130437298687488,
0.130437298687488,
0.130437298687488,
0.130437298687488,
0.150923085635624,
0.1092906420717,
0.1092906420717,
0.0781735959970571};
return faceToArea[fIdx];
}
else if (geomType.isPrism())
{
static const Scalar faceToArea[] = {0.166666666666667,
0.166666666666667,
0.166666666666667,
0.186338998124982,
0.186338998124982,
0.117851130197758,
0.186338998124982,
0.186338998124982,
0.117851130197758};
return faceToArea[fIdx];
}
else {
DUNE_THROW(Dune::NotImplemented, "scvf area for unknown GeometryType");
}
}
protected:
const Problem& problem_;
bool velocityOutput_;
......
test/references/1p2cniccconvection-reference.vtu
+ 8
− 8
View file @ 418b9c4f
......@@ -15,20 +15,20 @@
<DataArray type="Float32" Name="delp" NumberOfComponents="1" format="ascii">
21447.2 21183 20918.7 20654.5 20390.2 20126 19861.7 19597.4 19333.1 19068.6 18804 18539.1
18274 18008.5 17742.6 17476.1 17209.1 16941.5 16673.3 16404.6 16135.5 15865.9 15596 15325.7
15055.2 14784.5 14513.6 14242.6 13971.6 13700.4 13429.2 13158 12886.7 12615.5 12344.2 12072.9
15055.2 14784.5 14513.6 14242.6 13971.5 13700.4 13429.2 13158 12886.7 12615.5 12344.2 12072.9
11801.6 11530.3 11259 10987.7 10716.4 10445.1 10173.8 9902.53 9631.23 9359.92 9088.62 8817.32
8546.02 8274.72 8003.42 7732.12 7460.82 7189.52 6918.21 6646.91 6375.61 6104.31 5833.01 5561.71
5290.4 5019.1 4747.8 4476.5 4205.2 3933.89 3662.59 3391.29 3119.99 2848.68 2577.38 2306.08
2034.77 1763.47 1492.17 1220.87 949.562 678.259 406.955 135.652
</DataArray>
<DataArray type="Float32" Name="velocity" NumberOfComponents="3" format="ascii">
5e-05 0 0 0.0001 0 0 0.0001 0 0 0.0001 0 0
0.0001 0 0 0.0001 0 0 0.0001 0 0 0.0001 0 0
0.0001 0 0 0.0001 0 0 0.0001 0 0 0.0001 0 0
9.99999e-05 0 0 9.99997e-05 0 0 9.99994e-05 0 0 9.99991e-05 0 0
9.99987e-05 0 0 9.99982e-05 0 0 9.99977e-05 0 0 9.99972e-05 0 0
9.99967e-05 0 0 9.99961e-05 0 0 9.99956e-05 0 0 9.99952e-05 0 0
9.99948e-05 0 0 9.99944e-05 0 0 9.99941e-05 0 0 9.99938e-05 0 0
9.99936e-05 0 0 9.99934e-05 0 0 9.99933e-05 0 0 9.99932e-05 0 0
9.99936e-05 0 0 9.99935e-05 0 0 9.99933e-05 0 0 9.99932e-05 0 0
9.99932e-05 0 0 9.99931e-05 0 0 9.99931e-05 0 0 9.99931e-05 0 0
9.99931e-05 0 0 9.99931e-05 0 0 9.99931e-05 0 0 9.99931e-05 0 0
9.99931e-05 0 0 9.99931e-05 0 0 9.99931e-05 0 0 9.99931e-05 0 0
......@@ -89,9 +89,9 @@
998.79 998.79 998.79 998.79 998.79 998.789 998.789 998.789
</DataArray>
<DataArray type="Float32" Name="mu" NumberOfComponents="1" format="ascii">
0.00105696 0.00105696 0.00105697 0.00105697 0.001057 0.00105706 0.0010572 0.00105744 0.00105785 0.00105847 0.00105933 0.00106047
0.00106196 0.00106384 0.00106594 0.00106817 0.00107044 0.00107266 0.00107476 0.00107668 0.00107839 0.00107986 0.00108109 0.0010821
0.00108291 0.00108354 0.00108402 0.00108438 0.00108465 0.00108484 0.00108498 0.00108508 0.00108514 0.00108519 0.00108522 0.00108524
0.00105696 0.00105696 0.00105697 0.00105697 0.001057 0.00105706 0.0010572 0.00105744 0.00105785 0.00105847 0.00105934 0.00106047
0.00106196 0.00106385 0.00106595 0.00106817 0.00107044 0.00107266 0.00107476 0.00107668 0.00107839 0.00107985 0.00108109 0.0010821
0.0010829 0.00108354 0.00108402 0.00108438 0.00108465 0.00108484 0.00108498 0.00108508 0.00108514 0.00108519 0.00108522 0.00108524
0.00108525 0.00108526 0.00108527 0.00108527 0.00108527 0.00108527 0.00108528 0.00108528 0.00108528 0.00108528 0.00108528 0.00108528
0.00108528 0.00108528 0.00108528 0.00108528 0.00108528 0.00108528 0.00108528 0.00108528 0.00108528 0.00108528 0.00108528 0.00108528
0.00108528 0.00108528 0.00108528 0.00108528 0.00108528 0.00108528 0.00108528 0.00108528 0.00108528 0.00108528 0.00108528 0.00108528
......@@ -107,8 +107,8 @@
0 0 0 0 0 0 0 0
</DataArray>
<DataArray type="Float32" Name="temperature" NumberOfComponents="1" format="ascii">
291 291 291 291 290.999 290.996 290.991 290.981 290.965 290.94 290.906 290.861
290.805 290.74 290.667 290.591 290.512 290.436 290.363 290.297 290.238 290.188 290.145 290.111
291 291 291 291 290.999 290.996 290.991 290.981 290.965 290.94 290.906 290.86
290.805 290.74 290.667 290.59 290.512 290.436 290.363 290.297 290.239 290.188 290.145 290.111
290.083 290.061 290.044 290.032 290.023 290.016 290.011 290.008 290.006 290.004 290.003 290.002
290.002 290.002 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001
290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001
......
test/references/1pniccconvection-reference.vtu
+ 2
− 2
View file @ 418b9c4f
......@@ -22,7 +22,7 @@
1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
</DataArray>
<DataArray type="Float32" Name="velocity" NumberOfComponents="3" format="ascii">
5e-05 0 0 0.0001 0 0 0.0001 0 0 0.0001 0 0
0.0001 0 0 0.0001 0 0 0.0001 0 0 0.0001 0 0
0.0001 0 0 0.0001 0 0 0.0001 0 0 0.0001 0 0
0.0001 0 0 0.0001 0 0 0.0001 0 0 0.0001 0 0
9.99996e-05 0 0 9.9999e-05 0 0 9.99982e-05 0 0 9.99973e-05 0 0
......@@ -54,7 +54,7 @@
</DataArray>
<DataArray type="Float32" Name="temperature" NumberOfComponents="1" format="ascii">
291 291 291 291.001 291.002 291 290.997 290.998 291.005 291.014 291.015 290.994
290.944 290.862 290.754 290.632 290.508 290.391 290.29 290.208 290.143 290.096 290.063 290.04
290.943 290.861 290.754 290.632 290.508 290.391 290.29 290.208 290.144 290.096 290.063 290.04
290.025 290.015 290.009 290.006 290.004 290.003 290.002 290.002 290.001 290.001 290.001 290.001
290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001
290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001
......
test/references/3pniconvectioncc-reference.vtu
+ 1
− 1
View file @ 418b9c4f
......@@ -112,7 +112,7 @@
290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001
</DataArray>
<DataArray type="Float32" Name="velocityW" NumberOfComponents="3" format="ascii">
4.3676e-05 0 0 8.73521e-05 0 0 8.73521e-05 0 0 8.73521e-05 0 0
8.73521e-05 0 0 8.73521e-05 0 0 8.73521e-05 0 0 8.73521e-05 0 0
8.73521e-05 0 0 8.73521e-05 0 0 8.7352e-05 0 0 8.73519e-05 0 0
8.73517e-05 0 0 8.73514e-05 0 0 8.7351e-05 0 0 8.73505e-05 0 0
8.73501e-05 0 0 8.73496e-05 0 0 8.73491e-05 0 0 8.73486e-05 0 0
......
test/references/richardsniconvectioncc-reference.vtu
+ 6
− 6
View file @ 418b9c4f
......@@ -23,7 +23,7 @@
</DataArray>
<DataArray type="Float32" Name="pn" NumberOfComponents="1" format="ascii">
121447 121183 120919 120654 120390 120126 119862 119597 119333 119069 118804 118539
118274 118009 117743 117476 117209 116941 116673 116405 116136 115866 115596 115326
118274 118009 117743 117476 117209 116942 116673 116405 116136 115866 115596 115326
115055 114784 114514 114243 113972 113700 113429 113158 112887 112615 112344 112073
111802 111530 111259 110988 110716 110445 110174 109903 109631 109360 109089 108817
108546 108275 108003 107732 107461 107190 106918 106647 106376 106104 105833 105562
......@@ -32,7 +32,7 @@
</DataArray>
<DataArray type="Float32" Name="pw" NumberOfComponents="1" format="ascii">
121447 121183 120919 120654 120390 120126 119862 119597 119333 119069 118804 118539
118274 118009 117743 117476 117209 116941 116673 116405 116136 115866 115596 115326
118274 118009 117743 117476 117209 116942 116673 116405 116136 115866 115596 115326
115055 114784 114514 114243 113972 113700 113429 113158 112887 112615 112344 112073
111802 111530 111259 110988 110716 110445 110174 109903 109631 109360 109089 108817
108546 108275 108003 107732 107461 107190 106918 106647 106376 106104 105833 105562
......@@ -67,9 +67,9 @@
1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
</DataArray>
<DataArray type="Float32" Name="mobW" NumberOfComponents="1" format="ascii">
946.105 946.106 946.105 946.098 946.075 946.018 945.9 945.683 945.321 944.771 943.999 942.988
941.668 939.993 938.14 936.18 934.194 932.256 930.433 928.771 927.303 926.041 924.985 924.123
923.435 922.899 922.488 922.181 921.954 921.79 921.674 921.592 921.536 921.498 921.472 921.455
946.105 946.106 946.105 946.098 946.075 946.018 945.9 945.683 945.321 944.771 944 942.989
941.669 939.994 938.14 936.18 934.194 932.256 930.432 928.771 927.302 926.041 924.985 924.123
923.435 922.898 922.488 922.18 921.954 921.79 921.674 921.592 921.536 921.498 921.472 921.455
921.444 921.437 921.432 921.429 921.427 921.426 921.425 921.425 921.424 921.424 921.424 921.424
921.423 921.423 921.423 921.423 921.423 921.423 921.423 921.422 921.422 921.422 921.422 921.422
921.422 921.422 921.421 921.421 921.421 921.421 921.421 921.421 921.42 921.42 921.42 921.42
......@@ -103,7 +103,7 @@
290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001
</DataArray>
<DataArray type="Float32" Name="velocity" NumberOfComponents="3" format="ascii">
5e-05 0 0 0.0001 0 0 0.0001 0 0 0.0001 0 0
0.0001 0 0 0.0001 0 0 0.0001 0 0 0.0001 0 0
0.0001 0 0 0.0001 0 0 0.0001 0 0 0.0001 0 0
9.99999e-05 0 0 9.99997e-05 0 0 9.99995e-05 0 0 9.99991e-05 0 0
9.99987e-05 0 0 9.99982e-05 0 0 9.99977e-05 0 0 9.99972e-05 0 0
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment