diff --git a/dumux/implicit/1p/1pmodel.hh b/dumux/implicit/1p/1pmodel.hh
index 40afc0384823a7d9f9d6f9d678f8e83df9fc301c..40b419d33f8f7507715b649a72fe68ef61965928 100644
--- a/dumux/implicit/1p/1pmodel.hh
+++ b/dumux/implicit/1p/1pmodel.hh
@@ -59,6 +59,8 @@ class OnePBoxModel : public GET_PROP_TYPE(TypeTag, BaseModel)
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
enum { dim = GridView::dimension };
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+
public:
/*!
* \brief \copybrief Dumux::BoxModel::addOutputVtkFields
@@ -73,9 +75,9 @@ public:
typedef Dune::BlockVector<Dune::FieldVector<double, 1> > ScalarField;
// create the required scalar fields
- unsigned numVertices = this->problem_().gridView().size(dim);
- ScalarField *p = writer.allocateManagedBuffer(numVertices);
- ScalarField *K = writer.allocateManagedBuffer(numVertices);
+ unsigned numDofs = this->numDofs();
+ ScalarField *p = writer.allocateManagedBuffer(numDofs);
+ ScalarField *K = writer.allocateManagedBuffer(numDofs);
unsigned numElements = this->gridView_().size(0);
ScalarField *rank = writer.allocateManagedBuffer(numElements);
@@ -94,27 +96,40 @@ public:
fvGeometry.update(this->gridView_(), *elemIt);
elemBcTypes.update(this->problem_(), *elemIt, fvGeometry);
- int numVerts = elemIt->template count<dim> ();
- for (int i = 0; i < numVerts; ++i)
+ for (int scvIdx = 0; scvIdx < fvGeometry.numSCV; ++scvIdx)
{
- int globalIdx = this->vertexMapper().map(*elemIt, i, dim);
+ int globalIdx;
+ if (isBox)
+ globalIdx = this->vertexMapper().map(*elemIt, scvIdx, dim);
+ else
+ globalIdx = this->elementMapper().map(*elemIt);
+
volVars.update(sol[globalIdx],
this->problem_(),
*elemIt,
fvGeometry,
- i,
+ scvIdx,
false);
const SpatialParams &spatialParams = this->problem_().spatialParams();
(*p)[globalIdx] = volVars.pressure();
(*K)[globalIdx]= spatialParams.intrinsicPermeability(*elemIt,
fvGeometry,
- i);
+ scvIdx);
}
}
- writer.attachVertexData(*p, "p");
- writer.attachVertexData(*K, "K");
+ if (isBox)
+ {
+ writer.attachVertexData(*p, "p");
+ writer.attachVertexData(*K, "K");
+ }
+ else
+ {
+ writer.attachCellData(*p, "p");
+ writer.attachCellData(*K, "K");
+ }
+
writer.attachCellData(*rank, "process rank");
}
};
diff --git a/dumux/implicit/1p2c/1p2cfluxvariables.hh b/dumux/implicit/1p2c/1p2cfluxvariables.hh
index eca03684d0d2abb5d61c84759bb105e23b8ad148..ce408eee06c6449c307c009bec85277becd0beea 100644
--- a/dumux/implicit/1p2c/1p2cfluxvariables.hh
+++ b/dumux/implicit/1p2c/1p2cfluxvariables.hh
@@ -455,7 +455,7 @@ protected:
dispersionTensor_ = 0;
for (int i=0; i<dim; i++)
for (int j = 0; j<dim; j++)
- dispersionTensor_[i][j]=velocity[i]*velocity[j];
+ dispersionTensor_[i][j] = velocity[i]*velocity[j];
//normalize velocity product --> vv^T/||v||, [m/s]
Scalar vNorm = velocity.two_norm();
diff --git a/dumux/implicit/1p2c/1p2clocalresidual.hh b/dumux/implicit/1p2c/1p2clocalresidual.hh
index 9dbcde20fec33fcf660c9209ff9221a438ab00cb..815dd0a15678d6c2a7f1c0f1c045c34d4d67319b 100644
--- a/dumux/implicit/1p2c/1p2clocalresidual.hh
+++ b/dumux/implicit/1p2c/1p2clocalresidual.hh
@@ -25,7 +25,6 @@
#ifndef DUMUX_ONEP_TWOC_LOCAL_RESIDUAL_HH
#define DUMUX_ONEP_TWOC_LOCAL_RESIDUAL_HH
-#define VELOCITY_OUTPUT 1 //1 turns velocity output on, 0 turns it off
#include <dumux/boxmodels/common/boxmodel.hh>
diff --git a/dumux/implicit/1p2c/1p2cmodel.hh b/dumux/implicit/1p2c/1p2cmodel.hh
index 3b70302f8581ff7ad664753fdd8c5b17f78c403e..d565a1fa1a9fd1d29125bbb52be510b02bf11591 100644
--- a/dumux/implicit/1p2c/1p2cmodel.hh
+++ b/dumux/implicit/1p2c/1p2cmodel.hh
@@ -90,6 +90,8 @@ class OnePTwoCBoxModel : public GET_PROP_TYPE(TypeTag, BaseModel)
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef Dune::FieldVector<Scalar, dim> DimVector;
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+
public:
/*!
* \brief Constructor. Sets the upwind weight.
@@ -112,40 +114,47 @@ public:
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;
// create the required scalar fields
- unsigned numVertices = this->problem_().gridView().size(dim);
- ScalarField &pressure = *writer.allocateManagedBuffer(numVertices);
- ScalarField &delp = *writer.allocateManagedBuffer(numVertices);
- ScalarField &moleFraction0 = *writer.allocateManagedBuffer(numVertices);
- ScalarField &moleFraction1 = *writer.allocateManagedBuffer(numVertices);
- ScalarField &massFraction0 = *writer.allocateManagedBuffer(numVertices);
- ScalarField &massFraction1 = *writer.allocateManagedBuffer(numVertices);
- ScalarField &rho = *writer.allocateManagedBuffer(numVertices);
- ScalarField &mu = *writer.allocateManagedBuffer(numVertices);
-#ifdef VELOCITY_OUTPUT // check if velocity output is demanded
- ScalarField &velocityX = *writer.allocateManagedBuffer(numVertices);
- ScalarField &velocityY = *writer.allocateManagedBuffer(numVertices);
- ScalarField &velocityZ = *writer.allocateManagedBuffer(numVertices);
+ unsigned numDofs = this->numDofs();
+ ScalarField &pressure = *writer.allocateManagedBuffer(numDofs);
+ ScalarField &delp = *writer.allocateManagedBuffer(numDofs);
+ ScalarField &moleFraction0 = *writer.allocateManagedBuffer(numDofs);
+ ScalarField &moleFraction1 = *writer.allocateManagedBuffer(numDofs);
+ ScalarField &massFraction0 = *writer.allocateManagedBuffer(numDofs);
+ ScalarField &massFraction1 = *writer.allocateManagedBuffer(numDofs);
+ ScalarField &rho = *writer.allocateManagedBuffer(numDofs);
+ ScalarField &mu = *writer.allocateManagedBuffer(numDofs);
+ ScalarField &velocityX = *writer.allocateManagedBuffer(numDofs);
+ ScalarField &velocityY = *writer.allocateManagedBuffer(numDofs);
+ ScalarField &velocityZ = *writer.allocateManagedBuffer(numDofs);
//use vertical faces for vx and horizontal faces for vy calculation
- std::vector<DimVector> boxSurface(numVertices);
- // initialize velocity fields
- for (unsigned int i = 0; i < numVertices; ++i)
- {
-
- velocityX[i] = 0;
- if (dim > 1)
- {
- velocityY[i] = 0;
- }
- if (dim > 2)
- {
- velocityZ[i] = 0;
- }
- boxSurface[i] = Scalar(0.0); // initialize the boundary surface of the fv-boxes
- }
-#endif
+ std::vector<DimVector> boxSurface(numDofs);
+
+ // velocity output currently only works for the box discretization
+ if (!isBox)
+ velocityOutput = false;
+
+ if (velocityOutput)
+ {
+ // initialize velocity fields
+ for (unsigned int i = 0; i < numDofs; ++i)
+ {
+ velocityX[i] = 0;
+ if (dim > 1)
+ {
+ velocityY[i] = 0;
+ }
+ if (dim > 2)
+ {
+ velocityZ[i] = 0;
+ }
+ boxSurface[i] = Scalar(0.0); // initialize the boundary surface of the fv-boxes
+ }
+ }
+
unsigned numElements = this->gridView_().size(0);
ScalarField &rank =
*writer.allocateManagedBuffer(numElements);
@@ -164,15 +173,19 @@ public:
fvGeometry.update(this->gridView_(), *elemIt);
elemBcTypes.update(this->problem_(), *elemIt, fvGeometry);
- int numVerts = elemIt->template count<dim> ();
- for (int i = 0; i < numVerts; ++i)
+ for (int scvIdx = 0; scvIdx < fvGeometry.numSCV; ++scvIdx)
{
- int globalIdx = this->vertexMapper().map(*elemIt, i, dim);
+ int globalIdx;
+ if (isBox)
+ globalIdx = this->vertexMapper().map(*elemIt, scvIdx, dim);
+ else
+ globalIdx = this->elementMapper().map(*elemIt);
+
volVars.update(sol[globalIdx],
this->problem_(),
*elemIt,
fvGeometry,
- i,
+ scvIdx,
false);
pressure[globalIdx] = volVars.pressure();
@@ -185,134 +198,161 @@ public:
mu[globalIdx] = volVars.viscosity();
}
-#ifdef VELOCITY_OUTPUT // check if velocity output is demanded
- // In the box method, the velocity is evaluated on the FE-Grid. However, to get an
- // average apparent velocity at the vertex, all contributing velocities have to be interpolated.
- DimVector velocity;
-
- ElementVolumeVariables elemVolVars;
- elemVolVars.update(this->problem_(),
- *elemIt,
- fvGeometry,
- false /* isOldSol? */);
- // loop over the phases
- for (int faceIdx = 0; faceIdx < fvGeometry.numEdges; faceIdx++)
+ if (velocityOutput)
{
- velocity = 0.0;
- //prepare the flux calculations (set up and prepare geometry, FE gradients)
- FluxVariables fluxVars(this->problem_(),
- *elemIt,
- fvGeometry,
- faceIdx,
- elemVolVars);
+ // In the box method, the velocity is evaluated on the FE-Grid. However, to get an
+ // average apparent velocity at the vertex, all contributing velocities have to be interpolated.
+ DimVector velocity;
+
+ ElementVolumeVariables elemVolVars;
+ elemVolVars.update(this->problem_(),
+ *elemIt,
+ fvGeometry,
+ false /* isOldSol? */);
+ // loop over the phases
+ for (int faceIdx = 0; faceIdx < fvGeometry.numEdges; faceIdx++)
+ {
+ velocity = 0.0;
+ //prepare the flux calculations (set up and prepare geometry, FE gradients)
+ FluxVariables fluxVars(this->problem_(),
+ *elemIt,
+ fvGeometry,
+ faceIdx,
+ elemVolVars);
+
+ //use vertical faces for vx and horizontal faces for vy calculation
+ DimVector xVector(0), yVector(0);
+ xVector[0] = 1; yVector[1] = 1;
+
+ Dune::SeqScalarProduct<DimVector> sp;
+
+ Scalar xDir = std::abs(sp.dot(fluxVars.face().normal, xVector));
+ Scalar yDir = std::abs(sp.dot(fluxVars.face().normal, yVector));
+
+ // up+downstream node
+ const VolumeVariables &up =
+ elemVolVars[fluxVars.upstreamIdx()];
+ const VolumeVariables &dn =
+ elemVolVars[fluxVars.downstreamIdx()];
+
+ //get surface area to weight velocity at the IP with the surface area
+ Scalar scvfArea = fluxVars.face().normal.two_norm();
+
+ int vertIIdx = this->problem_().vertexMapper().map(
+ *elemIt, fluxVars.face().i, dim);
+ int vertJIdx = this->problem_().vertexMapper().map(
+ *elemIt, fluxVars.face().j, dim);
+
+ //use vertical faces (horizontal noraml vector) to calculate vx
+ //in case of heterogeneities it seams to be better to define intrinisc permeability elementwise
+ if(xDir > yDir)//(fluxVars.face().normal[0] > 1e-10 || fluxVars.face().normal[0] < -1e-10)// (xDir > yDir)
+ {
+ // get darcy velocity
+ //calculate (v n) n/A
+ Scalar tmp = fluxVars.KmvpNormal();
+ velocity = fluxVars.face().normal;
+ velocity *= tmp;
+ velocity /= scvfArea;
+ velocity *= (upwindWeight_ / up.viscosity() +
+ (1 - upwindWeight_)/ dn.viscosity());
+
+ // add surface area for weighting purposes
+ boxSurface[vertIIdx][0] += scvfArea;
+ boxSurface[vertJIdx][0] += scvfArea;
+
+ velocityX[vertJIdx] += velocity[0];
+ velocityX[vertIIdx] += velocity[0];
+
+ }
+ if (yDir > xDir)//(fluxVars.face().normal[1] > 1e-10 || fluxVars.face().normal[1] < -1e-10)// (yDir > xDir)
+ {
+ // get darcy velocity
+ //calculate (v n) n/A
+ Scalar tmp = fluxVars.KmvpNormal();
+ velocity = fluxVars.face().normal;
+ velocity *= tmp;
+ velocity /= scvfArea;
+ velocity *= (upwindWeight_ / up.viscosity() +
+ (1 - upwindWeight_)/ dn.viscosity());
+
+ // add surface area for weighting purposes
+ boxSurface[vertIIdx][1] += scvfArea;
+ boxSurface[vertJIdx][1] += scvfArea;
+
+ velocityY[vertJIdx] += velocity[1];
+ velocityY[vertIIdx] += velocity[1];
+ }
+ }
+ }
+ }
+
+ if (velocityOutput)
+ {
+ // normalize the velocities at the vertices
+ // calculate the bounding box of the grid view
+ VertexIterator vIt = this->gridView_().template begin<dim>();
+ const VertexIterator vEndIt = this->gridView_().template end<dim>();
+ for (; vIt!=vEndIt; ++vIt)
+ {
+ int i = this->problem_().vertexMapper().map(*vIt);
//use vertical faces for vx and horizontal faces for vy calculation
- DimVector xVector(0), yVector(0);
- xVector[0] = 1; yVector[1] = 1;
-
- Dune::SeqScalarProduct<DimVector> sp;
-
- Scalar xDir = std::abs(sp.dot(fluxVars.face().normal, xVector));
- Scalar yDir = std::abs(sp.dot(fluxVars.face().normal, yVector));
-
- // up+downstream node
- const VolumeVariables &up =
- elemVolVars[fluxVars.upstreamIdx()];
- const VolumeVariables &dn =
- elemVolVars[fluxVars.downstreamIdx()];
-
- //get surface area to weight velocity at the IP with the surface area
- Scalar scvfArea = fluxVars.face().normal.two_norm();
-
- int vertIIdx = this->problem_().vertexMapper().map(
- *elemIt, fluxVars.face().i, dim);
- int vertJIdx = this->problem_().vertexMapper().map(
- *elemIt, fluxVars.face().j, dim);
-
- //use vertical faces (horizontal noraml vector) to calculate vx
- //in case of heterogeneities it seams to be better to define intrinisc permeability elementwise
- if(xDir > yDir)//(fluxVars.face().normal[0] > 1e-10 || fluxVars.face().normal[0] < -1e-10)// (xDir > yDir)
+ velocityX[i] /= boxSurface[i][0];
+ if (dim >= 2)
{
- // get darcy velocity
- //calculate (v n) n/A
- Scalar tmp = fluxVars.KmvpNormal();
- velocity = fluxVars.face().normal;
- velocity *= tmp;
- velocity /= scvfArea;
- velocity *= (upwindWeight_ / up.viscosity() +
- (1 - upwindWeight_)/ dn.viscosity());
-
- // add surface area for weighting purposes
- boxSurface[vertIIdx][0] += scvfArea;
- boxSurface[vertJIdx][0] += scvfArea;
-
- velocityX[vertJIdx] += velocity[0];
- velocityX[vertIIdx] += velocity[0];
-
+ velocityY[i] /= boxSurface[i][1];
}
- if (yDir > xDir)//(fluxVars.face().normal[1] > 1e-10 || fluxVars.face().normal[1] < -1e-10)// (yDir > xDir)
+ if (dim == 3)
{
- // get darcy velocity
- //calculate (v n) n/A
- Scalar tmp = fluxVars.KmvpNormal();
- velocity = fluxVars.face().normal;
- velocity *= tmp;
- velocity /= scvfArea;
- velocity *= (upwindWeight_ / up.viscosity() +
- (1 - upwindWeight_)/ dn.viscosity());
-
- // add surface area for weighting purposes
- boxSurface[vertIIdx][1] += scvfArea;
- boxSurface[vertJIdx][1] += scvfArea;
-
- velocityY[vertJIdx] += velocity[1];
- velocityY[vertIIdx] += velocity[1];
+ velocityZ[i] /= boxSurface[i][2];
}
}
-#endif
}
-#ifdef VELOCITY_OUTPUT
- // normalize the velocities at the vertices
- // calculate the bounding box of the grid view
- VertexIterator vIt = this->gridView_().template begin<dim>();
- const VertexIterator vEndIt = this->gridView_().template end<dim>();
- for (; vIt!=vEndIt; ++vIt)
- {
- int i = this->problem_().vertexMapper().map(*vIt);
- //use vertical faces for vx and horizontal faces for vy calculation
- velocityX[i] /= boxSurface[i][0];
- if (dim >= 2)
- {
- velocityY[i] /= boxSurface[i][1];
- }
- if (dim == 3)
+ if (isBox)
+ {
+ writer.attachVertexData(pressure, "P");
+ writer.attachVertexData(delp, "delp");
+ if (velocityOutput)
{
- velocityZ[i] /= boxSurface[i][2];
+ writer.attachVertexData(velocityX, "Vx");
+ writer.attachVertexData(velocityY, "Vy");
+ if (dim > 2)
+ writer.attachVertexData(velocityZ, "Vz");
}
+ char nameMoleFraction0[42], nameMoleFraction1[42];
+ snprintf(nameMoleFraction0, 42, "x_%s", FluidSystem::componentName(0));
+ snprintf(nameMoleFraction1, 42, "x_%s", FluidSystem::componentName(1));
+ writer.attachVertexData(moleFraction0, nameMoleFraction0);
+ writer.attachVertexData(moleFraction1, nameMoleFraction1);
+
+ char nameMassFraction0[42], nameMassFraction1[42];
+ snprintf(nameMassFraction0, 42, "X_%s", FluidSystem::componentName(0));
+ snprintf(nameMassFraction1, 42, "X_%s", FluidSystem::componentName(1));
+ writer.attachVertexData(massFraction0, nameMassFraction0);
+ writer.attachVertexData(massFraction1, nameMassFraction1);
+ writer.attachVertexData(rho, "rho");
+ writer.attachVertexData(mu, "mu");
}
-#endif
- writer.attachVertexData(pressure, "P");
- writer.attachVertexData(delp, "delp");
-#ifdef VELOCITY_OUTPUT // check if velocity output is demanded
- writer.attachVertexData(velocityX, "Vx");
- writer.attachVertexData(velocityY, "Vy");
- if (dim > 2)
- writer.attachVertexData(velocityZ, "Vz");
-#endif
- char nameMoleFraction0[42], nameMoleFraction1[42];
- snprintf(nameMoleFraction0, 42, "x_%s", FluidSystem::componentName(0));
- snprintf(nameMoleFraction1, 42, "x_%s", FluidSystem::componentName(1));
- writer.attachVertexData(moleFraction0, nameMoleFraction0);
- writer.attachVertexData(moleFraction1, nameMoleFraction1);
-
- char nameMassFraction0[42], nameMassFraction1[42];
- snprintf(nameMassFraction0, 42, "X_%s", FluidSystem::componentName(0));
- snprintf(nameMassFraction1, 42, "X_%s", FluidSystem::componentName(1));
- writer.attachVertexData(massFraction0, nameMassFraction0);
- writer.attachVertexData(massFraction1, nameMassFraction1);
- writer.attachVertexData(rho, "rho");
- writer.attachVertexData(mu, "mu");
+ else
+ {
+ writer.attachCellData(pressure, "P");
+ writer.attachCellData(delp, "delp");
+ char nameMoleFraction0[42], nameMoleFraction1[42];
+ snprintf(nameMoleFraction0, 42, "x_%s", FluidSystem::componentName(0));
+ snprintf(nameMoleFraction1, 42, "x_%s", FluidSystem::componentName(1));
+ writer.attachCellData(moleFraction0, nameMoleFraction0);
+ writer.attachCellData(moleFraction1, nameMoleFraction1);
+
+ char nameMassFraction0[42], nameMassFraction1[42];
+ snprintf(nameMassFraction0, 42, "X_%s", FluidSystem::componentName(0));
+ snprintf(nameMassFraction1, 42, "X_%s", FluidSystem::componentName(1));
+ writer.attachCellData(massFraction0, nameMassFraction0);
+ writer.attachCellData(massFraction1, nameMassFraction1);
+ writer.attachCellData(rho, "rho");
+ writer.attachCellData(mu, "mu");
+ }
+
writer.attachCellData(rank, "process rank");
}
diff --git a/dumux/implicit/1p2c/1p2cproperties.hh b/dumux/implicit/1p2c/1p2cproperties.hh
index 37aa1eb857f4fac6c2e5c9a6ec038b9fe20bbcd5..856141eb584d7a6f2e61997c7348a5ad1e938b6b 100644
--- a/dumux/implicit/1p2c/1p2cproperties.hh
+++ b/dumux/implicit/1p2c/1p2cproperties.hh
@@ -63,6 +63,7 @@ NEW_PROP_TAG(ProblemEnableGravity); //!< Returns whether gravity is considered i
NEW_PROP_TAG(UseMoles); //!Defines whether mole (true) or mass (false) fractions are used
NEW_PROP_TAG(Scaling); //!Defines Scaling of the model
NEW_PROP_TAG(SpatialParamsForchCoeff); //!< Property for the forchheimer coefficient
+NEW_PROP_TAG(VtkAddVelocity); //!< Returns whether velocity vectors are written into the vtk output
}
// \}
}
diff --git a/dumux/implicit/1p2c/1p2cpropertydefaults.hh b/dumux/implicit/1p2c/1p2cpropertydefaults.hh
index 6e6340dcc551fee88f1afc6e5e4cf0bdf725225e..081a8b1d7dace3f143d23df836b01b3366a546fa 100644
--- a/dumux/implicit/1p2c/1p2cpropertydefaults.hh
+++ b/dumux/implicit/1p2c/1p2cpropertydefaults.hh
@@ -79,6 +79,9 @@ SET_TYPE_PROP(OnePTwoC, SpatialParams, ImplicitSpatialParamsOneP<TypeTag>);
//! Set the phaseIndex per default to zero (important for two-phase fluidsystems).
SET_INT_PROP(OnePTwoC, PhaseIdx, 0);
+// disable velocity output by default
+SET_BOOL_PROP(OnePTwoC, VtkAddVelocity, false);
+
// enable gravity by default
SET_BOOL_PROP(OnePTwoC, ProblemEnableGravity, true);
diff --git a/test/implicit/1p/1ptestproblem.hh b/test/implicit/1p/1ptestproblem.hh
index 1ec22776767b502fe97197cc5b79f637a5c4a2fe..799ab619a5e9cab957e8c1ac6adb8df35436ed25 100644
--- a/test/implicit/1p/1ptestproblem.hh
+++ b/test/implicit/1p/1ptestproblem.hh
@@ -49,7 +49,9 @@ class OnePTestProblem;
namespace Properties
{
-NEW_TYPE_TAG(OnePTestProblem, INHERITS_FROM(BoxOneP));
+NEW_TYPE_TAG(OnePTestProblem, INHERITS_FROM(OneP));
+NEW_TYPE_TAG(OnePTestBoxProblem, INHERITS_FROM(BoxModel, OnePTestProblem));
+NEW_TYPE_TAG(OnePTestCCProblem, INHERITS_FROM(CCModel, OnePTestProblem));
SET_PROP(OnePTestProblem, Fluid)
{
@@ -141,8 +143,12 @@ class OnePTestProblem : public ImplicitPorousMediaProblem<TypeTag>
public:
OnePTestProblem(TimeManager &timeManager, const GridView &gridView)
- : ParentType(timeManager, gridView)
+ : ParentType(timeManager, gridView)
{
+ name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag,
+ std::string,
+ Problem,
+ Name);
}
/*!
@@ -156,7 +162,9 @@ public:
* This is used as a prefix for files generated by the simulation.
*/
const char *name() const
- { return "1ptest"; }
+ {
+ return name_.c_str();
+ }
/*!
* \brief Return the temperature within the domain.
@@ -179,13 +187,15 @@ public:
// \{
/*!
- * \brief Specify which kind of boundary condition should be
- * used for which equation on a given boundary segment.
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary control volume.
+ *
+ * \param values The boundary types for the conservation equations
+ * \param globalPos The position of the center of the finite volume
*/
- void boundaryTypes(BoundaryTypes &values, const Vertex &vertex) const
+ void boundaryTypesAtPos(BoundaryTypes &values,
+ const GlobalPosition &globalPos) const
{
- const GlobalPosition globalPos = vertex.geometry().center();
-
double eps = 1.0e-3;
if (globalPos[dim-1] < eps || globalPos[dim-1] > this->bboxMax()[dim-1] - eps)
values.setAllDirichlet();
@@ -194,22 +204,18 @@ public:
}
/*!
- * \brief Evaluate the boundary conditions for a Dirichlet
- * boundary segment.
+ * \brief Evaluate the boundary conditions for a dirichlet
+ * control volume.
*
- * For this method, the \a priVars parameter stores primary variables.
+ * \param values The dirichlet values for the primary variables
+ * \param globalPos The center of the finite volume which ought to be set.
+ *
+ * For this method, the \a values parameter stores primary variables.
*/
- void dirichlet(PrimaryVariables &priVars, const Vertex &vertex) const
+ void dirichletAtPos(PrimaryVariables &values,
+ const GlobalPosition &globalPos) const
{
- double eps = 1.0e-3;
- const GlobalPosition globalPos = vertex.geometry().center();
-
- if (globalPos[dim-1] < eps) {
- priVars[pressureIdx] = 2.0e+5;
- }
- else if (globalPos[dim-1] > this->bboxMax()[dim-1] - eps) {
- priVars[pressureIdx] = 1.0e+5;
- }
+ values[pressureIdx] = 1.0e+5*(2.0 - globalPos[dim-1]);
}
/*!
@@ -228,8 +234,6 @@ public:
const int scvIdx,
const int boundaryFaceIdx) const
{
- // const GlobalPosition &globalPos = fvGeometry.boundaryFace[boundaryFaceIdx].ipGlobal;
-
priVars[conti0EqIdx] = 0;
}
@@ -251,11 +255,13 @@ public:
const FVElementGeometry &fvGeometry,
const int scvIdx) const
{
- //const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
- priVars[pressureIdx] = 1.0e+5;// + 9.81*1.23*(20-globalPos[dim-1]);
+ priVars[pressureIdx] = 1.0e+5;
}
// \}
+
+private:
+ std::string name_;
};
} //end namespace
diff --git a/test/implicit/1p/Makefile.am b/test/implicit/1p/Makefile.am
index 37b54a153ed5b865fa2e76bb48a7ae95335e9377..a50f7eb4f1bb50229333c3650e3db71268282622 100644
--- a/test/implicit/1p/Makefile.am
+++ b/test/implicit/1p/Makefile.am
@@ -1,10 +1,14 @@
-check_PROGRAMS = test_1p
+check_PROGRAMS = test_box1p test_cc1p
noinst_HEADERS = *.hh
EXTRA_DIST=*reference.vtu *.input grids/*.dgf CMakeLists.txt
-test_1p_SOURCES = test_1p.cc
-test_1p_CXXFLAGS = $(DUNEMPICPPFLAGS)
-test_1p_LDADD = $(DUNEMPILDFLAGS)
+test_box1p_SOURCES = test_box1p.cc
+test_box1p_CXXFLAGS = $(DUNEMPICPPFLAGS)
+test_box1p_LDADD = $(DUNEMPILDFLAGS)
+
+test_cc1p_SOURCES = test_cc1p.cc
+test_cc1p_CXXFLAGS = $(DUNEMPICPPFLAGS)
+test_cc1p_LDADD = $(DUNEMPILDFLAGS)
include $(top_srcdir)/am/global-rules
diff --git a/test/implicit/1p/test_box1p.cc b/test/implicit/1p/test_box1p.cc
new file mode 100644
index 0000000000000000000000000000000000000000..5436a789cabe623fb30c4b82ebd5595b9343e331
--- /dev/null
+++ b/test/implicit/1p/test_box1p.cc
@@ -0,0 +1,64 @@
+// -*- 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 test for the one-phase box model
+ */
+#include "config.h"
+#include "1ptestproblem.hh"
+#include <dumux/common/start.hh>
+
+/*!
+ * \brief Provides an interface for customizing error messages associated with
+ * reading in parameters.
+ *
+ * \param progName The name of the program, that was tried to be started.
+ * \param errorMsg The error message that was issued by the start function.
+ * Comprises the thing that went wrong and a general help message.
+ */
+void usage(const char *progName, const std::string &errorMsg)
+{
+ if (errorMsg.size() > 0) {
+ std::string errorMessageOut = "\nUsage: ";
+ errorMessageOut += progName;
+ errorMessageOut += " [options]\n";
+ errorMessageOut += errorMsg;
+ errorMessageOut += "\n\nThe list of mandatory arguments for this program is:\n"
+ "\t-TimeManager.TEnd End of the simulation [s] \n"
+ "\t-TimeManager.DtInitial Initial timestep size [s] \n"
+ "\t-Grid.File Name of the file containing the grid \n"
+ "\t definition in DGF format\n"
+ "\t-SpatialParams.LensLowerLeftX x-coordinate of the lower left corner of the lens [m] \n"
+ "\t-SpatialParams.LensLowerLeftY y-coordinate of the lower left corner of the lens [m] \n"
+ "\t-SpatialParams.LensUpperRightX x-coordinate of the upper right corner of the lens [m] \n"
+ "\t-SpatialParams.LensUpperRightY y-coordinate of the upper right corner of the lens [m] \n"
+ "\t-SpatialParams.Permeability Permeability of the domain [m^2] \n"
+ "\t-SpatialParams.PermeabilityLens Permeability of the lens [m^2] \n";
+
+ std::cout << errorMessageOut
+ << "\n";
+ }
+}
+
+int main(int argc, char** argv)
+{
+ typedef TTAG(OnePTestBoxProblem) ProblemTypeTag;
+ return Dumux::start<ProblemTypeTag>(argc, argv, usage);
+}
diff --git a/test/implicit/1p/test_box1p.input b/test/implicit/1p/test_box1p.input
new file mode 100644
index 0000000000000000000000000000000000000000..2245b7effb1dfee08667235d0b3fe22137142e84
--- /dev/null
+++ b/test/implicit/1p/test_box1p.input
@@ -0,0 +1,45 @@
+###############################################################
+# Parameter file for test_1p.
+# Everything behind a '#' is a comment.
+# Type "./test_1p --help" for more information.
+###############################################################
+
+###############################################################
+# Mandatory arguments
+###############################################################
+
+[TimeManager]
+DtInitial = 1 # [s]
+TEnd = 1 # [s]
+
+[Grid]
+File = ./grids/test_1p_2d.dgf
+
+[Problem]
+Name = 1ptestbox # name passed to the output routines
+
+[SpatialParams]
+# Lens from (0.25|0.25) to (0.75|0.75)
+# Box-wise definition of permability
+LensLowerLeftX = 0.25 # [m]
+LensLowerLeftY = 0.25 # [m]
+LensLowerLeftZ = 0.25 # [m]
+LensUpperRightX = 0.75 # [m]
+LensUpperRightY = 0.75 # [m]
+LensUpperRightZ = 0.75 # [m]
+
+# Permeabilities
+Permeability = 1e-10 # [m^2]
+PermeabilityLens = 1e-12 # [m^2]
+
+###############################################################
+# Simulation restart
+#
+# DuMux simulations can be restarted from *.drs files
+# Set Restart to the value of a specific file,
+# e.g.: 'Restart = 27184.1' for the restart file
+# name_time=27184.1_rank=0.drs
+# Please comment in the two lines below, if restart is desired.
+###############################################################
+# [TimeManager]
+# Restart = ...
diff --git a/test/implicit/1p/test_1p.cc b/test/implicit/1p/test_cc1p.cc
similarity index 98%
rename from test/implicit/1p/test_1p.cc
rename to test/implicit/1p/test_cc1p.cc
index 81bb87c96a3801e59279a670647eee67169859d1..488461f0c858b2822644661f72c8d453a94fd9ab 100644
--- a/test/implicit/1p/test_1p.cc
+++ b/test/implicit/1p/test_cc1p.cc
@@ -59,6 +59,6 @@ void usage(const char *progName, const std::string &errorMsg)
int main(int argc, char** argv)
{
- typedef TTAG(OnePTestProblem) ProblemTypeTag;
+ typedef TTAG(OnePTestCCProblem) ProblemTypeTag;
return Dumux::start<ProblemTypeTag>(argc, argv, usage);
}
diff --git a/test/implicit/1p/test_1p.input b/test/implicit/1p/test_cc1p.input
similarity index 95%
rename from test/implicit/1p/test_1p.input
rename to test/implicit/1p/test_cc1p.input
index 8e0bfa3599f284364b0e657e89bd22dfedc693e6..6c9fb69a43165fa4da0370eec9fdd5832fa836cf 100644
--- a/test/implicit/1p/test_1p.input
+++ b/test/implicit/1p/test_cc1p.input
@@ -15,6 +15,9 @@ TEnd = 1 # [s]
[Grid]
File = ./grids/test_1p_2d.dgf
+[Problem]
+Name = 1ptestcc # name passed to the output routines
+
[SpatialParams]
# Lens from (0.25|0.25) to (0.75|0.75)
# Box-wise definition of permability
diff --git a/test/implicit/1p2c/1p2coutflowproblem.hh b/test/implicit/1p2c/1p2coutflowproblem.hh
index 6291ed94c48b2abdf5c155c789f781b0f71aa370..195b6ebb5bfe953fe55a221b26cbec8b636c0c94 100644
--- a/test/implicit/1p2c/1p2coutflowproblem.hh
+++ b/test/implicit/1p2c/1p2coutflowproblem.hh
@@ -44,12 +44,14 @@ class OnePTwoCOutflowProblem;
namespace Properties
{
-NEW_TYPE_TAG(OnePTwoCOutflowProblem, INHERITS_FROM(BoxOnePTwoC));
+NEW_TYPE_TAG(OnePTwoCOutflowProblem, INHERITS_FROM(OnePTwoC));
+NEW_TYPE_TAG(OnePTwoCOutflowBoxProblem, INHERITS_FROM(BoxModel, OnePTwoCOutflowProblem));
+NEW_TYPE_TAG(OnePTwoCOutflowCCProblem, INHERITS_FROM(CCModel, OnePTwoCOutflowProblem));
// Set the grid type
SET_PROP(OnePTwoCOutflowProblem, Grid)
{
-#if HAVE_UG
+#if 0 //HAVE_UG
typedef Dune::UGGrid<2> type;
#else
typedef Dune::SGrid<2, 2> type;
@@ -153,6 +155,11 @@ public:
{
//initialize fluid system
FluidSystem::init();
+
+ name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag,
+ std::string,
+ Problem,
+ Name);
}
/*!
@@ -166,7 +173,9 @@ public:
* This is used as a prefix for files generated by the simulation.
*/
const char *name() const
- { return "outflow"; }
+ {
+ return name_.c_str();
+ }
/*!
* \brief Returns the temperature within the domain.
@@ -186,35 +195,46 @@ public:
/*!
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary segment.
+ *
+ * \param values The boundary types for the conservation equations
+ * \param globalPos The position for which the bc type should be evaluated
*/
- void boundaryTypes(BoundaryTypes &values, const Vertex &vertex) const
+ void boundaryTypesAtPos(BoundaryTypes &values,
+ const GlobalPosition &globalPos) const
{
- const GlobalPosition globalPos = vertex.geometry().center();
-
if(globalPos[0] < eps_ || globalPos[0] > this->bboxMax()[0] - eps_)
+ {
values.setAllDirichlet();
+ std::cout << "set Dirichlet values at pos " << globalPos << std::endl;
+ }
else
+ {
values.setAllNeumann();
-
+ std::cout << "set Neumann values at pos " << globalPos << std::endl;
+ }
+
//outflow condition for the transport equation at right boundary
- if(globalPos[0] > this->bboxMax()[0] - eps_)
- values.setOutflow(transportEqIdx);
+ //if(globalPos[0] > this->bboxMax()[0] - eps_)
+ // values.setOutflow(transportEqIdx);
}
/*!
- * \brief Evaluate the boundary conditions for a Dirichlet
+ * \brief Evaluate the boundary conditions for a dirichlet
* boundary segment.
*
- * For this method, the \a priVars parameter stores primary variables.
+ * \param values The dirichlet values for the primary variables
+ * \param globalPos The position for which the bc type should be evaluated
+ *
+ * For this method, the \a values parameter stores primary variables.
*/
- void dirichlet(PrimaryVariables &priVars, const Vertex &vertex) const
+ void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
{
- const GlobalPosition globalPos = vertex.geometry().center();
-
- initial_(priVars, globalPos);
+ std::cout << "set Dirichlet values " << values << " at pos " << globalPos << std::endl;
+ initial_(values, globalPos);
+ std::cout << "set Dirichlet values " << values << " at pos " << globalPos << std::endl;
//condition for the N2 molefraction at left boundary
- if(globalPos[0] < eps_)
- priVars[massOrMoleFracIdx] = 2.0e-5;
+ //if (globalPos[0] < 0.5*(this->bboxMin()[0] + this->bboxMax()[0]))
+ // values[massOrMoleFracIdx] = 2.0e-5;
}
/*!
@@ -232,7 +252,6 @@ public:
const int scvIdx,
const int boundaryFaceIdx) const
{
- //const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
priVars = 0;
}
@@ -261,18 +280,15 @@ public:
/*!
* \brief Evaluate the initial value for a control volume.
*
- * For this method, the \a priVars parameter stores primary
+ * \param values The initial values for the primary variables
+ * \param globalPos The position for which the initial condition should be evaluated
+ *
+ * For this method, the \a values parameter stores primary
* variables.
*/
- void initial(PrimaryVariables &priVars,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- const int scvIdx) const
+ void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
{
- const GlobalPosition &globalPos
- = element.geometry().corner(scvIdx);
-
- initial_(priVars, globalPos);
+ initial_(values, globalPos);
}
// \}
@@ -282,11 +298,12 @@ private:
void initial_(PrimaryVariables &priVars,
const GlobalPosition &globalPos) const
{
- priVars[pressureIdx] = 2e5 - 1e5*globalPos[0];//0.0; //initial condition for the pressure
- priVars[massOrMoleFracIdx] = 0.0; //initial condition for the N2 molefraction
+ priVars[pressureIdx] = 2e5 - 1e5*globalPos[0]; // initial condition for the pressure
+ priVars[massOrMoleFracIdx] = 0.0; // initial condition for the N2 molefraction
}
const Scalar eps_;
+ std::string name_;
};
} //end namespace
diff --git a/test/implicit/1p2c/Makefile.am b/test/implicit/1p2c/Makefile.am
index c8d0376f870cb901581b203cad25c713916f13da..60d4b1539c6c9682f095ac2fc9ae35e6ed9b926d 100644
--- a/test/implicit/1p2c/Makefile.am
+++ b/test/implicit/1p2c/Makefile.am
@@ -1,8 +1,9 @@
-check_PROGRAMS = test_1p2c
+check_PROGRAMS = test_box1p2c test_cc1p2c
noinst_HEADERS = *.hh
EXTRA_DIST=*reference.vtu *.input grids/*.dgf CMakeLists.txt
-test_1p2c_SOURCES = test_1p2c.cc
+test_box1p2c_SOURCES = test_box1p2c.cc
+test_cc1p2c_SOURCES = test_cc1p2c.cc
include $(top_srcdir)/am/global-rules
diff --git a/test/implicit/1p2c/test_box1p2c.cc b/test/implicit/1p2c/test_box1p2c.cc
new file mode 100644
index 0000000000000000000000000000000000000000..a6a454645c98de27a89046a9b537481bcbc81fe9
--- /dev/null
+++ b/test/implicit/1p2c/test_box1p2c.cc
@@ -0,0 +1,57 @@
+// -*- 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 test for the 1p2c box model
+ */
+#include "config.h"
+#include "1p2coutflowproblem.hh"
+#include <dumux/common/start.hh>
+
+/*!
+ * \brief Provides an interface for customizing error messages associated with
+ * reading in parameters.
+ *
+ * \param progName The name of the program, that was tried to be started.
+ * \param errorMsg The error message that was issued by the start function.
+ * Comprises the thing that went wrong and a general help message.
+ */
+void usage(const char *progName, const std::string &errorMsg)
+{
+ if (errorMsg.size() > 0) {
+ std::string errorMessageOut = "\nUsage: ";
+ errorMessageOut += progName;
+ errorMessageOut += " [options]\n";
+ errorMessageOut += errorMsg;
+ errorMessageOut += "\n\nThe list of mandatory options for this program is:\n"
+ "\t-TimeManager.TEnd End of the simulation [s] \n"
+ "\t-TimeManager.DtInitial Initial timestep size [s] \n"
+ "\t-Grid.File Name of the file containing the grid \n"
+ "\t definition in DGF format\n";
+ std::cout << errorMessageOut
+ << "\n";
+ }
+}
+
+int main(int argc, char** argv)
+{
+ typedef TTAG(OnePTwoCOutflowBoxProblem) ProblemTypeTag;
+ return Dumux::start<ProblemTypeTag>(argc, argv, usage);
+}
diff --git a/test/implicit/1p2c/test_box1p2c.input b/test/implicit/1p2c/test_box1p2c.input
new file mode 100644
index 0000000000000000000000000000000000000000..e8346519c3f77cf3944eb190a5888808f9be265c
--- /dev/null
+++ b/test/implicit/1p2c/test_box1p2c.input
@@ -0,0 +1,31 @@
+###############################################################
+# Parameter file for test_1p2c.
+# Everything behind a '#' is a comment.
+# Type "./test_1p2c --help" for more information.
+###############################################################
+
+###############################################################
+# Mandatory arguments
+###############################################################
+
+[TimeManager]
+DtInitial = 1 # [s]
+TEnd = 100 # [s]
+
+[Grid]
+File = ./grids/test_1p2c.dgf
+
+[Problem]
+Name = outflowbox # name passed to the output routines
+
+###############################################################
+# Simulation restart
+#
+# DuMux simulations can be restarted from *.drs files
+# Set Restart to the value of a specific file,
+# e.g.: 'Restart = 27184.1' for the restart file
+# name_time=27184.1_rank=0.drs
+# Please comment in the two lines below, if restart is desired.
+###############################################################
+# [TimeManager]
+# Restart = ...
diff --git a/test/implicit/1p2c/test_1p2c.cc b/test/implicit/1p2c/test_cc1p2c.cc
similarity index 97%
rename from test/implicit/1p2c/test_1p2c.cc
rename to test/implicit/1p2c/test_cc1p2c.cc
index d0a4d6a9432a6b1731b9cba59dd6d35df0f3a0aa..50ef574f77e3137f7597fbb19ce2a61a65e08b68 100644
--- a/test/implicit/1p2c/test_1p2c.cc
+++ b/test/implicit/1p2c/test_cc1p2c.cc
@@ -52,6 +52,6 @@ void usage(const char *progName, const std::string &errorMsg)
int main(int argc, char** argv)
{
- typedef TTAG(OnePTwoCOutflowProblem) ProblemTypeTag;
+ typedef TTAG(OnePTwoCOutflowCCProblem) ProblemTypeTag;
return Dumux::start<ProblemTypeTag>(argc, argv, usage);
}
diff --git a/test/implicit/1p2c/test_1p2c.input b/test/implicit/1p2c/test_cc1p2c.input
similarity index 93%
rename from test/implicit/1p2c/test_1p2c.input
rename to test/implicit/1p2c/test_cc1p2c.input
index 1c0887cac9a1382a27ce0ee989116b9fff19f12b..7578c20cb5ead014431ca7742e032f0c5a7ef45e 100644
--- a/test/implicit/1p2c/test_1p2c.input
+++ b/test/implicit/1p2c/test_cc1p2c.input
@@ -15,6 +15,9 @@ TEnd = 100 # [s]
[Grid]
File = ./grids/test_1p2c.dgf
+[Problem]
+Name = outflowcc # name passed to the output routines
+
###############################################################
# Simulation restart
#
diff --git a/test/implicit/2p/test_box2p.input b/test/implicit/2p/test_box2p.input
index 6a8c23673a65996fda3a0e8b69db31c24e3f0f70..f92213bafa6708d69a71e1e7e2c5dab0361e4efa 100644
--- a/test/implicit/2p/test_box2p.input
+++ b/test/implicit/2p/test_box2p.input
@@ -22,7 +22,7 @@ LensUpperRightX = 4.0 # [m] x-coordinate of the upper right lens corner
LensUpperRightY = 3.0 # [m] y-coordinate of the upper right lens corner
[Problem]
-Name = lensbox # for output files
+Name = lensbox # name passed to the output routines
###############################################################
# Simulation restart
diff --git a/test/implicit/2p/test_cc2p.input b/test/implicit/2p/test_cc2p.input
index a63fa6c64ec03bf90e441ec23f724dc9064d9bec..0d17b547e1c16569d6d77747fdc774d3f86a5ee9 100644
--- a/test/implicit/2p/test_cc2p.input
+++ b/test/implicit/2p/test_cc2p.input
@@ -22,7 +22,7 @@ LensUpperRightX = 4.0 # [m] x-coordinate of the upper right lens corner
LensUpperRightY = 3.0 # [m] y-coordinate of the upper right lens corner
[Problem]
-Name = lenscc # for output files
+Name = lenscc # name passed to the output routines
###############################################################
# Simulation restart
diff --git a/test/implicit/2p2c/Makefile.am b/test/implicit/2p2c/Makefile.am
index 5acdf0abf9cbda548062f7869dbaf6e6091b02f8..64ecc401981d6115285230d2373567089399c419 100644
--- a/test/implicit/2p2c/Makefile.am
+++ b/test/implicit/2p2c/Makefile.am
@@ -1,8 +1,9 @@
-check_PROGRAMS = test_2p2c
+check_PROGRAMS = test_box2p2c test_cc2p2c
noinst_HEADERS = *.hh
EXTRA_DIST=*reference.vtu *.input grids/*.dgf CMakeLists.txt
-test_2p2c_SOURCES = test_2p2c.cc
+test_box2p2c_SOURCES = test_box2p2c.cc
+test_cc2p2c_SOURCES = test_cc2p2c.cc
include $(top_srcdir)/am/global-rules
diff --git a/test/implicit/2p2c/injectionproblem.hh b/test/implicit/2p2c/injectionproblem.hh
index 6faef0455e261febffb53e5e8049cfa25665b9f3..1ad5ae4d16567929ff1ddd8b53aaab270afe23cf 100644
--- a/test/implicit/2p2c/injectionproblem.hh
+++ b/test/implicit/2p2c/injectionproblem.hh
@@ -40,8 +40,10 @@ class InjectionProblem;
namespace Properties
{
-NEW_TYPE_TAG(InjectionProblem, INHERITS_FROM(BoxTwoPTwoC, InjectionSpatialParams));
-
+NEW_TYPE_TAG(InjectionProblem, INHERITS_FROM(TwoPTwoC, InjectionSpatialParams));
+NEW_TYPE_TAG(InjectionBoxProblem, INHERITS_FROM(BoxModel, InjectionProblem));
+NEW_TYPE_TAG(InjectionCCProblem, INHERITS_FROM(CCModel, InjectionProblem));
+
// Set the grid type
SET_PROP(InjectionProblem, Grid)
{
@@ -130,6 +132,8 @@ class InjectionProblem : public ImplicitPorousMediaProblem<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, GridCreator) GridCreator;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
public:
/*!
@@ -243,12 +247,11 @@ public:
* used for which equation on a given boundary segment.
*
* \param values The boundary types for the conservation equations
- * \param vertex The vertex for which the boundary type is set
+ * \param globalPos The position for which the bc type should be evaluated
*/
- void boundaryTypes(BoundaryTypes &values, const Vertex &vertex) const
+ void boundaryTypesAtPos(BoundaryTypes &values,
+ const GlobalPosition &globalPos) const
{
- const GlobalPosition globalPos = vertex.geometry().center();
-
if (globalPos[0] < eps_)
values.setAllDirichlet();
else
@@ -260,14 +263,12 @@ public:
* boundary segment.
*
* \param values The dirichlet values for the primary variables
- * \param vertex The vertex for which the boundary type is set
+ * \param globalPos The position for which the bc type should be evaluated
*
* For this method, the \a values parameter stores primary variables.
*/
- void dirichlet(PrimaryVariables &values, const Vertex &vertex) const
+ void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
{
- const GlobalPosition globalPos = vertex.geometry().center();
-
initial_(values, globalPos);
}
@@ -292,10 +293,15 @@ public:
int scvIdx,
int boundaryFaceIdx) const
{
- const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
-
values = 0;
- if (globalPos[1] < 15 && globalPos[1] > 5) {
+
+ GlobalPosition globalPos;
+ if (isBox)
+ globalPos = element.geometry().corner(scvIdx);
+ else
+ globalPos = is.geometry().center();
+
+ if (globalPos[1] < 15 && globalPos[1] > 7) {
values[contiN2EqIdx] = -1e-3; // kg/(s*m^2)
}
}
@@ -311,20 +317,13 @@ public:
* \brief Evaluate the initial value for a control volume.
*
* \param values The initial values for the primary variables
- * \param element The finite element
- * \param fvGeometry The finite-volume geometry in the box scheme
- * \param scvIdx The local vertex index
+ * \param globalPos The position for which the initial condition should be evaluated
*
* For this method, the \a values parameter stores primary
* variables.
*/
- void initial(PrimaryVariables &values,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- int scvIdx) const
+ void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
{
- const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
-
initial_(values, globalPos);
}
diff --git a/test/implicit/2p2c/test_box2p2c.cc b/test/implicit/2p2c/test_box2p2c.cc
new file mode 100644
index 0000000000000000000000000000000000000000..03dd3a88b5caa1c24b589a51f5ea1e5fa2a099fc
--- /dev/null
+++ b/test/implicit/2p2c/test_box2p2c.cc
@@ -0,0 +1,67 @@
+// -*- 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 Test for the 2p2c box model.
+ */
+#include "config.h"
+#include "injectionproblem.hh"
+#include <dumux/common/start.hh>
+
+/*!
+ * \brief Provides an interface for customizing error messages associated with
+ * reading in parameters.
+ *
+ * \param progName The name of the program, that was tried to be started.
+ * \param errorMsg The error message that was issued by the start function.
+ * Comprises the thing that went wrong and a general help message.
+ */
+void usage(const char *progName, const std::string &errorMsg)
+{
+ if (errorMsg.size() > 0) {
+ std::string errorMessageOut = "\nUsage: ";
+ errorMessageOut += progName;
+ errorMessageOut += " [options]\n";
+ errorMessageOut += errorMsg;
+ errorMessageOut += "\n\nThe list of mandatory options for this program is:\n"
+ "\t-TimeManager.TEnd End of the simulation [s] \n"
+ "\t-TimeManager.DtInitial Initial timestep size [s] \n"
+ "\t-Grid.File Name of the file containing the grid \n"
+ "\t definition in DGF format\n"
+ "\t-FluidSystem.NTemperature Number of tabularization entries [-] \n"
+ "\t-FluidSystem.NPressure Number of tabularization entries [-] \n"
+ "\t-FluidSystem.PressureLow Low end for tabularization of fluid properties [Pa] \n"
+ "\t-FluidSystem.PressureHigh High end for tabularization of fluid properties [Pa] \n"
+ "\t-FluidSystem.TemperatureLow Low end for tabularization of fluid properties [Pa] \n"
+ "\t-FluidSystem.TemperatureHigh High end for tabularization of fluid properties [Pa] \n"
+ "\t-SimulationControl.Name The name of the output files [-] \n"
+ "\t-InitialConditions.Temperature Initial temperature in the reservoir [K] \n"
+ "\t-InitialConditions.DepthBOR Depth below ground surface [m] \n";
+
+ std::cout << errorMessageOut
+ << "\n";
+ }
+}
+
+int main(int argc, char** argv)
+{
+ typedef TTAG(InjectionBoxProblem) ProblemTypeTag;
+ return Dumux::start<ProblemTypeTag>(argc, argv, usage);
+}
diff --git a/test/implicit/2p2c/test_box2p2c.input b/test/implicit/2p2c/test_box2p2c.input
new file mode 100644
index 0000000000000000000000000000000000000000..103b865b9f4f6f4ecafae5299f68e1b7d02497e1
--- /dev/null
+++ b/test/implicit/2p2c/test_box2p2c.input
@@ -0,0 +1,42 @@
+###############################################################
+# Parameter file for test_2p2c.
+# Everything behind a '#' is a comment.
+# Type "./test_2p2c --help" for more information.
+###############################################################
+
+###############################################################
+# Mandatory arguments
+###############################################################
+
+[TimeManager]
+DtInitial = 250 # [s]
+TEnd = 1e4 # [s]
+
+[Grid]
+File = ./grids/test_2p2c.dgf
+
+[Problem]
+Name = injectionbox # [-] the name of the output files
+
+InitialTemperature = 313.15 # [K] 40°C initial temperature in the reservoir
+DepthBOR = 2700 # [m] depth below ground surface
+
+# tabularization
+NTemperature = 3 # [-] number of temperature table entries
+NPressure = 200 # [-] number of pressure table entries
+PressureLow = 1e5 # [Pa] lower pressure limit for tabularization
+PressureHigh = 3e7 # [Pa] upper pressure limit for tabularization
+TemperatureLow = 312.15 # [Pa] lower temperature limit for tabularization
+TemperatureHigh = 314.15 # [Pa] upper temperature limit for tabularization
+
+###############################################################
+# Simulation restart
+#
+# DuMux simulations can be restarted from *.drs files
+# Set Restart to the value of a specific file,
+# e.g.: 'Restart = 27184.1' for the restart file
+# name_time=27184.1_rank=0.drs
+# Please comment in the two lines below, if restart is desired.
+###############################################################
+# [TimeManager]
+# Restart = ...
diff --git a/test/implicit/2p2c/test_2p2c.cc b/test/implicit/2p2c/test_cc2p2c.cc
similarity index 98%
rename from test/implicit/2p2c/test_2p2c.cc
rename to test/implicit/2p2c/test_cc2p2c.cc
index 0c1daa929ae1eee6011d1fcb3f8d075a9fbe8bf6..48fb4ba1faa41100678c56838c3dad0dd177a1e5 100644
--- a/test/implicit/2p2c/test_2p2c.cc
+++ b/test/implicit/2p2c/test_cc2p2c.cc
@@ -62,6 +62,6 @@ void usage(const char *progName, const std::string &errorMsg)
int main(int argc, char** argv)
{
- typedef TTAG(InjectionProblem) ProblemTypeTag;
+ typedef TTAG(InjectionCCProblem) ProblemTypeTag;
return Dumux::start<ProblemTypeTag>(argc, argv, usage);
}
diff --git a/test/implicit/2p2c/test_2p2c.input b/test/implicit/2p2c/test_cc2p2c.input
similarity index 96%
rename from test/implicit/2p2c/test_2p2c.input
rename to test/implicit/2p2c/test_cc2p2c.input
index e5d98a4cd0f5a7539f449800dd098bef95c1d939..e1aed4b66009fd6e905e22abc72920e290f2f60e 100644
--- a/test/implicit/2p2c/test_2p2c.input
+++ b/test/implicit/2p2c/test_cc2p2c.input
@@ -16,7 +16,7 @@ TEnd = 1e4 # [s]
File = ./grids/test_2p2c.dgf
[Problem]
-Name = injection # [-] the name of the output files
+Name = injectioncc # [-] the name of the output files
InitialTemperature = 313.15 # [K] 40°C initial temperature in the reservoir
DepthBOR = 2700 # [m] depth below ground surface