main.cc 7.61 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
// -*- 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
21
22
 * \ingroup ShallowWaterTests
 * \brief A test for the shallow water model (wet dam break).
23
24
25
26
27
28
29
 */
#include <config.h>


#include <ctime>
#include <iostream>

30
#include <dune/common/version.hh>
31
32
33
34
35
36
37
38
39
40
41
42
#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 <dumux/io/vtkoutputmodule.hh>
#include <dune/istl/io.hh>

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

#include <dumux/io/grid/gridmanager.hh>
43
44
45

#include <dumux/linear/linearsolvertraits.hh>

46
#if DUNE_VERSION_GT_REV(DUNE_ISTL,2,7,0)
47
48
49
50
#include <dumux/linear/istlsolverfactorybackend.hh>
#else
#include <dumux/linear/amgbackend.hh>
#endif
51

52
53
54
55
#include <dumux/nonlinear/newtonsolver.hh>

#include <dumux/assembly/fvassembler.hh>

56
#include "problem.hh"
57
58
59
60
61
62
63
64
65

////////////////////////
// the main function
////////////////////////
int main(int argc, char** argv) try
{
    using namespace Dumux;

    // define the type tag for this problem
66
    using TypeTag = Properties::TTag::DamBreakWet;
67
68
69
70
71
72
73
74
75

    // 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
76
    Parameters::init(argc, argv);
77
78
79
80
81
82
83
84
85
86
87
88
89

    // try to create a grid (from the given grid file or the input file)
    GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
    gridManager.init();

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

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

    // create the finite volume grid geometry
90
    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
91
92
    auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
    gridGeometry->update();
93
94
95

    // the problem (initial and boundary conditions)
    using Problem = GetPropType<TypeTag, Properties::Problem>;
96
    auto problem = std::make_shared<Problem>(gridGeometry);
97
98
99

    // the solution vector
    using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
100
    SolutionVector x(gridGeometry->numDofs());
101
102
103
104
105
    problem->applyInitialSolution(x);
    auto xOld = x;

    // the grid variables
    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
106
    auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
107
108
109
110
111
112
113
114
115
116
117
118
    gridVariables->init(x);

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

    // intialize the vtk output module
    using IOFields = GetPropType<TypeTag, Properties::IOFields>;

    VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables,x, problem->name());
119
120
121
    vtkWriter.addField(problem->getExactWaterDepth(), "exactWaterDepth");
    vtkWriter.addField(problem->getExactVelocityX(), "exactVelocityX");
    problem->updateAnalyticalSolution(x,*gridVariables,0.0);
122
123
124
125
126
127
128
129
130
    IOFields::initOutputModule(vtkWriter);
    vtkWriter.write(0.0);

    // instantiate time loop
    auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0, dt, tEnd);
    timeLoop->setMaxTimeStepSize(maxDt);

    // the assembler with time loop for instationary problem
    using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
131
    auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, xOld);
132
133

    // the linear solver
134
#if DUNE_VERSION_GT_REV(DUNE_ISTL,2,7,0)
135
    using LinearSolver = IstlSolverFactoryBackend<LinearSolverTraits<GridGeometry>>;
136
137
138
#else
    using LinearSolver = AMGBiCGSTABBackend<LinearSolverTraits<GridGeometry>>;
#endif
139
    auto linearSolver = std::make_shared<LinearSolver>(leafGridView, gridGeometry->dofMapper());
140
141
142
143
144
145
146
147
148
149
150
151
152

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

    //! set some check point at the end of the time loop
    timeLoop->setCheckPoint(tEnd);

    // time loop
    timeLoop->start(); do
    {
        nonLinearSolver.solve(x,*timeLoop);

153
154
155
        // update the analytical solution
        problem->updateAnalyticalSolution(x,*gridVariables,timeLoop->time()+timeLoop->timeStepSize());

156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
        // 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
        if (timeLoop->isCheckPoint())
            vtkWriter.write(timeLoop->time());

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

        // set new dt as suggested by newton controller
        timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));


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

192
catch (const Dumux::ParameterException &e)
193
194
195
196
{
    std::cerr << std::endl << e << " ---> Abort!" << std::endl;
    return 1;
}
197
catch (const Dune::DGFException & e)
198
199
200
201
202
203
204
205
{
    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;
}
206
catch (const Dune::Exception &e)
207
208
209
210
211
212
213
214
215
{
    std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
    return 3;
}
catch (...)
{
    std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
    return 4;
}