main.cc 9.2 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
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh>
#include <dumux/common/dumuxmessage.hh>
36
#include <dumux/common/parameters.hh>
37
#include <dumux/common/properties.hh>
38

39
40
41
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager_yasp.hh>

42
43
#include <dumux/linear/seqsolverbackend.hh>

44
45
46
47
48
#include <dumux/multidomain/fvassembler.hh>
#include <dumux/multidomain/traits.hh>
#include <dumux/multidomain/newtonsolver.hh>

#include <dumux/freeflow/navierstokes/velocityoutput.hh>
49
#include <test/freeflow/navierstokes/analyticalsolutionvectors.hh>
50
51
52
#include <test/freeflow/navierstokes/errors.hh>

#include "properties.hh"
53

54
int main(int argc, char** argv)
55
56
57
{
    using namespace Dumux;

58
59
60
61
     // define the type tag for this problem
    using MomentumTypeTag = Properties::TTag::AngeliTestMomentum;
    using MassTypeTag = Properties::TTag::AngeliTestMass;

62
63
64
65
66
67
68
69
70

    // 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
71
    Parameters::init(argc, argv);
72
73

    // try to create a grid (from the given grid file or the input file)
74
    GridManager<GetPropType<MomentumTypeTag, Properties::Grid>> gridManager;
75
    gridManager.init();
76
77
78
79
80
81

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

    // we compute on the leaf grid view
82
    const auto& leafGridView = gridManager.grid().leafGridView();
83
84

    // create the finite volume grid geometry
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
    using MomentumGridGeometry = GetPropType<MomentumTypeTag, Properties::GridGeometry>;
    auto momentumGridGeometry = std::make_shared<MomentumGridGeometry>(leafGridView);

    using MassGridGeometry = GetPropType<MassTypeTag, Properties::GridGeometry>;
    auto massGridGeometry = std::make_shared<MassGridGeometry>(leafGridView);

    // the coupling manager
    using Traits = MultiDomainTraits<MomentumTypeTag, MassTypeTag>;
    using CouplingManager = StaggeredFreeFlowCouplingManager<Traits>;

    auto couplingManager = std::make_shared<CouplingManager>();

    // the problem (boundary conditions)
    using MomentumProblem = GetPropType<MomentumTypeTag, Properties::Problem>;
    auto momentumProblem = std::make_shared<MomentumProblem>(momentumGridGeometry, couplingManager);

    using MassProblem = GetPropType<MassTypeTag, Properties::Problem>;
    auto massProblem = std::make_shared<MassProblem>(massGridGeometry, couplingManager);
103
104

    // get some time loop parameters
105
    using Scalar = typename Traits::Scalar;
106
    const auto tStart = getParam<Scalar>("TimeLoop.TStart", 0.0);
107
108
    const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
    const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
109
    auto dt = getParam<Scalar>("TimeLoop.DtInitial");
110
111

    // instantiate time loop
112
    auto timeLoop = std::make_shared<TimeLoop<Scalar>>(tStart, dt, tEnd);
113
114
115
    timeLoop->setMaxTimeStepSize(maxDt);

    // the solution vector
116
117
    constexpr auto momentumIdx = CouplingManager::freeFlowMomentumIndex;
    constexpr auto massIdx = CouplingManager::freeFlowMassIndex;
118
    using SolutionVector = typename Traits::SolutionVector;
119
    SolutionVector x;
120
121
    momentumProblem->applyInitialSolution(x[momentumIdx]);
    massProblem->applyInitialSolution(x[massIdx]);
122
123
124
    auto xOld = x;

    // the grid variables
125
126
127
128
129
130
131
    using MomentumGridVariables = GetPropType<MomentumTypeTag, Properties::GridVariables>;
    auto momentumGridVariables = std::make_shared<MomentumGridVariables>(momentumProblem, momentumGridGeometry);

    using MassGridVariables = GetPropType<MassTypeTag, Properties::GridVariables>;
    auto massGridVariables = std::make_shared<MassGridVariables>(massProblem, massGridGeometry);

    couplingManager->init(momentumProblem, massProblem, std::make_tuple(momentumGridVariables, massGridVariables), x, xOld);
132

133
134
135
136
137
138
    massGridVariables->init(x[massIdx]);
    momentumGridVariables->init(x[momentumIdx]);

    // intialize the vtk output module
    using IOFields = GetPropType<MassTypeTag, Properties::IOFields>;
    VtkOutputModule vtkWriter(*massGridVariables, x[massIdx], massProblem->name());
139
    IOFields::initOutputModule(vtkWriter); // Add model specific output fields
140
141
142
143
144

    vtkWriter.addVelocityOutput(std::make_shared<NavierStokesVelocityOutput<MassGridVariables>>());
    NavierStokesTest::AnalyticalSolutionVectors analyticalSolVectors(momentumProblem, massProblem);
    vtkWriter.addField(analyticalSolVectors.analyticalPressureSolution(), "pressureExact");
    vtkWriter.addField(analyticalSolVectors.analyticalVelocitySolution(), "velocityExact");
145
146
147
    vtkWriter.write(0.0);

    // the assembler with time loop for instationary problem
148
149
150
151
152
153
    using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
    auto assembler = std::make_shared<Assembler>(std::make_tuple(momentumProblem, massProblem),
                                                 std::make_tuple(momentumGridGeometry, massGridGeometry),
                                                 std::make_tuple(momentumGridVariables, massGridVariables),
                                                 couplingManager, timeLoop, xOld);

154
    // the linear solver
155
    using LinearSolver = Dumux::UMFPackBackend;
156
157
158
    auto linearSolver = std::make_shared<LinearSolver>();

    // the non-linear solver
159
160
161
    // the non-linear solver
    using NewtonSolver = Dumux::MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
    NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
162

163
164
    // the discrete L2 and Linfity errors
    const bool printErrors = getParam<bool>("Problem.PrintErrors", false);
165
166
    NavierStokesTest::Errors errors(momentumProblem, massProblem, x);
    NavierStokesTest::ErrorCSVWriter errorCSVWriter(momentumProblem, massProblem);
167

168
169
170
    // time loop
    timeLoop->start(); do
    {
171
172
        massProblem->updateTime(timeLoop->time() + timeLoop->timeStepSize());
        momentumProblem->updateTime(timeLoop->time() + timeLoop->timeStepSize());
173

174
175
        // solve the non-linear system with time step control
        nonLinearSolver.solve(x, *timeLoop);
176
177
178

        // make the new solution the old solution
        xOld = x;
179
180
        momentumGridVariables->advanceTimeStep();
        massGridVariables->advanceTimeStep();
181

182
183
        // print discrete L2 and Linfity errors
        if (printErrors)
184
        {
185
186
            errors.update(x, timeLoop->time() + timeLoop->timeStepSize());
            errorCSVWriter.printErrors(errors);
187
        }
188

189
        // update the analytical solution
190
        analyticalSolVectors.update(timeLoop->time() + timeLoop->timeStepSize());
191

192
193
        // advance to the time loop to the next step
        timeLoop->advanceTimeStep();
194
        analyticalSolVectors.update(timeLoop->time());
195
196
197
198
199
200
201

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

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

202
203
        // set new dt as suggested by newton solver
        timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
204
205
206
207
208

    } while (!timeLoop->finished());

    timeLoop->finalize(leafGridView.comm());

209

210
211
212
213
214
215
216
217
218
219
220
221
222
    ////////////////////////////////////////////////////////////
    // 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