Commit 278351e9 authored by Timo Koch's avatar Timo Koch
Browse files

[test][ff] Port sincos convergence tests to new staggered implementation

parent 3d1bc8d4
......@@ -8,7 +8,7 @@ dumux_add_test(NAME test_ff_navierstokes_sincos
CMAKE_GUARD HAVE_UMFPACK
COMPILE_DEFINITIONS LINEARSOLVER=UMFPackBackend
COMMAND ./convergencetest.py
CMD_ARGS test_ff_navierstokes_sincos params.input
CMD_ARGS test_ff_navierstokes_sincos params.input
-Grid.UpperRight "6.28 6.28"
-Grid.Cells "150 150"
-Problem.Name test_ff_navierstokes_sincos_stationary_convergence
......@@ -38,7 +38,7 @@ dumux_add_test(NAME test_ff_navierstokes_sincos_uzawapreconditioner_factory
LABELS freeflow
TIMEOUT 5000
CMAKE_GUARD "( HAVE_UMFPACK AND ( ( DUNE_ISTL_VERSION VERSION_GREATER 2.7 ) OR ( DUNE_ISTL_VERSION VERSION_EQUAL 2.7 ) ) )"
COMPILE_DEFINITIONS LINEARSOLVER=IstlSolverFactoryBackend<LinearSolverTraits<GridGeometry>>
COMPILE_DEFINITIONS LINEARSOLVER=IstlSolverFactoryBackend<LinearSolverTraits<MassGridGeometry>>
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_ff_navierstokes_sincos_instationary-reference.vtu
......
......@@ -32,59 +32,47 @@
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh>
#include <dune/grid/io/file/vtk.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/io/grid/gridmanager_yasp.hh>
#include <dumux/io/staggeredvtkoutputmodule.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/linear/linearsolvertraits.hh>
#include <dumux/linear/istlsolverfactorybackend.hh>
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/multidomain/fvassembler.hh>
#include <dumux/multidomain/traits.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"
/*!
* \brief Creates analytical solution.
* Returns a tuple of the analytical solution for the pressure, the velocity and the velocity at the faces
* \param time the time at which to evaluate the analytical solution
* \param problem the problem for which to evaluate the analytical solution
*/
template<class Problem>
auto createSource(const Problem& problem)
template<class MomentumProblem>
auto createSource(const MomentumProblem& momentumProblem)
{
using Scalar = double;
using Indices = typename Problem::Indices;
using Indices = typename MomentumProblem::Indices;
const auto& gridGeometry = problem.gridGeometry();
std::array<std::vector<Scalar>, Problem::ModelTraits::numEq()> source;
const auto& gridGeometry = momentumProblem.gridGeometry();
std::array<std::vector<Scalar>, 2> source;
for (auto& component : source)
{
component.resize(gridGeometry.numCellCenterDofs());
}
component.resize(gridGeometry.gridView().size(0));
auto fvGeometry = localView(gridGeometry);
for (const auto& element : elements(gridGeometry.gridView()))
{
fvGeometry.bindElement(element);
for (auto&& scv : scvs(fvGeometry))
{
const auto ccDofIdx = scv.dofIndex();
const auto ccDofPosition = scv.dofPosition();
const auto sourceAtPosVal = problem.sourceAtPos(ccDofPosition);
source[Indices::momentumXBalanceIdx][ccDofIdx] = sourceAtPosVal[Indices::momentumXBalanceIdx];
source[Indices::momentumYBalanceIdx][ccDofIdx] = sourceAtPosVal[Indices::momentumYBalanceIdx];
}
const auto& center = element.geometry().center();
const auto eIdx = gridGeometry.elementMapper().index(element);
auto sourceAtPosVal = momentumProblem.sourceAtPos(center);
source[Indices::momentumXBalanceIdx][eIdx] = sourceAtPosVal[Indices::momentumXBalanceIdx];
source[Indices::momentumYBalanceIdx][eIdx] = sourceAtPosVal[Indices::momentumYBalanceIdx];
}
return source;
......@@ -95,7 +83,8 @@ int main(int argc, char** argv)
using namespace Dumux;
// define the type tag for this problem
using TypeTag = Properties::TTag::SincosTest;
using MomentumTypeTag = Properties::TTag::SincosTestMomentum;
using MassTypeTag = Properties::TTag::SincosTestMass;
// initialize MPI, finalize is done automatically on exit
const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
......@@ -108,7 +97,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<MomentumTypeTag, Properties::Grid>> gridManager;
gridManager.init();
////////////////////////////////////////////////////////////
......@@ -119,11 +108,18 @@ 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>();
// get some time loop parameters
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Scalar = GetPropType<MomentumTypeTag, Properties::Scalar>;
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
auto dt = getParam<Scalar>("TimeLoop.DtInitial");
......@@ -133,61 +129,80 @@ int main(int argc, char** argv)
timeLoop->setMaxTimeStepSize(maxDt);
// the problem (initial and boundary conditions)
using Problem = GetPropType<TypeTag, Properties::Problem>;
auto problem = std::make_shared<Problem>(gridGeometry);
problem->updateTimeStepSize(timeLoop->timeStepSize());
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>;
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);
momentumProblem->applyInitialSolution(x[momentumIdx]);
massProblem->applyInitialSolution(x[massIdx]);
auto xOld = x;
// 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);
const bool isStationary = getParam<bool>("Problem.IsStationary");
if (isStationary)
couplingManager->init(momentumProblem, massProblem, std::make_tuple(momentumGridVariables, massGridVariables), x);
else
couplingManager->init(momentumProblem, massProblem, std::make_tuple(momentumGridVariables, massGridVariables), x, xOld);
massGridVariables->init(x[massIdx]);
momentumGridVariables->init(x[momentumIdx]);
// initialize the vtk output module
StaggeredVtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
using IOFields = GetPropType<TypeTag, Properties::IOFields>;
using IOFields = GetPropType<MassTypeTag, Properties::IOFields>;
VtkOutputModule vtkWriter(*massGridVariables, x[massIdx], massProblem->name());
IOFields::initOutputModule(vtkWriter); // Add model specific output fields
auto source = createSource(*problem);
auto sourceX = source[Problem::Indices::momentumXBalanceIdx];
auto sourceY = source[Problem::Indices::momentumYBalanceIdx];
vtkWriter.addVelocityOutput(std::make_shared<NavierStokesVelocityOutput<MassGridVariables>>());
NavierStokesTest::AnalyticalSolutionVectors analyticalSolVectors(momentumProblem, massProblem);
vtkWriter.addField(analyticalSolVectors.analyticalPressureSolution(), "pressureExact");
vtkWriter.addField(analyticalSolVectors.analyticalVelocitySolution(), "velocityExact");
const auto source = createSource(*momentumProblem);
const auto& sourceX = source[MomentumProblem::Indices::momentumXBalanceIdx];
const auto& sourceY = source[MomentumProblem::Indices::momentumYBalanceIdx];
vtkWriter.addField(sourceX, "sourceX");
vtkWriter.addField(sourceY, "sourceY");
NavierStokesAnalyticalSolutionVectors analyticalSolVectors(problem, 0.0);
vtkWriter.addField(analyticalSolVectors.getAnalyticalPressureSolution(), "pressureExact");
vtkWriter.addField(analyticalSolVectors.getAnalyticalVelocitySolution(), "velocityExact");
vtkWriter.addFaceField(analyticalSolVectors.getAnalyticalVelocitySolutionOnFace(), "faceVelocityExact");
vtkWriter.write(0.0);
const bool isStationary = getParam<bool>("Problem.IsStationary");
// 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);
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 = LINEARSOLVER;
auto linearSolver = std::make_shared<LinearSolver>();
// the non-linear solver
using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);
using NewtonSolver = Dumux::MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
// the discrete L2 and Linfity errors
const bool printErrors = getParam<bool>("Problem.PrintErrors", false);
const bool printConvergenceTestFile = getParam<bool>("Problem.PrintConvergenceTestFile", false);
NavierStokesErrors errors(problem, x);
NavierStokesErrorCSVWriter errorCSVWriter(problem);
NavierStokesTest::Errors errors(momentumProblem, massProblem, x);
NavierStokesTest::ErrorCSVWriter errorCSVWriter(momentumProblem, massProblem);
if (isStationary)
{
......@@ -202,7 +217,7 @@ int main(int argc, char** argv)
errorCSVWriter.printErrors(errors);
if (printConvergenceTestFile)
convergenceTestAppendErrors(problem, errors);
convergenceTestAppendErrors(momentumProblem, massProblem, errors);
}
// write vtk output
......@@ -221,12 +236,17 @@ int main(int argc, char** argv)
// time loop
timeLoop->start(); do
{
// set the correct time level for the problem's boundary conditions
momentumProblem->updateTime(timeLoop->time() + timeLoop->timeStepSize());
massProblem->updateTime(timeLoop->time() + timeLoop->timeStepSize());
// solve the non-linear system with time step control
nonLinearSolver.solve(x, *timeLoop);
// make the new solution the old solution
xOld = x;
gridVariables->advanceTimeStep();
momentumGridVariables->advanceTimeStep();
massGridVariables->advanceTimeStep();
// print discrete L2 and Linfity errors
if (printErrors)
......@@ -237,7 +257,6 @@ int main(int argc, char** argv)
// advance to the time loop to the next step
timeLoop->advanceTimeStep();
problem->updateTime(timeLoop->time());
analyticalSolVectors.update(timeLoop->time());
// write vtk output
......@@ -248,7 +267,6 @@ int main(int argc, char** argv)
// set new dt as suggested by newton solver
timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
problem->updateTimeStepSize(timeLoop->timeStepSize());
} while (!timeLoop->finished());
......
......@@ -47,32 +47,33 @@ class SincosTestProblem : 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 PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
using NumEqVector = typename ParentType::NumEqVector;
using PrimaryVariables = typename ParentType::PrimaryVariables;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
using SubControlVolume = typename GridGeometry::SubControlVolume;
using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
using FVElementGeometry = typename GridGeometry::LocalView;
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 CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
public:
using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
SincosTestProblem(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry), time_(0.0), timeStepSize_(0.0)
SincosTestProblem(std::shared_ptr<const GridGeometry> gridGeometry, std::shared_ptr<CouplingManager> couplingManager)
: ParentType(gridGeometry, couplingManager), time_(0.0)
{
isStationary_ = getParam<bool>("Problem.IsStationary");
enableInertiaTerms_ = getParam<bool>("Problem.EnableInertiaTerms");
rho_ = getParam<Scalar>("Component.LiquidDensity");
//kinematic
Scalar nu = getParam<Scalar>("Component.LiquidKinematicViscosity", 1.0);
//dynamic
mu_ = rho_*nu;
mu_ = rho_*nu; // dynamic viscosity
useNeumann_ = getParam<bool>("Problem.UseNeumann", false);
}
/*!
......@@ -83,7 +84,7 @@ public:
Scalar temperature() const
{ return 298.0; }
/*!
/*!
* \brief Return the sources within the domain.
*
* \param globalPos The global position
......@@ -91,20 +92,20 @@ public:
NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
{
NumEqVector source(0.0);
const Scalar x = globalPos[0];
const Scalar y = globalPos[1];
const Scalar t = time_ + timeStepSize_;
using std::cos;
using std::sin;
source[Indices::momentumXBalanceIdx] = rho_*dtU_(x,y,t) - 2.*mu_*dxxU_(x,y,t) - mu_*dyyU_(x,y,t) - mu_*dxyV_(x,y,t) + dxP_(x,y,t);
source[Indices::momentumYBalanceIdx] = rho_*dtV_(x,y,t) - 2.*mu_*dyyV_(x,y,t) - mu_*dxyU_(x,y,t) - mu_*dxxV_(x,y,t) + dyP_(x,y,t);
if (enableInertiaTerms_)
if constexpr (ParentType::isMomentumProblem())
{
source[Indices::momentumXBalanceIdx] += rho_*dxUU_(x,y,t) + rho_*dyUV_(x,y,t);
source[Indices::momentumYBalanceIdx] += rho_*dxUV_(x,y,t) + rho_*dyVV_(x,y,t);
const Scalar x = globalPos[0];
const Scalar y = globalPos[1];
const Scalar t = time_;
source[Indices::momentumXBalanceIdx] = rho_*dtU_(x,y,t) - 2.0*mu_*dxxU_(x,y,t) - mu_*dyyU_(x,y,t) - mu_*dxyV_(x,y,t) + dxP_(x,y,t);
source[Indices::momentumYBalanceIdx] = rho_*dtV_(x,y,t) - 2.0*mu_*dyyV_(x,y,t) - mu_*dxyU_(x,y,t) - mu_*dxxV_(x,y,t) + dyP_(x,y,t);
if (this->enableInertiaTerms())
{
source[Indices::momentumXBalanceIdx] += rho_*dxUU_(x,y,t) + rho_*dyUV_(x,y,t);
source[Indices::momentumYBalanceIdx] += rho_*dxUV_(x,y,t) + rho_*dyVV_(x,y,t);
}
}
return source;
......@@ -116,49 +117,139 @@ public:
*/
// \{
/*!
/*!
* \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 boundaryTypesAtPos(const GlobalPosition& globalPos) const
{
BoundaryTypes values;
// set Dirichlet values for the velocity everywhere
values.setDirichlet(Indices::velocityXIdx);
values.setDirichlet(Indices::velocityYIdx);
if constexpr (ParentType::isMomentumProblem())
{
if (useNeumann_)
{
if (globalPos[1] > this->gridGeometry().bBoxMax()[1] - 1e-6
|| globalPos[0] > this->gridGeometry().bBoxMax()[0] - 1e-6
|| globalPos[1] < this->gridGeometry().bBoxMin()[1] + 1e-6)
{
values.setNeumann(Indices::velocityXIdx);
values.setNeumann(Indices::velocityYIdx);
}
else
{
values.setDirichlet(Indices::velocityXIdx);
values.setDirichlet(Indices::velocityYIdx);
}
}
else
{
// set Dirichlet values for the velocity everywhere
values.setAllDirichlet();
}
}
else
values.setAllNeumann();
return values;
}
//! Enable internal Dirichlet constraints
static constexpr bool enableInternalDirichletConstraints()
{ return !ParentType::isMomentumProblem(); }
/*!
* \brief Returns whether a fixed Dirichlet value shall be used at a given cell.
* \brief Tag a degree of freedom to carry internal Dirichlet constraints.
* If true is returned for a dof, the equation for this dof is replaced
* by the constraint that its primary variable values must match the
* user-defined values obtained from the function internalDirichlet(),
* which must be defined in the problem.
*
* \param element The finite element
* \param fvGeometry The finite-volume geometry
* \param scv The sub control volume
* \param pvIdx The primary variable index in the solution vector
* \param scv The sub-control volume
*/
bool isDirichletCell(const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolume& scv,
int pvIdx) const
std::bitset<PrimaryVariables::dimension> hasInternalDirichletConstraint(const Element& element, const SubControlVolume& scv) const
{
// set fixed pressure in one cell
return (scv.dofIndex() == 0) && pvIdx == Indices::pressureIdx;
std::bitset<PrimaryVariables::dimension> values;
if (!useNeumann_ && scv.dofIndex() == 0)
values.set(Indices::pressureIdx);
return values;
}
/*!
/*!
* \brief Define the values of internal Dirichlet constraints for a degree of freedom.
* \param element The finite element
* \param scv The sub-control volume
*/
PrimaryVariables internalDirichlet(const Element& element, const SubControlVolume& scv) const
{ return PrimaryVariables(analyticalSolution(scv.center())[Indices::pressureIdx]); }
/*!
* \brief Returns Dirichlet boundary values at a given position.
*
* \param globalPos The global position
*/
PrimaryVariables dirichletAtPos(const GlobalPosition & globalPos) const
PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
{
// use the values of the analytical solution
return analyticalSolution(globalPos, time_+timeStepSize_);
return analyticalSolution(globalPos, time_);
}
/*!
* \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);
if constexpr (ParentType::isMomentumProblem())
{
const auto flux = [&](Scalar x, Scalar y)
{
Dune::FieldMatrix<Scalar, dimWorld, dimWorld> momentumFlux(0.0);
const Scalar t = time_;
momentumFlux[0][0] = -2*mu_*dxU_(x,y,t) + p_(x,y,t);
momentumFlux[0][1] = -mu_*(dyU_(x,y,t) + dxV_(x,y,t));
momentumFlux[1][0] = -mu_*(dyU_(x,y,t) + dxV_(x,y,t));
momentumFlux[1][1] = -2*mu_*dyV_(x,y,t) + p_(x,y,t);
if (this->enableInertiaTerms())
{
momentumFlux[0][0] += rho_*u_(x,y,t)*u_(x,y,t);
momentumFlux[0][1] += rho_*u_(x,y,t)*v_(x,y,t);
momentumFlux[1][0] += rho_*v_(x,y,t)*u_(x,y,t);
momentumFlux[1][1] += rho_*v_(x,y,t)*v_(x,y,t);
}
return momentumFlux;
};
flux(scvf.ipGlobal()[0], scvf.ipGlobal()[1]).mv(scvf.unitOuterNormal(), values);
}
else
{
const auto insideDensity = elemVolVars[scvf.insideScvIdx()].density();
values[Indices::conti0EqIdx] = insideDensity * (this->faceVelocity(element, fvGeometry, scvf) * scvf.unitOuterNormal());
}
return values;
}
/*!
......@@ -171,14 +262,16 @@ public:
{
const Scalar x = globalPos[0];
const Scalar y = globalPos[1];
const Scalar t = time;
PrimaryVariables values;
values[Indices::pressureIdx] = (f1_(x) + f1_(y)) * f_(t) * f_(t);
values[Indices::velocityXIdx] = u_(x,y,t);
values[Indices::velocityYIdx] = v_(x,y,t);
if constexpr (ParentType::isMomentumProblem())
{
values[Indices::velocityXIdx] = u_(x,y,t);
values[Indices::velocityYIdx] = v_(x,y,t);
}
else
values[Indices::pressureIdx] = p_(x,y,t);
return values;
}
......@@ -190,7 +283,7 @@ public:
*/
// \{
/*!
/*!
* \brief Evaluates the initial value for a control volume.
*
* \param globalPos The global position
......@@ -198,70 +291,69 @@ public:
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
if (isStationary_)
{
PrimaryVariables values;
values[Indices::pressureIdx] = 0.0;
values[Indices::velocityXIdx] = 0.0;
values[Indices::velocityYIdx] = 0.0;
return values;
}
return PrimaryVariables(0.0);
else
{
return analyticalSolution(globalPos, 0.0);
}
}
/*!
* \brief Updates the time
* \brief Returns the analytical solution of the problem at a given position.