Commit 3d1bc8d4 authored by Timo Koch's avatar Timo Koch
Browse files

[test][ff] Port channel 2d tests to new staggered implementation

parent c5b3c40c
......@@ -34,7 +34,7 @@ dumux_add_test(NAME test_ff_stokes_channel_neumann_x_neumann_y
${CMAKE_CURRENT_BINARY_DIR}/test_ff_stokes_channel_nxny-00001.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_ff_channel params.input -Grid.UpperRight \"2 1\" -Grid.Cells \"50 25\"
-Problem.Name test_ff_stokes_channel_nxny -Problem.OutletCondition NeumannX_NeumannY
-Problem.UseVelocityProfile true -Problem.OutletPressure 0 -Problem.IsStationary 1")
-Problem.UseVelocityProfile true -Problem.OutletPressure 0 -Problem.IsStationary true")
dumux_add_test(NAME test_ff_stokes_channel_do_nothing
TARGET test_ff_channel
......
......@@ -24,24 +24,28 @@
#include <config.h>
#ifndef ISOTHERMAL
#define ISOTHERMAL 0
#endif
#include <ctime>
#include <iostream>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh>
#include <dune/grid/io/file/vtk.hh>
#include <dune/istl/io.hh>
#include <dumux/assembly/staggeredfvassembler.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/properties.hh>
#include <dumux/freeflow/navierstokes/staggered/fluxoversurface.hh>
#include <dumux/io/grid/gridmanager.hh>
#include <dumux/io/staggeredvtkoutputmodule.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/multidomain/newtonsolver.hh>
#include <dumux/multidomain/fvassembler.hh>
#include <dumux/multidomain/traits.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/freeflow/navierstokes/velocityoutput.hh>
#include <test/freeflow/navierstokes/analyticalsolutionvectors.hh>
#include <test/freeflow/navierstokes/errors.hh>
......@@ -53,7 +57,8 @@ int main(int argc, char** argv)
using namespace Dumux;
// define the type tag for this problem
using TypeTag = Properties::TTag::ChannelTest;
using MomentumTypeTag = Properties::TTag::ChannelTestMomentum;
using MassTypeTag = Properties::TTag::ChannelTestMass;
// initialize MPI, finalize is done automatically on exit
const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
......@@ -66,7 +71,7 @@ int main(int argc, char** argv)
Parameters::init(argc, argv);
// try to create a grid (from the given grid file or the input file)
GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
GridManager<GetPropType<MassTypeTag, Properties::Grid>> gridManager;
gridManager.init();
////////////////////////////////////////////////////////////
......@@ -77,15 +82,26 @@ int main(int argc, char** argv)
const auto& leafGridView = gridManager.grid().leafGridView();
// create the finite volume grid geometry
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
using MomentumGridGeometry = GetPropType<MomentumTypeTag, Properties::GridGeometry>;
auto momentumGridGeometry = std::make_shared<MomentumGridGeometry>(leafGridView);
using MassGridGeometry = GetPropType<MassTypeTag, Properties::GridGeometry>;
auto massGridGeometry = std::make_shared<MassGridGeometry>(leafGridView);
// the coupling manager
using Traits = MultiDomainTraits<MomentumTypeTag, MassTypeTag>;
using CouplingManager = StaggeredFreeFlowCouplingManager<Traits>;
auto couplingManager = std::make_shared<CouplingManager>();
// the problem (initial and boundary conditions)
using Problem = GetPropType<TypeTag, Properties::Problem>;
auto problem = std::make_shared<Problem>(gridGeometry);
// 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);
// get some time loop parameters
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Scalar = GetPropType<MassTypeTag, Properties::Scalar>;
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
auto dt = getParam<Scalar>("TimeLoop.DtInitial");
......@@ -94,100 +110,85 @@ int main(int argc, char** argv)
Scalar restartTime = getParam<Scalar>("Restart.Time", 0);
// the solution vector
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
constexpr auto momentumIdx = CouplingManager::freeFlowMomentumIndex;
constexpr auto massIdx = CouplingManager::freeFlowMassIndex;
using SolutionVector = typename Traits::SolutionVector;
SolutionVector x;
x[GridGeometry::cellCenterIdx()].resize(gridGeometry->numCellCenterDofs());
x[GridGeometry::faceIdx()].resize(gridGeometry->numFaceDofs());
problem->applyInitialSolution(x);
x[momentumIdx].resize(momentumGridGeometry->numDofs());
x[massIdx].resize(massGridGeometry->numDofs());
momentumProblem->applyInitialSolution(x[momentumIdx]);
massProblem->applyInitialSolution(x[massIdx]);
auto xOld = x;
// instantiate time loop
auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(restartTime, dt, tEnd);
timeLoop->setMaxTimeStepSize(maxDt);
problem->setTimeLoop(timeLoop);
massProblem->setTime(timeLoop->time());
momentumProblem->setTime(timeLoop->time());
// the grid variables
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
gridVariables->init(x);
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);
// initialize the coupling stencils
couplingManager->init(momentumProblem, massProblem, std::make_tuple(momentumGridVariables, massGridVariables), x, xOld);
// intializing the gridvariables requires the coupling manager to be set up
momentumGridVariables->init(x[momentumIdx]);
massGridVariables->init(x[massIdx]);
// initialize the vtk output module
using IOFields = GetPropType<TypeTag, Properties::IOFields>;
StaggeredVtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, 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>>());
const bool isStationary = getParam<bool>("Problem.IsStationary", false);
NavierStokesAnalyticalSolutionVectors analyticalSolVectors(problem);
if (problem->hasAnalyticalSolution())
using NSAnalyticSol = NavierStokesTest::AnalyticalSolutionVectors<MomentumProblem, MassProblem>;
std::unique_ptr<NSAnalyticSol> analyticalSolVectors;
if (momentumProblem->hasAnalyticalSolution())
{
vtkWriter.addField(analyticalSolVectors.getAnalyticalPressureSolution(), "pressureExact");
vtkWriter.addField(analyticalSolVectors.getAnalyticalVelocitySolution(), "velocityExact");
vtkWriter.addFaceField(analyticalSolVectors.getAnalyticalVelocitySolutionOnFace(), "faceVelocityExact");
analyticalSolVectors = std::make_unique<NSAnalyticSol>(momentumProblem, massProblem);
vtkWriter.addField(analyticalSolVectors->analyticalPressureSolution(), "pressureExact");
vtkWriter.addField(analyticalSolVectors->analyticalVelocitySolution(), "velocityExact");
// vtkWriter.addFaceField(analyticalSolVectors.analyticalVelocitySolutionOnFace(), "faceVelocityExact");
}
vtkWriter.write(restartTime);
// the assembler with time loop for instationary problem
using Assembler = StaggeredFVAssembler<TypeTag, DiffMethod::numeric>;
auto assembler = isStationary ? std::make_shared<Assembler>(problem, gridGeometry, gridVariables)
: std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, xOld);
const bool isStationary = getParam<bool>("Problem.IsStationary", false);
using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
auto assembler = isStationary ? std::make_shared<Assembler>(
std::make_tuple(momentumProblem, massProblem),
std::make_tuple(momentumGridGeometry, massGridGeometry),
std::make_tuple(momentumGridVariables, massGridVariables),
couplingManager
) : std::make_shared<Assembler>(
std::make_tuple(momentumProblem, massProblem),
std::make_tuple(momentumGridGeometry, massGridGeometry),
std::make_tuple(momentumGridVariables, massGridVariables),
couplingManager,
timeLoop, xOld
);
// 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);
// set up two surfaces over which fluxes are calculated
FluxOverSurface<GridVariables,
SolutionVector,
GetPropType<TypeTag, Properties::ModelTraits>,
GetPropType<TypeTag, Properties::LocalResidual>> flux(*gridVariables, x);
using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
const Scalar xMin = gridGeometry->bBoxMin()[0];
const Scalar xMax = gridGeometry->bBoxMax()[0];
const Scalar yMin = gridGeometry->bBoxMin()[1];
const Scalar yMax = gridGeometry->bBoxMax()[1];
// The first surface 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.
const Scalar planePosMiddleX = xMin + 0.5*(xMax - xMin);
int numCellsX = getParam<std::vector<int>>("Grid.Cells")[0];
const unsigned int refinement = getParam<unsigned int>("Grid.Refinement", 0);
numCellsX *= (1<<refinement);
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.addSurface("middle", p0middle, p1middle);
// The second surface is placed at the outlet of the channel.
const auto p0outlet = GlobalPosition{xMax, yMin};
const auto p1outlet = GlobalPosition{xMax, yMax};
flux.addSurface("outlet", p0outlet, p1outlet);
using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
if (isStationary)
{
Dune::Timer timer;
// solve the non-linear system with time step control
nonLinearSolver.solve(x);
// print discrete L2 and Linfity errors
if (problem->hasAnalyticalSolution())
if (momentumProblem->hasAnalyticalSolution())
{
// print discrete L2 and Linfity errors
const bool printErrors = getParam<bool>("Problem.PrintErrors", false);
......@@ -195,41 +196,24 @@ int main(int argc, char** argv)
if (printErrors || printConvergenceTestFile)
{
NavierStokesErrors errors(problem, x);
NavierStokesErrorCSVWriter(
problem, std::to_string(x[GridGeometry::cellCenterIdx()].size())
NavierStokesTest::Errors errors(momentumProblem, massProblem, x);
NavierStokesTest::ErrorCSVWriter(
momentumProblem, massProblem, std::to_string(x[massIdx].size())
).printErrors(errors);
if (printConvergenceTestFile)
convergenceTestAppendErrors(problem, errors);
NavierStokesTest::convergenceTestAppendErrors(momentumProblem, massProblem, errors);
}
}
// write vtk output
vtkWriter.write(1.0);
// calculate and print mass fluxes over the planes
flux.calculateMassOrMoleFluxes();
if(GetPropType<TypeTag, Properties::ModelTraits>::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;
// write vtk output
vtkWriter.write(1.0);
}
// 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;
timer.stop();
}
else
{
if (getParam<Scalar>("Problem.InletVelocity") > 1e-6)
timeLoop->setCheckPoint({200.0, 210.0});
// time loop
timeLoop->start(); do
{
......@@ -238,7 +222,8 @@ int main(int argc, char** argv)
// make the new solution the old solution
xOld = x;
gridVariables->advanceTimeStep();
momentumGridVariables->advanceTimeStep();
massGridVariables->advanceTimeStep();
// advance to the time loop to the next step
timeLoop->advanceTimeStep();
......@@ -246,30 +231,16 @@ int main(int argc, char** argv)
// write vtk output
vtkWriter.write(timeLoop->time());
// calculate and print mass fluxes over the planes
flux.calculateMassOrMoleFluxes();
if(GetPropType<TypeTag, Properties::ModelTraits>::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 solver
timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
// update time
massProblem->setTime(timeLoop->time());
momentumProblem->setTime(timeLoop->time());
} while (!timeLoop->finished());
timeLoop->finalize(leafGridView.comm());
......@@ -287,4 +258,4 @@ int main(int argc, char** argv)
}
return 0;
} // end main
}
......@@ -33,6 +33,9 @@
#include <dumux/freeflow/navierstokes/boundarytypes.hh>
#include <dumux/freeflow/navierstokes/problem.hh>
#include <dumux/freeflow/navierstokes/momentum/fluxhelper.hh>
#include <dumux/freeflow/navierstokes/scalarfluxhelper.hh>
namespace Dumux {
/*!
......@@ -51,20 +54,21 @@ class ChannelTestProblem : public NavierStokesProblem<TypeTag>
{
using ParentType = NavierStokesProblem<TypeTag>;
using BoundaryTypes = Dumux::NavierStokesBoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
using BoundaryTypes = typename ParentType::BoundaryTypes;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
using PrimaryVariables = typename ParentType::PrimaryVariables;
using NumEqVector = typename ParentType::NumEqVector;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
static constexpr auto dimWorld = GridGeometry::GridView::dimensionworld;
using Element = typename FVElementGeometry::Element;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using TimeLoopPtr = std::shared_ptr<CheckPointTimeLoop<Scalar>>;
using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
// the types of outlet boundary conditions
enum class OutletCondition
......@@ -75,30 +79,27 @@ class ChannelTestProblem : public NavierStokesProblem<TypeTag>
public:
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
ChannelTestProblem(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry)
ChannelTestProblem(std::shared_ptr<const GridGeometry> gridGeometry, std::shared_ptr<CouplingManager> couplingManager)
: ParentType(gridGeometry, couplingManager)
{
inletVelocity_ = getParam<Scalar>("Problem.InletVelocity");
dynamicViscosity_ = getParam<Scalar>("Component.LiquidKinematicViscosity", 1.0) * getParam<Scalar>("Component.LiquidDensity", 1.0);
dynamicViscosity_ = getParam<Scalar>("Component.LiquidKinematicViscosity", 1.0)
* getParam<Scalar>("Component.LiquidDensity", 1.0);
const auto tmp = getParam<std::string>("Problem.OutletCondition", "Outflow");
if (tmp == "Outflow")
const auto outletBC = getParam<std::string>("Problem.OutletCondition", "Outflow");
if (outletBC == "Outflow")
outletCondition_ = OutletCondition::outflow;
else if (tmp == "DoNothing")
else if (outletBC == "DoNothing")
outletCondition_ = OutletCondition::doNothing;
else if (tmp == "NeumannX_DirichletY")
else if (outletBC == "NeumannX_DirichletY")
outletCondition_ = OutletCondition::neumannXdirichletY;
else if (tmp == "NeumannX_NeumannY")
else if (outletBC == "NeumannX_NeumannY")
outletCondition_ = OutletCondition::neumannXneumannY;
else
DUNE_THROW(Dune::InvalidStateException, tmp + " is not a valid outlet boundary condition");
DUNE_THROW(Dune::InvalidStateException, outletBC + " is not a valid outlet boundary condition");
useVelocityProfile_ = getParam<bool>("Problem.UseVelocityProfile", false);
outletPressure_ = getParam<Scalar>("Problem.OutletPressure", 1.1e5);
const bool isStationary = getParam<bool>("Problem.IsStationary", false);
hasAnalyticalSolution_ = isStationary && useVelocityProfile_ && (outletCondition_ != OutletCondition::doNothing);
}
/*!
......@@ -121,70 +122,71 @@ public:
*
* \param globalPos The position of the center of the finite volume
*/
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
BoundaryTypes boundaryTypesAtPos(const GlobalPosition& globalPos) const
{
BoundaryTypes values;
if (isInlet_(globalPos))
if constexpr (ParentType::isMomentumProblem())
{
values.setDirichlet(Indices::velocityXIdx);
values.setDirichlet(Indices::velocityYIdx);
#if NONISOTHERMAL
values.setDirichlet(Indices::temperatureIdx);
#endif
}
else if (isOutlet_(globalPos))
{
if (outletCondition_ == OutletCondition::outflow)
values.setDirichlet(Indices::pressureIdx);
else if (outletCondition_ == OutletCondition::doNothing ||
outletCondition_ == OutletCondition::neumannXneumannY)
{
values.setNeumann(Indices::momentumXBalanceIdx);
values.setNeumann(Indices::momentumYBalanceIdx);
}
else // (outletCondition_ == OutletCondition::neumannXdirichletY)
if (isOutlet_(globalPos))
{
values.setDirichlet(Indices::velocityYIdx);
values.setNeumann(Indices::momentumXBalanceIdx);
if (outletCondition_ == OutletCondition::neumannXdirichletY)
{
values.setNeumann(Indices::momentumXBalanceIdx);
values.setDirichlet(Indices::velocityYIdx);
}
else
values.setAllNeumann();
}
#if NONISOTHERMAL
values.setOutflow(Indices::energyEqIdx);
#endif
}
else
{
values.setDirichlet(Indices::velocityXIdx);
values.setDirichlet(Indices::velocityYIdx);
values.setAllNeumann();
if (isInlet_(globalPos))
{
values.setDirichlet(Indices::pressureIdx);
#if NONISOTHERMAL
values.setNeumann(Indices::energyEqIdx);
values.setDirichlet(Indices::temperatureIdx);
#endif
}
}
return values;
}
/*!
* \brief Evaluates 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
* \brief Evaluates the boundary conditions for a Dirichlet control volume.
*
* \param globalPos The center of the finite volume which ought to be set.
*/
PrimaryVariables dirichlet(const Element& element, const SubControlVolumeFace& scvf) const
{
const auto& globalPos = scvf.ipGlobal();
PrimaryVariables values = initialAtPos(globalPos);
if(isInlet_(globalPos))
if constexpr (ParentType::isMomentumProblem())
{
values[Indices::velocityXIdx] = parabolicProfile(globalPos[1], inletVelocity_);
if (!useVelocityProfile_ && isInlet_(globalPos))
values[Indices::velocityXIdx] = inletVelocity_;
}
else
{
values[Indices::pressureIdx] = this->couplingManager().cellPressure(element, scvf);
#if NONISOTHERMAL
// give the system some time so that the pressure can equilibrate, then start the injection of the hot liquid
if(time() >= 200.0)
values[Indices::temperatureIdx] = 293.15;
// give the system some time so that the pressure can equilibrate, then start the injection of the hot liquid
if (time() >= 200.0)
values[Indices::temperatureIdx] = 293.15;
#endif
}
else
values[Indices::velocityXIdx] = 0.0;
return values;
return values;
}
/*!
......@@ -196,23 +198,45 @@ public:
* \param elemFaceVars The element face variables
* \param scvf The boundary sub control volume face
*/
template<class ElementVolumeVariables, class ElementFaceVariables>
template<class ElementVolumeVariables, class ElementFluxVariablesCache>
NumEqVector neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFaceVariables& elemFaceVars,
const ElementFluxVariablesCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
NumEqVector values(0.0);
values[scvf.directionIndex()] = initialAtPos(scvf.center())[Indices::pressureIdx];
if constexpr (ParentType::isMomentumProblem())
{
using FluxHelper = NavierStokesMomentumBoundaryFluxHelper;
// make sure to normalize the pressure if the property is set true
if (getPropValue<TypeTag, Properties::NormalizePressure>())
values[scvf.directionIndex()] -= initialAtPos(scvf.center())[Indices::pressureIdx];
if (outletCondition_ == OutletCondition::doNothing)
values = FluxHelper::fixedPressureMomentumFlux(*this, fvGeometry, scvf, elemVolVars, elemFluxVarsCache, outletPressure_, false /*zeroNormalVelocityGradient*/);
else if (outletCondition_ == OutletCondition::outflow)
values = FluxHelper::fixedPressureMomentumFlux(*this, fvGeometry, scvf, elemVolVars, elemFluxVarsCache, outletPressure_, true /*zeroNormalVelocityGradient*/);
else
{
assert(outletCondition_ == OutletCondition::neumannXneumannY || outletCondition_ == OutletCondition::neumannXdirichletY);
values = FluxHelper::fixedPressureMomentumFlux(*this, fvGeometry, scvf, elemVolVars, elemFluxVarsCache, outletPressure_, false /*zeroNormalVelocityGradient*/);
if (outletCondition_ != OutletCondition::doNothing)
values[1] = -dudy(scvf.center()[1], inletVelocity_) * elemVolVars[scvf.insideScvIdx()].viscosity();
Dune::FieldMatrix<Scalar, dimWorld, dimWorld> shearRate(0.0); // gradV + gradV^T
shearRate[0][1] = dudy(scvf.ipGlobal()[1], inletVelocity_); // we assume du/dx = dv/dy = dv/dx = 0
shearRate[1][0] = shearRate[0][1];
const auto normal = scvf.unitOuterNormal();
NumEqVector normalGradient(0.0);
shearRate.mv(normal, normalGradient);
values -= this->effectiveViscosity(element, fvGeometry, scvf)*normalGradient;
}
}
else
{
using FluxHelper = NavierStokesScalarBoundaryFluxHelper<AdvectiveFlux<ModelTraits>>;
if (isOutlet_(scvf.ipGlobal()))
values = FluxHelper::scalarOutflowFlux(*this, element, fvGeometry, scvf, elemVolVars);
}
return values;
}
......@@ -248,24 +272,25 @@ public: