main_momentum.cc 7.41 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
// -*- 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 3 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
 * \ingroup NavierStokesTests
 * \brief Test for the staggered grid Navier-Stokes model (Donea 2003, \cite Donea2003).
 */

#include <config.h>

#include <iostream>
#include <dune/common/parallel/mpihelper.hh>

#include <dumux/common/dumuxmessage.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/properties.hh>

#include <dumux/assembly/fvassembler.hh>

#include <dumux/io/grid/gridmanager.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/vtk/function.hh>
#include <dumux/io/vtk/intersectionwriter.hh>

#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/nonlinear/newtonsolver.hh>

#include <dumux/linear/linearsolvertraits.hh>
45
#include <dumux/linear/algebratraits.hh>
46
47
48
49
#include <dumux/linear/istlsolverfactorybackend.hh>

#include "properties_momentum.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
template<class GridGeometry, class GridVariables, class SolutionVector>
void updateVelocities(
    std::vector<Dune::FieldVector<double, 2>>& velocity,
    std::vector<Dune::FieldVector<double, 2>>& faceVelocity,
    const GridGeometry& gridGeometry,
    const GridVariables& gridVariables,
    const SolutionVector& x
){
    auto fvGeometry = localView(gridGeometry);
    auto elemVolVars = localView(gridVariables.curGridVolVars());
    for (const auto& element : elements(gridGeometry.gridView()))
    {
        fvGeometry.bind(element);
        elemVolVars.bind(element, fvGeometry, x);

        for (const auto& scv : scvs(fvGeometry))
        {
            const auto& vars = elemVolVars[scv];
            velocity[scv.elementIndex()][scv.dofAxis()] += 0.5*vars.velocity();
            faceVelocity[scv.dofIndex()][scv.dofAxis()] = vars.velocity();
        }
    }
}

template<class GridGeometry>
void updateRank(
    std::vector<int>& rank,
    const GridGeometry& gridGeometry
){
    for (const auto& element : elements(gridGeometry.gridView(), Dune::Partitions::interior))
    {
        const auto eIdxGlobal = gridGeometry.elementMapper().index(element);
        rank[eIdxGlobal] = gridGeometry.gridView().comm().rank();
    }
}

86
87
88
89
90
91
92
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
int main(int argc, char** argv)
{
    using namespace Dumux;

    using TypeTag = Properties::TTag::DoneaTestMomentum;

    const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);

    if (mpiHelper.rank() == 0)
        DumuxMessage::print(/*firstCall=*/true);

    Parameters::init(argc, argv);

    using GridManager = Dumux::GridManager<GetPropType<TypeTag, Properties::Grid>>;
    GridManager gridManager;
    gridManager.init();

    ////////////////////////////////////////////////////////////
    // solve Stokes problem on this grid
    ////////////////////////////////////////////////////////////

    const auto& leafGridView = gridManager.grid().leafGridView();

    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
    auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);

    using Problem = GetPropType<TypeTag, Properties::Problem>;
    auto problem = std::make_shared<Problem>(gridGeometry);

    using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
    SolutionVector x(gridGeometry->numDofs());

    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
    auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
    gridVariables->init(x);

122
123
124
    std::vector<Dune::FieldVector<double, 2>> velocity(gridGeometry->gridView().size(0));
    std::vector<Dune::FieldVector<double, 2>> faceVelocity(x.size());
    updateVelocities(velocity, faceVelocity, *gridGeometry, *gridVariables, x);
125

126
127
    std::vector<int> rank(gridGeometry->gridView().size(0));
    updateRank(rank, *gridGeometry);
128

129
130
    std::vector<std::size_t> dofIdx(x.size());
    for (const auto& facet : facets(gridGeometry->gridView()))
131
    {
132
133
        const auto idx = gridGeometry->gridView().indexSet().index(facet);
        dofIdx[idx] = idx;
134
135
    }

136
137
    Dune::VTKWriter<typename GridGeometry::GridView> writer(gridGeometry->gridView());
    using Field = Vtk::template Field<typename GridGeometry::GridView>;
138
139
140
141
142
    writer.addCellData(Field(
        gridGeometry->gridView(), gridGeometry->elementMapper(), velocity,
        "velocity", /*numComp*/2, /*codim*/0
    ).get());
    writer.addCellData(rank, "rank");
143
    writer.write("donea_momentum_0");
144
145
146

    ConformingIntersectionWriter faceVtk(gridGeometry->gridView());
    faceVtk.addField(dofIdx, "dofIdx");
147
148
    faceVtk.addField(faceVelocity, "velocityVector");
    faceVtk.addField([&](const auto& is, const auto idx) {
149
150
        const auto& facet = is.inside().template subEntity <1> (is.indexInInside());
        return facet.partitionType();
151
152
153
154
155
156
    }, "partitionType");
    faceVtk.write("donea_momentum_face_0_" + std::to_string(gridGeometry->gridView().comm().rank()), Dune::VTK::ascii);

    using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
    auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables);

157
158
    using LinearSolver = IstlSolverFactoryBackend<LinearSolverTraits<GridGeometry>,
                                                  LinearAlgebraTraitsFromAssembler<Assembler>>;
159
160
    const auto dofMapper = LinearSolverTraits<GridGeometry>::DofMapper(gridGeometry->gridView());
    auto linearSolver = std::make_shared<LinearSolver>(gridGeometry->gridView(), dofMapper);
161

162
163
164
165
166
167
168
169
170
171
172
    using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
    NewtonSolver nonLinearSolver(assembler, linearSolver);

    nonLinearSolver.solve(x);

    ////////////////////////////////////////////////////////////
    // write VTK output
    ////////////////////////////////////////////////////////////
    updateVelocities(velocity, faceVelocity, *gridGeometry, *gridVariables, x);
    writer.write("donea_momentum_1");
    faceVtk.write("donea_momentum_face_1_" + std::to_string(gridGeometry->gridView().comm().rank()), Dune::VTK::ascii);
173
174
175
176
177
178
179
180
181
182
183
184
185
186

    ////////////////////////////////////////////////////////////
    // finalize, print parameters and Dumux message to say goodbye
    ////////////////////////////////////////////////////////////

    // print dumux end message
    if (mpiHelper.rank() == 0)
    {
        Parameters::print();
        DumuxMessage::print(/*firstCall=*/false);
    }

    return 0;
}