### Merge branch 'feature/freeflowchannel_example' into 'master'

Feature/freeflowchannel example

See merge request !1682
 ... ... @@ -33,3 +33,10 @@ You learn how to * solve a shallow water flow problem including bottom friction * computate and output (VTK) an analytical reference solution ### [:open_file_folder: Example 4: Freeflow channel](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/tree/master/examples/freeflowchannel) In this example, we simulate a free flow between two plates in two dimensions. You learn how to * solve a free flow problem * set outflow boundary conditions in the free-flow context
 doc/intro.md problem.hh main.cc doc/results.md
 dune_symlink_to_source_files(FILES "params.input") dumux_add_test(NAME example_freeflow_channel_navierstokes LABELS freeflow navierstokes SOURCES main.cc CMAKE_GUARD HAVE_UMFPACK COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py CMD_ARGS --script fuzzy --files${CMAKE_SOURCE_DIR}/test/references/example_ff_channel-reference.vtu ${CMAKE_CURRENT_BINARY_DIR}/example_ff_channel-00001.vtu --command "${CMAKE_CURRENT_BINARY_DIR}/example_freeflow_channel_navierstokes params.input -Problem.Name example_ff_channel")
This diff is collapsed.
 This example is based on dumux/test/freeflow/navierstokes/channel/2d. # Freeflow through a channel ## Problem set-up This example contains a stationary free flow of a fluid through two parallel solid plates in two dimensions from left to right. The figure below shows the simulation set-up. The fluid flows into the system at the left with a constant velocity of $ v = 1~\frac{\text{m}}{\text{s}} $. The inflow velocity profile is a block profile. Due to the no-slip, no-flow boundary conditions at the top and bottom plate, the velocity profile gradually assumes a parabolic shape along the channel. At the outlet, the pressure is fixed and a zero velocity gradient in x-direction is assumed. The physical domain, which is modeled is the rectangular domain $x\in[0,10],~y\in[0,1]$. ![](./img/setup.png) ## Model description The Stokes model without gravitation and without sources or sinks for a stationary, incompressible, laminar, single phase, one-component, isothermal ($T=10^\circ C$) flow is considered assuming a Newtonian fluid of constant density $ \varrho = 1~\frac{\text{kg}}{\text{m}^3} $ and constant kinematic viscosity $ \nu = 1~\frac{\text{m}^2}{\text{s}} $. The momentum balance math - \nabla\cdot\left(\mu\left(\nabla\boldsymbol{u}+\nabla\boldsymbol{u}^{\text{T}}\right)\right)+ \nabla p = 0  with density $\varrho$, velocity $\boldsymbol{u}$, dynamic viscosity $\mu=\varrho\nu$ and pressure $p$ and the mass balance math \nabla \cdot \left(\varrho\boldsymbol{u}\right) =0  are discretized using a staggered-grid finite-volume scheme as spatial discretization with pressures and velocity components as primary variables. For details on the discretization scheme, have a look at the dumux [handbook](https://dumux.org/handbook). In the following, we take a close look at the files containing the set-up: At first, boundary conditions are set in problem.hh for the Navier-Stokes model. Afterwards, we show the different steps for solving the model in the source file main.cc. At the end, we show some simulation results.
 ## Results This example computes the following stationary velocity profile: ![](./img/velocity.png) and stationary pressure profile: ![](./img/pressure.png)

11.6 KB

40.3 KB

This source diff could not be displayed because it is too large. You can view the blob instead.

16.8 KB

 // -*- 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 . * *****************************************************************************/ // We look now at the main file for the channel problem. // ### Includes // Necessary files are included. //
