Skip to content
GitLab
Menu
Projects
Groups
Snippets
Loading...
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
hyteg
hyteg
Commits
56c5fe8b
Commit
56c5fe8b
authored
May 20, 2022
by
wagnandr
Browse files
applies clang-format to GMRES files
parent
d6f649cc
Pipeline
#40117
failed with stages
in 164 minutes and 25 seconds
Changes
3
Pipelines
1
Expand all
Hide whitespace changes
Inline
Side-by-side
apps/GMRESTestApp/GMRESApp.cpp
View file @
56c5fe8b
#include <cmath>
#include "core/DataTypes.h"
#include "core/Environment.h"
#include "core/mpi/MPIManager.h"
#include "hyteg/dataexport/VTKOutput.hpp"
#include "hyteg/gridtransferoperators/P1toP1LinearProlongation.hpp"
#include "hyteg/gridtransferoperators/P1toP1LinearRestriction.hpp"
...
...
@@ -12,69 +12,70 @@
#include "hyteg/primitivestorage/PrimitiveStorage.hpp"
#include "hyteg/primitivestorage/loadbalancing/SimpleBalancer.hpp"
#include "hyteg/solvers/CGSolver.hpp"
#include "hyteg/solvers/GMRESSolver.hpp"
#include "hyteg/solvers/GaussSeidelSmoother.hpp"
#include "hyteg/solvers/GeometricMultigridSolver.hpp"
#include "hyteg/solvers/GMRESSolver.hpp"
#include "hyteg/solvers/preconditioners/JacobiPreconditioner.hpp"
using
walberla
::
real_t
;
using
walberla
::
uint_c
;
using
walberla
::
uint_t
;
int
main
(
int
argc
,
char
**
argv
)
{
walberla
::
Environment
env
(
argc
,
argv
);
walberla
::
mpi
::
MPIManager
::
instance
()
->
useWorldComm
();
uint_t
numProcesses
=
uint_c
(
walberla
::
mpi
::
MPIManager
::
instance
()
->
numProcesses
()
);
walberla
::
shared_ptr
<
walberla
::
config
::
Config
>
cfg
(
new
walberla
::
config
::
Config
);
cfg
->
readParameterFile
(
"./GMRESparam.prm"
);
walberla
::
Config
::
BlockHandle
parameters
=
cfg
->
getOneBlock
(
"Parameters"
);
parameters
.
listParameters
();
const
uint_t
minLevel
=
parameters
.
getParameter
<
uint_t
>
(
"minLevel"
);
const
uint_t
maxLevel
=
parameters
.
getParameter
<
uint_t
>
(
"maxLevel"
);
const
uint_t
maxKrylowDim
=
parameters
.
getParameter
<
uint_t
>
(
"maxKrylowDim"
);
const
uint_t
restartLength
=
parameters
.
getParameter
<
uint_t
>
(
"restartLength"
);
const
real_t
arnoldiTOL
=
parameters
.
getParameter
<
real_t
>
(
"arnoldiTOL"
);
const
real_t
approxTOL
=
parameters
.
getParameter
<
real_t
>
(
"approxTOL"
);
const
real_t
doubleOrthoTOL
=
parameters
.
getParameter
<
real_t
>
(
"doubleOrthoTOL"
);
hyteg
::
MeshInfo
meshInfo
=
hyteg
::
MeshInfo
::
meshRectangle
(
hyteg
::
Point2D
(
{
0.0
,
0.0
}
),
hyteg
::
Point2D
(
{
1.0
,
1.0
}
),
hyteg
::
MeshInfo
::
CRISS
,
20
,
20
);
const
uint_t
minLevel
=
parameters
.
getParameter
<
uint_t
>
(
"minLevel"
);
const
uint_t
maxLevel
=
parameters
.
getParameter
<
uint_t
>
(
"maxLevel"
);
const
uint_t
maxKrylowDim
=
parameters
.
getParameter
<
uint_t
>
(
"maxKrylowDim"
);
const
uint_t
restartLength
=
parameters
.
getParameter
<
uint_t
>
(
"restartLength"
);
const
real_t
arnoldiTOL
=
parameters
.
getParameter
<
real_t
>
(
"arnoldiTOL"
);
const
real_t
approxTOL
=
parameters
.
getParameter
<
real_t
>
(
"approxTOL"
);
const
real_t
doubleOrthoTOL
=
parameters
.
getParameter
<
real_t
>
(
"doubleOrthoTOL"
);
hyteg
::
MeshInfo
meshInfo
=
hyteg
::
MeshInfo
::
meshRectangle
(
hyteg
::
Point2D
(
{
0.0
,
0.0
}
),
hyteg
::
Point2D
(
{
1.0
,
1.0
}
),
hyteg
::
MeshInfo
::
CRISS
,
20
,
20
);
hyteg
::
SetupPrimitiveStorage
setupStorage
(
meshInfo
,
numProcesses
);
hyteg
::
loadbalancing
::
roundRobin
(
setupStorage
,
numProcesses
);
std
::
shared_ptr
<
hyteg
::
PrimitiveStorage
>
storage
=
std
::
make_shared
<
hyteg
::
PrimitiveStorage
>
(
setupStorage
);
std
::
shared_ptr
<
hyteg
::
PrimitiveStorage
>
storage
=
std
::
make_shared
<
hyteg
::
PrimitiveStorage
>
(
setupStorage
);
hyteg
::
P1Function
<
real_t
>
residual
(
"residual"
,
storage
,
minLevel
,
maxLevel
);
hyteg
::
P1Function
<
real_t
>
f
(
"f"
,
storage
,
minLevel
,
maxLevel
);
hyteg
::
P1Function
<
real_t
>
u
(
"u"
,
storage
,
minLevel
,
maxLevel
);
hyteg
::
P1Function
<
real_t
>
laplaceTimesFunction
(
"laplaceTimesFunction"
,
storage
,
minLevel
,
maxLevel
);
std
::
function
<
real_t
(
const
hyteg
::
Point3D
&
)
>
boundaryConditions
=
[](
const
hyteg
::
Point3D
&
x
)
{
if
(
x
[
0
]
<=
0.1
||
x
[
1
]
<=
0.1
||
x
[
0
]
>=
0.9
||
x
[
1
]
>=
0.9
)
if
(
x
[
0
]
<=
0.1
||
x
[
1
]
<=
0.1
||
x
[
0
]
>=
0.9
||
x
[
1
]
>=
0.9
)
{
return
(
0.25
-
std
::
pow
(
x
[
0
]
-
0.5
,
2
)
);
return
(
0.25
-
std
::
pow
(
x
[
0
]
-
0.5
,
2
)
);
}
WALBERLA_ABORT
(
"point is not on the boundary"
);
};
u
.
interpolate
(
boundaryConditions
,
maxLevel
,
hyteg
::
DirichletBoundary
);
hyteg
::
P1ConstantLaplaceOperator
L
(
storage
,
minLevel
,
maxLevel
);
L
.
computeInverseDiagonalOperatorValues
();
auto
preCondi
=
std
::
make_shared
<
hyteg
::
JacobiPreconditioner
<
hyteg
::
P1ConstantLaplaceOperator
>
>
(
storage
,
minLevel
,
maxLevel
,
2
);
hyteg
::
GMRESSolver
gmresSolver
=
hyteg
::
GMRESSolver
<
hyteg
::
P1ConstantLaplaceOperator
>
(
storage
,
minLevel
,
maxLevel
,
maxKrylowDim
,
restartLength
,
arnoldiTOL
,
approxTOL
,
doubleOrthoTOL
,
preCondi
);
auto
preCondi
=
std
::
make_shared
<
hyteg
::
JacobiPreconditioner
<
hyteg
::
P1ConstantLaplaceOperator
>
>
(
storage
,
minLevel
,
maxLevel
,
2
);
hyteg
::
GMRESSolver
gmresSolver
=
hyteg
::
GMRESSolver
<
hyteg
::
P1ConstantLaplaceOperator
>
(
storage
,
minLevel
,
maxLevel
,
maxKrylowDim
,
restartLength
,
arnoldiTOL
,
approxTOL
,
doubleOrthoTOL
,
preCondi
);
gmresSolver
.
solve
(
L
,
u
,
f
,
maxLevel
);
L
.
apply
(
u
,
laplaceTimesFunction
,
maxLevel
,
hyteg
::
Inner
);
residual
.
assign
(
{
1.0
,
-
1.0
},
{
f
,
laplaceTimesFunction
},
maxLevel
,
hyteg
::
Inner
);
residual
.
assign
(
{
1.0
,
-
1.0
},
{
f
,
laplaceTimesFunction
},
maxLevel
,
hyteg
::
Inner
);
real_t
residualEuclideanNorm
=
std
::
sqrt
(
residual
.
dotGlobal
(
residual
,
maxLevel
,
hyteg
::
Inner
)
);
WALBERLA_LOG_INFO_ON_ROOT
(
"Euclidean norm of residual: "
<<
residualEuclideanNorm
);
if
(
parameters
.
getParameter
<
bool
>
(
"vtkOutput"
)
)
if
(
parameters
.
getParameter
<
bool
>
(
"vtkOutput"
)
)
{
hyteg
::
VTKOutput
vtkOutput
(
"."
,
"GMRESApp"
,
storage
);
vtkOutput
.
add
(
u
);
...
...
apps/GMRESTestApp/StokesGMRES.cpp
View file @
56c5fe8b
...
...
@@ -24,8 +24,8 @@
#include "core/config/Config.h"
#include "core/mpi/MPIManager.h"
#include "hyteg/composites/P1StokesFunction.hpp"
#include "hyteg/composites/P1P1StokesOperator.hpp"
#include "hyteg/composites/P1StokesFunction.hpp"
#include "hyteg/dataexport/VTKOutput.hpp"
#include "hyteg/functions/FunctionProperties.hpp"
#include "hyteg/gridtransferoperators/P1P1StokesToP1P1StokesProlongation.hpp"
...
...
@@ -42,6 +42,7 @@
#include "hyteg/primitivestorage/loadbalancing/DistributedBalancer.hpp"
#include "hyteg/primitivestorage/loadbalancing/SimpleBalancer.hpp"
#include "hyteg/solvers/CGSolver.hpp"
#include "hyteg/solvers/GMRESSolver.hpp"
#include "hyteg/solvers/GaussSeidelSmoother.hpp"
#include "hyteg/solvers/GeometricMultigridSolver.hpp"
#include "hyteg/solvers/MinresSolver.hpp"
...
...
@@ -51,8 +52,6 @@
#include "hyteg/solvers/preconditioners/stokes/StokesPressureBlockPreconditioner.hpp"
#include "hyteg/solvers/preconditioners/stokes/StokesVelocityBlockBlockDiagonalPreconditioner.hpp"
#include "hyteg/solvers/GMRESSolver.hpp"
using
walberla
::
real_c
;
using
walberla
::
real_t
;
using
namespace
hyteg
;
...
...
@@ -63,12 +62,13 @@ int main( int argc, char* argv[] )
walberla
::
MPIManager
::
instance
()
->
useWorldComm
();
auto
cfg
=
std
::
make_shared
<
walberla
::
config
::
Config
>
();
if
(
env
.
config
()
==
nullptr
)
if
(
env
.
config
()
==
nullptr
)
{
auto
defaultFile
=
"./StokesGMRES.prm"
;
WALBERLA_LOG_INFO_ON_ROOT
(
"No Parameter file given loading default parameter file: "
<<
defaultFile
);
cfg
->
readParameterFile
(
defaultFile
);
}
else
}
else
{
cfg
=
env
.
config
();
}
...
...
@@ -78,7 +78,7 @@ int main( int argc, char* argv[] )
const
uint_t
ntan
=
mainConf
.
getParameter
<
uint_t
>
(
"ntan"
);
std
::
vector
<
double
>
layers
;
for
(
auto
it
:
layersParam
)
for
(
auto
it
:
layersParam
)
{
layers
.
push_back
(
layersParam
.
getParameter
<
double
>
(
it
.
first
)
);
}
...
...
@@ -86,20 +86,21 @@ int main( int argc, char* argv[] )
const
double
rmin
=
layers
.
front
();
const
double
rmax
=
layers
.
back
();
const
Point3D
sourcePoint
=
Point3D
(
{
rmin
,
0
,
0
}
)
+
0.5
*
Point3D
(
{
rmax
-
rmin
,
0
,
0
}
);
const
Point3D
sourcePoint
=
Point3D
(
{
rmin
,
0
,
0
}
)
+
0.5
*
Point3D
(
{
rmax
-
rmin
,
0
,
0
}
);
const
real_t
sourceRadius
=
0.5
;
const
uint_t
minLevel
=
mainConf
.
getParameter
<
uint_t
>
(
"minLevel"
);
const
uint_t
maxLevel
=
mainConf
.
getParameter
<
uint_t
>
(
"maxLevel"
);
const
uint_t
numVCycle
=
mainConf
.
getParameter
<
uint_t
>
(
"numVCycle"
);
const
uint_t
maxKrylowDim
=
mainConf
.
getParameter
<
uint_t
>
(
"maxKrylowDim"
);
const
uint_t
restartLength
=
mainConf
.
getParameter
<
uint_t
>
(
"restartLength"
);
const
real_t
arnoldiTOL
=
mainConf
.
getParameter
<
real_t
>
(
"arnoldiTOL"
);
const
real_t
approxTOL
=
mainConf
.
getParameter
<
real_t
>
(
"approxTOL"
);
const
real_t
doubleOrthoTOL
=
mainConf
.
getParameter
<
real_t
>
(
"doubleOrthoTOL"
);
const
uint_t
minLevel
=
mainConf
.
getParameter
<
uint_t
>
(
"minLevel"
);
const
uint_t
maxLevel
=
mainConf
.
getParameter
<
uint_t
>
(
"maxLevel"
);
const
uint_t
numVCycle
=
mainConf
.
getParameter
<
uint_t
>
(
"numVCycle"
);
const
uint_t
maxKrylowDim
=
mainConf
.
getParameter
<
uint_t
>
(
"maxKrylowDim"
);
const
uint_t
restartLength
=
mainConf
.
getParameter
<
uint_t
>
(
"restartLength"
);
const
real_t
arnoldiTOL
=
mainConf
.
getParameter
<
real_t
>
(
"arnoldiTOL"
);
const
real_t
approxTOL
=
mainConf
.
getParameter
<
real_t
>
(
"approxTOL"
);
const
real_t
doubleOrthoTOL
=
mainConf
.
getParameter
<
real_t
>
(
"doubleOrthoTOL"
);
hyteg
::
MeshInfo
meshInfo
=
hyteg
::
MeshInfo
::
meshSphericalShell
(
ntan
,
layers
);
hyteg
::
SetupPrimitiveStorage
setupStorage
(
meshInfo
,
walberla
::
uint_c
(
walberla
::
mpi
::
MPIManager
::
instance
()
->
numProcesses
()
)
);
hyteg
::
SetupPrimitiveStorage
setupStorage
(
meshInfo
,
walberla
::
uint_c
(
walberla
::
mpi
::
MPIManager
::
instance
()
->
numProcesses
()
)
);
hyteg
::
loadbalancing
::
roundRobin
(
setupStorage
);
setupStorage
.
setMeshBoundaryFlagsOnBoundary
(
1
,
0
,
true
);
...
...
@@ -112,10 +113,10 @@ int main( int argc, char* argv[] )
hyteg
::
P1StokesFunction
<
real_t
>
f
(
"f"
,
storage
,
minLevel
,
maxLevel
);
hyteg
::
P1StokesFunction
<
real_t
>
u
(
"u"
,
storage
,
minLevel
,
maxLevel
);
if
(
mainConf
.
getParameter
<
bool
>
(
"printDoFCount"
)
)
if
(
mainConf
.
getParameter
<
bool
>
(
"printDoFCount"
)
)
{
uint_t
totalGlobalDofsStokes
=
0
;
for
(
uint_t
lvl
=
minLevel
;
lvl
<=
maxLevel
;
++
lvl
)
for
(
uint_t
lvl
=
minLevel
;
lvl
<=
maxLevel
;
++
lvl
)
{
uint_t
tmpDofStokes
=
numberOfGlobalDoFs
<
hyteg
::
P1StokesFunctionTag
>
(
*
storage
,
lvl
);
WALBERLA_LOG_INFO_ON_ROOT
(
"Stokes DoFs on level "
<<
lvl
<<
" : "
<<
tmpDofStokes
);
...
...
@@ -124,8 +125,8 @@ int main( int argc, char* argv[] )
WALBERLA_LOG_INFO_ON_ROOT
(
"Total Stokes DoFs on all level :"
<<
totalGlobalDofsStokes
);
}
hyteg
::
VTKOutput
vtkOutput
(
"./output"
,
"StokesGMRES"
,
storage
);
if
(
mainConf
.
getParameter
<
bool
>
(
"VTKOutput"
)
)
hyteg
::
VTKOutput
vtkOutput
(
"./output"
,
"StokesGMRES"
,
storage
);
if
(
mainConf
.
getParameter
<
bool
>
(
"VTKOutput"
)
)
{
vtkOutput
.
add
(
u
.
uvw
()
);
vtkOutput
.
add
(
u
.
p
()
);
...
...
@@ -137,7 +138,7 @@ int main( int argc, char* argv[] )
std
::
function
<
real_t
(
const
hyteg
::
Point3D
&
)
>
rhsPlumeX
=
[
sourcePoint
,
sourceRadius
](
const
hyteg
::
Point3D
&
x
)
{
const
real_t
distToSourcePoint
=
(
x
-
sourcePoint
).
norm
();
if
(
distToSourcePoint
<
sourceRadius
)
if
(
distToSourcePoint
<
sourceRadius
)
return
x
[
0
]
*
(
sourceRadius
-
distToSourcePoint
);
else
return
0.0
;
...
...
@@ -145,7 +146,7 @@ int main( int argc, char* argv[] )
std
::
function
<
real_t
(
const
hyteg
::
Point3D
&
)
>
rhsPlumeY
=
[
sourcePoint
,
sourceRadius
](
const
hyteg
::
Point3D
&
x
)
{
const
real_t
distToSourcePoint
=
(
x
-
sourcePoint
).
norm
();
if
(
distToSourcePoint
<
sourceRadius
)
if
(
distToSourcePoint
<
sourceRadius
)
return
x
[
1
]
*
(
sourceRadius
-
distToSourcePoint
);
else
return
0.0
;
...
...
@@ -153,7 +154,7 @@ int main( int argc, char* argv[] )
std
::
function
<
real_t
(
const
hyteg
::
Point3D
&
)
>
rhsPlumeZ
=
[
sourcePoint
,
sourceRadius
](
const
hyteg
::
Point3D
&
x
)
{
const
real_t
distToSourcePoint
=
(
x
-
sourcePoint
).
norm
();
if
(
distToSourcePoint
<
sourceRadius
)
if
(
distToSourcePoint
<
sourceRadius
)
return
x
[
2
]
*
(
sourceRadius
-
distToSourcePoint
);
else
return
0.0
;
...
...
@@ -162,117 +163,138 @@ int main( int argc, char* argv[] )
std
::
function
<
real_t
(
const
hyteg
::
Point3D
&
)
>
zero
=
[](
const
hyteg
::
Point3D
&
)
{
return
0.0
;
};
std
::
function
<
real_t
(
const
hyteg
::
Point3D
&
)
>
ones
=
[](
const
hyteg
::
Point3D
&
)
{
return
1.0
;
};
f
.
uvw
().
interpolate
(
{
rhsPlumeX
,
rhsPlumeY
,
rhsPlumeZ
},
maxLevel
);
f
.
uvw
().
interpolate
(
{
rhsPlumeX
,
rhsPlumeY
,
rhsPlumeZ
},
maxLevel
);
if
(
mainConf
.
getParameter
<
bool
>
(
"VTKOutput"
)
)
if
(
mainConf
.
getParameter
<
bool
>
(
"VTKOutput"
)
)
{
vtkOutput
.
write
(
maxLevel
,
0
);
}
std
::
string
solverType
=
mainConf
.
getParameter
<
std
::
string
>
(
"solver"
);
if
(
solverType
==
"plainGMRES"
)
if
(
solverType
==
"plainGMRES"
)
{
typedef
hyteg
::
GMRESSolver
<
hyteg
::
P1P1StokesOperator
>
GMRESsolv
;
auto
KleinerFeinerGMRES
=
GMRESsolv
(
storage
,
minLevel
,
maxLevel
,
maxKrylowDim
,
restartLength
,
arnoldiTOL
,
approxTOL
,
doubleOrthoTOL
);
typedef
hyteg
::
GMRESSolver
<
hyteg
::
P1P1StokesOperator
>
GMRESsolv
;
auto
KleinerFeinerGMRES
=
GMRESsolv
(
storage
,
minLevel
,
maxLevel
,
maxKrylowDim
,
restartLength
,
arnoldiTOL
,
approxTOL
,
doubleOrthoTOL
);
KleinerFeinerGMRES
.
solve
(
L
,
u
,
f
,
maxLevel
);
L
.
apply
(
u
,
r
,
maxLevel
,
hyteg
::
Inner
|
hyteg
::
NeumannBoundary
);
r
.
assign
(
{
1.0
,
-
1.0
},
{
f
,
r
},
maxLevel
,
hyteg
::
Inner
|
hyteg
::
NeumannBoundary
);
r
.
assign
(
{
1.0
,
-
1.0
},
{
f
,
r
},
maxLevel
,
hyteg
::
Inner
|
hyteg
::
NeumannBoundary
);
real_t
currentResidualL2
=
sqrt
(
r
.
dotGlobal
(
r
,
maxLevel
,
hyteg
::
Inner
)
);
WALBERLA_LOG_INFO_ON_ROOT
(
"currentResidualL2 : "
<<
currentResidualL2
);
}
else
if
(
solverType
==
"simplePrecGMRES"
)
WALBERLA_LOG_INFO_ON_ROOT
(
"currentResidualL2 : "
<<
currentResidualL2
);
}
else
if
(
solverType
==
"simplePrecGMRES"
)
{
/// A block Preconditioner for GMRES /////
typedef
StokesBlockDiagonalPreconditioner
<
hyteg
::
P1P1StokesOperator
,
hyteg
::
P1LumpedInvMassOperator
>
Preconditioner_T
;
typedef
StokesBlockDiagonalPreconditioner
<
hyteg
::
P1P1StokesOperator
,
hyteg
::
P1LumpedInvMassOperator
>
Preconditioner_T
;
auto
prec
=
std
::
make_shared
<
Preconditioner_T
>
(
storage
,
minLevel
,
maxLevel
,
5
);
typedef
hyteg
::
GMRESSolver
<
hyteg
::
P1P1StokesOperator
>
GMRESsolv
;
auto
KleinerFeinerGMRES
=
GMRESsolv
(
storage
,
minLevel
,
maxLevel
,
maxKrylowDim
,
restartLength
,
arnoldiTOL
,
approxTOL
,
doubleOrthoTOL
,
prec
);
typedef
hyteg
::
GMRESSolver
<
hyteg
::
P1P1StokesOperator
>
GMRESsolv
;
auto
KleinerFeinerGMRES
=
GMRESsolv
(
storage
,
minLevel
,
maxLevel
,
maxKrylowDim
,
restartLength
,
arnoldiTOL
,
approxTOL
,
doubleOrthoTOL
,
prec
);
KleinerFeinerGMRES
.
solve
(
L
,
u
,
f
,
maxLevel
);
L
.
apply
(
u
,
r
,
maxLevel
,
hyteg
::
Inner
|
hyteg
::
NeumannBoundary
);
r
.
assign
(
{
1.0
,
-
1.0
},
{
f
,
r
},
maxLevel
,
hyteg
::
Inner
|
hyteg
::
NeumannBoundary
);
r
.
assign
(
{
1.0
,
-
1.0
},
{
f
,
r
},
maxLevel
,
hyteg
::
Inner
|
hyteg
::
NeumannBoundary
);
real_t
currentResidualL2
=
sqrt
(
r
.
dotGlobal
(
r
,
maxLevel
,
hyteg
::
Inner
)
);
WALBERLA_LOG_INFO_ON_ROOT
(
"currentResidualL2 : "
<<
currentResidualL2
);
}
else
if
(
solverType
==
"plainMINRES"
)
WALBERLA_LOG_INFO_ON_ROOT
(
"currentResidualL2 : "
<<
currentResidualL2
);
}
else
if
(
solverType
==
"plainMINRES"
)
{
typedef
hyteg
::
MinResSolver
<
hyteg
::
P1P1StokesOperator
>
MinResSolv
;
typedef
hyteg
::
MinResSolver
<
hyteg
::
P1P1StokesOperator
>
MinResSolv
;
auto
KleinerFeinerMinRes
=
MinResSolv
(
storage
,
minLevel
,
maxLevel
,
maxKrylowDim
,
approxTOL
);
KleinerFeinerMinRes
.
solve
(
L
,
u
,
f
,
maxLevel
);
L
.
apply
(
u
,
r
,
maxLevel
,
hyteg
::
Inner
|
hyteg
::
NeumannBoundary
);
r
.
assign
(
{
1.0
,
-
1.0
},
{
f
,
r
},
maxLevel
,
hyteg
::
Inner
|
hyteg
::
NeumannBoundary
);
r
.
assign
(
{
1.0
,
-
1.0
},
{
f
,
r
},
maxLevel
,
hyteg
::
Inner
|
hyteg
::
NeumannBoundary
);
real_t
currentResidualL2
=
sqrt
(
r
.
dotGlobal
(
r
,
maxLevel
,
hyteg
::
Inner
)
);
WALBERLA_LOG_INFO_ON_ROOT
(
"currentResidualL2 : "
<<
currentResidualL2
);
}
else
if
(
solverType
==
"simplePrecMINRES"
)
WALBERLA_LOG_INFO_ON_ROOT
(
"currentResidualL2 : "
<<
currentResidualL2
);
}
else
if
(
solverType
==
"simplePrecMINRES"
)
{
/// A block Preconditioner for MinRes /////
typedef
StokesBlockDiagonalPreconditioner
<
hyteg
::
P1P1StokesOperator
,
hyteg
::
P1LumpedInvMassOperator
>
Preconditioner_T
;
typedef
StokesBlockDiagonalPreconditioner
<
hyteg
::
P1P1StokesOperator
,
hyteg
::
P1LumpedInvMassOperator
>
Preconditioner_T
;
auto
prec
=
std
::
make_shared
<
Preconditioner_T
>
(
storage
,
minLevel
,
maxLevel
,
5
);
typedef
hyteg
::
MinResSolver
<
hyteg
::
P1P1StokesOperator
>
MinResSolv
;
typedef
hyteg
::
MinResSolver
<
hyteg
::
P1P1StokesOperator
>
MinResSolv
;
auto
KleinerFeinerMinRes
=
MinResSolv
(
storage
,
minLevel
,
maxLevel
,
maxKrylowDim
,
approxTOL
,
prec
);
KleinerFeinerMinRes
.
solve
(
L
,
u
,
f
,
maxLevel
);
L
.
apply
(
u
,
r
,
maxLevel
,
hyteg
::
Inner
|
hyteg
::
NeumannBoundary
);
r
.
assign
(
{
1.0
,
-
1.0
},
{
f
,
r
},
maxLevel
,
hyteg
::
Inner
|
hyteg
::
NeumannBoundary
);
r
.
assign
(
{
1.0
,
-
1.0
},
{
f
,
r
},
maxLevel
,
hyteg
::
Inner
|
hyteg
::
NeumannBoundary
);
real_t
currentResidualL2
=
sqrt
(
r
.
dotGlobal
(
r
,
maxLevel
,
hyteg
::
Inner
)
);
WALBERLA_LOG_INFO_ON_ROOT
(
"currentResidualL2 : "
<<
currentResidualL2
);
}
else
if
(
solverType
==
"multigrid"
)
WALBERLA_LOG_INFO_ON_ROOT
(
"currentResidualL2 : "
<<
currentResidualL2
);
}
else
if
(
solverType
==
"multigrid"
)
{
typedef
StokesPressureBlockPreconditioner
<
hyteg
::
P1P1StokesOperator
,
hyteg
::
P1LumpedInvMassOperator
>
PressurePreconditioner_T
;
auto
pressurePrec
=
std
::
make_shared
<
PressurePreconditioner_T
>
(
storage
,
minLevel
,
minLevel
);
typedef
StokesPressureBlockPreconditioner
<
hyteg
::
P1P1StokesOperator
,
hyteg
::
P1LumpedInvMassOperator
>
PressurePreconditioner_T
;
auto
pressurePrec
=
std
::
make_shared
<
PressurePreconditioner_T
>
(
storage
,
minLevel
,
minLevel
);
typedef
hyteg
::
GMRESSolver
<
hyteg
::
P1P1StokesOperator
>
PressurePreconditionedGMRES_T
;
typedef
hyteg
::
GMRESSolver
<
hyteg
::
P1P1StokesOperator
>
PressurePreconditionedGMRES_T
;
auto
pressurePreconditionedGMRESSolver
=
std
::
make_shared
<
PressurePreconditionedGMRES_T
>
(
storage
,
minLevel
,
maxLevel
,
maxKrylowDim
,
restartLength
,
arnoldiTOL
,
approxTOL
,
doubleOrthoTOL
,
pressurePrec
);
typedef
GeometricMultigridSolver
<
hyteg
::
P1P1StokesOperator
>
GMRESMultigrid_T
;
typedef
GeometricMultigridSolver
<
hyteg
::
P1P1StokesOperator
>
GMRESMultigrid_T
;
auto
stokesRestriction
=
std
::
make_shared
<
hyteg
::
P1P1StokesToP1P1StokesRestriction
>
();
auto
stokesRestriction
=
std
::
make_shared
<
hyteg
::
P1P1StokesToP1P1StokesRestriction
>
();
auto
stokesProlongation
=
std
::
make_shared
<
hyteg
::
P1P1StokesToP1P1StokesProlongation
>
();
auto
gaussSeidel
=
std
::
make_shared
<
hyteg
::
GaussSeidelSmoother
<
hyteg
::
P1P1StokesOperator
::
VelocityOperator_T
>
>
();
auto
multigridVelocityPreconditioner
=
std
::
make_shared
<
hyteg
::
StokesVelocityBlockBlockDiagonalPreconditioner
<
hyteg
::
P1P1StokesOperator
>
>
(
storage
,
gaussSeidel
);
auto
GMRESMultigridSmoother
=
std
::
make_shared
<
hyteg
::
UzawaSmoother
<
P1P1StokesOperator
>
>
(
storage
,
multigridVelocityPreconditioner
,
minLevel
,
maxLevel
,
0.3
);
GMRESMultigrid_T
GMRESsolv
(
storage
,
GMRESMultigridSmoother
,
pressurePreconditionedGMRESSolver
,
stokesRestriction
,
stokesProlongation
,
minLevel
,
maxLevel
,
2
,
2
,
2
);
auto
gaussSeidel
=
std
::
make_shared
<
hyteg
::
GaussSeidelSmoother
<
hyteg
::
P1P1StokesOperator
::
VelocityOperator_T
>
>
();
auto
multigridVelocityPreconditioner
=
std
::
make_shared
<
hyteg
::
StokesVelocityBlockBlockDiagonalPreconditioner
<
hyteg
::
P1P1StokesOperator
>
>
(
storage
,
gaussSeidel
);
auto
GMRESMultigridSmoother
=
std
::
make_shared
<
hyteg
::
UzawaSmoother
<
P1P1StokesOperator
>
>
(
storage
,
multigridVelocityPreconditioner
,
minLevel
,
maxLevel
,
0.3
);
GMRESMultigrid_T
GMRESsolv
(
storage
,
GMRESMultigridSmoother
,
pressurePreconditionedGMRESSolver
,
stokesRestriction
,
stokesProlongation
,
minLevel
,
maxLevel
,
2
,
2
,
2
);
auto
count
=
hyteg
::
Function
<
hyteg
::
vertexdof
::
VertexDoFFunction
<
real_t
>
>::
getLevelWiseFunctionCounter
();
if
(
mainConf
.
getParameter
<
bool
>
(
"printFunctionCount"
)
)
{
for
(
uint_t
i
=
minLevel
;
i
<=
maxLevel
;
++
i
)
{
WALBERLA_LOG_INFO_ON_ROOT
(
"Total number of P1 Functions on "
<<
i
<<
" : "
<<
count
[
i
]);
if
(
mainConf
.
getParameter
<
bool
>
(
"printFunctionCount"
)
)
{
for
(
uint_t
i
=
minLevel
;
i
<=
maxLevel
;
++
i
)
{
WALBERLA_LOG_INFO_ON_ROOT
(
"Total number of P1 Functions on "
<<
i
<<
" : "
<<
count
[
i
]
);
}
}
L
.
apply
(
u
,
r
,
maxLevel
,
hyteg
::
Inner
|
hyteg
::
NeumannBoundary
);
r
.
assign
(
{
1.0
,
-
1.0
},
{
f
,
r
},
maxLevel
,
hyteg
::
Inner
|
hyteg
::
NeumannBoundary
);
r
.
assign
(
{
1.0
,
-
1.0
},
{
f
,
r
},
maxLevel
,
hyteg
::
Inner
|
hyteg
::
NeumannBoundary
);
real_t
currentResidualL2
=
sqrt
(
r
.
dotGlobal
(
r
,
maxLevel
,
hyteg
::
Inner
)
)
/
real_c
(
hyteg
::
numberOfGlobalDoFs
<
hyteg
::
P1StokesFunctionTag
>
(
*
storage
,
maxLevel
)
);
real_t
lastResidualL2
=
currentResidualL2
;
WALBERLA_LOG_INFO_ON_ROOT
(
"[StokesSphere] iteration | residual (L2) | convergence rate "
);
WALBERLA_LOG_INFO_ON_ROOT
(
"[StokesSphere] ----------+---------------+------------------"
);
WALBERLA_LOG_INFO_ON_ROOT
(
"[StokesSphere] "
<<
std
::
setw
(
9
)
<<
0
<<
" | "
<<
std
::
setw
(
1
3
)
<<
std
::
scientific
<<
currentResidualL2
<<
" | "
<<
std
::
setw
(
16
)
<<
std
::
scientific
<<
currentResidualL2
/
lastResidualL2
);
for
(
uint_t
i
=
0
;
i
<
numVCycle
;
i
++
)
WALBERLA_LOG_INFO_ON_ROOT
(
"[StokesSphere] "
<<
std
::
setw
(
9
)
<<
0
<<
" | "
<<
std
::
setw
(
13
)
<<
std
::
scientific
<<
currentResidualL2
<<
" | "
<<
std
::
setw
(
1
6
)
<<
std
::
scientific
<<
currentResidualL2
/
lastResidualL2
);
for
(
uint_t
i
=
0
;
i
<
numVCycle
;
i
++
)
{
GMRESsolv
.
solve
(
L
,
u
,
f
,
maxLevel
);
lastResidualL2
=
currentResidualL2
;
L
.
apply
(
u
,
r
,
maxLevel
,
hyteg
::
Inner
|
hyteg
::
NeumannBoundary
);
r
.
assign
(
{
1.0
,
-
1.0
},
{
f
,
r
},
maxLevel
,
hyteg
::
Inner
|
hyteg
::
NeumannBoundary
);
r
.
assign
(
{
1.0
,
-
1.0
},
{
f
,
r
},
maxLevel
,
hyteg
::
Inner
|
hyteg
::
NeumannBoundary
);
currentResidualL2
=
sqrt
(
r
.
dotGlobal
(
r
,
maxLevel
,
hyteg
::
Inner
)
)
/
real_c
(
hyteg
::
numberOfGlobalDoFs
<
hyteg
::
P1StokesFunctionTag
>
(
*
storage
,
maxLevel
)
);
WALBERLA_LOG_INFO_ON_ROOT
(
"[StokesSphere] "
<<
std
::
setw
(
9
)
<<
i
+
1
<<
" | "
<<
std
::
setw
(
1
3
)
<<
std
::
scientific
<<
currentResidualL2
<<
" | "
<<
std
::
setw
(
16
)
<<
std
::
scientific
<<
currentResidualL2
/
lastResidualL2
)
WALBERLA_LOG_INFO_ON_ROOT
(
"[StokesSphere] "
<<
std
::
setw
(
9
)
<<
i
+
1
<<
" | "
<<
std
::
setw
(
13
)
<<
std
::
scientific
<<
currentResidualL2
<<
" | "
<<
std
::
setw
(
1
6
)
<<
std
::
scientific
<<
currentResidualL2
/
lastResidualL2
)
//WALBERLA_LOG_INFO_ON_ROOT( "after it " << i << ": " << std::scientific << residualMG );
}
}
else
}
else
{
WALBERLA_LOG_INFO_ON_ROOT
(
"Unkown solver type"
);
}
...
...
@@ -288,16 +310,17 @@ int main( int argc, char* argv[] )
f.w.assign( {1.0}, {&u.w}, level, DirichletBoundary );
petScLUSolver.solve( L, u, f, r, level, uzawaTolerance, maxIterations, Inner | NeumannBoundary );
#endif
if
(
mainConf
.
getParameter
<
bool
>
(
"VTKOutput"
)
)
if
(
mainConf
.
getParameter
<
bool
>
(
"VTKOutput"
)
)
{
vtkOutput
.
write
(
maxLevel
,
1
);
}
if
(
mainConf
.
getParameter
<
bool
>
(
"PrintTiming"
)
)
{
if
(
mainConf
.
getParameter
<
bool
>
(
"PrintTiming"
)
)
{
auto
tt
=
timingTree
->
getReduced
();
//19.07.2018 this is not in walberla master yet
//auto tt = timingTree->getCopyWithRemainder();
WALBERLA_LOG_INFO_ON_ROOT
(
tt
);
WALBERLA_LOG_INFO_ON_ROOT
(
tt
);
}
return
EXIT_SUCCESS
;
...
...
src/hyteg/solvers/GMRESSolver.hpp
View file @
56c5fe8b
This diff is collapsed.
Click to expand it.
Write
Preview
Markdown
is supported
0%
Try again
or
attach a new file
.
Attach a file
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment