Skip to content
GitLab
Commit bde6eb7b authored by moirap's avatar moirap
Browse files

[2pTracer] final version of 2p tracer implementation

parent 4174fa94
Hide whitespace changes
Inline Side-by-side
dumux/porousmediumflow/2p/vtkoutputfields.hh
View file @ bde6eb7b
......@@ -50,6 +50,7 @@ public:
vtk.addVolumeVariable([](const auto& v){ return v.mobility(FluidSystem::phase1Idx); },"mob_"+ FluidSystem::phaseName(FluidSystem::phase1Idx));
vtk.addVolumeVariable([](const auto& v){ return v.porosity(); }, "porosity");
vtk.addVolumeVariable([](const auto& v){ return v.permeability(); }, "permeability");
}
};
......
test/porousmediumflow/tracer/2ptracer/2pncimmiscible.hh
View file @ bde6eb7b
......@@ -81,6 +81,15 @@ public:
//TODO: for each tracercomponent, define a unique component index for all added tracers
static constexpr int comp2Idx = 2; //!< index of the first tracer
// static constexpr int comp3Idx = 3; //!< index of the second tracer
// static constexpr int comp4Idx = 4; //!< index of the ... tracer
// static constexpr int comp5Idx = 5;
// static constexpr int comp6Idx = 6;
// static constexpr int comp7Idx = 7;
// static constexpr int comp8Idx = 8;
// static constexpr int comp9Idx = 9;
// static constexpr int comp10Idx = 10;
// static constexpr int comp11Idx = 11; // index of the 10th tracer
/****************************************
* Fluid phase related static parameters
......@@ -125,7 +134,7 @@ public:
if (phaseIdx == phase0Idx)
return Fluid0::isGas();
return Fluid1::isGas();
return Fluid1::isGas();
}
/*!
......@@ -181,7 +190,7 @@ public:
// let the fluids decide
if (phaseIdx == phase0Idx)
return Fluid0::isCompressible();
return Fluid1::isCompressible();
return Fluid1::isCompressible();
}
/*!
......@@ -196,7 +205,7 @@ public:
// let the fluids decide
if (phaseIdx == phase0Idx)
return Fluid0::viscosityIsConstant();
return Fluid1::viscosityIsConstant();
return Fluid1::viscosityIsConstant();
}
/*!
......@@ -212,7 +221,7 @@ public:
// let the fluids decide
if (phaseIdx == phase0Idx)
return Fluid0::isIdealFluid1();
return Fluid1::isIdealFluid1();
return Fluid1::isIdealFluid1();
}
/****************************************
......@@ -234,8 +243,27 @@ public:
//TODO: for each tracercomponent, define a unique name! Otherwise Paraview will crash.
else if (compIdx == comp2Idx)
return "NaCl1";
// else if (compIdx == comp3Idx)
// return "NaCl2";
// else if (compIdx == comp4Idx)
// return "NaCl3";
// else if (compIdx == comp5Idx)
// return "NaCl4";
// else if (compIdx == comp5Idx)
// return "NaCl5";
// else if (compIdx == comp5Idx)
// return "NaCl6";
// else if (compIdx == comp5Idx)
// return "NaCl7";
// else if (compIdx == comp5Idx)
// return "NaCl8";
// else if (compIdx == comp5Idx)
// return "NaCl9";
// else if (compIdx == comp5Idx)
// return "NaCl10";
}
/*!
* \brief Return the molar mass of a component in \f$\mathrm{[kg/mol]}\f$.
* \param compIdx index of the component
......@@ -409,15 +437,17 @@ public:
assert(0 <= phaseIdx && phaseIdx < numPhases);
assert(0 <= compIdx && compIdx < numComponents);
if (phaseIdx == phase1Idx)
if (phaseIdx == phase1Idx && compIdx == comp1Idx)
return 1e-2;
// We could calculate the real fugacity coefficient of
// the component in the fluid. Probably that's not worth
// the effort, since the fugacity coefficient of the other
// component is infinite anyway...
else if (phaseIdx == phase1Idx && compIdx != comp1Idx)
return 1.0;
else if (phaseIdx == phase0Idx && compIdx == comp0Idx)
//water likes to stay in the water phase
return 1e-10;
return 1e-2;
else if (phaseIdx == phase0Idx && compIdx == comp1Idx)
//NAPL doesn't want to be in the water phase
return fluidState.pressure(phase0Idx);
......
test/porousmediumflow/tracer/2ptracer/2pnctestproblem.hh
View file @ bde6eb7b
......@@ -32,14 +32,14 @@
#include <dumux/material/components/trichloroethene.hh>
#include <dumux/material/components/simpleh2o.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
#include "2pncimmiscible.hh"
#include <dumux/porousmediumflow/2pnc/model.hh>
#include <dumux/porousmediumflow/problem.hh>
#include "2ptestspatialparams.hh"
// #include "2ptestspatialparams_randomfield.hh"
// #include "2pnctestlocalresidual.hh"
#include "2pncimmiscible.hh"
//#include "2ptestspatialparams.hh"
#include "2ptestspatialparams_randomfield.hh"
//#include "2pnctestlocalresidual.hh"
#ifndef ENABLEINTERFACESOLVER
#define ENABLEINTERFACESOLVER 0
......@@ -120,6 +120,15 @@ class TwoPNCTestProblem : public PorousMediumFlowProblem<TypeTag>
//the tracercomponent indices
contiTracer1EqIdx = Indices::conti0EqIdx + FluidSystem::comp2Idx,
tracer1Idx = Indices::conti0EqIdx + FluidSystem::comp2Idx,
// tracer2Idx = Indices::conti0EqIdx + FluidSystem::comp3Idx,
// tracer3Idx = Indices::conti0EqIdx + FluidSystem::comp4Idx,
// tracer4Idx = Indices::conti0EqIdx + FluidSystem::comp5Idx,
// tracer5Idx = Indices::conti0EqIdx + FluidSystem::comp6Idx,
// tracer6Idx = Indices::conti0EqIdx + FluidSystem::comp7Idx,
// tracer7Idx = Indices::conti0EqIdx + FluidSystem::comp8Idx,
// tracer8Idx = Indices::conti0EqIdx + FluidSystem::comp9Idx,
// tracer9Idx = Indices::conti0EqIdx + FluidSystem::comp10Idx,
// tracer10Idx = Indices::conti0EqIdx + FluidSystem::comp11Idx,
};
// static constexpr int dimWorld = GridView::dimensionworld;
......@@ -141,7 +150,8 @@ public:
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
BoundaryTypes values;
if (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos))
// if (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos))
if (onLowerBoundary_(globalPos))
values.setAllDirichlet();
else
values.setAllNeumann();
......@@ -179,19 +189,32 @@ public:
values[saturationDNAPLIdx] = 0.0;
//the tracer component's Dirichlet BC
// VERSION 1
// if (onUpperBoundary_(globalPos))
// VERSION 2
if (onUpperBoundary_(globalPos) || onStripe1_(globalPos) || onStripe2_(globalPos) || onStripe3_(globalPos))
// VERSION 3
// if (onUpperTrajectoryPoint_(globalPos) || onLeftTrajectroyPoint_(globalPos) || onRightTrajectoryPoint_(globalPos))
if (onUpperBoundary_(globalPos))
{
if (useMoles)
if (useMoles){
values[tracer1Idx] = 1e-9;
else
// values[tracer2Idx] = 1e-9;
// values[tracer3Idx] = 1e-9;
// values[tracer4Idx] = 1e-9;
// values[tracer5Idx] = 1e-9;
// values[tracer6Idx] = 1e-9;
// values[tracer7Idx] = 1e-9;
// values[tracer8Idx] = 1e-9;
// values[tracer9Idx] = 1e-9;
// values[tracer10Idx] = 1e-9;
}
else{
values[tracer1Idx] = 1e-9*FluidSystem::molarMass(tracer1Idx)/FluidSystem::molarMass(0);
// values[tracer2Idx] = 1e-9*FluidSystem::molarMass(tracer2Idx)/FluidSystem::molarMass(0);
// values[tracer3Idx] = 1e-9*FluidSystem::molarMass(tracer3Idx)/FluidSystem::molarMass(0);
// values[tracer4Idx] = 1e-9*FluidSystem::molarMass(tracer4Idx)/FluidSystem::molarMass(0);
// values[tracer5Idx] = 1e-9*FluidSystem::molarMass(tracer5Idx)/FluidSystem::molarMass(0);
// values[tracer6Idx] = 1e-9*FluidSystem::molarMass(tracer6Idx)/FluidSystem::molarMass(0);
// values[tracer7Idx] = 1e-9*FluidSystem::molarMass(tracer7Idx)/FluidSystem::molarMass(0);
// values[tracer8Idx] = 1e-9*FluidSystem::molarMass(tracer8Idx)/FluidSystem::molarMass(0);
/*values[tracer9Idx] = 1e-9*FluidSystem::molarMass(tracer9Idx)/FluidSystem::molarMass(0);
values[tracer10Idx] = 1e-9*FluidSystem::molarMass(tracer10Idx)/FluidSystem::molarMass(0);
*/}
}
return values;
......@@ -213,8 +236,8 @@ public:
NumEqVector values(0.0);
if (onInlet_(globalPos))
{
values[contiDNAPLEqIdx] = -0.04; // kg / (m * s)
values[Indices::conti0EqIdx] = -0.04;
values[contiDNAPLEqIdx] = -0.05/FluidSystem::molarMass(contiDNAPLEqIdx); // kg / (m * s)
values[Indices::conti0EqIdx] = -0.05/FluidSystem::molarMass(Indices::conti0EqIdx);
}
// in the test with the oil wet lens, use higher injection rate
......@@ -253,19 +276,32 @@ public:
values[saturationDNAPLIdx] = 0;
//the tracer component's initial values
// VERSION 1
// if (onUpperBoundary_(globalPos))
// VERSION 2
if (onUpperBoundary_(globalPos) || onStripe1_(globalPos) || onStripe2_(globalPos) || onStripe3_(globalPos))
// VERSION 3
// if (onUpperTrajectoryPoint_(globalPos) || onLeftTrajectroyPoint_(globalPos) || onRightTrajectoryPoint_(globalPos))
if (onUpperBoundary_(globalPos))
{
if (useMoles)
if (useMoles){
values[tracer1Idx] = 1e-9;
else
// values[tracer2Idx] = 1e-9;
// values[tracer3Idx] = 1e-9;
// values[tracer4Idx] = 1e-9;
// values[tracer5Idx] = 1e-9;
// values[tracer6Idx] = 1e-9;
// values[tracer7Idx] = 1e-9;
// values[tracer8Idx] = 1e-9;
// values[tracer9Idx] = 1e-9;
// values[tracer10Idx] = 1e-9;
}
else{
values[tracer1Idx] = 1e-9*FluidSystem::molarMass(tracer1Idx)/FluidSystem::molarMass(0);
// values[tracer2Idx] = 1e-9*FluidSystem::molarMass(tracer2Idx)/FluidSystem::molarMass(0);
// values[tracer3Idx] = 1e-9*FluidSystem::molarMass(tracer3Idx)/FluidSystem::molarMass(0);
// values[tracer4Idx] = 1e-9*FluidSystem::molarMass(tracer4Idx)/FluidSystem::molarMass(0);
// values[tracer5Idx] = 1e-9*FluidSystem::molarMass(tracer5Idx)/FluidSystem::molarMass(0);
// values[tracer6Idx] = 1e-9*FluidSystem::molarMass(tracer6Idx)/FluidSystem::molarMass(0);
// values[tracer7Idx] = 1e-9*FluidSystem::molarMass(tracer7Idx)/FluidSystem::molarMass(0);
// values[tracer8Idx] = 1e-9*FluidSystem::molarMass(tracer8Idx)/FluidSystem::molarMass(0);
// values[tracer9Idx] = 1e-9*FluidSystem::molarMass(tracer9Idx)/FluidSystem::molarMass(0);
// values[tracer10Idx] = 1e-9*FluidSystem::molarMass(tracer10Idx)/FluidSystem::molarMass(0);
}
}
return values;
......@@ -322,29 +358,20 @@ private:
Scalar cellWidth_ = xMax_/xNumCells_;
Scalar width_ = xMax_ - this->fvGridGeometry().bBoxMin()[0];
// VERSION 1
bool onUpperBoundary_(const GlobalPosition &globalPos) const
{
return globalPos[1] > yMax_ - 0.1 - eps_;
}
bool onUpperBoundary_(const GlobalPosition &globalPos) const
{
return globalPos[1] > yMax_ - 0.1 - eps_;
}
// VERSION 2
bool onStripe1_(const GlobalPosition &globalPos) const
{
// std::cout << "Stripe1 is: " << (( (yMax_ /4.0 - cellHeight_*0.5) <= globalPos[1] ) &&
// ( (yMax_/4.0 + cellHeight_*0.5) > globalPos[1] ));
// std::cout << " at globalPos:" << globalPos[1] << "\n";
// std::cout << "yMax is: " << yMax_ << "\n" ;
return (
( (yMax_ /4.0 - cellHeight_*0.5) <= globalPos[1] ) &&
( (yMax_/4.0 + cellHeight_*0.5) > globalPos[1] )
);
bool onStripe1_(const GlobalPosition &globalPos) const
{
return (
( (yMax_ /4.0 - cellHeight_*0.5) <= globalPos[1] ) &&
( (yMax_/4.0 + cellHeight_*0.5) > globalPos[1] )
);
}
bool onStripe2_(const GlobalPosition &globalPos) const
{
return (
......@@ -360,44 +387,6 @@ private:
( (3.0 * yMax_/4.0 + cellHeight_*0.5) > globalPos[1] )
);
}
// VERSION 3
bool onUpperTrajectoryPoint_(const GlobalPosition &globalPos) const
{
// std::cout << "onUpperTrajectoryPoint is: "
// << ( globalPos[1] > (yMax_ - 2.0*cellHeight_) &&
// globalPos[0] > (width_/2.0 - cellWidth_) &&
// globalPos[0] < (width_/2.0 + cellWidth_) );
// std::cout << " at x- globalPos:" << globalPos[0] << "\n";
// std::cout << " at y- globalPos:" << globalPos[1] << "\n";
return (
globalPos[1] > (yMax_ - 2.0*cellHeight_) &&
globalPos[0] > (width_/2.0 - cellWidth_) &&
globalPos[0] < (width_/2.0 + cellWidth_)
);
}
bool onLeftTrajectroyPoint_(const GlobalPosition &globalPos) const
{
return (
( globalPos[1] >= (3.0 * yMax_ /4.0 - cellHeight_) ) &&
( globalPos[1] < (3.0 * yMax_/4.0 + cellHeight_) ) &&
( globalPos[0] > (xMax_ /6.0 - cellHeight_) ) &&
( globalPos[0] < (xMax_ /6.0 + cellHeight_) )
);
}
bool onRightTrajectoryPoint_(const GlobalPosition &globalPos) const
{
return (
( globalPos[1] >= (3.0 * yMax_ /4.0 - cellHeight_) ) &&
( globalPos[1] < (3.0 * yMax_/4.0 + cellHeight_) ) &&
( globalPos[0] > (5.0 * xMax_/6.0 - cellHeight_) ) &&
( globalPos[0] < (5.0 * xMax_/6.0 + cellHeight_) )
);
}
};
} // end namespace Dumux
......
test/porousmediumflow/tracer/2ptracer/2pnctestproblem_3stripes.hh 0 → 100644
View file @ bde6eb7b
// -*- 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
* \ingroup TracerTests
* \brief The properties for the incompressible test
*/
#ifndef DUMUX_INCOMPRESSIBLE_TWOP_TEST_PROBLEM_HH
#define DUMUX_INCOMPRESSIBLE_TWOP_TEST_PROBLEM_HH
#include <dune/grid/yaspgrid.hh>
#include <dumux/discretization/box/properties.hh>
#include <dumux/discretization/cellcentered/tpfa/properties.hh>
#include <dumux/material/components/trichloroethene.hh>
#include <dumux/material/components/simpleh2o.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
#include <dumux/porousmediumflow/2pnc/model.hh>
#include <dumux/porousmediumflow/problem.hh>
#include "2pncimmiscible.hh"
#include "2ptestspatialparams.hh"
//#include "2ptestspatialparams_randomfield.hh"
//#include "2pnctestlocalresidual.hh"
#ifndef ENABLEINTERFACESOLVER
#define ENABLEINTERFACESOLVER 0
#endif
namespace Dumux
{
/*!
* \ingroup TracerTests
* \brief The properties for the incompressible 2p test
*/
// forward declarations
template<class TypeTag>
class TwoPNCTestProblem;
namespace Properties
{
NEW_TYPE_TAG(TwoPNCTestProblem, INHERITS_FROM(TwoPNC));
NEW_TYPE_TAG(TwoPNCTestProblemTpfa, INHERITS_FROM(CCTpfaModel, TwoPNCTestProblem, SpatialParams));
NEW_TYPE_TAG(TwoPNCTestProblemBox, INHERITS_FROM(BoxModel, TwoPNCTestProblem, SpatialParams));
// Set the grid type
SET_TYPE_PROP(TwoPNCTestProblem, Grid, Dune::YaspGrid<2>);
// Set the problem type
SET_TYPE_PROP(TwoPNCTestProblem, Problem, TwoPNCTestProblem<TypeTag>);
// the fluid system
SET_PROP(TwoPNCTestProblem, FluidSystem)
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using WettingPhase = FluidSystems::OnePLiquid<Scalar, Components::SimpleH2O<Scalar> >;
using NonwettingPhase = FluidSystems::OnePLiquid<Scalar, Components::Trichloroethene<Scalar> >;
using type = FluidSystems::TwoPNCImmiscible<Scalar, WettingPhase, NonwettingPhase>;
};
// Enable caching
SET_BOOL_PROP(TwoPNCTestProblem, EnableGridVolumeVariablesCache, false);
SET_BOOL_PROP(TwoPNCTestProblem, EnableGridFluxVariablesCache, false);
SET_BOOL_PROP(TwoPNCTestProblem, EnableFVGridGeometryCache, false);
// Maybe enable the box-interface solver
SET_BOOL_PROP(TwoPNCTestProblem, EnableBoxInterfaceSolver, ENABLEINTERFACESOLVER);
} // end namespace Properties
/*!
* \ingroup TracerTests
* \brief The incompressible 2p test problem.
*/
template<class TypeTag>
class TwoPNCTestProblem : public PorousMediumFlowProblem<TypeTag>
{
using ParentType = PorousMediumFlowProblem<TypeTag>;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
using GlobalPosition = Dune::FieldVector<Scalar, GridView::dimensionworld>;
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
//! property that defines whether mole or mass fractions are used
static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
enum {
pressureH2OIdx = Indices::pressureIdx,
saturationDNAPLIdx = Indices::switchIdx,
contiDNAPLEqIdx = Indices::conti0EqIdx + FluidSystem::comp1Idx,
waterPhaseIdx = FluidSystem::phase0Idx,
dnaplPhaseIdx = FluidSystem::phase1Idx,
//the tracercomponent indices
contiTracer1EqIdx = Indices::conti0EqIdx + FluidSystem::comp2Idx,
tracer1Idx = Indices::conti0EqIdx + FluidSystem::comp2Idx,
// tracer2Idx = Indices::conti0EqIdx + FluidSystem::comp3Idx,
// tracer3Idx = Indices::conti0EqIdx + FluidSystem::comp4Idx,
// tracer4Idx = Indices::conti0EqIdx + FluidSystem::comp5Idx,
// tracer5Idx = Indices::conti0EqIdx + FluidSystem::comp6Idx,
// tracer6Idx = Indices::conti0EqIdx + FluidSystem::comp7Idx,
// tracer7Idx = Indices::conti0EqIdx + FluidSystem::comp8Idx,
// tracer8Idx = Indices::conti0EqIdx + FluidSystem::comp9Idx,
// tracer9Idx = Indices::conti0EqIdx + FluidSystem::comp10Idx,
// tracer10Idx = Indices::conti0EqIdx + FluidSystem::comp11Idx,
};
// static constexpr int dimWorld = GridView::dimensionworld;
public:
TwoPNCTestProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
{
Dune::FMatrixPrecision<>::set_singular_limit(1e-35);
}
/*!
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary segment
*
* \param values Stores the value of the boundary type
* \param globalPos The global position
*/
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
BoundaryTypes values;
// if (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos))
if (onLowerBoundary_(globalPos))
values.setAllDirichlet();
else
values.setAllNeumann();
return values;
}
/*!
* \brief Evaluates the boundary conditions for a Dirichlet
* boundary segment
*
* \param values Stores the Dirichlet values for the conservation equations in
* \f$ [ \textnormal{unit of primary variable} ] \f$
* \param globalPos The global position
*/
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables values;
values.setState(Indices::firstPhaseOnly);
typename GET_PROP_TYPE(TypeTag, FluidState) fluidState;
fluidState.setTemperature(temperature());
fluidState.setPressure(waterPhaseIdx, /*pressure=*/1e5);
fluidState.setPressure(dnaplPhaseIdx, /*pressure=*/1e5);
Scalar densityW = FluidSystem::density(fluidState, waterPhaseIdx);
// Scalar height = this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1];
Scalar depth = this->fvGridGeometry().bBoxMax()[1] - globalPos[1];
// Scalar alpha = 1 + 1.5/height;
// Scalar width = this->fvGridGeometry().bBoxMax()[0] - this->fvGridGeometry().bBoxMin()[0];
Scalar factor = 1;
// hydrostatic pressure scaled by alpha
values[pressureH2OIdx] = 1e5 - factor*densityW*this->gravity()[1]*depth;
values[saturationDNAPLIdx] = 0.0;
//the tracer component's Dirichlet BC
// VERSION "3stripes"
if (onUpperBoundary_(globalPos) || onStripe1_(globalPos) || onStripe2_(globalPos) || onStripe3_(globalPos))
{
if (useMoles){
values[tracer1Idx] = 1e-9;
// values[tracer2Idx] = 1e-9;
// values[tracer3Idx] = 1e-9;
// values[tracer4Idx] = 1e-9;
// values[tracer5Idx] = 1e-9;
// values[tracer6Idx] = 1e-9;
// values[tracer7Idx] = 1e-9;
// values[tracer8Idx] = 1e-9;
// values[tracer9Idx] = 1e-9;
// values[tracer10Idx] = 1e-9;
}
else{
values[tracer1Idx] = 1e-9*FluidSystem::molarMass(tracer1Idx)/FluidSystem::molarMass(0);
// values[tracer2Idx] = 1e-9*FluidSystem::molarMass(tracer2Idx)/FluidSystem::molarMass(0);
// values[tracer3Idx] = 1e-9*FluidSystem::molarMass(tracer3Idx)/FluidSystem::molarMass(0);
// values[tracer4Idx] = 1e-9*FluidSystem::molarMass(tracer4Idx)/FluidSystem::molarMass(0);
// values[tracer5Idx] = 1e-9*FluidSystem::molarMass(tracer5Idx)/FluidSystem::molarMass(0);
// values[tracer6Idx] = 1e-9*FluidSystem::molarMass(tracer6Idx)/FluidSystem::molarMass(0);
// values[tracer7Idx] = 1e-9*FluidSystem::molarMass(tracer7Idx)/FluidSystem::molarMass(0);
// values[tracer8Idx] = 1e-9*FluidSystem::molarMass(tracer8Idx)/FluidSystem::molarMass(0);
/*values[tracer9Idx] = 1e-9*FluidSystem::molarMass(tracer9Idx)/FluidSystem::molarMass(0);
values[tracer10Idx] = 1e-9*FluidSystem::molarMass(tracer10Idx)/FluidSystem::molarMass(0);
*/}
}
return values;
}
/*!
* \brief Evaluate the boundary conditions for a neumann
* boundary segment.
*
* \param values Stores the Neumann values for the conservation equations in
* \f$ [ \textnormal{unit of conserved quantity} / (m^(dim-1) \cdot s )] \f$
* \param globalPos The position of the integration point of the boundary segment.
*
* For this method, the \a values parameter stores the mass flux
* in normal direction of each phase. Negative values mean influx.
*/
NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
{
NumEqVector values(0.0);
if (onInlet_(globalPos))
{
values[contiDNAPLEqIdx] = -0.05/FluidSystem::molarMass(contiDNAPLEqIdx); // kg / (m * s)
values[Indices::conti0EqIdx] = -0.05/FluidSystem::molarMass(Indices::conti0EqIdx);
}
// in the test with the oil wet lens, use higher injection rate
if (this->spatialParams().lensIsOilWet())
values[contiDNAPLEqIdx] *= 10;
//no tracer is injected in the tracerproblem, so we do not need to set a Neumann BC for it different from 0!
return values;
}
/*!
* \brief Evaluates the initial values for a control volume
*
* \param values Stores the initial values for the conservation equations in
* \f$ [ \textnormal{unit of primary variables} ] \f$
* \param globalPos The global position
*/
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables values;