// Toggle to expand details // // First, the Dune configuration file is include, the standard header file for C++ to get time and date information // and another standard header for in- and output. // //
// Toggle to expand code (includes of problem file and of standard headers) // #include #include #include //
// // Dumux is based on DUNE, the Distributed and Unified Numerics Environment, which provides several grid managers and // linear solvers. So we need some includes from that. //
// Toggle to expand code (dune includes) // #include #include #include #include #include //
// // In Dumux, a property system is used to specify the model. For this, different properties are defined containing // type definitions, values and methods. All properties are declared in the file properties.hh. // The file parameters.hh contains the parameter class, which manages the definition of input parameters by a default // value, the inputfile or the command line. // The file dumuxmessage.hh contains the class defining the start and end message of the simulation. // The file valgrind.hh contains debugging funcionality. // // The file seqsolverbackend.hh contains the class, which defines the sequential linear solver backends. // The nonlinear Newton's method is included, as well as the assembler, which assembles the linear systems for staggered-grid finite volume schemes. // The containing class in the file diffmethod.hh defines the different differentiation methods used to compute the derivatives of the residual. // // We need the class in staggeredvtkoutputmodule.hh to simplify the writing of dumux simulation data to VTK format. // The gridmanager constructs a grid from the information in the input or grid file. There is a specification for the // different supported grid managers. // // The following class contains functionality for additional flux output to the console. //
// Toggle to expand code (dumux includes) // #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "problem.hh" //
//
// // // ### Setup basic properties for our simulation // We setup the DuMux properties for our simulation (click [here](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux-course/blob/master/slides/dumux-course-properties.pdf) for DuMux course slides on the property system) within the namespace Properties, which is a sub-namespace of Dumux. // 1. For every test problem, a new TypeTag has to be created, which is done within the namespace TTag (subnamespace of Properties). It inherits from the Navier-Stokes flow model and the staggered-grid discretization scheme. // 2. The grid is chosen to be a two-dimensional YASP grid. // 3. We set the FluidSystem to be a one-phase liquid with a single component. The class Component::Constant refers to a component with constant fluid properties (density, viscosity, ...) that can be set via the input file in the group [0.Component] where the number is the identifier given as template argument to the class template Component::Constant. // 4. The problem class ChannelExampleProblem, which is forward declared before we enter namespace Dumux and defined later in this file, is defined to be the problem used in this test problem (charaterized by the TypeTag ChannelExample). The fluid system, which contains information about the properties such as density, viscosity or diffusion coefficient of the fluid we're simulating, is set to a constant one phase liquid. // 5. We enable caching for the following classes (which stores values that were already calculated for later usage and thus results in higher memory usage but improved CPU speed): the grid volume variables, the grid flux variables, the finite volume grid geometry. // namespace Dumux::Properties { namespace TTag { struct ChannelExample { using InheritsFrom = std::tuple; }; } // namespace TTag template struct Grid { using type = Dune::YaspGrid<2>; }; template struct Problem { using type = Dumux::ChannelExampleProblem ; }; template struct FluidSystem { using Scalar = GetPropType; using type = FluidSystems::OnePLiquid >; }; template struct EnableGridVolumeVariablesCache { static constexpr bool value = true; }; template struct EnableGridFluxVariablesCache { static constexpr bool value = true; }; template struct EnableGridGeometryCache { static constexpr bool value = true; }; } // end namespace Dumux::Properties // // // ### Beginning of the main function // We begin the main function by making the type tag ChannelExample, that we defined in problem.hh for this test problem available here. // Then we initializing the message passing interface (MPI), even if we do not plan to run the application in parallel. Finalizing of the MPI is done automatically on exit. // We continue by printing the dumux start message and parsing the command line arguments and runtimeparameters from the input file in the init function. // int main(int argc, char** argv) try { using namespace Dumux; using TypeTag = Properties::TTag::ChannelExample; const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv); if (mpiHelper.rank() == 0) DumuxMessage::print(/*firstCall=*/true); Parameters::init(argc, argv); // // ### Set-up and solving of the problem // // A gridmanager tries to create the grid either from a grid file or the input file. You can learn more about grids in // Dumux in the [grid exercise](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux-course/-/tree/master/exercises/exercise-grids). // Then, we compute the leaf grid view, which is the overlap of all grid levels in hierarchical grids. // We create (std::make_shared) and initialize (update) the finite volume geometry to build up the subcontrolvolumes // and subcontrolvolume faces for each element of the grid partition. // In the problem, we define the boundary and initial conditions. // We set a solution vector which has a part (indexed by cellCenterIdx) for degrees of freedom (dofs) living // in grid cell centers - pressures - and a part (indexed by faceIdx) for degrees of freedom livin on grid cell faces. // We initialize the solution vector by what was defined as the initial solution in problem.hh (applyInitialSolution) // and then use the solution vector to intialize the gridVariables. Grid variables in general contain the s // primary variables (velocities, pressures) as well as secondary variables (density, viscosity, ...). // We then initialize the vtkoutput. Each model has a predefined model-specific output with relevant parameters // for that model. Here, it is pressure, velocity, density and process rank (relevant in the case of parallelisation). // GridManager> gridManager; gridManager.init(); const auto& leafGridView = gridManager.grid().leafGridView(); using GridGeometry = GetPropType; auto gridGeometry = std::make_shared(leafGridView); gridGeometry->update(); using Problem = GetPropType; auto problem = std::make_shared(gridGeometry); using SolutionVector = GetPropType; SolutionVector x; x[GridGeometry::cellCenterIdx()].resize(gridGeometry->numCellCenterDofs()); x[GridGeometry::faceIdx()].resize(gridGeometry->numFaceDofs()); problem->applyInitialSolution(x); using GridVariables = GetPropType; auto gridVariables = std::make_shared(problem, gridGeometry); gridVariables->init(x); using IOFields = GetPropType; StaggeredVtkOutputModule vtkWriter(*gridVariables, x, problem->name()); IOFields::initOutputModule(vtkWriter); // Add model specific output fields vtkWriter.write(0.0); // // We set up two surfaces over which fluxes are calculated. // We determine the extensions [xMin,xMax]x[yMin,yMax] of the physical domain. // The first surface (added by the first call of addSurface) 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 position of the surface (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 surface. This assumes a regular cartesian grid. // The second surface (second call of addSurface) is placed at the outlet of the channel. // FluxOverSurface, GetPropType> flux(*gridVariables, x); using Scalar = GetPropType; const Scalar xMin = gridGeometry->bBoxMin(); const Scalar xMax = gridGeometry->bBoxMax(); const Scalar yMin = gridGeometry->bBoxMin(); const Scalar yMax = gridGeometry->bBoxMax(); const Scalar planePosMiddleX = xMin + 0.5*(xMax - xMin); int numCellsX = getParam>("Grid.Cells"); const unsigned int refinement = getParam("Grid.Refinement", 0); numCellsX *= (1<; using Element = typename GridView::template Codim<0>::Entity; using GlobalPosition = typename Element::Geometry::GlobalCoordinate; const auto p0middle = GlobalPosition{planePosMiddleX + offsetX, yMin}; const auto p1middle = GlobalPosition{planePosMiddleX + offsetX, yMax}; flux.addSurface("middle", p0middle, p1middle); const auto p0outlet = GlobalPosition{xMax, yMin}; const auto p1outlet = GlobalPosition{xMax, yMax}; flux.addSurface("outlet", p0outlet, p1outlet); // // The incompressible Stokes equation depends only linearly on the velocity, so we could use a linear solver to solve the problem. // Here, we use the show the more general case which would also work for incompressible fluids or the // Navier-Stokes equation. We use non-linear Newton solver for the solution. // In each step of the Newton solver, we assemble the respective linear system, including the jacobian matrix and the residual by the // assembler. The linear systems are here solved by the direct linear solver UMFPack. // As a postprocessing, we calculate mass and volume fluxes over the planes specified above. //
// Toggle to expand code (assembly, solution process, postprocessing) // using Assembler = StaggeredFVAssembler; auto assembler = std::make_shared(problem, gridGeometry, gridVariables); using LinearSolver = Dumux::UMFPackBackend; auto linearSolver = std::make_shared(); using NewtonSolver = Dumux::NewtonSolver; NewtonSolver nonLinearSolver(assembler, linearSolver); nonLinearSolver.solve(x); flux.calculateMassOrMoleFluxes(); flux.calculateVolumeFluxes(); //
// // ### Final Output // We write the vtk output and print the mass/energy/volume fluxes over the planes. // We conclude by printing the dumux end message. After the end of the main function, // possibly catched error messages are printed. //
// Toggle to expand code (final output) // vtkWriter.write(1.0); if(GetPropType::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; } std::cout << "volume flux at middle is: " << flux.netFlux("middle") << std::endl; std::cout << "volume flux at outlet is: " << flux.netFlux("outlet") << std::endl; 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; } //
//
 [Grid] UpperRight = 10 1 # [m] Cells = 100 50 # [m] [Problem] Name = example_ff_channel # name passed to the output routines InletVelocity = 1 # [m/s] EnableGravity = false EnableInertiaTerms = false # Stokes, not Navier Stokes, in this example [Component] LiquidDensity = 1.0 # [kg/m^3] LiquidKinematicViscosity = 1.0 # [m^2/s] [ Newton ] MaxSteps = 10 MaxRelativeShift = 1e-8 [Vtk] WriteFaceData = false
 // -*- 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 . * *****************************************************************************/ // ### Header guard #ifndef DUMUX_CHANNEL_TEST_PROBLEM_HH #define DUMUX_CHANNEL_TEST_PROBLEM_HH // // ### The problem class // We enter the problem class where all necessary initial and boundary conditions are set for our simulation. // // As this is a Stokes problem, we inherit from the basic NavierStokesProblem. //
