Skip to content
GitLab
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
15129c16
Commit
15129c16
authored
Jul 19, 2018
by
Thomas Fetzer
Browse files
[rans] Fix ID vs. Idx confusion
parent
cbc39ecb
Changes
15
Hide whitespace changes
Inline
Side-by-side
dumux/freeflow/rans/oneeq/problem.hh
View file @
15129c16
...
...
@@ -96,7 +96,7 @@ public:
for
(
const
auto
&
element
:
elements
(
this
->
fvGridGeometry
().
gridView
()))
{
unsigned
int
elementI
D
=
this
->
fvGridGeometry
().
elementMapper
().
index
(
element
);
unsigned
int
elementI
dx
=
this
->
fvGridGeometry
().
elementMapper
().
index
(
element
);
auto
fvGeometry
=
localView
(
this
->
fvGridGeometry
());
fvGeometry
.
bindElement
(
element
);
...
...
@@ -107,26 +107,26 @@ public:
PrimaryVariables
priVars
=
makePriVarsFromCellCenterPriVars
<
PrimaryVariables
>
(
cellCenterPriVars
);
auto
elemSol
=
elementSolution
<
typename
FVGridGeometry
::
LocalView
>
(
std
::
move
(
priVars
));
// NOTE: first update the turbulence quantities
storedViscosityTilde_
[
elementI
D
]
=
elemSol
[
0
][
Indices
::
viscosityTildeIdx
];
storedViscosityTilde_
[
elementI
dx
]
=
elemSol
[
0
][
Indices
::
viscosityTildeIdx
];
// NOTE: then update the volVars
VolumeVariables
volVars
;
volVars
.
update
(
elemSol
,
asImp_
(),
element
,
scv
);
storedDynamicEddyViscosity_
[
elementI
D
]
=
volVars
.
calculateEddyViscosity
();
storedDynamicEddyViscosity_
[
elementI
dx
]
=
volVars
.
calculateEddyViscosity
();
}
}
// calculate cell-center-averaged velocity gradients, maximum, and minimum values
for
(
const
auto
&
element
:
elements
(
this
->
fvGridGeometry
().
gridView
()))
{
unsigned
int
elementI
D
=
this
->
fvGridGeometry
().
elementMapper
().
index
(
element
);
unsigned
int
elementI
dx
=
this
->
fvGridGeometry
().
elementMapper
().
index
(
element
);
for
(
unsigned
int
dimIdx
=
0
;
dimIdx
<
Grid
::
dimension
;
++
dimIdx
)
{
storedViscosityTildeGradient_
[
elementI
D
][
dimIdx
]
=
(
storedViscosityTilde_
[
ParentType
::
neighborI
D
_
[
elementI
D
][
dimIdx
][
1
]]
-
storedViscosityTilde_
[
ParentType
::
neighborI
D
_
[
elementI
D
][
dimIdx
][
0
]])
/
(
ParentType
::
cellCenter_
[
ParentType
::
neighborI
D
_
[
elementI
D
][
dimIdx
][
1
]][
dimIdx
]
-
ParentType
::
cellCenter_
[
ParentType
::
neighborI
D
_
[
elementI
D
][
dimIdx
][
0
]][
dimIdx
]);
storedViscosityTildeGradient_
[
elementI
dx
][
dimIdx
]
=
(
storedViscosityTilde_
[
ParentType
::
neighborI
dx
_
[
elementI
dx
][
dimIdx
][
1
]]
-
storedViscosityTilde_
[
ParentType
::
neighborI
dx
_
[
elementI
dx
][
dimIdx
][
0
]])
/
(
ParentType
::
cellCenter_
[
ParentType
::
neighborI
dx
_
[
elementI
dx
][
dimIdx
][
1
]][
dimIdx
]
-
ParentType
::
cellCenter_
[
ParentType
::
neighborI
dx
_
[
elementI
dx
][
dimIdx
][
0
]][
dimIdx
]);
}
auto
fvGeometry
=
localView
(
this
->
fvGridGeometry
());
...
...
@@ -139,13 +139,13 @@ public:
// face Value
Scalar
dirichletViscosityTilde
=
asImp_
().
dirichlet
(
element
,
scvf
)[
Indices
::
viscosityTildeIdx
];
unsigned
int
neighborI
D
=
ParentType
::
neighborI
D
_
[
elementI
D
][
normDim
][
0
];
if
(
scvf
.
center
()[
normDim
]
<
ParentType
::
cellCenter_
[
elementI
D
][
normDim
])
neighborI
D
=
ParentType
::
neighborI
D
_
[
elementI
D
][
normDim
][
1
];
unsigned
int
neighborI
dx
=
ParentType
::
neighborI
dx
_
[
elementI
dx
][
normDim
][
0
];
if
(
scvf
.
center
()[
normDim
]
<
ParentType
::
cellCenter_
[
elementI
dx
][
normDim
])
neighborI
dx
=
ParentType
::
neighborI
dx
_
[
elementI
dx
][
normDim
][
1
];
storedViscosityTildeGradient_
[
elementI
D
][
normDim
]
=
(
storedViscosityTilde_
[
neighborI
D
]
-
dirichletViscosityTilde
)
/
(
ParentType
::
cellCenter_
[
neighborI
D
][
normDim
]
-
scvf
.
center
()[
normDim
]);
storedViscosityTildeGradient_
[
elementI
dx
][
normDim
]
=
(
storedViscosityTilde_
[
neighborI
dx
]
-
dirichletViscosityTilde
)
/
(
ParentType
::
cellCenter_
[
neighborI
dx
][
normDim
]
-
scvf
.
center
()[
normDim
]);
}
}
}
...
...
dumux/freeflow/rans/oneeq/volumevariables.hh
View file @
15129c16
...
...
@@ -86,12 +86,12 @@ public:
{
RANSParentType
::
updateRANSProperties
(
elemSol
,
problem
,
element
,
scv
);
viscosityTilde_
=
elemSol
[
0
][
Indices
::
viscosityTildeIdx
];
storedViscosityTilde_
=
problem
.
storedViscosityTilde_
[
RANSParentType
::
elementI
D
()];
storedViscosityTildeGradient_
=
problem
.
storedViscosityTildeGradient_
[
RANSParentType
::
elementI
D
()];
stressTensorScalarProduct_
=
problem
.
stressTensorScalarProduct_
[
RANSParentType
::
elementI
D
()];
vorticityTensorScalarProduct_
=
problem
.
vorticityTensorScalarProduct_
[
RANSParentType
::
elementI
D
()];
storedViscosityTilde_
=
problem
.
storedViscosityTilde_
[
RANSParentType
::
elementI
dx
()];
storedViscosityTildeGradient_
=
problem
.
storedViscosityTildeGradient_
[
RANSParentType
::
elementI
dx
()];
stressTensorScalarProduct_
=
problem
.
stressTensorScalarProduct_
[
RANSParentType
::
elementI
dx
()];
vorticityTensorScalarProduct_
=
problem
.
vorticityTensorScalarProduct_
[
RANSParentType
::
elementI
dx
()];
if
(
problem
.
useStoredEddyViscosity_
)
RANSParentType
::
setDynamicEddyViscosity_
(
problem
.
storedDynamicEddyViscosity_
[
RANSParentType
::
elementI
D
()]);
RANSParentType
::
setDynamicEddyViscosity_
(
problem
.
storedDynamicEddyViscosity_
[
RANSParentType
::
elementI
dx
()]);
else
RANSParentType
::
setDynamicEddyViscosity_
(
calculateEddyViscosity
());
RANSParentType
::
calculateEddyDiffusivity
(
problem
);
...
...
dumux/freeflow/rans/problem.hh
View file @
15129c16
...
...
@@ -92,9 +92,9 @@ public:
calledUpdateStaticWallProperties
=
true
;
// update size and initial values of the global vectors
wallElementI
D
_
.
resize
(
this
->
fvGridGeometry
().
elementMapper
().
size
());
wallElementI
dx
_
.
resize
(
this
->
fvGridGeometry
().
elementMapper
().
size
());
wallDistance_
.
resize
(
this
->
fvGridGeometry
().
elementMapper
().
size
(),
std
::
numeric_limits
<
Scalar
>::
max
());
neighborI
D
_
.
resize
(
this
->
fvGridGeometry
().
elementMapper
().
size
());
neighborI
dx
_
.
resize
(
this
->
fvGridGeometry
().
elementMapper
().
size
());
cellCenter_
.
resize
(
this
->
fvGridGeometry
().
elementMapper
().
size
(),
GlobalPosition
(
0.0
));
velocity_
.
resize
(
this
->
fvGridGeometry
().
elementMapper
().
size
(),
DimVector
(
0.0
));
velocityMaximum_
.
resize
(
this
->
fvGridGeometry
().
elementMapper
().
size
(),
DimVector
(
0.0
));
...
...
@@ -141,8 +141,8 @@ public:
// search for shortest distance to wall for each element
for
(
const
auto
&
element
:
elements
(
gridView
))
{
unsigned
int
elementI
D
=
this
->
fvGridGeometry
().
elementMapper
().
index
(
element
);
cellCenter_
[
elementI
D
]
=
element
.
geometry
().
center
();
unsigned
int
elementI
dx
=
this
->
fvGridGeometry
().
elementMapper
().
index
(
element
);
cellCenter_
[
elementI
dx
]
=
element
.
geometry
().
center
();
for
(
unsigned
int
i
=
0
;
i
<
wallPositions
.
size
();
++
i
)
{
static
const
int
problemWallNormalAxis
...
...
@@ -158,26 +158,26 @@ public:
GlobalPosition
global
=
element
.
geometry
().
center
();
global
-=
wallPositions
[
i
];
// second and argument ensures to use only aligned elements
if
(
abs
(
global
[
searchAxis
])
<
wallDistance_
[
elementI
D
]
if
(
abs
(
global
[
searchAxis
])
<
wallDistance_
[
elementI
dx
]
&&
abs
(
global
[
searchAxis
])
<
global
.
two_norm
()
+
1e-8
&&
abs
(
global
[
searchAxis
])
>
global
.
two_norm
()
-
1e-8
)
{
wallDistance_
[
elementI
D
]
=
abs
(
global
[
searchAxis
]);
wallElementI
D
_
[
elementI
D
]
=
wallElements
[
i
];
wallNormalAxis_
[
elementI
D
]
=
searchAxis
;
sandGrainRoughness_
[
elementI
D
]
=
asImp_
().
sandGrainRoughnessAtPos
(
wallPositions
[
i
]);
wallDistance_
[
elementI
dx
]
=
abs
(
global
[
searchAxis
]);
wallElementI
dx
_
[
elementI
dx
]
=
wallElements
[
i
];
wallNormalAxis_
[
elementI
dx
]
=
searchAxis
;
sandGrainRoughness_
[
elementI
dx
]
=
asImp_
().
sandGrainRoughnessAtPos
(
wallPositions
[
i
]);
}
}
}
// search for neighbor I
D
s
// search for neighbor I
dx
s
for
(
const
auto
&
element
:
elements
(
gridView
))
{
unsigned
int
elementI
D
=
this
->
fvGridGeometry
().
elementMapper
().
index
(
element
);
unsigned
int
elementI
dx
=
this
->
fvGridGeometry
().
elementMapper
().
index
(
element
);
for
(
unsigned
int
dimIdx
=
0
;
dimIdx
<
dim
;
++
dimIdx
)
{
neighborI
D
_
[
elementI
D
][
dimIdx
][
0
]
=
elementI
D
;
neighborI
D
_
[
elementI
D
][
dimIdx
][
1
]
=
elementI
D
;
neighborI
dx
_
[
elementI
dx
][
dimIdx
][
0
]
=
elementI
dx
;
neighborI
dx
_
[
elementI
dx
][
dimIdx
][
1
]
=
elementI
dx
;
}
for
(
const
auto
&
intersection
:
intersections
(
gridView
,
element
))
...
...
@@ -185,18 +185,18 @@ public:
if
(
intersection
.
boundary
())
continue
;
unsigned
int
neighborI
D
=
this
->
fvGridGeometry
().
elementMapper
().
index
(
intersection
.
outside
());
unsigned
int
neighborI
dx
=
this
->
fvGridGeometry
().
elementMapper
().
index
(
intersection
.
outside
());
for
(
unsigned
int
dimIdx
=
0
;
dimIdx
<
dim
;
++
dimIdx
)
{
if
(
abs
(
cellCenter_
[
elementI
D
][
dimIdx
]
-
cellCenter_
[
neighborI
D
][
dimIdx
])
>
1e-8
)
if
(
abs
(
cellCenter_
[
elementI
dx
][
dimIdx
]
-
cellCenter_
[
neighborI
dx
][
dimIdx
])
>
1e-8
)
{
if
(
cellCenter_
[
elementI
D
][
dimIdx
]
>
cellCenter_
[
neighborI
D
][
dimIdx
])
if
(
cellCenter_
[
elementI
dx
][
dimIdx
]
>
cellCenter_
[
neighborI
dx
][
dimIdx
])
{
neighborI
D
_
[
elementI
D
][
dimIdx
][
0
]
=
neighborI
D
;
neighborI
dx
_
[
elementI
dx
][
dimIdx
][
0
]
=
neighborI
dx
;
}
if
(
cellCenter_
[
elementI
D
][
dimIdx
]
<
cellCenter_
[
neighborI
D
][
dimIdx
])
if
(
cellCenter_
[
elementI
dx
][
dimIdx
]
<
cellCenter_
[
neighborI
dx
][
dimIdx
])
{
neighborI
D
_
[
elementI
D
][
dimIdx
][
1
]
=
neighborI
D
;
neighborI
dx
_
[
elementI
dx
][
dimIdx
][
1
]
=
neighborI
dx
;
}
}
}
...
...
@@ -235,7 +235,7 @@ public:
{
auto
fvGeometry
=
localView
(
this
->
fvGridGeometry
());
fvGeometry
.
bindElement
(
element
);
unsigned
int
elementI
D
=
this
->
fvGridGeometry
().
elementMapper
().
index
(
element
);
unsigned
int
elementI
dx
=
this
->
fvGridGeometry
().
elementMapper
().
index
(
element
);
// calculate velocities
DimVector
velocityTemp
(
0.0
);
...
...
@@ -246,44 +246,44 @@ public:
velocityTemp
[
scvf
.
directionIndex
()]
+=
numericalSolutionFace
;
}
for
(
unsigned
int
dimIdx
=
0
;
dimIdx
<
dim
;
++
dimIdx
)
velocity_
[
elementI
D
][
dimIdx
]
=
velocityTemp
[
dimIdx
]
*
0.5
;
// faces are equidistant to cell center
velocity_
[
elementI
dx
][
dimIdx
]
=
velocityTemp
[
dimIdx
]
*
0.5
;
// faces are equidistant to cell center
}
// calculate cell-center-averaged velocity gradients, maximum, and minimum values
for
(
const
auto
&
element
:
elements
(
this
->
fvGridGeometry
().
gridView
()))
{
unsigned
int
elementI
D
=
this
->
fvGridGeometry
().
elementMapper
().
index
(
element
);
unsigned
int
wallElementI
D
=
wallElementI
D
_
[
elementI
D
];
unsigned
int
elementI
dx
=
this
->
fvGridGeometry
().
elementMapper
().
index
(
element
);
unsigned
int
wallElementI
dx
=
wallElementI
dx
_
[
elementI
dx
];
Scalar
maxVelocity
=
0.0
;
for
(
unsigned
int
dimIdx
=
0
;
dimIdx
<
dim
;
++
dimIdx
)
{
for
(
unsigned
int
velIdx
=
0
;
velIdx
<
dim
;
++
velIdx
)
{
velocityGradients_
[
elementI
D
][
velIdx
][
dimIdx
]
=
(
velocity_
[
neighborI
D
_
[
elementI
D
][
dimIdx
][
1
]][
velIdx
]
-
velocity_
[
neighborI
D
_
[
elementI
D
][
dimIdx
][
0
]][
velIdx
])
/
(
cellCenter_
[
neighborI
D
_
[
elementI
D
][
dimIdx
][
1
]][
dimIdx
]
-
cellCenter_
[
neighborI
D
_
[
elementI
D
][
dimIdx
][
0
]][
dimIdx
]);
velocityGradients_
[
elementI
dx
][
velIdx
][
dimIdx
]
=
(
velocity_
[
neighborI
dx
_
[
elementI
dx
][
dimIdx
][
1
]][
velIdx
]
-
velocity_
[
neighborI
dx
_
[
elementI
dx
][
dimIdx
][
0
]][
velIdx
])
/
(
cellCenter_
[
neighborI
dx
_
[
elementI
dx
][
dimIdx
][
1
]][
dimIdx
]
-
cellCenter_
[
neighborI
dx
_
[
elementI
dx
][
dimIdx
][
0
]][
dimIdx
]);
}
if
(
abs
(
velocity_
[
elementI
D
][
dimIdx
])
>
abs
(
velocityMaximum_
[
wallElementI
D
][
dimIdx
]))
if
(
abs
(
velocity_
[
elementI
dx
][
dimIdx
])
>
abs
(
velocityMaximum_
[
wallElementI
dx
][
dimIdx
]))
{
velocityMaximum_
[
wallElementI
D
][
dimIdx
]
=
velocity_
[
elementI
D
][
dimIdx
];
velocityMaximum_
[
wallElementI
dx
][
dimIdx
]
=
velocity_
[
elementI
dx
][
dimIdx
];
}
if
(
abs
(
velocity_
[
elementI
D
][
dimIdx
])
<
abs
(
velocityMinimum_
[
wallElementI
D
][
dimIdx
]))
if
(
abs
(
velocity_
[
elementI
dx
][
dimIdx
])
<
abs
(
velocityMinimum_
[
wallElementI
dx
][
dimIdx
]))
{
velocityMinimum_
[
wallElementI
D
][
dimIdx
]
=
velocity_
[
elementI
D
][
dimIdx
];
velocityMinimum_
[
wallElementI
dx
][
dimIdx
]
=
velocity_
[
elementI
dx
][
dimIdx
];
}
if
(
0
<=
flowNormalAxis
&&
flowNormalAxis
<
dim
)
{
flowNormalAxis_
[
elementI
D
]
=
flowNormalAxis
;
flowNormalAxis_
[
elementI
dx
]
=
flowNormalAxis
;
}
else
if
(
abs
(
maxVelocity
)
<
abs
(
velocity_
[
elementI
D
][
dimIdx
]))
else
if
(
abs
(
maxVelocity
)
<
abs
(
velocity_
[
elementI
dx
][
dimIdx
]))
{
maxVelocity
=
abs
(
velocity_
[
elementI
D
][
dimIdx
]);
flowNormalAxis_
[
elementI
D
]
=
dimIdx
;
maxVelocity
=
abs
(
velocity_
[
elementI
dx
][
dimIdx
]);
flowNormalAxis_
[
elementI
dx
]
=
dimIdx
;
}
}
...
...
@@ -302,13 +302,13 @@ public:
Scalar
dirichletVelocity
=
asImp_
().
dirichlet
(
element
,
scvf
)[
Indices
::
velocity
(
velIdx
)];
unsigned
int
neighborI
D
=
neighborI
D
_
[
elementI
D
][
scvfNormDim
][
0
];
if
(
scvf
.
center
()[
scvfNormDim
]
<
cellCenter_
[
elementI
D
][
scvfNormDim
])
neighborI
D
=
neighborI
D
_
[
elementI
D
][
scvfNormDim
][
1
];
unsigned
int
neighborI
dx
=
neighborI
dx
_
[
elementI
dx
][
scvfNormDim
][
0
];
if
(
scvf
.
center
()[
scvfNormDim
]
<
cellCenter_
[
elementI
dx
][
scvfNormDim
])
neighborI
dx
=
neighborI
dx
_
[
elementI
dx
][
scvfNormDim
][
1
];
velocityGradients_
[
elementI
D
][
velIdx
][
scvfNormDim
]
=
(
velocity_
[
neighborI
D
][
velIdx
]
-
dirichletVelocity
)
/
(
cellCenter_
[
neighborI
D
][
scvfNormDim
]
-
scvf
.
center
()[
scvfNormDim
]);
velocityGradients_
[
elementI
dx
][
velIdx
][
scvfNormDim
]
=
(
velocity_
[
neighborI
dx
][
velIdx
]
-
dirichletVelocity
)
/
(
cellCenter_
[
neighborI
dx
][
scvfNormDim
]
-
scvf
.
center
()[
scvfNormDim
]);
}
}
...
...
@@ -327,14 +327,14 @@ public:
unsigned
int
normalNormDim
=
normalFace
.
directionIndex
();
if
(
normalFace
.
boundary
()
&&
(
asImp_
().
boundaryTypes
(
element
,
normalFace
).
isBJS
(
Indices
::
velocity
(
velIdx
))))
{
unsigned
int
neighborI
D
=
neighborI
D
_
[
elementI
D
][
normalNormDim
][
0
];
if
(
normalFace
.
center
()[
normalNormDim
]
<
cellCenter_
[
elementI
D
][
normalNormDim
])
neighborI
D
=
neighborI
D
_
[
elementI
D
][
normalNormDim
][
1
];
bjsVelocityAverage
[
normalNormDim
]
+=
ParentType
::
bjsVelocity
(
scvf
,
normalFace
,
localSubFaceIdx
,
velocity_
[
elementI
D
][
velIdx
]);
if
(
bjsNumFaces
[
normalNormDim
]
>
0
&&
neighborI
D
!=
bjsNeighbor
[
normalNormDim
])
DUNE_THROW
(
Dune
::
InvalidStateException
,
"Two different neighborI
D
should not occur"
);
bjsNeighbor
[
normalNormDim
]
=
neighborI
D
;
unsigned
int
neighborI
dx
=
neighborI
dx
_
[
elementI
dx
][
normalNormDim
][
0
];
if
(
normalFace
.
center
()[
normalNormDim
]
<
cellCenter_
[
elementI
dx
][
normalNormDim
])
neighborI
dx
=
neighborI
dx
_
[
elementI
dx
][
normalNormDim
][
1
];
bjsVelocityAverage
[
normalNormDim
]
+=
ParentType
::
bjsVelocity
(
scvf
,
normalFace
,
localSubFaceIdx
,
velocity_
[
elementI
dx
][
velIdx
]);
if
(
bjsNumFaces
[
normalNormDim
]
>
0
&&
neighborI
dx
!=
bjsNeighbor
[
normalNormDim
])
DUNE_THROW
(
Dune
::
InvalidStateException
,
"Two different neighborI
dx
should not occur"
);
bjsNeighbor
[
normalNormDim
]
=
neighborI
dx
;
normalNormCoordinate
[
normalNormDim
]
=
normalFace
.
center
()[
normalNormDim
];
bjsNumFaces
[
normalNormDim
]
++
;
}
...
...
@@ -344,12 +344,12 @@ public:
if
(
bjsNumFaces
[
dirIdx
]
==
0
)
continue
;
unsigned
int
neighborI
D
=
bjsNeighbor
[
dirIdx
];
unsigned
int
neighborI
dx
=
bjsNeighbor
[
dirIdx
];
bjsVelocityAverage
[
dirIdx
]
/=
bjsNumFaces
[
dirIdx
];
velocityGradients_
[
elementI
D
][
velIdx
][
dirIdx
]
=
(
velocity_
[
neighborI
D
][
velIdx
]
-
bjsVelocityAverage
[
dirIdx
])
/
(
cellCenter_
[
neighborI
D
][
dirIdx
]
-
normalNormCoordinate
[
dirIdx
]);
velocityGradients_
[
elementI
dx
][
velIdx
][
dirIdx
]
=
(
velocity_
[
neighborI
dx
][
velIdx
]
-
bjsVelocityAverage
[
dirIdx
])
/
(
cellCenter_
[
neighborI
dx
][
dirIdx
]
-
normalNormCoordinate
[
dirIdx
]);
}
}
...
...
@@ -358,23 +358,23 @@ public:
// calculate or call all secondary variables
for
(
const
auto
&
element
:
elements
(
this
->
fvGridGeometry
().
gridView
()))
{
unsigned
int
elementI
D
=
this
->
fvGridGeometry
().
elementMapper
().
index
(
element
);
unsigned
int
elementI
dx
=
this
->
fvGridGeometry
().
elementMapper
().
index
(
element
);
Dune
::
FieldMatrix
<
Scalar
,
GridView
::
dimension
,
GridView
::
dimension
>
stressTensor
(
0.0
);
for
(
unsigned
int
dimIdx
=
0
;
dimIdx
<
dim
;
++
dimIdx
)
{
for
(
unsigned
int
velIdx
=
0
;
velIdx
<
dim
;
++
velIdx
)
{
stressTensor
[
dimIdx
][
velIdx
]
=
0.5
*
velocityGradients_
[
elementI
D
][
dimIdx
][
velIdx
]
+
0.5
*
velocityGradients_
[
elementI
D
][
velIdx
][
dimIdx
];
stressTensor
[
dimIdx
][
velIdx
]
=
0.5
*
velocityGradients_
[
elementI
dx
][
dimIdx
][
velIdx
]
+
0.5
*
velocityGradients_
[
elementI
dx
][
velIdx
][
dimIdx
];
}
}
stressTensorScalarProduct_
[
elementI
D
]
=
0.0
;
stressTensorScalarProduct_
[
elementI
dx
]
=
0.0
;
for
(
unsigned
int
dimIdx
=
0
;
dimIdx
<
dim
;
++
dimIdx
)
{
for
(
unsigned
int
velIdx
=
0
;
velIdx
<
dim
;
++
velIdx
)
{
stressTensorScalarProduct_
[
elementI
D
]
+=
stressTensor
[
dimIdx
][
velIdx
]
*
stressTensor
[
dimIdx
][
velIdx
];
stressTensorScalarProduct_
[
elementI
dx
]
+=
stressTensor
[
dimIdx
][
velIdx
]
*
stressTensor
[
dimIdx
][
velIdx
];
}
}
...
...
@@ -383,16 +383,16 @@ public:
{
for
(
unsigned
int
velIdx
=
0
;
velIdx
<
dim
;
++
velIdx
)
{
vorticityTensor
[
dimIdx
][
velIdx
]
=
0.5
*
velocityGradients_
[
elementI
D
][
dimIdx
][
velIdx
]
-
0.5
*
velocityGradients_
[
elementI
D
][
velIdx
][
dimIdx
];
vorticityTensor
[
dimIdx
][
velIdx
]
=
0.5
*
velocityGradients_
[
elementI
dx
][
dimIdx
][
velIdx
]
-
0.5
*
velocityGradients_
[
elementI
dx
][
velIdx
][
dimIdx
];
}
}
vorticityTensorScalarProduct_
[
elementI
D
]
=
0.0
;
vorticityTensorScalarProduct_
[
elementI
dx
]
=
0.0
;
for
(
unsigned
int
dimIdx
=
0
;
dimIdx
<
dim
;
++
dimIdx
)
{
for
(
unsigned
int
velIdx
=
0
;
velIdx
<
dim
;
++
velIdx
)
{
vorticityTensorScalarProduct_
[
elementI
D
]
+=
vorticityTensor
[
dimIdx
][
velIdx
]
*
vorticityTensor
[
dimIdx
][
velIdx
];
vorticityTensorScalarProduct_
[
elementI
dx
]
+=
vorticityTensor
[
dimIdx
][
velIdx
]
*
vorticityTensor
[
dimIdx
][
velIdx
];
}
}
...
...
@@ -409,7 +409,7 @@ public:
VolumeVariables
volVars
;
volVars
.
update
(
elemSol
,
asImp_
(),
element
,
scv
);
kinematicViscosity_
[
elementI
D
]
=
volVars
.
viscosity
()
/
volVars
.
density
();
kinematicViscosity_
[
elementI
dx
]
=
volVars
.
viscosity
()
/
volVars
.
density
();
}
}
}
...
...
@@ -466,9 +466,9 @@ public:
public:
bool
calledUpdateStaticWallProperties
=
false
;
std
::
vector
<
unsigned
int
>
wallElementI
D
_
;
std
::
vector
<
unsigned
int
>
wallElementI
dx
_
;
std
::
vector
<
Scalar
>
wallDistance_
;
std
::
vector
<
std
::
array
<
std
::
array
<
unsigned
int
,
2
>
,
dim
>>
neighborI
D
_
;
std
::
vector
<
std
::
array
<
std
::
array
<
unsigned
int
,
2
>
,
dim
>>
neighborI
dx
_
;
std
::
vector
<
GlobalPosition
>
cellCenter_
;
std
::
vector
<
DimVector
>
velocity_
;
std
::
vector
<
DimVector
>
velocityMaximum_
;
...
...
dumux/freeflow/rans/twoeq/kepsilon/problem.hh
View file @
15129c16
...
...
@@ -98,7 +98,7 @@ public:
ParentType
::
updateStaticWallProperties
();
// update size and initial values of the global vectors
matchingPointI
D
_
.
resize
(
this
->
fvGridGeometry
().
elementMapper
().
size
(),
0
);
matchingPointI
dx
_
.
resize
(
this
->
fvGridGeometry
().
elementMapper
().
size
(),
0
);
storedDensity_
.
resize
(
this
->
fvGridGeometry
().
elementMapper
().
size
(),
0.0
);
storedDissipation_
.
resize
(
this
->
fvGridGeometry
().
elementMapper
().
size
(),
0.0
);
storedTurbulentKineticEnergy_
.
resize
(
this
->
fvGridGeometry
().
elementMapper
().
size
(),
0.0
);
...
...
@@ -118,7 +118,7 @@ public:
// update the stored eddy viscosities
for
(
const
auto
&
element
:
elements
(
this
->
fvGridGeometry
().
gridView
()))
{
unsigned
int
elementI
D
=
this
->
fvGridGeometry
().
elementMapper
().
index
(
element
);
unsigned
int
elementI
dx
=
this
->
fvGridGeometry
().
elementMapper
().
index
(
element
);
auto
fvGeometry
=
localView
(
this
->
fvGridGeometry
());
fvGeometry
.
bindElement
(
element
);
...
...
@@ -129,13 +129,13 @@ public:
PrimaryVariables
priVars
=
makePriVarsFromCellCenterPriVars
<
PrimaryVariables
>
(
cellCenterPriVars
);
auto
elemSol
=
elementSolution
<
typename
FVGridGeometry
::
LocalView
>
(
std
::
move
(
priVars
));
// NOTE: first update the turbulence quantities
storedDissipation_
[
elementI
D
]
=
elemSol
[
0
][
Indices
::
dissipationEqIdx
];
storedTurbulentKineticEnergy_
[
elementI
D
]
=
elemSol
[
0
][
Indices
::
turbulentKineticEnergyEqIdx
];
storedDissipation_
[
elementI
dx
]
=
elemSol
[
0
][
Indices
::
dissipationEqIdx
];
storedTurbulentKineticEnergy_
[
elementI
dx
]
=
elemSol
[
0
][
Indices
::
turbulentKineticEnergyEqIdx
];
// NOTE: then update the volVars
VolumeVariables
volVars
;
volVars
.
update
(
elemSol
,
asImp_
(),
element
,
scv
);
storedDynamicEddyViscosity_
[
elementI
D
]
=
volVars
.
calculateEddyViscosity
();
storedDensity_
[
elementI
D
]
=
volVars
.
density
();
storedDynamicEddyViscosity_
[
elementI
dx
]
=
volVars
.
calculateEddyViscosity
();
storedDensity_
[
elementI
dx
]
=
volVars
.
density
();
}
}
...
...
@@ -143,18 +143,18 @@ public:
unsigned
int
numElementsInNearWallRegion
=
0
;
for
(
const
auto
&
element
:
elements
(
this
->
fvGridGeometry
().
gridView
()))
{
unsigned
int
elementI
D
=
this
->
fvGridGeometry
().
elementMapper
().
index
(
element
);
unsigned
int
wallNormalAxis
=
asImp_
().
wallNormalAxis_
[
elementI
D
];
unsigned
int
neighborI
D
0
=
asImp_
().
neighborI
D
_
[
elementI
D
][
wallNormalAxis
][
0
];
unsigned
int
neighborI
D
1
=
asImp_
().
neighborI
D
_
[
elementI
D
][
wallNormalAxis
][
1
];
numElementsInNearWallRegion
=
inNearWallRegion
(
elementI
D
)
unsigned
int
elementI
dx
=
this
->
fvGridGeometry
().
elementMapper
().
index
(
element
);
unsigned
int
wallNormalAxis
=
asImp_
().
wallNormalAxis_
[
elementI
dx
];
unsigned
int
neighborI
dx
0
=
asImp_
().
neighborI
dx
_
[
elementI
dx
][
wallNormalAxis
][
0
];
unsigned
int
neighborI
dx
1
=
asImp_
().
neighborI
dx
_
[
elementI
dx
][
wallNormalAxis
][
1
];
numElementsInNearWallRegion
=
inNearWallRegion
(
elementI
dx
)
?
numElementsInNearWallRegion
+
1
:
numElementsInNearWallRegion
+
0
;
if
((
!
inNearWallRegion
(
elementI
D
)
&&
(
inNearWallRegion
(
neighborI
D
0
)
||
inNearWallRegion
(
neighborI
D
1
)))
||
(
!
inNearWallRegion
(
elementI
D
)
&&
elementI
D
==
asImp_
().
wallElementI
D
_
[
elementI
D
])
||
(
inNearWallRegion
(
elementI
D
)
&&
(
asImp_
().
wallElementI
D
_
[
neighborI
D
0
]
!=
asImp_
().
wallElementI
D
_
[
neighborI
D
1
])))
if
((
!
inNearWallRegion
(
elementI
dx
)
&&
(
inNearWallRegion
(
neighborI
dx
0
)
||
inNearWallRegion
(
neighborI
dx
1
)))
||
(
!
inNearWallRegion
(
elementI
dx
)
&&
elementI
dx
==
asImp_
().
wallElementI
dx
_
[
elementI
dx
])
||
(
inNearWallRegion
(
elementI
dx
)
&&
(
asImp_
().
wallElementI
dx
_
[
neighborI
dx
0
]
!=
asImp_
().
wallElementI
dx
_
[
neighborI
dx
1
])))
{
matchingPointI
D
_
[
asImp_
().
wallElementI
D
_
[
elementI
D
]]
=
elementI
D
;
matchingPointI
dx
_
[
asImp_
().
wallElementI
dx
_
[
elementI
dx
]]
=
elementI
dx
;
}
}
std
::
cout
<<
"numElementsInNearWallRegion: "
<<
numElementsInNearWallRegion
<<
std
::
endl
;
...
...
@@ -162,8 +162,8 @@ public:
// calculate the potential zeroeq eddy viscosities for two-layer model
for
(
const
auto
&
element
:
elements
(
this
->
fvGridGeometry
().
gridView
()))
{
unsigned
int
elementI
D
=
this
->
fvGridGeometry
().
elementMapper
().
index
(
element
);
zeroEqDynamicEddyViscosity_
[
elementI
D
]
=
zeroEqEddyViscosityModel
(
elementI
D
);
unsigned
int
elementI
dx
=
this
->
fvGridGeometry
().
elementMapper
().
index
(
element
);
zeroEqDynamicEddyViscosity_
[
elementI
dx
]
=
zeroEqEddyViscosityModel
(
elementI
dx
);
}
// then make them match at the matching point
...
...
@@ -173,24 +173,24 @@ public:
{
for
(
const
auto
&
element
:
elements
(
this
->
fvGridGeometry
().
gridView
()))
{
unsigned
int
elementI
D
=
this
->
fvGridGeometry
().
elementMapper
().
index
(
element
);
unsigned
int
matchingPointI
D
=
matchingPointI
D
_
[
asImp_
().
wallElementI
D
_
[
elementI
D
]];
unsigned
int
elementI
dx
=
this
->
fvGridGeometry
().
elementMapper
().
index
(
element
);
unsigned
int
matchingPointI
dx
=
matchingPointI
dx
_
[
asImp_
().
wallElementI
dx
_
[
elementI
dx
]];
Scalar
scalingFactor
=
storedDynamicEddyViscosity_
[
matchingPointI
D
]
/
zeroEqDynamicEddyViscosity_
[
matchingPointI
D
];
if
(
!
isMatchingPoint
(
elementI
D
)
Scalar
scalingFactor
=
storedDynamicEddyViscosity_
[
matchingPointI
dx
]
/
zeroEqDynamicEddyViscosity_
[
matchingPointI
dx
];
if
(
!
isMatchingPoint
(
elementI
dx
)
&&
!
std
::
isnan
(
scalingFactor
)
&&
!
std
::
isinf
(
scalingFactor
))
{
zeroEqDynamicEddyViscosity_
[
elementI
D
]
*=
scalingFactor
;
zeroEqDynamicEddyViscosity_
[
elementI
dx
]
*=
scalingFactor
;
}
}
for
(
const
auto
&
element
:
elements
(
this
->
fvGridGeometry
().
gridView
()))
{
unsigned
int
elementI
D
=
this
->
fvGridGeometry
().
elementMapper
().
index
(
element
);
unsigned
int
matchingPointI
D
=
matchingPointI
D
_
[
asImp_
().
wallElementI
D
_
[
elementI
D
]];
if
(
isMatchingPoint
(
elementI
D
))
unsigned
int
elementI
dx
=
this
->
fvGridGeometry
().
elementMapper
().
index
(
element
);
unsigned
int
matchingPointI
dx
=
matchingPointI
dx
_
[
asImp_
().
wallElementI
dx
_
[
elementI
dx
]];
if
(
isMatchingPoint
(
elementI
dx
))
{
zeroEqDynamicEddyViscosity_
[
matchingPointI
D
]
=
storedDynamicEddyViscosity_
[
matchingPointI
D
];
zeroEqDynamicEddyViscosity_
[
matchingPointI
dx
]
=
storedDynamicEddyViscosity_
[
matchingPointI
dx
];
}
}
}
...
...
@@ -199,109 +199,109 @@ public:
/*!
* \brief Returns if an element is located in the near-wall region
*/
const
bool
inNearWallRegion
(
unsigned
int
elementI
D
)
const
const
bool
inNearWallRegion
(
unsigned
int
elementI
dx
)
const
{
unsigned
int
wallElementI
D
=
asImp_
().
wallElementI
D
_
[
elementI
D
];
unsigned
int
matchingPointI
D
=
matchingPointI
D
_
[
wallElementI
D
];
return
wallElementI
D
==
matchingPointI
D
?
yPlusNominal
(
elementI
D
)
<
yPlusThreshold_
:
yPlus
(
elementI
D
)
<
yPlusThreshold_
;
unsigned
int
wallElementI
dx
=
asImp_
().
wallElementI
dx
_
[
elementI
dx
];
unsigned
int
matchingPointI
dx
=
matchingPointI
dx
_
[
wallElementI
dx
];
return
wallElementI
dx
==
matchingPointI
dx
?
yPlusNominal
(
elementI
dx
)
<
yPlusThreshold_
:
yPlus
(
elementI
dx
)
<
yPlusThreshold_
;
}
/*!
* \brief Returns if an element is the matching point
*/
const
bool
isMatchingPoint
(
unsigned
int
elementI
D
)
const
{
return
matchingPointI
D
_
[
asImp_
().
wallElementI
D
_
[
elementI
D
]]
==
elementI
D
;
}
const
bool
isMatchingPoint
(
unsigned
int
elementI
dx
)
const
{
return
matchingPointI
dx
_
[
asImp_
().
wallElementI
dx
_
[
elementI
dx
]]
==
elementI
dx
;
}
/*!
* \brief Returns the \f$ y^+ \f$ value at an element center
*/
const
Scalar
yPlus
(
unsigned
int
elementI
D
)
const
const
Scalar
yPlus
(
unsigned
int
elementI
dx
)
const
{
return
asImp_
().
wallDistance_
[
elementI
D
]
*
uStar
(
elementI
D
)
/
asImp_
().
kinematicViscosity_
[
elementI
D
];
return
asImp_
().
wallDistance_
[
elementI
dx
]
*
uStar
(
elementI
dx
)
/
asImp_
().
kinematicViscosity_
[
elementI
dx
];
}
/*!
* \brief Returns the nominal \f$ y^+ \f$ value at an element center
*/
const
Scalar
yPlusNominal
(
unsigned
int
elementI
D
)
const
const
Scalar
yPlusNominal
(
unsigned
int
elementI
dx
)
const
{
return
asImp_
().
wallDistance_
[
elementI
D
]
*
uStarNominal
(
elementI
D
)
/
asImp_
().
kinematicViscosity_
[
elementI
D
];
return
asImp_
().
wallDistance_
[
elementI
dx
]
*
uStarNominal
(
elementI
dx
)
/
asImp_
().
kinematicViscosity_
[
elementI
dx
];
}
/*!
* \brief Returns the kinematic eddy viscosity of a 0-Eq. model
*/
const
Scalar
zeroEqEddyViscosityModel
(
unsigned
int
elementI
D
)
const
const
Scalar
zeroEqEddyViscosityModel
(
unsigned
int
elementI
dx
)
const
{
using
std
::
abs
;
using
std
::
exp
;
using
std
::
sqrt
;
// use VanDriest's model
Scalar
yPlusValue
=
yPlus
(
elementI
D
);
Scalar
yPlusValue
=
yPlus
(
elementI
dx
);
Scalar
mixingLength
=
0.0
;
if
(
yPlusValue
>
0.0
)
{
mixingLength
=
asImp_
().
karmanConstant
()
*
asImp_
().
wallDistance_
[
elementI
D
]
mixingLength
=
asImp_
().
karmanConstant
()
*
asImp_
().
wallDistance_
[
elementI
dx
]
*
(
1.0
-
exp
(
-
yPlusValue
/
26.0
))
/
sqrt
(
1.0
-
exp
(
-
0.26
*
yPlusValue
));
}
unsigned
int
wallNormalAxis
=
asImp_
().
wallNormalAxis_
[
elementI
D
];
unsigned
int
flowNormalAxis
=
asImp_
().
flowNormalAxis_
[
elementI
D
];
Scalar
velocityGradient
=
asImp_
().
velocityGradients_
[
elementI
D
][
flowNormalAxis
][
wallNormalAxis
];
return
mixingLength
*
mixingLength
*
abs
(
velocityGradient
)
*
storedDensity_
[
elementI
D
];
unsigned
int
wallNormalAxis
=
asImp_
().
wallNormalAxis_
[
elementI
dx
];
unsigned
int
flowNormalAxis
=
asImp_
().
flowNormalAxis_
[
elementI
dx
];
Scalar
velocityGradient
=
asImp_
().
velocityGradients_
[
elementI
dx
][
flowNormalAxis
][
wallNormalAxis
];
return
mixingLength
*
mixingLength
*
abs
(
velocityGradient
)
*
storedDensity_
[
elementI
dx
];
}
//! \brief Returns the wall shear stress velocity
const
Scalar
uStar
(
unsigned
int
elementI
D
)
const
const
Scalar
uStar
(
unsigned
int
elementI
dx
)
const
{
using
std
::
abs
;
using
std
::
sqrt
;
unsigned
int
wallElementI
D
=
asImp_
().
wallElementI
D
_
[
elementI
D
];
unsigned
int
wallNormalAxis
=
asImp_
().
wallNormalAxis_
[
elementI
D
];
unsigned
int
flowNormalAxis
=
asImp_
().
flowNormalAxis_
[
elementI
D
];
return
sqrt
(
asImp_
().
kinematicViscosity_
[
wallElementI
D
]
*
abs
(
asImp_
().
velocityGradients_
[
wallElementI
D
][
flowNormalAxis
][
wallNormalAxis
]));
unsigned
int
wallElementI
dx
=
asImp_
().
wallElementI
dx
_
[
elementI
dx
];