Commit bf041ac6 authored by Timo Koch's avatar Timo Koch
Browse files

Merge branch 'fix/upstream-sequential-remove' into 'master'

Fix/upstream sequential remove

See merge request !164
parents ca80243e e83ed623
Pipeline #8253 passed with stage
// -*- 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/>. *
*****************************************************************************/
#ifndef DUMUX_FVVELOCITY_HH
#define DUMUX_FVVELOCITY_HH
// dumux environment
#include <dumux/common/math.hh>
#include <dumux/porousmediumflow/sequential/pressureproperties.hh>
#include "velocitydefault.hh"
/**
* @file
* @brief Finite volume velocity reconstruction
*/
namespace Dumux
{
/*! \ingroup IMPET
*
* \brief Base class for finite volume velocity reconstruction
*
* Provides a basic frame for calculating a global velocity field.
* The definition of the local velocity calculation as well as the storage or other postprocessing
* has to be provided by the local velocity implementation.
* This local implementation has to have the form of VelocityDefault.
*
* \tparam TypeTag The Type Tag
* \tparam Velocity The implementation of the local velocity calculation
*/
template<class TypeTag, class Velocity> class FVVelocity
{
using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
using Problem = GetPropType<TypeTag, Properties::Problem>;
using CellData = GetPropType<TypeTag, Properties::CellData>;
public:
//!Initialize velocity implementation
void initialize()
{
velocity_.initialize();
}
//function which iterates through the grid and calculates the global velocity field
void calculateVelocity();
/*! \brief Adds velocity output to the output file
*
* \tparam MultiWriter Class defining the output writer
* \param writer The output writer (usually a <tt>VTKMultiWriter</tt> object)
*
*/
template<class MultiWriter>
void addOutputVtkFields(MultiWriter &writer)
{
velocity_.addOutputVtkFields(writer);
}
//! Constructs a FVVelocity object
/**
* \param problem A problem class object
*/
FVVelocity(Problem& problem) :
problem_(problem), velocity_(problem)
{}
private:
Problem& problem_;
Velocity velocity_;
};
/*! \brief Function which reconstructs a global velocity field
*
* Iterates through the grid and calls the local calculateVelocity(...) or calculateVelocityOnBoundary(...)
* functions which have to be provided by the local velocity implementation (see e.g. VelocityDefault )
*/
template<class TypeTag, class Velocity>
void FVVelocity<TypeTag, Velocity>::calculateVelocity()
{
for (const auto& element : elements(problem_.gridView()))
{
// cell information
int globalIdxI = problem_.variables().index(element);
CellData& cellDataI = problem_.variables().cellData(globalIdxI);
/***** flux term ***********/
// iterate over all faces of the cell
for (const auto& intersection : intersections(problem_.gridView(), element))
{
/************* handle interior face *****************/
if (intersection.neighbor())
{
int isIndex = intersection.indexInInside();
if (!cellDataI.fluxData().haveVelocity(isIndex))
velocity_.calculateVelocity(intersection, cellDataI);
} // end neighbor
/************* boundary face ************************/
else
{
velocity_.calculateVelocityOnBoundary(intersection, cellDataI);
}
} //end interfaces loop
} // end grid traversal
return;
}
}//end namespace Dumux
#endif
// -*- 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/>. *
*****************************************************************************/
#ifndef DUMUX_ONE_MODEL_PROBLEM_HH
#define DUMUX_ONE_MODEL_PROBLEM_HH
#include <dune/common/shared_ptr.hh>
#include <dumux/common/properties.hh>
#include <dumux/porousmediumflow/sequential/properties.hh>
#include <dumux/io/vtkmultiwriter.hh>
#include <dumux/io/restart.hh>
/**
* @file
* @brief Base class for definition of an sequential diffusion (pressure) or transport problem
*/
namespace Dumux
{
/*! \ingroup IMPET
*
* @brief Base class for definition of an sequential diffusion (pressure) or transport problem
*
* @tparam TypeTag The Type Tag
*/
template<class TypeTag>
class OneModelProblem
{
private:
using Implementation = GetPropType<TypeTag, Properties::Problem>;
using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
using Grid = typename GridView::Grid;
using TimeManager = GetPropType<TypeTag, Properties::TimeManager>;
using VtkMultiWriter = Dumux::VtkMultiWriter<GridView>;
using Variables = GetPropType<TypeTag, Properties::Variables>;
using Model = GetPropType<TypeTag, Properties::Model>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
using VertexMapper = typename SolutionTypes::VertexMapper;
using ElementMapper = typename SolutionTypes::ElementMapper;
enum
{
dim = GridView::dimension,
dimWorld = GridView::dimensionworld
};
enum
{
wetting = 0, nonwetting = 1
};
using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
using Element = typename GridView::Traits::template Codim<0>::Entity;
using Intersection = typename GridView::Intersection;
using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
using BoundaryTypes = GetPropType<TypeTag, Properties::SequentialBoundaryTypes>;
// private!! copy constructor
OneModelProblem(const OneModelProblem&)
{}
public:
//! Constructs an object of type OneModelProblemProblem
/*!
* \tparam TypeTag The TypeTag
* \tparam verbose Output level for TimeManager
*/
OneModelProblem(Grid& grid, bool verbose = true)
: gridView_(grid.leafGridView()),
bBoxMin_(std::numeric_limits<double>::max()),
bBoxMax_(-std::numeric_limits<double>::max()),
variables_(grid.leafGridView()),
outputInterval_(1),
outputTimeInterval_(0)
{
// calculate the bounding box of the grid view
using std::max;
using std::min;
for (const auto& vertex : vertices(grid.leafGridView())) {
for (int i=0; i<dim; i++) {
bBoxMin_[i] = min(bBoxMin_[i], vertex.geometry().center()[i]);
bBoxMax_[i] = max(bBoxMax_[i], vertex.geometry().center()[i]);
}
}
timeManager_ = std::make_shared<TimeManager>(verbose);
model_ = std::make_shared<Model>(asImp_()) ;
maxTimeStepSize_ = getParam<Scalar>("TimeManager.MaxTimeStepSize", std::numeric_limits<Scalar>::max());
}
//! Constructs an object of type OneModelProblemProblem
/*!
* \tparam TypeTag The TypeTag
* \tparam verbose Output level for TimeManager
*/
OneModelProblem(TimeManager& timeManager, Grid& grid)
: gridView_(grid.leafGridView()),
bBoxMin_(std::numeric_limits<double>::max()),
bBoxMax_(-std::numeric_limits<double>::max()),
variables_(grid.leafGridView()),
outputInterval_(1),
outputTimeInterval_(0)
{
// calculate the bounding box of the grid view
using std::max;
using std::min;
for (const auto& vertex : vertices(grid.leafGridView())) {
for (int i=0; i<dim; i++) {
bBoxMin_[i] = min(bBoxMin_[i], vertex.geometry().center()[i]);
bBoxMax_[i] = max(bBoxMax_[i], vertex.geometry().center()[i]);
}
}
timeManager_ = Dune::stackobject_to_shared_ptr<TimeManager>(timeManager);
model_ = std::make_shared<Model>(asImp_()) ;
maxTimeStepSize_ = getParam<Scalar>("TimeManager.MaxTimeStepSize", std::numeric_limits<Scalar>::max());
}
/*!
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary segment.
*
* \param bcTypes The boundary types for the conservation equations
* \param intersection The intersection for which the boundary type is set
*/
void boundaryTypes(BoundaryTypes &bcTypes,
const Intersection &intersection) const
{
// forward it to the method which only takes the global coordinate
asImp_().boundaryTypesAtPos(bcTypes, intersection.geometry().center());
}
/*!
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary segment.
*
* \param bcTypes The boundary types for the conservation equations
* \param globalPos The position of the center of the boundary intersection
*/
void boundaryTypesAtPos(BoundaryTypes &bcTypes,
const GlobalPosition &globalPos) const
{
// Throw an exception (there is no reasonable default value
// for Dirichlet conditions)
DUNE_THROW(Dune::InvalidStateException,
"The problem does not provide "
"a boundaryTypesAtPos() method.");
}
/*!
* \brief Evaluate the boundary conditions for a dirichlet
* control volume.
*
* \param values The dirichlet values for the primary variables
* \param intersection The boundary intersection
*
* For this method, the \a values parameter stores primary variables.
*/
void dirichlet(PrimaryVariables &values,
const Intersection &intersection) const
{
// forward it to the method which only takes the global coordinate
asImp_().dirichletAtPos(values, intersection.geometry().center());
}
/*!
* \brief Evaluate the boundary conditions for a dirichlet
* control volume.
*
* \param values The dirichlet values for the primary variables
* \param globalPos The position of the center of the boundary intersection
*
* For this method, the \a values parameter stores primary variables.
*/
void dirichletAtPos(PrimaryVariables &values,
const GlobalPosition &globalPos) const
{
// Throw an exception (there is no reasonable default value
// for Dirichlet conditions)
DUNE_THROW(Dune::InvalidStateException,
"The problem specifies that some boundary "
"segments are dirichlet, but does not provide "
"a dirichletAtPos() method.");
}
/*!
* \brief Evaluate the boundary conditions for a neumann
* boundary segment.
*
* \param values The neumann values for the conservation equations [kg / (m^2 *s )]
* \param intersection The boundary intersection
*
* For this method, the \a values parameter stores the mass flux
* in normal direction of each phase. Negative values mean influx.
*/
void neumann(PrimaryVariables &values,
const Intersection &intersection) const
{
// forward it to the interface with only the global position
asImp_().neumannAtPos(values, intersection.geometry().center());
}
/*!
* \brief Evaluate the boundary conditions for a neumann
* boundary segment.
*
* \param values The neumann values for the conservation equations [kg / (m^2 *s )]
* \param globalPos The position of the center of the boundary intersection
*
* For this method, the \a values parameter stores the mass flux
* in normal direction of each phase. Negative values mean influx.
*/
void neumannAtPos(PrimaryVariables &values,
const GlobalPosition &globalPos) const
{
// Throw an exception (there is no reasonable default value
// for Neumann conditions)
DUNE_THROW(Dune::InvalidStateException,
"The problem specifies that some boundary "
"segments are neumann, but does not provide "
"a neumannAtPos() method.");
}
/*!
* \brief Evaluate the source term
*
* \param values The source and sink values for the conservation equations
* \param element The element
*
* For this method, the \a values parameter stores the rate mass
* generated or annihilate per volume unit. Positive values mean
* that mass is created, negative ones mean that it vanishes.
*/
void source(PrimaryVariables &values,
const Element &element) const
{
// forward to generic interface
asImp_().sourceAtPos(values, element.geometry().center());
}
/*!
* \brief Evaluate the source term for all phases within a given
* sub-control-volume.
*
* \param values The source and sink values for the conservation equations
* \param globalPos The position of the center of the finite volume
* for which the source term ought to be
* specified in global coordinates
*
* For this method, the \a values parameter stores the rate mass
* generated or annihilate per volume unit. Positive values mean
* that mass is created, negative ones mean that it vanishes.
*/
void sourceAtPos(PrimaryVariables &values,
const GlobalPosition &globalPos) const
{ // Throw an exception (there is no initial condition)
DUNE_THROW(Dune::InvalidStateException,
"The problem does not provide "
"a sourceAtPos() method.");
}
/*!
* \brief Evaluate the initial value for a control volume.
*
* \param values The initial values for the primary variables
* \param element The element
*
* For this method, the \a values parameter stores primary
* variables.
*/
void initial(PrimaryVariables &values,
const Element &element) const
{
// forward to generic interface
asImp_().initialAtPos(values, element.geometry().center());
}
/*!
* \brief Evaluate the initial value for a control volume.
*
* \param values The dirichlet values for the primary variables
* \param globalPos The position of the center of the finite volume
* for which the initial values ought to be
* set (in global coordinates)
*
* For this method, the \a values parameter stores primary variables.
*/
void initialAtPos(PrimaryVariables &values,
const GlobalPosition &globalPos) const
{
// initialize with 0 by default
values = 0;
}
/*!
* \brief Called by the TimeManager in order to
* initialize the problem.
*/
void init()
{
// set the initial condition of the model
variables_.initialize();
model().initialize();
}
/*!
* \brief Called by TimeManager just before the time
* integration.
*/
void preTimeStep()
{}
/*!
* \brief Called by TimeManager in order to do a time
* integration on the model.
*/
void timeIntegration()
{}
/*!
* \brief Called by TimeManager whenever a solution for a
* timestep has been computed and the simulation time has
* been updated.
*
* This is used to do some janitorial tasks like writing the
* current solution to disk.
*/
void postTimeStep()
{}
/*!
* \brief Called by the time manager after everything which can be
* done about the current time step is finished and the
* model should be prepared to do the next time integration.
*/
void advanceTimeLevel()
{}
/*!
* \brief Returns the user specified maximum time step size
*
* Overload in problem for custom needs.
*/
Scalar maxTimeStepSize() const
{ return maxTimeStepSize_; }
/*!
* \brief Returns the current time step size [seconds].
*/
Scalar timeStepSize() const
{ return timeManager().timeStepSize(); }
/*!
* \brief Sets the current time step size [seconds].
*/
void setTimeStepSize(Scalar dt)
{ timeManager().setTimeStepSize(dt); }
/*!
* \brief Called by TimeManager whenever a solution for a
* timestep has been computed and the simulation time has
* been updated.
*/
Scalar nextTimeStepSize(Scalar dt)
{ return timeManager().timeStepSize();}
/*!
* \brief Returns true if a restart file should be written to
* disk.
*
* The default behaviour is to write one restart file every 5 time
* steps. This file is intented to be overwritten by the
* implementation.
*/
bool shouldWriteRestartFile() const
{
return
timeManager().timeStepIndex() > 0 &&
(timeManager().timeStepIndex() % 5 == 0);
}
/*!
* \brief Sets a time interval for Output
*
* The default is 0.0 -> Output determined by output number interval (<tt>setOutputInterval(int)</tt>)
*/
void setOutputTimeInterval(const Scalar timeInterval)
{
outputTimeInterval_ = (timeInterval > 0.0) ? timeInterval : 1e100;
timeManager().startNextEpisode(outputTimeInterval_);
}
/*!
* \brief Sets the interval for Output
*
* The default is 1 -> Output every time step
*/
void setOutputInterval(int interval)
{ outputInterval_ = interval; }
/*!
* \brief Returns true if the current solution should be written to
* disk (i.e. as a VTK file)
*
* The default behaviour is to write out every the solution for
* very time step. This file is intented to be overwritten by the
* implementation.
*/
bool shouldWriteOutput() const
{
if (outputInterval_ > 0)
{
if (timeManager().timeStepIndex() % outputInterval_ == 0
|| timeManager().willBeFinished()
|| timeManager().episodeWillBeFinished())
{
return true;
}
}
else if (timeManager().willBeFinished()
|| timeManager().episodeWillBeFinished() || timeManager().timeStepIndex() == 0)
{
return true;
}
return false;
}
void addOutputVtkFields()
{}
//! Write the fields current solution into an VTK output file.
void writeOutput(bool verbose = true)
{
if (verbose && gridView().comm().rank() == 0)
std::cout << "Writing result file for current time step\n";
if (!resultWriter_)
resultWriter_ = std::make_shared<VtkMultiWriter>(gridView(), asImp_().name());
resultWriter_->beginWrite(timeManager().time() + timeManager().timeStepSize());
model().addOutputVtkFields(*resultWriter_);
asImp_().addOutputVtkFields();
resultWriter_->endWrite();
}
/*!
* \brief Called when the end of an simulation episode is reached.