From 102e1491761e7e97a73ad4d69beccc38903877cb Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Sat, 6 Jan 2018 19:23:37 +0100
Subject: [PATCH] [test][stokes] Include flux calculation in channel test
* calculate the mass and volume fluxes over two different planes
---
test/freeflow/navierstokes/test_channel.cc | 40 ++++++++++++++++++++++
1 file changed, 40 insertions(+)
diff --git a/test/freeflow/navierstokes/test_channel.cc b/test/freeflow/navierstokes/test_channel.cc
index edee5e1a0a..ec1804236c 100644
--- a/test/freeflow/navierstokes/test_channel.cc
+++ b/test/freeflow/navierstokes/test_channel.cc
@@ -52,6 +52,8 @@
#include <dumux/io/staggeredvtkoutputmodule.hh>
+#include <dumux/freeflow/navierstokes/staggered/fluxoverplane.hh>
+
/*!
* \brief Provides an interface for customizing error messages associated with
* reading in parameters.
@@ -177,6 +179,34 @@ int main(int argc, char** argv) try
auto newtonController = std::make_shared<NewtonController>(leafGridView.comm(), timeLoop);
NewtonMethod nonLinearSolver(newtonController, assembler, linearSolver);
+ // set up two planes over which fluxes are calculated
+ FluxOverPlane<TypeTag> flux(*assembler, x);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using GlobalPosition = Dune::FieldVector<Scalar, GridView::dimensionworld>;
+
+ const Scalar xMin = fvGridGeometry->bBoxMin()[0];
+ const Scalar xMax = fvGridGeometry->bBoxMax()[0];
+ const Scalar yMin = fvGridGeometry->bBoxMin()[1];
+ const Scalar yMax = fvGridGeometry->bBoxMax()[1];
+
+ // The first plane 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 postion of the plane (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 plane. This assumes a regular cartesian grid.
+ const Scalar planePosMiddleX = xMin + 0.5*(xMax - xMin);
+ const int numCellsX = getParam<std::vector<int>>("Grid.Cells")[0];
+ 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.addPlane("middle", p0middle, p1middle);
+
+ // The second plane is placed at the outlet of the channel.
+ const auto p0outlet = GlobalPosition{xMax, yMin};
+ const auto p1outlet = GlobalPosition{xMax, yMax};
+ flux.addPlane("outlet", p0outlet, p1outlet);
+
// time loop
timeLoop->start(); do
{
@@ -212,6 +242,16 @@ int main(int argc, char** argv) try
// write vtk output
vtkWriter.write(timeLoop->time());
+ // calculate and print mass fluxes over the planes
+ flux.calculateMassOrMoleFluxes();
+ 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") << std::endl;
+ std::cout << "volume flux at outlet is: " << flux.netFlux("outlet") << std::endl;
+
// report statistics of this time step
timeLoop->reportTimeStep();
--
GitLab