Commit 3ce8271d authored by Leopold Stadler's avatar Leopold Stadler Committed by Timo Koch
Browse files

[shallowwater] added dam break example

parent f68ab65b
......@@ -2,3 +2,4 @@ add_subdirectory(navierstokes)
add_subdirectory(navierstokesnc)
add_subdirectory(rans)
add_subdirectory(ransnc)
add_subdirectory(shallowwater)
dune_symlink_to_source_files(FILES "params.input" references)
dune_add_test(NAME test_shallowwater
SOURCES main.cc
COMPILE_DEFINITIONS TYPETAG=DamBreakWet
COMMAND ${CMAKE_CURRENT_BINARY_DIR}/test_shallowwater)
// -*- 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
* \ingroup ShallowWaterTest
* \brief Test for the shallow water model (wet dam break).
*/
#include <config.h>
#include "problem.hh"
#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 <dumux/io/vtkoutputmodule.hh>
#include <dune/istl/io.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/valgrind.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/common/defaultusagemessage.hh>
#include <dumux/io/grid/gridmanager.hh>
#include <dumux/linear/amgbackend.hh>
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/assembly/fvassembler.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 options 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";
std::cout << errorMessageOut
<< "\n";
}
}
////////////////////////
// the main function
////////////////////////
int main(int argc, char** argv) try
{
using namespace Dumux;
// define the type tag for this problem
using TypeTag = Properties::TTag::TYPETAG;
// 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)
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
using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
fvGridGeometry->update();
// the problem (initial and boundary conditions)
using Problem = GetPropType<TypeTag, Properties::Problem>;
auto problem = std::make_shared<Problem>(fvGridGeometry);
// the solution vector
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
SolutionVector x(fvGridGeometry->numDofs());
problem->applyInitialSolution(x);
auto xOld = x;
// the grid variables
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
gridVariables->init(x);
// get some time loop parameters
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
const auto maxDivisions = getParam<int>("TimeLoop.MaxTimeStepDivisions");
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());
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>;
auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
// the linear solver
using LinearSolver = Dumux::AMGBackend<TypeTag>;
auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
// 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
{
// set previous solution for storage evaluations
assembler->setPreviousSolution(xOld);
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
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;
}
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;
}
[Problem]
Name = damBreak
[TimeLoop]
TEnd = 1.5 # [s]
MaxTimeStepSize = 0.01 # [s]
DtInitial = 0.01 # [s]
[Grid]
Positions0 = 0.0 20.0
Positions1 = 0.0 5.0
Cells0 = 100
Cells1 = 25
[Newton]
EnablePartialReassembly = true
// -*- 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
* \ingroup ShallowWaterTest
* \brief Test for the Shallow water model (wet dam break).
*/
#ifndef DUMUX_DAM_BREAK_TEST_PROBLEM_HH
#define DUMUX_DAM_BREAK_TEST_PROBLEM_HH
#include <dune/grid/yaspgrid.hh>
#include <dumux/discretization/cctpfa.hh>
#include "spatialparams.hh"
#include <dumux/freeflow/shallowwater/model.hh>
#include <dumux/freeflow/shallowwater/problem.hh>
#include <dumux/flux/shallowwater/riemannproblem.hh>
namespace Dumux
{
/*!
* \ingroup ShallowWaterTests
* \brief A simple dam break test for the shallow water equations
*/
template <class TypeTag>
class DamBreakProblem;
// Specify the properties for the problem
namespace Properties{
// Create new type tags
namespace TTag
{
struct ShallowWaterModel{ using InheritsFrom = std::tuple<ShallowWater>; };
struct DamBreakWet{ using InheritsFrom = std::tuple<ShallowWaterModel, CCTpfaModel>; };
} // end namespace TTag
template<class TypeTag>
struct Grid<TypeTag, TTag::ShallowWaterModel>
{ using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
// Set the problem property
template<class TypeTag>
struct Problem<TypeTag, TTag::ShallowWaterModel>{ using type = Dumux::DamBreakProblem<TypeTag>; };
// Set the spatial parameters
template<class TypeTag>
struct SpatialParams<TypeTag, TTag::DamBreakWet>
{
private:
using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
public:
using type = DamBreakSpatialParams<FVGridGeometry, Scalar>;
};
template<class TypeTag>
struct EnableFVGridGeometryCache<TypeTag, TTag::ShallowWaterModel> { static constexpr bool value = true; };
template<class TypeTag>
struct EnableGridVolumeVariablesCache<TypeTag, TTag::ShallowWaterModel> { static constexpr bool value = false; };
template<class TypeTag>
struct EnableGridFluxVariablesCache<TypeTag, TTag::ShallowWaterModel> { static constexpr bool value = false; };
} // end namespace Properties
/*!
* \ingroup Shallow water equations model
* \ingroup ImplicitTestProblems
*
* \brief A simple dam break test
*
* This problem uses the \ref ShallowWaterModel
*
* To run the simulation execute the following line in shell:
* <tt>./test_shallowwater -parameterFile test_shallowwater.input -TimeManager.TEnd 10</tt>
*
* where the initial time step is 0.01 seconds, and the end of the
* simulation time is 10 seconds
*/
template <class TypeTag>
class DamBreakProblem : public ShallowWaterProblem<TypeTag>
{
using ParentType = ShallowWaterProblem<TypeTag>;
using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
using PrimaryVariables = typename VolumeVariables::PrimaryVariables;
using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
using NeumannFluxes = GetPropType<TypeTag, Properties::NumEqVector>;
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using Element = typename GridView::template Codim<0>::Entity;
enum {
// copy some indices for convenience
massBalanceIdx = Indices::massBalanceIdx,
velocityXIdx = Indices::velocityXIdx,
velocityYIdx = Indices::velocityYIdx,
// Grid and world dimension
dimWorld = FVGridGeometry::GridView::dimensionworld,
};
using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
public:
DamBreakProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
{
name_ = getParam<std::string>("Problem.Name");
}
/*!
* \name Problem parameters
*/
// \{
/*!
* \brief The problem name
*
* This is used as a prefix for files generated by the simulation.
*/
const std::string& name() const
{ return name_; }
// \}
/*!
* \name Boundary conditions
*/
// \{
/*!
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary segment.
*
* \param globalPos The position for which the boundary type is set
*/
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
BoundaryTypes bcTypes;
bcTypes.setAllNeumann();
return bcTypes;
}
/*!
* \brief Specifies the neumann bounday
* \param element
* \param fvGeometry
* \param elemVolVars
* \param scvf
*/
NeumannFluxes neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf) const
{
NeumannFluxes values(0.0);
//we need the Riemann invariants to compute the values depending of the boundary type
//since we use a weak imposition we do not have a dirichlet value. We impose fluxes
//based on q,h, etc. computed with the Riemann invariants
const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
const auto& insideVolVars = elemVolVars[insideScv];
const auto& nxy = scvf.unitOuterNormal();
const auto& ip = scvf.ipGlobal();
auto riemannFlux = ShallowWater::riemannProblem(insideVolVars.waterDepth(),
insideVolVars.waterDepth(),
insideVolVars.velocity(0),
-insideVolVars.velocity(0),
insideVolVars.velocity(1),
-insideVolVars.velocity(1),
insideVolVars.bedSurface(),
insideVolVars.bedSurface(),
insideVolVars.gravity(),
nxy);
values[massBalanceIdx] = riemannFlux[0];
values[velocityXIdx] = riemannFlux[1];
values[velocityYIdx] = riemannFlux[2];
return values;
}
/*!
* \brief Evaluate the boundary conditions for a dirichlet boundary
* segment. For the shallow water equations we do a weak
* imposition of boundary conditions. You may use this for
* supercritical flow.
*
* \param globalPos The position for which the Dirichlet value is set
*
* For this method, the \a values parameter stores primary variables.
*/
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables values(0.0);
return values;
}
// \}
/*!
* \name Volume terms
*/
// \{
/*!
* \brief Evaluate the initial values for a control volume.
*
* For this method, the \a values parameter stores primary
* variables.
*
* \param globalPos The position for which the boundary type is set
*/
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables values(0.0);
values[0] = 1.0;
values[1] = 0.0;
values[2] = 0.0;
// water level on the left side of the gate
if (globalPos[0] < 10.0 + eps_)
{
values[0] = 4.0;
}
//water level on the right side of the gate
else
{
values[0] = 1.0;
}
return values;
};
// \}
private:
static constexpr Scalar eps_ = 1.0e-6;
std::string name_;
};
} //end namespace Dumux
#endif
This diff is collapsed.
// -*- 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
* \ingroup ShallowWaterTest
* \brief The spatial parameters for the dam break problem.
*/
#ifndef DUMUX_DAM_BREAK_SPATIAL_PARAMETERS_HH
#define DUMUX_DAM_BREAK_SPATIAL_PARAMETERS_HH
#include <dumux/material/spatialparams/fv.hh>
namespace Dumux{
/*!
* \ingroup ShallowWaterTests
* \brief The spatial parameters class for the dam break test.
*
*/
template<class FVGridGeometry, class Scalar>
class DamBreakSpatialParams
: public FVSpatialParams<FVGridGeometry, Scalar,
DamBreakSpatialParams<FVGridGeometry, Scalar>>
{
using GridView = typename FVGridGeometry::GridView;
using FVElementGeometry = typename FVGridGeometry::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using IndexSet = typename GridView::IndexSet;
using Element = typename GridView::template Codim<0>::Entity;
using ThisType = DamBreakSpatialParams<FVGridGeometry, Scalar>;
using ParentType = FVSpatialParams<FVGridGeometry, Scalar, ThisType>;
enum {
dim=GridView::dimension,
dimWorld=GridView::dimensionworld
};
using GlobalPosition = Dune::FieldVector<Scalar,dimWorld>;
public:
DamBreakSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
{}