Commit c78aa6ef authored by Timo Koch's avatar Timo Koch
Browse files

[test][co2] Port co2 test

parent a8df4839
add_input_file_links()
dune_symlink_to_source_files(FILES grids)
# build target for the CO2 test problem
# Ignore the porosity for all box models since it is defined element-wise in these test
# but the default 2p2c implementation outputs porosity per vertex.
# Depending on the order of the elements, the porosity would differ in these cases.
dune_add_test(NAME test_co2_box
SOURCES test_co2_fv.cc
COMPILE_DEFINITIONS TYPETAG=HeterogeneousBoxTypeTag
CMAKE_GUARD "( dune-alugrid_FOUND AND DUNE_GRID_EXPERIMENTAL_GRID_EXTENSIONS )"
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/co2box-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/co2_box-00019.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_co2_box test_co2_fv.input -Problem.Name co2_box"
--zeroThreshold {"porosity":1})
add_dumux_test(test_boxco2 test_boxco2 test_boxco2.cc
python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
--script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/co2box-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/heterogeneousbox-00019.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_boxco2"
--zeroThreshold {"porosity":1})
dune_add_test(NAME test_co2_tpfa
SOURCES test_co2_fv.cc
COMPILE_DEFINITIONS TYPETAG=HeterogeneousCCTpfaTypeTag
CMAKE_GUARD "( dune-alugrid_FOUND AND DUNE_GRID_EXPERIMENTAL_GRID_EXTENSIONS )"
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/co2cc-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/co2_tpfa-00018.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_co2_tpfa test_co2_fv.input -Problem.Name co2_tpfa")
add_dumux_test(test_ccco2 test_ccco2 test_ccco2.cc
python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
--script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/co2cc-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/heterogeneouscc-00018.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_ccco2")
add_dumux_test(test_restartco2 test_boxco2 test_boxco2.cc
python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
--script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/co2box-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/heterogeneousbox-00019.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_boxco2 -ParameterFile ${CMAKE_CURRENT_SOURCE_DIR}/test_restartco2.input"
--zeroThreshold {"porosity":1})
# # restart test (currently disbled because restart is not yet implemented)
# dune_add_test(NAME test_co2_box_restart
# TARGET test_co2_box
# COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
# CMD_ARGS --script fuzzy
# --files ${CMAKE_SOURCE_DIR}/test/references/co2box-reference.vtu
# ${CMAKE_CURRENT_BINARY_DIR}/heterogeneousbox-00019.vtu
# --command "${CMAKE_CURRENT_BINARY_DIR}/test_co2_box -ParameterFile ${CMAKE_CURRENT_SOURCE_DIR}/test_restartco2.input"
# --zeroThreshold {"porosity":1})
# the restart test has to run after the test that produces the restart file
set_tests_properties(test_restartco2 PROPERTIES DEPENDS test_boxco2)
# set_tests_properties(test_restartco2 PROPERTIES DEPENDS test_co2_box)
# build target for the CO2 non-isothermal test problem
add_dumux_test(test_boxco2ni test_boxco2ni test_boxco2ni.cc
python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
--script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/co2nibox-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/heterogeneousboxni-00019.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_boxco2ni"
--zeroThreshold {"porosity":1})
add_dumux_test(test_ccco2ni test_ccco2ni test_ccco2ni.cc
python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
--script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/co2nicc-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/heterogeneousccni-00018.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_ccco2ni")
dune_add_test(NAME test_co2ni_box
SOURCES test_co2_fv.cc
COMPILE_DEFINITIONS TYPETAG=HeterogeneousNIBoxTypeTag
COMPILE_DEFINITIONS ISOTHERMAL=0
CMAKE_GUARD "( dune-alugrid_FOUND AND DUNE_GRID_EXPERIMENTAL_GRID_EXTENSIONS )"
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/co2nibox-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/co2ni_box-00019.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_co2ni_box test_co2ni_fv.input -Problem.Name co2ni_box"
--zeroThreshold {"porosity":1})
dune_add_test(NAME test_co2ni_tpfa
SOURCES test_co2_fv.cc
COMPILE_DEFINITIONS TYPETAG=HeterogeneousNICCTpfaTypeTag
COMPILE_DEFINITIONS ISOTHERMAL=0
CMAKE_GUARD "( dune-alugrid_FOUND AND DUNE_GRID_EXPERIMENTAL_GRID_EXTENSIONS )"
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/co2nicc-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/co2ni_tpfa-00018.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_co2ni_tpfa test_co2ni_fv.input -Problem.Name co2ni_tpfa")
#install sources
install(FILES
......@@ -55,8 +69,5 @@ heterogeneousco2tables.hh
heterogeneousproblem.hh
heterogeneousproblemni.hh
heterogeneousspatialparameters.hh
test_boxco2.cc
test_boxco2ni.cc
test_ccco2.cc
test_ccco2ni.cc
test_co2_fv.cc
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/test/implicit/co2)
......@@ -25,17 +25,16 @@
#ifndef DUMUX_HETEROGENEOUS_CO2TABLES_HH
#define DUMUX_HETEROGENEOUS_CO2TABLES_HH
#include <cassert>
#include <dumux/material/components/co2tablereader.hh>
namespace Dumux
{
namespace HeterogeneousCO2Tables
{
namespace Dumux {
namespace HeterogeneousCO2Tables {
// the real work is done by some external program which provides
// ready-to-use tables.
#include "co2values.inc"
}
}
} // end namespace HeterogeneousCO2Tables
} // end namespace Dumux
#endif
......@@ -54,11 +54,10 @@ SET_PROP(HeterogeneousSpatialParams, MaterialLaw)
private:
// define the material law which is parameterized by effective
// saturations
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef RegularizedBrooksCorey<Scalar> EffMaterialLaw;
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
public:
// define the material law parameterized by absolute saturations
typedef EffToAbsLaw<EffMaterialLaw> type;
typedef EffToAbsLaw<RegularizedBrooksCorey<Scalar>> type;
};
}
......@@ -72,37 +71,33 @@ public:
template<class TypeTag>
class HeterogeneousSpatialParams : public ImplicitSpatialParams<TypeTag>
{
typedef ImplicitSpatialParams<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename GET_PROP_TYPE(TypeTag, GridCreator) GridCreator;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
enum { dim=GridView::dimension };
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GridView::template Codim<0>::Entity Element;
using ParentType = ImplicitSpatialParams<TypeTag>;
using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
enum { dimWorld = GridView::dimensionworld };
using GlobalPosition = Dune::FieldVector<typename GridView::ctype, dimWorld>;
public:
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
typedef typename MaterialLaw::Params MaterialLawParams;
using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
using MaterialLawParams = typename MaterialLaw::Params;
using PermeabilityType = Scalar;
/*!
* \brief The constructor
*
* \param gridView The grid view
*/
HeterogeneousSpatialParams(const GridView &gridView)
: ParentType(gridView), gridView_(gridView)
HeterogeneousSpatialParams(const Problem& problem)
: ParentType(problem)
{
/*
* Layer Index Setup:
*/
barrierTop_ = 1;
barrierMiddle_ = 2;
reservoir_ = 3;
// heat conductivity of granite
lambdaSolid_ = 2.8;
......@@ -123,40 +118,48 @@ public:
materialParams_.setPe(1e4);
}
~HeterogeneousSpatialParams()
{}
/*!
* \brief Reads layer information from the grid
*
*/
void setParams()
void getParamsFromGrid()
{
int numElements = gridView_.size(0);
paramIdx_.resize(numElements);
const auto& gridView = this->problem().fvGridGeometry().gridView();
paramIdx_.resize(gridView.size(0));
for (const auto& element : elements(gridView_))
for (const auto& element : elements(gridView))
{
int eIdx = gridView_.indexSet().index(element);
int param = GridCreator::parameters(element)[0];
paramIdx_[eIdx] = param;
const auto eIdx = this->problem().fvGridGeometry().elementMapper().index(element);
paramIdx_[eIdx] = GridCreator::parameters(element)[0];
}
}
/*!
* \brief Returns the scalar intrinsic permeability \f$[m^2]\f$
* \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$
* \note It is possibly solution dependent.
*
* \param element The finite element
* \param fvGeometry The finite volume geometry of the element
* \param scvIdx The local index of the sub-control volume
* \param element The current element
* \param scv The sub-control volume inside the element.
* \param elemSol The solution at the dofs connected to the element.
* \return instrinsic permeability
*/
const Scalar intrinsicPermeability(const Element &element,
const FVElementGeometry &fvGeometry,
int scvIdx) const
PermeabilityType permeability(const Element& element,
const SubControlVolume& scv,
const ElementSolutionVector& elemSol) const
{
//Get the global index of the element
int eIdx = gridView_.indexSet().index(element);
// Get the global index of the element
const auto eIdx = this->problem().fvGridGeometry().elementMapper().index(element);
return permeability(eIdx);
}
/*!
* \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$
* \note It is possibly solution dependent.
*
* \param eIdx the element index
*/
PermeabilityType permeability(std::size_t eIdx) const
{
if (paramIdx_[eIdx] == barrierTop_)
return barrierTopK_;
else if (paramIdx_[eIdx] == barrierMiddle_)
......@@ -168,36 +171,44 @@ public:
/*!
* \brief Returns the porosity \f$[-]\f$
*
* \param element The finite element
* \param fvGeometry The finite volume geometry of the element
* \param scvIdx The local index of the sub-control volume
* \param element The current element
* \param scv The sub-control volume inside the element.
* \param elemSol The solution at the dofs connected to the element.
* \return porosity
*/
Scalar porosity(const Element &element,
const FVElementGeometry &fvGeometry,
int scvIdx) const
Scalar porosity(const Element& element,
const SubControlVolume& scv,
const ElementSolutionVector& elemSol) const
{
//Get the global index of the element
int eIdx = gridView_.indexSet().index(element);
// Get the global index of the element
const auto eIdx = this->problem().fvGridGeometry().elementMapper().index(element);
return porosity(eIdx);
}
/*!
* \brief Returns the porosity \f$[-]\f$
*
* \param eIdx The element index
*/
Scalar porosity(std::size_t eIdx) const
{
if (paramIdx_[eIdx] == barrierTop_)
return barrierTopPorosity_;
else if (paramIdx_[eIdx] == barrierMiddle_)
return barrierMiddlePorosity_;
else
return reservoirPorosity_;
}
/*!
* \brief Returns the parameter object for the Brooks-Corey material law
* \brief Function for defining the parameters needed by constitutive relationships (kr-sw, pc-sw, etc.).
*
* \param element The finite element
* \param fvGeometry The finite volume geometry of the element
* \param scvIdx The local index of the sub-control volume
* \return the material parameters object
* \param globalPos The position of the center of the element
*/
const MaterialLawParams& materialLawParams(const Element &element,
const FVElementGeometry &fvGeometry,
int scvIdx) const
const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const
{
return materialParams_;
}
......@@ -207,13 +218,9 @@ public:
*
* This is only required for non-isothermal models.
*
* \param element The finite element
* \param fvGeometry The finite volume geometry
* \param scvIdx The local index of the sub-control volume
* \param globalPos The position of the center of the element
*/
Scalar solidHeatCapacity(const Element &element,
const FVElementGeometry &fvGeometry,
const int scvIdx) const
Scalar solidHeatCapacityAtPos(const GlobalPosition& globalPos) const
{
return 790; // specific heat capacity of granite [J / (kg K)]
}
......@@ -223,41 +230,27 @@ public:
*
* This is only required for non-isothermal models.
*
* \param element The finite element
* \param fvGeometry The finite volume geometry
* \param scvIdx The local index of the sub-control volume
* \param globalPos The position of the center of the element
*/
Scalar solidDensity(const Element &element,
const FVElementGeometry &fvGeometry,
const int scvIdx) const
Scalar solidDensityAtPos(const GlobalPosition& globalPos) const
{
return 2700; // density of granite [kg/m^3]
}
/*!
* \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the solid
*
* This is only required for non-isothermal models.
* \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the porous material.
*
* \param element The finite element
* \param fvGeometry The finite volume geometry of the element
* \param scvIdx The local index of the sub-control volume
* \param globalPos The position of the center of the element
*/
Scalar solidThermalConductivity(const Element &element,
const FVElementGeometry &fvGeometry,
const int scvIdx) const
Scalar solidThermalConductivityAtPos(const GlobalPosition& globalPos) const
{
return lambdaSolid_;
}
private:
int barrierTop_;
int barrierMiddle_;
int reservoir_;
int barrierTop_ = 1;
int barrierMiddle_ = 2;
int reservoir_ = 3;
Scalar barrierTopPorosity_;
Scalar barrierMiddlePorosity_;
......@@ -269,11 +262,9 @@ private:
Scalar lambdaSolid_;
MaterialLawParams materialParams_;
const GridView gridView_;
std::vector<int> paramIdx_;
};
}
} // 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
*
* \brief Test for the CO2 box model.
*/
#include <config.h>
#include "heterogeneousproblem.hh"
#include <dumux/common/start.hh>
/*!
* \brief Provides an interface for customizing error messages associated with
* reading in parameters. This test only works if the flags for experimental
* grid extensions and ALUGrid are set. See the commented part in optim.opts and
* debug.opts for reference.
*
* \param progName The name of the program, that was tried to be started.
* \param errorMsg The error message that was issued by the start function.
* Comprises the thing that went wrong and a general help message.
*/
void usage(const char *progName, const std::string &errorMsg)
{
if (errorMsg.size() > 0) {
std::string errorMessageOut = "\nUsage: ";
errorMessageOut += progName;
errorMessageOut += " [options]\n";
errorMessageOut += errorMsg;
errorMessageOut += "\n\nThe List of Mandatory arguments for this program is:\n"
"\t-TimeManager.DtInitial The initial timestep size. [s] \n"
"\t-TimeManager.TEnd The end of the simulation. [s] \n"
"\t-Grid.GridFile The file name of the file containing the grid \n"
"\t definition in DGF format\n"
"\t-FluidSystem.nTemperature Number of tabularization entries [-] \n"
"\t-FluidSystem.nPressure Number of tabularization entries [-] \n"
"\t-FluidSystem.pressureLow Low end for tabularization of fluid properties [Pa] \n"
"\t-FluidSystem.pressureHigh High end for tabularization of fluid properties [Pa] \n"
"\t-FluidSystem.temperatureLow Low end for tabularization of fluid properties [Pa] \n"
"\t-FluidSystem.temperatureHigh High end for tabularization of fluid properties [Pa] \n"
"\t-Problem.name The name of the output files [-] \n"
"\t-Problem.depthBOR Depth below ground surface [m] \n"
"\t-Problem.injectionRate Injection rate in kg/(s*m^2) or mole/(m^2*s) \n"
"\t-LinearSolver.ResidualReduction [-] \n";
std::cout << errorMessageOut
<< "\n";
}
}
int main(int argc, char** argv)
{
#if !DUNE_GRID_EXPERIMENTAL_GRID_EXTENSIONS
#warning Test needs experimental grid extensions.
std::cerr << "Test skipped, it needs experimental grid extensions, see optim.opts or debug.opts." << std::endl;
return 77;
#elif !HAVE_DUNE_ALUGRID
#warning You need to have dune-ALUGrid installed to run this test.
std::cerr << "You need to have dune-ALUGrid installed to run this test." << std::endl;
return 77;
#else
typedef TTAG(HeterogeneousBoxProblem) ProblemTypeTag;
return Dumux::start<ProblemTypeTag>(argc, argv, usage);
#endif
}
[TimeManager]
DtInitial = 250# [s]
TEnd = 1e6# [s]
[Grid]
File = ./grids/heterogeneousSmall.dgf # relative path to the grid file
[FluidSystem]
NTemperature = 100 # [-] number of tabularization entries
NPressure = 100 # [-] number of tabularization entries
PressureLow = 1e5# [Pa] low end for tabularization of fluid properties
PressureHigh = 3e7# [Pa] high end for tabularization of fluid properties
TemperatureLow = 290.00 # [Pa] low end for tabularization of fluid properties
TemperatureHigh = 331.00 # [Pa] high end for tabularization of fluid properties
[Problem]
Name = heterogeneousbox # [-] the name of the output files
EnableGravity = TRUE
DepthBOR = 1200# [m] depth below ground surface
InjectionRate = 1e-4 # [kg/sq/s]
[Vtk]
AddVelocity = FALSE
[LinearSolver]
ResidualReduction = 1e-10
// -*- 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 CO2 non-isothermal box model.
*/
#include <config.h>
#include "heterogeneousproblemni.hh"
#include <dumux/common/start.hh>
/*!
* \brief Provides an interface for customizing error messages associated with
* reading in parameters. This test only works if the flags for experimental
* grid extensions and ALUGrid are set. See the commented part in optim.opts and
* debug.opts for reference.
*
* \param progName The name of the program, that was tried to be started.
* \param errorMsg The error message that was issued by the start function.
* Comprises the thing that went wrong and a general help message.
*/
void usage(const char *progName, const std::string &errorMsg)
{
if (errorMsg.size() > 0) {
std::string errorMessageOut = "\nUsage: ";
errorMessageOut += progName;
errorMessageOut += " [options]\n";
errorMessageOut += errorMsg;
errorMessageOut += "\n\nThe List of Mandatory arguments for this program is:\n"
"\t-TimeManager.DtInitial The initial timestep size. [s] \n"
"\t-TimeManager.TEnd The end of the simulation. [s] \n"
"\t-Grid.GridFile The file name of the file containing the grid \n"
"\t definition in DGF format\n"
"\t-FluidSystem.nTemperature Number of tabularization entries [-] \n"
"\t-FluidSystem.nPressure Number of tabularization entries [-] \n"
"\t-FluidSystem.pressureLow Low end for tabularization of fluid properties [Pa] \n"
"\t-FluidSystem.pressureHigh High end for tabularization of fluid properties [Pa] \n"
"\t-FluidSystem.temperatureLow Low end for tabularization of fluid properties [Pa] \n"
"\t-FluidSystem.temperatureHigh High end for tabularization of fluid properties [Pa] \n"
"\t-Problem.name The name of the output files [-] \n"
"\t-Problem.depthBOR Depth below ground surface [m] \n"
"\t-Problem.injectionRate Injection rate in kg/(s*m^2) or mole/(m^2*s) \n"
"\t-Problem.injectionPressure Injection pressure [Pa] \n"
"\t-Problem.injectionTemperature Injection temperature in [K] \n"
"\t-LinearSolver.ResidualReduction [-] \n";