Skip to content
GitLab
Commit b3c7e76f authored by Timo Koch's avatar Timo Koch Committed by Dennis Gläser
Browse files

[next] Make first example compile testing only the fvgeometry

parent 85fdc156
Expand all Hide whitespace changes
Inline Side-by-side
dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh
View file @ b3c7e76f
......@@ -137,7 +137,7 @@ public:
void bindElement(const Element& element)
{
elementPtr_ = &element;
scvIndices_ = std::vector<IndexType>({fvGridGeometry().problem_().elementMapper().index(*elementPtr_)});
scvIndices_ = std::vector<IndexType>({fvGridGeometry().elementMapper().index(*elementPtr_)});
}
//! The global finite volume geometry we are a restriction of
......@@ -299,21 +299,20 @@ public:
}
}
//! Check if user added additional DOF dependencies, i.e. the residual of DOF globalI depends
//! TODO Check if user added additional DOF dependencies, i.e. the residual of DOF globalI depends
//! on additional DOFs not included in the discretization schemes' occupation pattern
const auto& problem = fvGridGeometry().problem_();
const auto globalI = problem.elementMapper().index(element);
const auto& additionalDofDependencies = problem.getAdditionalDofDependencies(globalI);
if (!additionalDofDependencies.empty())
{
neighborScvs_.reserve(neighborScvs_.size() + additionalDofDependencies.size());
neighborScvIndices_.reserve(neighborScvIndices_.size() + additionalDofDependencies.size());
for (auto globalJ : additionalDofDependencies)
{
neighborScvs_.emplace_back(fvGridGeometry().element(globalJ).geometry(), globalJ);
neighborScvIndices_.emplace_back(globalJ);
}
}
// const auto globalI = fvGridGeometry().elementMapper().index(element);
// const auto& additionalDofDependencies = problem.getAdditionalDofDependencies(globalI);
// if (!additionalDofDependencies.empty())
// {
// neighborScvs_.reserve(neighborScvs_.size() + additionalDofDependencies.size());
// neighborScvIndices_.reserve(neighborScvIndices_.size() + additionalDofDependencies.size());
// for (auto globalJ : additionalDofDependencies)
// {
// neighborScvs_.emplace_back(fvGridGeometry().element(globalJ).geometry(), globalJ);
// neighborScvIndices_.emplace_back(globalJ);
// }
// }
}
//! Binding of an element preparing the geometries only inside the element
......@@ -361,7 +360,7 @@ private:
//! create scvs and scvfs of the bound element
void makeElementGeometries(const Element& element)
{
auto eIdx = fvGridGeometry().problem_().elementMapper().index(element);
const auto eIdx = fvGridGeometry().elementMapper().index(element);
scvs_.emplace_back(element.geometry(), eIdx);
scvIndices_.emplace_back(eIdx);
......@@ -378,8 +377,8 @@ private:
int scvfCounter = 0;
for (const auto& intersection : intersections(fvGridGeometry().gridView(), element))
{
// check if intersection is on interior boundary
const auto isInteriorBoundary = fvGridGeometry().problem_().isInteriorBoundary(element, intersection);
// TODO check if intersection is on interior boundary
const auto isInteriorBoundary = false;
if (dim < dimWorld)
if (handledScvf[intersection.indexInInside()])
......@@ -408,7 +407,7 @@ private:
void makeNeighborGeometries(const Element& element)
{
// create the neighbor scv
auto eIdx = fvGridGeometry().problem_().elementMapper().index(element);
const auto eIdx = fvGridGeometry().elementMapper().index(element);
neighborScvs_.emplace_back(element.geometry(), eIdx);
neighborScvIndices_.push_back(eIdx);
......@@ -425,8 +424,8 @@ private:
int scvfCounter = 0;
for (const auto& intersection : intersections(fvGridGeometry().gridView(), element))
{
// check if intersection is on interior boundary
const auto isInteriorBoundary = fvGridGeometry().problem_().isInteriorBoundary(element, intersection);
// TODO check if intersection is on interior boundary
const auto isInteriorBoundary = false;
if (dim < dimWorld)
if (handledScvf[intersection.indexInInside()])
......@@ -458,7 +457,7 @@ private:
{
for (unsigned outsideScvIdx = 0; outsideScvIdx < neighborVolVarIndices[scvfCounter].size(); ++outsideScvIdx)
{
if (neighborVolVarIndices[scvfCounter][outsideScvIdx] == fvGridGeometry().problem_().elementMapper().index(*elementPtr_))
if (neighborVolVarIndices[scvfCounter][outsideScvIdx] == fvGridGeometry().elementMapper().index(*elementPtr_))
{
std::vector<IndexType> scvIndices({eIdx});
scvIndices.insert(scvIndices.end(), neighborVolVarIndices[scvfCounter].begin(), neighborVolVarIndices[scvfCounter].end());
......
dumux/implicit/cellcentered/assembler.hh
View file @ b3c7e76f
......@@ -32,15 +32,21 @@
namespace Dumux {
template<class TypeTag, DiscretizationMethods DM>
class ImplicitAssemblerImplementation;
template<class TypeTag>
class CCImplicitAssembler;
//! TPFA and MPFA use the same assembler class
template<class TypeTag>
class ImplicitAssemblerImplementation<TypeTag, DiscretizationMethods::CCTpfa> : CCImplicitAssembler<TypeTag>
class ImplicitAssemblerImplementation<TypeTag, DiscretizationMethods::CCTpfa> : public CCImplicitAssembler<TypeTag>
{
using CCImplicitAssembler<TypeTag>::CCImplicitAssembler;
};
template<class TypeTag>
class ImplicitAssemblerImplementation<TypeTag, DiscretizationMethods::CCMpfa> : CCImplicitAssembler<TypeTag>
class ImplicitAssemblerImplementation<TypeTag, DiscretizationMethods::CCMpfa> : public CCImplicitAssembler<TypeTag>
{
using CCImplicitAssembler<TypeTag>::CCImplicitAssembler;
};
......@@ -54,21 +60,23 @@ template<class TypeTag>
class CCImplicitAssembler
{
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using AssemblyMap = typename GET_PROP_TYPE(TypeTag, AssemblyMap);
using LocalResidual = typename GET_PROP_TYPE(TypeTag, LocalResidual);
using GridFvGeometry = typename GET_PROP_TYPE(TypeTag, GridFvGeometry);
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using LocalAssembler = typename GET_PROP_TYPE(TypeTag, LocalAssembler);
// using LocalAssembler = typename GET_PROP_TYPE(TypeTag, LocalAssembler);
using LocalAssembler = CCImplicitLocalAssembler<TypeTag, DifferentiationMethods::numeric>;
public:
using ResidualType = SolutionVector;
//! The constructor
CCImplicitAssembler(std::shared_ptr<const Problem> problem,
std::shared_ptr<const GridFvGeometry> gridFvGeometry,
std::shared_ptr<const FVGridGeometry> gridFvGeometry,
std::shared_ptr<const SolutionVector> prevSol,
std::shared_ptr<const SolutionVector> curSol,
std::shared_ptr<GridVariables> gridVariables)
......@@ -79,7 +87,7 @@ public:
, gridVariables_(gridVariables)
{
// build assembly map
assemblyMap_.init(globalFvGeometry);
assemblyMap_.init(gridFvGeometry);
}
/*!
......@@ -100,18 +108,18 @@ public:
// if we get here, everything worked well
succeeded = true;
if (gridView_().comm().size() > 1)
succeeded = gridView_().comm().min(succeeded);
if (gridView().comm().size() > 1)
succeeded = gridView().comm().min(succeeded);
}
// throw exception if a problem ocurred
catch (NumericalProblem &e)
{
std::cout << "rank " << problem_().gridView().comm().rank()
std::cout << "rank " << gridView().comm().rank()
<< " caught an exception while assembling:" << e.what()
<< "\n";
succeeded = false;
if (gridView_().comm().size() > 1)
succeeded = gridView_().comm().min(succeeded);
if (gridView().comm().size() > 1)
succeeded = gridView().comm().min(succeeded);
}
if (!succeeded)
DUNE_THROW(NumericalProblem, "A process did not succeed in linearizing the system");
......@@ -134,8 +142,8 @@ public:
// if we get here, everything worked well
succeeded = true;
if (gridView_().comm().size() > 1)
succeeded = gridView_().comm().min(succeeded);
if (gridView().comm().size() > 1)
succeeded = gridView().comm().min(succeeded);
}
// throw exception if a problem ocurred
catch (NumericalProblem &e)
......@@ -144,8 +152,8 @@ public:
<< " caught an exception while assembling:" << e.what()
<< "\n";
succeeded = false;
if (gridView_().comm().size() > 1)
succeeded = gridView_().comm().min(succeeded);
if (gridView().comm().size() > 1)
succeeded = gridView().comm().min(succeeded);
}
if (!succeeded)
DUNE_THROW(NumericalProblem, "A process did not succeed in linearizing the system");
......@@ -229,7 +237,7 @@ public:
// calculate the square norm of the residual
Scalar result2 = residual.two_norm2();
if (gridView().comm().size() > 1)
result2 = gridView_().comm().sum(result2);
result2 = gridView().comm().sum(result2);
using std::sqrt;
return sqrt(result2);
......@@ -238,7 +246,7 @@ public:
const Problem& problem() const
{ return *problem; }
const GridFvGeometry& gridFvGeometry() const
const FVGridGeometry& gridFvGeometry() const
{ return *gridFvGeometry_; }
const GridView& gridView() const
......@@ -302,7 +310,7 @@ private:
std::shared_ptr<const Problem> problem_;
// the finite volume geometry of the grid
std::shared_ptr<const GridFvGeometry> gridFvGeometry_;
std::shared_ptr<const FVGridGeometry> gridFvGeometry_;
// previous and current solution to the problem
std::shared_ptr<const SolutionVector> prevSol_;
......
dumux/implicit/cellcentered/assemblymap.hh
View file @ b3c7e76f
......@@ -47,7 +47,7 @@ class CCSimpleAssemblyMap
{
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using IndexType = typename GridView::IndexSet::IndexType;
struct DataJ
......
dumux/implicit/cellcentered/localassembler.hh
View file @ b3c7e76f
......@@ -30,6 +30,11 @@
namespace Dumux {
enum DifferentiationMethods
{
numeric
};
/*!
* \ingroup ImplicitModel
* \brief An assembler for the local contributions (per element) to the global
......@@ -47,6 +52,10 @@ class CCImplicitLocalAssembler<TypeTag, DifferentiationMethods::numeric>
using Implementation = typename GET_PROP_TYPE(TypeTag, LocalAssembler);
using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
using Element = typename GET_PROP_TYPE(TypeTag, GridView)::template Codim<0>::Entity;
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
public:
......@@ -57,6 +66,7 @@ public:
template<class Assembler>
static void assemble(Assembler& assembler, SolutionVector& res, const Element& element)
{
const auto globalI = assembler.fvGridGeometry().elementMapper().index(element);
res[globalI] = Implementation::assemble_(assembler, element);
}
......@@ -76,7 +86,7 @@ public:
*
* \param priVar The value of the primary variable
*/
static Scalar numericEpsilon(const Scalar priVar) const
static Scalar numericEpsilon(const Scalar priVar)
{
// define the base epsilon as the geometric mean of 1 and the
// resolution of the scalar type. E.g. for standard 64 bit
......@@ -115,7 +125,7 @@ private:
auto& A = assembler.matrix();
// prepare the local views
auto fvGeometry = localView(globalFvGeometry());
auto fvGeometry = localView(assembler.fvGridGeometry());
fvGeometry.bind(element);
auto curElemVolVars = localView(gridVariables.curGlobalVolVars());
......@@ -251,7 +261,7 @@ private:
neighborDeriv = origFlux;
}
if (numericDifferenceMethod_ <= 0)
if (numericDifferenceMethod <= 0)
{
// we are not using forward differences, i.e. we
// need to calculate f(x - \epsilon)
......@@ -266,7 +276,8 @@ private:
// calculate the residual with the deflected primary variables and subtract it
if (!isGhost)
partialDeriv -= localResidual.eval(element,
partialDeriv -= localResidual.eval(problem,
element,
fvGeometry,
prevElemVolVars,
curElemVolVars,
......@@ -275,12 +286,13 @@ private:
// calculate the fluxes into element with the deflected primary variables
for (std::size_t k = 0; k < numNeighbors; ++k)
for (auto scvfIdx : assemblyMap_[globalI][k].scvfsJ)
neighborDeriv[k] -= this->localResidual().evalFlux_(neighborElements[k],
fvGeometry,
curElemVolVars,
fvGeometry.scvf(scvfIdx),
elemFluxVarsCache);
for (auto scvfIdx : assemblyMap[globalI][k].scvfsJ)
neighborDeriv[k] -= localResidual.evalFlux(problem,
neighborElements[k],
fvGeometry,
curElemVolVars,
fvGeometry.scvf(scvfIdx),
elemFluxVarsCache);
}
else
{
......@@ -312,7 +324,7 @@ private:
// off-diagonal entries
j = 0;
for (const auto& dataJ : assemblyMap_[globalI])
for (const auto& dataJ : assemblyMap[globalI])
A[dataJ.globalJ][globalI][eqIdx][pvIdx] += neighborDeriv[j++][pvIdx];
}
}
......@@ -325,101 +337,103 @@ private:
// //
//////////////////////////////////////////////////////////////////////////////////////////////
const auto& additionalDofDepedencies = problem.getAdditionalDofDependencies(globalI);
if (!additionalDofDepedencies.empty() && !isGhost)
{
// compute the source in the undeflected state
auto source = localResidual.computeSource(element, fvGeometry, curElemVolVars, scv);
source *= -scv.volume()*curVolVarsI.extrusionFactor();
// deflect solution at given dofs and recalculate the source
for (auto globalJ : additionalDofDependencies)
{
const auto& scvJ = fvGeometry.scv(globalJ);
auto& curVolVarsJ = curElemVolVars[scv];
const auto& elementJ = gridFvGeometry.element(globalJ);
// save a copy of the original privars and volvars
// to restore original solution after deflection
const auto origPriVars = curSol[globalJ];
const auto origVolVarsJ = curVolVarsJ;
// derivatives with repect to the additional DOF we depend on
for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
{
// derivatives of element dof with respect to itself
NumEqVector partialDeriv(0.0);
const auto eps = Implementation::numericEpsilon(curVolVarsJ.priVar(pvIdx));
Scalar delta = 0;
if (numericDifferenceMethod_ >= 0)
{
// we are not using backward differences, i.e. we need to
// calculate f(x + \epsilon)
// deflect primary variables
curSol[globalJ][pvIdx] += eps;
delta += eps;
// update the volume variables and the flux var cache
curVolVarsJ.update(gridVariables.elementSolution(elementJ, curSol), problem, elementJ, scvJ);
// calculate the source with the deflected primary variables
auto deflSource = localResidual.computeSource(element, fvGeometry, curElemVolVars, scv);
deflSource *= -scv.volume()*curVolVarsI.extrusionFactor();
partialDeriv = std::move(deflSource);
}
else
{
// we are using backward differences, i.e. we don't need
// to calculate f(x + \epsilon) and we can recycle the
// (already calculated) source f(x)
partialDeriv = source;
}
if (numericDifferenceMethod_ <= 0)
{
// we are not using forward differences, i.e. we
// need to calculate f(x - \epsilon)
// deflect the primary variables
curSol[globalJ][pvIdx] -= delta + eps;
delta += eps;
// update the volume variables and the flux var cache
curVolVarsJ.update(gridVariables.elementSolution(elementJ, curSol), problem, elementJ, scvJ);
// calculate the source with the deflected primary variables and subtract
auto deflSource = localResidual.computeSource(element, fvGeometry, curElemVolVars, scv);
deflSource *= -scv.volume()*curVolVarsI.extrusionFactor();
partialDeriv -= std::move(deflSource);
}
else
{
// we are using forward differences, i.e. we don't need to
// calculate f(x - \epsilon) and we can recycle the
// (already calculated) source f(x)
partialDeriv -= source;
}
// divide difference in residuals by the magnitude of the
// deflections between the two function evaluation
partialDeriv /= delta;
// restore the original state of the dofs privars and the volume variables
curSol[globalJ] = origPriVars;
curVolVarsJ = origVolVarsJ;
// add the current partial derivatives to the global jacobian matrix
for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
A[globalI][globalJ][eqIdx][pvIdx] += partialDeriv[eqIdx];
}
}
}
// const auto& additionalDofDepedencies = problem.getAdditionalDofDependencies(globalI);
// if (!additionalDofDepedencies.empty() && !isGhost)
// {
// // compute the source in the undeflected state
// auto source = localResidual.computeSource(element, fvGeometry, curElemVolVars, scv);
// source *= -scv.volume()*curVolVarsI.extrusionFactor();
// // deflect solution at given dofs and recalculate the source
// for (auto globalJ : additionalDofDependencies)
// {
// const auto& scvJ = fvGeometry.scv(globalJ);
// auto& curVolVarsJ = curElemVolVars[scv];
// const auto& elementJ = gridFvGeometry.element(globalJ);
// // save a copy of the original privars and volvars
// // to restore original solution after deflection
// const auto origPriVars = curSol[globalJ];
// const auto origVolVarsJ = curVolVarsJ;
// // derivatives with repect to the additional DOF we depend on
// for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
// {
// // derivatives of element dof with respect to itself
// NumEqVector partialDeriv(0.0);
// const auto eps = Implementation::numericEpsilon(curVolVarsJ.priVar(pvIdx));
// Scalar delta = 0;
// if (numericDifferenceMethod >= 0)
// {
// // we are not using backward differences, i.e. we need to
// // calculate f(x + \epsilon)
// // deflect primary variables
// curSol[globalJ][pvIdx] += eps;
// delta += eps;
// // update the volume variables and the flux var cache
// curVolVarsJ.update(gridVariables.elementSolution(elementJ, curSol), problem, elementJ, scvJ);
// // calculate the source with the deflected primary variables
// auto deflSource = localResidual.computeSource(element, fvGeometry, curElemVolVars, scv);
// deflSource *= -scv.volume()*curVolVarsI.extrusionFactor();
// partialDeriv = std::move(deflSource);
// }
// else
// {
// // we are using backward differences, i.e. we don't need
// // to calculate f(x + \epsilon) and we can recycle the
// // (already calculated) source f(x)
// partialDeriv = source;
// }
// if (numericDifferenceMethod <= 0)
// {
// // we are not using forward differences, i.e. we
// // need to calculate f(x - \epsilon)
// // deflect the primary variables
// curSol[globalJ][pvIdx] -= delta + eps;
// delta += eps;
// // update the volume variables and the flux var cache
// curVolVarsJ.update(gridVariables.elementSolution(elementJ, curSol), problem, elementJ, scvJ);
// // calculate the source with the deflected primary variables and subtract
// auto deflSource = localResidual.computeSource(element, fvGeometry, curElemVolVars, scv);
// deflSource *= -scv.volume()*curVolVarsI.extrusionFactor();
// partialDeriv -= std::move(deflSource);
// }
// else
// {
// // we are using forward differences, i.e. we don't need to
// // calculate f(x - \epsilon) and we can recycle the
// // (already calculated) source f(x)
// partialDeriv -= source;
// }
// // divide difference in residuals by the magnitude of the
// // deflections between the two function evaluation
// partialDeriv /= delta;
// // restore the original state of the dofs privars and the volume variables
// curSol[globalJ] = origPriVars;
// curVolVarsJ = origVolVarsJ;
// // add the current partial derivatives to the global jacobian matrix
// for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
// A[globalI][globalJ][eqIdx][pvIdx] += partialDeriv[eqIdx];
// }
// }
// }
// return the original residual
return residual;
}
};
}
\ No newline at end of file
} // end namespace Dumux
#endif
dumux/implicit/cellcentered/localresidual.hh
View file @ b3c7e76f
......@@ -44,262 +44,62 @@ template<class TypeTag>
class CCLocalResidual : public ImplicitLocalResidual<TypeTag>
{
using ParentType = ImplicitLocalResidual<TypeTag>;
friend class ImplicitLocalResidual<TypeTag>;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
using Element = typename GridView::template Codim<0>::Entity;
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using Element = typename GET_PROP_TYPE(TypeTag, GridView)::template Codim<0>::Entity;
using ElementResidualVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
public:
// copying the local residual class is not a good idea
CCLocalResidual(const CCLocalResidual &) = delete;
CCLocalResidual() : ParentType() {}
protected:
void evalFluxes_(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementBoundaryTypes& bcTypes,
const ElementFluxVariablesCache& elemFluxVarsCache)
void evalFlux(ElementResidualVector& residual,
const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementBoundaryTypes& elemBcTypes,
const ElementFluxVariablesCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
// calculate the mass flux over the scv faces
for (auto&& scvf : scvfs(fvGeometry))
{
this->residual_[0] += this->asImp_().computeFlux_(element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
}
}
const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
const auto localScvIdx = scv.indexInElement();
PrimaryVariables computeFlux_(const Element &element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace &scvf,
const ElementFluxVariablesCache& elemFluxVarsCache)
{
// inner faces
if (!scvf.boundary())
return this->asImp_().computeFlux(element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
else
return PrimaryVariables(0.0);
}
void evalBoundary_(const Element &element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementBoundaryTypes& elemBcTypes,
const ElementFluxVariablesCache& elemFluxVarsCache)
{