From 09571e74590b066e6861c385677e2a11d9b1bb46 Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Tue, 21 May 2013 13:51:05 +0000
Subject: [PATCH] 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
---
dumux/implicit/2pdfm/2pdfmmodel.hh | 293 +----------------------------
1 file changed, 1 insertion(+), 292 deletions(-)
diff --git a/dumux/implicit/2pdfm/2pdfmmodel.hh b/dumux/implicit/2pdfm/2pdfmmodel.hh
index e6ab7565d5..208067e147 100644
--- a/dumux/implicit/2pdfm/2pdfmmodel.hh
+++ b/dumux/implicit/2pdfm/2pdfmmodel.hh
@@ -69,299 +69,8 @@ namespace Dumux
*/
template<class 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
#include "2pdfmpropertydefaults.hh"
--
GitLab