main.cc 11.5 KB
Newer Older
1 2 3 4 5 6 7
// -*- 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    *
8
 *   the Free Software Foundation, either version 3 of the License, or       *
9 10 11 12 13 14 15 16 17 18 19 20
 *   (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
21
 * \ingroup NavierStokesTests
22 23
 * \brief Test for the instationary staggered grid Navier-Stokes model
 *        with analytical solution (Angeli et al. 2017, \cite Angeli2017).
24
 */
25

26
#include <config.h>
27

28 29
#include <ctime>
#include <iostream>
30 31
#include <type_traits>
#include <tuple>
32

33 34 35 36
#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>
37

38 39 40
#include <dumux/assembly/staggeredfvassembler.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/common/dumuxmessage.hh>
41
#include <dumux/common/parameters.hh>
42
#include <dumux/common/properties.hh>
43
#include <dumux/common/valgrind.hh>
44 45
#include <dumux/io/grid/gridmanager.hh>
#include <dumux/io/staggeredvtkoutputmodule.hh>
46
#include <dumux/linear/seqsolverbackend.hh>
47
#include <dumux/nonlinear/newtonsolver.hh>
48

49
#include "problem.hh"
50

51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106
/*!
* \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 Scalar, class Problem>
auto createAnalyticalSolution(const Scalar time, const Problem& problem)
{
    const auto& fvGridGeometry = problem.fvGridGeometry();
    using GridView = typename std::decay_t<decltype(fvGridGeometry)>::GridView;

    static constexpr auto dim = GridView::dimension;
    static constexpr auto dimWorld = GridView::dimensionworld;

    using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;

    std::vector<Scalar> analyticalPressure;
    std::vector<VelocityVector> analyticalVelocity;
    std::vector<VelocityVector> analyticalVelocityOnFace;

    analyticalPressure.resize(fvGridGeometry.numCellCenterDofs());
    analyticalVelocity.resize(fvGridGeometry.numCellCenterDofs());
    analyticalVelocityOnFace.resize(fvGridGeometry.numFaceDofs());

    using Indices = typename Problem::Indices;
    for (const auto& element : elements(fvGridGeometry.gridView()))
    {
        auto fvGeometry = localView(fvGridGeometry);
        fvGeometry.bindElement(element);
        for (auto&& scv : scvs(fvGeometry))
        {
            auto ccDofIdx = scv.dofIndex();
            auto ccDofPosition = scv.dofPosition();
            auto analyticalSolutionAtCc = problem.analyticalSolution(ccDofPosition, time);

            // velocities on faces
            for (auto&& scvf : scvfs(fvGeometry))
            {
                const auto faceDofIdx = scvf.dofIndex();
                const auto faceDofPosition = scvf.center();
                const auto dirIdx = scvf.directionIndex();
                const auto analyticalSolutionAtFace = problem.analyticalSolution(faceDofPosition, time);
                analyticalVelocityOnFace[faceDofIdx][dirIdx] = analyticalSolutionAtFace[Indices::velocity(dirIdx)];
            }

            analyticalPressure[ccDofIdx] = analyticalSolutionAtCc[Indices::pressureIdx];

            for(int dirIdx = 0; dirIdx < dim; ++dirIdx)
                analyticalVelocity[ccDofIdx][dirIdx] = analyticalSolutionAtCc[Indices::velocity(dirIdx)];
        }
    }

    return std::make_tuple(analyticalPressure, analyticalVelocity, analyticalVelocityOnFace);
}

107 108 109 110 111
int main(int argc, char** argv) try
{
    using namespace Dumux;

    // define the type tag for this problem
112
    using TypeTag = Properties::TTag::AngeliTest;
113 114 115 116 117 118 119 120 121

    // 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
122
    Parameters::init(argc, argv);
123 124

    // try to create a grid (from the given grid file or the input file)
125
    GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
126
    gridManager.init();
127 128 129 130 131 132

    ////////////////////////////////////////////////////////////
    // run instationary non-linear problem on this grid
    ////////////////////////////////////////////////////////////

    // we compute on the leaf grid view
133
    const auto& leafGridView = gridManager.grid().leafGridView();
134 135

    // create the finite volume grid geometry
136
    using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
137 138 139 140
    auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
    fvGridGeometry->update();

    // get some time loop parameters
141
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
142 143
    const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
    const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
144
    auto dt = getParam<Scalar>("TimeLoop.DtInitial");
145 146

    // instantiate time loop
147
    auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0, dt, tEnd);
148 149 150
    timeLoop->setMaxTimeStepSize(maxDt);

    // the problem (initial and boundary conditions)
151
    using Problem = GetPropType<TypeTag, Properties::Problem>;
152
    auto problem = std::make_shared<Problem>(fvGridGeometry);
153
    problem->updateTimeStepSize(timeLoop->timeStepSize());
154 155

    // the solution vector
156
    using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
157
    SolutionVector x;
158 159
    x[FVGridGeometry::cellCenterIdx()].resize(fvGridGeometry->numCellCenterDofs());
    x[FVGridGeometry::faceIdx()].resize(fvGridGeometry->numFaceDofs());
160 161 162 163
    problem->applyInitialSolution(x);
    auto xOld = x;

    // the grid variables
164
    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
165
    auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
166
    gridVariables->init(x);
167

168
    // initialize the vtk output module
169
    StaggeredVtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
170
    using IOFields = GetPropType<TypeTag, Properties::IOFields>;
171
    IOFields::initOutputModule(vtkWriter); // Add model specific output fields
172 173 174 175 176

    auto analyticalSolution = createAnalyticalSolution(timeLoop->time(), *problem);
    vtkWriter.addField(std::get<0>(analyticalSolution), "pressureExact");
    vtkWriter.addField(std::get<1>(analyticalSolution), "velocityExact");
    vtkWriter.addFaceField(std::get<2>(analyticalSolution), "faceVelocityExact");
177 178 179 180 181 182 183
    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
184
    using LinearSolver = Dumux::UMFPackBackend;
185 186 187
    auto linearSolver = std::make_shared<LinearSolver>();

    // the non-linear solver
188 189
    using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
    NewtonSolver nonLinearSolver(assembler, linearSolver);
190 191

    // time loop
192
    const bool printL2Error = getParam<bool>("Problem.PrintL2Error");
193 194 195 196 197
    timeLoop->start(); do
    {
        // set previous solution for storage evaluations
        assembler->setPreviousSolution(xOld);

198 199
        // solve the non-linear system with time step control
        nonLinearSolver.solve(x, *timeLoop);
200 201 202 203 204

        // make the new solution the old solution
        xOld = x;
        gridVariables->advanceTimeStep();

205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222
        if (printL2Error)
        {
            using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
            using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
            using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;

            using L2Error = NavierStokesTestL2Error<Scalar, ModelTraits, PrimaryVariables>;
            const auto l2error = L2Error::calculateL2Error(*problem, x);
            const int numCellCenterDofs = fvGridGeometry->numCellCenterDofs();
            const int numFaceDofs = fvGridGeometry->numFaceDofs();
            std::cout << std::setprecision(8) << "** L2 error (abs/rel) for "
                    << std::setw(6) << numCellCenterDofs << " cc dofs and " << numFaceDofs << " face dofs (total: " << numCellCenterDofs + numFaceDofs << "): "
                    << std::scientific
                    << "L2(p) = " << l2error.first[Indices::pressureIdx] << " / " << l2error.second[Indices::pressureIdx]
                    << ", L2(vx) = " << l2error.first[Indices::velocityXIdx] << " / " << l2error.second[Indices::velocityXIdx]
                    << ", L2(vy) = " << l2error.first[Indices::velocityYIdx] << " / " << l2error.second[Indices::velocityYIdx]
                    << std::endl;
        }
223

224 225
        // advance to the time loop to the next step
        timeLoop->advanceTimeStep();
226
        problem->updateTime(timeLoop->time());
227
        analyticalSolution = createAnalyticalSolution(timeLoop->time(), *problem);
228 229 230 231 232 233 234

        // write vtk output
        vtkWriter.write(timeLoop->time());

        // report statistics of this time step
        timeLoop->reportTimeStep();

235 236
        // set new dt as suggested by newton solver
        timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
237
        problem->updateTimeStepSize(timeLoop->timeStepSize());
238 239 240 241 242 243 244 245 246 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

    } 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;
}