Skip to content
Commits on Source (8)
Hide whitespace changes
Inline Side-by-side
src/petsc4Foam/Make/options
View file @ 8ebbdb63
......@@ -13,6 +13,7 @@ EXE_INC = \
-Wno-old-style-cast \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/Pstream/mpi \
$(foreach dir,$(PETSC_INC_DIR),-I$(dir))
LIB_LIBS = \
......
src/petsc4Foam/solvers/petscSolver.C
View file @ 8ebbdb63
......@@ -33,6 +33,7 @@ License
#include "cyclicLduInterface.H"
#include "cyclicAMILduInterface.H"
#include "addToRunTimeSelectionTable.H"
#include "PstreamGlobals.H"
#include "petscSolver.H"
#include "petscControls.H"
......@@ -165,14 +166,24 @@ Foam::solverPerformance Foam::petscSolver::scalarSolve
ctx.initialized() = true;
needsCacheUpdate = false;
PetscLogStageRegister
(
("foam_" + eqName_ + "_mat").c_str(),
&ctx.matstage
);
PetscLogStageRegister
(
("foam_" + eqName_ + "_ksp").c_str(),
&ctx.kspstage
);
computeMatAllocation(Amat, lowNonZero, maxLowNonZeroPerRow);
buildMat(Amat, lowNonZero, maxLowNonZeroPerRow);
buildMat(Amat, lowNonZero, maxLowNonZeroPerRow, ctx.matstage);
buildKsp(Amat, ksp);
}
if (ctx.caching.needsMatrixUpdate() && needsCacheUpdate)
{
buildMat(Amat, lowNonZero, maxLowNonZeroPerRow);
buildMat(Amat, lowNonZero, maxLowNonZeroPerRow, ctx.matstage);
updateKsp(ksp, Amat, !ctx.caching.needsPrecondUpdate());
}
......@@ -219,8 +230,6 @@ Foam::solverPerformance Foam::petscSolver::scalarSolve
fieldName_
);
ctx.solverPerf = &solverPerf;
// Create solution and rhs vectors for PETSc
// We could create these once, and call PlaceArray/ResetArray instead
PetscWrappedVector petsc_psi(psi, Amat);
......@@ -270,6 +279,8 @@ Foam::solverPerformance Foam::petscSolver::scalarSolve
);
}
PetscLogStagePush(ctx.kspstage);
// Setup KSP (this is not explicitly needed, but we call it to separate
// PCSetUp from KSPSolve timings when requesting -log_view from PETSc)
KSPSetUp(ksp);
......@@ -277,6 +288,8 @@ Foam::solverPerformance Foam::petscSolver::scalarSolve
// Solve A x = b
KSPSolve(ksp, petsc_source, petsc_psi);
PetscLogStagePop();
// Set nIterations and final residual
PetscInt nIters;
KSPGetIterationNumber(ksp, &nIters);
......@@ -289,6 +302,9 @@ Foam::solverPerformance Foam::petscSolver::scalarSolve
solverPerf.finalResidual() = rnorm;
}
// Retain copy of solverPerformance
ctx.performance = solverPerf;
ctx.caching.eventEnd();
// Return solver performance to OpenFOAM
......@@ -354,7 +370,7 @@ void Foam::petscSolver::buildKsp
) const
{
// Create parallel solver context
KSPCreate(PETSC_COMM_WORLD, &ksp);
KSPCreate(PetscObjectComm((PetscObject)Amat), &ksp);
// Set the prefix for the options db (e.g. -eqn_p_)
KSPSetOptionsPrefix(ksp, prefix_.c_str());
......@@ -384,6 +400,17 @@ void Foam::petscSolver::computeMatAllocation
label& maxLowNonZeroPerRow
) const
{
// communicator for the matrix
// PETSc internally duplicates the communicator to
// avoid tag/attributes clashing
const auto comm =
(
PstreamGlobals::MPICommunicators_.empty()
? PETSC_COMM_SELF
: PstreamGlobals::MPICommunicators_[matrix_.mesh().comm()]
);
// traversal-based algorithm
typedef List<PetscInt> integerList;
......@@ -398,7 +425,7 @@ void Foam::petscSolver::computeMatAllocation
const label nrows_ = lduAddr.size();
// Create a petsc matrix
MatCreate(PETSC_COMM_WORLD, &Amat);
MatCreate(comm, &Amat);
MatSetSizes(Amat, nrows_, nrows_, PETSC_DECIDE, PETSC_DECIDE);
// Set the prefix for the options db (e.g. -eqn_p_)
......@@ -514,9 +541,12 @@ void Foam::petscSolver::buildMat
(
Mat& Amat,
List<label>& lowNonZero,
label& maxLowNonZeroPerRow
label& maxLowNonZeroPerRow,
PetscLogStage stage
) const
{
PetscLogStagePush(stage);
const lduAddressing& lduAddr = matrix_.mesh().lduAddr();
const lduInterfacePtrsList interfaces(matrix_.mesh().interfaces());
......@@ -697,6 +727,7 @@ void Foam::petscSolver::buildMat
MatAssemblyBegin(Amat, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(Amat, MAT_FINAL_ASSEMBLY);
PetscLogStagePop();
}
......
src/petsc4Foam/solvers/petscSolver.H
View file @ 8ebbdb63
......@@ -36,19 +36,16 @@ Description
For example,
\verbatim
solver petsc;
preconditioner bjacobi;
solver petsc;
petsc
{
solver cg;
sub_preconditioner ilu;
use_petsc_residual_norm true;
cgCoeffs
options
{
single_reduction true;
ksp_type cg;
ksp_cg_single_reduction true;
}
}
\endverbatim
......@@ -108,7 +105,8 @@ class petscSolver
(
Mat& Amat,
List<label>& lowNonZero,
label& maxLowNonZeroPerRow
label& maxLowNonZeroPerRow,
PetscLogStage stage
) const;
//- Build the ksp_ krylov solver context
......
src/petsc4Foam/utils/petscCacheManager.H
View file @ 8ebbdb63
......@@ -296,8 +296,8 @@ private:
}
}
// Default: Never
return false;
// Default: Always update
return true;
}
};
......
src/petsc4Foam/utils/petscLinearSolverContext.H
View file @ 8ebbdb63
......@@ -69,7 +69,10 @@ public:
Mat Amat;
KSP ksp;
solverPerformance* solverPerf;
PetscLogStage matstage;
PetscLogStage kspstage;
solverPerformance performance;
petscCacheManager caching;
List<PetscInt> lowNonZero;
......@@ -143,7 +146,8 @@ public:
}
}
//- Compute the normFactor used in convergence testing, assumes ArowsSum has been already computed
//- Compute the normFactor used in convergence testing,
//- assumes ArowsSum has been already computed
void computeNormFactor(Mat A, Vec psi, Vec source)
{
if (!init_aux_vectors_) return;
......
src/petsc4Foam/utils/petscUtils.C
View file @ 8ebbdb63
......@@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2019 OpenCFD Ltd.
Copyright (C) 2019-2020 OpenCFD Ltd.
Copyright (C) 2019 Simone Bna
Copyright (C) 2020 Stefano Zampini
-------------------------------------------------------------------------------
......@@ -31,11 +31,10 @@ License
#include "error.H"
#include "petscUtils.H"
#include "petscLinearSolverContext.H"
#include <algorithm>
// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
Foam::scalar Foam::gAverage(Vec input)
Foam::solveScalar Foam::gAverage(Vec input)
{
PetscInt len;
VecGetSize(input, &len);
......@@ -51,11 +50,11 @@ Foam::scalar Foam::gAverage(Vec input)
WarningInFunction
<< "Empty PETSc Vec, returning zero" << endl;
return Zero;
return 0;
}
Foam::scalar Foam::gSum(Vec input)
Foam::solveScalar Foam::gSum(Vec input)
{
PetscScalar val;
VecSum(input, &val);
......@@ -64,63 +63,6 @@ Foam::scalar Foam::gSum(Vec input)
}
void Foam::PetscUtils::copyVector
(
Vec input,
List<scalar>& output,
const bool syncPar /* unused */
)
{
const PetscScalar *data;
VecGetArrayRead(input, &data);
PetscInt len;
VecGetLocalSize(input, &len);
output.resize(len);
if (sizeof(PetscScalar) == sizeof(scalar))
{
// Slight optimization
std::memcpy(output.data(), data, output.byteSize());
}
else
{
std::copy(data, (data + len), output.begin());
}
VecRestoreArrayRead(input, &data);
}
void Foam::PetscUtils::copyVector
(
const UList<scalar>& input,
Vec output,
const bool syncPar
)
{
if (syncPar)
{
VecSetSizes(output, input.size(), PETSC_DECIDE);
}
PetscScalar *data;
VecGetArrayWrite(output, &data);
if (sizeof(PetscScalar) == sizeof(scalar))
{
// Slight optimization
std::memcpy(data, input.cdata(), input.byteSize());
}
else
{
std::copy(input.cbegin(), input.cend(), data);
}
VecRestoreArrayWrite(output, &data);
}
Foam::scalar Foam::PetscUtils::normFactor
(
Vec AdotPsi,
......@@ -155,7 +97,7 @@ Foam::scalar Foam::PetscUtils::normFactor
const PetscScalar* sourceValues;
VecGetArrayRead(source, &sourceValues);
scalar normFactor = Zero;
scalar normFactor{0};
PetscInt len;
VecGetLocalSize(psi, &len);
......@@ -183,6 +125,7 @@ Foam::scalar Foam::PetscUtils::normFactor
);
}
PetscErrorCode Foam::PetscUtils::foamKSPMonitorFoam
(
KSP ksp,
......@@ -200,6 +143,7 @@ PetscErrorCode Foam::PetscUtils::foamKSPMonitorFoam
PetscFunctionBeginUser;
auto* ctx = static_cast<petscLinearSolverContext*>(cctx);
// compute L1 norm and rescale by normFactor
ierr = KSPBuildResidual(ksp, ctx->res_l1_w[0], ctx->res_l1_w[1], &ctx->res_l1_w[1]);CHKERRQ(ierr);
ierr = VecNorm(ctx->res_l1_w[1], NORM_1, &fnorm);CHKERRQ(ierr);
......@@ -238,6 +182,7 @@ PetscErrorCode Foam::PetscUtils::foamKSPMonitorFoam
PetscFunctionReturn(0);
}
PetscErrorCode Foam::PetscUtils::foamKSPMonitorRecordInit
(
KSP ksp,
......@@ -248,12 +193,15 @@ PetscErrorCode Foam::PetscUtils::foamKSPMonitorRecordInit
{
PetscFunctionBeginUser;
if (it) PetscFunctionReturn(0);
auto* ctx = static_cast<petscLinearSolverContext*>(cctx);
solverPerformance* solverPerf = ctx->solverPerf;
solverPerf->initialResidual() = rnorm;
solverPerformance& solverPerf = ctx->performance;
solverPerf.initialResidual() = rnorm;
PetscFunctionReturn(0);
}
PetscErrorCode Foam::PetscUtils::foamKSPConverge
(
KSP ksp,
......@@ -270,7 +218,7 @@ PetscErrorCode Foam::PetscUtils::foamKSPConverge
// Equivalent to the OpenFOAM checkConvergence function
//
auto* ctx = static_cast<petscLinearSolverContext*>(cctx);
solverPerformance* solverPerf = ctx->solverPerf;
solverPerformance& solverPerf = ctx->performance;
PetscReal rtol;
PetscReal abstol;
......@@ -292,33 +240,30 @@ PetscErrorCode Foam::PetscUtils::foamKSPConverge
if (it == 0)
{
solverPerf->initialResidual() = residual;
solverPerf.initialResidual() = residual;
}
if (residual < abstol)
{
*reason = KSP_CONVERGED_ATOL;
solverPerf->finalResidual() = residual;
}
else if
(
rtol > Foam::SMALL /** pTraits<Type>::one */
&& residual < rtol * solverPerf->initialResidual()
&& residual < rtol * solverPerf.initialResidual()
)
{
*reason = KSP_CONVERGED_RTOL;
solverPerf->finalResidual() = residual;
}
else if (it >= maxits)
{
*reason = KSP_DIVERGED_ITS;
solverPerf->finalResidual() = residual;
}
else if (residual > divtol)
{
*reason = KSP_DIVERGED_DTOL;
solverPerf->finalResidual() = residual;
}
solverPerf.finalResidual() = residual;
PetscFunctionReturn(0);
}
......
src/petsc4Foam/utils/petscUtils.H
View file @ 8ebbdb63
......@@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2019 OpenCFD Ltd.
Copyright (C) 2019-2020 OpenCFD Ltd.
Copyright (C) 2019 Simone Bna
-------------------------------------------------------------------------------
License
......@@ -51,10 +51,10 @@ namespace Foam
{
//- Global average value in PETSc vector
scalar gAverage(Vec input);
solveScalar gAverage(Vec input);
//- Global sum of PETSc vector
scalar gSum(Vec input);
solveScalar gSum(Vec input);
namespace PetscUtils
......@@ -64,22 +64,6 @@ namespace PetscUtils
Namespace PetscUtils
\*---------------------------------------------------------------------------*/
//- Copy a PETSc vector to an OpenFOAM List
void copyVector
(
Vec input,
List<scalar>& output,
const bool syncPar = false
);
//- Copy an OpenFOAM List to a PETSc vector
void copyVector
(
const UList<scalar>& input,
Vec output,
const bool syncPar = false
);
//- Calculate the OpenFOAM normFactor (akin to an L1 norm)
scalar normFactor
(
......
tutorials/basic/laplacianFoam/flange/Allrun.init
View file @ 8ebbdb63
......@@ -3,6 +3,7 @@ cd "${0%/*}" || exit # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions # Tutorial run functions
#------------------------------------------------------------------------------
tutorial="basic/laplacianFoam/flange"
solution="fvSolution-petsc"
src="$FOAM_TUTORIALS/$tutorial"
dst=run
......@@ -57,13 +58,15 @@ then
echo "libs (petscFoam);" >> "$dst"/system/controlDict
fi
# Use fvSolution-petsc
# Use PETSc fvSolution
if [ -f "$dst"/system/fvSolution ] \
&& [ ! -f "$dst"/system/fvSolution-petsc ]
&& [ -f "$solution" ] \
&& [ ! -f "$dst"/system/"$solution" ]
then
echo "Rename fvSolution and relink to fvSolution-petsc"
echo "Rename fvSolution and relink to $solution"
mv "$dst"/system/fvSolution "$dst"/system/fvSolution-foam
(cd "$dst"/system && ln -sf ../../fvSolution-petsc fvSolution)
cp -f "$solution" "$dst"/system/"$solution"
(cd "$dst"/system && ln -sf "$solution" fvSolution)
fi
#------------------------------------------------------------------------------
tutorials/basic/laplacianFoam/flange/fvSolution-petsc
View file @ 8ebbdb63
......@@ -19,8 +19,7 @@ solvers
{
T
{
solver petsc;
preconditioner petsc;
solver petsc;
petsc
{
......
tutorials/basic/laplacianFoam/pipeOneD/system/fvSolution-petsc
View file @ 8ebbdb63
......@@ -19,8 +19,7 @@ solvers
{
T
{
solver petsc;
preconditioner petsc;
solver petsc;
petsc
{
......
tutorials/incompressible/simpleFoam/motorBike/Allrun.init
View file @ 8ebbdb63
......@@ -3,6 +3,7 @@ cd "${0%/*}" || exit # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions # Tutorial run functions
#------------------------------------------------------------------------------
tutorial="incompressible/simpleFoam/motorBike"
solution="fvSolution-petsc"
src="$FOAM_TUTORIALS/$tutorial"
dst=run
......@@ -57,13 +58,15 @@ then
echo "libs (petscFoam);" >> "$dst"/system/controlDict
fi
# Use fvSolution-petsc
# Use PETSc fvSolution
if [ -f "$dst"/system/fvSolution ] \
&& [ ! -f "$dst"/system/fvSolution-petsc ]
&& [ -f "$solution" ] \
&& [ ! -f "$dst"/system/"$solution" ]
then
echo "Rename fvSolution and relink to fvSolution-petsc"
echo "Rename fvSolution and relink to $solution"
mv "$dst"/system/fvSolution "$dst"/system/fvSolution-foam
(cd "$dst"/system && ln -sf ../../fvSolution-petsc fvSolution)
cp -f "$solution" "$dst"/system/"$solution"
(cd "$dst"/system && ln -sf "$solution" fvSolution)
fi
#------------------------------------------------------------------------------
tutorials/incompressible/simpleFoam/motorBike/fvSolution-petsc
View file @ 8ebbdb63
......@@ -18,8 +18,7 @@ solvers
{
p
{
solver petsc;
preconditioner petsc;
solver petsc;
petsc
{
......@@ -75,8 +74,7 @@ solvers
U
{
solver petsc;
preconditioner petsc;
solver petsc;
petsc
{
......@@ -85,7 +83,6 @@ solvers
ksp_type bicg;
pc_type bjacobi;
sub_pc_type ilu;
ksp_cg_single_reduction true;
}
caching
......@@ -108,8 +105,7 @@ solvers
k
{
solver petsc;
preconditioner petsc;
solver petsc;
petsc
{
......@@ -118,7 +114,6 @@ solvers
ksp_type bicg;
pc_type bjacobi;
sub_pc_type ilu;
ksp_cg_single_reduction true;
}
caching
......@@ -141,8 +136,7 @@ solvers
omega
{
solver petsc;
preconditioner petsc;
solver petsc;
petsc
{
......@@ -151,7 +145,6 @@ solvers
ksp_type bicg;
pc_type bjacobi;
sub_pc_type ilu;
ksp_cg_single_reduction true;
}
caching
......
tutorials/incompressible/simpleFoam/pitzDaily/Allrun.init
View file @ 8ebbdb63
......@@ -3,6 +3,7 @@ cd "${0%/*}" || exit # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions # Tutorial run functions
#------------------------------------------------------------------------------
tutorial="incompressible/simpleFoam/pitzDaily"
solution="fvSolution-petsc"
src="$FOAM_TUTORIALS/$tutorial"
dst=run
......@@ -57,13 +58,15 @@ then
echo "libs (petscFoam);" >> "$dst"/system/controlDict
fi
# Use fvSolution-petsc
# Use PETSc fvSolution
if [ -f "$dst"/system/fvSolution ] \
&& [ ! -f "$dst"/system/fvSolution-petsc ]
&& [ -f "$solution" ] \
&& [ ! -f "$dst"/system/"$solution" ]
then
echo "Rename fvSolution and relink to fvSolution-petsc"
echo "Rename fvSolution and relink to $solution"
mv "$dst"/system/fvSolution "$dst"/system/fvSolution-foam
(cd "$dst"/system && ln -sf ../../fvSolution-petsc fvSolution)
cp -f "$solution" "$dst"/system/"$solution"
(cd "$dst"/system && ln -sf "$solution" fvSolution)
fi
#------------------------------------------------------------------------------
tutorials/incompressible/simpleFoam/pitzDaily/fvSolution-petsc
View file @ 8ebbdb63
......@@ -19,8 +19,7 @@ solvers
{
p
{
solver petsc;
preconditioner petsc;
solver petsc;
petsc
{
......@@ -71,8 +70,7 @@ solvers
"(U|k|epsilon|omega|f|v2)"
{
solver petsc;
preconditioner petsc;
solver petsc;
petsc
{
......@@ -81,7 +79,6 @@ solvers
ksp_type bicg;
pc_type bjacobi;
sub_pc_type ilu;
ksp_cg_single_reduction true;
}
caching
......