Commit 1687727b authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

[test][tracer] restructure and improve tracer conservation test

Outsource properties to a separate header. Move function to check
conservation from problem to main. Improve variable names. Switch
to one space dimension, adapt simulation time interval. Throw an
exception if tracer is not conserved.
parent 07862579
Pipeline #7152 waiting for manual action with stages
......@@ -29,6 +29,7 @@
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh>
#include <dune/common/float_cmp.hh>
#include <dune/common/exceptions.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
......@@ -43,8 +44,7 @@
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager.hh>
#include "problem_2p.hh"
#include "problem_tracer.hh"
#include "properties.hh"
template <class MoleFractionVector, class SaturationVector>
void equilibrateTracer(MoleFractionVector& moleFracVec,
......@@ -59,6 +59,70 @@ void equilibrateTracer(MoleFractionVector& moleFracVec,
}
}
template <class Scalar, class Problem, class SolutionVector, class GridVariables>
Scalar checkConservation(const Problem& problem,
const SolutionVector& sol,
const GridVariables& gridVars,
Scalar dt,
Scalar accumulatedInflow)
{
Scalar storedAmount = 0.0;
const auto& gridGeometry = problem.gridGeometry();
for (const auto& element : elements(gridGeometry.gridView()))
{
auto fvGeometry = localView(gridGeometry);
fvGeometry.bind(element);
auto elemVolVars = localView(gridVars.curGridVolVars());
elemVolVars.bind(element, fvGeometry, sol);
for (auto&& scv : scvs(fvGeometry))
{
const auto& volVars = elemVolVars[scv];
if(problem.useMoles())
{
storedAmount += volVars.moleFraction(/*phaseIdx*/0, /*compIdx*/0) * volVars.molarDensity(/*phaseIdx*/0)
* scv.volume() * volVars.saturation(/*phaseIdx*/0) * volVars.porosity() * volVars.extrusionFactor();
}
else
storedAmount += volVars.massFraction(/*phaseIdx*/0, /*compIdx*/0) * volVars.density(/*phaseIdx*/0)
* scv.volume() * volVars.saturation(/*phaseIdx*/0) * volVars.porosity() * volVars.extrusionFactor();
}
for (auto&& scvf : scvfs(fvGeometry))
{
if (scvf.boundary())
{
const auto& volVars = elemVolVars[scvf.insideScvIdx()];
if (problem.useDirichlet())
{
const auto moleFraction = problem.dirichletAtPos(scvf.ipGlobal());
accumulatedInflow -= problem.spatialParams().volumeFlux(element, fvGeometry, elemVolVars, scvf)
* moleFraction
* volVars.molarDensity(/*phaseIdx*/0)
* scvf.area() * volVars.extrusionFactor()
* dt;
}
else
accumulatedInflow -= problem.neumannAtPos(scvf.ipGlobal())[0]
* scvf.area() * volVars.extrusionFactor()
* dt;
}
}
}
const std::string unit = problem.useMoles() ? " moles " : " kg ";
std::cout << "\033[1;31m" << storedAmount << unit << "tracer are inside the domain. \033[0m" << '\n';
std::cout << "\033[1;31m" << accumulatedInflow << unit << "should have flown into the domain. \033[0m" << '\n';
std::cout << "\033[1;31m" << accumulatedInflow - storedAmount << unit << "are missing. \033[0m" << '\n';
if (std::abs(accumulatedInflow - storedAmount) > std::numeric_limits<Scalar>::epsilon())
DUNE_THROW(Dune::InvalidStateException, "tracer is not conserved");
return accumulatedInflow;
}
int main(int argc, char** argv)
{
using namespace Dumux;
......@@ -220,7 +284,7 @@ int main(int argc, char** argv)
////////////////////////////////////////////////////////////
// time loop
double inflowMass = 0.0;
Scalar accumulatedInflow = 0.0;
timeLoop->start(); do
{
// solve the non-linear system with time step control
......@@ -228,6 +292,7 @@ int main(int argc, char** argv)
// make the new solution the old solution
pOld = p;
oldSaturation_ = saturation_;
twoPGridVariables->advanceTimeStep();
// loop over elements to compute fluxes, saturations, densities for tracer
......@@ -277,7 +342,6 @@ int main(int argc, char** argv)
}
}
oldSaturation_ = saturation_;
for (const auto& scv : scvs(fvGeometry))
{
const auto& volVars = elemVolVars[scv];
......@@ -325,7 +389,8 @@ int main(int argc, char** argv)
<< updateTimer.elapsed() << "(" << 100*updateTimer.elapsed()/elapsedTot << "%)"
<< std::endl;
inflowMass = tracerProblem->printTracerMass(x, *tracerGridVariables, timeLoop->timeStepSize(), inflowMass);
accumulatedInflow = checkConservation(*tracerProblem, x, *tracerGridVariables,
timeLoop->timeStepSize(), accumulatedInflow);
// make the new solution the old solution
xOld = x;
......
[TimeLoop]
DtInitial = 1e4 # [s]
TEnd = 1e5 # [s]
DtInitial = 1e5 # [s]
TEnd = 3e6 # [s]
[Grid]
UpperRight = 1 10
Cells = 1 10
UpperRight = 10
Cells = 10
[SpatialParams]
LinearPcEntry = 0.0
......
......@@ -18,79 +18,20 @@
*****************************************************************************/
/*!
* \ingroup TwoPTests
* \brief The properties for the flow problem of a radioactive 2p test.
* \brief initial and boundary conditions for the flow problem of the tracer conservation test
*/
#ifndef DUMUX_TWOP_FLOW_TEST_PROBLEM_HH
#define DUMUX_TWOP_FLOW_TEST_PROBLEM_HH
#include <iostream>
#include <dune/grid/yaspgrid.hh>
#include <dumux/common/boundarytypes.hh>
#include <dumux/common/numeqvector.hh>
#include <dumux/discretization/box.hh>
#include <dumux/discretization/cctpfa.hh>
#include <dumux/material/components/simpleh2o.hh>
#include <dumux/material/components/constant.hh>
#include <dumux/material/fluidsystems/h2oair.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
#include <dumux/material/fluidsystems/1pgas.hh>
#include <dumux/material/fluidsystems/2pimmiscible.hh>
#include <dumux/porousmediumflow/problem.hh>
#include <dumux/porousmediumflow/2p/model.hh>
#include <dumux/porousmediumflow/2p/incompressiblelocalresidual.hh>
#include "spatialparams_2p.hh"
namespace Dumux {
// forward declarations
template<class TypeTag> class TwoPFlowTestProblem;
namespace Properties {
// Create new type tags
namespace TTag {
struct TwoPFlow { using InheritsFrom = std::tuple<TwoP>; };
struct TwoPFlowTpfa { using InheritsFrom = std::tuple<TwoPFlow, CCTpfaModel>; };
} // end namespace TTag
// Set the grid type
template<class TypeTag>
struct Grid<TypeTag, TTag::TwoPFlow> { using type = Dune::YaspGrid<2>; };
// Set the problem type
template<class TypeTag>
struct Problem<TypeTag, TTag::TwoPFlow> { using type = TwoPFlowTestProblem<TypeTag>; };
// Set the fluid system
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::TwoPFlow>
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using WettingPhase = FluidSystems::OnePLiquid<Scalar, Components::Constant<0, Scalar> >;
using NonwettingPhase = FluidSystems::OnePGas<Scalar, Components::Constant<1, Scalar> >;
using type = FluidSystems::TwoPImmiscible<Scalar, WettingPhase, NonwettingPhase>;
};
// Set the spatial parameters
template<class TypeTag>
struct SpatialParams<TypeTag, TTag::TwoPFlow>
{
private:
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
public:
using type = TwoPFlowTestSpatialParams<GridGeometry, Scalar>;
};
} // end namespace Properties
/*!
* \ingroup TwoPTests
* \brief The properties for the flow problem of a radioactive 2p test.
* \brief initial and boundary conditions for the flow problem of the tracer conservation test
*/
template<class TypeTag>
class TwoPFlowTestProblem : public PorousMediumFlowProblem<TypeTag>
......@@ -120,7 +61,7 @@ public:
{
BoundaryTypes values;
if (onLowerBoundary_(globalPos))
if (onRightBoundary_(globalPos))
values.setAllDirichlet();
else
values.setAllNeumann();
......@@ -132,7 +73,7 @@ public:
{
NumEqVector values(0.0);
if (onUpperBoundary_(globalPos))
if (onLeftBoundary_(globalPos))
{
values[conti0EqIdx] = -injectionMass_;
}
......@@ -161,19 +102,17 @@ public:
{ return 293.15; }
private:
bool onLowerBoundary_(const GlobalPosition &globalPos) const
bool onLeftBoundary_(const GlobalPosition &globalPos) const
{
return globalPos[1] < this->gridGeometry().bBoxMin()[1] + eps_;
return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_;
}
bool onUpperBoundary_(const GlobalPosition &globalPos) const
bool onRightBoundary_(const GlobalPosition &globalPos) const
{
return globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_;
return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_;
}
static constexpr Scalar eps_ = 1e-6;
Scalar injectionMass_;
};
......
......@@ -21,103 +21,16 @@
* \ingroup TracerTests
* \brief problem for a tracer conservation test
*/
#ifndef DUMUX_TEST_PROBLEM_TRACER_HH
#define DUMUX_TEST_PROBLEM_TRACER_HH
#ifndef DUMUX_TEST_TRACER_CONSERVATION_PROBLEM_HH
#define DUMUX_TEST_TRACER_CONSERVATION_PROBLEM_HH
#include <dune/grid/yaspgrid.hh>
#include <iostream>
#include <dumux/common/boundarytypes.hh>
#include <dumux/common/numeqvector.hh>
#include <dumux/discretization/cctpfa.hh>
#include <dumux/porousmediumflow/tracer/model.hh>
#include <dumux/porousmediumflow/problem.hh>
#include <dumux/material/fluidsystems/base.hh>
#include "spatialparams_tracer.hh"
namespace Dumux {
/*!
* \ingroup TracerTests
* \brief problem for a tracer conservation test
*/
template <class TypeTag>
class TracerConservationTestProblem;
namespace Properties {
//Create new type tags
namespace TTag {
struct TracerConservationTest { using InheritsFrom = std::tuple<Tracer>; };
struct TracerConservationTestTpfa { using InheritsFrom = std::tuple<TracerConservationTest, CCTpfaModel>; };
} // end namespace TTag
// Set the grid type
template<class TypeTag>
struct Grid<TypeTag, TTag::TracerConservationTest> { using type = Dune::YaspGrid<2>; };
// Set the problem property
template<class TypeTag>
struct Problem<TypeTag, TTag::TracerConservationTest> { using type = TracerConservationTestProblem<TypeTag>; };
// Set the spatial parameters
template<class TypeTag>
struct SpatialParams<TypeTag, TTag::TracerConservationTest>
{
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using type = TracerConservationTestSpatialParams<GridGeometry, Scalar>;
};
//! A simple fluid system with one tracer component
template<class TypeTag>
class TracerFluidSystem : public FluidSystems::Base<GetPropType<TypeTag, Properties::Scalar>,
TracerFluidSystem<TypeTag>>
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Problem = GetPropType<TypeTag, Properties::Problem>;
using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
using Element = typename GridView::template Codim<0>::Entity;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
public:
//! If the fluid system only contains tracer components
static constexpr bool isTracerFluidSystem()
{ return true; }
//! No component is the main component
static constexpr int getMainComponent(int phaseIdx)
{ return -1; }
//! The number of components
static constexpr int numComponents = 1;
//! Human readable component name (for vtk output)
static std::string componentName(int compIdx)
{ return "Cs137"; }
//! Molar mass in kg/mol
static Scalar molarMass(int compIdx)
{ return 0.1; }
//! Binary diffusion coefficient
//! (might depend on spatial parameters like pressure / temperature)
static Scalar binaryDiffusionCoefficient(unsigned int compIdx,
const Problem& problem,
const Element& element,
const SubControlVolume& scv)
{
static const Scalar D = getParam<Scalar>("Problem.BinaryDiffusionCoefficient");
return D;
}
};
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::TracerConservationTest> { using type = TracerFluidSystem<TypeTag>; };
} // end namespace Properties
template <class TypeTag>
class TracerConservationTestProblem : public PorousMediumFlowProblem<TypeTag>
......@@ -132,14 +45,12 @@ class TracerConservationTestProblem : public PorousMediumFlowProblem<TypeTag>
using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
static constexpr bool useMoles = getPropValue<TypeTag, Properties::UseMoles>();
public:
TracerConservationTestProblem(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry)
{
// stating in the console whether mole or mass fractions are used
if(useMoles)
if(useMoles_)
std::cout<<"problem uses mole fractions" << '\n';
else
std::cout<<"problem uses mass fractions" << '\n';
......@@ -151,7 +62,7 @@ public:
{
BoundaryTypes values;
if (useDirichlet_ && onUpperBoundary_(globalPos))
if (useDirichlet_ && onLeftBoundary_(globalPos))
values.setAllDirichlet();
else
values.setAllNeumann();
......@@ -163,9 +74,9 @@ public:
{
NumEqVector values(0.0);
if (onUpperBoundary_(globalPos))
if (onLeftBoundary_(globalPos))
{
if (useMoles)
if (useMoles_)
values = -1e-10;
else
values = -1e-10*FluidSystem::molarMass(0)/this->spatialParams().fluidMolarMass(globalPos);
......@@ -178,7 +89,7 @@ public:
{
PrimaryVariables values;
if (onUpperBoundary_(globalPos))
if (onLeftBoundary_(globalPos))
values = 1e-9;
else
values = initialAtPos(globalPos);
......@@ -193,83 +104,19 @@ public:
return initialValues;
}
//! Helper function for mass balance evaluations
template <class TracerSolutionVector, class GridVariables>
Scalar printTracerMass(const TracerSolutionVector& sol,
const GridVariables& gridVars,
Scalar dt,
Scalar inflowMass)
{
// compute the mass in the entire domain to make sure the tracer is conserved
Scalar tracerMass = 0.0;
const auto& gg = this->gridGeometry();
for (const auto& element : elements(gg.gridView()))
{
auto fvGeometry = localView(gg);
fvGeometry.bind(element);
auto elemVolVars = localView(gridVars.curGridVolVars());
elemVolVars.bind(element, fvGeometry, sol);
auto elemFluxVarsCache = localView(gridVars.gridFluxVarsCache());
elemFluxVarsCache.bind(element, fvGeometry, elemVolVars);
for (auto&& scv : scvs(fvGeometry))
{
const auto& volVars = elemVolVars[scv];
if(useMoles)
{
tracerMass += volVars.moleFraction(/*phaseIdx*/0, /*compIdx*/0) * volVars.molarDensity(/*phaseIdx*/0)
* scv.volume() * volVars.saturation(/*phaseIdx*/0) * volVars.porosity() * volVars.extrusionFactor();
}
else
tracerMass += volVars.massFraction(/*phaseIdx*/0, /*compIdx*/0) * volVars.density(/*phaseIdx*/0)
* scv.volume() * volVars.saturation(/*phaseIdx*/0) * volVars.porosity() * volVars.extrusionFactor();
}
for (auto&& scvf : scvfs(fvGeometry))
{
if (scvf.boundary())
{
const auto& volVars = elemVolVars[scvf.insideScvIdx()];
if (useDirichlet_)
{
const auto moleFraction = dirichletAtPos(scvf.ipGlobal());
inflowMass -= this->spatialParams().volumeFlux(element, fvGeometry, elemVolVars, scvf)
* moleFraction
* volVars.molarDensity(/*phaseIdx*/0)
* scvf.area() * volVars.extrusionFactor()
* dt;
}
else
inflowMass -= neumannAtPos(scvf.ipGlobal())[0]
* scvf.area() * volVars.extrusionFactor()
* dt;
}
}
}
std::cout << "\033[1;31m" << tracerMass << " moles tracer are inside the domain. \033[0m" << '\n';
std::cout << "\033[1;31m" << inflowMass << " moles should have flown into the domain. \033[0m" << '\n';
std::cout << "\033[1;31m" << inflowMass - tracerMass << " moles are missing. \033[0m" << '\n';
return inflowMass;
}
bool useMoles() const
{ return useMoles_; }
bool useDirichlet() const
{ return useDirichlet_; }
private:
static constexpr Scalar eps_ = 1e-6;
bool onLowerBoundary_(const GlobalPosition &globalPos) const
{
return globalPos[1] < this->gridGeometry().bBoxMin()[1] + eps_;
}
bool onUpperBoundary_(const GlobalPosition &globalPos) const
bool onLeftBoundary_(const GlobalPosition &globalPos) const
{
return globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_;
return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_;
}
static constexpr Scalar eps_ = 1e-6;
static constexpr bool useMoles_ = getPropValue<TypeTag, Properties::UseMoles>();
bool useDirichlet_;
};
......
// -*- 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 3 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 properties for the tracer conservation test
*/
#ifndef DUMUX_TEST_TRACER_CONSERVATION_PROPERTIES_HH
#define DUMUX_TEST_TRACER_CONSERVATION_PROPERTIES_HH
#include <dune/grid/yaspgrid.hh>
#include <dumux/discretization/cctpfa.hh>
#include <dumux/material/components/constant.hh>
#include <dumux/material/fluidsystems/base.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
#include <dumux/material/fluidsystems/1pgas.hh>
#include <dumux/material/fluidsystems/2pimmiscible.hh>
#include <dumux/porousmediumflow/tracer/model.hh>
#include <dumux/porousmediumflow/2p/model.hh>
#include "spatialparams_tracer.hh"
#include "problem_tracer.hh"
#include "spatialparams_2p.hh"
#include "problem_2p.hh"
namespace Dumux {
namespace Properties {
////////////////////////////////////////
// tracer properties
////////////////////////////////////////
//Create new type tags
namespace TTag {
struct TracerConservationTest { using InheritsFrom = std::tuple<Tracer>; };
struct TracerConservationTestTpfa { using InheritsFrom = std::tuple<TracerConservationTest, CCTpfaModel>; };
} // end namespace TTag
// Set the grid type
template<class TypeTag>
struct Grid<TypeTag, TTag::TracerConservationTest> { using type = Dune::YaspGrid<1>; };
// Set the problem property
template<class TypeTag>
struct Problem<TypeTag, TTag::TracerConservationTest> { using type = TracerConservationTestProblem<TypeTag>; };
// Set the spatial parameters
template<class TypeTag>
struct SpatialParams<TypeTag, TTag::TracerConservationTest>
{
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using type = TracerConservationTestSpatialParams<GridGeometry, Scalar>;
};
//! A simple fluid system with one tracer component
template<class TypeTag>
class TracerFluidSystem : public FluidSystems::Base<GetPropType<TypeTag, Properties::Scalar>,
TracerFluidSystem<TypeTag>>
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Problem = GetPropType<TypeTag, Properties::Problem>;
using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
using Element = typename GridView::template Codim<0>::Entity;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
public: