Commit af8308dd authored by Kilian Weishaupt's avatar Kilian Weishaupt Committed by Timo Koch
Browse files

[navierstokes][test] Adapt test_kovasznay to changes

* make stationary
parent b007df5a
......@@ -77,7 +77,7 @@ dune_add_test(NAME test_navierstokes_kovasznay
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_kovasznay-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_kovasznay-00002.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_kovasznay-00001.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_navierstokes_kovasznay")
dune_add_test(NAME test_navierstokes_angeli
......
......@@ -32,22 +32,11 @@
#include <dumux/discretization/staggered/freeflow/properties.hh>
#include <dumux/freeflow/navierstokes/model.hh>
// solve Navier-Stokes equations
#define ENABLE_NAVIERSTOKES 1
namespace Dumux
{
template <class TypeTag>
class KovasznayTestProblem;
namespace Capabilities
{
template<class TypeTag>
struct isStationary<KovasznayTestProblem<TypeTag>>
{ static const bool value = true; };
}
namespace Properties
{
NEW_TYPE_TAG(KovasznayTestProblem, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokes));
......@@ -70,11 +59,7 @@ SET_BOOL_PROP(KovasznayTestProblem, EnableFVGridGeometryCache, true);
SET_BOOL_PROP(KovasznayTestProblem, EnableGridFluxVariablesCache, true);
SET_BOOL_PROP(KovasznayTestProblem, EnableGridVolumeVariablesCache, true);
#if ENABLE_NAVIERSTOKES
SET_BOOL_PROP(KovasznayTestProblem, EnableInertiaTerms, true);
#else
SET_BOOL_PROP(KovasznayTestProblem, EnableInertiaTerms, false);
#endif
}
/*!
......@@ -299,10 +284,10 @@ public:
// treat cell-center dofs
const auto dofIdxCellCenter = scv.dofIndex();
const auto& posCellCenter = scv.dofPosition();
const auto analyticalSolutionCellCenter = dirichletAtPos(posCellCenter)[cellCenterIdx];
const auto numericalSolutionCellCenter = curSol[cellCenterIdx][dofIdxCellCenter];
sumError[cellCenterIdx] += squaredDiff_(analyticalSolutionCellCenter, numericalSolutionCellCenter) * scv.volume();
sumReference[cellCenterIdx] += analyticalSolutionCellCenter * analyticalSolutionCellCenter * scv.volume();
const auto analyticalSolutionCellCenter = dirichletAtPos(posCellCenter)[pressureIdx];
const auto numericalSolutionCellCenter = curSol[cellCenterIdx][dofIdxCellCenter][pressureIdx];
sumError[pressureIdx] += squaredDiff_(analyticalSolutionCellCenter, numericalSolutionCellCenter) * scv.volume();
sumReference[pressureIdx] += analyticalSolutionCellCenter * analyticalSolutionCellCenter * scv.volume();
totalVolume += scv.volume();
// treat face dofs
......@@ -310,7 +295,7 @@ public:
{
const int dofIdxFace = scvf.dofIndex();
const int dirIdx = scvf.directionIndex();
const auto analyticalSolutionFace = dirichletAtPos(scvf.center())[faceIdx][dirIdx];
const auto analyticalSolutionFace = dirichletAtPos(scvf.center())[Indices::velocity(dirIdx)];
const auto numericalSolutionFace = curSol[faceIdx][dofIdxFace][momentumBalanceIdx];
directionIndex[dofIdxFace] = dirIdx;
errorVelocity[dofIdxFace] = squaredDiff_(analyticalSolutionFace, numericalSolutionFace);
......@@ -322,8 +307,8 @@ public:
}
// get the absolute and relative discrete L2-error for cell-center dofs
l2NormAbs[cellCenterIdx] = std::sqrt(sumError[cellCenterIdx] / totalVolume);
l2NormRel[cellCenterIdx] = std::sqrt(sumError[cellCenterIdx] / sumReference[cellCenterIdx]);
l2NormAbs[pressureIdx] = std::sqrt(sumError[pressureIdx] / totalVolume);
l2NormRel[pressureIdx] = std::sqrt(sumError[pressureIdx] / sumReference[pressureIdx]);
// get the absolute and relative discrete L2-error for face dofs
for(int i = 0; i < numFaceDofs; ++i)
......@@ -332,14 +317,14 @@ public:
const auto error = errorVelocity[i];
const auto ref = velocityReference[i];
const auto volume = staggeredVolume[i];
sumError[faceIdx][dirIdx] += error * volume;
sumReference[faceIdx][dirIdx] += ref * volume;
sumError[Indices::velocity(dirIdx)] += error * volume;
sumReference[Indices::velocity(dirIdx)] += ref * volume;
}
for(int dirIdx = 0; dirIdx < dimWorld; ++dirIdx)
{
l2NormAbs[faceIdx][dirIdx] = std::sqrt(sumError[faceIdx][dirIdx] / totalVolume);
l2NormRel[faceIdx][dirIdx] = std::sqrt(sumError[faceIdx][dirIdx] / sumReference[faceIdx][dirIdx]);
l2NormAbs[Indices::velocity(dirIdx)] = std::sqrt(sumError[Indices::velocity(dirIdx)] / totalVolume);
l2NormRel[Indices::velocity(dirIdx)] = std::sqrt(sumError[Indices::velocity(dirIdx)] / sumReference[Indices::velocity(dirIdx)]);
}
return std::make_pair(l2NormAbs, l2NormRel);
}
......@@ -396,11 +381,13 @@ private:
const auto faceDofPosition = scvf.center();
const auto dirIdx = scvf.directionIndex();
const auto analyticalSolutionAtFace = analyticalSolution(faceDofPosition);
analyticalVelocityOnFace_[faceDofIdx][dirIdx] = analyticalSolutionAtFace[faceIdx][dirIdx];
analyticalVelocityOnFace_[faceDofIdx][dirIdx] = analyticalSolutionAtFace[Indices::velocity(dirIdx)];
}
analyticalPressure_[ccDofIdx] = analyticalSolutionAtCc[pressureIdx];
analyticalVelocity_[ccDofIdx] = analyticalSolutionAtCc[faceIdx];
for(int dirIdx = 0; dirIdx < dim; ++dirIdx)
analyticalVelocity_[ccDofIdx][dirIdx] = analyticalSolutionAtCc[Indices::velocity(dirIdx)];
}
}
}
......
......@@ -19,7 +19,7 @@
/*!
* \file
*
* \brief Test for the staggered grid Stokes model (Donea et al., 2003)
* \brief Stationary test for the staggered grid Navier-Stokes model (Kovasznay 1947)
*/
#include <config.h>
......@@ -119,7 +119,7 @@ int main(int argc, char** argv) try
auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
fvGridGeometry->update();
// the problem (initial and boundary conditions)
// the problem (boundary conditions)
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
auto problem = std::make_shared<Problem>(fvGridGeometry);
......@@ -133,25 +133,11 @@ int main(int argc, char** argv) try
SolutionVector x;
x[cellCenterIdx].resize(numDofsCellCenter);
x[faceIdx].resize(numDofsFace);
problem->applyInitialSolution(x);
auto xOld = x;
// the grid variables
using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
gridVariables->init(x, xOld);
// get some time loop parameters
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
const auto maxDivisions = getParam<int>("TimeLoop.MaxTimeStepDivisions");
const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
auto dt = getParam<Scalar>("TimeLoop.DtInitial");
// check if we are about to restart a previously interrupted simulation
Scalar restartTime = 0;
if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
restartTime = getParam<Scalar>("TimeLoop.Restart");
gridVariables->init(x);
// intialize the vtk output module
using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
......@@ -162,13 +148,9 @@ int main(int argc, char** argv) try
vtkWriter.addFaceField(problem->getAnalyticalVelocitySolutionOnFace(), "faceVelocityExact");
vtkWriter.write(0.0);
// instantiate time loop
auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, dt, tEnd);
timeLoop->setMaxTimeStepSize(maxDt);
// the assembler with time loop for instationary problem
using Assembler = StaggeredFVAssembler<TypeTag, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables);
// the linear solver
using LinearSolver = Dumux::UMFPackBackend<TypeTag>;
......@@ -177,55 +159,23 @@ int main(int argc, char** argv) try
// the non-linear solver
using NewtonController = StaggeredNewtonController<TypeTag>;
using NewtonMethod = Dumux::NewtonMethod<NewtonController, Assembler, LinearSolver>;
auto newtonController = std::make_shared<NewtonController>(leafGridView.comm(), timeLoop);
auto newtonController = std::make_shared<NewtonController>(leafGridView.comm());
NewtonMethod nonLinearSolver(newtonController, assembler, linearSolver);
// time loop
timeLoop->start(); do
{
// set previous solution for storage evaluations
assembler->setPreviousSolution(xOld);
// try solving the non-linear system
for (int i = 0; i < maxDivisions; ++i)
{
// linearize & solve
auto converged = nonLinearSolver.solve(x);
if (converged)
break;
if (!converged && i == maxDivisions-1)
DUNE_THROW(Dune::MathError,
"Newton solver didn't converge after "
<< maxDivisions
<< " time-step divisions. dt="
<< timeLoop->timeStepSize()
<< ".\nThe solutions of the current and the previous time steps "
<< "have been saved to restart files.");
}
// make the new solution the old solution
xOld = x;
gridVariables->advanceTimeStep();
// advance to the time loop to the next step
timeLoop->advanceTimeStep();
problem->postTimeStep(x);
// write vtk output
vtkWriter.write(timeLoop->time());
// linearize & solve
Dune::Timer timer;
nonLinearSolver.solve(x);
// report statistics of this time step
timeLoop->reportTimeStep();
// write vtk output
vtkWriter.write(1.0);
// set new dt as suggested by newton controller
timeLoop->setTimeStepSize(newtonController->suggestTimeStepSize(timeLoop->timeStepSize()));
timer.stop();
} while (!timeLoop->finished());
const auto& comm = Dune::MPIHelper::getCollectiveCommunication();
std::cout << "Simulation took " << timer.elapsed() << " seconds on "
<< comm.size() << " processes.\n"
<< "The cumulative CPU time was " << timer.elapsed()*comm.size() << " seconds.\n";
timeLoop->finalize(leafGridView.comm());
////////////////////////////////////////////////////////////
// finalize, print dumux message to say goodbye
......
[TimeLoop]
DtInitial = 1 # [s]
TEnd = 2 # [s]
[Grid]
LowerLeft = -0.5 -0.5
UpperRight = 2 1.5
......
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