Commit d885ac02 authored by Edward 'Ned' Coltman's avatar Edward 'Ned' Coltman
Browse files

Exercise-runtimeparams exercise files

parent 6d1f17c2
# Exercise Runtime Parameters (DuMuX course)
<br>
## Problem set-up
Here we will expand on what we've covered in the mainfile exercise, and the problem set up will remain the same.
## Preparing the exercise
* Navigate to the directory `dumux-course/exercises/exercise-runtimeparams/`
<br><br>
### Task 1: Understanding Input Parameters
<hr>
For this task we will edit the following files:
* The shared __problem file__: `injection2pproblem.hh`
* And the shared __input file__: `exercise_runtimeparams.input`
Parameters can either be directly defined within your program, or specified
via the input file. Within every main file, (*.cc), the following function
is called, which will read in the imput file parameters
```c++
// parse command line arguments and input file
Parameters::init(argc, argv);
```
This input file should either be named the same as the executable file, with a
trailing ".input", as is standard in our CMake system, or be explicitly written
as the first shell argument after the executable file is called.
```bash
./exercise_runtimeparams
(Calls a file (exercise_runtimeparams.input) as the default input file.)
```
```bash
./exercise_runtimeparams exercise1.input
(Calls the input file provided (exercise1.input) as the input file.)
```
In the input file `exercise_runtimeparams.input` you can find the following section
```ini
[SpatialParams]
PermeabilityAquitard = 1e-15 # m^2
EntryPressureAquitard = 4.5e4 # Pa
PermeabilityAquifer = 1e-12 # m^2
EntryPressureAquifer = 1e4 # Pa
```
When a parameter is defined directly within your program, you'll need to recompile your program everytime you change the value. When a parameter is passed via the imput file, this is not the case. If we decided to vary the entry pressure in our geologic units a few times via the parameters listed above, there would be no need to recompile between simulation runs.
* > __Task 1__: Change the aquitard's 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.
<br><br>
### Task 2: Setting a variable to collect a runtime parameter
<hr>
Let's free one of the variables in this exercise. Within the class `injection2p2cproblem.hh`, in the first line of the neumann boundary condition definition function neumannAtPos(-), the injection rate is defined as $`1e-4 kg/(s m^2)`$.
```c++
// inject nitrogen. negative values mean injection
// units kg/(s*m^2)
values[Indices::conti0EqIdx + FluidSystem::N2Idx]= -1e-4/FluidSystem::molarMass(FluidSystem::N2Idx);
values[Indices::conti0EqIdx + FluidSystem::H2OIdx] = 0.0;
```
This parameter may need to change, and if we choose to always change this within the class,
we will need to recompile every time.
Instead of hard defining this parameter within the funciton, we can set a variable to read into our parameter tree via the input file and use this in our function instead.
To do this, there are two functions defined in `dumux/dumux/common/parameters.hh`, `getParam()` and `getParamFromGroup()`. They use the following format:
```c++
variable_ = getParam<TYPE>("GROUPNAME.PARAMNAME");
```
or
```c++
variable_ = getParamFromGroup<TYPE>("GROUPNAME", "PARAMNAME");
```
`<TYPE>`,`<GROUPNAME>`,`<PARAMNAME>` should be appropriatly defined for your variable:
* `<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
An example of this is already performed in the problem constructor. The Injection Duration (injectionDuration_) is defined via the imput file, and can then be used later in the problem.
```c++
// depth of the aquifer, units: m
aquiferDepth_ = getParam<Scalar>("Problem.AquiferDepth");
// the duration of the injection, units: second
injectionDuration_ = getParamFromGroup<Scalar>("Problem","InjectionDuration");
```
The injection duration parameter is located in the "Problem" group, is named "InjectionDuration", and has the type "Scalar".
This variable should then be defined as Scalar at the bottom of this problem class in the private section.
```c++
Scalar injectionDuration_;
```
In the input file, within the group [Problem], a value is set to the Parameter name called in the class.
```ini
[Problem]
InjectionDuration = 2.628e6 # in seconds, i.e. one month
```
* > __Task 2__: The goal is to replace the value '-1e-4' in
```c++
values[Indices::conti0EqIdx + FluidSystem::N2Idx]= -1e-4/FluidSystem::molarMass(FluidSystem::N2Idx);
```
with a runtime variable. (2a) Develop a new variable called totalAreaSpecificInflow_, (2b) record a value to this variable from a path in the input file, and (2c) incorporate this variable into the injection boundary condition.
When your problem file and the input file are edited, compile and run the exercise.
```bash
make exercise_runtimeparams && ./exercise_runtimeparams exercise_runtimeparams.input
```
<br><br>
### Task 3: Default Values for Runtime Parameters
<hr>
In the case that no path to a required runtime parameter is provided in the input file, the program will no longer run. You can try this if you like by removing the TotalAreaSpecificInflow from the input file. The resulting error message should look like this:
```bash
loggingparametertree.hh:316:
Key Problem.TotalAreaSpecificInflow not found in the parameter tree ---> Abort!
```
To avoid this, we can place a default value in the variable defintion. Whenever the parameter is specified in the input file, this default value in your class will be overwritten, but in the case that no value is provided in the input file, this default value will be used. In order to do this, follow the following template.
```c++
variable_ = getParam<TYPE>("GROUPNAME.PARAMNAME", DEFAULTVALUE);
```
```c++
variable_ = getParamFromGroup<TYPE>("GROUPNAME","PARAMNAME", DEFAULTVALUE);
```
* > __Task 3__: Set up the totalAreaSpecificInflow_ variable to record a default value of -1e-4, and run this with and without a provided value of -1e-3 in the input file.
<br><br>
### Task 4: Other Runtime Parameter Functions
<hr>
Setting default values for variables defined with runtime parameters can also lead to problems. If one runtime parameter from the input file is set to multiple different variables, each with a different Default value, changing the variable in one location can lead to unexpected changes elsewhere. On top of this, in DUMUX, there are a few base variables that are set with default values for all dumux simulations. These can be found in the header file `dumux/common/parameters.hh` in the function `globalDefaultParameters`.
One way to check this is to use either the `hasParam()` or the `hasParamInGroup()` function. These bool functions will check to see if a parameter is read in via the input file. These functions are also both defined in the `dumux/dumux/common/parameters.hh` class, and follow a similar format to that of `getParam` and `getParamFromGroup`
An example of this would look like this:
```c++
if (hasParam("GROUPNAME.PARAMNAME"))
std::cout << "Parameter value is read from the input file." << std::endl;
else
std::cout << "Using the default parameter value." << std::endl;
```
or,
```c++
if (hasParamInGroup("GROUPNAME","PARAMNAME"))
std::cout << "Parameter value is read from the input file." << std::endl;
else
std::cout << "Using the default parameter value." << std::endl;
```
Using these functions we can better check which parameter values are being included in our program.
* > __Task 4__: Using one of the bool hasParam functions, place output in the problem file to alert the user to where the parameter value comes from.
// -*- 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 runtime parameters
*/
#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;
}
[TimeLoop]
DtInitial = 3600 # in seconds
TEnd = 3.154e9 # in seconds, i.e ten years
[Grid]
UpperRight = 60 40
Cells = 24 16
[Problem]
Name = runtimeparams_exercise
OnlyPlotMaterialLaws = true
AquiferDepth = 2700.0 # m
InjectionDuration = 2.628e6 # in seconds, i.e. one month
# TODO: Task 2: Create a parameter called "TotalAreaSpecificInflow"
[SpatialParams]
PermeabilityAquitard = 1e-15 # m^2
# TODO: Task 1: Change the Aquitard's Entry Pressure
EntryPressureAquitard = 4.5e4 # Pa
PermeabilityAquifer = 1e-12 # m^2
EntryPressureAquifer = 1e4 # Pa
// -*- 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 two-phase porousmediumflow problem for exercise runtime parameters
*/
#ifndef DUMUX_EXRUNTIMEPARAMS_INJECTION_PROBLEM_2P_HH
#define DUMUX_EXRUNTIMEPARAMS_INJECTION_PROBLEM_2P_HH
#include <dumux/discretization/cellcentered/tpfa/properties.hh>
#include <dumux/porousmediumflow/2p/model.hh>
#include <dumux/porousmediumflow/problem.hh>
#include <dumux/material/fluidsystems/h2on2.hh>
#include "injection2pspatialparams.hh"
namespace Dumux {
// forward declare problem
template <class TypeTag>
class InjectionProblem2P;
namespace Properties {
// define the TypeTag for this problem with a cell-centered two-point flux approximation spatial discretization.
NEW_TYPE_TAG(Injection2pTypeTag, INHERITS_FROM(TwoP));
NEW_TYPE_TAG(Injection2pCCTypeTag, INHERITS_FROM(CCTpfaModel, Injection2pTypeTag));
// Set the grid type
SET_TYPE_PROP(Injection2pTypeTag, Grid, Dune::YaspGrid<2>);
// Set the problem property
SET_TYPE_PROP(Injection2pTypeTag, Problem, InjectionProblem2P<TypeTag>);
// Set the spatial parameters
SET_TYPE_PROP(Injection2pTypeTag, SpatialParams,
InjectionSpatialParams<typename GET_PROP_TYPE(TypeTag, FVGridGeometry),
typename GET_PROP_TYPE(TypeTag, Scalar)>);
// Set fluid configuration
SET_TYPE_PROP(Injection2pTypeTag, FluidSystem, FluidSystems::H2ON2<typename GET_PROP_TYPE(TypeTag, Scalar), /*useComplexRelations=*/ false>);
} // end namespace Properties
/*!
* \ingroup TwoPModel
* \ingroup ImplicitTestProblems
* \brief Gas injection problem where a gas (here nitrogen) is injected into a fully
* water saturated medium. During buoyancy driven upward migration the gas
* passes a high temperature area.
*
* The domain is sized 60 m times 40 m.
*
* For the mass conservation equation neumann boundary conditions are used on
* the top, on the bottom and on the right of the domain, while dirichlet conditions
* apply on the left boundary.
*
* Gas is injected at the right boundary from 7 m to 15 m at a rate of
* 0.001 kg/(s m), the remaining neumann boundaries are no-flow
* boundaries.
*
* At the dirichlet boundaries a hydrostatic pressure and a gas saturation of zero a
*
* This problem uses the \ref TwoPModel model.
*/
template<class TypeTag>
class InjectionProblem2P : public PorousMediumFlowProblem<TypeTag>
{
using ParentType = PorousMediumFlowProblem<TypeTag>;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
enum { dimWorld = GridView::dimensionworld };
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
public:
InjectionProblem2P(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
{
// initialize the tables of the fluid system
FluidSystem::init(/*tempMin=*/273.15,
/*tempMax=*/423.15,
/*numTemp=*/50,
/*pMin=*/0.0,
/*pMax=*/30e6,
/*numP=*/300);
// name of the problem and output file
// getParam<TYPE>("GROUPNAME.PARAMNAME") reads and sets parameter PARAMNAME
// of type TYPE given in the group GROUPNAME from the input file
name_ = getParam<std::string>("Problem.Name");
// depth of the aquifer, units: m
aquiferDepth_ = getParam<Scalar>("Problem.AquiferDepth");
// the duration of the injection, units: second
injectionDuration_ = getParamFromGroup<Scalar>("Problem","InjectionDuration");
//TODO: Task 2: Set a variable "TotalAreaSpecificInflow" to read in a value from the parameter tree via the input file
//TODO: Task 3: Set a default value for the above parameter.
//TODO: Task 4: Provide output describing where the parameter value comes from using parameter bool functions.
}
/*!
* \name Problem parameters
*/
// \{
/*!
* \brief Returns the problem name
*
* This is used as a prefix for files generated by the simulation.
*/
std::string name() const
{ return name_+"-2p"; }
/*!
* \brief Returns the temperature \f$ K \f$
*/
Scalar temperature() const
{
return 273.15 + 30; // [K]
}
// \}
/*!
* \name Boundary conditions
*/
// \{
/*!
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary segment.
*
* \param bcTypes The boundary types for the conservation equations
* \param globalPos The position for which the bc type should be evaluated
*/
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
BoundaryTypes bcTypes;
// set the left of the domain (with the global position in "0 = x" direction as a Dirichlet boundary
if (globalPos[0] < eps_)
bcTypes.setAllDirichlet();
// set all other as Neumann boundaries
else
bcTypes.setAllNeumann();
return bcTypes;
}
/*!
* \brief Evaluates the boundary conditions for a Dirichlet
* boundary segment
*
* \param globalPos The global position
*/
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
return initialAtPos(globalPos);
}
/*!
* \brief Evaluate the boundary conditions for a neumann
* boundary segment.
*
* \param values Stores the Neumann values for the conservation equations in
* \f$ [ \textnormal{unit of conserved quantity} / (m^(dim-1) \cdot s )] \f$
* \param globalPos The position of the integration point of the boundary segment.
*
* For this method, the \a values parameter stores the mass flux
* in normal direction of each phase. Negative values mean influx.
*/
PrimaryVariables neumannAtPos(const GlobalPosition &globalPos) const
{
// initialize values to zero, i.e. no-flow Neumann boundary conditions
PrimaryVariables values(0.0);
// if we are inside the injection zone set inflow Neumann boundary conditions
// using < boundary + eps_ or > boundary - eps_ is safer for floating point comparisons
// than using <= or >= as it is robust with regard to imprecision introduced by rounding errors.
if (time_ < injectionDuration_
&& globalPos[1] < 15 + eps_ && globalPos[1] > 7 - eps_ && globalPos[0] > 0.9*this->fvGridGeometry().bBoxMax()[0])
{
// inject nitrogen. negative values mean injection
// units kg/(s*m^2)
//TODO: Task 2: incorporate "totalAreaSpecificInflow_" into the injection boundary condition
values[Indices::conti0EqIdx + FluidSystem::N2Idx]= -1e-4/FluidSystem::molarMass(FluidSystem::N2Idx);
values[Indices::conti0EqIdx + FluidSystem::H2OIdx] = 0.0;
}