Skip to content
GitLab
Menu
Projects
Groups
Snippets
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in / Register
Toggle navigation
Menu
Open sidebar
dumux-repositories
dumux
Commits
e3fa9e4d
Commit
e3fa9e4d
authored
Dec 21, 2017
by
Dennis Gläser
Browse files
[discretization] improve docu in folder mpfa
parent
351d69a5
Changes
12
Hide whitespace changes
Inline
Side-by-side
dumux/discretization/cellcentered/mpfa/darcyslaw.hh
View file @
e3fa9e4d
...
...
@@ -228,15 +228,15 @@ public:
const
auto
&
tij
=
fluxVarsCache
.
advectionTij
();
const
auto
&
pj
=
fluxVarsCache
.
pressures
(
phaseIdx
);
//
!
compute t_ij*p_j
// compute t_ij*p_j
Scalar
scvfFlux
=
tij
*
pj
;
//
!
maybe add gravitational acceleration
// maybe add gravitational acceleration
static
const
bool
enableGravity
=
getParamFromGroup
<
bool
>
(
GET_PROP_VALUE
(
TypeTag
,
ModelParameterGroup
),
"Problem.EnableGravity"
);
if
(
enableGravity
)
scvfFlux
+=
fluxVarsCache
.
gravity
(
phaseIdx
);
//
!
switch the sign if necessary
// switch the sign if necessary
if
(
fluxVarsCache
.
advectionSwitchFluxSign
())
scvfFlux
*=
-
1.0
;
...
...
dumux/discretization/cellcentered/mpfa/dualgridindexset.hh
View file @
e3fa9e4d
...
...
@@ -21,8 +21,8 @@
* \ingroup CCMpfaDiscretization
* \brief Class for the index sets of the dual grid in mpfa schemes.
*/
#ifndef DUMUX_DISCRETIZATION_MPFA_DUALGRID_INDEX
SET_BA
SE_HH
#define DUMUX_DISCRETIZATION_MPFA_DUALGRID_INDEX
SET_BA
SE_HH
#ifndef DUMUX_DISCRETIZATION_MPFA_DUALGRID_INDEX
_
SE
T
_HH
#define DUMUX_DISCRETIZATION_MPFA_DUALGRID_INDEX
_
SE
T
_HH
namespace
Dumux
{
...
...
dumux/discretization/cellcentered/mpfa/fickslaw.hh
View file @
e3fa9e4d
...
...
@@ -151,7 +151,7 @@ class FicksLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
{
switchFluxSign_
[
phaseIdx
][
compIdx
]
=
localFaceData
.
isOutside
();
//
!
store pointer to the mole fraction vector of this iv
// store pointer to the mole fraction vector of this iv
xj_
[
phaseIdx
][
compIdx
]
=
&
dataHandle
.
moleFractions
(
phaseIdx
,
compIdx
);
const
auto
ivLocalIdx
=
localFaceData
.
ivLocalScvfIndex
();
...
...
dumux/discretization/cellcentered/mpfa/fourierslaw.hh
View file @
e3fa9e4d
...
...
@@ -141,7 +141,7 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
{
switchFluxSign_
=
localFaceData
.
isOutside
();
//
!
store pointer to the temperature vector of this iv
// store pointer to the temperature vector of this iv
Tj_
=
&
dataHandle
.
temperatures
();
const
auto
ivLocalIdx
=
localFaceData
.
ivLocalScvfIndex
();
...
...
@@ -189,7 +189,7 @@ public:
const
auto
&
tij
=
fluxVarsCache
.
heatConductionTij
();
const
auto
&
Tj
=
fluxVarsCache
.
temperatures
();
//
!
compute Tij*tj
// compute Tij*tj
Scalar
flux
=
tij
*
Tj
;
if
(
fluxVarsCache
.
heatConductionSwitchFluxSign
())
flux
*=
-
1.0
;
...
...
dumux/discretization/cellcentered/mpfa/helper.hh
View file @
e3fa9e4d
...
...
@@ -475,9 +475,9 @@ public:
{
case
3
:
// triangle
{
//
!
Only build the maps the first time we encounter a triangle
static
const
std
::
uint8_t
vo
=
1
;
//
!<
vertex offset in point vector p
static
const
std
::
uint8_t
eo
=
4
;
//
!<
edge offset in point vector p
// Only build the maps the first time we encounter a triangle
static
const
std
::
uint8_t
vo
=
1
;
// vertex offset in point vector p
static
const
std
::
uint8_t
eo
=
4
;
// edge offset in point vector p
static
const
std
::
uint8_t
map
[
3
][
4
]
=
{
{
0
,
eo
+
1
,
eo
+
0
,
vo
+
0
},
...
...
@@ -492,9 +492,9 @@ public:
}
case
4
:
// quadrilateral
{
//
!
Only build the maps the first time we encounter a quadrilateral
static
const
std
::
uint8_t
vo
=
1
;
//
!<
vertex offset in point vector p
static
const
std
::
uint8_t
eo
=
5
;
//
!<
face offset in point vector p
// Only build the maps the first time we encounter a quadrilateral
static
const
std
::
uint8_t
vo
=
1
;
// vertex offset in point vector p
static
const
std
::
uint8_t
eo
=
5
;
// face offset in point vector p
static
const
std
::
uint8_t
map
[
4
][
4
]
=
{
{
0
,
eo
+
0
,
eo
+
2
,
vo
+
0
},
...
...
dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh
View file @
e3fa9e4d
...
...
@@ -80,7 +80,7 @@ class CCMpfaInteractionVolumeBase
using
Element
=
typename
GridView
::
template
Codim
<
0
>
::
Entity
;
public:
//
!
Types required in the assembly of the local eq system
// Types required in the assembly of the local eq system
using
Problem
=
typename
GET_PROP_TYPE
(
TypeTag
,
Problem
);
using
FVElementGeometry
=
typename
GET_PROP_TYPE
(
TypeTag
,
FVElementGeometry
);
using
ElementVolumeVariables
=
typename
GET_PROP_TYPE
(
TypeTag
,
ElementVolumeVariables
);
...
...
dumux/discretization/cellcentered/mpfa/interactionvolumedatahandle.hh
View file @
e3fa9e4d
...
...
@@ -43,7 +43,7 @@ public:
template
<
class
TypeTag
,
class
InteractionVolume
,
bool
EnableAdvection
>
class
AdvectionDataHandle
{
//
!
export matrix & vector types from interaction volume
// export matrix & vector types from interaction volume
using
Matrix
=
typename
InteractionVolume
::
Traits
::
Matrix
;
using
Vector
=
typename
InteractionVolume
::
Traits
::
Vector
;
using
Scalar
=
typename
Vector
::
value_type
;
...
...
@@ -172,14 +172,14 @@ private:
std
::
array
<
Vector
,
numPhases
>
p_
;
//!< The interaction volume-wide phase pressures
std
::
array
<
Vector
,
numPhases
>
g_
;
//!< The gravitational acceleration at each scvf (only for enabled gravity)
std
::
vector
<
std
::
vector
<
Vector
>
>
outsideT_
;
//!< The transmissibilities for "outside" faces (only on surface grids)
std
::
array
<
std
::
vector
<
Vector
>
,
numPhases
>
outsideG_
;
//!< The gravitational acceleration on "outside" faces (only on surface grids)
std
::
array
<
std
::
vector
<
Vector
>
,
numPhases
>
outsideG_
;
//!< The gravitational acceleration on "outside" faces (only on surface grids)
};
//! Data handle for quantities related to diffusion
template
<
class
TypeTag
,
class
InteractionVolume
,
bool
EnableDiffusion
>
class
DiffusionDataHandle
{
//
!
export matrix & vector types from interaction volume
// export matrix & vector types from interaction volume
using
Matrix
=
typename
InteractionVolume
::
Traits
::
Matrix
;
using
Vector
=
typename
InteractionVolume
::
Traits
::
Vector
;
using
OutsideTContainer
=
std
::
vector
<
std
::
vector
<
Vector
>
>
;
...
...
@@ -205,11 +205,11 @@ public:
if
(
pIdx
==
cIdx
)
continue
;
//
!
resize transmissibility matrix & mole fraction vector
// resize transmissibility matrix & mole fraction vector
T_
[
pIdx
][
cIdx
].
resize
(
iv
.
numFaces
(),
iv
.
numKnowns
());
xj_
[
pIdx
][
cIdx
].
resize
(
iv
.
numKnowns
());
//
!
resize outsideTij on surface grids
// resize outsideTij on surface grids
if
(
dim
<
dimWorld
)
{
outsideT_
[
pIdx
][
cIdx
].
resize
(
iv
.
numFaces
());
...
...
@@ -236,8 +236,8 @@ public:
private:
//! diffusion-related variables
unsigned
int
phaseIdx_
;
//!< The phase index set for the context
unsigned
int
compIdx_
;
//!< The component index set for the context
unsigned
int
phaseIdx_
;
//!< The phase index set for the context
unsigned
int
compIdx_
;
//!< The component index set for the context
std
::
array
<
std
::
array
<
Matrix
,
numComponents
>
,
numPhases
>
T_
;
//!< The transmissibilities such that f_i = T_ij*x_j
std
::
array
<
std
::
array
<
Vector
,
numComponents
>
,
numPhases
>
xj_
;
//!< The interaction volume-wide mole fractions
std
::
array
<
std
::
array
<
OutsideTContainer
,
numComponents
>
,
numPhases
>
outsideT_
;
...
...
dumux/discretization/cellcentered/mpfa/localassembler.hh
View file @
e3fa9e4d
...
...
@@ -77,7 +77,7 @@ class InteractionVolumeAssemblerBase
elemVolVarsPtr_
=
&
elemVolVars
;
}
//
!
return functions to the local views & problem
// return functions to the local views & problem
const
Problem
&
problem
()
const
{
return
*
problemPtr_
;
}
const
FVElementGeometry
&
fvGeometry
()
const
{
return
*
fvGeometryPtr_
;
}
const
ElementVolumeVariables
&
elemVolVars
()
const
{
return
*
elemVolVarsPtr_
;
}
...
...
dumux/discretization/cellcentered/mpfa/localfacedata.hh
View file @
e3fa9e4d
...
...
@@ -67,7 +67,7 @@ public:
,
isOutside_
(
true
)
{}
//
!
Functions to return stored data
// Functions to return stored data
LocalIndexType
ivLocalScvfIndex
()
const
{
return
ivLocalScvfIndex_
;
}
LocalIndexType
ivLocalInsideScvIndex
()
const
{
return
ivLocalInsideScvIndex_
;
}
LocalIndexType
scvfLocalOutsideScvfIndex
()
const
{
assert
(
isOutside_
);
return
scvfLocalOutsideScvfIndex_
;
}
...
...
dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
View file @
e3fa9e4d
...
...
@@ -19,7 +19,7 @@
/*!
* \file
* \ingroup CCMpfaDiscretization
* \brief Class
es
for the interaction volume of the mpfa-o scheme.
* \brief Class for the interaction volume of the mpfa-o scheme.
*/
#ifndef DUMUX_DISCRETIZATION_CC_MPFA_O_INTERACTIONVOLUME_HH
#define DUMUX_DISCRETIZATION_CC_MPFA_O_INTERACTIONVOLUME_HH
...
...
@@ -135,7 +135,6 @@ public:
//! publicly state the mpfa-scheme this interaction volume is associated with
static
constexpr
MpfaMethods
MpfaMethod
=
MpfaMethods
::
oMethod
;
//! Types required in the assembly of the local eq system
using
Problem
=
typename
GET_PROP_TYPE
(
TypeTag
,
Problem
);
using
FVElementGeometry
=
typename
GET_PROP_TYPE
(
TypeTag
,
FVElementGeometry
);
using
ElementVolumeVariables
=
typename
GET_PROP_TYPE
(
TypeTag
,
ElementVolumeVariables
);
...
...
@@ -274,16 +273,19 @@ public:
//! returns a reference to the information container on Dirichlet BCs within this iv
const
std
::
vector
<
DirichletData
>&
dirichletData
()
const
{
return
dirichletData_
;
}
//! return
functions for matrices involved in local system of
equation
s
//! return
s the matrix associated with face unknowns in local
equation
system
const
Matrix
&
A
()
const
{
return
A_
;
}
Matrix
&
A
()
{
return
A_
;
}
//! returns the matrix associated with cell unknowns in local equation system
const
Matrix
&
B
()
const
{
return
B_
;
}
Matrix
&
B
()
{
return
B_
;
}
//! returns the matrix associated with face unknowns in flux expressions
const
Matrix
&
C
()
const
{
return
C_
;
}
Matrix
&
C
()
{
return
C_
;
}
//! returns container storing the transmissibilities for each face & coordinate
const
std
::
vector
<
std
::
vector
<
DimVector
>
>&
omegas
()
const
{
return
wijk_
;
}
std
::
vector
<
std
::
vector
<
DimVector
>
>&
omegas
()
{
return
wijk_
;
}
...
...
dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh
View file @
e3fa9e4d
...
...
@@ -74,13 +74,13 @@ public:
template
<
class
GetTensorFunction
>
void
assemble
(
Matrix
&
T
,
InteractionVolume
&
iv
,
const
GetTensorFunction
&
getTensor
)
{
//
!
assemble D into T directly
// assemble D into T directly
assembleLocalMatrices_
(
iv
.
A
(),
iv
.
B
(),
iv
.
C
(),
T
,
iv
,
getTensor
);
//
!
maybe solve the local system
// maybe solve the local system
if
(
iv
.
numUnknowns
()
>
0
)
{
//
!
T = C*A^-1*B + D
// T = C*A^-1*B + D
iv
.
A
().
invert
();
iv
.
C
().
rightmultiply
(
iv
.
A
());
T
+=
multiplyMatrices
(
iv
.
C
(),
iv
.
B
());
...
...
@@ -100,10 +100,10 @@ public:
template
<
class
OutsideTijContainer
,
class
GetTensorFunction
>
void
assemble
(
OutsideTijContainer
&
outsideTij
,
Matrix
&
T
,
InteractionVolume
&
iv
,
const
GetTensorFunction
&
getTensor
)
{
//
!
assemble D into T directly
// assemble D into T directly
assembleLocalMatrices_
(
iv
.
A
(),
iv
.
B
(),
iv
.
C
(),
T
,
iv
,
getTensor
);
//
!
maybe solve the local system
// maybe solve the local system
if
(
iv
.
numUnknowns
()
>
0
)
{
// T = C*A^-1*B + D
...
...
@@ -111,7 +111,7 @@ public:
iv
.
B
().
leftmultiply
(
iv
.
A
());
T
+=
multiplyMatrices
(
iv
.
C
(),
iv
.
B
());
//
!
compute outside transmissibilities
// compute outside transmissibilities
for
(
const
auto
&
localFaceData
:
iv
.
localFaceData
())
{
// continue only for "outside" faces
...
...
@@ -170,19 +170,19 @@ public:
InteractionVolume
&
iv
,
const
GetTensorFunction
&
getTensor
)
{
//
!
assemble D into T & C into CA directly
// assemble D into T & C into CA directly
assembleLocalMatrices_
(
iv
.
A
(),
iv
.
B
(),
CA
,
T
,
iv
,
getTensor
);
//
!
maybe solve the local system
// maybe solve the local system
if
(
iv
.
numUnknowns
()
>
0
)
{
//
!
T = C*A^-1*B + D
// T = C*A^-1*B + D
iv
.
A
().
invert
();
CA
.
rightmultiply
(
iv
.
A
());
T
+=
multiplyMatrices
(
CA
,
iv
.
B
());
}
//
!
assemble gravitational acceleration container (enforce usage of mpfa-o type version)
// assemble gravitational acceleration container (enforce usage of mpfa-o type version)
assembleGravity
(
g
,
iv
,
CA
,
getTensor
);
}
...
...
@@ -215,10 +215,10 @@ public:
InteractionVolume
&
iv
,
const
GetTensorFunction
&
getTensor
)
{
//
!
assemble D into T directly
// assemble D into T directly
assembleLocalMatrices_
(
iv
.
A
(),
iv
.
B
(),
iv
.
C
(),
T
,
iv
,
getTensor
);
//
!
maybe solve the local system
// maybe solve the local system
if
(
iv
.
numUnknowns
()
>
0
)
{
// T = C*A^-1*B + D
...
...
@@ -228,7 +228,7 @@ public:
A
=
iv
.
A
();
CA
=
iv
.
C
().
rightmultiply
(
A
);
//
!
compute outside transmissibilities
// compute outside transmissibilities
for
(
const
auto
&
localFaceData
:
iv
.
localFaceData
())
{
// continue only for "outside" faces
...
...
@@ -281,15 +281,15 @@ public:
template
<
class
GetU
>
void
assemble
(
Vector
&
u
,
const
InteractionVolume
&
iv
,
const
GetU
&
getU
)
{
//
!
resize given container
// resize given container
u
.
resize
(
iv
.
numKnowns
());
//
!
put the cell pressures first
// put the cell pressures first
for
(
LocalIndexType
i
=
0
;
i
<
iv
.
numScvs
();
++
i
)
u
[
i
]
=
getU
(
iv
.
localScv
(
i
).
globalScvIndex
()
);
//
!
Dirichlet BCs come afterwards
// Dirichlet BCs come afterwards
unsigned
int
i
=
iv
.
numScvs
();
for
(
const
auto
&
data
:
iv
.
dirichletData
())
u
[
i
++
]
=
getU
(
data
.
volVarIndex
()
);
...
...
@@ -315,15 +315,15 @@ public:
const
Matrix
&
CA
,
const
GetTensorFunction
&
getTensor
)
{
//
!
we require the CA matrix and the g vector to have the correct size already
// we require the CA matrix and the g vector to have the correct size already
assert
(
g
.
size
()
==
numPhases
&&
"Provided gravity container does not have NumPhases entries"
);
assert
(
g
[
0
].
size
()
==
iv
.
numFaces
()
&&
"Gravitation vector g does not have the correct size"
);
assert
(
CA
.
rows
()
==
iv
.
numFaces
()
&&
CA
.
cols
()
==
iv
.
numUnknowns
()
&&
"Matrix CA does not have the correct size"
);
//! For each face, we...
//! - arithmetically average the phase densities
//! - compute the term alpha:= A
*
rho
*n^T*K*g
in each neighboring cell
//! - compute
sum_
alpha
s
= alpha_outside - alpha_inside
//! - compute the term
\f$ \
alpha
:= A
\
rho
\ \mathbf{n}^T \mathbf{K} \mathbf{g} \f$
in each neighboring cell
//! - compute
\f$ \
alpha
^*
=
\
alpha_
{
outside
}
-
\
alpha_
{
inside
} \f$
using
Scalar
=
typename
Vector
::
value_type
;
std
::
array
<
Vector
,
numPhases
>
sum_alphas
;
...
...
@@ -335,12 +335,12 @@ public:
for
(
LocalIndexType
faceIdx
=
0
;
faceIdx
<
iv
.
numFaces
();
++
faceIdx
)
{
//
!
gravitational acceleration on this face
// gravitational acceleration on this face
const
auto
&
curLocalScvf
=
iv
.
localScvf
(
faceIdx
);
const
auto
&
curGlobalScvf
=
this
->
fvGeometry
().
scvf
(
curLocalScvf
.
globalScvfIndex
());
const
auto
gravity
=
this
->
problem
().
gravityAtPos
(
curGlobalScvf
.
ipGlobal
());
//
!
get permeability tensor in "positive" sub volume
// get permeability tensor in "positive" sub volume
const
auto
&
neighborScvIndices
=
curLocalScvf
.
neighboringLocalScvIndices
();
const
auto
&
posLocalScv
=
iv
.
localScv
(
neighborScvIndices
[
0
]);
const
auto
&
posGlobalScv
=
this
->
fvGeometry
().
scv
(
posLocalScv
.
globalScvIndex
());
...
...
@@ -348,8 +348,8 @@ public:
const
auto
&
posElement
=
iv
.
element
(
neighborScvIndices
[
0
]);
const
auto
tensor
=
getTensor
(
this
->
problem
(),
posElement
,
posVolVars
,
this
->
fvGeometry
(),
posGlobalScv
);
//
!
This function should never be called for surface grids,
//
!
for which there is the specialization of this function below
// This function should never be called for surface grids,
// for which there is the specialization of this function below
assert
(
neighborScvIndices
.
size
()
<=
2
&&
"Scvf seems to have more than one outside scv!"
);
std
::
array
<
Scalar
,
numPhases
>
rho
;
...
...
@@ -361,7 +361,7 @@ public:
if
(
!
curGlobalScvf
.
boundary
())
{
//
!
obtain outside tensor
// obtain outside tensor
const
auto
&
negLocalScv
=
iv
.
localScv
(
neighborScvIndices
[
1
]
);
const
auto
&
negGlobalScv
=
this
->
fvGeometry
().
scv
(
negLocalScv
.
globalScvIndex
());
const
auto
&
negVolVars
=
this
->
elemVolVars
()[
negGlobalScv
];
...
...
@@ -386,7 +386,7 @@ public:
sum_alphas
[
pIdx
][
localDofIdx
]
-=
alpha_inside
*
rho
[
pIdx
]
*
curGlobalScvf
.
area
();
}
}
//
!
use Dirichlet BC densities
// use Dirichlet BC densities
else
{
const
auto
&
dirichletVolVars
=
this
->
elemVolVars
()[
curGlobalScvf
.
outsideScvIdx
()];
...
...
@@ -394,12 +394,12 @@ public:
rho
[
pIdx
]
=
dirichletVolVars
.
density
(
pIdx
);
}
//
!
add "inside" alpha to gravity container
// add "inside" alpha to gravity container
for
(
unsigned
int
pIdx
=
0
;
pIdx
<
numPhases
;
++
pIdx
)
g
[
pIdx
][
faceIdx
]
+=
alpha_inside
*
rho
[
pIdx
]
*
curGlobalScvf
.
area
();
}
//
!
g += CA*sum_alphas
// g += CA*sum_alphas
for
(
unsigned
int
pIdx
=
0
;
pIdx
<
numPhases
;
++
pIdx
)
CA
.
umv
(
sum_alphas
[
pIdx
],
g
[
pIdx
]);
}
...
...
@@ -431,7 +431,7 @@ public:
const
Matrix
&
A
,
const
GetTensorFunction
&
getTensor
)
{
//
!
we require the CA matrix and the gravity containers to have the correct size already
// we require the CA matrix and the gravity containers to have the correct size already
assert
(
CA
.
rows
()
==
iv
.
numFaces
()
&&
CA
.
cols
()
==
iv
.
numUnknowns
()
&&
"Matrix CA does not have the correct size"
);
assert
(
g
.
size
()
==
numPhases
&&
"Provided gravity container does not have NumPhases entries"
);
assert
(
outsideG
.
size
()
==
numPhases
&&
"Provided outside gravity container does not have NumPhases entries"
);
...
...
@@ -442,8 +442,8 @@ public:
//! For each face, we...
//! - arithmetically average the phase densities
//! - compute the term alpha:=
A*rho*n^T*K*g
in each neighboring cell
//! - compute
sum_
alpha
s
= alpha_outside - alpha_inside
//! - compute the term
\f$ \
alpha
:=
\mathbf{A} \rho \ \mathbf{n}^T \mathbf{K} \mathbf{g} \f$
in each neighboring cell
//! - compute
\f$ \
alpha
^*
=
\sum{\
alpha_
{
outside
, i}}
-
\
alpha_
{
inside
} \f$
using
Scalar
=
typename
Vector
::
value_type
;
// reset everything to zero
...
...
@@ -457,12 +457,12 @@ public:
for
(
LocalIndexType
faceIdx
=
0
;
faceIdx
<
iv
.
numFaces
();
++
faceIdx
)
{
//
!
gravitational acceleration on this face
// gravitational acceleration on this face
const
auto
&
curLocalScvf
=
iv
.
localScvf
(
faceIdx
);
const
auto
&
curGlobalScvf
=
this
->
fvGeometry
().
scvf
(
curLocalScvf
.
globalScvfIndex
());
const
auto
gravity
=
this
->
problem
().
gravityAtPos
(
curGlobalScvf
.
ipGlobal
());
//
!
get permeability tensor in "positive" sub volume
// get permeability tensor in "positive" sub volume
const
auto
&
neighborScvIndices
=
curLocalScvf
.
neighboringLocalScvIndices
();
const
auto
&
posLocalScv
=
iv
.
localScv
(
neighborScvIndices
[
0
]);
const
auto
&
posGlobalScv
=
this
->
fvGeometry
().
scv
(
posLocalScv
.
globalScvIndex
());
...
...
@@ -486,7 +486,7 @@ public:
{
for
(
unsigned
int
idxInOutside
=
0
;
idxInOutside
<
curGlobalScvf
.
numOutsideScvs
();
++
idxInOutside
)
{
//
!
obtain outside tensor
// obtain outside tensor
const
auto
&
negLocalScv
=
iv
.
localScv
(
neighborScvIndices
[
idxInOutside
]
);
const
auto
&
negGlobalScv
=
this
->
fvGeometry
().
scv
(
negLocalScv
.
globalScvIndex
());
const
auto
&
negVolVars
=
this
->
elemVolVars
()[
negGlobalScv
];
...
...
@@ -512,7 +512,7 @@ public:
sum_alphas
[
pIdx
][
localDofIdx
]
*=
rho
[
pIdx
]
*
curGlobalScvf
.
area
();
}
}
//
!
use Dirichlet BC densities
// use Dirichlet BC densities
else
{
const
auto
&
dirichletVolVars
=
this
->
elemVolVars
()[
curGlobalScvf
.
outsideScvIdx
()];
...
...
@@ -520,7 +520,7 @@ public:
rho
[
pIdx
]
=
dirichletVolVars
.
density
(
pIdx
);
}
//
!
add "inside" & "outside" alphas to gravity containers
// add "inside" & "outside" alphas to gravity containers
for
(
unsigned
int
pIdx
=
0
;
pIdx
<
numPhases
;
++
pIdx
)
{
g
[
pIdx
][
faceIdx
]
+=
alpha_inside
*
rho
[
pIdx
]
*
curGlobalScvf
.
area
();
...
...
@@ -530,8 +530,8 @@ public:
}
}
//
!
g += CA*sum_alphas
//
!
outsideG = wikj*A^-1*sum_alphas + outsideG
// g += CA*sum_alphas
// outsideG = wikj*A^-1*sum_alphas + outsideG
for
(
unsigned
int
pIdx
=
0
;
pIdx
<
numPhases
;
++
pIdx
)
{
CA
.
umv
(
sum_alphas
[
pIdx
],
g
[
pIdx
]);
...
...
@@ -539,7 +539,7 @@ public:
Vector
AG
(
iv
.
numUnknowns
());
A
.
mv
(
sum_alphas
[
pIdx
],
AG
);
//
!
compute gravitational accelerations
// compute gravitational accelerations
for
(
const
auto
&
localFaceData
:
iv
.
localFaceData
())
{
// continue only for "outside" faces
...
...
@@ -573,10 +573,9 @@ private:
* in the interaction volume-local system of equations resulting from flux
* and solution continuity across the scvfs.
*
* Flux expressions:
* \f$\mathbf{f} = \mathbf{C} \bar{matbf{u}} + \mathbf{D} matbf{u}\f$.
* Continuity equations
* \f$\mathbf{A} \, \bar{matbf{u}} = \mathbf{B} \, matbf{u}\f$,
* Flux expressions: \f$\mathbf{f} = \mathbf{C} \bar{\mathbf{u}} + \mathbf{D} \mathbf{u}\f$.
*
* Continuity equations: \f$\mathbf{A} \, \bar{\mathbf{u}} = \mathbf{B} \, \mathbf{u}\f$.
*
* \note The matrices are expected to have been resized beforehand.
*
...
...
@@ -592,31 +591,31 @@ private:
InteractionVolume
&
iv
,
const
GetTensorFunction
&
getTensor
)
{
//
!
Matrix D is assumed to have the right size already
// Matrix D is assumed to have the right size already
assert
(
D
.
rows
()
==
iv
.
numFaces
()
&&
D
.
cols
()
==
iv
.
numKnowns
()
&&
"Matrix D does not have the correct size"
);
//
!
if only Dirichlet faces are present in the iv,
//
!
the matrices A, B & C are undefined and D = T
// if only Dirichlet faces are present in the iv,
// the matrices A, B & C are undefined and D = T
if
(
iv
.
numUnknowns
()
==
0
)
{
//
!
reset matrix beforehand
// reset matrix beforehand
D
=
0.0
;
//
!
Loop over all the faces, in this case these are all dirichlet boundaries
// Loop over all the faces, in this case these are all dirichlet boundaries
for
(
LocalIndexType
faceIdx
=
0
;
faceIdx
<
iv
.
numFaces
();
++
faceIdx
)
{
const
auto
&
curLocalScvf
=
iv
.
localScvf
(
faceIdx
);
const
auto
&
curGlobalScvf
=
this
->
fvGeometry
().
scvf
(
curLocalScvf
.
globalScvfIndex
());
const
auto
&
neighborScvIndices
=
curLocalScvf
.
neighboringLocalScvIndices
();
//
!
get tensor in "positive" sub volume
// get tensor in "positive" sub volume
const
auto
&
posLocalScv
=
iv
.
localScv
(
neighborScvIndices
[
0
]);
const
auto
&
posGlobalScv
=
this
->
fvGeometry
().
scv
(
posLocalScv
.
globalScvIndex
());
const
auto
&
posVolVars
=
this
->
elemVolVars
()[
posGlobalScv
];
const
auto
&
posElement
=
iv
.
element
(
neighborScvIndices
[
0
]);
const
auto
tensor
=
getTensor
(
this
->
problem
(),
posElement
,
posVolVars
,
this
->
fvGeometry
(),
posGlobalScv
);
//
!
the omega factors of the "positive" sub volume
// the omega factors of the "positive" sub volume
const
auto
wijk
=
computeMpfaTransmissibility
(
posLocalScv
,
curGlobalScvf
,
tensor
,
posVolVars
.
extrusionFactor
());
const
auto
posScvLocalDofIdx
=
posLocalScv
.
localDofIndex
();
...
...
@@ -631,12 +630,12 @@ private:
}
else
{
//
!
we require the matrices A,B,C to have the correct size already
// we require the matrices A,B,C to have the correct size already
assert
(
A
.
rows
()
==
iv
.
numUnknowns
()
&&
A
.
cols
()
==
iv
.
numUnknowns
()
&&
"Matrix A does not have the correct size"
);
assert
(
B
.
rows
()
==
iv
.
numUnknowns
()
&&
B
.
cols
()
==
iv
.
numKnowns
()
&&
"Matrix B does not have the correct size"
);
assert
(
C
.
rows
()
==
iv
.
numFaces
()
&&
C
.
cols
()
==
iv
.
numKnowns
()
&&
"Matrix C does not have the correct size"
);
//
!
reset matrices
// reset matrices
A
=
0.0
;
B
=
0.0
;
C
=
0.0
;
...
...
@@ -650,7 +649,7 @@ private:
const
auto
curIsDirichlet
=
curLocalScvf
.
isDirichlet
();
const
auto
curLocalDofIdx
=
curLocalScvf
.
localDofIndex
();
//
!
get tensor in "positive" sub volume
// get tensor in "positive" sub volume
const
auto
&
neighborScvIndices
=
curLocalScvf
.
neighboringLocalScvIndices
();
const
auto
&
posLocalScv
=
iv
.
localScv
(
neighborScvIndices
[
0
]);
const
auto
&
posGlobalScv
=
this
->
fvGeometry
().
scv
(
posLocalScv
.
globalScvIndex
());
...
...
@@ -658,24 +657,24 @@ private:
const
auto
&
posElement
=
iv
.
element
(
neighborScvIndices
[
0
]);
const
auto
tensor
=
getTensor
(
this
->
problem
(),
posElement
,
posVolVars
,
this
->
fvGeometry
(),
posGlobalScv
);
//
!
the omega factors of the "positive" sub volume
// the omega factors of the "positive" sub volume
wijk
[
faceIdx
][
0
]
=
computeMpfaTransmissibility
(
posLocalScv
,
curGlobalScvf
,
tensor
,
posVolVars
.
extrusionFactor
());
//
!
go over the coordinate directions in the positive sub volume
// go over the coordinate directions in the positive sub volume
for
(
unsigned
int
localDir
=
0
;
localDir
<
dim
;
localDir
++
)
{
const
auto
&
otherLocalScvf
=
iv
.
localScvf
(
posLocalScv
.
scvfIdxLocal
(
localDir
)
);
const
auto
otherLocalDofIdx
=
otherLocalScvf
.
localDofIndex
();
//
!
if we are not on a Dirichlet face, add entries associated with unknown face pressures
//
!
i.e. in matrix C and maybe A (if current face is not a Dirichlet face)
// if we are not on a Dirichlet face, add entries associated with unknown face pressures
// i.e. in matrix C and maybe A (if current face is not a Dirichlet face)
if
(
!
otherLocalScvf
.
isDirichlet
())
{
C
[
faceIdx
][
otherLocalDofIdx
]
-=
wijk
[
faceIdx
][
0
][
localDir
];
if
(
!
curIsDirichlet
)
A
[
curLocalDofIdx
][
otherLocalDofIdx
]
-=
wijk
[
faceIdx
][
0
][
localDir
];
}
//
!
the current face is a Dirichlet face and creates entries in D & maybe B
// the current face is a Dirichlet face and creates entries in D & maybe B
else
{
D
[
faceIdx
][
otherLocalDofIdx
]
-=
wijk
[
faceIdx
][
0
][
localDir
];
...
...
@@ -683,7 +682,7 @@ private:
B
[
curLocalDofIdx
][
otherLocalDofIdx
]
+=
wijk
[
faceIdx
][
0
][
localDir
];
}
//
!
add entries related to pressures at the scv centers (dofs)
// add entries related to pressures at the scv centers (dofs)
const
auto
posScvLocalDofIdx
=
posLocalScv
.
localDofIndex
();
D
[
faceIdx
][
posScvLocalDofIdx
]
+=
wijk
[
faceIdx
][
0
][
localDir
];
...
...
dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh
View file @
e3fa9e4d
...
...
@@ -44,7 +44,7 @@ class CCMpfaOInteractionVolumeLocalScv
{
public:
//
! publicly state
some types
//
export
some types
using
GridIndexType
=
typename
IvIndexSet
::
GridIndexType
;
using
LocalIndexType
=
typename
IvIndexSet
::
LocalIndexType
;
using
GlobalCoordinate
=
GC
;
...
...
@@ -134,7 +134,7 @@ struct CCMpfaOInteractionVolumeLocalScvf
using
LocalIndexContainer
=
typename
IvIndexSet
::
LocalIndexContainer
;