Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
hyteg
hyteg
Commits
537bdc9b
Commit
537bdc9b
authored
Apr 07, 2022
by
Benjamin Mann
Browse files
Merge branch 'master' into mogli/adaptiveRefinement_optimization
parents
31e7f702
474976d2
Pipeline
#39025
passed with stages
in 129 minutes and 1 second
Changes
173
Pipelines
1
Expand all
Hide whitespace changes
Inline
Side-by-side
apps/Chorin.cpp
View file @
537bdc9b
...
...
@@ -19,7 +19,7 @@
*/
#include
<core/math/Constants.h>
#include
"hyteg/dgfunctionspace/DGUpwindOperator.hpp"
#include
"hyteg/dgfunctionspace
_old
/DGUpwindOperator.hpp"
#include
"hyteg/mesh/MeshInfo.hpp"
#include
"hyteg/p1functionspace/P1Function.hpp"
#include
"hyteg/p1functionspace/P1ConstantOperator.hpp"
...
...
@@ -121,11 +121,11 @@ int main( int argc, char* argv[] )
hyteg
::
P1Function
<
real_t
>
res
(
"res"
,
storage
,
minLevel
,
maxLevel
);
hyteg
::
P1Function
<
real_t
>
ones
(
"ones"
,
storage
,
minLevel
,
maxLevel
);
auto
u_dg
=
std
::
make_shared
<
hyteg
::
DGFunction
<
real_t
>
>
(
"u_dg"
,
storage
,
minLevel
,
maxLevel
);
auto
v_dg
=
std
::
make_shared
<
hyteg
::
DGFunction
<
real_t
>
>
(
"v_dg"
,
storage
,
minLevel
,
maxLevel
);
auto
u_dg
=
std
::
make_shared
<
hyteg
::
DGFunction
_old
<
real_t
>
>
(
"u_dg"
,
storage
,
minLevel
,
maxLevel
);
auto
v_dg
=
std
::
make_shared
<
hyteg
::
DGFunction
_old
<
real_t
>
>
(
"v_dg"
,
storage
,
minLevel
,
maxLevel
);
auto
u_dg_old
=
std
::
make_shared
<
hyteg
::
DGFunction
<
real_t
>
>
(
"u_dg"
,
storage
,
minLevel
,
maxLevel
);
auto
v_dg_old
=
std
::
make_shared
<
hyteg
::
DGFunction
<
real_t
>
>
(
"v_dg"
,
storage
,
minLevel
,
maxLevel
);
auto
u_dg_old
=
std
::
make_shared
<
hyteg
::
DGFunction
_old
<
real_t
>
>
(
"u_dg"
,
storage
,
minLevel
,
maxLevel
);
auto
v_dg_old
=
std
::
make_shared
<
hyteg
::
DGFunction
_old
<
real_t
>
>
(
"v_dg"
,
storage
,
minLevel
,
maxLevel
);
hyteg
::
P1ConstantLaplaceOperator
A
(
storage
,
minLevel
,
maxLevel
);
hyteg
::
P1ConstantLaplaceOperator
Ascaled
(
storage
,
minLevel
,
maxLevel
);
...
...
apps/Surrogates/VariableCoeffLaplace_Surrogates.cpp
View file @
537bdc9b
...
...
@@ -384,11 +384,11 @@ void showStencilFunction(std::shared_ptr<PrimitiveStorage> storage, const uint_t
P2Form_divKgrad
form
;
form
.
setGeometryMap
(
face
.
getGeometryMap
());
Point3D
x0
(
face
.
coords
[
0
]),
x
;
Point3D
x0
(
face
.
getCoordinates
()
[
0
]),
x
;
real_t
h
=
1.0
/
(
walberla
::
real_c
(
levelinfo
::
num_microvertices_per_edge
(
level
)
-
1
));
Point3D
d0
=
h
*
(
face
.
coords
[
1
]
-
face
.
coords
[
0
]);
Point3D
d2
=
h
*
(
face
.
coords
[
2
]
-
face
.
coords
[
0
]);
Point3D
d0
=
h
*
(
face
.
getCoordinates
()[
1
]
-
face
.
getCoordinates
()
[
0
]);
Point3D
d2
=
h
*
(
face
.
getCoordinates
()[
2
]
-
face
.
getCoordinates
()
[
0
]);
// directions
const
Point3D
dirS
=
-
d2
;
...
...
@@ -412,8 +412,8 @@ void showStencilFunction(std::shared_ptr<PrimitiveStorage> storage, const uint_t
P2
::
variablestencil
::
macroface
::
assembleStencil
(
form
,
x
,
dirS
,
dirSE
,
dirE
,
dirN
,
dirNW
,
dirW
,
dirNE
,
VtVStencil
,
EtVStencil
,
VtEStencil
,
EtEStencil
);
uint
_t
i
=
it
.
col
();
uint
_t
j
=
it
.
row
();
idx
_t
i
=
it
.
col
();
idx
_t
j
=
it
.
row
();
// VERTEX DoF
if
(
!
vertexdof
::
macroface
::
isVertexOnBoundary
(
level
,
it
))
...
...
apps/benchmarks/ApplyPerformanceAnalysis-2D-P2/AbstractApply.cpp
View file @
537bdc9b
...
...
@@ -37,9 +37,9 @@ void apply_2d_vertex_to_vertex( double* WALBERLA_RESTRICT dstPtr,
const
uint_t
level
)
{
uint_t
rowsize
=
levelinfo
::
num_microvertices_per_edge
(
level
);
for
(
uint
_t
j
=
1
;
j
<
rowsize
-
2
;
++
j
)
for
(
idx
_t
j
=
1
;
j
<
idx_t
(
rowsize
)
-
2
;
++
j
)
{
for
(
uint
_t
i
=
1
;
i
<
rowsize
-
j
-
1
;
++
i
)
for
(
idx
_t
i
=
1
;
i
<
idx_t
(
rowsize
)
-
j
-
1
;
++
i
)
{
auto
tmp
=
real_t
(
0
);
for
(
const
auto
direction
:
vertexdof
::
macroface
::
neighborsWithCenter
)
...
...
apps/benchmarks/CompareInterpolate.cpp
View file @
537bdc9b
...
...
@@ -46,17 +46,17 @@ inline void interpolateStdFunction( Face&
uint_t
rowsize
=
levelinfo
::
num_microvertices_per_edge
(
Level
);
Point3D
x
,
x0
;
auto
dstPtr
=
faceMemory
->
getPointer
(
Level
);
x0
=
face
.
coords
[
0
];
Point3D
d0
=
(
face
.
coords
[
1
]
-
face
.
coords
[
0
]
)
/
(
walberla
::
real_c
(
rowsize
-
1
)
);
Point3D
d2
=
(
face
.
coords
[
2
]
-
face
.
coords
[
0
]
)
/
(
walberla
::
real_c
(
rowsize
-
1
)
);
x0
=
face
.
getCoordinates
()
[
0
];
Point3D
d0
=
(
face
.
getCoordinates
()[
1
]
-
face
.
getCoordinates
()
[
0
]
)
/
(
walberla
::
real_c
(
rowsize
-
1
)
);
Point3D
d2
=
(
face
.
getCoordinates
()[
2
]
-
face
.
getCoordinates
()
[
0
]
)
/
(
walberla
::
real_c
(
rowsize
-
1
)
);
uint_t
inner_rowsize
=
rowsize
;
for
(
uint
_t
i
=
1
;
i
<
rowsize
-
2
;
++
i
)
for
(
idx
_t
i
=
1
;
i
<
idx_t
(
rowsize
)
-
2
;
++
i
)
{
x
=
x0
;
x
+=
real_c
(
i
)
*
d2
+
d0
;
for
(
uint
_t
j
=
1
;
j
<
inner_rowsize
-
2
;
++
j
)
for
(
idx
_t
j
=
1
;
j
<
idx_t
(
inner_rowsize
)
-
2
;
++
j
)
{
dstPtr
[
vertexdof
::
macroface
::
indexFromVertex
(
Level
,
j
,
i
,
stencilDirection
::
VERTEX_C
)]
=
expr
(
x
);
x
+=
d0
;
...
...
@@ -74,17 +74,17 @@ inline void
uint_t
rowsize
=
levelinfo
::
num_microvertices_per_edge
(
Level
);
Point3D
x
,
x0
;
auto
dstPtr
=
faceMemory
->
getPointer
(
Level
);
x0
=
face
.
coords
[
0
];
Point3D
d0
=
(
face
.
coords
[
1
]
-
face
.
coords
[
0
]
)
/
(
walberla
::
real_c
(
rowsize
-
1
)
);
Point3D
d2
=
(
face
.
coords
[
2
]
-
face
.
coords
[
0
]
)
/
(
walberla
::
real_c
(
rowsize
-
1
)
);
x0
=
face
.
getCoordinates
()
[
0
];
Point3D
d0
=
(
face
.
getCoordinates
()[
1
]
-
face
.
getCoordinates
()
[
0
]
)
/
(
walberla
::
real_c
(
rowsize
-
1
)
);
Point3D
d2
=
(
face
.
getCoordinates
()[
2
]
-
face
.
getCoordinates
()
[
0
]
)
/
(
walberla
::
real_c
(
rowsize
-
1
)
);
uint_t
inner_rowsize
=
rowsize
;
for
(
uint
_t
i
=
1
;
i
<
rowsize
-
2
;
++
i
)
for
(
idx
_t
i
=
1
;
i
<
idx_t
(
rowsize
)
-
2
;
++
i
)
{
x
=
x0
;
x
+=
real_c
(
i
)
*
d2
+
d0
;
for
(
uint
_t
j
=
1
;
j
<
inner_rowsize
-
2
;
++
j
)
for
(
idx
_t
j
=
1
;
j
<
idx_t
(
inner_rowsize
)
-
2
;
++
j
)
{
dstPtr
[
vertexdof
::
macroface
::
indexFromVertex
(
Level
,
j
,
i
,
stencilDirection
::
VERTEX_C
)]
=
expr
(
x
);
x
+=
d0
;
...
...
@@ -103,17 +103,17 @@ inline void interpolateFunctor( Face&
uint_t
rowsize
=
levelinfo
::
num_microvertices_per_edge
(
Level
);
Point3D
x
,
x0
;
auto
dstPtr
=
faceMemory
->
getPointer
(
Level
);
x0
=
face
.
coords
[
0
];
Point3D
d0
=
(
face
.
coords
[
1
]
-
face
.
coords
[
0
]
)
/
(
walberla
::
real_c
(
rowsize
-
1
)
);
Point3D
d2
=
(
face
.
coords
[
2
]
-
face
.
coords
[
0
]
)
/
(
walberla
::
real_c
(
rowsize
-
1
)
);
x0
=
face
.
getCoordinates
()
[
0
];
Point3D
d0
=
(
face
.
getCoordinates
()[
1
]
-
face
.
getCoordinates
()
[
0
]
)
/
(
walberla
::
real_c
(
rowsize
-
1
)
);
Point3D
d2
=
(
face
.
getCoordinates
()[
2
]
-
face
.
getCoordinates
()
[
0
]
)
/
(
walberla
::
real_c
(
rowsize
-
1
)
);
uint_t
inner_rowsize
=
rowsize
;
for
(
uint
_t
i
=
1
;
i
<
rowsize
-
2
;
++
i
)
for
(
idx
_t
i
=
1
;
i
<
idx_t
(
rowsize
)
-
2
;
++
i
)
{
x
=
x0
;
x
+=
real_c
(
i
)
*
d2
+
d0
;
for
(
uint
_t
j
=
1
;
j
<
inner_rowsize
-
2
;
++
j
)
for
(
idx
_t
j
=
1
;
j
<
idx_t
(
inner_rowsize
)
-
2
;
++
j
)
{
dstPtr
[
vertexdof
::
macroface
::
indexFromVertex
(
Level
,
j
,
i
,
stencilDirection
::
VERTEX_C
)]
=
exprFunctor
(
x
);
x
+=
d0
;
...
...
@@ -130,18 +130,18 @@ inline void interpolateWithoutFunction( Face& face, const PrimitiveDataID< Funct
uint_t
rowsize
=
levelinfo
::
num_microvertices_per_edge
(
Level
);
Point3D
x
,
x0
;
auto
dstPtr
=
faceMemory
->
getPointer
(
Level
);
x0
=
face
.
coords
[
0
];
Point3D
d0
=
(
face
.
coords
[
1
]
-
face
.
coords
[
0
]
)
/
(
walberla
::
real_c
(
rowsize
-
1
)
);
Point3D
d2
=
(
face
.
coords
[
2
]
-
face
.
coords
[
0
]
)
/
(
walberla
::
real_c
(
rowsize
-
1
)
);
x0
=
face
.
getCoordinates
()
[
0
];
Point3D
d0
=
(
face
.
getCoordinates
()[
1
]
-
face
.
getCoordinates
()
[
0
]
)
/
(
walberla
::
real_c
(
rowsize
-
1
)
);
Point3D
d2
=
(
face
.
getCoordinates
()[
2
]
-
face
.
getCoordinates
()
[
0
]
)
/
(
walberla
::
real_c
(
rowsize
-
1
)
);
uint_t
inner_rowsize
=
rowsize
;
for
(
uint
_t
i
=
1
;
i
<
rowsize
-
2
;
++
i
)
for
(
idx
_t
i
=
1
;
i
<
idx_t
(
rowsize
)
-
2
;
++
i
)
{
x
=
x0
;
x
+=
real_c
(
i
)
*
d2
+
d0
;
for
(
uint
_t
j
=
1
;
j
<
inner_rowsize
-
2
;
++
j
)
for
(
idx
_t
j
=
1
;
j
<
idx_t
(
inner_rowsize
)
-
2
;
++
j
)
{
dstPtr
[
vertexdof
::
macroface
::
indexFromVertex
(
Level
,
j
,
i
,
stencilDirection
::
VERTEX_C
)]
=
sqrt
(
x
[
0
]
*
x
[
0
]
+
x
[
1
]
*
x
[
1
]
);
...
...
apps/dg0_transport.cpp
View file @
537bdc9b
...
...
@@ -21,8 +21,8 @@
#include
<core/timing/Timer.h>
#include
"hyteg/dataexport/VTKOutput.hpp"
#include
"hyteg/dgfunctionspace/DGUpwindOperator.hpp"
#include
"hyteg/dgfunctionspace/DGFunction.hpp"
#include
"hyteg/dgfunctionspace
_old
/DGUpwindOperator.hpp"
#include
"hyteg/dgfunctionspace
_old
/DGFunction.hpp"
#include
"hyteg/mesh/MeshInfo.hpp"
#include
"hyteg/p1functionspace/P1Function.hpp"
#include
"hyteg/primitivestorage/PrimitiveStorage.hpp"
...
...
@@ -74,8 +74,8 @@ int main(int argc, char* argv[])
std
::
shared_ptr
<
hyteg
::
PrimitiveStorage
>
storage
=
std
::
make_shared
<
hyteg
::
PrimitiveStorage
>
(
setupStorage
);
std
::
shared_ptr
<
hyteg
::
DGFunction
<
real_t
>>
c_old
=
std
::
make_shared
<
hyteg
::
DGFunction
<
real_t
>>
(
"c_old"
,
storage
,
minLevel
,
maxLevel
);
std
::
shared_ptr
<
hyteg
::
DGFunction
<
real_t
>>
c
=
std
::
make_shared
<
hyteg
::
DGFunction
<
real_t
>>
(
"c"
,
storage
,
minLevel
,
maxLevel
);
std
::
shared_ptr
<
hyteg
::
DGFunction
_old
<
real_t
>>
c_old
=
std
::
make_shared
<
hyteg
::
DGFunction
_old
<
real_t
>>
(
"c_old"
,
storage
,
minLevel
,
maxLevel
);
std
::
shared_ptr
<
hyteg
::
DGFunction
_old
<
real_t
>>
c
=
std
::
make_shared
<
hyteg
::
DGFunction
_old
<
real_t
>>
(
"c"
,
storage
,
minLevel
,
maxLevel
);
std
::
shared_ptr
<
hyteg
::
P1Function
<
real_t
>>
u
=
std
::
make_shared
<
hyteg
::
P1Function
<
real_t
>>
(
"u"
,
storage
,
minLevel
,
maxLevel
);
std
::
shared_ptr
<
hyteg
::
P1Function
<
real_t
>>
v
=
std
::
make_shared
<
hyteg
::
P1Function
<
real_t
>>
(
"v"
,
storage
,
minLevel
,
maxLevel
);
...
...
apps/geophysics_plume.cpp
View file @
537bdc9b
...
...
@@ -24,8 +24,8 @@
#include
"hyteg/composites/P1StokesFunction.hpp"
#include
"hyteg/composites/P1P1StokesOperator.hpp"
#include
"hyteg/dataexport/VTKOutput.hpp"
#include
"hyteg/dgfunctionspace/DGFunction.hpp"
#include
"hyteg/dgfunctionspace/DGUpwindOperator.hpp"
#include
"hyteg/dgfunctionspace
_old
/DGFunction.hpp"
#include
"hyteg/dgfunctionspace
_old
/DGUpwindOperator.hpp"
#include
"hyteg/functions/FunctionProperties.hpp"
#include
"hyteg/gridtransferoperators/P1P1StokesToP1P1StokesProlongation.hpp"
#include
"hyteg/gridtransferoperators/P1P1StokesToP1P1StokesRestriction.hpp"
...
...
@@ -101,10 +101,10 @@ int main( int argc, char* argv[] )
#endif
// Setting up Functions
auto
c_old
=
std
::
make_shared
<
hyteg
::
DGFunction
<
real_t
>
>
(
"c"
,
storage
,
minLevel
,
maxLevel
);
auto
c
=
std
::
make_shared
<
hyteg
::
DGFunction
<
real_t
>
>
(
"c"
,
storage
,
minLevel
,
maxLevel
);
auto
c_old
=
std
::
make_shared
<
hyteg
::
DGFunction
_old
<
real_t
>
>
(
"c"
,
storage
,
minLevel
,
maxLevel
);
auto
c
=
std
::
make_shared
<
hyteg
::
DGFunction
_old
<
real_t
>
>
(
"c"
,
storage
,
minLevel
,
maxLevel
);
auto
f_dg
=
std
::
make_shared
<
hyteg
::
DGFunction
<
real_t
>
>
(
"f_dg"
,
storage
,
minLevel
,
maxLevel
);
auto
f_dg
=
std
::
make_shared
<
hyteg
::
DGFunction
_old
<
real_t
>
>
(
"f_dg"
,
storage
,
minLevel
,
maxLevel
);
auto
r
=
std
::
make_shared
<
hyteg
::
P1StokesFunction
<
real_t
>
>
(
"r"
,
storage
,
minLevel
,
maxLevel
);
auto
f
=
std
::
make_shared
<
hyteg
::
P1StokesFunction
<
real_t
>
>
(
"f"
,
storage
,
minLevel
,
maxLevel
);
...
...
@@ -165,11 +165,11 @@ int main( int argc, char* argv[] )
for
(
uint_t
lvl
=
minLevel
;
lvl
<=
maxLevel
;
++
lvl
)
{
uint_t
tmpDofStokes
=
numberOfGlobalDoFs
<
hyteg
::
P1StokesFunctionTag
>
(
*
storage
,
lvl
);
uint_t
tmpDofTemp
=
numberOfGlobalDoFs
<
hyteg
::
DGFunctionTag
>
(
*
storage
,
lvl
);
//
uint_t tmpDofTemp = numberOfGlobalDoFs< hyteg::DGFunctionTag >( *storage, lvl );
WALBERLA_LOG_DETAIL_ON_ROOT
(
"Stokes DoFs on level "
<<
lvl
<<
" : "
<<
tmpDofStokes
);
WALBERLA_LOG_DETAIL_ON_ROOT
(
"Temperature DoFs on level "
<<
lvl
<<
" : "
<<
tmpDofTemp
);
//
WALBERLA_LOG_DETAIL_ON_ROOT( "Temperature DoFs on level " << lvl << " : " << tmpDofTemp );
totalGlobalDofsStokes
+=
tmpDofStokes
;
totalGlobalDofsTemp
+=
tmpDofTemp
;
//
totalGlobalDofsTemp += tmpDofTemp;
}
WALBERLA_LOG_INFO_ON_ROOT
(
"Total Stokes DoFs on all level :"
<<
totalGlobalDofsStokes
);
WALBERLA_LOG_INFO_ON_ROOT
(
"Total Temperature DoFs on all level :"
<<
totalGlobalDofsTemp
);
...
...
src/coupling_hyteg_convection_particles/MMOCTransport.hpp
View file @
537bdc9b
...
...
@@ -58,39 +58,40 @@ enum class TimeSteppingScheme
RK4
,
// fourth order
};
const
static
std
::
vector
<
std
::
vector
<
real_t
>
>
RK_A_ExplicitEuler
=
{
{{
0
}}
};
const
static
std
::
vector
<
real_t
>
RK_b_ExplicitEuler
=
{
1
};
const
static
std
::
vector
<
real_t
>
RK_c_ExplicitEuler
=
{
0
};
const
static
std
::
vector
<
std
::
vector
<
real_t
>
>
RK_A_ExplicitEuler
=
{
{
{
0
}
}
};
const
static
std
::
vector
<
real_t
>
RK_b_ExplicitEuler
=
{
1
};
const
static
std
::
vector
<
real_t
>
RK_c_ExplicitEuler
=
{
0
};
const
static
std
::
vector
<
std
::
vector
<
real_t
>
>
RK_A_RK3
=
{
{{
0
,
0
,
0
},
{
0.5
,
0
,
0
},
{
-
1.
,
2.
,
0
}}
};
const
static
std
::
vector
<
real_t
>
RK_b_RK3
=
{
1.
/
6.
,
2.
/
3.
,
1.
/
6.
};
const
static
std
::
vector
<
real_t
>
RK_c_RK3
=
{
0
,
0.5
,
1.
};
const
static
std
::
vector
<
std
::
vector
<
real_t
>
>
RK_A_RK3
=
{
{
{
0
,
0
,
0
},
{
0.5
,
0
,
0
},
{
-
1.
,
2.
,
0
}
}
};
const
static
std
::
vector
<
real_t
>
RK_b_RK3
=
{
1.
/
6.
,
2.
/
3.
,
1.
/
6.
};
const
static
std
::
vector
<
real_t
>
RK_c_RK3
=
{
0
,
0.5
,
1.
};
const
static
std
::
vector
<
std
::
vector
<
real_t
>
>
RK_A_Ralston
=
{
{{
0
,
0
,
0
},
{
0.5
,
0
,
0
},
{
0
,
0.75
,
0
}}
};
const
static
std
::
vector
<
real_t
>
RK_b_Ralston
=
{
2.
/
9.
,
1.
/
3.
,
4.
/
9.
};
const
static
std
::
vector
<
real_t
>
RK_c_Ralston
=
{
0
,
0.5
,
0.75
};
const
static
std
::
vector
<
std
::
vector
<
real_t
>
>
RK_A_Ralston
=
{
{
{
0
,
0
,
0
},
{
0.5
,
0
,
0
},
{
0
,
0.75
,
0
}
}
};
const
static
std
::
vector
<
real_t
>
RK_b_Ralston
=
{
2.
/
9.
,
1.
/
3.
,
4.
/
9.
};
const
static
std
::
vector
<
real_t
>
RK_c_Ralston
=
{
0
,
0.5
,
0.75
};
const
static
std
::
vector
<
std
::
vector
<
real_t
>
>
RK_A_RK4
=
{{{
0
,
0
,
0
,
0
},
{
0.5
,
0
,
0
,
0
},
{
0
,
0.5
,
0
,
0
},
{
0
,
0
,
1
,
0
}}};
const
static
std
::
vector
<
real_t
>
RK_b_RK4
=
{
1.
/
6.
,
1.
/
3.
,
1.
/
3.
,
1.
/
6.
};
const
static
std
::
vector
<
real_t
>
RK_c_RK4
=
{
0
,
0.5
,
0.5
,
1.
};
const
static
std
::
vector
<
std
::
vector
<
real_t
>
>
RK_A_RK4
=
{
{
{
0
,
0
,
0
,
0
},
{
0.5
,
0
,
0
,
0
},
{
0
,
0.5
,
0
,
0
},
{
0
,
0
,
1
,
0
}
}
};
const
static
std
::
vector
<
real_t
>
RK_b_RK4
=
{
1.
/
6.
,
1.
/
3.
,
1.
/
3.
,
1.
/
6.
};
const
static
std
::
vector
<
real_t
>
RK_c_RK4
=
{
0
,
0.5
,
0.5
,
1.
};
const
static
std
::
map
<
TimeSteppingScheme
,
std
::
vector
<
std
::
vector
<
real_t
>
>
>
RK_A
=
{
{
TimeSteppingScheme
::
ExplicitEuler
,
RK_A_ExplicitEuler
},
{
TimeSteppingScheme
::
RK3
,
RK_A_RK3
},
{
TimeSteppingScheme
::
Ralston
,
RK_A_Ralston
},
{
TimeSteppingScheme
::
RK4
,
RK_A_RK4
}
};
{
TimeSteppingScheme
::
ExplicitEuler
,
RK_A_ExplicitEuler
},
{
TimeSteppingScheme
::
RK3
,
RK_A_RK3
},
{
TimeSteppingScheme
::
Ralston
,
RK_A_Ralston
},
{
TimeSteppingScheme
::
RK4
,
RK_A_RK4
}
};
const
static
std
::
map
<
TimeSteppingScheme
,
std
::
vector
<
real_t
>
>
RK_b
=
{
{
TimeSteppingScheme
::
ExplicitEuler
,
RK_b_ExplicitEuler
},
{
TimeSteppingScheme
::
RK3
,
RK_b_RK3
},
{
TimeSteppingScheme
::
Ralston
,
RK_b_Ralston
},
{
TimeSteppingScheme
::
RK4
,
RK_b_RK4
}
};
{
TimeSteppingScheme
::
ExplicitEuler
,
RK_b_ExplicitEuler
},
{
TimeSteppingScheme
::
RK3
,
RK_b_RK3
},
{
TimeSteppingScheme
::
Ralston
,
RK_b_Ralston
},
{
TimeSteppingScheme
::
RK4
,
RK_b_RK4
}
};
const
static
std
::
map
<
TimeSteppingScheme
,
std
::
vector
<
real_t
>
>
RK_c
=
{
{
TimeSteppingScheme
::
ExplicitEuler
,
RK_c_ExplicitEuler
},
{
TimeSteppingScheme
::
RK3
,
RK_c_RK3
},
{
TimeSteppingScheme
::
Ralston
,
RK_c_Ralston
},
{
TimeSteppingScheme
::
RK4
,
RK_c_RK4
}
};
{
TimeSteppingScheme
::
ExplicitEuler
,
RK_c_ExplicitEuler
},
{
TimeSteppingScheme
::
RK3
,
RK_c_RK3
},
{
TimeSteppingScheme
::
Ralston
,
RK_c_Ralston
},
{
TimeSteppingScheme
::
RK4
,
RK_c_RK4
}
};
inline
std
::
vector
<
PrimitiveID
>
getNeighboringPrimitives
(
const
PrimitiveID
&
primitiveID
,
const
PrimitiveStorage
&
storage
)
{
...
...
@@ -147,12 +148,12 @@ inline void updateParticlePosition( const PrimitiveStorage&
Point3D
computationalLocation
;
face
->
getGeometryMap
()
->
evalFinv
(
toPoint3D
(
p
->
getPosition
()
),
computationalLocation
);
Point2D
computationalLocation2D
(
{
computationalLocation
[
0
],
computationalLocation
[
1
]}
);
Point2D
computationalLocation2D
(
{
computationalLocation
[
0
],
computationalLocation
[
1
]
}
);
if
(
isPointInTriangle
(
computationalLocation2D
,
Point2D
(
{
face
->
getCoordinates
().
at
(
0
)[
0
],
face
->
getCoordinates
().
at
(
0
)[
1
]}
),
Point2D
(
{
face
->
getCoordinates
().
at
(
1
)[
0
],
face
->
getCoordinates
().
at
(
1
)[
1
]}
),
Point2D
(
{
face
->
getCoordinates
().
at
(
2
)[
0
],
face
->
getCoordinates
().
at
(
2
)[
1
]}
)
)
)
Point2D
(
{
face
->
getCoordinates
().
at
(
0
)[
0
],
face
->
getCoordinates
().
at
(
0
)[
1
]
}
),
Point2D
(
{
face
->
getCoordinates
().
at
(
1
)[
0
],
face
->
getCoordinates
().
at
(
1
)[
1
]
}
),
Point2D
(
{
face
->
getCoordinates
().
at
(
2
)[
0
],
face
->
getCoordinates
().
at
(
2
)[
1
]
}
)
)
)
{
p
->
setContainingPrimitive
(
faceID
);
foundByPointLocation
=
true
;
...
...
@@ -168,13 +169,13 @@ inline void updateParticlePosition( const PrimitiveStorage&
const
auto
neighborFace
=
storage
.
getFace
(
neighborFaceID
);
Point3D
computationalLocationNeighbor
;
neighborFace
->
getGeometryMap
()
->
evalFinv
(
toPoint3D
(
p
->
getPosition
()
),
computationalLocationNeighbor
);
Point2D
computationalLocationNeighbor2D
(
{
computationalLocationNeighbor
[
0
],
computationalLocationNeighbor
[
1
]}
);
Point2D
computationalLocationNeighbor2D
(
{
computationalLocationNeighbor
[
0
],
computationalLocationNeighbor
[
1
]
}
);
if
(
isPointInTriangle
(
computationalLocationNeighbor2D
,
Point2D
(
{
neighborFace
->
getCoordinates
().
at
(
0
)[
0
],
neighborFace
->
getCoordinates
().
at
(
0
)[
1
]}
),
Point2D
(
{
neighborFace
->
getCoordinates
().
at
(
1
)[
0
],
neighborFace
->
getCoordinates
().
at
(
1
)[
1
]}
),
Point2D
(
{
neighborFace
->
getCoordinates
().
at
(
2
)[
0
],
neighborFace
->
getCoordinates
().
at
(
2
)[
1
]}
)
)
)
Point2D
(
{
neighborFace
->
getCoordinates
().
at
(
0
)[
0
],
neighborFace
->
getCoordinates
().
at
(
0
)[
1
]
}
),
Point2D
(
{
neighborFace
->
getCoordinates
().
at
(
1
)[
0
],
neighborFace
->
getCoordinates
().
at
(
1
)[
1
]
}
),
Point2D
(
{
neighborFace
->
getCoordinates
().
at
(
2
)[
0
],
neighborFace
->
getCoordinates
().
at
(
2
)[
1
]
}
)
)
)
{
// set it to the first neighbor we found to contain the particle
p
->
setContainingPrimitive
(
neighborFaceID
);
...
...
@@ -262,8 +263,7 @@ inline void updateParticlePosition( const PrimitiveStorage&
else
{
// check for neighbor cells if we did not find it in its previous cell
const
auto
&
neighboringCells
=
cell
->
getIndirectNeighborCellIDs
();
for
(
const
auto
&
neighborCellID
:
neighboringCells
)
for
(
const
auto
&
neighborCellID
:
cell
->
getIndirectNeighborCellIDsOverVertices
()
)
{
WALBERLA_ASSERT
(
storage
.
cellExistsLocally
(
neighborCellID
)
||
storage
.
cellExistsInNeighborhood
(
neighborCellID
)
);
...
...
@@ -308,8 +308,7 @@ inline void updateParticlePosition( const PrimitiveStorage&
}
else
{
const
auto
&
neighboringCells
=
cell
->
getIndirectNeighborCellIDs
();
for
(
const
auto
&
neighborCellID
:
neighboringCells
)
for
(
const
auto
&
neighborCellID
:
cell
->
getIndirectNeighborCellIDsOverVertices
()
)
{
WALBERLA_ASSERT
(
storage
.
cellExistsLocally
(
neighborCellID
)
||
storage
.
cellExistsInNeighborhood
(
neighborCellID
)
);
...
...
@@ -524,10 +523,10 @@ inline uint_t initializeParticles( walberla::convection_particles::data::Particl
particleStorage
.
clear
();
const
uint_t
rank
=
uint_c
(
walberla
::
mpi
::
MPIManager
::
instance
()
->
rank
()
);
const
uint_t
rank
=
uint_c
(
walberla
::
mpi
::
MPIManager
::
instance
()
->
rank
()
);
// const std::vector< std::vector< real_t > >& A = RK_A.at( timeSteppingScheme ); //this seems to be unused?
const
std
::
vector
<
real_t
>&
b
=
RK_b
.
at
(
timeSteppingScheme
);
const
uint_t
rkStages
=
b
.
size
();
const
std
::
vector
<
real_t
>&
b
=
RK_b
.
at
(
timeSteppingScheme
);
const
uint_t
rkStages
=
b
.
size
();
if
constexpr
(
std
::
is_same
<
FunctionType
,
P1Function
<
real_t
>
>::
value
)
{
...
...
@@ -806,10 +805,10 @@ inline void particleIntegration( walberla::convection_particles::data::ParticleS
storage
.
getTimingTree
()
->
start
(
"Evaluate at particle position"
);
std
::
vector
<
real_t
>
results
(
{
0
,
0
}
);
std
::
vector
<
real_t
>
resultsLastTimeStep
(
{
0
,
0
}
);
std
::
vector
<
FunctionType
>
functions
=
{
ux
,
uy
};
std
::
vector
<
FunctionType
>
functionsLastTimeStep
=
{
uxLastTimeStep
,
uyLastTimeStep
};
std
::
vector
<
real_t
>
results
(
{
0
,
0
}
);
std
::
vector
<
real_t
>
resultsLastTimeStep
(
{
0
,
0
}
);
std
::
vector
<
FunctionType
>
functions
=
{
ux
,
uy
};
std
::
vector
<
FunctionType
>
functionsLastTimeStep
=
{
uxLastTimeStep
,
uyLastTimeStep
};
if
(
storage
.
hasGlobalCells
()
)
{
results
.
push_back
(
0
);
...
...
@@ -944,15 +943,22 @@ inline void evaluateTemperature( walberla::convection_particles::data::ParticleS
const
int
TAG
=
98234
;
#ifdef _MSC_VER
//need a first receive to avoid blocking on windows while communicating with itself
int
selfCommMessage
=
0
;
int
selfCommMessage
=
0
;
MPI_Request
selfCommRequest
;
MPI_Irecv
(
&
selfCommMessage
,
1
,
MPI_INT
,
walberla
::
mpi
::
MPIManager
::
instance
()
->
rank
(),
TAG
,
walberla
::
mpi
::
MPIManager
::
instance
()
->
comm
(),
&
selfCommRequest
);
MPI_Irecv
(
&
selfCommMessage
,
1
,
MPI_INT
,
walberla
::
mpi
::
MPIManager
::
instance
()
->
rank
(),
TAG
,
walberla
::
mpi
::
MPIManager
::
instance
()
->
comm
(),
&
selfCommRequest
);
#endif
for
(
const
auto
&
p
:
particleStorage
)
{
if
(
numParticlesToSendToRank
.
count
(
p
.
getStartProcess
()
)
==
0
){
if
(
numParticlesToSendToRank
.
count
(
p
.
getStartProcess
()
)
==
0
)
{
numParticlesToSendToRank
[
p
.
getStartProcess
()]
=
0
;
sendRequests
[
p
.
getStartProcess
()]
=
MPI_Request
();
sendRequests
[
p
.
getStartProcess
()]
=
MPI_Request
();
}
numParticlesToSendToRank
[
p
.
getStartProcess
()]
++
;
}
...
...
@@ -977,21 +983,23 @@ inline void evaluateTemperature( walberla::convection_particles::data::ParticleS
}
int
numReceivedParticleLocations
=
0
;
#ifdef _MSC_VER
#ifdef _MSC_VER
//get self communication
{
int
ready
;
int
ready
;
MPI_Status
status
;
MPI_Test
(
&
selfCommRequest
,
&
ready
,
&
status
);
if
(
ready
){
numReceivedParticleLocations
=
selfCommMessage
;
MPI_Test
(
&
selfCommRequest
,
&
ready
,
&
status
);
if
(
ready
)
{
numReceivedParticleLocations
=
selfCommMessage
;
numParticlesToReceiveFromRank
[
uint_c
(
status
.
MPI_SOURCE
)]
=
numReceivedParticleLocations
;
}
else
{
WALBERLA_LOG_INFO
(
"somethings very wrong here"
);
}
else
{
WALBERLA_LOG_INFO
(
"somethings very wrong here"
);
}
}
#endif
#endif
while
(
numReceivedParticleLocations
<
numberOfCreatedParticles
)
{
MPI_Status
status
;
...
...
@@ -1184,16 +1192,16 @@ class MMOCTransport
particleLocationRadius_
=
0.1
*
MeshQuality
::
getMinimalEdgeLength
(
storage
,
maxLevel
);
}
void
step
(
const
FunctionType
&
c
,
const
vecfun_t
&
u
,
const
vecfun_t
&
uLastTimeStep
,
const
uint_t
&
level
,
const
DoFType
&
flag
,
const
real_t
&
dt
,
const
uint_t
&
innerSteps
,
const
bool
&
resetParticles
=
true
,
const
bool
&
globalMaxLimiter
=
true
,
const
bool
&
setParticlesOutsideDomainToZero
=
false
)
void
step
(
const
FunctionType
&
c
,
const
vecfun_t
&
u
,
const
vecfun_t
&
uLastTimeStep
,
const
uint_t
&
level
,
const
DoFType
&
flag
,
const
real_t
&
dt
,
const
uint_t
&
innerSteps
,
const
bool
&
resetParticles
=
true
,
const
bool
&
globalMaxLimiter
=
true
,
const
bool
&
setParticlesOutsideDomainToZero
=
false
)
{
uint_t
aux
=
u
.
getDimension
()
==
3
?
2
:
0
;
step
(
c
,
...
...
@@ -1213,18 +1221,18 @@ class MMOCTransport
}
template
<
typename
MassOperator
>
void
step
(
const
FunctionType
&
c
,
const
vecfun_t
&
u
,
const
vecfun_t
&
uLastTimeStep
,
const
uint_t
&
level
,
const
DoFType
&
flag
,
const
real_t
&
dt
,
const
uint_t
&
innerSteps
,
const
MassOperator
&
massOperator
,
const
real_t
&
allowedRelativeMassDifference
,
const
real_t
&
dtPertubationAdjustedAdvection
,
const
bool
&
globalMaxLimiter
=
true
,
const
bool
&
setParticlesOutsideDomainToZero
=
false
)
void
step
(
const
FunctionType
&
c
,
const
vecfun_t
&
u
,
const
vecfun_t
&
uLastTimeStep
,
const
uint_t
&
level
,
const
DoFType
&
flag
,
const
real_t
&
dt
,
const
uint_t
&
innerSteps
,
const
MassOperator
&
massOperator
,
const
real_t
&
allowedRelativeMassDifference
,
const
real_t
&
dtPertubationAdjustedAdvection
,
const
bool
&
globalMaxLimiter
=
true
,
const
bool
&
setParticlesOutsideDomainToZero
=
false
)
{
uint_t
aux
=
u
.
getDimension
()
==
3
?
2
:
0
;
step
(
c
,
...
...
@@ -1271,7 +1279,7 @@ class MMOCTransport
if
(
resetParticles
)
{
cOld_
.
assign
(
{
1.0
},
{
c
},
level
,
All
);
cOld_
.
assign
(
{
1.0
},
{
c
},
level
,
All
);
storage_
->
getTimingTree
()
->
start
(
"Particle initialization"
);
numberOfCreatedParticles_
=
initializeParticles
(
particleStorage_
,
*
storage_
,
c
,
ux
,
uy
,
uz
,
level
,
Inner
,
timeSteppingSchemeConvection_
,
0
);
...
...
@@ -1335,8 +1343,8 @@ class MMOCTransport
cMinus_
.
copyBoundaryConditionFromFunction
(
c
);
cAdjusted_
.
copyBoundaryConditionFromFunction
(
c
);
cPlus_
.
assign
(
{
1.0
},
{
c
},
level
,
All
);
cMinus_
.
assign
(
{
1.0
},
{
c
},
level
,
All
);
cPlus_
.
assign
(
{
1.0
},
{
c
},
level
,
All
);
cMinus_
.
assign
(
{
1.0
},
{
c
},
level
,
All
);
// calculate old mass
massOperator
.
apply
(
c
,
cTmp_
,
level
,
flag
);
...
...
@@ -1407,11 +1415,11 @@ class MMOCTransport
if
(
massAfter
<=
massBefore
)
{
cAdjusted_
.
interpolate
(
maxAssignment
,
{
cPlus_
,
cMinus_
},
level
);
cAdjusted_
.
interpolate
(
maxAssignment
,
{
cPlus_
,
cMinus_
},
level
);
}
else
{
cAdjusted_
.
interpolate
(
minAssignment
,
{
cPlus_
,
cMinus_
},
level
);
cAdjusted_
.
interpolate
(
minAssignment
,
{
cPlus_
,
cMinus_
},
level
);
}
// calculate adjustment mass
...
...
@@ -1420,7 +1428,7 @@ class MMOCTransport
auto
theta
=
(
massBefore
-
massAdjusted
)
/
(
massAfter
-
massAdjusted
);
c
.
assign
(
{
theta
,
1
-
theta
},
{
c
,
cAdjusted_
},
level
,
flag
);
c
.
assign
(
{
theta
,
1
-
theta
},
{
c
,
cAdjusted_
},
level
,
flag
);
}
const
std
::
shared_ptr
<
PrimitiveStorage
>
storage_
;
...
...
src/hyteg/CMakeLists.txt
View file @
537bdc9b
...
...
@@ -35,6 +35,8 @@ add_subdirectory( boundary )
add_subdirectory
(
composites
)
add_subdirectory
(
solvers
)
add_subdirectory
(
dgfunctionspace
)
add_subdirectory
(
dgfunctionspace_old
)
add_subdirectory
(
volumedofspace
)
add_subdirectory
(
mixedoperators
)
add_subdirectory
(
eigen
)
add_subdirectory
(
elementwiseoperators
)
...
...
@@ -50,6 +52,7 @@ add_subdirectory( petsc )
add_subdirectory
(
geometry
)
add_subdirectory
(
indexing
)
add_subdirectory
(
communication
)
add_subdirectory
(
p0functionspace
)
add_subdirectory
(
p1functionspace
)
add_subdirectory
(
facedofspace
)
add_subdirectory
(
facedofspace
_old
)
add_subdirectory
(
polynomial
)
src/hyteg/Levelinfo.hpp
View file @
537bdc9b
/*
* Copyright (c) 2017-20
19
Dominik Thoennes, Nils Kohl.
* Copyright (c) 2017-20
22
Dominik Thoennes, Nils Kohl.