Commit 3159cda5 authored by Timo Koch's avatar Timo Koch
Browse files

Merge branch 'port/ff-channel-pipe-with-staggered-grid' into 'master'

[ff/channel/pipe] Port ff-channel-pipe case to staggered grid

See merge request !2841
parents 9aaefc8f f7d910b3
Pipeline #9403 failed with stages
in 0 seconds
......@@ -18,7 +18,8 @@ print("Removed old log file ({})!".format(testname + '.log'))
# do the runs with different refinement
for i in [0, 1, 2]:
subprocess.call(['./' + testname] + testargs + ['-Grid.Refinement', str(i),
'-Problem.Name', testname])
'-Problem.Name', testname,
'-Problem.PrintConvergenceTestFile', 'true'])
# check the rates and append them to the log file
logfile = open(testname + '.log', "r+")
......
......@@ -26,15 +26,21 @@
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/partial.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/staggeredvtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager_yasp.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/assembly/staggeredfvassembler.hh>
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/multidomain/fvassembler.hh>
#include <dumux/multidomain/traits.hh>
#include <dumux/multidomain/staggeredfreeflow/couplingmanager.hh>
#include <dumux/multidomain/newtonsolver.hh>
#include <dumux/freeflow/navierstokes/velocityoutput.hh>
#include <test/freeflow/navierstokes/analyticalsolutionvectors.hh>
#include <test/freeflow/navierstokes/errors.hh>
#include "properties.hh"
......@@ -43,76 +49,102 @@ int main(int argc, char** argv)
{
using namespace Dumux;
// define the type tag for this problem
using MomentumTypeTag = Properties::TTag::PipeFlowMomentum;
using MassTypeTag = Properties::TTag::PipeFlowMass;
// initialize MPI, finalize is done automatically on exit
const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
// parse command line arguments and input file
Parameters::init(argc, argv);
// Define the sub problem type tags
using TypeTag = Properties::TTag::PipeFlow;
// try to create a grid (from the given grid file or the input file)
Dumux::GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
Dumux::GridManager<GetPropType<MomentumTypeTag, Properties::Grid>> gridManager;
gridManager.init();
// we compute on the leaf grid view
const auto& gridView = gridManager.grid().leafGridView();
// create the finite volume grid geometry
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
auto gridGeometry = std::make_shared<GridGeometry>(gridView);
using MomentumGridGeometry = GetPropType<MomentumTypeTag, Properties::GridGeometry>;
auto momentumGridGeometry = std::make_shared<MomentumGridGeometry>(gridView);
using MassGridGeometry = GetPropType<MassTypeTag, Properties::GridGeometry>;
auto massGridGeometry = std::make_shared<MassGridGeometry>(gridView);
// the problem (initial and boundary conditions)
using Problem = GetPropType<TypeTag, Properties::Problem>;
auto problem = std::make_shared<Problem>(gridGeometry);
// the coupling manager
using Traits = MultiDomainTraits<MomentumTypeTag, MassTypeTag>;
using CouplingManager = StaggeredFreeFlowCouplingManager<Traits>;
auto couplingManager = std::make_shared<CouplingManager>();
// the problem (boundary conditions)
using MomentumProblem = GetPropType<MomentumTypeTag, Properties::Problem>;
auto momentumProblem = std::make_shared<MomentumProblem>(momentumGridGeometry, couplingManager);
using MassProblem = GetPropType<MassTypeTag, Properties::Problem>;
auto massProblem = std::make_shared<MassProblem>(massGridGeometry, couplingManager);
// the solution vector
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
SolutionVector sol;
sol[GridGeometry::cellCenterIdx()].resize(gridGeometry->numCellCenterDofs());
sol[GridGeometry::faceIdx()].resize(gridGeometry->numFaceDofs());
constexpr auto momentumIdx = CouplingManager::freeFlowMomentumIndex;
constexpr auto massIdx = CouplingManager::freeFlowMassIndex;
using SolutionVector = typename Traits::SolutionVector;
SolutionVector x;
momentumProblem->applyInitialSolution(x[momentumIdx]);
massProblem->applyInitialSolution(x[massIdx]);
// the grid variables
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
gridVariables->init(sol);
using MomentumGridVariables = GetPropType<MomentumTypeTag, Properties::GridVariables>;
auto momentumGridVariables = std::make_shared<MomentumGridVariables>(momentumProblem, momentumGridGeometry);
using MassGridVariables = GetPropType<MassTypeTag, Properties::GridVariables>;
auto massGridVariables = std::make_shared<MassGridVariables>(massProblem, massGridGeometry);
couplingManager->init(momentumProblem, massProblem, std::make_tuple(momentumGridVariables, massGridVariables), x);
massGridVariables->init(x[massIdx]);
momentumGridVariables->init(x[momentumIdx]);
using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(std::make_tuple(momentumProblem, massProblem),
std::make_tuple(momentumGridGeometry, massGridGeometry),
std::make_tuple(momentumGridVariables, massGridVariables),
couplingManager);
// intialize the vtk output module
using IOFields = GetPropType<TypeTag, Properties::IOFields>;
StaggeredVtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, sol, problem->name());
using IOFields = GetPropType<MassTypeTag, Properties::IOFields>;
VtkOutputModule vtkWriter(*massGridVariables, x[massIdx], massProblem->name());
IOFields::initOutputModule(vtkWriter); // Add model specific output fields
vtkWriter.addVelocityOutput(std::make_shared<NavierStokesVelocityOutput<MassGridVariables>>());
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, gridGeometry, gridVariables);
// 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);
using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
nonLinearSolver.solve(sol);
nonLinearSolver.solve(x);
vtkWriter.write(1.0);
// print discrete L2 and Linfity errors
// the discrete L2 and Linfity errors
const bool printErrors = getParam<bool>("Problem.PrintErrors", false);
const bool printConvergenceTestFile = getParam<bool>("Problem.PrintConvergenceTestFile", false);
if (printErrors || printConvergenceTestFile)
{
NavierStokesErrors errors(problem, sol);
NavierStokesErrorCSVWriter(
problem, std::to_string(sol[GridGeometry::cellCenterIdx()].size())
).printErrors(errors);
if (printConvergenceTestFile)
convergenceTestAppendErrors(problem, errors);
}
NavierStokesTest::Errors errors(momentumProblem, massProblem, x);
NavierStokesTest::ErrorCSVWriter errorCSVWriter(momentumProblem, massProblem);
if (printErrors || printConvergenceTestFile)
{
NavierStokesTest::Errors errors(momentumProblem, massProblem, x);
NavierStokesTest::ErrorCSVWriter errorCSVWriter(momentumProblem, massProblem);
errorCSVWriter.printErrors(errors);
if (printConvergenceTestFile)
convergenceTestAppendErrors(momentumProblem, massProblem, errors);
}
// print dumux end message
if (mpiHelper.rank() == 0)
......
......@@ -20,15 +20,17 @@
#ifndef DUMUX_TEST_FREEFLOW_PIPE_PROBLEM_HH
#define DUMUX_TEST_FREEFLOW_PIPE_PROBLEM_HH
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/numeqvector.hh>
#include <dumux/freeflow/navierstokes/problem.hh>
#include <dumux/freeflow/navierstokes/boundarytypes.hh>
#include <dumux/freeflow/navierstokes/momentum/fluxhelper.hh>
#include <dumux/freeflow/navierstokes/scalarfluxhelper.hh>
namespace Dumux {
/*!
* \ingroup NavierStokesTests
* \brief Freeflow problem for pipe flow
* Simulation of a radially-symmetric pipe flow with circular cross-section
*/
......@@ -37,19 +39,24 @@ class FreeFlowPipeProblem : public NavierStokesProblem<TypeTag>
{
using ParentType = NavierStokesProblem<TypeTag>;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using GridView = typename GridGeometry::GridView;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using BoundaryTypes = NavierStokesBoundaryTypes<PrimaryVariables::size()>;
using PrimaryVariables = typename ParentType::PrimaryVariables;
using NumEqVector = typename ParentType::NumEqVector;
using BoundaryTypes = typename ParentType::BoundaryTypes;
using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
public:
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
using Indices = typename ModelTraits::Indices;
FreeFlowPipeProblem(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry)
FreeFlowPipeProblem(std::shared_ptr<const GridGeometry> gridGeometry, std::shared_ptr<CouplingManager> couplingManager)
: ParentType(gridGeometry, couplingManager)
{
name_ = getParamFromGroup<std::string>(this->paramGroup(), "Problem.Name");
......@@ -76,62 +83,106 @@ public:
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary segment.
*/
BoundaryTypes boundaryTypes(const Element& element,
const SubControlVolumeFace& scvf) const
BoundaryTypes boundaryTypesAtPos(const GlobalPosition& globalPos) const
{
BoundaryTypes values;
const auto& globalPos = scvf.dofPosition();
// inlet
if (onLowerBoundary_(globalPos))
if constexpr (ParentType::isMomentumProblem())
{
values.setDirichlet(Indices::velocityXIdx);
values.setDirichlet(Indices::velocityYIdx);
if (onLowerBoundary_(globalPos) || onOuterBoundary_(globalPos))
values.setAllDirichlet();
else if (onInnerBoundary_(globalPos))
{
values.setDirichlet(Indices::velocityXIdx);
values.setNeumann(Indices::momentumYBalanceIdx);
}
else
values.setAllNeumann();
}
// outlet
else if (onUpperBoundary_(globalPos))
{
values.setDirichlet(Indices::pressureIdx);
}
// pipe centerline
else if (onInnerBoundary_(globalPos))
else
values.setAllNeumann();
return values;
}
/*!
* \brief Evaluates the boundary conditions for a Neumann control volume.
*
* \param element The element for which the Neumann boundary condition is set
* \param fvGeometry The fvGeometry
* \param elemVolVars The element volume variables
* \param elemFaceVars The element face variables
* \param scvf The boundary sub control volume face
*/
template<class ElementVolumeVariables, class ElementFluxVariablesCache>
NumEqVector neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFluxVariablesCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
NumEqVector values(0.0);
const auto& globalPos = scvf.ipGlobal();
if constexpr (ParentType::isMomentumProblem())
{
values.setAllSymmetry();
if (onInnerBoundary_(globalPos))
return values; // zero shear stress at symmetry axis
if (onUpperBoundary_(globalPos))
values = NavierStokesMomentumBoundaryFluxHelper::fixedPressureMomentumFlux(*this, fvGeometry, scvf, elemVolVars, elemFluxVarsCache, analyticalPressure(globalPos), true /*zeroNormalVelocityGradient*/);
}
// pipe wall
else if (onOuterBoundary_(globalPos))
else
{
values.setDirichlet(Indices::velocityXIdx);
values.setDirichlet(Indices::velocityYIdx);
using FluxHelper = NavierStokesScalarBoundaryFluxHelper<AdvectiveFlux<ModelTraits>>;
if (onLowerBoundary_(globalPos) || onUpperBoundary_(globalPos))
values = FluxHelper::scalarOutflowFlux(*this, element, fvGeometry, scvf, elemVolVars);
}
return values;
}
PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
{ return analyticalSolution(globalPos); }
/*!
* \brief Evaluates the initial value for a control volume.
*
* \param globalPos The global position
*/
PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
{
PrimaryVariables values(0.0);
values[Indices::pressureIdx] = initialPressure_;
if constexpr (ParentType::isMomentumProblem())
return values;
else
values[Indices::pressureIdx] = initialPressure_;
return values;
}
PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
{ return analyticalSolution(globalPos); }
PrimaryVariables analyticalSolution(const GlobalPosition& globalPos, Scalar time = 0.0) const
{
PrimaryVariables values(0.0);
// paraboloid velocity profile
const auto r = globalPos[0] - this->gridGeometry().bBoxMin()[0];
const auto y = globalPos[1] - this->gridGeometry().bBoxMin()[1];
values[Indices::velocityXIdx] = 0.0;
values[Indices::velocityYIdx] = 2.0*meanInletVelocity_*(1.0 - r*r/(pipeRadius_*pipeRadius_));
values[Indices::pressureIdx] = (pipeLength_-y)*meanInletVelocity_*8.0*mu_/(pipeRadius_*pipeRadius_);
if constexpr (ParentType::isMomentumProblem())
{
const auto r = globalPos[0] - this->gridGeometry().bBoxMin()[0];
values[Indices::velocityXIdx] = 0.0;
values[Indices::velocityYIdx] = 2.0*meanInletVelocity_*(1.0 - r*r/(pipeRadius_*pipeRadius_));
}
else
values[Indices::pressureIdx] = analyticalPressure(globalPos);
return values;
}
Scalar analyticalPressure(const GlobalPosition& globalPos) const
{
const auto y = globalPos[1];
return (pipeLength_-y)*meanInletVelocity_*8.0*mu_/(pipeRadius_*pipeRadius_);
}
private:
bool onInnerBoundary_(const GlobalPosition &globalPos) const
{ return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_; }
......
......@@ -27,19 +27,26 @@
#include <dune/grid/yaspgrid.hh>
#include <dumux/discretization/staggered/freeflow/properties.hh>
#include <dumux/discretization/extrusion.hh>
#include <dumux/freeflow/navierstokes/model.hh>
#include <dumux/freeflow/navierstokes/momentum/model.hh>
#include <dumux/freeflow/navierstokes/mass/1p/model.hh>
#include <dumux/discretization/fcstaggered.hh>
#include <dumux/discretization/cctpfa.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
#include <dumux/material/components/constant.hh>
#include <dumux/multidomain/staggeredfreeflow/couplingmanager.hh>
#include "problem.hh"
namespace Dumux::Properties {
// Create new type tags
namespace TTag {
struct PipeFlow { using InheritsFrom = std::tuple<NavierStokes, StaggeredFreeFlowModel>; };
struct PipeFlow {};
struct PipeFlowMomentum { using InheritsFrom = std::tuple<PipeFlow, NavierStokesMomentum, FaceCenteredStaggeredModel>; };
struct PipeFlowMass { using InheritsFrom = std::tuple<PipeFlow, NavierStokesMassOneP, CCTpfaModel>; };
} // end namespace TTag
// the fluid system
......@@ -69,16 +76,37 @@ struct EnableGridVolumeVariablesCache<TypeTag, TTag::PipeFlow> { static constexp
// rotation-symmetric grid geometry forming a cylinder channel
template<class TypeTag>
struct GridGeometry<TypeTag, TTag::PipeFlow>
struct GridGeometry<TypeTag, TTag::PipeFlowMomentum>
{
static constexpr auto upwindSchemeOrder = getPropValue<TypeTag, Properties::UpwindSchemeOrder>();
static constexpr bool enableCache = getPropValue<TypeTag, Properties::EnableGridGeometryCache>();
using GridView = typename GetPropType<TypeTag, Properties::Grid>::LeafGridView;
struct GGTraits : public StaggeredFreeFlowDefaultFVGridGeometryTraits<GridView, upwindSchemeOrder>
struct GGTraits : public FaceCenteredStaggeredDefaultGridGeometryTraits<GridView>
{ using Extrusion = RotationalExtrusion<0>; };
using type = StaggeredFVGridGeometry<GridView, enableCache, GGTraits>;
using type = FaceCenteredStaggeredFVGridGeometry<GridView, enableCache, GGTraits>;
};
// rotation-symmetric grid geometry forming a cylinder channel
template<class TypeTag>
struct GridGeometry<TypeTag, TTag::PipeFlowMass>
{
static constexpr bool enableCache = getPropValue<TypeTag, Properties::EnableGridGeometryCache>();
using GridView = typename GetPropType<TypeTag, Properties::Grid>::LeafGridView;
struct GGTraits : public CCTpfaDefaultGridGeometryTraits<GridView>
{ using Extrusion = RotationalExtrusion<0>; };
using type = CCTpfaFVGridGeometry<GridView, enableCache, GGTraits>;
};
template<class TypeTag>
struct CouplingManager<TypeTag, TTag::PipeFlow>
{
private:
using Traits = MultiDomainTraits<TTag::PipeFlowMomentum, TTag::PipeFlowMass>;
public:
using type = StaggeredFreeFlowCouplingManager<Traits>;
};
} // end namespace Dumux::Properties
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment