Commit f0098df5 authored by Dennis Gläser's avatar Dennis Gläser
Browse files

[assembly] add fv local operator

parent 38126113
Pipeline #19950 waiting for manual action with stages
in 41 seconds
// -*- 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 Assembly
* \brief Local operator for finite-volume schemes.
*/
#ifndef DUMUX_ASSEMBLY_FV_LOCAL_OPERATOR_HH
#define DUMUX_ASSEMBLY_FV_LOCAL_OPERATOR_HH
#include <utility>
#include <concepts>
#include <dune/common/fvector.hh>
#include <dune/common/reservedvector.hh>
#include <dumux/experimental/new_assembly/dumux/assembly/concepts.hh>
#include <dumux/experimental/new_assembly/dumux/common/concepts.hh>
#include <dumux/experimental/new_assembly/dumux/common/variables.hh>
namespace Dumux {
template<int numEq,
typename EG,
typename EV,
Concepts::ScvOperator<EG, EV> SourceOperator,
Concepts::ScvOperator<EG, EV> StorageOperator,
Concepts::ScvfOperator<EG, EV> FluxOperator>
class FVLocalOperator
{
using Scalar = Variables::ScalarType<typename EV::GridVariables>;
using NumEqVector = Dune::FieldVector<Scalar, numEq>;
public:
using ElementGeometry = EG;
using ElementVariables = EV;
using LocalResidual = Dune::ReservedVector<NumEqVector, ElementGeometry::maxNumElementScvs>;
FVLocalOperator(const ElementGeometry& elemGeom,
const ElementVariables& elemVars,
SourceOperator&& sourceOp,
StorageOperator&& storageOp,
FluxOperator&& fluxOp)
: elemGeom_(elemGeom)
, elemVars_(elemVars)
, sourceOperator_(std::move(sourceOp))
, storageOperator_(std::move(storageOp))
, fluxOperator_(std::move(fluxOp))
{}
void addStorage(LocalResidual& res) const
{
for (const auto& scv : scvs(elemGeom_))
addTo_(
res[scv.localDofIndex()],
storageOperator_(elemGeom_, elemVars_, scv)*scv.volume()
);
}
void addFluxes(LocalResidual& res) const
{
for (const auto& scvf : scvfs(elemGeom_))
addTo_(
res[insideLocalDofIndex_(scvf)],
fluxOperator_(elemGeom_, elemVars_, scvf)*scvf.area()
);
}
void addSources(LocalResidual& res) const
{
for (const auto& scv : scvs(elemGeom_))
addTo_(
res[scv.localDofIndex()],
sourceOperator_(elemGeom_, elemVars_, scv)*scv.volume()
);
}
void reset(LocalResidual& res) const
{
std::for_each(res.begin(), res.end(), [] (auto& v) {
v = 0.0;
});
}
LocalResidual getEmptyResidual() const
{
LocalResidual res;
res.resize(elemGeom_.numScv());
return res;
}
private:
template<typename SubControlVolumeFace>
auto insideLocalDofIndex_(const SubControlVolumeFace& scvf) const
{ return elemGeom_.scv(scvf.insideScvIdx()).localDofIndex(); }
template<typename Values> requires(
Concepts::AddAssignable<NumEqVector, Values>)
void addTo_(NumEqVector& result, const Values& values) const
{ result += values; }
const ElementGeometry& elemGeom_;
const ElementVariables& elemVars_;
SourceOperator sourceOperator_;
StorageOperator storageOperator_;
FluxOperator fluxOperator_;
};
template<int numEq, typename EG, typename EV, typename S, typename ST, typename F>
auto makeFVLocalOperator(const EG& eg, const EV& ev, S&& source, ST&& storage, F&& flux)
{
return FVLocalOperator<numEq, EG, EV, S, ST, F>{
eg, ev, std::move(source), std::move(storage), std::move(flux)
};
}
} // namespace Dumux
#endif
add_subdirectory(poisson)
dumux_add_test(SOURCES test_matrix_pattern.cc LABELS unit new_assembly)
dumux_add_test(SOURCES test_fv_local_operator.cc LABELS unit new_assembly)
dumux_add_test(SOURCES test_composite_operator.cc LABELS unit new_assembly)
#include <iostream>
#include <dune/common/exceptions.hh>
#include <dune/grid/yaspgrid.hh>
#include <dumux/common/initialize.hh>
#include <dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh>
#include <dumux/experimental/new_assembly/dumux/common/indexstrategies.hh>
#include <dumux/experimental/new_assembly/dumux/assembly/fvlocaloperator.hh>
#include <dumux/experimental/new_assembly/dumux/discretization/gridvariables.hh>
#include <dumux/experimental/new_assembly/dumux/linear/dune/labackend.hh>
template<typename GV>
class TestGridVariablesLocalView
{
public:
using GridVariables = GV;
TestGridVariablesLocalView(const GV& gv)
: gv_(gv)
{}
template<typename ElementGeometry>
TestGridVariablesLocalView bind(const ElementGeometry& eg) &&
{ return *this; }
template<typename ElementGeometry>
void bind(const ElementGeometry& eg) &
{}
private:
const GV& gv_;
};
template<typename GridGeometry, typename Dofs>
class TestGridVariables
: public Dumux::GridVariables<GridGeometry,
Dofs,
Dumux::BlockedIndexStrategy<1>>
{
using IS = Dumux::BlockedIndexStrategy<1>;
using ParentType = Dumux::GridVariables<GridGeometry, Dofs, IS>;
public:
using ParentType::ParentType;
using LocalView = TestGridVariablesLocalView<TestGridVariables>;
friend LocalView localView(const TestGridVariables& gv)
{ return {gv}; }
};
int main(int argc, char** argv)
{
Dumux::initialize(argc, argv);
static constexpr int dim = 2;
using Scalar = double;
// make grid
const Scalar domainSize = 0.5;
const Scalar cellVolume = domainSize*domainSize;
const Scalar boundaryArea = 4*domainSize;
using Grid = Dune::YaspGrid<dim>;
const Grid grid{
Dune::FieldVector<Scalar, dim>{{domainSize, domainSize}},
std::array<int, dim>{{1, 1}}
};
// make grid geometry of the scheme to be used
using GridGeometry = Dumux::CCTpfaFVGridGeometry<typename Grid::LeafGridView>;
auto gridGeometry = std::make_shared<GridGeometry>(grid.leafGridView());
// make grid variables
using Dofs = typename Dumux::DefaultBlockedLinearAlgebraBackend<Scalar, 0>::Vector;
using GridVariables = TestGridVariables<GridGeometry, Dofs>;
auto gridVariables = std::make_shared<GridVariables>(
gridGeometry,
std::make_shared<typename GridVariables::IndexStrategy>()
);
// build a local operator
static constexpr int numEq = 1;
const auto element = *elements(gridGeometry->gridView()).begin();
auto fvGeometry = localView(*gridGeometry).bind(element);
auto elemVars = localView(*gridVariables).bind(element);
auto localOperator = Dumux::makeFVLocalOperator<numEq>(
fvGeometry,
elemVars,
[] (auto... args) { return Dune::FieldVector<Scalar, numEq>{{1.5}}; },
[] (auto... args) { return Dune::FieldVector<Scalar, numEq>{{2.5}}; },
[] (auto... args) { return Dune::FieldVector<Scalar, numEq>{{3.5}}; }
);
auto result = localOperator.getEmptyResidual();
localOperator.addSources(result);
localOperator.addSources(result);
if (std::abs(result[0][0] - 3.0*cellVolume) > 1e-8)
DUNE_THROW(Dune::InvalidStateException, "Unexpected source operator");
localOperator.reset(result);
if (std::abs(result[0][0]) > 1e-8)
DUNE_THROW(Dune::InvalidStateException, "Reset failed");
localOperator.addStorage(result);
localOperator.addStorage(result);
if (std::abs(result[0][0] - 5.0*cellVolume) > 1e-8)
DUNE_THROW(Dune::InvalidStateException, "Unexpected storage operator");
localOperator.reset(result);
localOperator.addFluxes(result);
localOperator.addFluxes(result);
if (std::abs(result[0][0] - 7.0*boundaryArea) > 1e-8)
DUNE_THROW(Dune::InvalidStateException, "Unexpected flux operator");
std::cout << "Test passed" << std::endl;
return 0;
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment