Commit c993aecc authored by Katharina Heck's avatar Katharina Heck
Browse files

add the neumann option to the angeli test

parent 58008210
......@@ -79,6 +79,7 @@ public:
: ParentType(gridGeometry, couplingManager)
{
kinematicViscosity_ = getParam<Scalar>("Component.LiquidKinematicViscosity", 1.0);
useNeumann_ = getParam<bool>("Problem.UseNeumann", false);
useQuad_ = getParam<bool>("Problem.UseQuad", false);
}
......@@ -105,7 +106,14 @@ public:
BoundaryTypes values;
if constexpr (ParentType::isMomentumProblem())
{
static constexpr Scalar eps = 1e-8;
if (useNeumann_ && (globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps ||
globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps))
values.setAllNeumann();
else
values.setAllDirichlet();
}
else
values.setAllNeumann();
......@@ -134,6 +142,8 @@ public:
// We don't need internal Dirichlet conditions if a Neumann BC is set for the momentum balance (which accounts for the pressure).
// If only Dirichlet BCs are set for the momentum balance, fix the pressure at some cells such that the solution is fully defined.
if (!useNeumann_)
{
auto fvGeometry = localView(this->gridGeometry());
fvGeometry.bindElement(element);
......@@ -146,6 +156,7 @@ public:
if (isAtBoundary(fvGeometry))
values.set(Indices::pressureIdx);
}
return values;
}
......@@ -177,6 +188,51 @@ public:
return analyticalSolution(scvf.center(), 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 insideDensity = elemVolVars[scvf.insideScvIdx()].density();
values[Indices::conti0EqIdx] = this->faceVelocity(element, fvGeometry, scvf) * insideDensity * scvf.unitOuterNormal();
}
else
{
using std::sin;
using std::cos;
Dune::FieldMatrix<Scalar, 2, 2> momentumFlux(0.0);
const auto x = scvf.ipGlobal()[0];
const auto y = scvf.ipGlobal()[1];
const Scalar mu = kinematicViscosity_;
const Scalar t = time_;
momentumFlux[0][0] = M_PI*M_PI *(-4.0*mu*exp(10.0*M_PI*M_PI*mu*t)*sin(M_PI*x)*sin(2*M_PI*y) + (4.0*sin(2*M_PI*y)*sin(2*M_PI*y)*cos(M_PI*x)*cos(M_PI*x) - 1.0*cos(2*M_PI*x) - 0.25*cos(4*M_PI*y))*exp(5.0*M_PI*M_PI*mu*t))*exp(-15.0*M_PI*M_PI*mu*t);
momentumFlux[0][1] = M_PI*M_PI *( 3.0*mu*exp(10.0*M_PI*M_PI*mu*t) - 2.0*exp(5.0*M_PI*M_PI*mu*t)*sin(M_PI*x)*sin(2*M_PI*y))*exp(-15.0*M_PI*M_PI*mu*t)*cos(M_PI*x)*cos(2*M_PI*y);
momentumFlux[1][0] = M_PI*M_PI *( 3.0*mu*exp(10.0*M_PI*M_PI*mu*t) - 2.0*exp(5.0*M_PI*M_PI*mu*t)*sin(M_PI*x)*sin(2*M_PI*y))*exp(-15.0*M_PI*M_PI*mu*t)*cos(M_PI*x)*cos(2*M_PI*y);
momentumFlux[1][1] = M_PI*M_PI *( 4.0*mu*exp(10.0*M_PI*M_PI*mu*t)*sin(M_PI*x)*sin(2*M_PI*y) + (sin(M_PI*x)*sin(M_PI*x)*cos(2*M_PI*y)*cos(2*M_PI*y) - 1.0*cos(2*M_PI*x) - 0.25*cos(4*M_PI*y))*exp(5.0*M_PI*M_PI*mu*t))*exp(-15.0*M_PI*M_PI*mu*t);
const auto normal = scvf.unitOuterNormal();
momentumFlux.mv(normal, values);
}
return values;
}
/*!
* \brief Returns the analytical solution of the problem at a given time and position.
......@@ -282,6 +338,7 @@ private:
Scalar kinematicViscosity_;
Scalar time_ = 0.0;
bool useNeumann_;
bool useQuad_;
};
} // end namespace Dumux
......
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