Commit f0f325d4 authored by Timo Koch's avatar Timo Koch
Browse files

[test] Add tracer multiphase test

parent 37530a96
add_subdirectory(2ptracer)
add_subdirectory(constvel)
add_subdirectory(multiphase)
dune_symlink_to_source_files(FILES "params.input")
dune_add_test(NAME test_tracer_multiphase_tpfa
LABELS porousmediumflow tracer
SOURCES main.cc
COMPILE_DEFINITIONS TYPETAG=TracerTestTpfa
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_tracer_multiphase_tpfa-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_tracer_multiphase_tpfa-00010.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_tracer_multiphase_tpfa params.input -Problem.Name test_tracer_multiphase_tpfa")
dune_add_test(NAME test_tracer_multiphase_mpfa
LABELS porousmediumflow tracer
SOURCES main.cc
COMPILE_DEFINITIONS TYPETAG=TracerTestMpfa
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_tracer_multiphase_tpfa-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_tracer_multiphase_mpfa-00010.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_tracer_multiphase_mpfa params.input -Problem.Name test_tracer_multiphase_mpfa")
// -*- 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 The test fluidsystem
*/
#ifndef DUMUX_TRACER_MULTIPHASE_TEST_FLUIDSYSTEM_HH
#define DUMUX_TRACER_MULTIPHASE_TEST_FLUIDSYSTEM_HH
#include <dumux/material/fluidsystems/base.hh>
namespace Dumux::FluidSystems {
//! A simple fluid system with two independent tracer components
template<class Scalar>
class TracerTest : public Base<Scalar, TracerTest<Scalar>>
{
public:
static constexpr bool isTracerFluidSystem()
{ return true; }
//! The number of components
static constexpr int numComponents = 2;
//! Human readable component name (index compIdx) (for vtk output)
static std::string componentName(int compIdx)
{ return "tracer_" + std::to_string(compIdx); }
//! Human readable phase name (index phaseIdx) (for velocity vtk output)
static std::string phaseName(int phaseIdx = 0)
{ return "aq"; }
//! Molar mass in kg/mol of the component with index compIdx
static Scalar molarMass(unsigned int compIdx)
{ return 0.300; }
//! Binary diffusion coefficient
//! (might depend on spatial parameters like pressure / temperature)
template<class Problem, class Element, class SubControlVolume>
static Scalar binaryDiffusionCoefficient(unsigned int compIdx,
const Problem& problem,
const Element& element,
const SubControlVolume& scv)
{
static const Scalar D = getParam<Scalar>("Problem.D");
static const Scalar D2 = getParam<Scalar>("Problem.D2");
if (compIdx == 0)
return D;
else
return D2;
}
/*!
* \copydoc Dumux::FluidSystems::Base::isCompressible
*/
static constexpr bool isCompressible(int phaseIdx)
{ return false; }
/*!
* \copydoc Dumux::FluidSystems::Base::viscosityIsConstant
*/
static constexpr bool viscosityIsConstant(int phaseIdx)
{ return true; }
};
} // end namespace Dumux::FluidSystems
#endif
// -*- 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 Test for the tracer CC model.
*/
#include <config.h>
#include <ctime>
#include <iostream>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh>
#include <dune/grid/io/file/dgfparser/dgfexception.hh>
#include <dune/grid/io/file/vtk/vtksequencewriter.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/common/defaultusagemessage.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/linear/pdesolver.hh>
#include <dumux/assembly/fvassembler.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager.hh>
#include "properties.hh"
int main(int argc, char** argv) try
{
using namespace Dumux;
//! define the type tag for this problem
using TypeTag = Properties::TTag::TYPETAG;
//! initialize MPI, finalize is done automatically on exit
const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
//! print dumux start message
if (mpiHelper.rank() == 0)
DumuxMessage::print(/*firstCall=*/true);
// parse the command line arguments and input file
Parameters::init(argc, argv);
// try to create a grid (from the given grid file or the input file)
GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
gridManager.init();
////////////////////////////////////////////////////////////
// setup instationary non-linear problem on this grid
////////////////////////////////////////////////////////////
//! we compute on the leaf grid view
const auto& leafGridView = gridManager.grid().leafGridView();
//! create the finite volume grid geometry
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
gridGeometry->update();
//! the problem (initial and boundary conditions)
using Problem = GetPropType<TypeTag, Properties::Problem>;
auto problem = std::make_shared<Problem>(gridGeometry);
//! the solution vector
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
SolutionVector x(gridGeometry->numDofs());
problem->applyInitialSolution(x);
auto xOld = x;
//! the grid variables
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
gridVariables->init(x);
//! intialize the vtk output module
VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
using VelocityOutput = GetPropType<TypeTag, Properties::VelocityOutput>;
vtkWriter.addVelocityOutput(std::make_shared<VelocityOutput>(*gridVariables));
using IOFields = GetPropType<TypeTag, Properties::IOFields>;
IOFields::initOutputModule(vtkWriter); // Add model specific output fields
vtkWriter.addVolumeVariable([](const auto& v){ return v.saturation(); }, "S_aq");
vtkWriter.write(0.0);
//! get some time loop parameters
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
auto dt = getParam<Scalar>("TimeLoop.Dt");
//! instantiate time loop
auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0, dt, tEnd);
//! the assembler with time loop for instationary problem
using Assembler = FVAssembler<TypeTag, DiffMethod::numeric, /*implicit=*/false>;
auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop);
assembler->setPreviousSolution(xOld);
//! the linear solver
using LinearSolver = ExplicitDiagonalSolver;
auto linearSolver = std::make_shared<LinearSolver>();
//! the system solver
LinearPDESolver<Assembler, LinearSolver> solver(assembler, linearSolver);
/////////////////////////////////////////////////////////////////////////////////////////////////
// run instationary non-linear simulation
/////////////////////////////////////////////////////////////////////////////////////////////////
//! set some check points for the time loop
timeLoop->setPeriodicCheckPoint(tEnd/10.0);
//! start the time loop
timeLoop->start();
while (!timeLoop->finished())
{
// assemble & solve
solver.solve(x);
// make the new solution the old solution
xOld = x;
gridVariables->advanceTimeStep();
// advance to the time loop to the next step
timeLoop->advanceTimeStep();
// report statistics of this time step
timeLoop->reportTimeStep();
// write vtk output on check points
if (timeLoop->isCheckPoint() || timeLoop->finished())
vtkWriter.write(timeLoop->time());
// set new dt TODO implement CFL-criterion
timeLoop->setTimeStepSize(dt);
}
timeLoop->finalize(leafGridView.comm());
////////////////////////////////////////////////////////////
// finalize, print dumux message to say goodbye
////////////////////////////////////////////////////////////
//! print dumux end message
if (mpiHelper.rank() == 0)
Parameters::print();
return 0;
}
catch (const Dumux::ParameterException &e)
{
std::cerr << std::endl << e << " ---> Abort!" << std::endl;
return 1;
}
catch (const Dune::DGFException & e)
{
std::cerr << "DGF exception thrown (" << e <<
"). Most likely, the DGF file name is wrong "
"or the DGF file is corrupted, "
"e.g. missing hash at end of file or wrong number (dimensions) of entries."
<< " ---> Abort!" << std::endl;
return 2;
}
catch (const Dune::Exception &e)
{
std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
return 3;
}
catch (...)
{
std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
return 4;
}
[TimeLoop]
Dt = 10000 # [s]
TEnd = 1e7 # [s]
[Grid]
UpperRight = 1 1
Cells = 50 50
[Problem]
Name = tracer
D = 0
D2 = 1e-8
FrontVelocity = 1e-7
Porosity = 0.5
[Vtk]
AddVelocity = true
[Assembly]
NumericDifference.BaseEpsilon = 1
// -*- 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 Definition of a problem for the tracer problem:
* A rotating velocity field mixes a tracer band in a porous groundwater reservoir.
*/
#ifndef DUMUX_TRACER_MULTIPHASE_TEST_PROBLEM_HH
#define DUMUX_TRACER_MULTIPHASE_TEST_PROBLEM_HH
#include <dumux/porousmediumflow/problem.hh>
namespace Dumux {
/*!
* \ingroup TracerTests
* \brief A problem for the multiphase tracer model
*/
template <class TypeTag>
class TracerTest : public PorousMediumFlowProblem<TypeTag>
{
using ParentType = PorousMediumFlowProblem<TypeTag>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using GridView = typename GridGeometry::GridView;
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
//! property that defines whether mole or mass fractions are used
static constexpr bool useMoles = getPropValue<TypeTag, Properties::UseMoles>();
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
public:
TracerTest(std::shared_ptr<const GridGeometry> gg)
: ParentType(gg)
{
// stating in the console whether mole or mass fractions are used
if(useMoles)
std::cout<<"problem uses mole fractions" << '\n';
else
std::cout<<"problem uses mass fractions" << '\n';
}
/*!
* \brief The boundary types
*/
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
BoundaryTypes values;
values.setAllNeumann();
// tracer is only present in water
if (this->spatialParams().isWater(globalPos) && globalPos[0] < eps_)
values.setAllDirichlet();
return values;
}
/*!
* \brief The Dirichlet values
*/
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
if (globalPos[0] < eps_)
return PrimaryVariables({1e-8, 1e-8});
else
return PrimaryVariables({0.0, 0.0});
}
/*!
* \brief Evaluate the boundary conditions for a neumann
* boundary segment.
*
* This is the method for the case where the Neumann condition is
* potentially solution dependent
*
* \param element The finite element
* \param fvGeometry The finite-volume geometry
* \param elemVolVars All volume variables for the element
* \param scvf The sub control volume face
*
* Negative values mean influx.
* E.g. for the mass balance that would the mass flux in \f$ [ kg / (m^2 \cdot s)] \f$.
*/
template<class ElemFluxVarsCache>
NumEqVector neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElemFluxVarsCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
NumEqVector values(0.0);
const auto& globalPos = scvf.ipGlobal();
if (globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_)
{
const auto& volVars = elemVolVars[scvf.insideScvIdx()];
values = NumEqVector({volVars.massFraction(0,0), volVars.massFraction(0,1)});
values *= volVars.density()*(this->spatialParams().velocity(scvf)*scvf.unitOuterNormal());
}
return values;
}
/*!
* \brief The initial values
*/
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
return PrimaryVariables({0.0, 0.0});
}
private:
static constexpr Scalar eps_ = 1e-7;
};
} // end namespace Dumux
#endif
// -*- 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 The test properties
*/
#ifndef DUMUX_TRACER_MULTIPHASE_TEST_PROPERTIES_HH
#define DUMUX_TRACER_MULTIPHASE_TEST_PROPERTIES_HH
#include <dune/grid/yaspgrid.hh>
#include <dumux/discretization/box.hh>
#include <dumux/discretization/cctpfa.hh>
#include <dumux/discretization/ccmpfa.hh>
#include <dumux/porousmediumflow/tracer/model.hh>
#include "fluidsystem.hh"
#include "problem.hh"
#include "spatialparams.hh"
namespace Dumux::Properties {
// Create new type tags
namespace TTag {
struct TracerTest { using InheritsFrom = std::tuple<Tracer>; };
struct TracerTestTpfa { using InheritsFrom = std::tuple<TracerTest, CCTpfaModel>; };
struct TracerTestMpfa { using InheritsFrom = std::tuple<TracerTest, CCMpfaModel>; };
struct TracerTestBox { using InheritsFrom = std::tuple<TracerTest, BoxModel>; };
} // end namespace TTag
// enable caching
template<class TypeTag>
struct EnableGridVolumeVariablesCache<TypeTag, TTag::TracerTest> { static constexpr bool value = true; };
template<class TypeTag>
struct EnableGridFluxVariablesCache<TypeTag, TTag::TracerTest> { static constexpr bool value = true; };
template<class TypeTag>
struct EnableGridGeometryCache<TypeTag, TTag::TracerTest> { static constexpr bool value = true; };
// Set the grid type
template<class TypeTag>
struct Grid<TypeTag, TTag::TracerTest> { using type = Dune::YaspGrid<2>; };
// Set the problem property
template<class TypeTag>
struct Problem<TypeTag, TTag::TracerTest> { using type = TracerTest<TypeTag>; };
// Set the spatial parameters
template<class TypeTag>
struct SpatialParams<TypeTag, TTag::TracerTest>
{
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using type = TracerTestSpatialParams<GridGeometry, Scalar>;
};
// Define whether mole(true) or mass (false) fractions are used
template<class TypeTag>
struct UseMoles<TypeTag, TTag::TracerTest> { static constexpr bool value = false; };
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::TracerTest> { using type = FluidSystems::TracerTest<GetPropType<TypeTag, Properties::Scalar>>; };
} // end namespace Dumux::Properties
#endif
// -*- 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 Definition of the spatial parameters for the tracer problem.
*/
#ifndef DUMUX_TRACER_TEST_SPATIAL_PARAMS_HH
#define DUMUX_TRACER_TEST_SPATIAL_PARAMS_HH
#include <dumux/common/parameters.hh>
#include <dumux/material/spatialparams/fv1p.hh>
namespace Dumux {
/*!
* \ingroup TracerTests
* \brief Definition of the spatial parameters for the tracer problem.
*/
template<class GridGeometry, class Scalar>
class TracerTestSpatialParams
: public FVSpatialParamsOneP<GridGeometry, Scalar,
TracerTestSpatialParams<GridGeometry, Scalar>>