Commit 1ae9db04 authored by Felix Weinhardt's avatar Felix Weinhardt Committed by Katharina Heck
Browse files

added 1ptracer to examples

parent 81110be6
dune_symlink_to_source_files(FILES "params.input")
dune_add_test(NAME example_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}/example_1ptracer params.input -Problem.Name example_1ptracer")
This tutorial was copied from /home/felix/Dokumente/Aktuelles_U/19_03_27_DumuxTag/dumux/dumux/test/porousmediumflow/tracer/1ptracer.
# One-phase flow with with random permeability distribution (and a tracer model)
## Problem set-up
This example combines a stationary One-phase flow in a porous medium with randomly generated permeabilities with a tracer model. In the first step, the groundwater-velocity is evaluated under stationary conditions. Based on the velocity field the tracer model is running instationary.
## Random permeability generation
### spatialparams_1p.hh
The follwing code can be found in lines 64-72.
```C++
// generate random fields
std::mt19937 rand(0);
std::lognormal_distribution<Scalar> K(std::log(permeability_), std::log(permeability_)*0.1);
std::lognormal_distribution<Scalar> KLens(std::log(permeabilityLens_), std::log(permeabilityLens_)*0.1);
for (const auto& element : elements(fvGridGeometry->gridView()))
{
const auto eIdx = fvGridGeometry->elementMapper().index(element);
const auto globalPos = element.geometry().center();
K_[eIdx] = isInLens_(globalPos) ? KLens(rand) : K(rand);
}
```
Two lognormal distributions are defined: K and KLens with different means and standard deviations. In this case the mean is defined as the log of the permeability given from the input file and the standard deviation is 10% of the mean.
Within a loop through all elements the randomly generated permeabilities are assigned according to their position in the domain (inside or outside the lense).
## Tracer model
...
\ No newline at end of file
// -*- 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 "problem_1p.hh"
#include "problem_tracer.hh"
#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/assembly/fvassembler.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager.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
////////////////////////////////////////////////////////////
//! create the finite volume grid geometry
using FVGridGeometry = GetPropType<OnePTypeTag, Properties::FVGridGeometry>;
auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
fvGridGeometry->update();
//! the problem (boundary conditions)
using OnePProblem = GetPropType<OnePTypeTag, Properties::Problem>;
auto problemOneP = std::make_shared<OnePProblem>(fvGridGeometry);
//! the solution vector
using JacobianMatrix = GetPropType<OnePTypeTag, Properties::JacobianMatrix>;
using SolutionVector = GetPropType<OnePTypeTag, Properties::SolutionVector>;
SolutionVector p(leafGridView.size(0));
//! the linear system
auto A = std::make_shared<JacobianMatrix>();
auto r = std::make_shared<SolutionVector>();
//! the grid variables
using OnePGridVariables = GetPropType<OnePTypeTag, Properties::GridVariables>;
auto onePGridVariables = std::make_shared<OnePGridVariables>(problemOneP, fvGridGeometry);
onePGridVariables->init(p);
//! the assembler
using OnePAssembler = FVAssembler<OnePTypeTag, DiffMethod::analytic>;
auto assemblerOneP = std::make_shared<OnePAssembler>(problemOneP, fvGridGeometry, onePGridVariables);
assemblerOneP->setLinearSystem(A, r);
Dune::Timer timer;
// assemble the local jacobian and the residual
Dune::Timer assemblyTimer; std::cout << "Assembling linear system ..." << std::flush;
assemblerOneP->assembleJacobianAndResidual(p);
assemblyTimer.stop(); std::cout << " took " << assemblyTimer.elapsed() << " seconds." << std::endl;
// we solve Ax = -r
(*r) *= -1.0;
//! solve the 1p problem
using LinearSolver = UMFPackBackend;
Dune::Timer solverTimer; std::cout << "Solving linear system ..." << std::flush;
auto linearSolver = std::make_shared<LinearSolver>();
linearSolver->solve(*A, p, *r);
solverTimer.stop(); std::cout << " took " << solverTimer.elapsed() << " seconds." << std::endl;
//! update the grid variables
Dune::Timer updateTimer; std::cout << "Updating variables ..." << std::flush;
onePGridVariables->update(p);
updateTimer.elapsed(); std::cout << " took " << updateTimer.elapsed() << std::endl;
//! 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(fvGridGeometry->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(*fvGridGeometry);
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>(fvGridGeometry);
// set the flux from the 1p problem
tracerProblem->spatialParams().setVolumeFlux(volumeFlux);
//! the solution vector
SolutionVector x(leafGridView.size(0));
tracerProblem->applyInitialSolution(x);
auto xOld = x;
//! the grid variables
using GridVariables = GetPropType<TracerTypeTag, Properties::GridVariables>;
auto gridVariables = std::make_shared<GridVariables>(tracerProblem, fvGridGeometry);
gridVariables->init(x);
//! 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, fvGridGeometry, gridVariables, timeLoop);
assembler->setLinearSystem(A, r);
//! 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);
/////////////////////////////////////////////////////////////////////////////////////////////////
// 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
{
// set previous solution for storage evaluations
assembler->setPreviousSolution(xOld);
Dune::Timer assembleTimer;
assembler->assembleJacobianAndResidual(x);
assembleTimer.stop();
// solve the linear system A(xOld-xNew) = r
Dune::Timer solveTimer;
SolutionVector xDelta(x);
linearSolver->solve(*A, xDelta, *r);
solveTimer.stop();
// update solution and grid variables
updateTimer.reset();
x -= xDelta;
gridVariables->update(x);
updateTimer.stop();
// statistics
const auto elapsedTot = assembleTimer.elapsed() + solveTimer.elapsed() + updateTimer.elapsed();
std::cout << "Assemble/solve/update time: "
<< assembleTimer.elapsed() << "(" << 100*assembleTimer.elapsed()/elapsedTot << "%)/"
<< solveTimer.elapsed() << "(" << 100*solveTimer.elapsed()/elapsedTot << "%)/"
<< updateTimer.elapsed() << "(" << 100*updateTimer.elapsed()/elapsedTot << "%)"
<< std::endl;
// 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 (Dumux::ParameterException &e)
{
std::cerr << std::endl << e << " ---> Abort!" << std::endl;
return 1;
}
catch (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 (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;
}
[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
// -*- 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 FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using type = OnePTestSpatialParams<FVGridGeometry, 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 EnableFVGridGeometryCache<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::FVGridGeometry>::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
static constexpr int dimWorld = GridView::dimensionworld;
public:
OnePTestProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry) {}
/*!
* \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->fvGridGeometry().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.
*
* This is not specific to the discretization. By default it just
* throws an exception so it must be overloaded by the problem if
* no energy equation is used.
*/
Scalar temperature() const
{
return 283.15; // 10°C
}
};
} // 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 Definition of a problem, for the tracer problem:
* A rotating velocity field mixes a tracer band in a porous groundwater reservoir.
*/
#ifndef DUMUX_TRACER_TEST_PROBLEM_HH
#define DUMUX_TRACER_TEST_PROBLEM_HH
#include <dune/grid/yaspgrid.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 Definition of a problem, for the tracer problem:
* A rotating velocity field mixes a tracer band in a porous groundwater reservoir.
*/
template <class TypeTag>
class TracerTestProblem;
namespace Properties {
// Create new type tags
namespace TTag {
struct TracerTest { using InheritsFrom = std::tuple<Tracer>; };