diff --git a/.patches/exercise-basic/exercise-basic.patch b/.patches/exercise-basic/exercise-basic.patch index 7e0ee9a827d59e0d4dc4282ac65782bfed3c1e37..43e320fbdf6e6245526b91e26def9cf46879a7aa 100644 --- a/.patches/exercise-basic/exercise-basic.patch +++ b/.patches/exercise-basic/exercise-basic.patch @@ -1,163 +1,5 @@ -diff -ruN exercises/exercise-basic/2p2cmain.cc exercises/solution/exercise-basic/2p2cmain.cc ---- exercises/exercise-basic/2p2cmain.cc 2024-07-08 09:06:59.581159992 +0200 -+++ exercises/solution/exercise-basic/2p2cmain.cc 1970-01-01 01:00:00.000000000 +0100 -@@ -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 -ruN exercises/exercise-basic/2pmain.cc exercises/solution/exercise-basic/2pmain.cc ---- exercises/exercise-basic/2pmain.cc 2024-07-08 09:06:59.581159992 +0200 +--- exercises/exercise-basic/2pmain.cc 2023-03-31 13:48:56.158686713 +0200 +++ exercises/solution/exercise-basic/2pmain.cc 1970-01-01 01:00:00.000000000 +0100 @@ -1,157 +0,0 @@ -// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- @@ -319,7 +161,7 @@ diff -ruN exercises/exercise-basic/2pmain.cc exercises/solution/exercise-basic/2 -} // end main diff -ruN exercises/exercise-basic/2pnimain.cc exercises/solution/exercise-basic/2pnimain.cc --- exercises/exercise-basic/2pnimain.cc 1970-01-01 01:00:00.000000000 +0100 -+++ exercises/solution/exercise-basic/2pnimain.cc 2024-07-08 09:06:59.613159864 +0200 ++++ exercises/solution/exercise-basic/2pnimain.cc 2023-03-31 13:48:56.182686617 +0200 @@ -0,0 +1,150 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: @@ -472,17 +314,13 @@ diff -ruN exercises/exercise-basic/2pnimain.cc exercises/solution/exercise-basic + return 0; +} // end main diff -ruN exercises/exercise-basic/CMakeLists.txt exercises/solution/exercise-basic/CMakeLists.txt ---- exercises/exercise-basic/CMakeLists.txt 2024-07-08 09:06:59.581159992 +0200 -+++ exercises/solution/exercise-basic/CMakeLists.txt 2024-07-08 09:06:59.613159864 +0200 -@@ -1,13 +1,6 @@ +--- exercises/exercise-basic/CMakeLists.txt 2024-07-16 17:38:07.971854664 +0200 ++++ exercises/solution/exercise-basic/CMakeLists.txt 2023-03-31 13:48:56.182686617 +0200 +@@ -1,9 +1,6 @@ -# the immiscible two-phase simulation program -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 - +# the two-phase non-isothermal simulation program @@ -490,243 +328,12 @@ diff -ruN exercises/exercise-basic/CMakeLists.txt exercises/solution/exercise-ba + SOURCES 2pnimain.cc) # add a symlink for each input file - add_input_file_links() -diff -ruN exercises/exercise-basic/injection2p2cproblem.hh exercises/solution/exercise-basic/injection2p2cproblem.hh ---- exercises/exercise-basic/injection2p2cproblem.hh 2024-07-08 09:06:59.581159992 +0200 -+++ exercises/solution/exercise-basic/injection2p2cproblem.hh 1970-01-01 01:00:00.000000000 +0100 -@@ -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 +-add_input_file_links() +\ No newline at end of file ++add_input_file_links() diff -ruN exercises/exercise-basic/injection2pniproblem.hh exercises/solution/exercise-basic/injection2pniproblem.hh ---- exercises/exercise-basic/injection2pniproblem.hh 2024-07-16 16:55:12.734099365 +0200 -+++ exercises/solution/exercise-basic/injection2pniproblem.hh 2024-07-16 16:55:12.738099338 +0200 +--- exercises/exercise-basic/injection2pniproblem.hh 2024-07-16 17:38:07.971854664 +0200 ++++ exercises/solution/exercise-basic/injection2pniproblem.hh 2024-07-16 17:38:07.971854664 +0200 @@ -19,7 +19,7 @@ /*! * \file @@ -774,102 +381,100 @@ diff -ruN exercises/exercise-basic/injection2pniproblem.hh exercises/solution/ex static constexpr int dimWorld = GridView::dimensionworld; using Element = typename GridView::template Codim<0>::Entity; using GlobalPosition = typename Element::Geometry::GlobalCoordinate; -@@ -115,14 +116,6 @@ +@@ -83,13 +84,13 @@ + /*pMax=*/30e6, + /*numP=*/300); + +- // name of the problem and output file ++ // 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, unit: m ++ // Depth of the aquifer, unit: m + aquiferDepth_ = getParam<Scalar>("Problem.AquiferDepth"); +- // the duration of the injection, unit: seconds ++ // The duration of the injection, unit: seconds + injectionDuration_ = getParam<Scalar>("Problem.InjectionDuration"); + } + +@@ -115,13 +116,6 @@ else bcTypes.setAllNeumann(); - /*! -- * TODO:dumux-course-task: -- * dumux-course-task: -- * set Dirichlet conditions for the energy equation on the left boundary -- * and Neumann everywhere else -- * think about: is there anything necessary to do here? +- * TODO:dumux-course-task 4: +- * Set Dirichlet conditions for the energy equation on the left boundary +- * and Neumann everywhere else. +- * Think about: is there anything necessary to do here? - */ - return bcTypes; } -@@ -135,13 +128,6 @@ +@@ -134,12 +128,6 @@ PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const { return initialAtPos(globalPos); - - /*! -- * TODO:dumux-course-task: -- * dumux-course-task: -- * set Dirichlet conditions for the energy equation on the left boundary -- * think about: is there anything necessary to do here? +- * TODO:dumux-course-task 4: +- * Set Dirichlet conditions for the energy equation on the left boundary. +- * Think about: is there anything necessary to do here? - */ } /*! -@@ -156,19 +142,17 @@ +@@ -154,19 +142,20 @@ NumEqVector values(0.0); // if we are inside the injection zone set inflow Neumann boundary conditions - if (injectionActive() && onInjectionBoundary(globalPos)) + if (injectionActive() && onInjectionBoundary(globalPos)) { -+ - // inject nitrogen. negative values mean injection - // units kg/(s*m^2) - values[Indices::conti0EqIdx + FluidSystem::N2Idx] = -1e-4; +- // inject nitrogen. negative values mean injection ++ const Scalar injectionRate = -1e-4; ++ ++ // inject nitrogen. Negative values mean injection + // unit: kg/(s*m^2) +- values[Indices::conti0EqIdx + FluidSystem::N2Idx] = -1e-4; ++ values[Indices::conti0EqIdx + FluidSystem::N2Idx] = injectionRate; values[Indices::conti0EqIdx + FluidSystem::H2OIdx] = 0.0; - /*! -- * TODO:dumux-course-task: +- * TODO:dumux-course-task 4: - * Set Neumann noflow conditions for the energy equation everywhere else except the left boundary. -- * Additionally, consider the energy flux at the injection point which is equal to the product of the respective mass flux and the matching enthalpy. Use the function `gasEnthalpy(...)` from the N2 component to access the necessary enthalpy. +- * Additionally, consider the energy flux at the injection point which is equal to the product of the respective mass flux and the matching enthalpy. Use the function `gasEnthalpy(temperature,pressure)` from the N2 component to access the necessary enthalpy. - * hint: use `Indices::energyEqIdx` to access the entry belonging to the energy flux. - */ + // energy fluxes are always mass specific -+ //units W/(m^2) -+ values[Indices::energyEqIdx] = -1e-4 /*kg/(m^2 s)*/ *N2::gasEnthalpy(injectionTemperature(), 1e5 /*Pa*/)/*J/kg*/; ++ // unit: W/(m^2) ++ const Scalar temperatureAtInjection = initialAtPos(globalPos)[Indices::temperatureIdx];/*K*/ ++ const Scalar pressureAtInjection = initialAtPos(globalPos)[Indices::pressureIdx];/*Pa*/ ++ values[Indices::energyEqIdx] = injectionRate /*kg/(m^2 s)*/ *N2::gasEnthalpy(temperatureAtInjection, pressureAtInjection)/*J/kg*/; } return values; -@@ -204,20 +188,17 @@ +@@ -202,13 +191,10 @@ values[Indices::pressureIdx] = pw; values[Indices::saturationIdx] = 0.0; - /*! -- * TODO:dumux-course-task: -- * set a temperature gradient of 0.03 K per m beginning at 283 K here. -- * Hint: you can use aquiferDepth_ and the globalPos similar to the pressure gradient -- * use globalPos[0] and globalPos[1] to implement the high temperature lens with 380 K +- * TODO:dumux-course-task 4: +- * Set a temperature gradient of 0.03 K per m beginning at 283 K here. +- * Hint: you can use aquiferDepth_ and the globalPos similar to the pressure gradient. +- * Use globalPos[0] and globalPos[1] to implement the high temperature lens with 380 K - * Hint : use Indices::temperatureIdx to address the initial values for temperature - */ + values[Indices::temperatureIdx] = 283.0 + (aquiferDepth_ - globalPos[1])*0.03; + if (globalPos[0] > 20 - eps_ && globalPos[0] < 30 + eps_ && globalPos[1] > 5 - eps_ && globalPos[1] < 35 + eps_) -+ values[Indices::temperatureIdx] = injectionTemperature(); ++ values[Indices::temperatureIdx] = 380.0; + return values; } -- //! set the time for the time dependent boundary conditions (called from main) -+ //! 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_; } -@@ -230,6 +211,12 @@ - && globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; - } - -+ //! Set injection temperature for the domain -+ Scalar injectionTemperature() const -+ { -+ return 380.0; -+ } -+ - private: - static constexpr Scalar eps_ = 1e-6; - std::string name_; //! Problem name diff -ruN exercises/exercise-basic/injection2pproblem.hh exercises/solution/exercise-basic/injection2pproblem.hh ---- exercises/exercise-basic/injection2pproblem.hh 2024-07-08 09:06:59.581159992 +0200 +--- exercises/exercise-basic/injection2pproblem.hh 2024-07-16 17:38:07.971854664 +0200 +++ exercises/solution/exercise-basic/injection2pproblem.hh 1970-01-01 01:00:00.000000000 +0100 @@ -1,223 +0,0 @@ -// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- @@ -949,7 +554,7 @@ diff -ruN exercises/exercise-basic/injection2pproblem.hh exercises/solution/exer - Injection2PProblem(std::shared_ptr<const GridGeometry> gridGeometry) - : ParentType(gridGeometry) - { -- // initialize the tables of the fluid system +- // Initialize the tables of the fluid system - FluidSystem::init(/*tempMin=*/273.15, - /*tempMax=*/423.15, - /*numTemp=*/50, @@ -957,13 +562,13 @@ diff -ruN exercises/exercise-basic/injection2pproblem.hh exercises/solution/exer - /*pMax=*/30e6, - /*numP=*/300); - -- // name of the problem and output file +- // 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 +- // Depth of the aquifer, unit: m - aquiferDepth_ = getParam<Scalar>("Problem.AquiferDepth"); -- // the duration of the injection, units: second +- // The duration of the injection, unit: seconds - injectionDuration_ = getParam<Scalar>("Problem.InjectionDuration"); - } - @@ -1023,8 +628,8 @@ diff -ruN exercises/exercise-basic/injection2pproblem.hh exercises/solution/exer - // than using <= or >= as it is robust with regard to imprecision introduced by rounding errors. - if (injectionActive() && onInjectionBoundary(globalPos)) - { -- // inject nitrogen. negative values mean injection -- // units kg/(s*m^2) +- // Inject nitrogen. Negative values mean injection +- // unit: kg/(s*m^2) - values[Indices::conti0EqIdx + FluidSystem::N2Idx] = -1e-4; - values[Indices::conti0EqIdx + FluidSystem::H2OIdx] = 0.0; - } @@ -1068,7 +673,7 @@ diff -ruN exercises/exercise-basic/injection2pproblem.hh exercises/solution/exer - - // \} - -- //! set the time for the time dependent boundary conditions (called from main) +- //! Set the time for the time dependent boundary conditions (called from main) - void setTime(Scalar time) - { time_ = time; } - @@ -1096,8 +701,8 @@ diff -ruN exercises/exercise-basic/injection2pproblem.hh exercises/solution/exer - -#endif diff -ruN exercises/exercise-basic/params.input exercises/solution/exercise-basic/params.input ---- exercises/exercise-basic/params.input 2024-07-08 09:06:59.581159992 +0200 -+++ exercises/solution/exercise-basic/params.input 2024-07-08 09:06:59.613159864 +0200 +--- exercises/exercise-basic/params.input 2023-03-31 13:48:56.158686713 +0200 ++++ exercises/solution/exercise-basic/params.input 2024-07-08 11:41:07.994134860 +0200 @@ -1,6 +1,6 @@ [TimeLoop] DtInitial = 3600 # in seconds @@ -1118,90 +723,8 @@ diff -ruN exercises/exercise-basic/params.input exercises/solution/exercise-basi +SolidDensity = 2700 # solid density of granite +SolidThermalConductivity = 2.8 # solid thermal conducitivity of granite +SolidHeatCapacity = 790 # solid heat capacity of granite -diff -ruN exercises/exercise-basic/properties2p2c.hh exercises/solution/exercise-basic/properties2p2c.hh ---- exercises/exercise-basic/properties2p2c.hh 2024-07-08 09:06:59.581159992 +0200 -+++ exercises/solution/exercise-basic/properties2p2c.hh 1970-01-01 01:00:00.000000000 +0100 -@@ -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 diff -ruN exercises/exercise-basic/properties2p.hh exercises/solution/exercise-basic/properties2p.hh ---- exercises/exercise-basic/properties2p.hh 2024-07-08 09:06:59.581159992 +0200 +--- exercises/exercise-basic/properties2p.hh 2024-07-08 11:41:07.986134888 +0200 +++ exercises/solution/exercise-basic/properties2p.hh 1970-01-01 01:00:00.000000000 +0100 @@ -1,75 +0,0 @@ -// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- @@ -1280,8 +803,8 @@ diff -ruN exercises/exercise-basic/properties2p.hh exercises/solution/exercise-b - -#endif diff -ruN exercises/exercise-basic/properties2pni.hh exercises/solution/exercise-basic/properties2pni.hh ---- exercises/exercise-basic/properties2pni.hh 2024-07-08 09:06:59.581159992 +0200 -+++ exercises/solution/exercise-basic/properties2pni.hh 2024-07-08 09:06:59.613159864 +0200 +--- exercises/exercise-basic/properties2pni.hh 2024-07-08 11:41:07.986134888 +0200 ++++ exercises/solution/exercise-basic/properties2pni.hh 2024-07-08 11:41:07.994134860 +0200 @@ -31,18 +31,14 @@ #include <dumux/porousmediumflow/2p/model.hh> #include <dumux/material/fluidsystems/h2on2.hh> @@ -1315,9 +838,9 @@ diff -ruN exercises/exercise-basic/properties2pni.hh exercises/solution/exercise } // end namespace Dumux::Properties diff -ruN exercises/exercise-basic/README.md exercises/solution/exercise-basic/README.md ---- exercises/exercise-basic/README.md 2024-07-16 16:55:12.734099365 +0200 +--- exercises/exercise-basic/README.md 2024-07-16 17:38:07.971854664 +0200 +++ exercises/solution/exercise-basic/README.md 1970-01-01 01:00:00.000000000 +0100 -@@ -1,105 +0,0 @@ +@@ -1,103 +0,0 @@ -# Exercise Basics (DuMuX course) -<br> - @@ -1332,7 +855,7 @@ diff -ruN exercises/exercise-basic/README.md exercises/solution/exercise-basic/R - -* 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 @@ -1340,11 +863,10 @@ diff -ruN exercises/exercise-basic/README.md exercises/solution/exercise-basic/R - -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` - @@ -1358,23 +880,22 @@ diff -ruN exercises/exercise-basic/README.md exercises/solution/exercise-basic/R -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> @@ -1392,7 +913,7 @@ diff -ruN exercises/exercise-basic/README.md exercises/solution/exercise-basic/R - SOURCES 2pnimain.cc) -``` - --* Test that everything compiles without error +-* In the respective build-cmake folder, test that everything compiles without error - -```bash -make # should rerun cmake @@ -1413,14 +934,14 @@ diff -ruN exercises/exercise-basic/README.md exercises/solution/exercise-basic/R - -* The following set-up should be realized: - -- __Boundary conditions:__ Dirichlet conditions at the left boundary. For the primary variable __temperature__ use a varying temperature of <br/> +- __Initial conditions:__ For the primary variable __temperature__ use a varying temperature of <br/> -$`\displaystyle T(y) = 283~\text{K} + 0.03~\frac{\text{K}}{\text{m}} \cdot \left( d_\text{aquifer} - y \right)`$, <br/> -with the aquifer depth --$\displaystyle d_\text{aquifer}=2700~\text{m}$. Add an energy flux at the injection point of N$_2$. Assign Neumann no-flow for the energy balance to the rest of the boundaries. +-$\displaystyle d_\text{aquifer}=2700~\text{m}$. Additionally, add a subdomain (20 < x < 30, 5 < y < 35), where you assign a constant initial temperature of 380 K. - -- __Initial conditions:__ The same temperature gradient as in the boundary conditions with an exception in the subdomain (20 < x < 30, 5 < y < 35), where you assign a constant initial temperature of 380 K. +- __Boundary conditions:__ Dirichlet boundary conditions at the left boundary with the same temperature gradient as in the initial conditions. For the Neumann conditions, assign an energy flux at the injection point of N$_2$ and no-flow conditions for the energy balance to the rest of the boundaries. - -<img src="https://git.iws.uni-stuttgart.de/dumux-repositories/dumux-course/raw/master/exercises/extradoc/exercise1_nonisothermal.png" width="800"> - -The non-isothermal model requires additional parameters like the thermal conductivity of the solid component. They are already implemented and set in `params.input`, you just need to _uncomment_ them. -\ Kein Zeilenumbruch am Dateiende. +\ No newline at end of file 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..ecd78ecfe82dccd2c101be4ba9619da1499e74f0 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> @@ -72,7 +70,7 @@ dumux_add_test(NAME exercise_basic_2pni SOURCES 2pnimain.cc) ``` -* Test that everything compiles without error +* In the respective build-cmake folder, test that everything compiles without error ```bash make # should rerun cmake @@ -93,12 +91,12 @@ Look for comments containing * The following set-up should be realized: - __Boundary conditions:__ Dirichlet conditions at the left boundary. For the primary variable __temperature__ use a varying temperature of <br/> + __Initial conditions:__ For the primary variable __temperature__ use a varying temperature of <br/> $`\displaystyle T(y) = 283~\text{K} + 0.03~\frac{\text{K}}{\text{m}} \cdot \left( d_\text{aquifer} - y \right)`$, <br/> with the aquifer depth -$\displaystyle d_\text{aquifer}=2700~\text{m}$. Add an energy flux at the injection point of N$_2$. Assign Neumann no-flow for the energy balance to the rest of the boundaries. +$\displaystyle d_\text{aquifer}=2700~\text{m}$. Additionally, add a subdomain (20 < x < 30, 5 < y < 35), where you assign a constant initial temperature of 380 K. - __Initial conditions:__ The same temperature gradient as in the boundary conditions with an exception in the subdomain (20 < x < 30, 5 < y < 35), where you assign a constant initial temperature of 380 K. + __Boundary conditions:__ Dirichlet boundary conditions at the left boundary with the same temperature gradient as in the initial conditions. For the Neumann conditions, assign an energy flux at the injection point of N$_2$ and no-flow conditions for the energy balance to the rest of the boundaries. <img src="https://git.iws.uni-stuttgart.de/dumux-repositories/dumux-course/raw/master/exercises/extradoc/exercise1_nonisothermal.png" width="800"> 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/injection2pniproblem.hh b/exercises/exercise-basic/injection2pniproblem.hh index 76e52cc1c2e0b069f60860cee76f7e4ebed86b6c..a08be51ef39587409027233569e376795fdd19c0 100644 --- a/exercises/exercise-basic/injection2pniproblem.hh +++ b/exercises/exercise-basic/injection2pniproblem.hh @@ -36,7 +36,7 @@ namespace Dumux { * \ingroup TwoPModel * \ingroup ImplicitTestProblems * \brief Gas injection problem where a gas (here nitrogen) is injected into a fully - * water saturated medium. During buoyancy driven upward migration the gas + * water saturated medium. During buoyancy-driven upward migration the gas * passes a high temperature area. * * The domain is sized 60 m times 40 m. @@ -75,7 +75,7 @@ public: Injection2PNIProblem(std::shared_ptr<const GridGeometry> gridGeometry) : ParentType(gridGeometry) { - // initialize the tables of the fluid system + // Initialize the tables of the fluid system FluidSystem::init(/*tempMin=*/273.15, /*tempMax=*/423.15, /*numTemp=*/50, @@ -87,9 +87,9 @@ public: // 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 + // depth of the aquifer, unit: m aquiferDepth_ = getParam<Scalar>("Problem.AquiferDepth"); - // the duration of the injection, units: second + // the duration of the injection, unit: seconds injectionDuration_ = getParam<Scalar>("Problem.InjectionDuration"); } @@ -116,11 +116,10 @@ public: bcTypes.setAllNeumann(); /*! - * TODO:dumux-course-task: - * dumux-course-task: - * set Dirichlet conditions for the energy equation on the left boundary - * and Neumann everywhere else - * think about: is there anything necessary to do here? + * TODO:dumux-course-task 4: + * Set Dirichlet conditions for the energy equation on the left boundary + * and Neumann everywhere else. + * Think about: is there anything necessary to do here? */ return bcTypes; @@ -137,10 +136,9 @@ public: return initialAtPos(globalPos); /*! - * TODO:dumux-course-task: - * dumux-course-task: - * set Dirichlet conditions for the energy equation on the left boundary - * think about: is there anything necessary to do here? + * TODO:dumux-course-task 4: + * Set Dirichlet conditions for the energy equation on the left boundary. + * Think about: is there anything necessary to do here? */ } @@ -159,14 +157,14 @@ public: if (injectionActive() && onInjectionBoundary(globalPos)) { // inject nitrogen. negative values mean injection - // units kg/(s*m^2) + // unit: kg/(s*m^2) values[Indices::conti0EqIdx + FluidSystem::N2Idx] = -1e-4; values[Indices::conti0EqIdx + FluidSystem::H2OIdx] = 0.0; /*! - * TODO:dumux-course-task: + * TODO:dumux-course-task 4: * Set Neumann noflow conditions for the energy equation everywhere else except the left boundary. - * Additionally, consider the energy flux at the injection point which is equal to the product of the respective mass flux and the matching enthalpy. Use the function `gasEnthalpy(...)` from the N2 component to access the necessary enthalpy. + * Additionally, consider the energy flux at the injection point which is equal to the product of the respective mass flux and the matching enthalpy. Use the function `gasEnthalpy(temperature,pressure)` from the N2 component to access the necessary enthalpy. * hint: use `Indices::energyEqIdx` to access the entry belonging to the energy flux. */ } @@ -205,16 +203,16 @@ public: values[Indices::saturationIdx] = 0.0; /*! - * TODO:dumux-course-task: - * set a temperature gradient of 0.03 K per m beginning at 283 K here. - * Hint: you can use aquiferDepth_ and the globalPos similar to the pressure gradient - * use globalPos[0] and globalPos[1] to implement the high temperature lens with 380 K + * TODO:dumux-course-task 4: + * Set a temperature gradient of 0.03 K per m beginning at 283 K here. + * Hint: you can use aquiferDepth_ and the globalPos similar to the pressure gradient. + * Use globalPos[0] and globalPos[1] to implement the high temperature lens with 380 K * Hint : use Indices::temperatureIdx to address the initial values for temperature */ return values; } - //! set the time for the time dependent boundary conditions (called from main) + //! Set the time for the time dependent boundary conditions (called from main) void setTime(Scalar time) { time_ = time; } diff --git a/exercises/exercise-basic/injection2pproblem.hh b/exercises/exercise-basic/injection2pproblem.hh index 01e5d0bb419bf751d3f5229f91305e9fa9ad2787..9ef1163c7d0cf63e2070f0d17b9e048b566fe55f 100644 --- a/exercises/exercise-basic/injection2pproblem.hh +++ b/exercises/exercise-basic/injection2pproblem.hh @@ -75,7 +75,7 @@ public: Injection2PProblem(std::shared_ptr<const GridGeometry> gridGeometry) : ParentType(gridGeometry) { - // initialize the tables of the fluid system + // Initialize the tables of the fluid system FluidSystem::init(/*tempMin=*/273.15, /*tempMax=*/423.15, /*numTemp=*/50, @@ -83,13 +83,13 @@ public: /*pMax=*/30e6, /*numP=*/300); - // name of the problem and output file + // 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 + // Depth of the aquifer, unit: m aquiferDepth_ = getParam<Scalar>("Problem.AquiferDepth"); - // the duration of the injection, units: second + // The duration of the injection, unit: seconds injectionDuration_ = getParam<Scalar>("Problem.InjectionDuration"); } @@ -149,8 +149,8 @@ public: // than using <= or >= as it is robust with regard to imprecision introduced by rounding errors. if (injectionActive() && onInjectionBoundary(globalPos)) { - // inject nitrogen. negative values mean injection - // units kg/(s*m^2) + // Inject nitrogen. Negative values mean injection + // unit: kg/(s*m^2) values[Indices::conti0EqIdx + FluidSystem::N2Idx] = -1e-4; values[Indices::conti0EqIdx + FluidSystem::H2OIdx] = 0.0; } @@ -194,7 +194,7 @@ public: // \} - //! set the time for the time dependent boundary conditions (called from main) + //! Set the time for the time dependent boundary conditions (called from main) void setTime(Scalar time) { time_ = time; } 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 diff --git a/exercises/solution/exercise-basic/injection2pniproblem.hh b/exercises/solution/exercise-basic/injection2pniproblem.hh index 5a61045bc1e83102dcf4fd7688ebca3b9b048624..7cb3c10bb66320268a87954b9e38329ae22a04e9 100644 --- a/exercises/solution/exercise-basic/injection2pniproblem.hh +++ b/exercises/solution/exercise-basic/injection2pniproblem.hh @@ -35,7 +35,7 @@ namespace Dumux { * \ingroup TwoPModel * \ingroup ImplicitTestProblems * \brief Gas injection problem where a gas (here nitrogen) is injected into a fully - * water saturated medium. During buoyancy driven upward migration the gas + * water saturated medium. During buoyancy-driven upward migration the gas * passes a high temperature area. * * The domain is sized 60 m times 40 m. @@ -76,7 +76,7 @@ public: Injection2PNIProblem(std::shared_ptr<const GridGeometry> gridGeometry) : ParentType(gridGeometry) { - // initialize the tables of the fluid system + // Initialize the tables of the fluid system FluidSystem::init(/*tempMin=*/273.15, /*tempMax=*/423.15, /*numTemp=*/50, @@ -84,13 +84,13 @@ public: /*pMax=*/30e6, /*numP=*/300); - // name of the problem and output file + // 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 + // Depth of the aquifer, unit: m aquiferDepth_ = getParam<Scalar>("Problem.AquiferDepth"); - // the duration of the injection, units: second + // The duration of the injection, unit: seconds injectionDuration_ = getParam<Scalar>("Problem.InjectionDuration"); } @@ -144,15 +144,18 @@ public: // if we are inside the injection zone set inflow Neumann boundary conditions if (injectionActive() && onInjectionBoundary(globalPos)) { - - // inject nitrogen. negative values mean injection - // units kg/(s*m^2) - values[Indices::conti0EqIdx + FluidSystem::N2Idx] = -1e-4; + const Scalar injectionRate = -1e-4; + + // inject nitrogen. Negative values mean injection + // unit: kg/(s*m^2) + values[Indices::conti0EqIdx + FluidSystem::N2Idx] = injectionRate; values[Indices::conti0EqIdx + FluidSystem::H2OIdx] = 0.0; // energy fluxes are always mass specific - //units W/(m^2) - values[Indices::energyEqIdx] = -1e-4 /*kg/(m^2 s)*/ *N2::gasEnthalpy(injectionTemperature(), 1e5 /*Pa*/)/*J/kg*/; + // unit: W/(m^2) + const Scalar temperatureAtInjection = initialAtPos(globalPos)[Indices::temperatureIdx];/*K*/ + const Scalar pressureAtInjection = initialAtPos(globalPos)[Indices::pressureIdx];/*Pa*/ + values[Indices::energyEqIdx] = injectionRate /*kg/(m^2 s)*/ *N2::gasEnthalpy(temperatureAtInjection, pressureAtInjection)/*J/kg*/; } return values; @@ -190,7 +193,7 @@ public: values[Indices::temperatureIdx] = 283.0 + (aquiferDepth_ - globalPos[1])*0.03; if (globalPos[0] > 20 - eps_ && globalPos[0] < 30 + eps_ && globalPos[1] > 5 - eps_ && globalPos[1] < 35 + eps_) - values[Indices::temperatureIdx] = injectionTemperature(); + values[Indices::temperatureIdx] = 380.0; return values; } @@ -198,7 +201,7 @@ public: //! 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_; } @@ -211,12 +214,6 @@ public: && globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; } - //! Set injection temperature for the domain - Scalar injectionTemperature() const - { - return 380.0; - } - private: static constexpr Scalar eps_ = 1e-6; std::string name_; //! Problem name