diff --git a/lecture/mm/heavyoil/sagd/problem.hh b/lecture/mm/heavyoil/sagd/problem.hh index 0bb354a4beb923ce9ebf148c8cde9a72b5225433..c47acf8967dd5e28b6b5a1abd8afa13f1a170b8b 100644 --- a/lecture/mm/heavyoil/sagd/problem.hh +++ b/lecture/mm/heavyoil/sagd/problem.hh @@ -28,7 +28,7 @@ #include -#include +#include #include #include @@ -48,36 +48,50 @@ template class SagdProblem; namespace Properties { -NEW_TYPE_TAG(SagdTypeTag, INHERITS_FROM(ThreePWaterOilNI)); -NEW_TYPE_TAG(ThreePWaterOilSagdBoxTypeTag, INHERITS_FROM(BoxModel, SagdTypeTag)); +// Create new type tags +namespace TTag { +struct SagdTypeTag { using InheritsFrom = std::tuple; }; +struct ThreePWaterOilSagdBoxTypeTag { using InheritsFrom = std::tuple; }; +} // end namespace TTag // Set the grid type -SET_TYPE_PROP(SagdTypeTag, Grid, Dune::YaspGrid<2>); +template +struct Grid { using type = Dune::YaspGrid<2>; }; // Set the problem property -SET_TYPE_PROP(SagdTypeTag, Problem, Dumux::SagdProblem); +template +struct Problem { using type = Dumux::SagdProblem; }; // Set the spatial parameters -SET_PROP(SagdTypeTag, SpatialParams) +template +struct SpatialParams { - using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); - using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using FVGridGeometry = GetPropType; + using Scalar = GetPropType; using type = SagdSpatialParams; }; -// Set the fluid system -SET_TYPE_PROP(SagdTypeTag, - FluidSystem, - Dumux::FluidSystems::H2OHeavyOil); +// Set fluid configuration +template +struct FluidSystem +{ +private: + using Scalar = GetPropType; +public: + using type = FluidSystems::H2OHeavyOil; +}; -SET_BOOL_PROP(SagdTypeTag, OnlyGasPhaseCanDisappear, true); +template +struct OnlyGasPhaseCanDisappear { static constexpr bool value = true; }; -SET_BOOL_PROP(SagdTypeTag, UseMoles, true); +template +struct UseMoles { static constexpr bool value = true; }; // Set the fluid system -SET_PROP(SagdTypeTag, SolidSystem) +template +struct SolidSystem { - using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using Scalar = GetPropType; using InertComponent = Components::Constant<1, Scalar>; using type = SolidSystems::InertSolidPhase; }; @@ -96,11 +110,11 @@ template class SagdProblem : public PorousMediumFlowProblem { using ParentType = PorousMediumFlowProblem; - using GridView = typename GET_PROP_TYPE(TypeTag, GridView); - using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); - using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); - using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem); - using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices; + using GridView = GetPropType; + using FVGridGeometry = GetPropType; + using Scalar = GetPropType; + using FluidSystem = GetPropType; + using Indices = typename GetPropType::Indices; enum { pressureIdx = Indices::pressureIdx, switch1Idx = Indices::switch1Idx, @@ -121,12 +135,12 @@ class SagdProblem : public PorousMediumFlowProblem dimWorld = GridView::dimensionworld }; - using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables); - using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector); - using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView; - using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes); + using PrimaryVariables = GetPropType; + using NumEqVector = GetPropType; + using ElementVolumeVariables = typename GetPropType::LocalView; + using BoundaryTypes = GetPropType; using Element = typename GridView::template Codim<0>::Entity; - using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView; + using FVElementGeometry = typename GetPropType::LocalView; using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; using GlobalPosition = typename SubControlVolumeFace::GlobalPosition; diff --git a/lecture/mm/heavyoil/sagd/sagd.cc b/lecture/mm/heavyoil/sagd/sagd.cc index d30ad927aa17c316c6d6bfe27d3ce51ab211feeb..f62e0cf75d5455e4654030d7e41ad633f8ebc2ac 100644 --- a/lecture/mm/heavyoil/sagd/sagd.cc +++ b/lecture/mm/heavyoil/sagd/sagd.cc @@ -79,7 +79,7 @@ int main(int argc, char** argv) try using namespace Dumux; // define the type tag for this problem - using TypeTag = TTAG(ThreePWaterOilSagdBoxTypeTag); + using TypeTag = Properties::TTag::ThreePWaterOilSagdBoxTypeTag; //////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////// @@ -95,7 +95,7 @@ int main(int argc, char** argv) try Parameters::init(argc, argv, usage); // try to create a grid (from the given grid file or the input file) - GridManager gridManager; + GridManager> gridManager; gridManager.init(); //////////////////////////////////////////////////////////// @@ -106,37 +106,37 @@ int main(int argc, char** argv) try const auto& leafGridView = gridManager.grid().leafGridView(); // create the finite volume grid geometry - using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); + using FVGridGeometry = GetPropType; auto fvGridGeometry = std::make_shared(leafGridView); fvGridGeometry->update(); // the problem (initial and boundary conditions) - using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + using Problem = GetPropType; auto problem = std::make_shared(fvGridGeometry); // the solution vector - using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector); + using SolutionVector = GetPropType; SolutionVector x(fvGridGeometry->numDofs()); problem->applyInitialSolution(x); auto xOld = x; // the grid variables - using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables); + using GridVariables = GetPropType; auto gridVariables = std::make_shared(problem, fvGridGeometry); - gridVariables->init(x, xOld); + gridVariables->init(x); // get some time loop parameters - using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using Scalar = GetPropType; const auto tEnd = getParam("TimeLoop.TEnd"); const auto maxDt = getParam("TimeLoop.MaxTimeStepSize"); auto dt = getParam("TimeLoop.DtInitial"); // intialize the vtk output module - using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields); + using IOFields = GetPropType; VtkOutputModule vtkWriter(*gridVariables, x, problem->name()); - using VelocityOutput = typename GET_PROP_TYPE(TypeTag, VelocityOutput); + using VelocityOutput = GetPropType; vtkWriter.addVelocityOutput(std::make_shared(*gridVariables)); - VtkOutputFields::init(vtkWriter); //! Add model specific output fields + IOFields::initOutputModule(vtkWriter); //!< Add model specific output fields vtkWriter.write(0.0); // instantiate time loop @@ -165,7 +165,7 @@ int main(int argc, char** argv) try // solve the non-linear system with time step control nonLinearSolver.solve(x, *timeLoop); - using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables); + using PrimaryVariables = GetPropType; PrimaryVariables storage(0); const auto& localResidual = assembler->localResidual(); for (const auto& element : elements(leafGridView, Dune::Partitions::interior)) diff --git a/lecture/mm/heavyoil/sagdcyclic/CMakeLists.txt b/lecture/mm/heavyoil/sagdcyclic/CMakeLists.txt index 6733f9d8e7b41bda679c27e258dcac7145246225..f296fb32e4f31fe64268a4a5937939d28fd7a050 100644 --- a/lecture/mm/heavyoil/sagdcyclic/CMakeLists.txt +++ b/lecture/mm/heavyoil/sagdcyclic/CMakeLists.txt @@ -1,7 +1,7 @@ dune_add_test(NAME sagd_cyclic SOURCES sagd_cyclic.cc COMMAND ${CMAKE_CURRENT_BINARY_DIR}/sagd_cyclic - CMD_ARGS -TimeLoop.TEnd 43200 -TimeLoop.MaxTimeStepSize21600) + CMD_ARGS -TimeLoop.TEnd 43200 -TimeLoop.MaxTimeStepSize 21600) install(FILES problem.hh diff --git a/lecture/mm/heavyoil/sagdcyclic/problem.hh b/lecture/mm/heavyoil/sagdcyclic/problem.hh index 6e25f04b144a9342fd79b91622601eef61c1b311..ec475896da046f40850e296c7674596e5b4a3fe2 100644 --- a/lecture/mm/heavyoil/sagdcyclic/problem.hh +++ b/lecture/mm/heavyoil/sagdcyclic/problem.hh @@ -28,7 +28,7 @@ #include -#include +#include #include #include @@ -48,36 +48,50 @@ template class SagdCyclicProblem; namespace Properties { -NEW_TYPE_TAG(SagdTypeTag, INHERITS_FROM(ThreePWaterOilNI)); -NEW_TYPE_TAG(ThreePWaterOilSagdCyclicBoxTypeTag, INHERITS_FROM(BoxModel, SagdTypeTag)); +// Create new type tags +namespace TTag { +struct SagdTypeTag { using InheritsFrom = std::tuple; }; +struct ThreePWaterOilSagdCyclicBoxTypeTag { using InheritsFrom = std::tuple; }; +} // end namespace TTag // Set the grid type -SET_TYPE_PROP(SagdTypeTag, Grid, Dune::YaspGrid<2>); +template +struct Grid { using type = Dune::YaspGrid<2>; }; // Set the problem property -SET_TYPE_PROP(SagdTypeTag, Problem, Dumux::SagdCyclicProblem); +template +struct Problem { using type = Dumux::SagdCyclicProblem; }; // Set the spatial parameters -SET_PROP(SagdTypeTag, SpatialParams) +template +struct SpatialParams { - using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); - using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using FVGridGeometry = GetPropType; + using Scalar = GetPropType; using type = SagdSpatialParams; }; -// Set the fluid system -SET_TYPE_PROP(SagdTypeTag, - FluidSystem, - Dumux::FluidSystems::H2OHeavyOil); +// Set fluid configuration +template +struct FluidSystem +{ +private: + using Scalar = GetPropType; +public: + using type = FluidSystems::H2OHeavyOil; +}; -SET_BOOL_PROP(SagdTypeTag, OnlyGasPhaseCanDisappear, true); +template +struct OnlyGasPhaseCanDisappear { static constexpr bool value = true; }; -SET_BOOL_PROP(SagdTypeTag, UseMoles, true); +template +struct UseMoles { static constexpr bool value = true; }; // Set the fluid system -SET_PROP(SagdTypeTag, SolidSystem) +template +struct SolidSystem { - using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using Scalar = GetPropType; using InertComponent = Components::Constant<1, Scalar>; using type = SolidSystems::InertSolidPhase; }; @@ -96,11 +110,11 @@ template class SagdCyclicProblem : public PorousMediumFlowProblem { using ParentType = PorousMediumFlowProblem; - using GridView = typename GET_PROP_TYPE(TypeTag, GridView); - using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); - using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); - using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem); - using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices; + using GridView = GetPropType; + using FVGridGeometry = GetPropType; + using Scalar = GetPropType; + using FluidSystem = GetPropType; + using Indices = typename GetPropType::Indices; enum { pressureIdx = Indices::pressureIdx, switch1Idx = Indices::switch1Idx, @@ -121,14 +135,14 @@ class SagdCyclicProblem : public PorousMediumFlowProblem dimWorld = GridView::dimensionworld }; - using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables); - using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector); - using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables); - using Problem = typename GET_PROP_TYPE(TypeTag, Problem); - using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView; - using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes); + using PrimaryVariables = GetPropType; + using NumEqVector = GetPropType; + using VolumeVariables = GetPropType; + using Problem = GetPropType; + using ElementVolumeVariables = typename GetPropType::LocalView; + using BoundaryTypes = GetPropType; using Element = typename GridView::template Codim<0>::Entity; - using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView; + using FVElementGeometry = typename GetPropType::LocalView; using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; using GlobalPosition = typename SubControlVolumeFace::GlobalPosition; diff --git a/lecture/mm/heavyoil/sagdcyclic/sagd_cyclic.cc b/lecture/mm/heavyoil/sagdcyclic/sagd_cyclic.cc index 427279c07b9e6571b6a2d1e9c8fd8e8552fa97c7..cadceda02c84044978b0960f6e1b772ef69839f9 100644 --- a/lecture/mm/heavyoil/sagdcyclic/sagd_cyclic.cc +++ b/lecture/mm/heavyoil/sagdcyclic/sagd_cyclic.cc @@ -79,7 +79,7 @@ int main(int argc, char** argv) try using namespace Dumux; // define the type tag for this problem - using TypeTag = TTAG(ThreePWaterOilSagdCyclicBoxTypeTag); + using TypeTag = Properties::TTag::ThreePWaterOilSagdCyclicBoxTypeTag; //////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////// @@ -95,7 +95,7 @@ int main(int argc, char** argv) try Parameters::init(argc, argv, usage); // try to create a grid (from the given grid file or the input file) - GridManager gridManager; + GridManager> gridManager; gridManager.init(); //////////////////////////////////////////////////////////// @@ -106,37 +106,37 @@ int main(int argc, char** argv) try const auto& leafGridView = gridManager.grid().leafGridView(); // create the finite volume grid geometry - using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); + using FVGridGeometry = GetPropType; auto fvGridGeometry = std::make_shared(leafGridView); fvGridGeometry->update(); // the problem (initial and boundary conditions) - using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + using Problem = GetPropType; auto problem = std::make_shared(fvGridGeometry); // the solution vector - using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector); + using SolutionVector = GetPropType; SolutionVector x(fvGridGeometry->numDofs()); problem->applyInitialSolution(x); auto xOld = x; // the grid variables - using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables); + using GridVariables = GetPropType; auto gridVariables = std::make_shared(problem, fvGridGeometry); - gridVariables->init(x, xOld); + gridVariables->init(x); // get some time loop parameters - using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using Scalar = GetPropType; const auto tEnd = getParam("TimeLoop.TEnd"); const auto maxDt = getParam("TimeLoop.MaxTimeStepSize"); auto dt = getParam("TimeLoop.DtInitial"); // intialize the vtk output module - using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields); + using IOFields = GetPropType; VtkOutputModule vtkWriter(*gridVariables, x, problem->name()); - using VelocityOutput = typename GET_PROP_TYPE(TypeTag, VelocityOutput); + using VelocityOutput = GetPropType; vtkWriter.addVelocityOutput(std::make_shared(*gridVariables)); - VtkOutputFields::init(vtkWriter); //! Add model specific output fields + IOFields::initOutputModule(vtkWriter); //!< Add model specific output fields vtkWriter.write(0.0); // instantiate time loop @@ -167,7 +167,7 @@ int main(int argc, char** argv) try // solve the non-linear system with time step control nonLinearSolver.solve(x, *timeLoop); - using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables); + using PrimaryVariables = GetPropType; PrimaryVariables storage(0); const auto& localResidual = assembler->localResidual(); for (const auto& element : elements(leafGridView, Dune::Partitions::interior)) diff --git a/lecture/mm/heavyoil/sagdcyclichyst/problem.hh b/lecture/mm/heavyoil/sagdcyclichyst/problem.hh index 0428c7adfd8faea956957df6f563b8b2dbfcb83f..a7ad02e4f011ec12f401db6c8f0645bb2e481361 100644 --- a/lecture/mm/heavyoil/sagdcyclichyst/problem.hh +++ b/lecture/mm/heavyoil/sagdcyclichyst/problem.hh @@ -28,7 +28,7 @@ #include -#include +#include #include #include @@ -48,36 +48,50 @@ template class SagdCyclicHystProblem; namespace Properties { -NEW_TYPE_TAG(SagdTypeTag, INHERITS_FROM(ThreePWaterOilNI)); -NEW_TYPE_TAG(ThreePWaterOilSagdCyclicHystBoxTypeTag, INHERITS_FROM(BoxModel, SagdTypeTag)); +// Create new type tags +namespace TTag { +struct SagdTypeTag { using InheritsFrom = std::tuple; }; +struct ThreePWaterOilSagdCyclicHystBoxTypeTag { using InheritsFrom = std::tuple; }; +} // end namespace TTag // Set the grid type -SET_TYPE_PROP(SagdTypeTag, Grid, Dune::YaspGrid<2>); +template +struct Grid { using type = Dune::YaspGrid<2>; }; // Set the problem property -SET_TYPE_PROP(SagdTypeTag, Problem, Dumux::SagdCyclicHystProblem); +template +struct Problem { using type = Dumux::SagdCyclicHystProblem; }; // Set the spatial parameters -SET_PROP(SagdTypeTag, SpatialParams) +template +struct SpatialParams { - using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); - using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using FVGridGeometry = GetPropType; + using Scalar = GetPropType; using type = SagdSpatialParams; }; -// Set the fluid system -SET_TYPE_PROP(SagdTypeTag, - FluidSystem, - Dumux::FluidSystems::H2OHeavyOil); +// Set fluid configuration +template +struct FluidSystem +{ +private: + using Scalar = GetPropType; +public: + using type = FluidSystems::H2OHeavyOil; +}; -SET_BOOL_PROP(SagdTypeTag, OnlyGasPhaseCanDisappear, true); +template +struct OnlyGasPhaseCanDisappear { static constexpr bool value = true; }; -SET_BOOL_PROP(SagdTypeTag, UseMoles, true); +template +struct UseMoles { static constexpr bool value = true; }; // Set the fluid system -SET_PROP(SagdTypeTag, SolidSystem) +template +struct SolidSystem { - using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using Scalar = GetPropType; using InertComponent = Components::Constant<1, Scalar>; using type = SolidSystems::InertSolidPhase; }; @@ -96,11 +110,11 @@ template class SagdCyclicHystProblem : public PorousMediumFlowProblem { using ParentType = PorousMediumFlowProblem; - using GridView = typename GET_PROP_TYPE(TypeTag, GridView); - using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); - using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); - using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem); - using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices; + using GridView = GetPropType; + using FVGridGeometry = GetPropType; + using Scalar = GetPropType; + using FluidSystem = GetPropType; + using Indices = typename GetPropType::Indices; enum { pressureIdx = Indices::pressureIdx, switch1Idx = Indices::switch1Idx, @@ -121,14 +135,14 @@ class SagdCyclicHystProblem : public PorousMediumFlowProblem dimWorld = GridView::dimensionworld }; - using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables); - using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector); - using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables); - using Problem = typename GET_PROP_TYPE(TypeTag, Problem); - using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView; - using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes); + using PrimaryVariables = GetPropType; + using NumEqVector = GetPropType; + using VolumeVariables = GetPropType; + using Problem = GetPropType; + using ElementVolumeVariables = typename GetPropType::LocalView; + using BoundaryTypes = GetPropType; using Element = typename GridView::template Codim<0>::Entity; - using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView; + using FVElementGeometry = typename GetPropType::LocalView; using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; using GlobalPosition = typename SubControlVolumeFace::GlobalPosition; diff --git a/lecture/mm/heavyoil/sagdcyclichyst/problem.hh.SAVE.hh b/lecture/mm/heavyoil/sagdcyclichyst/problem.hh.SAVE.hh deleted file mode 100644 index c8de471acf095746892d400a176075c1eda98d59..0000000000000000000000000000000000000000 --- a/lecture/mm/heavyoil/sagdcyclichyst/problem.hh.SAVE.hh +++ /dev/null @@ -1,467 +0,0 @@ -// -*- 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 vesion. * - * * - * 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 . * - *****************************************************************************/ -/*! - * \file - * - * \brief Non-isothermal SAGD problem - * - */ -#ifndef DUMUX_SAGDCYCLICHYSTPROBLEM_HH -#define DUMUX_SAGDCYCLICHYSTPROBLEM_HH - -#include -#include -#include -#include -#include -#include -#include "spatialparams.hh" - -#define ISOTHERMAL 0 - -namespace Dumux { - -template -class SagdCyclicHystProblem; - -namespace Properties { - -NEW_TYPE_TAG(SagdCyclicHystProblemTypeTag, INHERITS_FROM(ThreePWaterOilNI, SagdCyclicHystSpatialParamsTypeTag)); - -NEW_TYPE_TAG(SagdCyclicHystBoxProblem, INHERITS_FROM(BoxModel, SagdCyclicHystProblemTypeTag)); - -NEW_TYPE_TAG(SagdCyclicHystCCProblem, INHERITS_FROM(CCModel, SagdCyclicHystProblemTypeTag)); - -// Set the grid type -SET_TYPE_PROP(SagdCyclicHystProblemTypeTag, Grid, Dune::YaspGrid<2>); - -// Set the problem property -SET_TYPE_PROP(SagdCyclicHystProblemTypeTag, Problem, Dumux::SagdCyclicHystProblem); - -// Set the fluid system -SET_TYPE_PROP(SagdCyclicHystProblemTypeTag, - FluidSystem, - Dumux::FluidSystems::H2OHeavyOil); - -// Enable gravity -SET_BOOL_PROP(SagdCyclicHystProblemTypeTag, ProblemEnableGravity, true); - -// Use forward differences instead of central differences -SET_INT_PROP(SagdCyclicHystProblemTypeTag, ImplicitNumericDifferenceMethod, +1); - -// Write newton convergence -SET_BOOL_PROP(SagdCyclicHystProblemTypeTag, NewtonWriteConvergence, false); - -SET_BOOL_PROP(SagdCyclicHystProblemTypeTag, UseSimpleModel, true); - -SET_BOOL_PROP(SagdCyclicHystProblemTypeTag, UseMoles, true); - -} // end namespace Properties - -/*! - * \ingroup ThreePWaterOilBoxModel - * \ingroup ImplicitTestProblems - * \brief Non-isothermal problem where ... - * - * This problem uses the \ref ThreePWaterOilModel. - * - * */ -template -class SagdCyclicHystProblem : public PorousMediumFlowProblem -{ - using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar) - using GridView = typename GET_PROP_TYPE(TypeTag, GridView); - using Grid = typename GridView::Grid; - using ParentType = PorousMediumFlowProblem; - // copy some indices for convenience - using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices; - static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases); - static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents); - static constexpr int pressureIdx = Indices::pressureIdx; - static constexpr int switch1Idx = Indices::switch1Idx; - static constexpr int switch2Idx = Indices::switch2Idx; - static constexpr int energyEqIdx = Indices::energyEqIdx; - // phase and component indices - static constexpr int wPhaseIdx = Indices::Phase0Idx; - static constexpr int nPhaseIdx = Indices::Phase1Idx; - static constexpr int gPhaseIdx = Indices::Phase2Idx; - static constexpr int wCompIdx = Indices::Comp0Idx; - static constexpr int nCompIdx = Indices::Comp1Idx; - // Phase State - static constexpr int wPhaseOnly = Indices::wPhaseOnly; - static constexpr int wnPhaseOnly = Indices::wnPhaseOnly; - static constexpr int wgPhaseOnly = Indices::wgPhaseOnly; - static constexpr int threePhases = Indices::threePhases; - // Grid and world dimension - static constexpr int dim = GridView::dimension; - static constexpr int dimWorld = GridView::dimensionworld; - using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables); - using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables); - using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes); - using Element = typename GridView::template Codim<0>::Entity; - using Vertex = typename GridView::template Codim::Entity; - using Intersection = typename GridView::Intersection; - using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView; - using FluidSystem = typedef typename GET_PROP_TYPE(TypeTag, FluidSystem); - using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView; - using ElementIterator = typename GridView::template Codim<0>::Iterator; - using VertexIterator = typename GridView::template Codim::Iterator; - using GlobalPosition = typename Element::Geometry::GlobalCoordinate; - static constexpr bool isBox = FVGridGeometry::discMethod == DiscretizationMethod::box; - static constexpr int dofCodim = isBox ? dim : 0; - static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles); - -public: - /*! - * \brief The constructor - * - * \param timeLoop The time manager - * \param gridView The grid view - */ - SagdCyclicHystProblem(Tstd::shared_ptr fvGridGeometry) - : ParentType(fvGridGeometry), pOut_(4.0e6) - { - name_ = getParam("Problem.Name"); - maxDepth_ = 400.0; // [m] - FluidSystem::init(); - totalMassProducedOil_ =0; - totalMassProducedWater_ =0; - - this->timeLoop().startNextEpisode(3600*12.); - - //stateing in the console whether mole or mass fractions are used - if (!useMoles) - std::cout<<"Problem uses mass-fractions."<timeLoop().timeStepIndex() == 0 || - this->timeLoop().timeStepIndex() % 10000 == 0 || - this->timeLoop().episodeWillBeFinished() || - this->timeLoop().willBeFinished(); - } - - bool shouldWriteOutput() const - { - return - this->timeLoop().timeStepIndex() == 0 || - this->timeLoop().timeStepIndex() % 10000 == 0 || - this->timeLoop().episodeWillBeFinished() || - this->timeLoop().willBeFinished(); - } - - void postTimeStep() - { - const Scalar time = this->timeLoop().time(); - const Scalar dt = this->timeLoop().timeStepSize(); - - // Calculate storage terms - PrimaryVariables storage; - this->model().globalStorage(storage); - - // Write mass balance information for rank 0 - if (this->gridView().comm().rank() == 0) - { - std::cout<<"Storage: " << storage << " Time: " << time+dt << std::endl; - massBalance.open ("massBalanceCyclicHyst.txt", std::ios::out | std::ios::app ); - massBalance << " Storage " << storage - << " Time " << time+dt - << std::endl; - massBalance.close(); - } - - - this->spatialParams().getMaxSaturation(*this); - this->spatialParams().trappedSat(*this); - - //mass of Oil - const Scalar newMassProducedOil_ = massProducedOil_; - - totalMassProducedOil_ += newMassProducedOil_; - //mass of Water - Scalar newMassProducedWater_ = massProducedWater_; - - totalMassProducedWater_ += newMassProducedWater_; - - const int timeStepIndex = this->timeLoop().timeStepIndex(); - if (timeStepIndex == 0 || - timeStepIndex % 100 == 0 || //after every 1000000 secs - this->timeLoop().episodeWillBeFinished() || - this->timeLoop().willBeFinished()) - - { - std::cout<<" totalMassProducedOil_ : "<< totalMassProducedOil_ << " Time: " << time+dt << std::endl; - std::cout<<" totalMassProducedWater_ : "<< totalMassProducedWater_ << " Time: " << time+dt << std::endl; - } - } - - void episodeEnd() - { - int indexEpisode = this->timeLoop().episodeIndex(); - this->timeLoop().startNextEpisode(3600*12); //episode length sent to 12 hours - std::cout<<"Episode index is set to: "< 40 - eps_) - bcTypes.setAllNeumann(); - - else if (globalPos[0] > 60 - eps_) - bcTypes.setAllDirichlet(); - - // on left - else - bcTypes.setAllNeumann(); - - return bcTypes; - } - - /*! - * \brief Evaluate the boundary conditions for a dirichlet - * boundary segment. - * - * \param values The dirichlet values for the primary variables - * \param globalPos The position for which the bc type should be evaluated - * - * For this method, the \a values parameter stores primary variables. - */ - PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const - { - PrimaryVariables values; - initial_(values, globalPos); - - return values; - } - - /*! - * \brief Evaluate the boundary conditions for a neumann - * boundary segment. - * - * \param values The neumann values for the conservation equations - * \param element The finite element - * \param fvGeomtry The finite-volume geometry in the box scheme - * \param is The intersection between element and boundary - * \param scvIdx The local vertex index - * \param boundaryFaceIdx The index of the boundary face - * - * For this method, the \a values parameter stores the mass flux - * in normal direction of each phase. Negative values mean influx. - */ - void solDependentNeumann(PrimaryVariables &values, - const Element &element, - const FVElementGeometry &fvGeometry, - const Intersection &is, - const int scvIdx, - const int boundaryFaceIdx, - const ElementVolumeVariables &elemVolVars) const - { - values = 0.0; - const int indexEpisode = this->timeLoop().episodeIndex(); - - GlobalPosition globalPos; - if (isBox) - globalPos = element.geometry().corner(scvIdx); - - else - globalPos = is.geometry().center(); - - if (globalPos[1] > 8.5 + eps_ && globalPos[1] < 9.5 - eps_ ) - { - - if (indexEpisode % 2 == 0) // if episodeIndex is even >> noInjection Phase - { - values[Indices::contiNEqIdx] = 0.0; - values[Indices::contiWEqIdx] = 0.0; - values[Indices::energyEqIdx] = 0.0; - } - - else - {// if episodeIndex is odd >> Injection Phase - values[Indices::contiNEqIdx] = -0.0; - values[Indices::contiWEqIdx] = -0.193*2; // (55.5 mol*47.6)/3600 mol/s m = 0.733833333333 - values[Indices::energyEqIdx] = -9132*2; // J/sec m 34788.3611111 - } - } - - else if (globalPos[1] > 2.5 + eps_ && globalPos[1] < 3.5 - eps_) // production well - { - const Scalar elemPressW = elemVolVars[scvIdx].pressure(wPhaseIdx); //Pressures - // const Scalar elemPressG = elemVolVars[scvIdx].pressure(gPhaseIdx); - const Scalar elemPressN = elemVolVars[scvIdx].pressure(nPhaseIdx); - - const Scalar densityW = elemVolVars[scvIdx].fluidState().density(wPhaseIdx); //Densities - // const Scalar densityG = elemVolVars[scvIdx].fluidState().density(gPhaseIdx); - const Scalar densityN = elemVolVars[scvIdx].fluidState().density(nPhaseIdx); - - const Scalar elemMobW = elemVolVars[scvIdx].mobility(wPhaseIdx); //Mobilities - // const Scalar elemMobG = elemVolVars[scvIdx].mobility(gPhaseIdx); - const Scalar elemMobN = elemVolVars[scvIdx].mobility(nPhaseIdx); - - const Scalar enthW = elemVolVars[scvIdx].enthalpy(wPhaseIdx); //Mobilities - // const Scalar enthG = elemVolVars[scvIdx].enthalpy(gPhaseIdx); - const Scalar enthN = elemVolVars[scvIdx].enthalpy(nPhaseIdx); - - const Scalar wellRadius = 0.50 * 0.3048; // 0.50 ft as specified by SPE9 - // const Scalar wellArea = M_PI*std::pow(wellRadius,2); // [m^2] - - const Scalar gridHeight_ = 0.5; - const Scalar effectiveRadius_ = 0.208 * gridHeight_; //Peaceman's Well Model - - const Scalar qW = (((2*3.1415*0.5*4e-14)/(std::log(effectiveRadius_/wellRadius))) * //divided by molarMass() of water to convert from kg/m s to mol/m s - densityW * elemMobW * ( elemPressW-pOut_))/0.01801528; - - const Scalar qN = (((2*3.1415*0.5*4e-14)/(std::log(effectiveRadius_/wellRadius))) * - densityN * elemMobN * (elemPressN-pOut_))/0.35; //divided by molarMass() of HeavyOil to convert from kg/m s to mol/m s - - //without cooling: - // const Scalar qE = qW*0.018*enthW + qN*enthN*0.350; - - //with cooling: see Diplomarbeit Stefan Roll, Sept. 2015 - Scalar wT = elemVolVars[scvIdx].temperature(); // well temperature - Scalar qE; - if ( wT > 495. ) - { - qE = qW*0.018*enthW + qN*enthN*0.350 + (wT-495.)*5000.; // ~3x injected enthalpy - std::cout<< "Cooling now! Extracted enthalpy: " << qE << std::endl; - } - - else - { - qE = qW*0.018*enthW + qN*enthN*0.350; - } - - values[Indices::contiWEqIdx] = qW; - values[Indices::contiNEqIdx] = qN; - values[Indices::energyEqIdx] = qE; - - massProducedOil_ = qN; - massProducedWater_ = qW; - } - } - - // \} - - /*! - * \name Volume terms - */ - // \{ - - /*! - * \brief Evaluate the initial value for a control volume. - * - * \param values The initial values for the primary variables - * \param globalPos The position for which the initial condition should be evaluated - * - * For this method, the \a values parameter stores primary - * variables. - */ - void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const - { - initial_(values, globalPos); - } - - /*! - * \brief Return the initial phase state inside a control volume. - * - * \param globalPos The global position - */ - int initialPhasePresenceAtPos(const GlobalPosition &globalPos) const - { - return wnPhaseOnly; - } - -private: - // internal method for the initial condition (reused for the dirichlet conditions!) - void initial_(PrimaryVariables &values, const GlobalPosition &globalPos) const - { - Scalar densityW = 1000.0; - values[pressureIdx] = 101300. + (maxDepth_ - globalPos[1])*densityW*9.81; - - values[switch1Idx] = 295.13; // temperature - values[switch2Idx] = 0.3; //NAPL saturation - } - - Scalar maxDepth_; - static constexpr Scalar eps_ = 1e-6 - Scalar pIn_; - Scalar pOut_; - //int vectorSize; - Scalar episodeLength_; - std::string name_; - Scalar totalMassProducedOil_; - Scalar totalMassProducedWater_; - - // TODO very evil hack - mutable Scalar massProducedOil_; - mutable Scalar massProducedWater_; - std::ofstream massBalance; - -}; - -} //end namespace Dumux - -#endif // DUMUX_SAGDCYCLICHYSTPROBLEM_HH diff --git a/lecture/mm/heavyoil/sagdcyclichyst/sagd_cyclic_hyst.cc b/lecture/mm/heavyoil/sagdcyclichyst/sagd_cyclic_hyst.cc index 28a4415ff5a1ad12e6d9ebe8742adb0a043700ad..2a2e80bcf9d285809a2c12e674913f77bbe45fdd 100644 --- a/lecture/mm/heavyoil/sagdcyclichyst/sagd_cyclic_hyst.cc +++ b/lecture/mm/heavyoil/sagdcyclichyst/sagd_cyclic_hyst.cc @@ -123,7 +123,7 @@ int main(int argc, char** argv) try // the grid variables using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables); auto gridVariables = std::make_shared(problem, fvGridGeometry); - gridVariables->init(x, xOld); + gridVariables->init(x); // get some time loop parameters using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); @@ -132,11 +132,11 @@ int main(int argc, char** argv) try auto dt = getParam("TimeLoop.DtInitial"); // intialize the vtk output module - using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields); + using IOFields = GetPropType; VtkOutputModule vtkWriter(*gridVariables, x, problem->name()); using VelocityOutput = typename GET_PROP_TYPE(TypeTag, VelocityOutput); vtkWriter.addVelocityOutput(std::make_shared(*gridVariables)); - VtkOutputFields::init(vtkWriter); //! Add model specific output fields + IOFields::initOutputModule(vtkWriter); //!< Add model specific output fields vtkWriter.write(0.0); // instantiate time loop diff --git a/lecture/mm/heavyoil/sagdcyclichyst/sagd_cyclic_hyst.cc.SAVE.cc b/lecture/mm/heavyoil/sagdcyclichyst/sagd_cyclic_hyst.cc.SAVE.cc deleted file mode 100644 index 6ec113acca79920c200ef615fa73b10001fbd1f9..0000000000000000000000000000000000000000 --- a/lecture/mm/heavyoil/sagdcyclichyst/sagd_cyclic_hyst.cc.SAVE.cc +++ /dev/null @@ -1,231 +0,0 @@ -// -*- 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 . * - *****************************************************************************/ -/*! - * \file - * - * \brief test for the 3p3cni box model - */ -#include -#include -#include - -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include - -#include -#include - -#include -#include - -#include - -#include -#include - -#include "problem.hh" - - -/*! - * \brief Provides an interface for customizing error messages associated with - * reading in parameters. - * - * \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 options for this program is:\n" - "\t-TimeManager.TEnd End of the simulation [s] \n" - "\t-TimeManager.DtInitial Initial timestep size [s] \n" - "\t-Grid.File Name of the file containing the grid \n" - "\t definition in DGF format\n"; - - std::cout << errorMessageOut - << "\n"; - } -} - -//////////////////////// -// the main function -//////////////////////// -int main(int argc, char** argv) try -{ - using namespace Dumux; - - // define the type tag for this problem - using TypeTag = TTAG(SagdCyclicHystBoxProblemTypeTag); - - // 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, usage); - - // try to create a grid (from the given grid file or the input file) - GridManager 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(leafGridView); - fvGridGeometry->update(); - - // the problem (initial and boundary conditions) - using Problem = typename GET_PROP_TYPE(TypeTag, Problem); - auto problem = std::make_shared(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(problem, fvGridGeometry); - gridVariables->init(x, xOld); - - // get some time loop parameters - using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); - const auto tEnd = getParam("TimeLoop.TEnd"); - const auto maxDt = getParam("TimeLoop.MaxTimeStepSize"); - auto dt = getParam("TimeLoop.DtInitial"); - - // check if we are about to restart a previously interrupted simulation - Scalar restartTime = 0; - if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart")) // Fehlt ebenso - restartTime = getParam("TimeLoop.Restart"); - - // intialize the vtk output module - using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields); - VtkOutputModule vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name()); - VtkOutputFields::init(vtkWriter); //!< Add model specific output fields - vtkWriter.write(0.0); - - // instantiate time loop - auto timeLoop = std::make_shared>(restartTime, dt, tEnd); - timeLoop->setMaxTimeStepSize(maxDt); - - // the assembler with time loop for instationary problem - using Assembler = FVAssembler; - auto assembler = std::make_shared(problem, fvGridGeometry, gridVariables, timeLoop); - - // the linear solver - using LinearSolver = AMGBackend; - auto linearSolver = std::make_shared(leafGridView, fvGridGeometry->dofMapper()); - - // the non-linear solver - using NewtonSolver = Dumux::NewtonSolver; - NewtonSolver nonLinearSolver(assembler, linearSolver); - - // time loop - timeLoop->start(); do - { - // set previous solution for storage evaluations - assembler->setPreviousSolution(xOld); - - // 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(); - - // write vtk output - 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()); - - //////////////////////////////////////////////////////////// - // finalize, print dumux message to say goodbye - //////////////////////////////////////////////////////////// - - // print dumux end message - if (mpiHelper.rank() == 0) - { - Parameters::print(); - DumuxMessage::print(/*firstCall=*/false); - } - - - -} - -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; -} diff --git a/lecture/mm/heavyoil/sagdcyclichyst/spatialparams.hh.SAVE.hh b/lecture/mm/heavyoil/sagdcyclichyst/spatialparams.hh.SAVE.hh deleted file mode 100644 index f3d99d839dda73612f53f9e48056ac48a1cfb94e..0000000000000000000000000000000000000000 --- a/lecture/mm/heavyoil/sagdcyclichyst/spatialparams.hh.SAVE.hh +++ /dev/null @@ -1,350 +0,0 @@ -// -*- 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 . * - *****************************************************************************/ -/*! - * \file - * - * \brief Definition of the spatial parameters for the SAGD problem. - */ -#ifndef DUMUX_SAGDCYCLICHYST_SPATIAL_PARAMS_HH -#define DUMUX_SAGDCYCLICHYST_SPATIAL_PARAMS_HH - -#include - -#include -#include -#include "../3p/parkervangenuchtenzerohysteresis.hh" -#include "../3p/parkervangenuchtenzerohysteresisparams.hh" - -namespace Dumux { - -//forward declaration -template -class SagdCyclicHystSpatialParams; - -namespace Properties { - -// The spatial parameters TypeTag -NEW_TYPE_TAG(SagdCyclicHystSpatialParamsTypeTag); - -// Set the spatial parameters -SET_TYPE_PROP(SagdCyclicHystSpatialParamsTypeTag, SpatialParams, Dumux::SagdCyclicHystSpatialParams); - -} // end namespace Properties - -/*! - * \ingroup ThreePThreeCNIModel - * - * \brief Definition of the spatial parameters for the SAGD problem - */ -template -class SagdCyclicHystSpatialParams : public FVSpatialParams> -{ - using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); - using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); - using ParentType = FVSpatialParams >; - using Problem = typename GET_PROP_TYPE(TypeTag, Problem); - using Grid = typename GET_PROP_TYPE(TypeTag, Grid); - using GridView = typename GET_PROP_TYPE(TypeTag, GridView); - using CoordScalar = typename Grid::ctype; - static constexpr int dim = GridView::dimension; - static constexpr int dimWorld=GridView::dimensionworld; - using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices; - using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits); - static constexpr int numPhases = ModelTraits::numPhases; - static constexpr int numComponents = ModelTraits::NumComponents; - static constexpr int wPhaseIdx = Indices::Phase0Idx; - static constexpr int nPhaseIdx = Indices::Phase1Idx; - static constexpr int gPhaseIdx = Indices::Phase2Idx; - using Element = typename GridView::template Codim<0>::Entity; - using GlobalPosition = typename Element::Geometry::GlobalCoordiante; -// TODO: What is this used for? Is this correct? - using DimVetor = Dune::FieldVector; - using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector); - using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables); - using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView; - using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView; - using ElementIterator = typename GridView::template Codim<0>::Iterator; - using VertexIterator = typename GridView::template Codim::Iterator; - using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables); - static constexpr auto isBox = (FVGridGeometry::discMethod == DiscretizationMethod::box); - static constexpr auto dofCodim = isBox ? dim : 0; - -public: - using EffectiveLaw = ParkerVanGenZeroHyst3P; - using MaterialLaw = EffToAbsLaw; - using MaterialLawParams = typename MaterialLaw::Params; - using MaterialLawHystParamsVector = std::vector; - - /*! - * \brief The constructor - * - * \param gridView The grid view - */ - SagdCyclicHystSpatialParams(std::shared_ptr fvGridGeometry) : ParentType(fvGridGeometry) - { - layerBottom_ = 35.0; - // intrinsic permeabilities - fineK_ = 1e-16; - coarseK_ = 4e-14; - // porosities - finePorosity_ = 0.10; - coarsePorosity_ = 0.1; - // heat conductivity of granite - lambdaSolid_ = 2.8; - // specific heat capacities - fineHeatCap_ = 850.; - coarseHeatCap_ = 850; - // residual saturations - fineMaterialParams_.setSwr(0.1); - fineMaterialParams_.setSwrx(0.12); //Total liquid Residual Saturation - fineMaterialParams_.setSnr(0.09); //Residual of NAPL if there is no water - fineMaterialParams_.setSgr(0.01); - coarseMaterialParams_.setSwr(0.1); - coarseMaterialParams_.setSwrx(0.12); - coarseMaterialParams_.setSnr(0.09); - coarseMaterialParams_.setSgr(0.01); - - // parameters for the 3phase van Genuchten law - fineMaterialParams_.setVgn(4.0); - coarseMaterialParams_.setVgn(4.0); - fineMaterialParams_.setVgAlpha(0.0005); - coarseMaterialParams_.setVgAlpha(0.0015); - coarseMaterialParams_.setKrRegardsSnr(false); - fineMaterialParams_.setKrRegardsSnr(false); - } - - /*! - * \brief Apply the intrinsic permeability tensor to a pressure - * potential gradient. - * - * \param element The current finite element - * \param fvElemGeom The current finite volume geometry of the element - * \param scvIdx The index of the sub-control volume - */ - Scalar intrinsicPermeability(const Element &element, - const FVElementGeometry &fvElemGeom, - int scvIdx) const - { - const GlobalPosition &pos = fvElemGeom.subContVol[scvIdx].global; - if (isFineMaterial_(pos)) - return fineK_; - - return coarseK_; - } - - /*! - * \brief Define the porosity \f$[-]\f$ of the spatial parameters - * - * \param element The finite element - * \param fvGeometry The finite volume geometry - * \param scvIdx The local index of the sub-control volume where - * the porosity needs to be defined - */ - Scalar porosity(const Element &element, - const FVElementGeometry &fvGeometry, - const int scvIdx) const - { - const GlobalPosition &pos = fvGeometry.subContVol[scvIdx].global; - if (isFineMaterial_(pos)) - return finePorosity_; - - else - return coarsePorosity_; - } - - /*! - * \brief return the parameter object for the Brooks-Corey material law which depends on the position - * - * \param element The current finite element - * \param fvGeometry The current finite volume geometry of the element - * \param scvIdx The index of the sub-control volume - */ - const MaterialLawParams& materialLawParams(const Element &element, - const FVElementGeometry &fvGeometry, - const int scvIdx) const - { - const GlobalPosition &pos = fvGeometry.subContVol[scvIdx].global; - if (isFineMaterial_(pos)) - return fineMaterialParams_; - - else - return coarseMaterialParams_; - } - - /*! - * \brief Returns the heat capacity \f$[J / (kg K)]\f$ of the rock matrix. - * - * 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 - */ - Scalar solidHeatCapacity(const Element &element, - const FVElementGeometry &fvGeometry, - const int scvIdx) const - { - const GlobalPosition &pos = fvGeometry.subContVol[scvIdx].global; - if (isFineMaterial_(pos)) - return fineHeatCap_; - - else - return coarseHeatCap_; - } - - /*! - * \brief Returns the mass density \f$[kg / m^3]\f$ of the rock matrix. - * - * 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 - */ - Scalar solidDensity(const Element &element, - const FVElementGeometry &fvGeometry, - const int scvIdx) const - { - return 2650; // density of sand [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. - * - * \param element The finite element - * \param fvGeometry The finite volume geometry of the element - * \param scvIdx The local index of the sub-control volume - */ - Scalar solidThermalConductivity(const Element &element, - const FVElementGeometry &fvGeometry, - const int scvIdx) const - { - return lambdaSolid_; - } - - struct MaxSaturations - { - Scalar MaxSatW; - Scalar MaxSatN; - Scalar MaxSatG; - }; - - struct trappedSaturations - { - Scalar trappedSatN; - }; - - void getMaxSaturation(Problem &problem) - { - // get the number of degrees of freedom - FVElementGeometry fvGeometry; - VolumeVariables volVars; - maxSats_.resize(problem.model().numDofs()); - - for (const auto& element : elements(problem.gridView())) - { - fvGeometry.update(problem.gridView(), element); - - for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx) - { - const unsigned int dofIdxGlobal = problem.model().dofMapper().subIndex(element, scvIdx, dofCodim); - volVars.update(problem.model().curSol()[dofIdxGlobal], - problem, - element, - fvGeometry, - scvIdx, - false); - - if (volVars.saturation(wPhaseIdx)>maxSats_[dofIdxGlobal].MaxSatW) - maxSats_[dofIdxGlobal].MaxSatW = volVars.saturation(wPhaseIdx); - - if (volVars.saturation(wPhaseIdx)>maxSats_[dofIdxGlobal].MaxSatW) - maxSats_[dofIdxGlobal].MaxSatW = volVars.saturation(wPhaseIdx); - - if (volVars.saturation(nPhaseIdx)>maxSats_[dofIdxGlobal].MaxSatN) - maxSats_[dofIdxGlobal].MaxSatN = volVars.saturation(nPhaseIdx); - - if (volVars.saturation(gPhaseIdx)>maxSats_[dofIdxGlobal].MaxSatG) - maxSats_[dofIdxGlobal].MaxSatG = volVars.saturation(gPhaseIdx); - } - } - } - - void trappedSat(Problem &problem) - { - // get the number of degrees of freedom - FVElementGeometry fvGeometry; - VolumeVariables volVars; - trapSats_.resize(problem.model().numDofs()); - - for (const auto& element : elements(problem.gridView())) - { - fvGeometry.update(problem.gridView(), element); - for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx) - { - const unsigned int dofIdxGlobal = problem.model().dofMapper().subIndex(element, scvIdx, dofCodim); - volVars.update(problem.model().curSol()[dofIdxGlobal], - problem, - element, - fvGeometry, - scvIdx, - false); - trapSats_[dofIdxGlobal].trappedSatN = - (maxSats_[dofIdxGlobal].MaxSatN)/(1+9*maxSats_[dofIdxGlobal].MaxSatN); - - coarseMaterialParams_.setTrappedSatN(trapSats_[dofIdxGlobal].trappedSatN); - fineMaterialParams_.setTrappedSatN(trapSats_[dofIdxGlobal].trappedSatN); - } - } - } - -private: - bool isFineMaterial_(const GlobalPosition &pos) const - { - return pos[dim-1] > layerBottom_ - eps_; - } - - std::vector trapSats_; - - Scalar layerBottom_; - Scalar lambdaSolid_; - - Scalar fineK_; - Scalar coarseK_; - - Scalar finePorosity_; - Scalar coarsePorosity_; - - Scalar fineHeatCap_; - Scalar coarseHeatCap_; - - MaterialLawParams fineMaterialParams_; - MaterialLawParams coarseMaterialParams_; - std::vector maxSats_; - - MaterialLawHystParamsVector materialLawHystParams_; - - static constexpr Scalar eps_ = 1.5e-7; -}; - -} // end namespace Dumux - -#endif // DUMUX_SAGDCYCLIC_SPATIAL_PARAMS_HH