diff --git a/exercises/exercise-basic/2p2cmain.cc b/exercises/exercise-basic/2p2cmain.cc deleted file mode 100644 index b9a07f08e16dce4b177694c369dea916aac494cc..0000000000000000000000000000000000000000 --- a/exercises/exercise-basic/2p2cmain.cc +++ /dev/null @@ -1,154 +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 3 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 2p2c porousmediumflow problem in exercise-basic - */ -#include <config.h> - -#include <iostream> - -#include <dumux/common/initialize.hh> -#include <dumux/common/properties.hh> -#include <dumux/common/parameters.hh> - -#include <dumux/linear/istlsolvers.hh> -#include <dumux/linear/linearsolvertraits.hh> -#include <dumux/linear/linearalgebratraits.hh> -#include <dumux/nonlinear/newtonsolver.hh> - -#include <dumux/assembly/fvassembler.hh> -#include <dumux/assembly/diffmethod.hh> - -#include <dumux/io/vtkoutputmodule.hh> -#include <dumux/io/grid/gridmanager_yasp.hh> - -// The properties file, where the compile time options are defined -#include "properties2p2c.hh" - -//////////////////////// -// the main function -//////////////////////// -int main(int argc, char** argv) -{ - using namespace Dumux; - - // define the type tag for this problem - using TypeTag = Properties::TTag::Injection2p2cCC; - - // initialize MPI+X, finalize is done automatically on exit - Dumux::initialize(argc, argv); - - // 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<GetPropType<TypeTag, Properties::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 GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; - auto gridGeometry = std::make_shared<GridGeometry>(leafGridView); - - // the problem (initial and boundary conditions) - using Problem = GetPropType<TypeTag, Properties::Problem>; - auto problem = std::make_shared<Problem>(gridGeometry); - - // the solution vector - using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>; - SolutionVector x; - problem->applyInitialSolution(x); - auto xOld = x; - - // the grid variables - using GridVariables = GetPropType<TypeTag, Properties::GridVariables>; - auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry); - gridVariables->init(x); - - // initialize the vtk output module - using IOFields = GetPropType<TypeTag, Properties::IOFields>; - VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name()); - using VelocityOutput = GetPropType<TypeTag, Properties::VelocityOutput>; - vtkWriter.addVelocityOutput(std::make_shared<VelocityOutput>(*gridVariables)); - IOFields::initOutputModule(vtkWriter); //!< Add model specific output fields - vtkWriter.write(0.0); - - // instantiate time loop - using Scalar = GetPropType<TypeTag, Properties::Scalar>; - const auto tEnd = getParam<Scalar>("TimeLoop.TEnd"); - const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize"); - const auto dt = getParam<Scalar>("TimeLoop.DtInitial"); - 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, gridGeometry, gridVariables, timeLoop, xOld); - - // the linear solver - using LinearSolver = AMGBiCGSTABIstlSolver<LinearSolverTraits<GridGeometry>, LinearAlgebraTraitsFromAssembler<Assembler>>; - auto linearSolver = std::make_shared<LinearSolver>(gridGeometry->gridView(), gridGeometry->dofMapper()); - - // the non-linear solver -// using PrimaryVariableSwitch = GetPropType<TypeTag, Properties::PrimaryVariableSwitch>; - using NewtonSolver = NewtonSolver<Assembler, LinearSolver>; - NewtonSolver nonLinearSolver(assembler, linearSolver); - - // time loop - timeLoop->start(); - while (!timeLoop->finished()) - { - //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()); - - // print parameter report - if (leafGridView.comm().rank() == 0) - Parameters::print(); - - return 0; -} // end main diff --git a/exercises/exercise-basic/CMakeLists.txt b/exercises/exercise-basic/CMakeLists.txt index e70eefef42b098146db4cbe4455c8f01c163e278..95bff476c7291c7e8d6d5dd0975e593dde0cb403 100644 --- a/exercises/exercise-basic/CMakeLists.txt +++ b/exercises/exercise-basic/CMakeLists.txt @@ -2,12 +2,8 @@ dumux_add_test(NAME exercise_basic_2p SOURCES 2pmain.cc) -# the compositional two-phase two-component simulation program -dumux_add_test(NAME exercise_basic_2p2c - SOURCES 2p2cmain.cc) - # here, add the two-phase non-isothermal simulation program # add a symlink for each input file -add_input_file_links() +add_input_file_links() \ No newline at end of file diff --git a/exercises/exercise-basic/README.md b/exercises/exercise-basic/README.md index 5dd08af047defb938ca1016978eea869d699a7d1..839644a057ca142394ee2382976a854b7526d773 100644 --- a/exercises/exercise-basic/README.md +++ b/exercises/exercise-basic/README.md @@ -12,7 +12,7 @@ The aquifer is situated 2700 m below sea level and the domain size is 60 m x 40 * Navigate to the directory `dumux-course/exercises/exercise-basic` -This exercise deals with two problems: a two-phase immiscible problem (__2p__) and a two-phase compositional problem (__2p2c__). They both set up the same scenario with the difference that the 2p2c assumes a miscible fluid state for the two fluids (water and gaseous N$_2$) and the 2p model assumes an immiscible fluid state. +This exercise deals with two problems: a two-phase immiscible problem (__2p__) and a two-phase non-isothermal problem (__2pni__). They both set up the same scenario with the difference that the 2pni model introduces an extra energy equation. <br><br> ### Task 1: Getting familiar with the code @@ -20,11 +20,10 @@ This exercise deals with two problems: a two-phase immiscible problem (__2p__) a Locate all the files you will need for this exercise * The __main file__ for the __2p__ problem : `2pmain.cc` -* The __main file__ for the __2p2c__ problem : `2p2cmain.cc` * The __problem file__ for the __2p__ problem: `injection2pproblem.hh` -* The __problem file__ for the __2p2c__ problem: `injection2p2cproblem.hh` +* The __problem file__ for the __2pni__ problem: `injection2pniproblem.hh` * The __properties file__ for the __2p__ problem: `properties2p.hh` -* The __properties file__ for the __2p2c__ problem: `properties2p2c.hh` +* The __properties file__ for the __2pni__ problem: `properties2pni.hh` * The shared __spatial parameters file__: `injection2pspatialparams.hh` * The shared __input file__: `params.input` @@ -38,23 +37,22 @@ Locate all the files you will need for this exercise cd ../../build-cmake/exercises/exercise-basic ``` -* Compile both executables `exercise_basic_2p` and `exercise_basic_2p2c` +* Compile the executable `exercise_basic_2p` ```bash -make exercise_basic_2p exercise_basic_2p2c +make exercise_basic_2p ``` -* Execute the two problems and inspect the result +* Execute the problem and inspect the result ```bash ./exercise_basic_2p params.input -./exercise_basic_2p2c params.input ``` * you can look at the results with paraview ```bash -paraview injection-2p2c.pvd +paraview injection-2p.pvd ``` <br><br> diff --git a/exercises/exercise-basic/injection2p2cproblem.hh b/exercises/exercise-basic/injection2p2cproblem.hh deleted file mode 100644 index 556be438c6357f77723804163bdc44cc3fbed409..0000000000000000000000000000000000000000 --- a/exercises/exercise-basic/injection2p2cproblem.hh +++ /dev/null @@ -1,229 +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 3 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-basic - */ - -#ifndef DUMUX_EX_BASIC_PROBLEM_2P2C_HH -#define DUMUX_EX_BASIC_PROBLEM_2P2C_HH - -#include <dumux/common/properties.hh> -#include <dumux/common/boundarytypes.hh> -#include <dumux/common/numeqvector.hh> -#include <dumux/porousmediumflow/problem.hh> -#include <dumux/material/binarycoefficients/h2o_n2.hh> - -namespace Dumux { - -/*! - * \ingroup TwoPTwoCModel - * \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 Injection2p2cProblem : public PorousMediumFlowProblem<TypeTag> -{ - using ParentType = PorousMediumFlowProblem<TypeTag>; - using Scalar = GetPropType<TypeTag, Properties::Scalar>; - using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices; - using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>; - using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>; - using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; - using FVElementGeometry = typename GridGeometry::LocalView; - using GridView = typename GridGeometry::GridView; - using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>; - using NumEqVector = Dumux::NumEqVector<PrimaryVariables>; - - static constexpr int dimWorld = GridView::dimensionworld; - using Element = typename GridView::template Codim<0>::Entity; - using GlobalPosition = typename Element::Geometry::GlobalCoordinate; - -public: - Injection2p2cProblem(std::shared_ptr<const GridGeometry> gridGeometry) - : ParentType(gridGeometry) - { - // 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"); - } - - /*! - * \brief Returns the problem name - * - * This is used as a prefix for files generated by the simulation. - */ - std::string name() const - { return name_+"-2p2c"; } - - /*! - * \brief Specifies which kind of boundary condition should be - * used for which equation on a given boundary segment. - * - * \param globalPos The position for which the bc type should be evaluated - */ - BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const - { - BoundaryTypes bcTypes; - if (globalPos[0] < eps_) - bcTypes.setAllDirichlet(); - else - bcTypes.setAllNeumann(); - - return bcTypes; - } - - /*! - * \brief Evaluates the boundary conditions for a Dirichlet - * boundary segment - * - * \param globalPos The global position - */ - PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const - { - return initialAtPos(globalPos); - } - - /*! - * \brief Evaluate the boundary conditions for a Neumann - * boundary segment. - * - * \param globalPos The position of the integration point of the boundary segment. - */ - NumEqVector neumannAtPos(const GlobalPosition &globalPos) const - { - // initialize values to zero, i.e. no-flow Neumann boundary conditions - NumEqVector values(0.0); - - // if we are inside the injection zone set inflow Neumann boundary conditions - if (injectionActive() && onInjectionBoundary(globalPos)) - { - // TODO: dumux-course-task - //instead of setting -1e-4 here directly use totalAreaSpecificInflow_ in the computation - - // inject nitrogen. negative values mean injection - // convert from units kg/(s*m^2) to mole/(s*m^2) - values[Indices::conti0EqIdx + FluidSystem::N2Idx] = -1e-4/FluidSystem::molarMass(FluidSystem::N2Idx); - values[Indices::conti0EqIdx + FluidSystem::H2OIdx] = 0.0; - } - - return values; - } - - /*! - * \brief Evaluate the source term for all phases within a given - * sub-control-volume. - * - * \param globalPos The position for which the source term should be evaluated - */ - NumEqVector sourceAtPos(const GlobalPosition &globalPos) const - { - return NumEqVector(0.0); - } - - /*! - * \brief Evaluate the initial value for a control volume. - * - * \param globalPos The position for which the initial condition should be evaluated - */ - PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const - { - PrimaryVariables values(0.0); - values.setState(Indices::firstPhaseOnly); - // get the water density at atmospheric conditions - const Scalar densityW = FluidSystem::H2O::liquidDensity(this->spatialParams().temperatureAtPos(globalPos), 1.0e5); - - // assume an initially hydrostatic liquid pressure profile - // note: we subtract rho_w*g*h because g is defined negative - const Scalar pw = 1.0e5 - densityW*this->spatialParams().gravity(globalPos)[dimWorld-1]*(aquiferDepth_ - globalPos[dimWorld-1]); - - // initially we have some nitrogen dissolved - // saturation mole fraction would be - // moleFracLiquidN2 = (pw + pc + p_vap^sat)/henry; - const Scalar moleFracLiquidN2 = pw*0.95/BinaryCoeff::H2O_N2::henry(this->spatialParams().temperatureAtPos(globalPos)); - - // note that because we start with a single phase system the primary variables - // are pl and x^w_N2. This will switch as soon after we start injecting to a two - // phase system so the primary variables will be pl and Sn (nonwetting saturation). - values[Indices::pressureIdx] = pw; - values[Indices::switchIdx] = moleFracLiquidN2; - - return values; - } - - //! set the time for the time dependent boundary conditions (called from main) - void setTime(Scalar time) - { time_ = time; } - - //! Return true if the injection is currently active - bool injectionActive() const - { return time_ < injectionDuration_; } - - //! Return true if the given position is in the injection boundary region - bool onInjectionBoundary(const GlobalPosition& globalPos) const - { - return globalPos[1] < 15. + eps_ - && globalPos[1] > 7. - eps_ - && globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; - } - -private: - static constexpr Scalar eps_ = 1e-6; - std::string name_; //! Problem name - Scalar aquiferDepth_; //! Depth of the aquifer in m - Scalar injectionDuration_; //! Duration of the injection in seconds - Scalar time_; - //TODO: dumux-course-task - //define the Scalar totalAreaSpecificInflow_ here - -}; - -} //end namespace Dumux - -#endif diff --git a/exercises/exercise-basic/properties2p2c.hh b/exercises/exercise-basic/properties2p2c.hh deleted file mode 100644 index dab9384cc064f50a262e6ac13032ab21bc36b7c1..0000000000000000000000000000000000000000 --- a/exercises/exercise-basic/properties2p2c.hh +++ /dev/null @@ -1,78 +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 3 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 two-component porousmediumflow properties file for exercise-basic - */ - -#ifndef DUMUX_EX_BASIC_PROPERTIES_2P2C_HH -#define DUMUX_EX_BASIC_PROPERTIES_2P2C_HH - -#include <dune/grid/yaspgrid.hh> - -#include <dumux/discretization/cctpfa.hh> -#include <dumux/porousmediumflow/2p2c/model.hh> -#include <dumux/material/fluidsystems/h2on2.hh> - -#include "injection2p2cproblem.hh" -#include "injection2pspatialparams.hh" - -namespace Dumux::Properties { - -// Create new type tags -namespace TTag { -struct Injection2p2c { using InheritsFrom = std::tuple<TwoPTwoC>; }; -struct Injection2p2cCC { using InheritsFrom = std::tuple<Injection2p2c, CCTpfaModel>; }; -} // end namespace TTag - -// Set the grid type -template<class TypeTag> -struct Grid<TypeTag, TTag::Injection2p2c> { using type = Dune::YaspGrid<2>; }; - -// Set the problem property -template<class TypeTag> -struct Problem<TypeTag, TTag::Injection2p2c> { using type = Injection2p2cProblem<TypeTag>; }; - -// Set the spatial parameters -template<class TypeTag> -struct SpatialParams<TypeTag, TTag::Injection2p2c> -{ -private: - using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; - using Scalar = GetPropType<TypeTag, Properties::Scalar>; -public: - using type = InjectionSpatialParams<GridGeometry, Scalar>; -}; - -// Set fluid configuration -template<class TypeTag> -struct FluidSystem<TypeTag, TTag::Injection2p2c> -{ - using type = FluidSystems::H2ON2< GetPropType<TypeTag, Properties::Scalar>, - FluidSystems::H2ON2DefaultPolicy</*fastButSimplifiedRelations=*/ true> >; -}; - -// Define whether mole (true) or mass (false) fractions are used -template<class TypeTag> -struct UseMoles<TypeTag, TTag::Injection2p2c> { static constexpr bool value = true; }; - -} // end namespace Dumux::Properties - -#endif