Commit 8065f08d authored by Timo Koch's avatar Timo Koch
Browse files

[exercises] Add exercises from last year's course

parent 40999f56
# add a target that builds all exercise solutions
add_custom_target(test_exercises)
add_subdirectory(exercise-basic)
add_subdirectory(exercise-properties)
add_subdirectory(exercise-fluidsystem)
add_subdirectory(solution)
# the immiscible two-phase simulation program
dune_add_test(NAME exercise1_2p
SOURCES exercise1_2p.cc
CMD_ARGS exercise1.input)
# the compositional two-phase two-component simulation program
dune_add_test(NAME exercise1_2p2c
SOURCES exercise1_2p2c.cc
CMD_ARGS exercise1.input)
# here, add the two-phase non-isothermal simulation program
# add tutorial to the common target
add_dependencies(test_exercises exercise1_2p exercise1_2p2c)
# add a symlink for the input file
dune_symlink_to_source_files(FILES "exercise1.input")
# Exercise #1 (DuMuX course)
<br>
## Problem set-up
N$`_2`$ is injected in an aquifer previously saturated with water with an injection rate of 0.0001 kg/(s*m$`^2`$).
The aquifer is situated 2700 m below sea level and the domain size is 60 m x 40 m. It consists of two layers, a moderately permeable one ($`\Omega_1`$) and a lower permeable one ($`\Omega_2`$).
<img src="https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/raw/master/tutorial/extradoc/exercise1_setup.png" width="1000">
## Preparing the exercise
* Navigate to the directory `dumux/tutorial/ex1`
_Exercise 1_ deals with two problems: a two-phase immiscible problem (__2p__) and a two-phase compositional problem (__2p2c__). They both set up the same scenario with the difference that the 2p2c assumes a miscible fluid state for the two fluids (water and gaseous N$`_2`$) and the 2p model assumes an immiscible fluid state.
<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 __2p__ problem : `exercise1_2p.cc`
* The __main file__ for the __2p2c__ problem : `exercise1_2p2c.cc`
* The __problem file__ for the __2p__ problem: `injection2pproblem.hh`
* The __problem file__ for the __2p2c__ problem: `injection2p2cproblem.hh`
* The shared __spatial parameters file__: `injection2pspatialparams.hh`
* The shared __input file__: `exercise1.input`
<hr><br><br>
### Task 2: Compiling and running an executable
<hr>
* Change to the build-directory
```bash
cd ../../build-cmake/tutorial/ex1
```
* Compile both executables `exercise1_2p` and `exercise1_2p2c`
```bash
make exercise1_2p exercise1_2p2c
```
* Execute the two problems and inspect the result
```bash
./exercise1_2p exercise1.input
./exercise1_2p2c exercise1.input
```
* you can look at the results with paraview
```bash
paraview injection-2p2c.pvd
```
<hr><br><br>
### Task 3: Changing input parameters
<hr>
In the input file `exercise1.input` you can find the following section
```ini
[SpatialParams]
EntryPressureAquitard = 4.5e4
EntryPressureAquifer = 1e4
```
* Change the values for the aquitard entry pressure in the input file to a lower value and compare the results with the previous solution. You don't need to recompile the executable.
<hr><br><br>
### Task 4: Runtime parameters
<hr>
The injection rate is currently hard-coded in `injection2p2cproblem.hh` to $`1e-4 kg/(s m^2)`$.
```c++
// set the Neumann values for the Nitrogen component balance
// convert from units kg/(s*m^2) to mole/(s*m^2)
values[Indices::contiNEqIdx] = -1e-4/FluidSystem::molarMass(FluidSystem::nCompIdx);
values[Indices::contiWEqIdx] = 0.0;
```
We want to be able to set it at runtime. To this end,
* use the following DuMuX macro to read a runtime parameter from the input file
```c++
// read the injection rate from the input file at run time
totalAreaSpecificInflow_ = getParam<TYPE>("GROUPNAME.PARAMNAME");
```
* Replace
`<TYPE>`,`<GROUPNAME>`,`<PARAMNAME>` by what is appropriate for your case:
* `<TYPE>` is the type of the parameter to read
* `<GROUPNAME>` is the group in the input file
* `<PARAMNAME>` is the name of the parameter in the input file
Note that due to the way the macro works, the names are specified as plain text within the "quotation marks".`<GROUPNAME>` and `<PARAMNAME>` need to be separated by a dot (.).
Follow the instructions given as a
```c++
// TODO: dumux-course-task
```
in the `injection2p2cproblem.hh` file and also remember to also set the parameter totalAreaSpecificInflow in the input file.
* Check the influence of that parameter on the simulation result by rerunning the simulation with different injection rates. Remember to also set the parameter totalAreaSpecificInflow in the input file.
Since you have changed your header file, you have to recompile the program.
<hr><br><br>
### 5. Setting up a new executable (for a non-isothermal simulation)
<hr>
* Copy the main file `exercise1_2p.cc` and rename it to `exercise1_2pni.cc`
* In `exercise1_2pni.cc`, include the header `injection2pniproblem.hh` instead of `injection2pproblem.hh`.
* In `exercise1_2pni.cc`, change `Injection2pCCTypeTag` to `Injection2pNICCTypeTag` in the line `using TypeTag = TTAG(Injection2pCCTypeTag);`
* Add a new executable in `CMakeLists.txt` by adding the lines
```cmake
# the two-phase non-isothermal simulation program
dune_add_test(NAME exercise1_2pni
SOURCES exercise1_2pni.cc
CMD_ARGS exercise1.input)
```
* Test that everything compiles without error
```bash
make # should rerun cmake
make exercise1_2pni # builds new executable
```
<hr><br><br>
### 6. Setting up a non-isothermal __2pni__ test problem
<hr>
* Open the file `injection2pniproblem.hh`. It is a copy of the `injection2pproblem.hh` with some useful comments on how to implement a non-isothermal model. Look for comments containing
```c++
// TODO: dumux-course-task
```
* The following set-up should be realized:
__Boundary conditions:__ Dirichlet conditions for the temperature with a temperature gradient of 0.03 K/m and a starting temperature of 283 K.
__Initial conditions:__ The same temperature gradient as in the boundary conditions with an additional lens (with position: 20 < x < 30, 5 < y < 35), which has an initial temperature of 380 K.
<img src="https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/raw/master/tutorial/extradoc/exercise1_nonisothermal.png" width="800">
The non-isothermal model requires additional parameters like the thermal conductivity of the solid component. They are already implemented and set in `exercise1.input`, you just need to _uncomment_ them.
// -*- 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 main file for the two-phase porousmediumflow problem of exercise 1
*/
#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/common/defaultusagemessage.hh>
#include <dumux/linear/amgbackend.hh>
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/assembly/fvassembler.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/discretization/methods.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager.hh>
#include "injection2pproblem.hh"
////////////////////////
// the main function
////////////////////////
int main(int argc, char** argv) try
{
using namespace Dumux;
// define the type tag for this problem
using TypeTag = TTAG(Injection2pCCTypeTag);
// 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 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<typename GET_PROP_TYPE(TypeTag, Grid)> gridManager;
gridManager.init();
////////////////////////////////////////////////////////////
// run 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 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());
problem->applyInitialSolution(x);
auto xOld = x;
// the grid variables
using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
gridVariables->init(x, xOld);
// get some time loop parameters
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
auto dt = getParam<Scalar>("TimeLoop.DtInitial");
// 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
// instantiate time loop
auto timeLoop = std::make_shared<TimeLoop<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 = AMGBackend<TypeTag>;
auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
// the non-linear solver
using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);
// time loop
timeLoop->start();
while (!timeLoop->finished())
{
// set previous solution for storage evaluations
assembler->setPreviousSolution(xOld);
//set time in problem (is used in time-dependent Neumann boundary condition)
problem->setTime(timeLoop->time()+timeLoop->timeStepSize());
// solve the non-linear system with time step control
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();
// report statistics of this time step
timeLoop->reportTimeStep();
// set new dt as suggested by the newton solver
timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
// output to vtk
vtkWriter.write(timeLoop->time());
}
timeLoop->finalize(leafGridView.comm());
////////////////////////////////////////////////////////////
// finalize, print dumux message to say goodbye
////////////////////////////////////////////////////////////
// print dumux end message
if (mpiHelper.rank() == 0)
{
Parameters::print();
DumuxMessage::print(/*firstCall=*/false);
}
return 0;
} // end main
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;
}
[TimeLoop]
DtInitial = 3600 # in seconds
TEnd = 3.154e9 # in seconds, i.e ten years
[Grid]
LowerLeft = 0 0
UpperRight = 60 40
Cells = 24 16
[Problem]
Name = injection
OnlyPlotMaterialLaws = true
AquiferDepth = 2700.0 # m
InjectionDuration = 2.628e6 # in seconds, i.e. one month
#TODO: dumux-course-task:
#set totalAreaSpecificInflow
[SpatialParams]
PermeabilityAquitard = 1e-15 # m^2
EntryPressureAquitard = 4.5e4 # Pa
PermeabilityAquifer = 1e-12 # m^2
EntryPressureAquifer = 1e4 # Pa
# these parameters are only used in the nonisothermal model. Uncomment them for that
#[Component]
#SolidDensity = 2700 # solid density of granite
#SolidThermalConductivity = 2.8 # solid thermal conducitivity of granite
#SolidHeatCapacity = 790 # solid heat capacity of granite
// -*- 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 main file for the two-phase porousmediumflow problem of exercise 1
*/
#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/common/defaultusagemessage.hh>
#include <dumux/linear/amgbackend.hh>
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/assembly/fvassembler.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/discretization/methods.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager.hh>
// The problem file, where setup-specific boundary and initial conditions are defined.
#include "injection2pproblem.hh"
////////////////////////
// the main function
////////////////////////
int main(int argc, char** argv) try
{
using namespace Dumux;
// define the type tag for this problem
using TypeTag = TTAG(Injection2pCCTypeTag);
// 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 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<typename GET_PROP_TYPE(TypeTag, Grid)> gridManager;
gridManager.init();
////////////////////////////////////////////////////////////
// run 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 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());
problem->applyInitialSolution(x);
auto xOld = x;
// the grid variables
using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
gridVariables->init(x, xOld);
// get some time loop parameters
// getParam<TYPE>("GROUPNAME.PARAMNAME") reads and sets parameter PARAMNAME
// of type TYPE given in the group GROUPNAME from the input file
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
auto dt = getParam<Scalar>("TimeLoop.DtInitial");
// 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
// instantiate time loop
auto timeLoop = std::make_shared<TimeLoop<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 = AMGBackend<TypeTag>;
auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
// the non-linear solver
using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);
// time loop
timeLoop->start();
while (!timeLoop->finished())
{
// set previous solution for storage evaluations
assembler->setPreviousSolution(xOld);
//set time in problem (is used in time-dependent Neumann boundary condition)
problem->setTime(timeLoop->time()+timeLoop->timeStepSize());
// solve the non-linear system with time step control
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();
// report statistics of this time step
timeLoop->reportTimeStep();
// set new dt as suggested by the newton solver
timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
// output to vtk
vtkWriter.write(timeLoop->time());
}
timeLoop->finalize(leafGridView.comm());
////////////////////////////////////////////////////////////
// finalize, print dumux message to say goodbye
////////////////////////////////////////////////////////////
// print dumux end message
if (mpiHelper.rank() == 0)
{
Parameters::print();
DumuxMessage::print(/*firstCall=*/false);
}
return 0;
} // end main
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;
}
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4: