Newer
Older
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
*
* \brief Channel flow test for the staggered grid (Navier-)Stokes model
*/
#include <config.h>
#include <ctime>
#include <iostream>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh>
#include <dune/grid/io/file/dgfparser/dgfexception.hh>
#include <dune/grid/io/file/vtk.hh>
#include <dune/istl/io.hh>
#include "channeltestproblem.hh"
#include <dumux/common/parameters.hh>
#include <dumux/common/valgrind.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/common/defaultusagemessage.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/assembly/staggeredfvassembler.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/discretization/methods.hh>
#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.
*
* \param progName The name of the program, that was tried to be started.
* \param errorMsg The error message that was issued by the start function.
* Comprises the thing that went wrong and a general help message.
*/
void usage(const char *progName, const std::string &errorMsg)
{
if (errorMsg.size() > 0) {
std::string errorMessageOut = "\nUsage: ";
errorMessageOut += progName;
errorMessageOut += " [options]\n";
errorMessageOut += errorMsg;
errorMessageOut += "\n\nThe list of mandatory arguments for this program is:\n"
"\t-TimeManager.TEnd End of the simulation [s] \n"
"\t-TimeManager.DtInitial Initial timestep size [s] \n"
"\t-Grid.File Name of the file containing the grid \n"
"\t definition in DGF format\n"
"\t-SpatialParams.LensLowerLeftX x-coordinate of the lower left corner of the lens [m] \n"
"\t-SpatialParams.LensLowerLeftY y-coordinate of the lower left corner of the lens [m] \n"
"\t-SpatialParams.LensUpperRightX x-coordinate of the upper right corner of the lens [m] \n"
"\t-SpatialParams.LensUpperRightY y-coordinate of the upper right corner of the lens [m] \n"
"\t-SpatialParams.Permeability Permeability of the domain [m^2] \n"
"\t-SpatialParams.PermeabilityLens Permeability of the lens [m^2] \n";
std::cout << errorMessageOut
<< "\n";
}
}
int main(int argc, char** argv) try
{
using namespace Dumux;
// define the type tag for this problem
using TypeTag = TTAG(ChannelTestTypeTag);
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
// initialize MPI, finalize is done automatically on exit
const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
// print dumux start message
if (mpiHelper.rank() == 0)
DumuxMessage::print(/*firstCall=*/true);
// parse command line arguments and input file
Parameters::init(argc, argv, usage);
// try to create a grid (from the given grid file or the input file)
using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator);
GridCreator::makeGrid();
GridCreator::loadBalance();
////////////////////////////////////////////////////////////
// run instationary non-linear problem on this grid
////////////////////////////////////////////////////////////
// we compute on the leaf grid view
const auto& leafGridView = GridCreator::grid().leafGridView();
// create the finite volume grid geometry
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
fvGridGeometry->update();
// the problem (initial and boundary conditions)
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
auto problem = std::make_shared<Problem>(fvGridGeometry);
// get some time loop parameters
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
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");
// instantiate time loop
auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(restartTime, dt, tEnd);
timeLoop->setMaxTimeStepSize(maxDt);
problem->setTimeLoop(timeLoop);
// the solution vector
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
typename DofTypeIndices::CellCenterIdx cellCenterIdx;
typename DofTypeIndices::FaceIdx faceIdx;
const auto numDofsCellCenter = leafGridView.size(0);
const auto numDofsFace = leafGridView.size(1);
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);
// intialize the vtk output module
using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
StaggeredVtkOutputModule<TypeTag, GET_PROP_VALUE(TypeTag, PhaseIdx)> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
VtkOutputFields::init(vtkWriter); //!< Add model specific output fields
vtkWriter.write(0.0);
// the assembler with time loop for instationary problem
using Assembler = StaggeredFVAssembler<TypeTag, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
// 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);
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
// 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
{
// set previous solution for storage evaluations
assembler->setPreviousSolution(xOld);
// solve the non-linear system with time step control
nonLinearSolver.solve(x, *timeLoop);
// make the new solution the old solution
xOld = x;
gridVariables->advanceTimeStep();
// advance to the time loop to the next step
timeLoop->advanceTimeStep();
// write vtk output
vtkWriter.write(timeLoop->time());
// calculate and print mass fluxes over the planes
flux.calculateMassOrMoleFluxes();
if(GET_PROP_TYPE(TypeTag, 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()));
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
} while (!timeLoop->finished());
timeLoop->finalize(leafGridView.comm());
////////////////////////////////////////////////////////////
// finalize, print dumux message to say goodbye
////////////////////////////////////////////////////////////
// print dumux end message
if (mpiHelper.rank() == 0)
{
Parameters::print();
DumuxMessage::print(/*firstCall=*/false);
}
return 0;
} // end main
catch (Dumux::ParameterException &e)
{
std::cerr << std::endl << e << " ---> Abort!" << std::endl;
return 1;
}
catch (Dune::DGFException & e)
{
std::cerr << "DGF exception thrown (" << e <<
"). Most likely, the DGF file name is wrong "
"or the DGF file is corrupted, "
"e.g. missing hash at end of file or wrong number (dimensions) of entries."
<< " ---> Abort!" << std::endl;
return 2;
}
catch (Dune::Exception &e)
{
std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
return 3;
}
catch (...)
std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
return 4;