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

[next] First working stationary 1p test with homogeneous Dirichlet

parent 8b44dfe0
Hide whitespace changes
Inline Side-by-side
dumux/discretization/cellcentered/globalvolumevariables.hh
View file @ 9cd50607
......@@ -156,7 +156,7 @@ class CCGlobalVolumeVariables<TypeTag, /*enableGlobalVolVarsCache*/false>
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
public:
void update(Problem& problem, const FVGridGeometry& fvGridGeometry, const SolutionVector& sol)
void update(const Problem& problem, const FVGridGeometry& fvGridGeometry, const SolutionVector& sol)
{ problemPtr_ = &problem; }
/*!
......@@ -168,10 +168,10 @@ public:
{ return ElementVolumeVariables(global); }
private:
Problem& problem_() const
const Problem& problem_() const
{ return *problemPtr_;}
Problem* problemPtr_;
const Problem* problemPtr_;
};
} // end namespace
......
dumux/implicit/cellcentered/assembler.hh
View file @ 9cd50607
......@@ -27,11 +27,18 @@
#include <dune/istl/matrixindexset.hh>
#include <dumux/implicit/properties.hh>
#include <dumux/discretization/methods.hh>
#include "localassembler.hh"
namespace Dumux {
namespace Properties
{
NEW_PROP_TAG(AssemblyMap);
NEW_PROP_TAG(LocalAssembler);
}
template<class TypeTag, DiscretizationMethods DM>
class ImplicitAssemblerImplementation;
......@@ -76,18 +83,18 @@ public:
//! The constructor
CCImplicitAssembler(std::shared_ptr<const Problem> problem,
std::shared_ptr<const FVGridGeometry> gridFvGeometry,
std::shared_ptr<const FVGridGeometry> fvGridGeometry,
std::shared_ptr<const SolutionVector> prevSol,
std::shared_ptr<const SolutionVector> curSol,
std::shared_ptr<GridVariables> gridVariables)
: problem_(problem)
, gridFvGeometry_(gridFvGeometry)
, fvGridGeometry_(fvGridGeometry)
, prevSol_(prevSol)
, curSol_(curSol)
, gridVariables_(gridVariables)
{
// build assembly map
assemblyMap_.init(gridFvGeometry);
assemblyMap_.init(*fvGridGeometry);
}
/*!
......@@ -137,7 +144,7 @@ public:
try
{
// let the local assembler add the element contributions
for (const auto element : elements(gridView()))
for (const auto& element : elements(gridView()))
LocalAssembler::assemble(*this, element);
// if we get here, everything worked well
......@@ -148,7 +155,7 @@ public:
// 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;
......@@ -164,11 +171,11 @@ public:
* This also resizes the containers to the required sizes and sets the
* sparsity pattern of the matrix.
*/
void setLinearSystem(std::shared_ptr<JacobianMatrix> matrix,
std::shared_ptr<SolutionVector> residual)
void setLinearSystem(std::shared_ptr<JacobianMatrix>& A,
std::shared_ptr<SolutionVector>& r)
{
matrix_ = matrix;
residual_ = residual;
matrix_ = A;
residual_ = r;
// check and/or set the BCRS matrix's build mode
if (matrix().buildMode() == JacobianMatrix::BuildMode::unknown)
......@@ -185,7 +192,7 @@ public:
void updateLinearSystem()
{
// resize the matrix and the residual
const auto numDofs = numDofs();
const auto numDofs = this->numDofs();
matrix().setSize(numDofs, numDofs);
residual().resize(numDofs);
......@@ -200,9 +207,9 @@ public:
occupationPattern.add(dataJ.globalJ, globalI);
// reserve index for additional user defined DOF dependencies
const auto& additionalDofDependencies = problem().getAdditionalDofDependencies(globalI);
for (auto globalJ : additionalDofDependencies)
occupationPattern.add(globalI, globalJ);
// const auto& additionalDofDependencies = problem().getAdditionalDofDependencies(globalI);
// for (auto globalJ : additionalDofDependencies)
// occupationPattern.add(globalI, globalJ);
}
// export pattern to matrix
......@@ -217,10 +224,10 @@ public:
{
if (element.partitionType != Dune::GhostEntity)
{
const auto eIdx = gridFvGeometry().elementMapper().index(element);
const auto eIdx = fvGridGeometry().elementMapper().index(element);
r[eIdx] += localResidual().eval(element,
problem(),
gridFvGeometry(),
fvGridGeometry(),
gridVariables(),
curSol(),
prevSol());
......@@ -244,13 +251,13 @@ public:
}
const Problem& problem() const
{ return *problem; }
{ return *problem_; }
const FVGridGeometry& gridFvGeometry() const
{ return *gridFvGeometry_; }
const FVGridGeometry& fvGridGeometry() const
{ return *fvGridGeometry_; }
const GridView& gridView() const
{ return gridFvGeometry_().gridView(); }
{ return fvGridGeometry().gridView(); }
const SolutionVector& prevSol() const
{ return *prevSol_; }
......@@ -269,8 +276,8 @@ public:
SolutionVector& residual()
{
assert(residual_ && "The right hand side has been set yet!");
return residual_;
assert(residual_ && "The residual has not been set yet!");
return *residual_;
}
const AssemblyMap& assemblyMap() const
......@@ -289,28 +296,19 @@ private:
// that the jacobian matrix must only be erased partially!
void resetSystem_()
{
residual_() = 0.0;
residual() = 0.0;
resetMatrix_();
}
void resetMatrix_()
{
matrix_() = 0.0;
}
/*!
* \brief Assemble the global Jacobian of the residual
* and the residual for the current solution.
*/
void assemble_()
{
matrix() = 0.0;
}
// pointer to the problem to be solved
std::shared_ptr<const Problem> problem_;
// the finite volume geometry of the grid
std::shared_ptr<const FVGridGeometry> gridFvGeometry_;
std::shared_ptr<const FVGridGeometry> fvGridGeometry_;
// previous and current solution to the problem
std::shared_ptr<const SolutionVector> prevSol_;
......
dumux/implicit/cellcentered/elementboundarytypes.hh
View file @ 9cd50607
......@@ -89,8 +89,8 @@ public:
for (auto&& scv : scvs(fvGeometry))
{
if (!problem.model().onBoundary(scv))
return;
// if (!problem.model().onBoundary(scv))
// return;
for (auto&& scvFace : scvfs(fvGeometry))
{
......
dumux/implicit/cellcentered/localassembler.hh
View file @ 9cd50607
......@@ -25,6 +25,7 @@
#define DUMUX_CC_IMPLICIT_LOCAL_ASSEMBLER_HH
#include <dune/istl/matrixindexset.hh>
#include <dune/istl/bvector.hh>
#include <dumux/implicit/properties.hh>
......@@ -54,6 +55,10 @@ class CCImplicitLocalAssembler<TypeTag, DifferentiationMethods::numeric>
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);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using GridVolumeVariables = typename GET_PROP_TYPE(TypeTag, GlobalVolumeVariables);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
......@@ -116,7 +121,7 @@ private:
{
// get some references for convenience
const auto& problem = assembler.problem();
const auto& gridFvGeometry = assembler.gridFvGeometry();
const auto& fvGridGeometry = assembler.fvGridGeometry();
const auto& assemblyMap = assembler.assemblyMap();
const auto& prevSol = assembler.prevSol();
auto& curSol = assembler.curSol();
......@@ -128,17 +133,17 @@ private:
auto fvGeometry = localView(assembler.fvGridGeometry());
fvGeometry.bind(element);
auto curElemVolVars = localView(gridVariables.curGlobalVolVars());
auto curElemVolVars = localView(gridVariables.curGridVolVars());
curElemVolVars.bind(element, fvGeometry, curSol);
auto prevElemVolVars = localView(gridVariables.prevGlobalVolVars());
auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
prevElemVolVars.bindElement(element, fvGeometry, prevSol);
auto elemFluxVarsCache = localView(gridVariables.globalFluxVarsCache());
auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
// the global dof of the actual element
const auto globalI = gridFvGeometry.elementMapper().index(element);
const auto globalI = fvGridGeometry.elementMapper().index(element);
// check for boundaries on the element
ElementBoundaryTypes elemBcTypes;
......@@ -148,12 +153,15 @@ private:
const bool isGhost = (element.partitionType() == Dune::GhostEntity);
// the actual element's current residual
const NumEqVector residual = isGhost ? 0.0 : localResidual.eval(element,
fvGeometry,
prevElemVolVars,
curElemVolVars,
elemBcTypes,
elemFluxVarsCache);
NumEqVector residual(0.0);
if (!isGhost)
residual = localResidual.eval(problem,
element,
fvGeometry,
prevElemVolVars,
curElemVolVars,
elemBcTypes,
elemFluxVarsCache)[0];
// TODO Do we really need this??????????
// this->model_().updatePVWeights(fvGeometry);
......@@ -182,20 +190,26 @@ private:
unsigned int j = 0;
for (const auto& dataJ : assemblyMap[globalI])
{
neighborElements.emplace_back(gridFvGeometry().element(dataJ.globalJ));
neighborElements.emplace_back(fvGridGeometry.element(dataJ.globalJ));
for (const auto scvfIdx : dataJ.scvfsJ)
origFlux[j] += localResidual.evalFlux_(neighborElements.back(),
fvGeometry,
curElemVolVars,
fvGeometry.scvf(scvfIdx),
elemFluxVarsCache);
{
ElementSolutionVector flux(1);
localResidual.evalFlux(flux, problem,
neighborElements.back(),
fvGeometry,
curElemVolVars,
elemBcTypes,
elemFluxVarsCache,
fvGeometry.scvf(scvfIdx));
origFlux[j] += flux[0];
}
// increment neighbor counter
++j;
}
// reference to the element's scv (needed later) and corresponding vol vars
const auto& scv = fvGeometry.scv(globalI);
auto& curVolVars = curElemVolVars[scv];
auto& curVolVars = getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scv);
// save a copy of the original privars and vol vars in order
// to restore the original solution after deflection
......@@ -235,21 +249,28 @@ private:
// calculate the residual with the deflected primary variables
if (!isGhost)
partialDeriv = localResidual.eval(element,
partialDeriv = localResidual.eval(problem,
element,
fvGeometry,
prevElemVolVars,
curElemVolVars,
elemBcTypes,
elemFluxVarsCache);
elemFluxVarsCache)[0];
// calculate the fluxes in the neighbors with the deflected primary variables
for (std::size_t k = 0; k < numNeighbors; ++k)
for (auto scvfIdx : assemblyMap[globalI][k].scvfsJ)
neighborDeriv[k] += localResidual.evalFlux_(neighborElements[k],
fvGeometry,
curElemVolVars,
fvGeometry.scvf(scvfIdx),
elemFluxVarsCache);
{
ElementSolutionVector flux(1);
localResidual.evalFlux(flux, problem,
neighborElements[k],
fvGeometry,
curElemVolVars,
elemBcTypes,
elemFluxVarsCache,
fvGeometry.scvf(scvfIdx));
neighborDeriv[k] += flux[0];
}
}
else
{
......@@ -282,17 +303,22 @@ private:
prevElemVolVars,
curElemVolVars,
elemBcTypes,
elemFluxVarsCache);
elemFluxVarsCache)[0];
// 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] -= localResidual.evalFlux(problem,
neighborElements[k],
fvGeometry,
curElemVolVars,
fvGeometry.scvf(scvfIdx),
elemFluxVarsCache);
{
ElementSolutionVector flux(1);
localResidual.evalFlux(flux, problem,
neighborElements[k],
fvGeometry,
curElemVolVars,
elemBcTypes,
elemFluxVarsCache,
fvGeometry.scvf(scvfIdx));
neighborDeriv[k] += flux[0];
}
}
else
{
......@@ -349,7 +375,7 @@ private:
// {
// const auto& scvJ = fvGeometry.scv(globalJ);
// auto& curVolVarsJ = curElemVolVars[scv];
// const auto& elementJ = gridFvGeometry.element(globalJ);
// const auto& elementJ = fvGridGeometry.element(globalJ);
// // save a copy of the original privars and volvars
// // to restore original solution after deflection
......@@ -432,6 +458,16 @@ private:
// return the original residual
return residual;
}
private:
template<class T = TypeTag>
static typename std::enable_if<!GET_PROP_VALUE(T, EnableGlobalVolumeVariablesCache), VolumeVariables&>::type
getVolVarAccess(GridVolumeVariables& gridVolVars, ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
{ return elemVolVars[scv]; }
template<class T = TypeTag>
static typename std::enable_if<GET_PROP_VALUE(T, EnableGlobalVolumeVariablesCache), VolumeVariables&>::type
getVolVarAccess(GridVolumeVariables& gridVolVars, ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
{ return gridVolVars.volVars(scv); }
};
} // end namespace Dumux
......
dumux/implicit/cellcentered/localresidual.hh
View file @ 9cd50607
......@@ -70,7 +70,7 @@ public:
// inner faces
if (!scvf.boundary())
{
residual[localScvIdx] += this->asImp_().computeFlux(element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
residual[localScvIdx] += this->asImp_().computeFlux(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
}
// boundary faces
......@@ -80,12 +80,12 @@ public:
// Dirichlet boundaries
if (bcTypes.hasDirichlet() && !bcTypes.hasNeumann())
residual[localScvIdx] += this->asImp_().computeFlux(element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
residual[localScvIdx] += this->asImp_().computeFlux(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
// Neumann and Robin ("solution dependent Neumann") boundary conditions
else if (bcTypes.hasNeumann() && !bcTypes.hasDirichlet())
{
auto neumannFluxes = this->problem().neumann(element, fvGeometry, elemVolVars, scvf);
auto neumannFluxes = problem.neumann(element, fvGeometry, elemVolVars, scvf);
// multiply neumann fluxes with the area and the extrusion factor
neumannFluxes *= scvf.area()*elemVolVars[scv].extrusionFactor();
......
dumux/implicit/localresidual.hh
View file @ 9cd50607
......@@ -163,7 +163,7 @@ public:
const ElementFluxVariablesCache& elemFluxVarsCache) const
{
// initialize the residual vector for all scvs in this element
ElementResidualVector residual(fvGeometry.numScv(), 0.0);
ElementResidualVector residual(fvGeometry.numScv());
// evaluate the volume terms (storage + source terms)
for (auto&& scv : scvs(fvGeometry))
......@@ -264,28 +264,28 @@ public:
const ElementVolumeVariables& prevElemVolVars,
const SubControlVolume& scv) const
{
const auto& curVolVars = curElemVolVars[scv];
const auto& prevVolVars = prevElemVolVars[scv];
// const auto& curVolVars = curElemVolVars[scv];
// const auto& prevVolVars = prevElemVolVars[scv];
// mass balance within the element. this is the
// \f$\frac{m}{\partial t}\f$ term if using implicit
// euler as time discretization.
//
// We might need a more explicit way for
// doing the time discretization...
// // mass balance within the element. this is the
// // \f$\frac{m}{\partial t}\f$ term if using implicit
// // euler as time discretization.
// //
// // We might need a more explicit way for
// // doing the time discretization...
//! Compute storage with the model specific storage residual
ResidualVector prevStorage = asImp_().computeStorage(problem, scv, prevVolVars);
ResidualVector storage = asImp_().computeStorage(problem, scv, curVolVars);
// //! Compute storage with the model specific storage residual
// ResidualVector prevStorage = asImp_().computeStorage(problem, scv, prevVolVars);
// ResidualVector storage = asImp_().computeStorage(problem, scv, curVolVars);
prevStorage *= prevVolVars.extrusionFactor();
storage *= curVolVars.extrusionFactor();
// prevStorage *= prevVolVars.extrusionFactor();
// storage *= curVolVars.extrusionFactor();
storage -= prevStorage;
storage *= scv.volume();
storage /= problem.timeManager().timeStepSize();
// storage -= prevStorage;
// storage *= scv.volume();
// storage /= problem.timeManager().timeStepSize();
residual[scv.indexInElement()] += storage;
// residual[scv.indexInElement()] += storage;
}
void evalSource(ElementResidualVector& residual,
......
test/porousmediumflow/1p/implicit/incompressible/properties.hh → test/porousmediumflow/1p/implicit/incompressible/problem.hh
View file @ 9cd50607
......@@ -20,21 +20,66 @@
* \file
* \brief The properties for the incompressible test
*/
#ifndef DUMUX_INCOMPRESSIBLE_ONEP_TEST_PROPERTIES_HH
#define DUMUX_INCOMPRESSIBLE_ONEP_TEST_PROPERTIES_HH
#ifndef DUMUX_INCOMPRESSIBLE_ONEP_TEST_PROBLEM_HH
#define DUMUX_INCOMPRESSIBLE_ONEP_TEST_PROBLEM_HH
#include <dumux/material/components/simpleh2o.hh>
#include <dumux/material/fluidsystems/liquidphase.hh>
#include <dumux/implicit/cellcentered/tpfa/properties.hh>
#include <dumux/implicit/cellcentered/localassembler.hh>
#include <dumux/implicit/gridvariables.hh>
#include <dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh>
#include <dumux/porousmediumflow/1p/implicit/propertydefaults.hh>
namespace Dumux
{
// forward declarations
template<class TypeTag> class OnePTestProblem;
template<class TypeTag> class OnePTestSpatialParams;
namespace Properties
{
NEW_PROP_TAG(EnableFVGridGeometryCache);
NEW_PROP_TAG(FVGridGeometry);
NEW_TYPE_TAG(IncompressibleTestProblem, INHERITS_FROM(CCTpfaModel, OneP));
// Set the grid type
SET_TYPE_PROP(IncompressibleTestProblem, Grid, Dune::YaspGrid<2>);
// Set the finite volume grid geometry
SET_TYPE_PROP(IncompressibleTestProblem, FVGridGeometry, CCTpfaFVGridGeometry<TypeTag, true>);
// Set the problem type
SET_TYPE_PROP(IncompressibleTestProblem, Problem, OnePTestProblem<TypeTag>);
SET_TYPE_PROP(IncompressibleTestProblem, SpatialParams, OnePTestSpatialParams<TypeTag>);
// the grid variables
SET_TYPE_PROP(IncompressibleTestProblem, GridVariables, GridVariables<TypeTag>);
// the local assembler
SET_TYPE_PROP(IncompressibleTestProblem, LocalAssembler, CCImplicitLocalAssembler<TypeTag, DifferentiationMethods::numeric>);
// the fluid system
SET_PROP(IncompressibleTestProblem, Fluid)
{
private:
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
public:
using type = FluidSystems::LiquidPhase<Scalar, SimpleH2O<Scalar> >;
};
// Enable caching
SET_BOOL_PROP(IncompressibleTestProblem, EnableGlobalVolumeVariablesCache, true);
SET_BOOL_PROP(IncompressibleTestProblem, EnableGlobalFluxVariablesCache, true);
SET_BOOL_PROP(IncompressibleTestProblem, EnableGlobalFVGeometryCache, true);
SET_BOOL_PROP(IncompressibleTestProblem, EnableFVGridGeometryCache, true);
} // end namespace Properties
template<class TypeTag>
class MockSpatialParams
class OnePTestSpatialParams
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
......@@ -61,20 +106,23 @@ public:
template<class TypeTag>
class MockProblem
class OnePTestProblem
{
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
using GlobalPosition = Dune::FieldVector<Scalar, GridView::dimension>;
public:
MockProblem(const GridView& gridView) {}
OnePTestProblem(const GridView& gridView) {}
/*!
* \brief Specifies which kind of boundary condition should be
......@@ -87,7 +135,10 @@ public:
const SubControlVolumeFace &scvf) const
{
BoundaryTypes values;
// const auto& pos = scvf.ipGlobal();
values.setAllDirichlet();
// if (pos[0] > 1.0 - 1e-8 || pos[0] < 1e-8)
// values.setAllNeumann();
return values;
}
......@@ -103,9 +154,70 @@ public:
PrimaryVariables dirichlet(const Element &element,
const SubControlVolumeFace &scvf) const
{
// const auto& pos = scvf.ipGlobal();
// return PrimaryVariables(pos[1]*1e5 + 1e5);
return PrimaryVariables(0.0);
}
/*!
* \brief Evaluate the boundary conditions for a dirichlet
* control volume.
*
* \param values The dirichlet values for the primary variables
* \param globalPos The center of the finite volume which ought to be set.
*
* For this method, the \a values parameter stores primary variables.
*/
ResidualVector neumann(const Element& element,
const FVElementGeometry& fvGeometry,