From 39b4e961abcd645dade3d8301342e08fd971c315 Mon Sep 17 00:00:00 2001
From: Andreas Lauser <and@poware.org>
Date: Fri, 2 Sep 2011 11:54:48 +0000
Subject: [PATCH] 2p2c box model: make it compile with Scalar != double
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@6568 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/boxmodels/2p2c/2p2cfluxvariables.hh | 9 ++--
dumux/boxmodels/2p2c/2p2clocalresidual.hh | 22 ++++++----
dumux/boxmodels/2p2c/2p2cmodel.hh | 46 ++++++++++-----------
dumux/boxmodels/2p2c/2p2cproblem.hh | 12 +++---
dumux/material/binarycoefficients/h2o_n2.hh | 2 +-
test/boxmodels/2p2c/injectionproblem.hh | 9 ++--
6 files changed, 52 insertions(+), 48 deletions(-)
diff --git a/dumux/boxmodels/2p2c/2p2cfluxvariables.hh b/dumux/boxmodels/2p2c/2p2cfluxvariables.hh
index ad95409af8..156212be94 100644
--- a/dumux/boxmodels/2p2c/2p2cfluxvariables.hh
+++ b/dumux/boxmodels/2p2c/2p2cfluxvariables.hh
@@ -70,8 +70,8 @@ class TwoPTwoCFluxVariables
typedef typename FVElementGeometry::SubControlVolume SCV;
typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
- typedef Dune::FieldVector<CoordScalar, dimWorld> Vector;
- typedef Dune::FieldMatrix<CoordScalar, dimWorld, dimWorld> Tensor;
+ typedef Dune::FieldVector<Scalar, dim> Vector;
+ typedef Dune::FieldMatrix<Scalar, dim, dim> Tensor;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices)) Indices;
enum {
@@ -227,7 +227,10 @@ private:
fvElemGeom_,
face().j));
K.mv(potentialGrad_[phaseIdx], Kmvp_[phaseIdx]);
- KmvpNormal_[phaseIdx] = - (Kmvp_[phaseIdx] * face().normal);
+ KmvpNormal_[phaseIdx] = 0;
+ for (int i = 0; i < Vector::size; ++i)
+ KmvpNormal_[phaseIdx] += Kmvp_[phaseIdx][i] * face().normal[i];
+ KmvpNormal_[phaseIdx] *= -1;
}
// set the upstream and downstream vertices
diff --git a/dumux/boxmodels/2p2c/2p2clocalresidual.hh b/dumux/boxmodels/2p2c/2p2clocalresidual.hh
index 64b9c3c05e..19894552df 100644
--- a/dumux/boxmodels/2p2c/2p2clocalresidual.hh
+++ b/dumux/boxmodels/2p2c/2p2clocalresidual.hh
@@ -306,10 +306,13 @@ public:
void computeDiffusiveFlux(PrimaryVariables &flux, const FluxVariables &vars) const
{
// add diffusive flux of gas component in liquid phase
- Scalar tmp =
- - vars.porousDiffCoeff(lPhaseIdx) *
- vars.molarDensityAtIP(lPhaseIdx) *
- (vars.molarConcGrad(lPhaseIdx) * vars.face().normal);
+ Scalar tmp = 0;
+ for (int i = 0; i < dim; ++i)
+ tmp += vars.molarConcGrad(lPhaseIdx)[i] * vars.face().normal[i];
+ tmp *= -1;
+ tmp *=
+ vars.porousDiffCoeff(lPhaseIdx) *
+ vars.molarDensityAtIP(lPhaseIdx);
// add the diffusive fluxes only to the component mass balance
if (replaceCompEqIdx != contiGEqIdx)
flux[contiGEqIdx] += tmp * FluidSystem::molarMass(gCompIdx);
@@ -317,10 +320,13 @@ public:
flux[contiLEqIdx] -= tmp * FluidSystem::molarMass(lCompIdx);
// add diffusive flux of liquid component in gas phase
- tmp =
- - vars.porousDiffCoeff(gPhaseIdx) *
- vars.molarDensityAtIP(gPhaseIdx) *
- (vars.molarConcGrad(gPhaseIdx) * vars.face().normal);
+ tmp = 0;
+ for (int i = 0; i < dim; ++i)
+ tmp += vars.molarConcGrad(gPhaseIdx)[i] * vars.face().normal[i];
+ tmp *= -1;
+ tmp *=
+ vars.porousDiffCoeff(gPhaseIdx) *
+ vars.molarDensityAtIP(gPhaseIdx);
// add the diffusive fluxes only to the component mass balance
if (replaceCompEqIdx != contiLEqIdx)
flux[contiLEqIdx] += tmp * FluidSystem::molarMass(lCompIdx);
diff --git a/dumux/boxmodels/2p2c/2p2cmodel.hh b/dumux/boxmodels/2p2c/2p2cmodel.hh
index 4504f04105..8d06df6073 100644
--- a/dumux/boxmodels/2p2c/2p2cmodel.hh
+++ b/dumux/boxmodels/2p2c/2p2cmodel.hh
@@ -140,6 +140,7 @@ class TwoPTwoCModel: public BoxModel<TypeTag>
formulation = GET_PROP_VALUE(TypeTag, PTAG(Formulation))
};
+ typedef typename GridView::ctype CoordScalar;
typedef typename GridView::template Codim<dim>::Entity Vertex;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
@@ -285,36 +286,31 @@ public:
MultiWriter &writer)
{
bool velocityOutput = GET_PROP_VALUE(TypeTag, PTAG(EnableVelocityOutput));
- typedef Dune::BlockVector<Dune::FieldVector<Scalar, 1> > ScalarField;
- typedef Dune::BlockVector<Dune::FieldVector<Scalar, dim> > VectorField;
+ typedef Dune::BlockVector<Dune::FieldVector<double, 1> > ScalarField;
+ typedef Dune::BlockVector<Dune::FieldVector<double, dim> > VectorField;
// create the required scalar fields
unsigned numVertices = this->problem_().gridView().size(dim);
- ScalarField *Sg = writer.allocateManagedBuffer (numVertices);
- ScalarField *Sl = writer.allocateManagedBuffer (numVertices);
- ScalarField *pg = writer.allocateManagedBuffer (numVertices);
- ScalarField *pl = writer.allocateManagedBuffer (numVertices);
- ScalarField *pc = writer.allocateManagedBuffer (numVertices);
- ScalarField *rhoL =
- writer.allocateManagedBuffer (numVertices);
- ScalarField *rhoG =
- writer.allocateManagedBuffer (numVertices);
- ScalarField *mobL =
- writer.allocateManagedBuffer (numVertices);
- ScalarField *mobG =
- writer.allocateManagedBuffer (numVertices);
- ScalarField *phasePresence =
- writer.allocateManagedBuffer (numVertices);
+ ScalarField *Sg = writer.allocateManagedBuffer(numVertices);
+ ScalarField *Sl = writer.allocateManagedBuffer(numVertices);
+ ScalarField *pg = writer.allocateManagedBuffer(numVertices);
+ ScalarField *pl = writer.allocateManagedBuffer(numVertices);
+ ScalarField *pc = writer.allocateManagedBuffer(numVertices);
+ ScalarField *rhoL = writer.allocateManagedBuffer(numVertices);
+ ScalarField *rhoG = writer.allocateManagedBuffer(numVertices);
+ ScalarField *mobL = writer.allocateManagedBuffer(numVertices);
+ ScalarField *mobG = writer.allocateManagedBuffer(numVertices);
+ ScalarField *phasePresence = writer.allocateManagedBuffer(numVertices);
ScalarField *massFrac[numPhases][numComponents];
for (int i = 0; i < numPhases; ++i)
for (int j = 0; j < numComponents; ++j)
- massFrac[i][j] = writer.allocateManagedBuffer (numVertices);
+ massFrac[i][j] = writer.allocateManagedBuffer(numVertices);
ScalarField *temperature = writer.allocateManagedBuffer(numVertices);
ScalarField *poro = writer.allocateManagedBuffer(numVertices);
ScalarField *cellNum =writer.allocateManagedBuffer (numVertices);
- VectorField *velocityG = writer.template allocateManagedBuffer <Scalar, dim> (numVertices);
- VectorField *velocityL = writer.template allocateManagedBuffer <Scalar, dim> (numVertices);
+ VectorField *velocityG = writer.template allocateManagedBuffer<double, dim>(numVertices);
+ VectorField *velocityL = writer.template allocateManagedBuffer<double, dim>(numVertices);
if(velocityOutput) // check if velocity output is demanded
{
@@ -431,7 +427,7 @@ public:
const Dune::FieldVector<Scalar, dim>& localPosIP = fvElemGeom.subContVolFace[faceIdx].ipLocal;
// Transformation of the global normal vector to normal vector in the reference element
- const Dune::FieldMatrix<Scalar, dim, dim> jacobianT1 = elemIt->geometry().jacobianTransposed(localPosIP);
+ const Dune::FieldMatrix<CoordScalar, dim, dim> jacobianT1 = elemIt->geometry().jacobianTransposed(localPosIP);
const GlobalPosition globalNormal = fluxDat.face().normal;
GlobalPosition localNormal(0);
@@ -467,18 +463,18 @@ public:
}
typedef Dune::GenericReferenceElements<Scalar, dim> ReferenceElements;
- const Dune::FieldVector<Scalar, dim>& localPos = ReferenceElements::general(elemIt->geometry().type()).position(0,
- 0);
+ const Dune::FieldVector<Scalar, dim>& localPos =
+ ReferenceElements::general(elemIt->geometry().type()).position(0, 0);
// get the transposed Jacobian of the element mapping
- const Dune::FieldMatrix<Scalar, dim, dim> jacobianT2 = elemIt->geometry().jacobianTransposed(localPos);
+ const Dune::FieldMatrix<CoordScalar, dim, dim> 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<Scalar, dim> scvVelocity(0);
+ Dune::FieldVector<CoordScalar, dim> scvVelocity(0);
jacobianT2.mtv(scvVelocityL[i], scvVelocity);
scvVelocity /= elemIt->geometry().integrationElement(localPos);
diff --git a/dumux/boxmodels/2p2c/2p2cproblem.hh b/dumux/boxmodels/2p2c/2p2cproblem.hh
index 492bbc4d4f..bb1a910b15 100644
--- a/dumux/boxmodels/2p2c/2p2cproblem.hh
+++ b/dumux/boxmodels/2p2c/2p2cproblem.hh
@@ -54,7 +54,9 @@ class TwoPTwoCProblem : public BoxProblem<TypeTag>
dimWorld = GridView::dimensionworld
};
- typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ typedef typename GridView::ctype CoordScalar;
+ typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
+ typedef Dune::FieldVector<Scalar, dim> Vector;
typedef typename GridView::template Codim<0>::Entity Element;
public:
@@ -121,7 +123,7 @@ public:
* This is the box discretization specific interface. By default
* it just calls gravityAtPos().
*/
- const GlobalPosition &boxGravity(const Element &element,
+ const Vector &boxGravity(const Element &element,
const FVElementGeometry &fvGeom,
int scvIdx) const
{ return asImp_().gravityAtPos(fvGeom.subContVol[scvIdx].global); }
@@ -132,7 +134,7 @@ public:
* This is discretization independent interface. By default it
* just calls gravity().
*/
- const GlobalPosition &gravityAtPos(const GlobalPosition &pos) const
+ const Vector &gravityAtPos(const GlobalPosition &pos) const
{ return asImp_().gravity(); }
/*!
@@ -144,7 +146,7 @@ public:
* property is true, \f$\boldsymbol{g} = ( 0,\dots,\ -9.81)^T \f$ holds,
* else \f$\boldsymbol{g} = ( 0,\dots, 0)^T \f$.
*/
- const GlobalPosition &gravity() const
+ const Vector &gravity() const
{ return gravity_; }
/*!
@@ -170,7 +172,7 @@ private:
const Implementation &asImp_() const
{ return *static_cast<const Implementation *>(this); }
- GlobalPosition gravity_;
+ Vector gravity_;
// spatial parameters
SpatialParameters spatialParams_;
diff --git a/dumux/material/binarycoefficients/h2o_n2.hh b/dumux/material/binarycoefficients/h2o_n2.hh
index dadfb64abf..a30c7ef16f 100644
--- a/dumux/material/binarycoefficients/h2o_n2.hh
+++ b/dumux/material/binarycoefficients/h2o_n2.hh
@@ -75,7 +75,7 @@ public:
// atomic diffusion volumes
const Scalar SigmaNu[2] = { 13.1 /* H2O */, 18.5 /* N2 */ };
// molar masses [g/mol]
- const Scalar M[2] = { H2O::molarMass()*1e3, N2::molarMass()*1e3 };
+ const Scalar M[2] = { H2O::molarMass()*Scalar(1e3), N2::molarMass()*Scalar(1e3) };
return fullerMethod(M, SigmaNu, temperature, pressure);
};
diff --git a/test/boxmodels/2p2c/injectionproblem.hh b/test/boxmodels/2p2c/injectionproblem.hh
index 494bbc1405..f8e1f7d773 100644
--- a/test/boxmodels/2p2c/injectionproblem.hh
+++ b/test/boxmodels/2p2c/injectionproblem.hh
@@ -24,9 +24,8 @@
*
* \brief Definition of a problem, where air is injected under a low permeable layer.
*/
-
-#ifndef DUMUX_INJECTIONPROBLEM_HH
-#define DUMUX_INJECTIONPROBLEM_HH
+#ifndef DUMUX_INJECTION_PROBLEM_HH
+#define DUMUX_INJECTION_PROBLEM_HH
#include <dune/grid/io/file/dgfparser/dgfug.hh>
#include <dune/grid/io/file/dgfparser/dgfs.hh>
@@ -41,7 +40,6 @@
#include "injectionspatialparameters.hh"
-
namespace Dumux
{
@@ -65,8 +63,7 @@ SET_PROP(InjectionProblem, Problem)
};
// Set fluid configuration
-SET_PROP(InjectionProblem,
- FluidSystem)
+SET_PROP(InjectionProblem, FluidSystem)
{
//typedef Dumux::Brine_CO2_System<TypeTag, Dumux::IFP::CO2Tables> type;
typedef Dumux::H2O_N2_System<TypeTag> type;
--
GitLab