OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #821 from akva2/move_dunefem_adaptation_to_separate_class

Browse Source 
move the dynamic refinement part related to dune fem to separate class
This commit is contained in:
Atgeirr Flø Rasmussen 2023-08-11 09:31:19 +02:00 committed by GitHub
commit 67baba967e
6 changed files with 292 additions and 129 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

191
opm/models/discretization/common/fvbasediscretization.hh
View File

@ -39,7 +39,6 @@
#include "fvbaseboundarycontext.hh"
#include "fvbaseconstraintscontext.hh"
#include "fvbaseconstraints.hh"
#include "fvbasediscretization.hh"
#include "fvbasegradientcalculator.hh"
#include "fvbasenewtonmethod.hh"
#include "fvbaseprimaryvariables.hh"
@ -64,12 +63,6 @@
#include <dune/common/fmatrix.hh>
#include <dune/istl/bvector.hh>
#if HAVE_DUNE_FEM
#include <dune/fem/space/common/adaptationmanager.hh>
#include <dune/fem/space/common/restrictprolongtuple.hh>
#include <dune/fem/function/blockvectorfunction.hh>
#include <dune/fem/misc/capabilities.hh>
#endif
#include <algorithm>
#include <cstddef>
@ -78,10 +71,13 @@
#include <stdexcept>
#include <sstream>
#include <string>
#include <type_traits>
#include <vector>
namespace Opm {
template<class TypeTag>
class FvBaseDiscretizationNoAdapt;
template<class TypeTag>
class FvBaseDiscretization;
} // namespace Opm
@ -322,6 +318,23 @@ struct UseVolumetricResidual<TypeTag, TTag::FvBaseDiscretization> { static const
template<class TypeTag>
struct EnableExperiments<TypeTag, TTag::FvBaseDiscretization> { static constexpr bool value = true; };
template <class TypeTag, class MyTypeTag>
struct BaseDiscretizationType {
using type = UndefinedProperty;
};
#if !HAVE_DUNE_FEM
template<class TypeTag>
struct BaseDiscretizationType<TypeTag,TTag::FvBaseDiscretization> {
using type = FvBaseDiscretizationNoAdapt<TypeTag>;
};
template<class TypeTag>
struct DiscreteFunction<TypeTag, TTag::FvBaseDiscretization> {
using BaseDiscretization = FvBaseDiscretization<TypeTag>;
using type = typename BaseDiscretization::BlockVectorWrapper;
};
#endif
} // namespace Opm::Properties
namespace Opm {
@ -422,25 +435,8 @@ public:
};
private:
#if HAVE_DUNE_FEM
using DiscreteFunctionSpace = GetPropType<TypeTag, Properties::DiscreteFunctionSpace> ;
// discrete function storing solution data
using DiscreteFunction = Dune::Fem::ISTLBlockVectorDiscreteFunction<DiscreteFunctionSpace, PrimaryVariables>;
// problem restriction and prolongation operator for adaptation
using Problem = GetPropType<TypeTag, Properties::Problem> ;
using ProblemRestrictProlongOperator = typename Problem :: RestrictProlongOperator ;
// discrete function restriction and prolongation operator for adaptation
using DiscreteFunctionRestrictProlong = Dune::Fem::RestrictProlongDefault< DiscreteFunction >;
using RestrictProlong = Dune::Fem::RestrictProlongTuple< DiscreteFunctionRestrictProlong, ProblemRestrictProlongOperator >;
// adaptation classes
using AdaptationManager = Dune::Fem::AdaptationManager<Grid, RestrictProlong >;
#else
using DiscreteFunction = BlockVectorWrapper ;
using DiscreteFunctionSpace = size_t ;
#endif
using DiscreteFunctionSpace = GetPropType<TypeTag, Properties::DiscreteFunctionSpace>;
using DiscreteFunction = GetPropType<TypeTag, Properties::DiscreteFunction>;
// copying a discretization object is not a good idea
FvBaseDiscretization(const FvBaseDiscretization& );
@ -457,25 +453,11 @@ public:
, newtonMethod_(simulator)
, localLinearizer_(ThreadManager::maxThreads())
, linearizer_(new Linearizer())
#if HAVE_DUNE_FEM
, space_( simulator.vanguard().gridPart() )
#else
, space_( asImp_().numGridDof() )
#endif
, enableGridAdaptation_( EWOMS_GET_PARAM(TypeTag, bool, EnableGridAdaptation) )
, enableIntensiveQuantityCache_(EWOMS_GET_PARAM(TypeTag, bool, EnableIntensiveQuantityCache))
, enableStorageCache_(EWOMS_GET_PARAM(TypeTag, bool, EnableStorageCache))
, enableThermodynamicHints_(EWOMS_GET_PARAM(TypeTag, bool, EnableThermodynamicHints))
{
#if HAVE_DUNE_FEM
if (enableGridAdaptation_ && !Dune::Fem::Capabilities::isLocallyAdaptive<Grid>::v)
throw std::invalid_argument("Grid adaptation enabled, but chosen Grid is not capable"
" of adaptivity");
#else
if (enableGridAdaptation_)
throw std::invalid_argument("Grid adaptation currently requires the presence of the "
"dune-fem module");
#endif
bool isEcfv = std::is_same<Discretization, EcfvDiscretization<TypeTag> >::value;
if (enableGridAdaptation_ && !isEcfv)
throw std::invalid_argument("Grid adaptation currently only works for the "
@ -486,8 +468,6 @@ public:
size_t numDof = asImp_().numGridDof();
for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
solution_[timeIdx].reset(new DiscreteFunction("solution", space_));
if (storeIntensiveQuantities()) {
intensiveQuantityCache_[timeIdx].resize(numDof);
intensiveQuantityCacheUpToDate_[timeIdx].resize(numDof, /*value=*/false);
@ -1420,48 +1400,8 @@ public:
*/
void adaptGrid()
{
#if HAVE_DUNE_FEM
// adapt the grid if enabled and if all dependencies are available
// adaptation is only done if markForGridAdaptation returns true
if (enableGridAdaptation_)
{
// check if problem allows for adaptation and cells were marked
if( simulator_.problem().markForGridAdaptation() )
{
// adapt the grid and load balance if necessary
adaptationManager().adapt();
// if the grid has potentially changed, we need to re-create the
// supporting data structures.
#if DUNE_VERSION_NEWER(DUNE_GRID, 2, 8)
elementMapper_.update(gridView_);
vertexMapper_.update(gridView_);
#else
elementMapper_.update();
vertexMapper_.update();
#endif
resetLinearizer();
// this is a bit hacky because it supposes that Problem::finishInit()
// works fine multiple times in a row.
//
// TODO: move this to Problem::gridChanged()
finishInit();
// notify the problem that the grid has changed
//
// TODO: come up with a mechanism to access the unadapted data structures
// outside of the problem (i.e., grid, mappers, solutions)
simulator_.problem().gridChanged();
// notify the modules for visualization output
auto outIt = outputModules_.begin();
auto outEndIt = outputModules_.end();
for (; outIt != outEndIt; ++outIt)
(*outIt)->allocBuffers();
}
}
#endif
throw std::invalid_argument("Grid adaptation need to be implemented for "
"specific settings of grid and function spaces");
}
/*!
@ -1498,7 +1438,9 @@ public:
void advanceTimeLevel()
{
// at this point we can adapt the grid
asImp_().adaptGrid();
if (this->enableGridAdaptation_) {
asImp_().adaptGrid();
}
// make the current solution the previous one.
solution(/*timeIdx=*/1) = solution(/*timeIdx=*/0);
@ -1907,22 +1849,6 @@ public:
bool storeIntensiveQuantities() const
{ return enableIntensiveQuantityCache_ || enableThermodynamicHints_; }
#if HAVE_DUNE_FEM
AdaptationManager& adaptationManager()
{
if( ! adaptationManager_ )
{
// create adaptation objects here, because when doing so in constructor
// problem is not yet intialized, aka seg fault
restrictProlong_.reset(
new RestrictProlong( DiscreteFunctionRestrictProlong(*(solution_[/*timeIdx=*/ 0] )),
simulator_.problem().restrictProlongOperator() ) );
adaptationManager_.reset( new AdaptationManager( simulator_.vanguard().grid(), *restrictProlong_ ) );
}
return *adaptationManager_;
}
#endif
const Timer& prePostProcessTimer() const
{ return prePostProcessTimer_; }
@ -1938,13 +1864,9 @@ public:
template<class Serializer>
void serializeOp(Serializer& serializer)
{
for (auto& sol : solution_) {
#if HAVE_DUNE_FEM
serializer(sol->blockVector());
#else
serializer(*sol);
#endif
}
using BaseDiscretization = GetPropType<TypeTag, Properties::BaseDiscretizationType>;
using Helper = typename BaseDiscretization::template SerializeHelper<Serializer>;
Helper::serializeOp(serializer, solution_);
}
bool operator==(const FvBaseDiscretization& rhs) const
@ -2049,15 +1971,8 @@ protected:
// while these are logically bools, concurrent writes to vector<bool> are not thread safe.
mutable std::vector<unsigned char> intensiveQuantityCacheUpToDate_[historySize];
DiscreteFunctionSpace space_;
mutable std::array< std::unique_ptr< DiscreteFunction >, historySize > solution_;
#if HAVE_DUNE_FEM
std::unique_ptr<RestrictProlong> restrictProlong_;
std::unique_ptr<AdaptationManager> adaptationManager_;
#endif
std::list<BaseOutputModule<TypeTag>*> outputModules_;
Scalar gridTotalVolume_;
@ -2071,6 +1986,50 @@ protected:
bool enableStorageCache_;
bool enableThermodynamicHints_;
};
/*!
* \ingroup FiniteVolumeDiscretizations
*
* \brief The base class for the finite volume discretization schemes without adaptation.
*/
template<class TypeTag>
class FvBaseDiscretizationNoAdapt : public FvBaseDiscretization<TypeTag>
{
using ParentType = FvBaseDiscretization<TypeTag>;
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using DiscreteFunction = GetPropType<TypeTag, Properties::DiscreteFunction>;
static constexpr unsigned historySize = getPropValue<TypeTag, Properties::TimeDiscHistorySize>();
public:
template<class Serializer>
struct SerializeHelper {
template<class SolutionType>
static void serializeOp(Serializer& serializer,
SolutionType& solution)
{
for (auto& sol : solution) {
serializer(*sol);
}
}
};
FvBaseDiscretizationNoAdapt(Simulator& simulator)
: ParentType(simulator)
{
if (this->enableGridAdaptation_) {
throw std::invalid_argument("Grid adaptation need to use"
" BaseDiscretization = FvBaseDiscretizationFemAdapt"
" which currently requires the presence of the"
" dune-fem module");
}
size_t numDof = this->asImp_().numGridDof();
for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
this->solution_[timeIdx] = std::make_unique<DiscreteFunction>("solution", numDof);
}
}
};
} // namespace Opm
#endif
#endif // EWOMS_FV_BASE_DISCRETIZATION_HH

