Commit 070423ab authored by Katharina Heck's avatar Katharina Heck
Browse files

[feature] add exercise main file with readme and solution

parent 0beaa300
......@@ -6,5 +6,6 @@ add_subdirectory(exercise-runtimeparams)
add_subdirectory(exercise-properties)
add_subdirectory(exercise-fluidsystem)
add_subdirectory(exercise-biomineralization)
add_subdirectory(exercise-mainfile)
add_subdirectory(solution)
// -*- 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
*
* \brief The one-phase porousmediumflow problem for exercise 1
*/
#ifndef DUMUX_EX1_ONEP_TEST_PROBLEM_HH
#define DUMUX_EX1_ONEP_TEST_PROBLEM_HH
#include <dune/grid/yaspgrid.hh>
#include <dumux/material/components/simpleh2o.hh>
#include <dumux/material/components/h2o.hh>
#include <dumux/material/components/tabulatedcomponent.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
#include <dumux/discretization/cellcentered/tpfa/properties.hh>
#include <dumux/discretization/cellcentered/mpfa/properties.hh>
#include <dumux/discretization/box/properties.hh>
// TODO: dumux-course-task
// uncomment the incompressiblelocalresidual which is a specialization of the standard immisible localresidual for one phase incompressible cases and provides an analytic jacobian.
//#include <dumux/porousmediumflow/1p/incompressiblelocalresidual.hh>
#include <dumux/porousmediumflow/problem.hh>
#include <dumux/porousmediumflow/1p/model.hh>
#include "1pspatialparams.hh"
namespace Dumux
{
// forward declarations
template<class TypeTag> class OnePTestProblem;
namespace Properties
{
// create the type tag nodes. Here we define the incompressible type tag as well as the compressible type tag. The incompressible uses a different fluidsystem than the compressible
NEW_TYPE_TAG(OnePBase, INHERITS_FROM(OneP));
NEW_TYPE_TAG(OnePIncompressible, INHERITS_FROM(CCTpfaModel, OnePBase));
NEW_TYPE_TAG(OnePCompressible, INHERITS_FROM(CCTpfaModel, OnePBase));
// Set the grid type
SET_TYPE_PROP(OnePBase, Grid, Dune::YaspGrid<2>);
// Set the problem type
SET_TYPE_PROP(OnePBase, Problem, OnePTestProblem<TypeTag>);
// set the spatial params
SET_TYPE_PROP(OnePBase, SpatialParams, OnePTestSpatialParams<TypeTag>);
// the fluid system for incompressible tests
SET_PROP(OnePIncompressible, FluidSystem)
{
private:
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
public:
using type = FluidSystems::OnePLiquid<Scalar, Components::SimpleH2O<Scalar> >;
};
// TODO: dumux-course-task
// set the OneP Incompressible local residual for the OnePIncompressible type tag. This provides an analytic jacobian to be used for the analytic solution. Change that by setting:
//SET_TYPE_PROP(OnePIncompressible, LocalResidual, OnePIncompressibleLocalResidual<TypeTag>);
// the fluid system for compressible tests
SET_PROP(OnePCompressible, FluidSystem)
{
private:
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
public:
using type = FluidSystems::OnePLiquid<Scalar, Components::TabulatedComponent<Components::H2O<Scalar>>>;
};
// Disable caching (for testing purposes)
SET_BOOL_PROP(OnePBase, EnableGridVolumeVariablesCache, false);
SET_BOOL_PROP(OnePBase, EnableGridFluxVariablesCache, false);
SET_BOOL_PROP(OnePBase, EnableFVGridGeometryCache, false);
} // end namespace Properties
/*!
* \ingroup OnePTests
* \brief Test problem for the compressible one-phase model:
* \todo doc me!
* <tt>./test_box1pfv</tt> or
* <tt>./test_cc1pfv</tt>
*/
template<class TypeTag>
class OnePTestProblem : 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 PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
static constexpr int dimWorld = GridView::dimensionworld;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
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 values The boundary types for the conservation equations
* \param globalPos The position of the center of the finite volume
*/
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
BoundaryTypes values;
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 Evaluate the boundary conditions for a dirichlet
* control volume.
*
* \param values The dirichlet values for the primary variables
* \param globalPos The center of the finite volume which ought to be set.
*
* For this method, the \a values parameter stores primary variables.
*/
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables values(0);
values[0] = 1.0e5*(2.0 - globalPos[dimWorld-1]);
return values;
}
/*!
* \brief Evaluate the initial conditions
*
* \param globalPos The center of the finite volume which ought to be set.
*/
PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
{
return PrimaryVariables(1.0e5);
}
/*!
* \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 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 OnePTests
* \brief The spatial params the incompressible test
*/
#ifndef DUMUX_EX1_ONEP_TEST_SPATIAL_PARAMS_HH
#define DUMUX_EX1_ONEP_TEST_SPATIAL_PARAMS_HH
#include <dumux/porousmediumflow/properties.hh>
#include <dumux/material/spatialparams/fv1p.hh>
namespace Dumux {
/*!
* \ingroup OnePTests
* \brief The spatial parameters class for the test problem using the
* compressible 1p model
*/
template<class TypeTag>
class OnePTestSpatialParams
: public FVSpatialParamsOneP<typename GET_PROP_TYPE(TypeTag, FVGridGeometry),
typename GET_PROP_TYPE(TypeTag, Scalar),
OnePTestSpatialParams<TypeTag>>
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using GridView = typename FVGridGeometry::GridView;
using Element = typename GridView::template Codim<0>::Entity;
using FVElementGeometry = typename FVGridGeometry::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using ParentType = FVSpatialParamsOneP<FVGridGeometry, Scalar, OnePTestSpatialParams<TypeTag>>;
static constexpr int dimWorld = GridView::dimensionworld;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
public:
using PermeabilityType = Scalar;
OnePTestSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
{
permeability_ = getParam<Scalar>("SpatialParams.Permeability");
permeabilityLens_ = getParam<Scalar>("SpatialParams.PermeabilityLens");
lensLowerLeft_ = getParam<GlobalPosition>("SpatialParams.LensLowerLeft");
lensUpperRight_ = getParam<GlobalPosition>("SpatialParams.LensUpperRight");
}
/*!
* \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$.
*
* \param element The element
* \param scv The sub control volume
* \param elemSol The element solution vector
* \return the intrinsic permeability
*/
template<class ElementSolution>
PermeabilityType permeability(const Element& element,
const SubControlVolume& scv,
const ElementSolution& elemSol) const
{
if (isInLens_(scv.dofPosition()))
return permeabilityLens_;
else
return permeability_;
}
/*!
* \brief Define the porosity \f$\mathrm{[-]}\f$.
*
* \param globalPos The global position
*/
Scalar porosityAtPos(const GlobalPosition& globalPos) const
{ return 0.4; }
private:
bool isInLens_(const GlobalPosition &globalPos) const
{
for (int i = 0; i < dimWorld; ++i) {
if (globalPos[i] < lensLowerLeft_[i] + eps_ || globalPos[i] > lensUpperRight_[i] - eps_)
return false;
}
return true;
}
GlobalPosition lensLowerLeft_;
GlobalPosition lensUpperRight_;
Scalar permeability_, permeabilityLens_;
static constexpr Scalar eps_ = 1.5e-7;
};
} // end namespace Dumux
#endif
# the one-phase simulation program
dune_add_test(NAME exercise1_1p_a
SOURCES exercise1_1p_a.cc
COMPILE_DEFINITIONS TYPETAG=OnePIncompressible
CMD_ARGS exercise1_1p_a.input)
dune_add_test(NAME exercise1_1p_b
SOURCES exercise1_1p_b.cc
COMPILE_DEFINITIONS TYPETAG=OnePCompressible
CMD_ARGS exercise1_1p_b.input)
dune_add_test(NAME exercise1_1p_c
SOURCES exercise1_1p_c.cc
COMPILE_DEFINITIONS TYPETAG=OnePCompressible
CMD_ARGS exercise1_1p_c.input)
# here, add the two-phase non-isothermal simulation program
# add tutorial to the common target
add_dependencies(test_exercises exercise1_1p_a exercise1_1p_b exercise1_1p_c)
# add a symlink for the input file
dune_symlink_to_source_files(FILES "exercise1_1p_a.input" "exercise1_1p_b.input" "exercise1_1p_c.input")
# Exercise #1 (DuMuX course)
<br>
## Problem set-up
This exercise will make you familiar the program sequence in DuMux and how different levels of complexity can be realized in the main file according to the complexity of your physical problem.
In order to do so, there are three examples of one phase flow problems. Two examples (a and b) are stationary problems and the third example (c) is an instationary problem.
The stationary examples differ in the fluidssystems they are using which means they differ in the fluid properties (e.g. density, thermal conductivity etc). The first problem (a) uses an incompressible fluid which means that the density does not change with pressure changes. This makes is possible to solve the system linearly. The second problem uses a compressible fluid, that means the density is a function of pressure and we need to use a nonlinear solver.
To summarize the problems differ in:
* exercise1_1p_a: a one-phase incompressible, stationary problem
* exercise1_1p_b: a one-phase compressible, stationary problem
* exercise1_1p_c: a one-phase compressible, instationary problem
The problem set-up for all three examples is always the same: It is a two dimensional problem and the domain is $`1m x 1m`$. It is a heterogeneous set-up with a lens in the middle of the domain which has a lower permeability ($`1e-12 m^2`$ compared to $`1e-10 m^2`$ in the rest of the domain).
<img src="https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/raw/master/tutorial/extradoc/exercise1_1p_setup.png" width="1000">
In the beginning there is a uniform pressure of $`1e5 Pa`$. On the top and the bottom border dirichlet boundary conditions are set with a pressure of $`1e5 Pa`$ on top and $`2e5 Pa`$ on the bottom. At the sides there is no in- or outflow and there are no source terms.
## Preparing the exercise
* Navigate to the directory `dumux/tutorial/ex1`
<br><br>
### Task 1: Getting familiar with the code
<hr>
Locate all the files you will need for this exercise
* The __main file__ for the __1p incompressible, stationary__ problem : `exercise1_1p_a.cc`
* The __main file__ for the __1p compressible, stationary__ problem : `exercise1_1p_b.cc`
* The __main file__ for the __1p compressible, instationary__ problem : `exercise1_1p_c.cc`
* The shared __problem file__: `1pproblem.hh`
* The shared __spatial parameters file__: `1pspatialparams.hh`
* The __input file__ for the __1p incompressible, stationary__ problem: `exercise1_1p_a.input`
* The __input file__ for the __1p compressible, stationary__ problem: `exercise1_1p_b.input`
* The __input file__ for the __1p compressible, instationary__ problem: `exercise1_1p_c.input`
Please pay special attention to the similarities and differences in the three main files. The first main file is solved linearly and does not need a newton solver or any other nonlinear solver method. The second problem is a nonlinear problem and uses newton's method to solve the system. The third problem is nonlinear and additionally instationary. Therefore a time loop needs to be included in the main file.
The general structure of any main file in DuMux is:
* the specific problem TypeTag is defined for that problem. Here this is done via a compile definition in the CMakeLists.txt. You will have a closer look at the CMakeList and how to set up your own problem later in the exercise.
```c++
// define the type tag for this problem
using TypeTag = TTAG(TYPETAG);
```
* a gridmanager tries to create the grid either from a grid file or the input file
```c++
GridManager<typename GET_PROP_TYPE(TypeTag, Grid)> gridManager;
gridManager.init();
```
* we create the finite volume grid geometry, the problem, solutionvector and the gridvariables and initialize them. Additionally we initialize the vtkoutput. Each model has a predefined model specific output with relevant parameters for that model.
```c++
// create the finite volume grid geometry
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
fvGridGeometry->update();
// the problem (initial and boundary conditions)
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
auto problem = std::make_shared<Problem>(fvGridGeometry);
// the solution vector
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
SolutionVector x(fvGridGeometry->numDofs());
// the grid variables
using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
gridVariables->init(x);
// intialize the vtk output module
using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
VtkOutputFields::init(vtkWriter); //!< Add model specific output fields
vtkWriter.write(0.0);
```
* then we need to assemble and solve the system. Depending on the problem this can be done with a linear solver or with a nonlinear solver. If the problem is time dependent we additionally need a time loop. An example for that is given in exercise1_1p_c:
```c++
// get some time loop parameters
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
auto dt = getParam<Scalar>("TimeLoop.DtInitial");
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 Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
// the linear solver
using LinearSolver = ILU0BiCGSTABBackend;
auto linearSolver = std::make_shared<LinearSolver>();
// the non-linear solver
using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);
// set some check points for the time loop
timeLoop->setPeriodicCheckPoint(tEnd/10.0);
// time loop
timeLoop->start(); do
{
// set previous solution for storage evaluations
assembler->setPreviousSolution(xOld);
// linearize & solve
nonLinearSolver.solve(x, *timeLoop);
// 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
if (timeLoop->isCheckPoint())
vtkWriter.write(timeLoop->time());
// report statistics of this time step
timeLoop->reportTimeStep();
// set new dt as suggested by the newton solver
timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
} while (!timeLoop->finished());
timeLoop->finalize(leafGridView.comm());
```
<hr><br><br>
### Task 2: Compiling and running an executable
<hr>
* Change to the build-directory
```bash
cd ../../build-cmake/tutorial/ex1
```
* Compile all three executables `exercise1_1p_a` and `exercise1_1p_b and exercise1_1p_c`
```bash
make exercise1_1p_a exercise1_1p_b exercise1_1p_c
```
* Execute the three problems and inspect the result
```bash
./exercise1_1p_a
./exercise1_1p_b
./exercise1_1p_c
```
* you can look at the results with paraview
```bash
paraview injection-1p_a.pvd
```
<hr><br><br>
### Task 3: Analytical differentiation
<hr>
For the incompressible one phase problem it is possible to also have an analytic solution method. To implement that follow the directions in the `exercise1_1p_a.cc` and the `1pproblem.hh` marked by:
```c++
// TODO: dumux-course-task
```
// -*- 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
*
* \brief test for the one-phase CC model
*/
#include <config.h>
#include "1pproblem.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 <dune/istl/io.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/assembly/fvassembler.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 tag for this problem
using TypeTag = 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);
// initialize parameter tree
Parameters::init(argc, argv);
//////////////////////////////////////////////////////////////////////
// try to create a grid (from the given grid file or the input file)
/////////////////////////////////////////////////////////////////////
GridManager<typename GET_PROP_TYPE(TypeTag, Grid)> gridManager;
gridManager.init();
////////////////////////////////////////////////////////////
// run stationary 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 FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
fvGridGeometry->update();
// the problem (initial and boundary conditions)
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
auto problem = std::make_shared<Problem>(fvGridGeometry);
// the solution vector
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);