From 6a38604e178ed7d317012a988adb9cd1cf58e16f Mon Sep 17 00:00:00 2001 From: Thomas Fetzer <thomas.fetzer@iws.uni-stuttgart.de> Date: Fri, 12 Jan 2018 14:44:16 +0100 Subject: [PATCH] [rans][doc] First basic test --- doc/handbook/dumux-handbook.bib | 12 + test/freeflow/CMakeLists.txt | 1 + test/freeflow/rans/CMakeLists.txt | 10 + test/freeflow/rans/pipelauferproblem.hh | 245 +++++++++++++++++ test/freeflow/rans/test_pipe_laufer.cc | 312 ++++++++++++++++++++++ test/freeflow/rans/test_pipe_laufer.input | 26 ++ 6 files changed, 606 insertions(+) create mode 100644 test/freeflow/rans/CMakeLists.txt create mode 100644 test/freeflow/rans/pipelauferproblem.hh create mode 100644 test/freeflow/rans/test_pipe_laufer.cc create mode 100644 test/freeflow/rans/test_pipe_laufer.input diff --git a/doc/handbook/dumux-handbook.bib b/doc/handbook/dumux-handbook.bib index 97600f24df..b37413ba0f 100644 --- a/doc/handbook/dumux-handbook.bib +++ b/doc/handbook/dumux-handbook.bib @@ -1571,3 +1571,15 @@ url={http://dx.doi.org/10.1007/s11242-015-0599-1} timestamp = {2017.05.11}, url = {http://dx.doi.org/10.1023/B:IJOT.0000022327.04529.f3} } + +@Article{Laufer1954a, + author = {John Laufer}, + title = {{The structure of turbulence in fully developed pipe flow}}, + journal = {NACA Report}, + year = {1954}, + volume = {1174}, + pages = {417--434}, + owner = {fetzer}, + timestamp = {2012.09.03}, + url = {https://ntrs.nasa.gov/search.jsp?R=19930092199} +} diff --git a/test/freeflow/CMakeLists.txt b/test/freeflow/CMakeLists.txt index a2cd437fb9..5fa27535b2 100644 --- a/test/freeflow/CMakeLists.txt +++ b/test/freeflow/CMakeLists.txt @@ -1,2 +1,3 @@ add_subdirectory("navierstokes") add_subdirectory("navierstokesnc") +add_subdirectory("rans") diff --git a/test/freeflow/rans/CMakeLists.txt b/test/freeflow/rans/CMakeLists.txt new file mode 100644 index 0000000000..a41d8ffb76 --- /dev/null +++ b/test/freeflow/rans/CMakeLists.txt @@ -0,0 +1,10 @@ +add_input_file_links() + +dune_add_test(NAME test_pipe_laufer + SOURCES test_pipe_laufer.cc + CMAKE_GUARD HAVE_UMFPACK + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMD_ARGS --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/pipe_laufer.vtu + ${CMAKE_CURRENT_BINARY_DIR}/pipe_laufer-00002.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_pipe_laufer") diff --git a/test/freeflow/rans/pipelauferproblem.hh b/test/freeflow/rans/pipelauferproblem.hh new file mode 100644 index 0000000000..933dc780bf --- /dev/null +++ b/test/freeflow/rans/pipelauferproblem.hh @@ -0,0 +1,245 @@ +// -*- 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 RANSTests + * \brief Pipe flow test for the staggered grid RANS model + * + * This test simulates is based on pipe flow experiments by + * John Laufers experiments in 1954 \cite Laufer1954a. + */ +#ifndef DUMUX_PIPE_LAUFER_PROBLEM_HH +#define DUMUX_PIPE_LAUFER_PROBLEM_HH + +#include <dumux/material/fluidsystems/gasphase.hh> +#include <dumux/material/components/air.hh> + +#include <dumux/freeflow/rans/model.hh> +#include <dumux/freeflow/rans/problem.hh> +#include <dumux/discretization/staggered/freeflow/properties.hh> + +namespace Dumux +{ +template <class TypeTag> +class PipeLauferProblem; + +namespace Properties +{ +NEW_TYPE_TAG(PipeLauferProblem, INHERITS_FROM(StaggeredFreeFlowModel, RANS)); + +// the fluid system +SET_PROP(PipeLauferProblem, FluidSystem) +{ + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using type = FluidSystems::GasPhase<Scalar, Air<Scalar> >; +}; + +// Set the grid type +SET_TYPE_PROP(PipeLauferProblem, Grid, + Dune::YaspGrid<2, Dune::TensorProductCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >); + +// Set the problem property +SET_TYPE_PROP(PipeLauferProblem, Problem, Dumux::PipeLauferProblem<TypeTag> ); + +SET_BOOL_PROP(PipeLauferProblem, EnableFVGridGeometryCache, true); + +SET_BOOL_PROP(PipeLauferProblem, EnableGridFluxVariablesCache, true); +SET_BOOL_PROP(PipeLauferProblem, EnableGridVolumeVariablesCache, true); +} + +/*! + * \ingroup NavierStokesTests + * \brief Test problem for the one-phase (Navier-) Stokes problem in a channel. + * + * This test simulates is based on pipe flow experiments by + * John Laufers experiments in 1954 \cite Laufer1954a. + */ +template <class TypeTag> +class PipeLauferProblem : public RANSProblem<TypeTag> +{ + using ParentType = RANSProblem<TypeTag>; + + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + + // copy some indices for convenience + using Indices = typename GET_PROP_TYPE(TypeTag, Indices); + enum { dimWorld = GridView::dimensionworld }; + enum { + massBalanceIdx = Indices::massBalanceIdx, + momentumBalanceIdx = Indices::momentumBalanceIdx, + pressureIdx = Indices::pressureIdx, + velocityXIdx = Indices::velocityXIdx, + velocityYIdx = Indices::velocityYIdx + }; + + 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 SourceValues = typename GET_PROP_TYPE(TypeTag, NumEqVector); + + using TimeLoopPtr = std::shared_ptr<CheckPointTimeLoop<Scalar>>; + +public: + PipeLauferProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry) + : ParentType(fvGridGeometry), eps_(1e-6) + { + inletVelocity_ = getParam<Scalar>("Problem.InletVelocity"); + } + + /*! + * \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 + */ + SourceValues sourceAtPos(const GlobalPosition &globalPos) const + { + return SourceValues(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; + + // set Dirichlet values for the velocity everywhere + values.setDirichlet(momentumBalanceIdx); + + // set a fixed pressure in one cell + if (isOutlet(globalPos)) + { + values.setDirichlet(massBalanceIdx); + values.setOutflow(momentumBalanceIdx); + } + else + values.setOutflow(massBalanceIdx); + + 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); + + if(isInlet(globalPos)) + { + values[velocityXIdx] = inletVelocity_; + } + + 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[velocityXIdx] = 0.0; + values[velocityYIdx] = 0.0; + + return values; + } + + // \} + void setTimeLoop(TimeLoopPtr timeLoop) + { + timeLoop_ = timeLoop; + if(inletVelocity_ > eps_) + timeLoop_->setCheckPoint({200.0, 210.0}); + } + + Scalar time() const + { + return timeLoop_->time(); + } + +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_; + } + + bool isWall(const GlobalPosition& globalPos) const + { + return globalPos[0] > eps_ || globalPos[0] < this->fvGridGeometry().bBoxMax()[0] - eps_; + } + + Scalar eps_; + Scalar inletVelocity_; + TimeLoopPtr timeLoop_; +}; +} //end namespace + +#endif diff --git a/test/freeflow/rans/test_pipe_laufer.cc b/test/freeflow/rans/test_pipe_laufer.cc new file mode 100644 index 0000000000..5c5f4d4a6b --- /dev/null +++ b/test/freeflow/rans/test_pipe_laufer.cc @@ -0,0 +1,312 @@ +// -*- 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 Pipe flow test for the staggered grid RANS model + * + * This test simulates is based on pipe flow experiments by + * John Laufers experiments in 1954 \cite Laufer1954a. + */ + #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 "pipelauferproblem.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/newtonmethod.hh> +#include <dumux/nonlinear/newtoncontroller.hh> + +#include <dumux/assembly/staggeredfvassembler.hh> +#include <dumux/assembly/diffmethod.hh> +#include <dumux/nonlinear/staggerednewtoncontroller.hh> + +#include <dumux/discretization/methods.hh> + +#include <dumux/io/staggeredvtkoutputmodule.hh> + +#include <dumux/freeflow/navierstokes/staggered/fluxoverplane.hh> + +/*! + * \brief Provides an interface for customizing error messages associated with + * reading in parameters. + * + * \param progName The name of the program, that was tried to be started. + * \param errorMsg The error message that was issued by the start function. + * Comprises the thing that went wrong and a general help message. + */ +void usage(const char *progName, const std::string &errorMsg) +{ +// if (errorMsg.size() > 0) { +// std::string errorMessageOut = "\nUsage: "; +// errorMessageOut += progName; +// errorMessageOut += " [options]\n"; +// errorMessageOut += errorMsg; +// errorMessageOut += "\n\nThe list of mandatory arguments for this program is:\n" +// "\t-TimeManager.TEnd End of the simulation [s] \n" +// "\t-TimeManager.DtInitial Initial timestep size [s] \n" +// "\t-Grid.File Name of the file containing the grid \n" +// "\t definition in DGF format\n" +// "\t-SpatialParams.LensLowerLeftX x-coordinate of the lower left corner of the lens [m] \n" +// "\t-SpatialParams.LensLowerLeftY y-coordinate of the lower left corner of the lens [m] \n" +// "\t-SpatialParams.LensUpperRightX x-coordinate of the upper right corner of the lens [m] \n" +// "\t-SpatialParams.LensUpperRightY y-coordinate of the upper right corner of the lens [m] \n" +// "\t-SpatialParams.Permeability Permeability of the domain [m^2] \n" +// "\t-SpatialParams.PermeabilityLens Permeability of the lens [m^2] \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(PipeLauferProblem); + + // 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 maxDivisions = getParam<int>("TimeLoop.MaxTimeStepDivisions"); + 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); + 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.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<TypeTag>; + auto linearSolver = std::make_shared<LinearSolver>(); + + // the non-linear solver + using NewtonController = StaggeredNewtonController<TypeTag>; + using NewtonMethod = Dumux::NewtonMethod<NewtonController, Assembler, LinearSolver>; + auto newtonController = std::make_shared<NewtonController>(leafGridView.comm(), timeLoop); + NewtonMethod nonLinearSolver(newtonController, assembler, linearSolver); + + // set up two planes over which fluxes are calculated + FluxOverPlane<TypeTag> flux(*assembler, x); + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using GlobalPosition = Dune::FieldVector<Scalar, GridView::dimensionworld>; + + const Scalar xMin = fvGridGeometry->bBoxMin()[0]; + const Scalar xMax = fvGridGeometry->bBoxMax()[0]; + const Scalar yMin = fvGridGeometry->bBoxMin()[1]; + const Scalar yMax = fvGridGeometry->bBoxMax()[1]; + + // The first plane shall be placed at the middle of the channel. + // If we have an odd number of cells in x-direction, there would not be any cell faces + // at the postion of the plane (which is required for the flux calculation). + // In this case, we add half a cell-width to the x-position in order to make sure that + // the cell faces lie on the plane. This assumes a regular cartesian grid. + const Scalar planePosMiddleX = xMin + 0.5*(xMax - xMin); + const int numCellsX = getParam<std::vector<int>>("Grid.Cells")[0]; + const Scalar offsetX = (numCellsX % 2 == 0) ? 0.0 : 0.5*((xMax - xMin) / numCellsX); + + const auto p0middle = GlobalPosition{planePosMiddleX + offsetX, yMin}; + const auto p1middle = GlobalPosition{planePosMiddleX + offsetX, yMax}; + flux.addPlane("middle", p0middle, p1middle); + + // The second plane is placed at the outlet of the channel. + const auto p0outlet = GlobalPosition{xMax, yMin}; + const auto p1outlet = GlobalPosition{xMax, yMax}; + flux.addPlane("outlet", p0outlet, p1outlet); + + // time loop + timeLoop->start(); do + { + // set previous solution for storage evaluations + assembler->setPreviousSolution(xOld); + + // try solving the non-linear system + for (int i = 0; i < maxDivisions; ++i) + { + // linearize & solve + auto converged = nonLinearSolver.solve(x); + + if (converged) + break; + + if (!converged && i == maxDivisions-1) + DUNE_THROW(Dune::MathError, + "Newton solver didn't converge after " + << maxDivisions + << " time-step divisions. dt=" + << timeLoop->timeStepSize() + << ".\nThe solutions of the current and the previous time steps " + << "have been saved to restart files."); + } + + // make the new solution the old solution + xOld = x; + gridVariables->advanceTimeStep(); + + // advance to the time loop to the next step + timeLoop->advanceTimeStep(); + + // write vtk output + vtkWriter.write(timeLoop->time()); + + // calculate and print mass fluxes over the planes + flux.calculateMassOrMoleFluxes(); + if(GET_PROP_VALUE(TypeTag, EnableEnergyBalance)) + { + std::cout << "mass / energy flux at middle is: " << flux.netFlux("middle") << std::endl; + std::cout << "mass / energy flux at outlet is: " << flux.netFlux("outlet") << std::endl; + } + else + { + std::cout << "mass flux at middle is: " << flux.netFlux("middle") << std::endl; + std::cout << "mass flux at outlet is: " << flux.netFlux("outlet") << std::endl; + } + + // calculate and print volume fluxes over the planes + flux.calculateVolumeFluxes(); + std::cout << "volume flux at middle is: " << flux.netFlux("middle")[0] << std::endl; + std::cout << "volume flux at outlet is: " << flux.netFlux("outlet")[0] << std::endl; + + // report statistics of this time step + timeLoop->reportTimeStep(); + + // set new dt as suggested by newton controller + timeLoop->setTimeStepSize(newtonController->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/rans/test_pipe_laufer.input b/test/freeflow/rans/test_pipe_laufer.input new file mode 100644 index 0000000000..e479e08600 --- /dev/null +++ b/test/freeflow/rans/test_pipe_laufer.input @@ -0,0 +1,26 @@ +[TimeLoop] +DtInitial = 1 # [s] +TEnd = 100 # [s] + +[Grid] +Verbosity = true +Positions0 = 0.0 10.0 +Positions1 = 0.0 0.12345 0.2469 +Cells0 = 16 +Cells1 = 8 8 +Grading1 = 1.4 -1.4 + +Cells = 0 # TODO Remove this requirement + +[Problem] +Name = test_pipe_laufer # name passed to the output routines +InletVelocity = 2.5 # [m/s] +EnableGravity = false + +[Newton] +MaxSteps = 10 +MaxRelativeShift = 1e-5 + +[Vtk] +AddVelocity = true +WriteFaceData = false -- GitLab