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

implicit 2pdfm: The TwoPDFMModel has always done the same as the TwoPModel.... 

implicit 2pdfm: The TwoPDFMModel has always done the same as the TwoPModel. The obsolete reimplementation of addOutputVtkFields has been removed. Reviewed by Christoph.

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@10720 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent bde044f1
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dumux/implicit/2pdfm/2pdfmmodel.hh
+ 1
− 292
View file @ 09571e74
...@@ -69,299 +69,8 @@ namespace Dumux ...@@ -69,299 +69,8 @@ namespace Dumux
*/ */
template<class TypeTag > template<class TypeTag >
class TwoPDFMModel : public TwoPModel<TypeTag> class TwoPDFMModel : public TwoPModel<TypeTag>
{ {};
typedef TwoPModel<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
enum {
nPhaseIdx = Indices::nPhaseIdx,
wPhaseIdx = Indices::wPhaseIdx,
pressureIdx = Indices::pressureIdx,
numPhases = GET_PROP_VALUE(TypeTag, NumPhases)
};
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::ctype CoordScalar;
enum {
dim = GridView::dimension,
dimWorld = GridView::dimensionworld
};
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef Dune::FieldVector<Scalar, numPhases> PhasesVector;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
public:
/*!
* \brief Append all quantities of interest which can be derived
* from the solution of the current time step to the VTK
* writer.
*
* \param sol The global solution vector
* \param writer The writer for multi-file VTK datasets
*/
template<class MultiWriter>
void addOutputVtkFields(const SolutionVector &sol, MultiWriter &writer)
{
bool velocityOutput = GET_PARAM_FROM_GROUP(TypeTag, bool, Vtk, AddVelocity);
typedef Dune::BlockVector<Dune::FieldVector<double, 1> > ScalarField;
typedef Dune::BlockVector<Dune::FieldVector<double, dim> > VectorField;
// get the number of degrees of freedom and the number of vertices
unsigned numDofs = this->numDofs();
unsigned numVertices = this->problem_().gridView().size(dim);
// velocity output currently only works for vertex data
if (numDofs != numVertices)
{
velocityOutput = false;
}
// create the required scalar fields
ScalarField *pW = writer.allocateManagedBuffer(numDofs);
ScalarField *pN = writer.allocateManagedBuffer(numDofs);
ScalarField *pC = writer.allocateManagedBuffer(numDofs);
ScalarField *Sw = writer.allocateManagedBuffer(numDofs);
ScalarField *Sn = writer.allocateManagedBuffer(numDofs);
ScalarField *rhoW = writer.allocateManagedBuffer(numDofs);
ScalarField *rhoN = writer.allocateManagedBuffer(numDofs);
ScalarField *mobW = writer.allocateManagedBuffer(numDofs);
ScalarField *mobN = writer.allocateManagedBuffer(numDofs);
ScalarField *poro = writer.allocateManagedBuffer(numDofs);
ScalarField *Te = writer.allocateManagedBuffer(numDofs);
ScalarField *cellNum =writer.allocateManagedBuffer (numDofs);
VectorField *velocityN = writer.template allocateManagedBuffer<double, dim>(numDofs);
VectorField *velocityW = writer.template allocateManagedBuffer<double, dim>(numDofs);
if(velocityOutput) // check if velocity output is demanded
{
// initialize velocity fields
for (unsigned int i = 0; i < numVertices; ++i)
{
(*velocityN)[i] = Scalar(0);
(*velocityW)[i] = Scalar(0);
(*cellNum)[i] = Scalar(0.0);
}
}
unsigned numElements = this->gridView_().size(0);
ScalarField *rank = writer.allocateManagedBuffer(numElements);
FVElementGeometry fvGeometry;
VolumeVariables volVars;
ElementVolumeVariables elemVolVars;
ElementIterator elemIt = this->gridView_().template begin<0>();
ElementIterator elemEndIt = this->gridView_().template end<0>();
for (; elemIt != elemEndIt; ++elemIt)
{
if(velocityOutput && !elemIt->geometry().type().isCube()){
DUNE_THROW(Dune::InvalidStateException,
"Currently, velocity output only works for cubes. "
"Please set the EnableVelocityOutput property to false!");
}
int idx = this->elementMapper().map(*elemIt);
(*rank)[idx] = this->gridView_().comm().rank();
fvGeometry.update(this->gridView_(), *elemIt);
if (numDofs == numElements) // element data
{
volVars.update(sol[idx],
this->problem_(),
*elemIt,
fvGeometry,
/*scvIdx=*/0,
false);
(*pW)[idx] = volVars.pressure(wPhaseIdx);
(*pN)[idx] = volVars.pressure(nPhaseIdx);
(*pC)[idx] = volVars.capillaryPressure();
(*Sw)[idx] = volVars.saturation(wPhaseIdx);
(*Sn)[idx] = volVars.saturation(nPhaseIdx);
(*rhoW)[idx] = volVars.density(wPhaseIdx);
(*rhoN)[idx] = volVars.density(nPhaseIdx);
(*mobW)[idx] = volVars.mobility(wPhaseIdx);
(*mobN)[idx] = volVars.mobility(nPhaseIdx);
(*poro)[idx] = volVars.porosity();
(*Te)[idx] = volVars.temperature();
}
else // vertex data
{
int numVerts = elemIt->template count<dim> ();
for (int scvIdx = 0; scvIdx < numVerts; ++scvIdx)
{
int globalIdx = this->vertexMapper().map(*elemIt, scvIdx, dim);
volVars.update(sol[globalIdx],
this->problem_(),
*elemIt,
fvGeometry,
scvIdx,
false);
(*pW)[globalIdx] = volVars.pressure(wPhaseIdx);
(*pN)[globalIdx] = volVars.pressure(nPhaseIdx);
(*pC)[globalIdx] = volVars.capillaryPressure();
(*Sw)[globalIdx] = volVars.saturation(wPhaseIdx);
(*Sn)[globalIdx] = volVars.saturation(nPhaseIdx);
(*rhoW)[globalIdx] = volVars.density(wPhaseIdx);
(*rhoN)[globalIdx] = volVars.density(nPhaseIdx);
(*mobW)[globalIdx] = volVars.mobility(wPhaseIdx);
(*mobN)[globalIdx] = volVars.mobility(nPhaseIdx);
(*poro)[globalIdx] = volVars.porosity();
(*Te)[globalIdx] = volVars.temperature();
if(velocityOutput)
{
(*cellNum)[globalIdx] += 1;
}
}
if(velocityOutput)
{
// calculate vertex velocities
GlobalPosition tmpVelocity[numPhases];
for(int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
tmpVelocity[phaseIdx] = Scalar(0.0);
}
typedef Dune::BlockVector<Dune::FieldVector<Scalar, dim> > SCVVelocities;
SCVVelocities scvVelocityW(8), scvVelocityN(8);
scvVelocityW = 0;
scvVelocityN = 0;
ElementVolumeVariables elemVolVars;
elemVolVars.update(this->problem_(),
*elemIt,
fvGeometry,
false /* oldSol? */);
for (int faceIdx = 0; faceIdx < fvGeometry.numEdges; faceIdx++)
{
FluxVariables fluxVars(this->problem_(),
*elemIt,
fvGeometry,
faceIdx,
elemVolVars);
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
// local position of integration point
const Dune::FieldVector<Scalar, dim>& localPosIP = fvGeometry.subContVolFace[faceIdx].ipLocal;
// Transformation of the global normal vector to normal vector in the reference element
const typename Element::Geometry::JacobianTransposed& jacobianT1 =
elemIt->geometry().jacobianTransposed(localPosIP);
const GlobalPosition globalNormal = fluxVars.face().normal;
GlobalPosition localNormal;
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 subcontrolvolumeface
// in the reference element.
PhasesVector q;
q[phaseIdx] = fluxVars.volumeFlux(phaseIdx) / localArea;
// transform the normal Darcy velocity into a vector
tmpVelocity[phaseIdx] = localNormal;
tmpVelocity[phaseIdx] *= q[phaseIdx];
if(phaseIdx == wPhaseIdx)
{
scvVelocityW[fluxVars.face().i] += tmpVelocity[phaseIdx];
scvVelocityW[fluxVars.face().j] += tmpVelocity[phaseIdx];
}
else if(phaseIdx == nPhaseIdx)
{
scvVelocityN[fluxVars.face().i] += tmpVelocity[phaseIdx];
scvVelocityN[fluxVars.face().j] += tmpVelocity[phaseIdx];
}
}
}
typedef Dune::GenericReferenceElements<Scalar, dim> ReferenceElements;
const Dune::FieldVector<Scalar, dim> &localPos
= ReferenceElements::general(elemIt->geometry().type()).position(0, 0);
// get the transposed Jacobian of the element mapping
const typename Element::Geometry::JacobianTransposed& jacobianT2
= elemIt->geometry().jacobianTransposed(localPos);
// transform vertex velocities from local to global coordinates
for (int i = 0; i < numVerts; ++i)
{
int globalIdx = this->vertexMapper().map(*elemIt, i, dim);
// calculate the subcontrolvolume velocity by the Piola transformation
Dune::FieldVector<CoordScalar, dim> scvVelocity(0);
jacobianT2.mtv(scvVelocityW[i], scvVelocity);
scvVelocity /= elemIt->geometry().integrationElement(localPos);
// add up the wetting phase subcontrolvolume velocities for each vertex
(*velocityW)[globalIdx] += scvVelocity;
jacobianT2.mtv(scvVelocityN[i], scvVelocity);
scvVelocity /= elemIt->geometry().integrationElement(localPos);
// add up the nonwetting phase subcontrolvolume velocities for each vertex
(*velocityN)[globalIdx] += scvVelocity;
}
}
}
}
if (numDofs == numElements) // element data
{
writer.attachCellData(*Sn, "Sn");
writer.attachCellData(*Sw, "Sw");
writer.attachCellData(*pN, "pn");
writer.attachCellData(*pW, "pw");
writer.attachCellData(*pC, "pc");
writer.attachCellData(*rhoW, "rhoW");
writer.attachCellData(*rhoN, "rhoN");
writer.attachCellData(*mobW, "mobW");
writer.attachCellData(*mobN, "mobN");
writer.attachCellData(*poro, "porosity");
writer.attachCellData(*Te, "temperature");
}
else // vertex data
{
writer.attachVertexData(*Sn, "Sn");
writer.attachVertexData(*Sw, "Sw");
writer.attachVertexData(*pN, "pn");
writer.attachVertexData(*pW, "pw");
writer.attachVertexData(*pC, "pc");
writer.attachVertexData(*rhoW, "rhoW");
writer.attachVertexData(*rhoN, "rhoN");
writer.attachVertexData(*mobW, "mobW");
writer.attachVertexData(*mobN, "mobN");
writer.attachVertexData(*poro, "porosity");
writer.attachVertexData(*Te, "temperature");
if(velocityOutput) // check if velocity output is demanded
{
// divide the vertex velocities by the number of adjacent scvs i.e. cells
for(unsigned int globalIdx = 0; globalIdx < numVertices; ++globalIdx)
{
(*velocityW)[globalIdx] /= (*cellNum)[globalIdx];
(*velocityN)[globalIdx] /= (*cellNum)[globalIdx];
}
writer.attachVertexData(*velocityW, "velocityW", dim);
writer.attachVertexData(*velocityN, "velocityN", dim);
}
}
writer.attachCellData(*rank, "process rank");
}
};
} // end namespace } // end namespace
#include "2pdfmpropertydefaults.hh" #include "2pdfmpropertydefaults.hh"
......
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