diff --git a/dumux/freeflow/CMakeLists.txt b/dumux/freeflow/CMakeLists.txt index 6f831c047117ec2d6b8df38d753709f959b2768f..344282fb7e130e514a6229f4d3a2f1b61b151db9 100644 --- a/dumux/freeflow/CMakeLists.txt +++ b/dumux/freeflow/CMakeLists.txt @@ -2,3 +2,4 @@ add_subdirectory("navierstokes") add_subdirectory("navierstokesnc") add_subdirectory("nonisothermal") add_subdirectory("rans") +add_subdirectory("ransnc") diff --git a/dumux/freeflow/ransnc/CMakeLists.txt b/dumux/freeflow/ransnc/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..c50a949f4e207d7ec35fbf1d465f324d1d85e68f --- /dev/null +++ b/dumux/freeflow/ransnc/CMakeLists.txt @@ -0,0 +1,6 @@ +#install headers +install(FILES +model.hh +volumevariables.hh +vtkoutputfields.hh +DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/freeflow/ransnc) diff --git a/dumux/freeflow/ransnc/model.hh b/dumux/freeflow/ransnc/model.hh new file mode 100644 index 0000000000000000000000000000000000000000..f3c83b064e20799ac9cac9375cb9835d5f95cdf3 --- /dev/null +++ b/dumux/freeflow/ransnc/model.hh @@ -0,0 +1,145 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see <http://www.gnu.org/licenses/>. * + *****************************************************************************/ +/*! + * \file + * \ingroup RANSModel + * + * \brief A single-phase, multi-component isothermal Navier-Stokes model + * + * \copydoc RANSModel + * * + * The system is closed by a <B> component mass/mole balance equation </B> for each component \f$\kappa\f$: + * \f[ + * \frac{\partial \left(\varrho X^\kappa\right)}{\partial t} + * + \nabla \cdot \left( \varrho {\boldsymbol{v}} X^\kappa + * - (D^\kappa + D_\text{t}) \varrho \frac{M^\kappa}{M} \textbf{grad}\, x^\kappa \right) + * - q^\kappa = 0 + * \f] + * + * Alternatively, one component balance equation can be replace by a <B> total mass/mole balance equation </B>: + * \f[ + * \frac{\partial \varrho_g}{\partial t} + * + \nabla \cdot \left( + * \varrho {\boldsymbol{v}} + * - \sum_\kappa (D^\kappa + D_\text{t}) \varrho \frac{M^\kappa}{M} \textbf{grad}\, x^\kappa + * \right) + * - q = 0 + * \f] + * + * The eddy diffusivity \$[ D_\text{t} \$] is related to the eddy viscosity by the turbulent + * Schmidt number: + * \f[ D_\text{t} = \frac{\nu_\text{t}}{\mathrm{Sc}_\text{t}} \f] + * + * So far, only the staggered grid spatial discretization (for structured grids) is available. + */ + +#ifndef DUMUX_RANS_NC_MODEL_HH +#define DUMUX_RANS_NC_MODEL_HH + + +// TODO: remove unused headers +#include <dumux/common/properties.hh> + +#include <dumux/freeflow/rans/model.hh> +#include <dumux/freeflow/rans/zeroeq/model.hh> +#include <dumux/freeflow/navierstokesnc/model.hh> +#include <dumux/freeflow/nonisothermal/model.hh> +#include <dumux/freeflow/nonisothermal/indices.hh> +#include <dumux/freeflow/nonisothermal/vtkoutputfields.hh> +#include <dumux/discretization/fickslaw.hh> +#include <dumux/discretization/fourierslaw.hh> + +#include "volumevariables.hh" +#include "vtkoutputfields.hh" + +#include <dumux/assembly/staggeredlocalresidual.hh> +#include <dumux/material/fluidsystems/gasphase.hh> +#include <dumux/material/fluidsystems/liquidphase.hh> + +#include <dumux/material/fluidstates/compositional.hh> + +namespace Dumux { + +/*! + * \ingroup RANSModel + * \brief Traits for the Reynolds-averaged Navier-Stokes multi-component model + */ +template<int dim, int nComp> +struct RANSNCModelTraits : NavierStokesNCModelTraits<dim, nComp> +{ }; + +/////////////////////////////////////////////////////////////////////////// +// properties for the single-phase, multi-component Reynolds-averaged Navier-Stokes model +/////////////////////////////////////////////////////////////////////////// +namespace Properties { + +////////////////////////////////////////////////////////////////// +// Type tags +////////////////////////////////////////////////////////////////// + +//! The type tag for the single-phase, multi-component isothermal Reynolds-averaged Navier-Stokes model +NEW_TYPE_TAG(RANSNC, INHERITS_FROM(RANS, NavierStokesNC)); +NEW_TYPE_TAG(ZeroEqNC, INHERITS_FROM(RANSNC, NavierStokesNC)); + +// //! The type tag for the single-phase, multi-component non-isothermal Reynolds-averaged Navier-Stokes model +// NEW_TYPE_TAG(RANSNCNI, INHERITS_FROM(RANS, NavierStokesNC)); + +/////////////////////////////////////////////////////////////////////////// +// default property values +/////////////////////////////////////////////////////////////////////////// + +//! The volume variables +SET_TYPE_PROP(RANSNC, VolumeVariables, RANSNCVolumeVariables<TypeTag>); + +//! The specific vtk output fields +SET_PROP(RANSNC, VtkOutputFields) +{ +private: + using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); + using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem); + static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx); +public: + using type = RANSNCVtkOutputFields<FVGridGeometry, FluidSystem, phaseIdx>; +}; + +////////////////////////////////////////////////////////////////////////// +// Property values for non-isothermal multi-component Reynolds-averaged Navier-Stokes model +////////////////////////////////////////////////////////////////////////// + +// //! The non-isothermal vtk output fields +// SET_PROP(RANSNCNI, VtkOutputFields) +// { +// private: +// using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); +// using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem); +// static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx); +// using IsothermalFields = RANSNCVtkOutputFields<FVGridGeometry, FluidSystem, phaseIdx>; +// public: +// using type = NavierStokesNonIsothermalVtkOutputFields<IsothermalFields>; +// }; +// +// //! Use Fourier's Law as default heat conduction type +// SET_TYPE_PROP(RANSNCNI, HeatConductionType, FouriersLaw<TypeTag>); + +// \} +} // end namespace Properties +} // end namespace Dumux + + +#endif diff --git a/dumux/freeflow/ransnc/volumevariables.hh b/dumux/freeflow/ransnc/volumevariables.hh new file mode 100644 index 0000000000000000000000000000000000000000..91b3570715c547d7c56617a244a70395cb0b12de --- /dev/null +++ b/dumux/freeflow/ransnc/volumevariables.hh @@ -0,0 +1,113 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see <http://www.gnu.org/licenses/>. * + *****************************************************************************/ +/*! + * \file + * \ingroup RANSNCModel + * + * \copydoc Dumux::RANSNCVolumeVariables + */ +#ifndef DUMUX_RANS_NC_VOLUMEVARIABLES_HH +#define DUMUX_RANS_NC_VOLUMEVARIABLES_HH + +#include <dumux/common/properties.hh> + +#include <dumux/freeflow/rans/volumevariables.hh> +#include <dumux/freeflow/navierstokesnc/volumevariables.hh> + +namespace Dumux { + +/*! + * \ingroup RANSNCModel + * \brief Volume variables for the single-phase, multi-component Reynolds-averaged Navier-Stokes model. + */ +template <class TypeTag> +class RANSNCVolumeVariables + : public RANSVolumeVariables<TypeTag>, + public NavierStokesNCVolumeVariables<TypeTag> +{ + using ParentTypeSinglePhase = RANSVolumeVariables<TypeTag>; + using ParentTypeCompositional = NavierStokesNCVolumeVariables<TypeTag>; + using Implementation = typename GET_PROP_TYPE(TypeTag, VolumeVariables); + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView; + using SubControlVolume = typename FVElementGeometry::SubControlVolume; + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using Element = typename GridView::template Codim<0>::Entity; + +public: + + /*! + * \brief Update all quantities for a given control volume + * + * \param elemSol A vector containing all primary variables connected to the element + * \param problem The object specifying the problem which ought to + * be simulated + * \param element An element which contains part of the control volume + * \param scv The sub-control volume + */ + template<class ElementSolution> + void update(const ElementSolution &elemSol, + const Problem &problem, + const Element &element, + const SubControlVolume& scv) + { + this->extrusionFactor_ = problem.extrusionFactor(element, scv, elemSol); + + ParentTypeSinglePhase::update(elemSol, problem, element, scv, this->fluidState_); + ParentTypeCompositional::update(elemSol, problem, element, scv, this->fluidState_); + + static const int turbulentSchmidtNumber + = getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), + "RANS.TurbulentSchmidtNumber", 1.0); + eddyDiffusivity_ = ParentTypeSinglePhase::kinematicViscosity() + / turbulentSchmidtNumber; + } + + /*! + * \brief Returns the eddy diffusivity \f$\mathrm{[m^2/s]}\f$ + */ + Scalar eddyDiffusivity() const + { + eddyDiffusivity_; + } + + /*! + * \brief Returns the effective diffusion coefficient \f$\mathrm{[m^2/s]}\f$ + */ + Scalar effectiveDiffusivity(int pIdx, int compIdx) const + { + return ParentTypeCompositional::diffusionCoefficient(pIdx, compIdx) + + eddyDiffusivity(); + } + +protected: + + Implementation &asImp_() + { return *static_cast<Implementation*>(this); } + + const Implementation &asImp_() const + { return *static_cast<const Implementation*>(this); } + + Scalar eddyDiffusivity_; +}; + +} // end namespace Dumux + +#endif diff --git a/dumux/freeflow/ransnc/vtkoutputfields.hh b/dumux/freeflow/ransnc/vtkoutputfields.hh new file mode 100644 index 0000000000000000000000000000000000000000..14280a2b3a3ef947f37d78dcd3b9effa2a635a33 --- /dev/null +++ b/dumux/freeflow/ransnc/vtkoutputfields.hh @@ -0,0 +1,54 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see <http://www.gnu.org/licenses/>. * + *****************************************************************************/ +/*! + * \file + * \ingroup RANSNCModel + * \copydoc Dumux::RANSNCVtkOutputFields + */ +#ifndef DUMUX_RANS_NC_VTK_OUTPUT_FIELDS_HH +#define DUMUX_RANS_NC_VTK_OUTPUT_FIELDS_HH + +#include <dumux/freeflow/rans/vtkoutputfields.hh> +#include <dumux/freeflow/navierstokesnc/vtkoutputfields.hh> + +namespace Dumux +{ + +/*! + * \ingroup RANSNCModel + * \brief Adds vtk output fields specific to the RANSNC model + */ +template<class FVGridGeometry, class FluidSystem, int phaseIdx> +class RANSNCVtkOutputFields +{ + +public: + //! Initialize the RANSNC specific vtk output fields. + template <class VtkOutputModule> + static void init(VtkOutputModule& vtk) + { + RANSVtkOutputFields<FVGridGeometry>::init(vtk); + NavierStokesNCVtkOutputFields<FVGridGeometry, FluidSystem, phaseIdx>::init(vtk); + vtk.addVolumeVariable([](const auto& v){ return v.eddyDiffusivity(); }, "D_t"); + } +}; + +} // end namespace Dumux + +#endif diff --git a/test/freeflow/CMakeLists.txt b/test/freeflow/CMakeLists.txt index 5fa27535b21e7a14de55bba0527bca6371711eda..ee2d61ab69e4cd74c942fbc4436f472594e821a1 100644 --- a/test/freeflow/CMakeLists.txt +++ b/test/freeflow/CMakeLists.txt @@ -1,3 +1,4 @@ add_subdirectory("navierstokes") add_subdirectory("navierstokesnc") add_subdirectory("rans") +add_subdirectory("ransnc") diff --git a/test/freeflow/ransnc/CMakeLists.txt b/test/freeflow/ransnc/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..abd0de20b7b70ed9a29228fea367ca7783a16314 --- /dev/null +++ b/test/freeflow/ransnc/CMakeLists.txt @@ -0,0 +1,10 @@ +add_input_file_links() + +dune_add_test(NAME test_channel_zeroeq2c + SOURCES test_channel_zeroeq.cc + CMAKE_GUARD HAVE_UMFPACK + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMD_ARGS --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/test_channel_zeroeq2c.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_channel_zeroeq2c-00013.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_channel_zeroeq2c test_channel_zeroeq.input") diff --git a/test/freeflow/ransnc/channeltestproblem.hh b/test/freeflow/ransnc/channeltestproblem.hh new file mode 100644 index 0000000000000000000000000000000000000000..79c7a7a1e8098320b676eb5d8e1fb05fa327a843 --- /dev/null +++ b/test/freeflow/ransnc/channeltestproblem.hh @@ -0,0 +1,324 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see <http://www.gnu.org/licenses/>. * + *****************************************************************************/ +/*! + * \file + * \ingroup RANSNCTests + * \brief Channel flow test for the multi-component staggered grid Reynolds-averaged Navier-Stokes model + */ +#ifndef DUMUX_RANS_NC_TEST_PROBLEM_HH +#define DUMUX_RANS_NC_TEST_PROBLEM_HH + +#include <dumux/material/fluidsystems/liquidphase.hh> +#include <dumux/material/components/simpleh2o.hh> +#include <dumux/material/components/constant.hh> +#include <dumux/material/fluidsystems/h2oair.hh> + +#include <dumux/discretization/staggered/freeflow/properties.hh> + +#include <dumux/freeflow/rans/problem.hh> +#include <dumux/freeflow/ransnc/model.hh> + +namespace Dumux +{ +template <class TypeTag> +class ChannelNCTestProblem; + +namespace Properties +{ + +#if !NONISOTHERMAL +NEW_TYPE_TAG(ChannelNCTestTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, ZeroEqNC)); +#else +NEW_TYPE_TAG(ChannelNCTestTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, RANSNCNI)); +#endif + +NEW_PROP_TAG(FluidSystem); + +// Select the fluid system +SET_TYPE_PROP(ChannelNCTestTypeTag, FluidSystem, + FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>); + +SET_INT_PROP(ChannelNCTestTypeTag, PhaseIdx, + GET_PROP_TYPE(TypeTag, FluidSystem)::nPhaseIdx); + +SET_INT_PROP(ChannelNCTestTypeTag, ReplaceCompEqIdx, 0); + +// Set the grid type +SET_TYPE_PROP(ChannelNCTestTypeTag, Grid, Dune::YaspGrid<2>); + +// Set the problem property +SET_TYPE_PROP(ChannelNCTestTypeTag, Problem, Dumux::ChannelNCTestProblem<TypeTag> ); + +SET_BOOL_PROP(ChannelNCTestTypeTag, EnableFVGridGeometryCache, true); + +SET_BOOL_PROP(ChannelNCTestTypeTag, EnableGridFluxVariablesCache, true); +SET_BOOL_PROP(ChannelNCTestTypeTag, EnableGridVolumeVariablesCache, true); + +// Enable gravity +SET_BOOL_PROP(ChannelNCTestTypeTag, UseMoles, true); +} + +/*! + * \ingroup RANSNCTests + * \brief Test problem for the one-phase model. + * \todo doc me! + */ +template <class TypeTag> +class ChannelNCTestProblem : public RANSProblem<TypeTag> +{ + using ParentType = RANSProblem<TypeTag>; + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem); + + // copy some indices for convenience + using Indices = typename GET_PROP_TYPE(TypeTag, Indices); + enum { dimWorld = GridView::dimensionworld }; + enum { + massBalanceIdx = Indices::massBalanceIdx, + transportEqIdx = 1, + momentumBalanceIdx = Indices::momentumBalanceIdx, + pressureIdx = Indices::pressureIdx, + velocityXIdx = Indices::velocityXIdx, + velocityYIdx = Indices::velocityYIdx, +#if NONISOTHERMAL + temperatureIdx = Indices::temperatureIdx, + energyBalanceIdx = Indices::energyBalanceIdx, +#endif + transportCompIdx = 1/*FluidSystem::wCompIdx*/ + }; + + using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes); + + using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); + + using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>; + + using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables); + using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector); + + using TimeLoopPtr = std::shared_ptr<CheckPointTimeLoop<Scalar>>; + using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables); + using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector); + +public: + ChannelNCTestProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry) + : ParentType(fvGridGeometry), eps_(1e-6) + { + inletVelocity_ = getParam<Scalar>("Problem.InletVelocity"); + FluidSystem::init(); + deltaP_.resize(this->fvGridGeometry().numCellCenterDofs()); + } + + /*! + * \name Problem parameters + */ + // \{ + + + bool shouldWriteRestartFile() const + { + return false; + } + + /*! + * \brief Return the temperature within the domain in [K]. + * + * This problem assumes a temperature of 10 degrees Celsius. + */ + Scalar temperature() const + { return 273.15 + 10; } // 10C + + /*! + * \brief Return the sources within the domain. + * + * \param globalPos The global position + */ + NumEqVector sourceAtPos(const GlobalPosition &globalPos) const + { + return NumEqVector(0.0); + } + // \} + /*! + * \name Boundary conditions + */ + // \{ + + /*! + * \brief Specifies which kind of boundary condition should be + * used for which equation on a given boundary control volume. + * + * \param globalPos The position of the center of the finite volume + */ + BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const + { + BoundaryTypes values; + + if(isInlet(globalPos)) + { + values.setDirichlet(momentumBalanceIdx); + values.setOutflow(massBalanceIdx); + values.setDirichlet(transportEqIdx); +#if NONISOTHERMAL + values.setDirichlet(energyBalanceIdx); +#endif + } + else if(isOutlet(globalPos)) + { + values.setOutflow(momentumBalanceIdx); + values.setDirichlet(massBalanceIdx); + values.setOutflow(transportEqIdx); +#if NONISOTHERMAL + values.setOutflow(energyBalanceIdx); +#endif + } + + else + { + // set Dirichlet values for the velocity everywhere + values.setDirichlet(momentumBalanceIdx); + values.setOutflow(massBalanceIdx); + values.setOutflow(transportEqIdx); +#if NONISOTHERMAL + values.setOutflow(energyBalanceIdx); +#endif + } + + return values; + } + + /*! + * \brief Evaluate the boundary conditions for a dirichlet + * control volume. + * + * \param globalPos The center of the finite volume which ought to be set. + */ + PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const + { + PrimaryVariables values = initialAtPos(globalPos); + + // give the system some time so that the pressure can equilibrate, then start the injection of the tracer + if(isInlet(globalPos)) + { + if(time() >= 10.0 || inletVelocity_ < eps_) + { + values[transportCompIdx] = 1e-3; +#if NONISOTHERMAL + values[temperatureIdx] = 293.15; +#endif + } + } + + return values; + } + + // \} + + /*! + * \name Volume terms + */ + // \{ + + /*! + * \brief Evaluate the initial value for a control volume. + * + * \param globalPos The global position + */ + PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const + { + PrimaryVariables values; + values[pressureIdx] = 1.1e+5; + values[transportCompIdx] = 0.0; +#if NONISOTHERMAL + values[temperatureIdx] = 283.15; +#endif + + // parabolic velocity profile + values[velocityXIdx] = inletVelocity_*(globalPos[1] - this->fvGridGeometry().bBoxMin()[1])*(this->fvGridGeometry().bBoxMax()[1] - globalPos[1]) + / (0.25*(this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1])*(this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1])); + + values[velocityYIdx] = 0.0; + + return values; + } + + /*! + * \brief Adds additional VTK output data to the VTKWriter. Function is called by the output module on every write. + * + * \param gridVariables The grid variables + * \param sol The solution vector + */ + void calculateDeltaP(const GridVariables& gridVariables, const SolutionVector& sol) + { + for (const auto& element : elements(this->fvGridGeometry().gridView())) + { + auto fvGeometry = localView(this->fvGridGeometry()); + fvGeometry.bindElement(element); + for (auto&& scv : scvs(fvGeometry)) + { + auto ccDofIdx = scv.dofIndex(); + + auto elemVolVars = localView(gridVariables.curGridVolVars()); + elemVolVars.bind(element, fvGeometry, sol); + + deltaP_[ccDofIdx] = elemVolVars[scv].pressure() - 1.1e5; + } + } + } + + auto& getDeltaP() const + { return deltaP_; } + + // \} + + void setTimeLoop(TimeLoopPtr timeLoop) + { + timeLoop_ = timeLoop; + } + + Scalar time() const + { + return timeLoop_->time(); + } + + bool isOnWall(const GlobalPosition& globalPos) const + { + return globalPos[0] > eps_ || globalPos[0] < this->fvGridGeometry().bBoxMax()[0] - eps_; + } + +private: + + bool isInlet(const GlobalPosition& globalPos) const + { + return globalPos[0] < eps_; + } + + bool isOutlet(const GlobalPosition& globalPos) const + { + return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_; + } + + const Scalar eps_; + Scalar inletVelocity_; + TimeLoopPtr timeLoop_; + std::vector<Scalar> deltaP_; +}; +} //end namespace + +#endif diff --git a/test/freeflow/ransnc/test_channel_zeroeq.cc b/test/freeflow/ransnc/test_channel_zeroeq.cc new file mode 100644 index 0000000000000000000000000000000000000000..eb7add37a4247571645aad701401a764eab04926 --- /dev/null +++ b/test/freeflow/ransnc/test_channel_zeroeq.cc @@ -0,0 +1,237 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see <http://www.gnu.org/licenses/>. * + *****************************************************************************/ +/*! + * \file + * + * \brief Test for the staggered grid multi-component (Navier-)Stokes model + */ + #include <config.h> + + #include <ctime> + #include <iostream> + + #include <dune/common/parallel/mpihelper.hh> + #include <dune/common/timer.hh> + #include <dune/grid/io/file/dgfparser/dgfexception.hh> + #include <dune/grid/io/file/vtk.hh> + #include <dune/istl/io.hh> + +#include "channeltestproblem.hh" + +#include <dumux/common/properties.hh> +#include <dumux/common/parameters.hh> +#include <dumux/common/valgrind.hh> +#include <dumux/common/dumuxmessage.hh> +#include <dumux/common/defaultusagemessage.hh> + +#include <dumux/linear/seqsolverbackend.hh> +#include <dumux/nonlinear/newtonsolver.hh> + +#include <dumux/assembly/staggeredfvassembler.hh> +#include <dumux/assembly/diffmethod.hh> + +#include <dumux/discretization/methods.hh> + +#include <dumux/io/staggeredvtkoutputmodule.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 += "\nPlease use the provided input files.\n"; + std::cout << errorMessageOut + << "\n"; + } +} + +int main(int argc, char** argv) try +{ + using namespace Dumux; + + // define the type tag for this problem + using TypeTag = TTAG(ChannelNCTestTypeTag); + + // 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) + using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator); + GridCreator::makeGrid(); + GridCreator::loadBalance(); + + //////////////////////////////////////////////////////////// + // run instationary non-linear problem on this grid + //////////////////////////////////////////////////////////// + + // we compute on the leaf grid view + const auto& leafGridView = GridCreator::grid().leafGridView(); + + // create the finite volume grid geometry + using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); + auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView); + fvGridGeometry->update(); + + // the problem (initial and boundary conditions) + using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + auto problem = std::make_shared<Problem>(fvGridGeometry); + + // get some time loop parameters + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + const auto tEnd = getParam<Scalar>("TimeLoop.TEnd"); + const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize"); + auto dt = getParam<Scalar>("TimeLoop.DtInitial"); + + // check if we are about to restart a previously interrupted simulation + Scalar restartTime = 0; + if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart")) + restartTime = getParam<Scalar>("TimeLoop.Restart"); + + // instantiate time loop + auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(restartTime, dt, tEnd); + timeLoop->setMaxTimeStepSize(maxDt); + problem->setTimeLoop(timeLoop); + + // the solution vector + using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector); + using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices); + typename DofTypeIndices::CellCenterIdx cellCenterIdx; + typename DofTypeIndices::FaceIdx faceIdx; + const auto numDofsCellCenter = leafGridView.size(0); + const auto numDofsFace = leafGridView.size(1); + SolutionVector x; + x[cellCenterIdx].resize(numDofsCellCenter); + x[faceIdx].resize(numDofsFace); + problem->applyInitialSolution(x); + problem->updateStaticWallProperties(); + problem->updateDynamicWallProperties(x); + auto xOld = x; + + // the grid variables + using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables); + auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry); + gridVariables->init(x, xOld); + + // intialize the vtk output module + using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields); + StaggeredVtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name()); + VtkOutputFields::init(vtkWriter); //!< Add model specific output fields + vtkWriter.addField(problem->getDeltaP(), "deltaP"); + vtkWriter.write(0.0); + + // the assembler with time loop for instationary problem + using Assembler = StaggeredFVAssembler<TypeTag, DiffMethod::numeric>; + auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop); + + // the linear solver + using LinearSolver = Dumux::UMFPackBackend; + auto linearSolver = std::make_shared<LinearSolver>(); + + // the non-linear solver + using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>; + 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(); + + problem->calculateDeltaP(*gridVariables, x); + + // update wall properties + problem->updateDynamicWallProperties(x); + + // write vtk output + vtkWriter.write(timeLoop->time()); + + // report statistics of this time step + timeLoop->reportTimeStep(); + + // set new dt as suggested by 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); + } + + return 0; +} // end main +catch (Dumux::ParameterException &e) +{ + std::cerr << std::endl << e << " ---> Abort!" << std::endl; + return 1; +} +catch (Dune::DGFException & e) +{ + std::cerr << "DGF exception thrown (" << e << + "). Most likely, the DGF file name is wrong " + "or the DGF file is corrupted, " + "e.g. missing hash at end of file or wrong number (dimensions) of entries." + << " ---> Abort!" << std::endl; + return 2; +} +catch (Dune::Exception &e) +{ + std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl; + return 3; +} +catch (...) +{ + std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl; + return 4; +} diff --git a/test/freeflow/ransnc/test_channel_zeroeq.input b/test/freeflow/ransnc/test_channel_zeroeq.input new file mode 100644 index 0000000000000000000000000000000000000000..bd58ed545c202399a01538b8ee2c5c0912c14e32 --- /dev/null +++ b/test/freeflow/ransnc/test_channel_zeroeq.input @@ -0,0 +1,42 @@ +[TimeLoop] +DtInitial = 1 # [s] +TEnd = 1000 # [s] + +[Grid] +UpperRight = 5 1 +Cells = 50 10 + +[Problem] +Name = test_channel_zeroeq # name passed to the output routines +InletVelocity = 0 +EnableGravity = false + +[ Newton ] +MaxSteps = 10 +MaxRelativeShift = 1e-8 + +[Vtk] +AddVelocity = true +WriteFaceData = false + +# ADVECTION +# [TimeLoop] +# DtInitial = 5 # [s] +# TEnd = 400 # [s] +# +# [Grid] +# UpperRight = 5 1 +# Cells = 50 25 +# +# [Problem] +# Name = test_advection # name passed to the output routines +# InletVelocity = 1e-2 +# EnableGravity = false +# +# [ Newton ] +# MaxSteps = 10 +# MaxRelativeShift = 1e-8 +# +# [Vtk] +# AddVelocity = true +# WriteFaceData = false