Commit 122525f6 authored by Timo Koch's avatar Timo Koch
Browse files

Merge branch 'feature/exercise-grids' into 'master'

Feature/exercise grids

See merge request !17
parents c2298f7f 7aa9e1ec
......@@ -3,6 +3,7 @@ add_custom_target(test_exercises)
add_subdirectory(exercise-basic)
add_subdirectory(exercise-runtimeparams)
add_subdirectory(exercise-grids)
add_subdirectory(exercise-properties)
add_subdirectory(exercise-fluidsystem)
add_subdirectory(exercise-biomineralization)
......
......@@ -4,3 +4,5 @@ Click on the exercise to go the description
* [01 Exercise on the main file](./01-mainfile/README.md)
* [Exercise on runtime parameters](./exercise-runtimeparams/README.md)
* [Exercise on Grid Development](./exercise-grids/README.md)
# the grid exercise simulation program
dune_add_test(NAME exercise_grids
SOURCES exercise_grids.cc
CMD_ARGS exercise_grids.input)
# add tutorial to the common target
add_dependencies(test_exercises exercise_grids)
# add a symlink for the input file
dune_symlink_to_source_files(FILES "exercise_grids.input")
dune_symlink_to_source_files(FILES grids)
# 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-grids/`
We'll be editing the following files during this exercise:
* The __input file__: `exercise_grids.input` ,
* and the __problem file__: `injection2pproblem.hh`
Currently our grid is defined using the following input parameters in the input file:
```bash
[Grid]
UpperRight = 60 40
Cells = 24 16
```
The upper right point of our domain is defined using the first parameter, and all directions are split into a number of equal sized cells defined by the second parameter.
As you may have noticed, the resolution of our grid isn't quite optimal, and a finer grid could help to show a more realistic picture of the physics.
<br><br>
### Task 1: Applying a Global Refinement
<hr>
A simple strategy for addressing this resolution issue would be a global refinement.
In the Grid section of your input file, you can add a parameter "Refinement" which will apply a global refinement to the grid you provided. A refinement of value 2 will split each cell into two cells per dimension.
* > __Task 1__: Apply a global refinement of value 2. Make note of the increase in simulation run time.
The simulation should take much more time. Luckily there are other ways to resolve your grid that will not be as computationally expensive.
<br><br>
### Task 2: Changing the Grid Type
<hr>
Up until now, we have used a basic uniform structured grid. This grid is provided in the Dune environment and is called the YASP grid (Yet Another Structured Parallel Grid).
In the Problem file, the grid properties are defined using the following command.
```c++
// Set the grid type
SET_TYPE_PROP(Injection2pTypeTag, Grid, Dune::YaspGrid<2>);
```
This sets the Grid, which belongs to the Injection2pTypeTag type tag, and calls the manager with a basic YaspGrid in the second dimension.
If we look in the grid manager header file, `dumux/dumux/io/grid/gridmanager.hh`,
we will see that there are a few grid development options. Until now, we have used the most basic definition.
Let's change our grid manager to a slightly more complicated version. We will still call a YaspGrid, but can use the Dune coordinates system named TensorProductCoordinates rather than the basic version. The following class can be found in `dumux/dumux/io/grid/gridmanager.hh`.
```c++
GridManager<Dune::YaspGrid<dim, Dune::TensorProductCoordinates<ctype, dim> >>
```
A YaspGrid, of dimension `dim` is used, and the coordinate system `Dune::TensorProductCoordinates`, with arguments coordinate type `ctype` and dimension `dim`, is provided as a secondary argument. This allows us to switch to a grid format that allows for more cells in certain zones, and a cell grading.
To use this format, we need to set the grid property in our problem file accordingly. Our problem dimension is 2 (`dim`), and we will provide the coordinates as doubles (`ctype`).
To switch the property we should replace the `Dune::YaspGrid<2>` argument with `Dune::YaspGrid<dim, Dune::TensorProductCoordinates<ctype, dim> >`, and the `dim` and `ctype` with `2` and `double`
In our input file, we can now specify more specific Grid sizing directions: Positions, Cells, and Grading, each of which can be defined per direction.
In order to redevelop the exact same grid that was build beforehand, we can enter the following parameters to the Grid section of the input file.
```bash
[Grid]
Positions0 = 0 60
Positions1 = 0 40
Cells0 = 24
Cells1 = 16
Grading0 = 1.0
Grading1 = 1.0
```
* > __Task 2__: Following the format shown above, change the grid manager. Make the required changes to the input file and run the simulation.
<br><br>
### Task 3: Grid Zoning and Grading
<hr>
The Positions0 and Positions1 parameters will segment the grid into zones split at the values provided. In this case, the X direction
will be made of one zone between 0 and 60, and the Y direction will be split into one zone between 0 and 40.
```bash
[Grid]
Positions0 = 0 60
Positions1 = 0 40
```
The Cells0 and Cells1 parameters will define a number of cells to each of the zones. The first number provided will be applied to the first zone.
In the first zone in the X direction, the domain is split into 24 cells. In the first zone in the Y direction, the domain is split into 16 cells.
```bash
[Grid]
Cells0 = 24
Cells1 = 16
```
The Grading0 and Grading1 will apply a cell size grading throughout the zone in which it is applied.
The value 1.0 will split the domain evenly such that each cell has the same size in one direction.
Any other values will distribute the cell size unevenly throughout the domain.
These sizes are defined such that the next cell in the domain will be sized as the grading value multiplied by the previous cell size.
```bash
[Grid]
Grading0 = 1.0
Grading1 = 1.0
```
Positive values will make the cells closer to the start of the zone small, and those further from the start of the zone larger.
Negative values will make the cells closer to the start of the zone large, and those further from the start of the zone smaller.
Considering that the majority of the action in this simulation will occur along the interface between the aquifer and the aquitard, and along the injection axis, we should provide more resolution in these areas, and less in the other areas.
* > __Task 3__: Using the Zone defintions described above, edit your grid to focus your simulation on the areas in greatest concern.
One example of this would be the following definition.
```bash
[Grid]
Positions0 = 0 45 60
Positions1 = 0 25 30 35 40
Cells0 = 9 10
Cells1 = 5 5 5 2
Grading0 = 1.0 -1.3
Grading1 = 1.0 -1.3 1.3 1.0
```
<br><br>
### Task 4: Reading in a Structured Grid (.dgf and .msh)
<hr>
DuMuX can also read in grids from external files. There are two supported file types: `.dgf` grids, and `.msh` grids.
__`.dgf`__ grids (DUNE Grid Format grids) are developed within the dune environment and are often used in dumux tests.
(see: `https://www.dune-project.org/doxygen/2.4.1/group__DuneGridFormatParser.html`)
__`.msh`__ grids are made from the open-source Gmsh software, and is also supported by dumux.
(see: `http://gmsh.info//` and `http://gmsh.info//doc/texinfo/gmsh.html`)
We can read in simple .dgf files using basic YaspGrids, but if we want to use more complicated external grids, we would need to use other grid formats, like UGGrid or ALUGrid.
For this example, we can rewrite our Grid Properties to develop an ALUGrid. The grid manager class takes the following form:
```c++
GridManager<Dune::ALUGrid<dim, dimworld, elType, refinementType>>
```
Our `dim` and `dimworld` will both be 2, we can use `dune:cube` as our `elType` (element type), and for the `refinementType` we can use `Dune::nonconforming`.
There are two external structured grid files located in the `grid/` folder (`grid_structured.msh` and `grid_structured.dgf`). A path to one of these grids should be included in the input file.
* > __Task 4__: Edit the grid properties to set up an ALUGrid made of cubes. Call one of the strucutred meshes via the input file.
<br><br>
### Task 5: Reading in a Unstructured Grid (.msh)
<hr>
Unstructured grids are also supported for some Dumux models. An example of an Unstructured grid is located in the `grid/` folder (`grid_unstructured.msh`)
This grid is made up of triangles, or simplexes. This element type can be set in the grid properties section.
* > __Task 5__: Change the grid property element type to include `dune::simplex` instead of `dune::cube`. Read in the unstructured grid via the input file and run the simulation.
// -*- 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-Grids
*/
#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
# TODO: Task 1: Globally refine your grid
# TODO: Task 2: Develop even grid input parameters that can be split into zones
# TODO: Task 3: Using the zoning and grading parameters, redevelop your grid to optimize your refinement in the areas that matter.
# TODO: Task 4: Run your simulation the provided structured grid file ./grids/grid_structured.msh
# TODO: Task 5: Run your simulation the provided structured grid file ./grids/grid_unstructured.msh
[Problem]
Name = grid_exercise
OnlyPlotMaterialLaws = true
AquiferDepth = 2700.0 # m
InjectionDuration = 2.628e6 # in seconds, i.e. one month
[SpatialParams]
PermeabilityAquitard = 1e-15 # m^2
EntryPressureAquitard = 4.5e4 # Pa
PermeabilityAquifer = 1e-12 # m^2
EntryPressureAquifer = 1e4 # Pa
DGF
Interval
0 0 % first corner
60 40 % second corner
30 20 % cells in x and y direction
#
This diff is collapsed.
X = 60;
Y = 40;
Res = 1;
Point(1) = {0,0,0,Res};
Point(2) = {X,0,0,Res};
Point(3) = {X,Y,0,Res};
Point(4) = {0,Y,0,Res};
Line(5) = {1,2};
Line(6) = {2,3};
Line(7) = {3,4};
Line(8) = {4,1};
Line Loop(9) = {5,6,7,8};
Plane Surface(10) = 9;
This diff is collapsed.
// -*- 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-runtimeparams
*/
#ifndef DUMUX_EXGRIDS_INJECTION_PROBLEM_2P_HH
#define DUMUX_EXGRIDS_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>);
// TODO: Task 2: Replace the above Grid Property definition with a more flexible grid (Use Dune::TensorProductCoordinates)
// TODO: Task 4: Replace the above Grid Property definition to read in a external structured grid via a .msh file (Use Dune::ALUGrid and Dune:cube)
// TODO: Task 5: Replace the above Grid Property definition to read in a external unstructured grid via a .msh file (Use Dune::ALUGrid and Dune::simplex)
// 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_ = getParam<Scalar>("Problem.InjectionDuration");
}
/*!
* \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
*/
// \{