Skip to content
GitLab
Commit 221952c5 authored by Timo Koch's avatar Timo Koch
Browse files

Merge branch 'feature/outflow-bc-tracer' into 'master' 

[tracer][reference] implements an outflow boundary for tracer problem

Closes #784

See merge request !1806
parents 641ded71 692c5a0d
Hide whitespace changes
Inline Side-by-side
examples/1ptracer/README.md
View file @ 221952c5
......@@ -665,10 +665,21 @@ We use convenient declarations that we derive from the property system.
using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
using ElementVolumeVariables = typename GridVariables::GridVolumeVariables::LocalView;
using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
```
We create a bool saying whether mole or mass fractions are used
```cpp
static constexpr bool useMoles = getPropValue<TypeTag, Properties::UseMoles>();
```
We create additional int to make dimWorld and numComponents available in the problem
```cpp
static constexpr int dimWorld = GridView::dimensionworld;
static const int numComponents = FluidSystem::numComponents;
public:
```
......@@ -716,8 +727,37 @@ We assign different values, depending on wether mole concentrations or mass conc
}
return initialValues;
}
```
We implement an outflow boundary on the top of the domain and prescribe zero-flux Neumann boundary conditions on all other boundaries.
```cpp
NumEqVector neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFluxVariablesCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
NumEqVector values(0.0);
const auto& volVars = elemVolVars[scvf.insideScvIdx()];
const auto& globalPos = scvf.center();
```
This is the outflow boundary, where tracer is transported by advection
with the given flux field and diffusive flux is enforced to be zero
```cpp
if (globalPos[dimWorld-1] > this->gridGeometry().bBoxMax()[dimWorld-1] - eps_)
{
values = this->spatialParams().volumeFlux(element, fvGeometry, elemVolVars, scvf)
* volVars.massFraction(0, 0) * volVars.density(0)
/ scvf.area();
assert(values>=0.0 && "Volume flux at outflow boundary is expected to have a positive sign");
}
```
prescribe zero-flux Neumann boundary conditions elsewhere
```cpp
else
values = 0.0;
return values;
}
private:
```
......
examples/1ptracer/problem_tracer.hh
View file @ 221952c5
......@@ -158,8 +158,17 @@ class TracerTestProblem : public PorousMediumFlowProblem<TypeTag>
using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
using ElementVolumeVariables = typename GridVariables::GridVolumeVariables::LocalView;
using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
//We create a bool saying whether mole or mass fractions are used
static constexpr bool useMoles = getPropValue<TypeTag, Properties::UseMoles>();
//We create additional int to make dimWorld and numComponents available in the problem
static constexpr int dimWorld = GridView::dimensionworld;
static const int numComponents = FluidSystem::numComponents;
public:
// This is the constructor of our problem class:
......@@ -198,7 +207,32 @@ public:
return initialValues;
}
//We implement an outflow boundary on the top of the domain and prescribe zero-flux Neumann boundary conditions on all other boundaries.
NumEqVector neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFluxVariablesCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
NumEqVector values(0.0);
const auto& volVars = elemVolVars[scvf.insideScvIdx()];
const auto& globalPos = scvf.center();
// This is the outflow boundary, where tracer is transported by advection
// with the given flux field and diffusive flux is enforced to be zero
if (globalPos[dimWorld-1] > this->gridGeometry().bBoxMax()[dimWorld-1] - eps_)
{
values = this->spatialParams().volumeFlux(element, fvGeometry, elemVolVars, scvf)
* volVars.massFraction(0, 0) * volVars.density(0)
/ scvf.area();
assert(values>=0.0 && "Volume flux at outflow boundary is expected to have a positive sign");
}
//prescribe zero-flux Neumann boundary conditions elsewhere
else
values = 0.0;
return values;
}
private:
// We assign a private global variable for the epsilon:
......
test/porousmediumflow/tracer/1ptracer/CMakeLists.txt deleted 100644 → 0
View file @ 641ded71
dune_symlink_to_source_files(FILES "params.input")
dumux_add_test(NAME test_1ptracer
LABELS porousmediumflow tracer
SOURCES main.cc
CMAKE_GUARD HAVE_UMFPACK
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_1ptracer_transport-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_1ptracer-00010.vtu
${CMAKE_SOURCE_DIR}/test/references/test_1ptracer_pressure-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/1p.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_1ptracer params.input -Problem.Name test_1ptracer")
test/porousmediumflow/tracer/1ptracer/main.cc deleted 100644 → 0
View file @ 641ded71
// -*- 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.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/linear/pdesolver.hh>
#include <dumux/assembly/fvassembler.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager_yasp.hh>
#include "problem_1p.hh"
#include "problem_tracer.hh"
int main(int argc, char** argv) try
{
using namespace Dumux;
//! define the type tags for this problem
using OnePTypeTag = Properties::TTag::IncompressibleTest;
using TracerTypeTag = Properties::TTag::TracerTestCC;
//! 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
////////////////////////////////////////////////////////////
//! parse command line arguments
Parameters::init(argc, argv);
//////////////////////////////////////////////////////////////////////
// try to create a grid (from the given grid file or the input file)
/////////////////////////////////////////////////////////////////////
// only create the grid once using the 1p type tag
GridManager<GetPropType<OnePTypeTag, Properties::Grid>> gridManager;
gridManager.init();
//! we compute on the leaf grid view
const auto& leafGridView = gridManager.grid().leafGridView();
////////////////////////////////////////////////////////////
// setup & solve 1p problem on this grid
////////////////////////////////////////////////////////////
Dune::Timer timer;
//! create the finite volume grid geometry
using GridGeometry = GetPropType<OnePTypeTag, Properties::GridGeometry>;
auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
gridGeometry->update();
//! the problem (boundary conditions)
using OnePProblem = GetPropType<OnePTypeTag, Properties::Problem>;
auto problemOneP = std::make_shared<OnePProblem>(gridGeometry);
//! the solution vector
using SolutionVector = GetPropType<OnePTypeTag, Properties::SolutionVector>;
SolutionVector p(leafGridView.size(0));
//! the grid variables
using OnePGridVariables = GetPropType<OnePTypeTag, Properties::GridVariables>;
auto onePGridVariables = std::make_shared<OnePGridVariables>(problemOneP, gridGeometry);
onePGridVariables->init(p);
//! the assembler
using OnePAssembler = FVAssembler<OnePTypeTag, DiffMethod::analytic>;
auto assemblerOneP = std::make_shared<OnePAssembler>(problemOneP, gridGeometry, onePGridVariables);
//! the linear solver
using OnePLinearSolver = UMFPackBackend;
auto onePLinearSolver = std::make_shared<OnePLinearSolver>();
//! the pde system solver
LinearPDESolver<OnePAssembler, OnePLinearSolver> onePSolver(assemblerOneP, onePLinearSolver);
onePSolver.solve(p);
//! write output to vtk
using GridView = GetPropType<OnePTypeTag, Properties::GridView>;
Dune::VTKWriter<GridView> onepWriter(leafGridView);
onepWriter.addCellData(p, "p");
const auto& k = problemOneP->spatialParams().getKField();
onepWriter.addCellData(k, "permeability");
onepWriter.write("1p");
timer.stop();
const auto& comm = Dune::MPIHelper::getCollectiveCommunication();
std::cout << "Simulation took " << timer.elapsed() << " seconds on "
<< comm.size() << " processes.\n"
<< "The cumulative CPU time was " << timer.elapsed()*comm.size() << " seconds.\n";
////////////////////////////////////////////////////////////
// compute volume fluxes for the tracer model
////////////////////////////////////////////////////////////
using Scalar = GetPropType<OnePTypeTag, Properties::Scalar>;
std::vector<Scalar> volumeFlux(gridGeometry->numScvf(), 0.0);
using FluxVariables = GetPropType<OnePTypeTag, Properties::FluxVariables>;
auto upwindTerm = [](const auto& volVars) { return volVars.mobility(0); };
for (const auto& element : elements(leafGridView))
{
auto fvGeometry = localView(*gridGeometry);
fvGeometry.bind(element);
auto elemVolVars = localView(onePGridVariables->curGridVolVars());
elemVolVars.bind(element, fvGeometry, p);
auto elemFluxVars = localView(onePGridVariables->gridFluxVarsCache());
elemFluxVars.bind(element, fvGeometry, elemVolVars);
for (const auto& scvf : scvfs(fvGeometry))
{
const auto idx = scvf.index();
if (!scvf.boundary())
{
FluxVariables fluxVars;
fluxVars.init(*problemOneP, element, fvGeometry, elemVolVars, scvf, elemFluxVars);
volumeFlux[idx] = fluxVars.advectiveFlux(0, upwindTerm);
}
else
{
const auto bcTypes = problemOneP->boundaryTypes(element, scvf);
if (bcTypes.hasOnlyDirichlet())
{
FluxVariables fluxVars;
fluxVars.init(*problemOneP, element, fvGeometry, elemVolVars, scvf, elemFluxVars);
volumeFlux[idx] = fluxVars.advectiveFlux(0, upwindTerm);
}
}
}
}
////////////////////////////////////////////////////////////
// setup & solve tracer problem on the same grid
////////////////////////////////////////////////////////////
//! the problem (initial and boundary conditions)
using TracerProblem = GetPropType<TracerTypeTag, Properties::Problem>;
auto tracerProblem = std::make_shared<TracerProblem>(gridGeometry);
// set the flux from the 1p problem
tracerProblem->spatialParams().setVolumeFlux(volumeFlux);
//! the solution vector
SolutionVector x;
tracerProblem->applyInitialSolution(x);
auto xOld = x;
//! the grid variables
using GridVariables = GetPropType<TracerTypeTag, Properties::GridVariables>;
auto gridVariables = std::make_shared<GridVariables>(tracerProblem, gridGeometry);
gridVariables->init(x);
//! intialize the vtk output module
VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, tracerProblem->name());
using IOFields = GetPropType<TracerTypeTag, Properties::IOFields>;
IOFields::initOutputModule(vtkWriter); // Add model specific output fields
using VelocityOutput = GetPropType<TracerTypeTag, Properties::VelocityOutput>;
vtkWriter.addVelocityOutput(std::make_shared<VelocityOutput>(*gridVariables));
vtkWriter.write(0.0);
//! get some time loop parameters
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
auto dt = getParam<Scalar>("TimeLoop.DtInitial");
const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
//! instantiate time loop
auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0.0, dt, tEnd);
timeLoop->setMaxTimeStepSize(maxDt);
//! the assembler with time loop for instationary problem
using TracerAssembler = FVAssembler<TracerTypeTag, DiffMethod::analytic, /*implicit=*/false>;
auto assembler = std::make_shared<TracerAssembler>(tracerProblem, gridGeometry, gridVariables, timeLoop, xOld);
//! the linear solver
using TracerLinearSolver = ExplicitDiagonalSolver;
auto linearSolver = std::make_shared<TracerLinearSolver>();
//! the pde system solver
LinearPDESolver<TracerAssembler, TracerLinearSolver> 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(); do
{
// assemble, solve, update
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();
// write vtk output on check points
if (timeLoop->isCheckPoint())
vtkWriter.write(timeLoop->time());
// report statistics of this time step
timeLoop->reportTimeStep();
// set new dt
timeLoop->setTimeStepSize(dt);
} while (!timeLoop->finished());
timeLoop->finalize(leafGridView.comm());
////////////////////////////////////////////////////////////
// finalize, print dumux message to say goodbye
////////////////////////////////////////////////////////////
//! print dumux end message
if (mpiHelper.rank() == 0)
DumuxMessage::print(/*firstCall=*/false);
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;
}
test/porousmediumflow/tracer/1ptracer/params.input deleted 100644 → 0
View file @ 641ded71
[Grid]
UpperRight = 1 1
Cells = 50 50
[SpatialParams]
LensLowerLeft = 0.2 0.2
LensUpperRight = 0.8 0.8
Permeability = 1e-10 # [m^2]
PermeabilityLens = 1e-11 # [m^2]
[TimeLoop]
DtInitial = 10 # [s]
MaxTimeStepSize = 10
TEnd = 5000 # [s]
[Problem]
Name = tracer # name passed to the output routines
[Vtk]
AddVelocity = true
test/porousmediumflow/tracer/1ptracer/problem_1p.hh deleted 100644 → 0
View file @ 641ded71
// -*- 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 properties for the incompressible test
*/
#ifndef DUMUX_ONEP_TRACER_TEST_PROBLEM_HH
#define DUMUX_ONEP_TRACER_TEST_PROBLEM_HH
#include <dune/grid/yaspgrid.hh>
#include <dumux/discretization/cctpfa.hh>
#include <dumux/porousmediumflow/1p/model.hh>
#include <dumux/porousmediumflow/problem.hh>
#include <dumux/material/components/simpleh2o.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
#include <dumux/porousmediumflow/1p/incompressiblelocalresidual.hh>
#include "spatialparams_1p.hh"
namespace Dumux {
/*!
* \ingroup TracerTests
* \brief The properties for the incompressible test
*/
// forward declarations
template<class TypeTag>
class OnePTestProblem;
namespace Properties {
// Create new type tags
namespace TTag {
struct IncompressibleTest { using InheritsFrom = std::tuple<OneP, CCTpfaModel>; };
} // end namespace TTag
// Set the grid type
template<class TypeTag>
struct Grid<TypeTag, TTag::IncompressibleTest> { using type = Dune::YaspGrid<2>; };
// Set the problem type
template<class TypeTag>
struct Problem<TypeTag, TTag::IncompressibleTest> { using type = OnePTestProblem<TypeTag>; };
template<class TypeTag>
struct SpatialParams<TypeTag, TTag::IncompressibleTest>
{
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using type = OnePTestSpatialParams<GridGeometry, Scalar>;
};
template<class TypeTag>
struct LocalResidual<TypeTag, TTag::IncompressibleTest> { using type = OnePIncompressibleLocalResidual<TypeTag>; };
// the fluid system
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::IncompressibleTest>
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using type = FluidSystems::OnePLiquid<Scalar, Components::SimpleH2O<Scalar> >;
};
// Enable caching
template<class TypeTag>
struct EnableGridVolumeVariablesCache<TypeTag, TTag::IncompressibleTest> { static constexpr bool value = true; };
template<class TypeTag>
struct EnableGridFluxVariablesCache<TypeTag, TTag::IncompressibleTest> { static constexpr bool value = true; };
template<class TypeTag>
struct EnableGridGeometryCache<TypeTag, TTag::IncompressibleTest> { static constexpr bool value = true; };
} // end namespace Properties
template<class TypeTag>
class OnePTestProblem : public PorousMediumFlowProblem<TypeTag>
{
using ParentType = PorousMediumFlowProblem<TypeTag>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using Element = typename GridView::template Codim<0>::Entity;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
static constexpr int dimWorld = GridView::dimensionworld;
public:
OnePTestProblem(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry) {}
/*!
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary control volume.
*
* \param element The finite element
* \param scvf The sub-control volume face
*/
BoundaryTypes boundaryTypes(const Element &element,
const SubControlVolumeFace &scvf) const
{
BoundaryTypes values;
const auto globalPos = scvf.ipGlobal();
Scalar eps = 1.0e-6;
if (globalPos[dimWorld-1] < eps || globalPos[dimWorld-1] > this->gridGeometry().bBoxMax()[dimWorld-1] - eps)
values.setAllDirichlet();
else
values.setAllNeumann();
return values;
}
/*!
* \brief Evaluates the boundary conditions for a Dirichlet control volume.
*
* \param element The finite element
* \param scvf The sub-control volume face
*
* For this method, the \a values parameter stores primary variables.
*/
PrimaryVariables dirichlet(const Element &element,
const SubControlVolumeFace &scvf) const
{
const auto& pos = scvf.ipGlobal();
PrimaryVariables values(0);
values[0] = 1.0e+5*(1.1 - pos[dimWorld-1]*0.1);
return values;
}
/*!
* \brief Returns the temperature \f$\mathrm{[K]}\f$ for an isothermal problem.