Toggle to expand code: #include namespace Dumux { template class ChannelExampleProblem : public NavierStokesProblem { //
// // We use convenient declarations that we derive from the property system. //
// Toggle to expand code (convenient declarations) // using ParentType = NavierStokesProblem; using BoundaryTypes = GetPropType; using GridGeometry = GetPropType; using FVElementGeometry = typename GridGeometry::LocalView; using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; using Indices = typename GetPropType::Indices; using NumEqVector = GetPropType; using PrimaryVariables = GetPropType; using Scalar = GetPropType; using Element = typename GridGeometry::GridView::template Codim<0>::Entity; using GlobalPosition = typename Element::Geometry::GlobalCoordinate; public: //
// // There follows the constructor of our problem class: // Within the constructor, we set the inlet velocity to a run-time specified value. // As no run-time value is specified, we set the outlet pressure to 1.1e5 Pa. //
// Toggle to expand code (constructor) // ChannelExampleProblem(std::shared_ptr gridGeometry) : ParentType(gridGeometry) { inletVelocity_ = getParam("Problem.InletVelocity"); outletPressure_ = getParam("Problem.OutletPressure", 1.1e5); } //
// // Now, we define the type of initial and boundary conditions depending on location. // Two types of boundary conditions can be specified: Dirichlet and Neumann. On a Dirichlet boundary, // the values of the primary variables need to be fixed. // On a Neumann boundary condition, values for derivatives need to be fixed. // When Dirichlet conditions are set for the pressure, the derivative of the velocity // vector with respect to the direction normal to the boundary is automatically set to // zero. This boundary condition is called in-/outflow boundary condition in Dumux. // In the following we specify Dirichlet boundaries for velocity on the left of our domain // if isInlet_ is true, Dirichlet boundaries for pressure on the right of our domain // if isOutlet_ is true and specify Dirichlet boundaries for velocity on the top and bottom // of our domain else. //
// Toggle to expand code (boundaryTypesAtPos) // BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const { BoundaryTypes values; if(isInlet_(globalPos)) { values.setDirichlet(Indices::velocityXIdx); values.setDirichlet(Indices::velocityYIdx); } else if(isOutlet_(globalPos)) { values.setDirichlet(Indices::pressureIdx); } else { values.setDirichlet(Indices::velocityXIdx); values.setDirichlet(Indices::velocityYIdx); } return values; } //
// // Second, we specify the values for the Dirichlet boundaries. We need to fix the values of our primary variables. // To ensure a no-slip boundary condition at the top and bottom of the channel, the Dirichlet velocity // in x-direction is set to zero if not at the inlet. //
// Toggle to expand code (dirichletAtPos) // PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const { PrimaryVariables values = initialAtPos(globalPos); if(!isInlet_(globalPos)) { values[Indices::velocityXIdx] = 0.0; } return values; } //
// // We specify the values for the initial conditions. // We assign constant values for pressure and velocity components. //
// Toggle to expand code (initialAtPos) // PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const { PrimaryVariables values; values[Indices::pressureIdx] = outletPressure_; values[Indices::velocityXIdx] = inletVelocity_; values[Indices::velocityYIdx] = 0.0; return values; } //
// // We need to specify a constant temperature for our isothermal problem. // We set it to 10°C. //
// Toggle to expand code (temperature) // Scalar temperature() const { return 273.15 + 10; } private: //
// // The inlet is at the left side of the physical domain. //
// Toggle to expand code (isInlet_) // bool isInlet_(const GlobalPosition& globalPos) const { return globalPos < eps_; } //
// // The outlet is at the right side of the physical domain. //
// Toggle to expand code (isOutlet_) // bool isOutlet_(const GlobalPosition& globalPos) const { return globalPos > this->gridGeometry().bBoxMax() - eps_; } //
// // Finally, private variables are declared: //
// Toggle to expand code (private variables) // static constexpr Scalar eps_=1e-6; Scalar inletVelocity_; Scalar outletPressure_; }; } #endif //
//
This diff is collapsed.