174
opm/models/discretization/common/fvbasediscretizationfemadapt.hh Normal file
View File

@ -0,0 +1,174 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \copydoc Opm::FvBaseDiscretization
*/
#ifndef EWOMS_FV_BASE_DISCRETIZATION_FEMADAPT_HH
#define EWOMS_FV_BASE_DISCRETIZATION_FEMADAPT_HH
#include <opm/models/discretization/common/fvbasediscretization.hh>
#include <dune/fem/space/common/adaptationmanager.hh>
#include <dune/fem/space/common/restrictprolongtuple.hh>
#include <dune/fem/function/blockvectorfunction.hh>
#include <dune/fem/misc/capabilities.hh>
namespace Opm {
template<class TypeTag>
class FvBaseDiscretizationFemAdapt;
namespace Properties {
template<class TypeTag>
struct BaseDiscretizationType<TypeTag,TTag::FvBaseDiscretization> {
using type = FvBaseDiscretizationFemAdapt<TypeTag>;
};
template<class TypeTag>
struct DiscreteFunction<TypeTag, TTag::FvBaseDiscretization> {
using DiscreteFunctionSpace = GetPropType<TypeTag, Properties::DiscreteFunctionSpace>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using type = Dune::Fem::ISTLBlockVectorDiscreteFunction<DiscreteFunctionSpace, PrimaryVariables>;
};
} // namespace Properties
/*!
* \ingroup FiniteVolumeDiscretizations
*
* \brief The base class for the finite volume discretization schemes.
*/
template <class TypeTag>
class FvBaseDiscretizationFemAdapt : public FvBaseDiscretization<TypeTag>
{
using Grid = GetPropType<TypeTag, Properties::Grid>;
using ParentType = FvBaseDiscretization<TypeTag>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using Problem = GetPropType<TypeTag, Properties::Problem>;
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
static constexpr unsigned historySize = getPropValue<TypeTag, Properties::TimeDiscHistorySize>();
using DiscreteFunctionSpace = GetPropType<TypeTag, Properties::DiscreteFunctionSpace>;
// discrete function storing solution data
using DiscreteFunction = Dune::Fem::ISTLBlockVectorDiscreteFunction<DiscreteFunctionSpace, PrimaryVariables>;
// problem restriction and prolongation operator for adaptation
using ProblemRestrictProlongOperator = typename Problem::RestrictProlongOperator;
// discrete function restriction and prolongation operator for adaptation
using DiscreteFunctionRestrictProlong = Dune::Fem::RestrictProlongDefault<DiscreteFunction>;
using RestrictProlong
= Dune::Fem::RestrictProlongTuple<DiscreteFunctionRestrictProlong, ProblemRestrictProlongOperator>;
// adaptation classes
using AdaptationManager = Dune::Fem::AdaptationManager<Grid, RestrictProlong>;
public:
template<class Serializer>
struct SerializeHelper {
template<class SolutionType>
static void serializeOp(Serializer& serializer,
SolutionType& solution)
{
for (auto& sol : solution) {
serializer(sol->blockVector());
}
}
};
FvBaseDiscretizationFemAdapt(Simulator& simulator)
: ParentType(simulator)
, space_(simulator.vanguard().gridPart())
{
if (this->enableGridAdaptation_ && !Dune::Fem::Capabilities::isLocallyAdaptive<Grid>::v) {
throw std::invalid_argument("Grid adaptation enabled, but chosen Grid is not capable"
" of adaptivity");
}
for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
this->solution_[timeIdx] = std::make_unique<DiscreteFunction>("solution", space_);
}
}
void adaptGrid()
{
// adapt the grid if enabled and if all dependencies are available
// adaptation is only done if markForGridAdaptation returns true
if (this->enableGridAdaptation_) {
// check if problem allows for adaptation and cells were marked
if (this->simulator_.problem().markForGridAdaptation()) {
// adapt the grid and load balance if necessary
adaptationManager().adapt();
// if the grid has potentially changed, we need to re-create the
// supporting data structures.
this->elementMapper_.update(this->gridView_);
this->vertexMapper_.update(this->gridView_);
this->resetLinearizer();
// this is a bit hacky because it supposes that Problem::finishInit()
// works fine multiple times in a row.
//
// TODO: move this to Problem::gridChanged()
this->finishInit();
// notify the problem that the grid has changed
//
// TODO: come up with a mechanism to access the unadapted data structures
// outside of the problem (i.e., grid, mappers, solutions)
this->simulator_.problem().gridChanged();
// notify the modules for visualization output
auto outIt = this->outputModules_.begin();
auto outEndIt = this->outputModules_.end();
for (; outIt != outEndIt; ++outIt)
(*outIt)->allocBuffers();
}
}
}
AdaptationManager& adaptationManager()
{
if (!adaptationManager_) {
// create adaptation objects here, because when doing so in constructor
// problem is not yet intialized, aka seg fault
restrictProlong_ = std::make_unique<RestrictProlong>(DiscreteFunctionRestrictProlong(*(this->solution_[/*timeIdx=*/0])),
this->simulator_.problem().restrictProlongOperator());
adaptationManager_ = std::make_unique<AdaptationManager>(this->simulator_.vanguard().grid(), *restrictProlong_);
}
return *adaptationManager_;
}
private:
DiscreteFunctionSpace space_;
std::unique_ptr<RestrictProlong> restrictProlong_;
std::unique_ptr<AdaptationManager> adaptationManager_;
};
} // namespace Opm
#endif

