diff --git a/dumux/boxmodels/1p2c/1p2cboundaryvariables.hh b/dumux/boxmodels/1p2c/1p2cboundaryvariables.hh index 340499a2a379be06b588526aeee6959042eb2113..d888a78cffe7b7f0b3c954ea63182b897287e5bf 100644 --- a/dumux/boxmodels/1p2c/1p2cboundaryvariables.hh +++ b/dumux/boxmodels/1p2c/1p2cboundaryvariables.hh @@ -264,7 +264,10 @@ protected: calculateK_(problem.spatialParameters().intrinsicPermeability(element, fvElemGeom_, scvIdx_)); ScalarGradient Kmvp; K_.mv(potentialGrad_, Kmvp); - KmvpNormal_ = - (Kmvp*boundaryFace_->normal); + KmvpNormal_ = 0; + for (int i = 0; i < dim; ++i) + KmvpNormal_ = Kmvp[i]*boundaryFace_->normal[i]; + KmvpNormal_ *= -1; const VolumeVariables &vertDat = elemDat[scvIdx_]; diff --git a/dumux/boxmodels/1p2c/1p2cfluxvariables.hh b/dumux/boxmodels/1p2c/1p2cfluxvariables.hh index d8a1db4bcc25aeef363e9ac52e7e9658b7964f13..770f206455df92a11663d00d14c64778f3679d47 100644 --- a/dumux/boxmodels/1p2c/1p2cfluxvariables.hh +++ b/dumux/boxmodels/1p2c/1p2cfluxvariables.hh @@ -62,7 +62,9 @@ class OnePTwoCFluxVariables enum { dim = GridView::dimension }; enum { dimWorld = GridView::dimensionworld }; - + + typedef typename GridView::ctype CoordScalar; + typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition; typedef Dune::FieldVector<Scalar, dimWorld> Vector; typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> Tensor; @@ -325,8 +327,10 @@ protected: } else { // use two-point gradients - tmp = element.geometry().corner(face().i); - tmp -= element.geometry().corner(face().j); + const GlobalPosition &posI = element.geometry().corner(face().i); + const GlobalPosition &posJ = element.geometry().corner(face().j); + for (int i = 0; i < Vector::size; ++ i) + tmp[i] = posI[i] - posJ[i]; Scalar dist = tmp.two_norm(); tmp = face().normal; @@ -398,11 +402,14 @@ protected: * \param elemDat The parameters stored in the considered element */ void calculateVelocities_(const Problem &problem, - const Element &element, - const ElementVolumeVariables &elemDat) + const Element &element, + const ElementVolumeVariables &elemDat) { K_.mv(potentialGrad_, Kmvp_); - KmvpNormal_ = - (Kmvp_ * face().normal); + KmvpNormal_ = 0; + for (int i = 0; i < Vector::size; ++i) + KmvpNormal_ += Kmvp_[i] * face().normal[i]; + KmvpNormal_ *= -1; // set the upstream and downstream vertices upstreamIdx_ = face().i; diff --git a/dumux/boxmodels/1p2c/1p2clocalresidual.hh b/dumux/boxmodels/1p2c/1p2clocalresidual.hh index ca2c88699f1a5ce08caeae0a6c6fb9627ed98928..a5789a8ced0fe7a0c37afd676ff06ac13dff6bed 100644 --- a/dumux/boxmodels/1p2c/1p2clocalresidual.hh +++ b/dumux/boxmodels/1p2c/1p2clocalresidual.hh @@ -247,18 +247,24 @@ public: if(!useMoles) { // diffusive flux of the second component - massfraction - tmp = - fluxVars.porousDiffCoeff() * fluxVars.densityAtIP()* - (fluxVars.massFracGrad(comp1Idx) * fluxVars.face().normal); + tmp = 0; + for (int i = 0; i < Vector::size; ++ i) + tmp += fluxVars.massFracGrad(comp1Idx)[i]*fluxVars.face().normal[i]; + tmp *= -1; + tmp *= fluxVars.porousDiffCoeff() * fluxVars.densityAtIP(); flux[transEqIdx] += tmp;// * FluidSystem::molarMass(comp1Idx); } else { // diffusive flux of the second component - molefraction - tmp = - fluxVars.porousDiffCoeff() * fluxVars.molarDensityAtIP()* - (fluxVars.moleFracGrad(comp1Idx) * fluxVars.face().normal); - - // dispersive flux of second component - molefraction + tmp = 0; + for (int i = 0; i < Vector::size; ++ i) + tmp += fluxVars.moleFracGrad(comp1Idx)[i]*fluxVars.face().normal[i]; + tmp *= -1; + tmp *= fluxVars.porousDiffCoeff() * fluxVars.molarDensityAtIP(); + + // dispersive flux of second component - molefraction // Vector normalDisp; // fluxVars.dispersionTensor().mv(fluxVars.face().normal, normalDisp); // tmp -= fluxVars.molarDensityAtIP()* @@ -437,8 +443,12 @@ protected: *vertVars.fluidState().massFrac(phaseIdx, comp1Idx); // diffusive flux of comp1 component in phase0 - Scalar tmp = -boundaryVars.porousDiffCoeff()*boundaryVars.densityAtIP() - *(boundaryVars.massFracGrad(comp1Idx)*boundaryVars.boundaryFace().normal); + Scalar tmp = 0; + for (int i = 0; i < Vector::size; ++ i) + tmp += boundaryVars.massFracGrad(comp1Idx)[i]*boundaryVars.boundaryFace().normal[i]; + tmp *= -1; + + tmp *= boundaryVars.porousDiffCoeff()*boundaryVars.densityAtIP(); flux[transEqIdx] += tmp;//* FluidSystem::molarMass(comp1Idx); } else //use molefractions @@ -448,8 +458,12 @@ protected: *vertVars.fluidState().moleFrac(phaseIdx, comp1Idx); // diffusive flux of comp1 component in phase0 - Scalar tmp = -boundaryVars.porousDiffCoeff()*boundaryVars.molarDensityAtIP() - *(boundaryVars.moleFracGrad(comp1Idx)*boundaryVars.boundaryFace().normal); + Scalar tmp = 0; + for (int i = 0; i < Vector::size; ++ i) + tmp += boundaryVars.moleFracGrad(comp1Idx)[i]*boundaryVars.boundaryFace().normal[i]; + tmp *= -1; + + tmp *= boundaryVars.porousDiffCoeff()*boundaryVars.molarDensityAtIP(); flux[transEqIdx] += tmp; } } diff --git a/dumux/boxmodels/1p2c/1p2cmodel.hh b/dumux/boxmodels/1p2c/1p2cmodel.hh index e646b563a9f447179fc6a968e552ea896d02bff4..cbbd4263949f223485c7781f68207533c5193fcc 100644 --- a/dumux/boxmodels/1p2c/1p2cmodel.hh +++ b/dumux/boxmodels/1p2c/1p2cmodel.hh @@ -130,7 +130,7 @@ public: void addOutputVtkFields(const SolutionVector &sol, MultiWriter &writer) { - typedef Dune::BlockVector<Dune::FieldVector<Scalar, 1> > ScalarField; + typedef Dune::BlockVector<Dune::FieldVector<double, 1> > ScalarField; // create the required scalar fields unsigned numVertices = this->problem_().gridView().size(dim); @@ -167,7 +167,7 @@ public: #endif unsigned numElements = this->gridView_().size(0); ScalarField &rank = - *writer.allocateManagedBuffer (numElements); + *writer.allocateManagedBuffer(numElements); FVElementGeometry fvElemGeom; VolumeVariables volVars; diff --git a/dumux/boxmodels/1p2c/1p2cproblem.hh b/dumux/boxmodels/1p2c/1p2cproblem.hh index a75340488cbf7430f12b1b9dc85463ca76ff796a..c5caf34a4c0a579a33f675667d7e9ab785be1bdc 100644 --- a/dumux/boxmodels/1p2c/1p2cproblem.hh +++ b/dumux/boxmodels/1p2c/1p2cproblem.hh @@ -55,8 +55,10 @@ class OnePTwoCBoxProblem : public BoxProblem<TypeTag> dim = GridView::dimension, 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: @@ -123,7 +125,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); } @@ -134,7 +136,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(); } /*! @@ -146,7 +148,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_; } /*! @@ -171,7 +173,7 @@ private: const Implementation &asImp_() const { return *static_cast<const Implementation *>(this); } - GlobalPosition gravity_; + Vector gravity_; // spatial parameters SpatialParameters spatialParams_;