3
opm/models/discretization/common/fvbaseproperties.hh
View File

@ -90,6 +90,9 @@ struct Stencil { using type = UndefinedProperty; };
template<class TypeTag, class MyTypeTag>
struct DiscreteFunctionSpace { using type = UndefinedProperty; };
template<class TypeTag, class MyTypeTag>
struct DiscreteFunction { using type = UndefinedProperty; };
//! The type of the problem
template<class TypeTag, class MyTypeTag>
struct Problem { using type = UndefinedProperty; };

18
opm/models/discretization/ecfv/ecfvdiscretization.hh
View File

@ -39,6 +39,7 @@
#include <opm/models/discretization/common/fvbasediscretization.hh>
#if HAVE_DUNE_FEM
#include <opm/models/discretization/common/fvbasediscretizationfemadapt.hh>
#include <dune/fem/space/common/functionspace.hh>
#include <dune/fem/space/finitevolume.hh>
#endif
@ -97,6 +98,18 @@ private:
public:
using type = Dune::Fem::FiniteVolumeSpace< FunctionSpace, GridPart, 0 >;
};
#else
template <class TypeTag>
struct DummySpaceEcfv {
using DiscreteFunctionSpace = GetPropType<TypeTag, Properties::DiscreteFunctionSpace>;
DummySpaceEcfv(const DiscreteFunctionSpace&) {};
DummySpaceEcfv(const int&) {};
};
template <class TypeTag>
struct DiscreteFunctionSpace<TypeTag, TTag::EcfvDiscretization> {
using type = DummySpaceEcfv<TypeTag>;
};
#endif
//! Set the border list creator for to the one of an element based
@ -130,10 +143,9 @@ namespace Opm {
* \brief The base class for the element-centered finite-volume discretization scheme.
*/
template<class TypeTag>
class EcfvDiscretization : public FvBaseDiscretization<TypeTag>
class EcfvDiscretization : public GetPropType<TypeTag, Properties::BaseDiscretizationType>
{
using ParentType = FvBaseDiscretization<TypeTag>;
using ParentType = GetPropType<TypeTag, Properties::BaseDiscretizationType>;
using Implementation = GetPropType<TypeTag, Properties::Model>;
using DofMapper = GetPropType<TypeTag, Properties::DofMapper>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;

17
opm/models/discretization/vcfv/vcfvdiscretization.hh
View File

@ -40,6 +40,7 @@
#include <opm/models/discretization/common/fvbasediscretization.hh>
#if HAVE_DUNE_FEM
#include <opm/models/discretization/common/fvbasediscretizationfemadapt.hh>
#include <dune/fem/space/common/functionspace.hh>
#include <dune/fem/space/lagrange.hh>
#endif
@ -109,6 +110,18 @@ public:
// Lagrange discrete function space with unknowns at the cell vertices
using type = Dune::Fem::LagrangeDiscreteFunctionSpace< FunctionSpace, GridPart, 1 >;
};
#else
template <class TypeTag>
struct DummySpaceVcfv {
using DiscreteFunctionSpace = GetPropType<TypeTag, Properties::DiscreteFunctionSpace>;
DummySpaceVcfv(const DiscreteFunctionSpace&) {};
DummySpaceVcfv(const int&) {};
};
template <class TypeTag>
struct DiscreteFunctionSpace<TypeTag, TTag::VcfvDiscretization> {
using type = DummySpaceVcfv<TypeTag>;
};
#endif
//! Set the border list creator for vertices
@ -138,9 +151,9 @@ namespace Opm {
* \brief The base class for the vertex centered finite volume discretization scheme.
*/
template<class TypeTag>
class VcfvDiscretization : public FvBaseDiscretization<TypeTag>
class VcfvDiscretization : public GetPropType<TypeTag, Properties::BaseDiscretizationType>
{
using ParentType = FvBaseDiscretization<TypeTag>;
using ParentType = GetPropType<TypeTag, Properties::BaseDiscretizationType>;
using Implementation = GetPropType<TypeTag, Properties::Model>;
using DofMapper = GetPropType<TypeTag, Properties::DofMapper>;
using GridView = GetPropType<TypeTag, Properties::GridView>;

18
opm/models/io/basevanguard.hh
View File

@ -120,15 +120,17 @@ protected:
void updateGridView_()
{
#if HAVE_DUNE_FEM
// first delete old grid part
// this is due to a bug in dune-fem (dangling reference)
gridPart_.reset();
gridPart_.reset(new GridPart(asImp_().grid()));
gridView_.reset(new GridView(static_cast<GridView>(*gridPart_)));
assert(gridView_->size(0) == asImp_().grid().leafGridView().size(0));
#else
gridView_.reset(new GridView(asImp_().grid().leafGridView()));
if constexpr (std::is_same_v<GridView,
typename GetPropType<TypeTag,
Properties::GridPart>::GridViewType>) {
gridPart_ = std::make_unique<GridPart>(asImp_().grid());
gridView_ = std::make_unique<GridView>(static_cast<GridView>(*gridPart_));
assert(gridView_->size(0) == asImp_().grid().leafGridView().size(0));
} else
#endif
{
gridView_ = std::make_unique<GridView>(asImp_().grid().leafGridView());
}
}
private: