zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge branch 'feature-single-precision' into 'develop'

Browse Source 
Feature single precision

See merge request Development/OpenFOAM-plus!234
This commit is contained in:
Andrew Heather
2019-06-05 14:55:59 +01:00
committed by Andrew Heather
parent 24371a1d84 262faabd7c
commit 4c25755202
132 changed files with 1484 additions and 2810 deletions
Download Patch File Download Diff File Expand all files Collapse all files

14
src/OpenFOAM/dimensionedTypes/dimensionedType/dimensionedType.C
View File

@ -614,9 +614,12 @@ Foam::sqr(const dimensioned<Type>& dt)
}
template<class Type>
Foam::dimensioned<Foam::scalar> Foam::magSqr(const dimensioned<Type>& dt)
Foam::dimensioned<typename Foam::typeOfMag<Type>::type>
Foam::magSqr(const dimensioned<Type>& dt)
{
return dimensioned<scalar>
typedef typename typeOfMag<Type>::type magType;
return dimensioned<magType>
(
"magSqr(" + dt.name() + ')',
magSqr(dt.dimensions()),
@ -625,9 +628,12 @@ Foam::dimensioned<Foam::scalar> Foam::magSqr(const dimensioned<Type>& dt)
}
template<class Type>
Foam::dimensioned<Foam::scalar> Foam::mag(const dimensioned<Type>& dt)
Foam::dimensioned<typename Foam::typeOfMag<Type>::type>
Foam::mag(const dimensioned<Type>& dt)
{
return dimensioned<scalar>
typedef typename typeOfMag<Type>::type magType;
return dimensioned<magType>
(
"mag(" + dt.name() + ')',
dt.dimensions(),

6
src/OpenFOAM/dimensionedTypes/dimensionedType/dimensionedType.H
View File

@ -437,10 +437,12 @@ dimensioned<typename outerProduct<Type, Type>::type>
sqr(const dimensioned<Type>&);
template<class Type>
dimensioned<scalar> magSqr(const dimensioned<Type>&);
dimensioned<typename typeOfMag<Type>::type>
magSqr(const dimensioned<Type>& dt);
template<class Type>
dimensioned<scalar> mag(const dimensioned<Type>&);
dimensioned<typename typeOfMag<Type>::type>
mag(const dimensioned<Type>& dt);
template<class Type>
dimensioned<Type> cmptMultiply

43
src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldFunctions.C
View File

@ -141,13 +141,13 @@ sqr(const tmp<DimensionedField<Type, GeoMesh>>& tdf)
template<class Type, class GeoMesh>
tmp<DimensionedField<scalar, GeoMesh>> magSqr
(
const DimensionedField<Type, GeoMesh>& df
)
tmp<DimensionedField<typename typeOfMag<Type>::type, GeoMesh>>
magSqr(const DimensionedField<Type, GeoMesh>& df)
{
typedef typename typeOfMag<Type>::type magType;
auto tres =
tmp<DimensionedField<scalar, GeoMesh>>::New
tmp<DimensionedField<magType, GeoMesh>>::New
(
IOobject
(
@ -165,15 +165,15 @@ tmp<DimensionedField<scalar, GeoMesh>> magSqr
}
template<class Type, class GeoMesh>
tmp<DimensionedField<scalar, GeoMesh>> magSqr
(
const tmp<DimensionedField<Type, GeoMesh>>& tdf
)
tmp<DimensionedField<typename typeOfMag<Type>::type, GeoMesh>>
magSqr(const tmp<DimensionedField<Type, GeoMesh>>& tdf)
{
typedef typename typeOfMag<Type>::type magType;
const DimensionedField<Type, GeoMesh>& df = tdf();
auto tres =
reuseTmpDimensionedField<scalar, Type, GeoMesh>::New
reuseTmpDimensionedField<magType, Type, GeoMesh>::New
(
tdf,
"magSqr(" + df.name() + ')',
@ -188,13 +188,13 @@ tmp<DimensionedField<scalar, GeoMesh>> magSqr
template<class Type, class GeoMesh>
tmp<DimensionedField<scalar, GeoMesh>> mag
(
const DimensionedField<Type, GeoMesh>& df
)
tmp<DimensionedField<typename typeOfMag<Type>::type, GeoMesh>>
mag(const DimensionedField<Type, GeoMesh>& df)
{
typedef typename typeOfMag<Type>::type magType;
auto tres =
tmp<DimensionedField<scalar, GeoMesh>>::New
tmp<DimensionedField<magType, GeoMesh>>::New
(
IOobject
(
@ -212,15 +212,15 @@ tmp<DimensionedField<scalar, GeoMesh>> mag
}
template<class Type, class GeoMesh>
tmp<DimensionedField<scalar, GeoMesh>> mag
(
const tmp<DimensionedField<Type, GeoMesh>>& tdf
)
tmp<DimensionedField<typename typeOfMag<Type>::type, GeoMesh>>
mag(const tmp<DimensionedField<Type, GeoMesh>>& tdf)
{
typedef typename typeOfMag<Type>::type magType;
const DimensionedField<Type, GeoMesh>& df = tdf();
auto tres =
reuseTmpDimensionedField<scalar, Type, GeoMesh>::New
reuseTmpDimensionedField<magType, Type, GeoMesh>::New
(
tdf,
"mag(" + df.name() + ')',
@ -325,11 +325,10 @@ UNARY_REDUCTION_FUNCTION(Type, max, gMax)
UNARY_REDUCTION_FUNCTION(Type, min, gMin)
UNARY_REDUCTION_FUNCTION(Type, sum, gSum)
UNARY_REDUCTION_FUNCTION(Type, average, gAverage)
UNARY_REDUCTION_FUNCTION(MinMax<Type>, minMax, gMinMax)
UNARY_REDUCTION_FUNCTION(scalarMinMax, minMaxMag, gMinMaxMag)
UNARY_REDUCTION_FUNCTION(scalar, sumMag, gSumMag)
UNARY_REDUCTION_FUNCTION(typename typeOfMag<Type>::type, sumMag, gSumMag)
#undef UNARY_REDUCTION_FUNCTION

36
src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldFunctions.H
View File

@ -62,34 +62,29 @@ tmp<DimensionedField<typename outerProduct<Type, Type>::type, GeoMesh>>
sqr(const tmp<DimensionedField<Type, GeoMesh>>& tdf);
template<class Type, class GeoMesh>
tmp<DimensionedField<scalar, GeoMesh>> magSqr
(
const DimensionedField<Type, GeoMesh>& df
);
tmp<DimensionedField<typename typeOfMag<Type>::type, GeoMesh>>
magSqr(const DimensionedField<Type, GeoMesh>& df);
template<class Type, class GeoMesh>
tmp<DimensionedField<scalar, GeoMesh>> magSqr
(
const tmp<DimensionedField<Type, GeoMesh>>& tdf
);
tmp<DimensionedField<typename typeOfMag<Type>::type, GeoMesh>>
magSqr(const tmp<DimensionedField<Type, GeoMesh>>& tdf);
template<class Type, class GeoMesh>
tmp<DimensionedField<scalar, GeoMesh>> mag
(
const DimensionedField<Type, GeoMesh>& df
);
tmp<DimensionedField<typename typeOfMag<Type>::type, GeoMesh>>
mag(const DimensionedField<Type, GeoMesh>& df);
template<class Type, class GeoMesh>
tmp<DimensionedField<scalar, GeoMesh>> mag
(
const tmp<DimensionedField<Type, GeoMesh>>& tdf
);
tmp<DimensionedField<typename typeOfMag<Type>::type, GeoMesh>>
mag(const tmp<DimensionedField<Type, GeoMesh>>& tdf);
template<class Type, class GeoMesh>
tmp
<
DimensionedField
<typename DimensionedField<Type, GeoMesh>::cmptType, GeoMesh>
<
typename DimensionedField<Type, GeoMesh>::cmptType,
GeoMesh
>
>
cmptAv(const DimensionedField<Type, GeoMesh>& df);
@ -97,7 +92,10 @@ template<class Type, class GeoMesh>
tmp
<
DimensionedField
<typename DimensionedField<Type, GeoMesh>::cmptType, GeoMesh>
<
typename DimensionedField<Type, GeoMesh>::cmptType,
GeoMesh
>
>
cmptAv(const tmp<DimensionedField<Type, GeoMesh>>& tdf);
@ -123,7 +121,7 @@ UNARY_REDUCTION_FUNCTION(Type, average, gAverage)
UNARY_REDUCTION_FUNCTION(MinMax<Type>, minMax, gMinMax)
UNARY_REDUCTION_FUNCTION(scalarMinMax, minMaxMag, gMinMaxMag)
UNARY_REDUCTION_FUNCTION(scalar, sumMag, gSumMag)
UNARY_REDUCTION_FUNCTION(typename typeOfMag<Type>::type, sumMag, gSumMag)
#undef UNARY_REDUCTION_FUNCTION

65
src/OpenFOAM/fields/FieldFields/FieldField/FieldFieldFunctions.C
View File

@ -158,7 +158,11 @@ sqr(const tmp<FieldField<Field, Type>>& tf)
template<template<class> class Field, class Type>
void magSqr(FieldField<Field, scalar>& sf, const FieldField<Field, Type>& f)
void magSqr
(
FieldField<Field, typename typeOfMag<Type>::type>& sf,
const FieldField<Field, Type>& f
)
{
forAll(sf, i)
{
@ -167,11 +171,14 @@ void magSqr(FieldField<Field, scalar>& sf, const FieldField<Field, Type>& f)
}
template<template<class> class Field, class Type>
tmp<FieldField<Field, scalar>> magSqr(const FieldField<Field, Type>& f)
tmp<FieldField<Field, typename typeOfMag<Type>::type>>
magSqr(const FieldField<Field, Type>& f)
{
typedef typename typeOfMag<Type>::type magType;
auto tres
(
FieldField<Field, scalar>::NewCalculatedType(f)
FieldField<Field, magType>::NewCalculatedType(f)
);
magSqr(tres.ref(), f);
@ -179,11 +186,14 @@ tmp<FieldField<Field, scalar>> magSqr(const FieldField<Field, Type>& f)
}
template<template<class> class Field, class Type>
tmp<FieldField<Field, scalar>> magSqr(const tmp<FieldField<Field, Type>>& tf)
tmp<FieldField<Field, typename typeOfMag<Type>::type>>
magSqr(const tmp<FieldField<Field, Type>>& tf)
{
typedef typename typeOfMag<Type>::type magType;
auto tres
(
reuseTmpFieldField<Field, scalar, Type>::New(tf)
reuseTmpFieldField<Field, magType, Type>::New(tf)
);
magSqr(tres.ref(), tf());
@ -193,7 +203,11 @@ tmp<FieldField<Field, scalar>> magSqr(const tmp<FieldField<Field, Type>>& tf)
template<template<class> class Field, class Type>
void mag(FieldField<Field, scalar>& sf, const FieldField<Field, Type>& f)
void mag
(
FieldField<Field, typename typeOfMag<Type>::type>& sf,
const FieldField<Field, Type>& f
)
{
forAll(sf, i)
{
@ -202,11 +216,14 @@ void mag(FieldField<Field, scalar>& sf, const FieldField<Field, Type>& f)
}
template<template<class> class Field, class Type>
tmp<FieldField<Field, scalar>> mag(const FieldField<Field, Type>& f)
tmp<FieldField<Field, typename typeOfMag<Type>::type>>
mag(const FieldField<Field, Type>& f)
{
typedef typename typeOfMag<Type>::type magType;
auto tres
(
FieldField<Field, scalar>::NewCalculatedType(f)
FieldField<Field, magType>::NewCalculatedType(f)
);
mag(tres.ref(), f);
@ -214,11 +231,14 @@ tmp<FieldField<Field, scalar>> mag(const FieldField<Field, Type>& f)
}
template<template<class> class Field, class Type>
tmp<FieldField<Field, scalar>> mag(const tmp<FieldField<Field, Type>>& tf)
tmp<FieldField<Field, typename typeOfMag<Type>::type>>
mag(const tmp<FieldField<Field, Type>>& tf)
{
typedef typename typeOfMag<Type>::type magType;
auto tres
(
reuseTmpFieldField<Field, scalar, Type>::New(tf)
reuseTmpFieldField<Field, magType, Type>::New(tf)
);
mag(tres.ref(), tf());
@ -480,19 +500,21 @@ Type sum(const FieldField<Field, Type>& f)
TMP_UNARY_FUNCTION(Type, sum)
template<template<class> class Field, class Type>
scalar sumMag(const FieldField<Field, Type>& f)
typename typeOfMag<Type>::type sumMag(const FieldField<Field, Type>& f)
{
scalar SumMag = 0.0;
typedef typename typeOfMag<Type>::type magType;
magType result = Zero;
forAll(f, i)
{
SumMag += sumMag(f[i]);
result += sumMag(f[i]);
}
return SumMag;
return result;
}
TMP_UNARY_FUNCTION(scalar, sumMag)
TMP_UNARY_FUNCTION(typename typeOfMag<Type>::type, sumMag)
template<template<class> class Field, class Type>
Type average(const FieldField<Field, Type>& f)
@ -555,25 +577,24 @@ TMP_UNARY_FUNCTION(scalarMinMax, minMaxMag)
// With reduction on ReturnType
#define G_UNARY_FUNCTION(returnType, gFunc, func, rFunc) \
#define G_UNARY_FUNCTION(ReturnType, gFunc, func, rFunc) \
\
template<template<class> class Field, class Type> \
returnType gFunc(const FieldField<Field, Type>& f) \
ReturnType gFunc(const FieldField<Field, Type>& f) \
{ \
returnType res = func(f); \
reduce(res, rFunc##Op<returnType>()); \
ReturnType res = func(f); \
reduce(res, rFunc##Op<ReturnType>()); \
return res; \
} \
TMP_UNARY_FUNCTION(returnType, gFunc)
TMP_UNARY_FUNCTION(ReturnType, gFunc)
G_UNARY_FUNCTION(Type, gMax, max, max)
G_UNARY_FUNCTION(Type, gMin, min, min)
G_UNARY_FUNCTION(Type, gSum, sum, sum)
G_UNARY_FUNCTION(MinMax<Type>, gMinMax, minMax, sum)
G_UNARY_FUNCTION(scalarMinMax, gMinMaxMag, minMaxMag, sum)
G_UNARY_FUNCTION(scalar, gSumMag, sumMag, sum)
G_UNARY_FUNCTION(typename typeOfMag<Type>::type, gSumMag, sumMag, sum)
#undef G_UNARY_FUNCTION

37
src/OpenFOAM/fields/FieldFields/FieldField/FieldFieldFunctions.H
View File

@ -90,23 +90,35 @@ sqr(const tmp<FieldField<Field, Type>>& tf);
template<template<class> class Field, class Type>
void magSqr(FieldField<Field, scalar>& sf, const FieldField<Field, Type>& f);
void magSqr
(
FieldField<Field, typename typeOfMag<Type>::type>& sf,
const FieldField<Field, Type>& f
);
template<template<class> class Field, class Type>
tmp<FieldField<Field, scalar>> magSqr(const FieldField<Field, Type>& f);
tmp<FieldField<Field, typename typeOfMag<Type>::type>>
magSqr(const FieldField<Field, Type>& f);
template<template<class> class Field, class Type>
tmp<FieldField<Field, scalar>> magSqr(const tmp<FieldField<Field, Type>>& tf);
tmp<FieldField<Field, typename typeOfMag<Type>::type>>
magSqr(const tmp<FieldField<Field, Type>>& tf);
template<template<class> class Field, class Type>
void mag(FieldField<Field, scalar>& sf, const FieldField<Field, Type>& f);
void mag
(
FieldField<Field, typename typeOfMag<Type>::type>& res,
const FieldField<Field, Type>& f
);
template<template<class> class Field, class Type>
tmp<FieldField<Field, scalar>> mag(const FieldField<Field, Type>& f);
tmp<FieldField<Field, typename typeOfMag<Type>::type>>
mag(const FieldField<Field, Type>& f);
template<template<class> class Field, class Type>
tmp<FieldField<Field, scalar>> mag(const tmp<FieldField<Field, Type>>& tf);
tmp<FieldField<Field, typename typeOfMag<Type>::type>>
mag(const tmp<FieldField<Field, Type>>& tf);
template<template<class> class Field, class Type>
@ -214,9 +226,9 @@ TMP_UNARY_FUNCTION(Type, sum)
template<template<class> class Field, class Type>
scalar sumMag(const FieldField<Field, Type>& f);
typename typeOfMag<Type>::type sumMag(const FieldField<Field, Type>& f);
TMP_UNARY_FUNCTION(scalar, sumMag)
TMP_UNARY_FUNCTION(typename typeOfMag<Type>::type, sumMag)
template<template<class> class Field, class Type>
@ -239,20 +251,19 @@ TMP_UNARY_FUNCTION(scalarMinMax, minMaxMag)
// With reduction on ReturnType
#define G_UNARY_FUNCTION(returnType, gFunc, func, rFunc) \
#define G_UNARY_FUNCTION(ReturnType, gFunc, func, rFunc) \
\
template<template<class> class Field, class Type> \
returnType gFunc(const FieldField<Field, Type>& f); \
TMP_UNARY_FUNCTION(returnType, gFunc)
ReturnType gFunc(const FieldField<Field, Type>& f); \
TMP_UNARY_FUNCTION(ReturnType, gFunc)
G_UNARY_FUNCTION(Type, gMax, max, max)
G_UNARY_FUNCTION(Type, gMin, min, min)
G_UNARY_FUNCTION(Type, gSum, sum, sum)
G_UNARY_FUNCTION(MinMax<Type>, gMinMax, minMax, sum)
G_UNARY_FUNCTION(scalarMinMax, gMinMaxMag, minMaxMag, sum)
G_UNARY_FUNCTION(scalar, gSumMag, sumMag, sum)
G_UNARY_FUNCTION(typename typeOfMag<Type>::type, gSumMag, sumMag, sum)
#undef G_UNARY_FUNCTION

121
src/OpenFOAM/fields/Fields/Field/FieldFunctions.C
View File

@ -146,23 +146,35 @@ sqr(const tmp<Field<Type>>& tf)
template<class Type>
void magSqr(Field<scalar>& res, const UList<Type>& f)
void magSqr
(
Field<typename typeOfMag<Type>::type>& res,
const UList<Type>& f
)
{
TFOR_ALL_F_OP_FUNC_F(scalar, res, =, magSqr, Type, f)
typedef typename typeOfMag<Type>::type magType;
TFOR_ALL_F_OP_FUNC_F(magType, res, =, magSqr, Type, f)
}
template<class Type>
tmp<Field<scalar>> magSqr(const UList<Type>& f)
tmp<Field<typename typeOfMag<Type>::type>>
magSqr(const UList<Type>& f)
{
auto tres = tmp<Field<scalar>>::New(f.size());
typedef typename typeOfMag<Type>::type magType;
auto tres = tmp<Field<magType>>::New(f.size());
magSqr(tres.ref(), f);
return tres;
}
template<class Type>
tmp<Field<scalar>> magSqr(const tmp<Field<Type>>& tf)
tmp<Field<typename typeOfMag<Type>::type>>
magSqr(const tmp<Field<Type>>& tf)
{
auto tres = reuseTmp<scalar, Type>::New(tf);
typedef typename typeOfMag<Type>::type magType;
auto tres = reuseTmp<magType, Type>::New(tf);
magSqr(tres.ref(), tf());
tf.clear();
return tres;
@ -170,23 +182,35 @@ tmp<Field<scalar>> magSqr(const tmp<Field<Type>>& tf)
template<class Type>
void mag(Field<scalar>& res, const UList<Type>& f)
void mag
(
Field<typename typeOfMag<Type>::type>& res,
const UList<Type>& f
)
{
TFOR_ALL_F_OP_FUNC_F(scalar, res, =, mag, Type, f)
typedef typename typeOfMag<Type>::type magType;
TFOR_ALL_F_OP_FUNC_F(magType, res, =, mag, Type, f)
}
template<class Type>
tmp<Field<scalar>> mag(const UList<Type>& f)
tmp<Field<typename typeOfMag<Type>::type>>
mag(const UList<Type>& f)
{
auto tres = tmp<Field<scalar>>::New(f.size());
typedef typename typeOfMag<Type>::type magType;
auto tres = tmp<Field<magType>>::New(f.size());
mag(tres.ref(), f);
return tres;
}
template<class Type>
tmp<Field<scalar>> mag(const tmp<Field<Type>>& tf)
tmp<Field<typename typeOfMag<Type>::type>>
mag(const tmp<Field<Type>>& tf)
{
auto tres = reuseTmp<scalar, Type>::New(tf);
typedef typename typeOfMag<Type>::type magType;
auto tres = reuseTmp<magType, Type>::New(tf);
mag(tres.ref(), tf());
tf.clear();
return tres;
@ -341,14 +365,14 @@ TMP_UNARY_FUNCTION(Type, min)
template<class Type>
Type sum(const UList<Type>& f)
{
Type Sum = Zero;
if (f.size())
{
Type Sum = Zero;
TFOR_ALL_S_OP_F(Type, Sum, +=, Type, f)
return Sum;
}
return Zero;
return Sum;
}
TMP_UNARY_FUNCTION(Type, sum)
@ -413,15 +437,15 @@ Type minMagSqr(const UList<Type>& f)
TMP_UNARY_FUNCTION(Type, minMagSqr)
template<class Type>
typename pTraits<Type>::cmptType
typename scalarProduct<Type, Type>::type
sumProd(const UList<Type>& f1, const UList<Type>& f2)
{
typedef typename pTraits<Type>::cmptType outType;
typedef typename scalarProduct<Type, Type>::type prodType;
outType result = Zero;
prodType result = Zero;
if (f1.size() && (f1.size() == f2.size()))
{
TFOR_ALL_S_OP_F_OP_F(outType, result, +=, Type, f1, &&, Type, f2)
TFOR_ALL_S_OP_F_OP_F(prodType, result, +=, Type, f1, &&, Type, f2)
}
return result;
}
@ -450,41 +474,54 @@ Type sumCmptProd(const UList<Type>& f1, const UList<Type>& f2)
template<class Type>
scalar sumSqr(const UList<Type>& f)
typename outerProduct1<Type>::type
sumSqr(const UList<Type>& f)
{
scalar SumSqr = 0;
typedef typename outerProduct1<Type>::type prodType;
prodType result = Zero;
if (f.size())
{
TFOR_ALL_S_OP_FUNC_F(scalar, SumSqr, +=, sqr, Type, f)
TFOR_ALL_S_OP_FUNC_F(prodType, result, +=, sqr, Type, f)
}
return SumSqr;
return result;
}
TMP_UNARY_FUNCTION(scalar, sumSqr)
template<class Type>
scalar sumMag(const UList<Type>& f)
typename outerProduct1<Type>::type
sumSqr(const tmp<Field<Type>>& tf)
{
scalar SumMag = 0;
if (f.size())
{
TFOR_ALL_S_OP_FUNC_F(scalar, SumMag, +=, mag, Type, f)
}
return SumMag;
typedef typename outerProduct1<Type>::type prodType;
prodType result = sumSqr(tf());
tf.clear();
return result;
}
TMP_UNARY_FUNCTION(scalar, sumMag)
template<class Type>
typename typeOfMag<Type>::type
sumMag(const UList<Type>& f)
{
typedef typename typeOfMag<Type>::type magType;
magType result = Zero;
if (f.size())
{
TFOR_ALL_S_OP_FUNC_F(magType, result, +=, mag, Type, f)
}
return result;
}
TMP_UNARY_FUNCTION(typename typeOfMag<Type>::type, sumMag)
template<class Type>
Type sumCmptMag(const UList<Type>& f)
{
Type SumMag = Zero;
Type result = Zero;
if (f.size())
{
TFOR_ALL_S_OP_FUNC_F(scalar, SumMag, +=, cmptMag, Type, f)
TFOR_ALL_S_OP_FUNC_F(Type, result, +=, cmptMag, Type, f)
}
return SumMag;
return result;
}
TMP_UNARY_FUNCTION(Type, sumCmptMag)
@ -530,24 +567,24 @@ G_UNARY_FUNCTION(Type, gSumCmptMag, sumCmptMag, sum)
G_UNARY_FUNCTION(MinMax<Type>, gMinMax, minMax, sum)
G_UNARY_FUNCTION(scalarMinMax, gMinMaxMag, minMaxMag, sum)
G_UNARY_FUNCTION(scalar, gSumSqr, sumSqr, sum)
G_UNARY_FUNCTION(scalar, gSumMag, sumMag, sum)
G_UNARY_FUNCTION(typename outerProduct1<Type>::type, gSumSqr, sumSqr, sum)
G_UNARY_FUNCTION(typename typeOfMag<Type>::type, gSumMag, sumMag, sum)
#undef G_UNARY_FUNCTION
template<class Type>
typename pTraits<Type>::cmptType gSumProd
typename scalarProduct<Type, Type>::type gSumProd
(
const UList<Type>& f1,
const UList<Type>& f2,
const label comm
)
{
typedef typename pTraits<Type>::cmptType outType;
typedef typename scalarProduct<Type, Type>::type prodType;
outType result = sumProd(f1, f2);
reduce(result, sumOp<outType>(), Pstream::msgType(), comm);
prodType result = sumProd(f1, f2);
reduce(result, sumOp<prodType>(), Pstream::msgType(), comm);
return result;
}

45
src/OpenFOAM/fields/Fields/Field/FieldFunctions.H
View File

@ -94,23 +94,35 @@ sqr(const tmp<Field<Type>>& tf);
template<class Type>
void magSqr(Field<scalar>& res, const UList<Type>& f);
void magSqr
(
Field<typename typeOfMag<Type>::type>& res,
const UList<Type>& f
);
template<class Type>
tmp<Field<scalar>> magSqr(const UList<Type>& f);
tmp<Field<typename typeOfMag<Type>::type>>
magSqr(const UList<Type>& f);
template<class Type>
tmp<Field<scalar>> magSqr(const tmp<Field<Type>>& tf);
tmp<Field<typename typeOfMag<Type>::type>>
magSqr(const tmp<Field<Type>>& tf);
template<class Type>
void mag(Field<scalar>& res, const UList<Type>& f);
void mag
(
Field<typename typeOfMag<Type>::type>& res,
const UList<Type>& f
);
template<class Type>
tmp<Field<scalar>> mag(const UList<Type>& f);
tmp<Field<typename typeOfMag<Type>::type>>
mag(const UList<Type>& f);
template<class Type>
tmp<Field<scalar>> mag(const tmp<Field<Type>>& tf);
tmp<Field<typename typeOfMag<Type>::type>>
mag(const tmp<Field<Type>>& tf);
template<class Type>
@ -197,7 +209,8 @@ TMP_UNARY_FUNCTION(Type, minMagSqr)
template<class Type>
typename pTraits<Type>::cmptType sumProd
typename scalarProduct<Type, Type>::type
sumProd
(
const UList<Type>& f1,
const UList<Type>& f2
@ -207,14 +220,15 @@ template<class Type>
Type sumCmptProd(const UList<Type>& f1, const UList<Type>& f2);
template<class Type>
scalar sumSqr(const UList<Type>& f);
TMP_UNARY_FUNCTION(scalar, sumSqr)
typename outerProduct1<Type>::type sumSqr(const UList<Type>& f);
template<class Type>
scalar sumMag(const UList<Type>& f);
typename outerProduct1<Type>::type sumSqr(const tmp<Field<Type>>& tf);
TMP_UNARY_FUNCTION(scalar, sumMag)
template<class Type>
typename typeOfMag<Type>::type sumMag(const UList<Type>& f);
TMP_UNARY_FUNCTION(typename typeOfMag<Type>::type, sumMag)
template<class Type>
Type sumCmptMag(const UList<Type>& f);
@ -244,13 +258,14 @@ G_UNARY_FUNCTION(Type, gSumCmptMag, sumCmptMag, sum)
G_UNARY_FUNCTION(MinMax<Type>, gMinMax, minMax, sum)
G_UNARY_FUNCTION(scalarMinMax, gMinMaxMag, minMaxMag, sum)
G_UNARY_FUNCTION(scalar, gSumSqr, sumSqr, sum)
G_UNARY_FUNCTION(scalar, gSumMag, sumMag, sum)
G_UNARY_FUNCTION(typename outerProduct1<Type>::type, gSumSqr, sumSqr, sum)
G_UNARY_FUNCTION(typename typeOfMag<Type>::type, gSumMag, sumMag, sum)
#undef G_UNARY_FUNCTION
template<class Type>
typename pTraits<Type>::cmptType gSumProd
typename scalarProduct<Type, Type>::type gSumProd
(
const UList<Type>& f1,
const UList<Type>& f2,

66
src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.C
View File

@ -257,7 +257,7 @@ sqr(const tmp<GeometricField<Type, PatchField, GeoMesh>>& tgf)
template<class Type, template<class> class PatchField, class GeoMesh>
void magSqr
(
GeometricField<scalar, PatchField, GeoMesh>& gsf,
GeometricField<typename typeOfMag<Type>::type, PatchField, GeoMesh>& gsf,
const GeometricField<Type, PatchField, GeoMesh>& gf
)
{
@ -268,13 +268,16 @@ void magSqr
template<class Type, template<class> class PatchField, class GeoMesh>
tmp<GeometricField<scalar, PatchField, GeoMesh>> magSqr
tmp<GeometricField<typename typeOfMag<Type>::type, PatchField, GeoMesh>>
magSqr
(
const GeometricField<Type, PatchField, GeoMesh>& gf
)
{
typedef typename typeOfMag<Type>::type magType;
auto tres =
tmp<GeometricField<scalar, PatchField, GeoMesh>>::New
tmp<GeometricField<magType, PatchField, GeoMesh>>::New
(
IOobject
(
@ -294,30 +297,13 @@ tmp<GeometricField<scalar, PatchField, GeoMesh>> magSqr
}
template<class Type, template<class> class PatchField, class GeoMesh>
tmp<GeometricField<scalar, PatchField, GeoMesh>> magSqr
tmp<GeometricField<typename typeOfMag<Type>::type, PatchField, GeoMesh>>
magSqr
(
const tmp<GeometricField<Type, PatchField, GeoMesh>>& tgf
)
{
const GeometricField<Type, PatchField, GeoMesh>& gf = tgf();
auto tres =
tmp<GeometricField<scalar, PatchField, GeoMesh>>::New
(
IOobject
(
"magSqr(" + gf.name() + ')',
gf.instance(),
gf.db(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
gf.mesh(),
sqr(gf.dimensions())
);
magSqr(tres.ref(), gf);
auto tres = magSqr(tgf.cref());
tgf.clear();
return tres;
@ -327,7 +313,7 @@ tmp<GeometricField<scalar, PatchField, GeoMesh>> magSqr
template<class Type, template<class> class PatchField, class GeoMesh>
void mag
(
GeometricField<scalar, PatchField, GeoMesh>& gsf,
GeometricField<typename typeOfMag<Type>::type, PatchField, GeoMesh>& gsf,
const GeometricField<Type, PatchField, GeoMesh>& gf
)
{
@ -338,13 +324,16 @@ void mag
template<class Type, template<class> class PatchField, class GeoMesh>
tmp<GeometricField<scalar, PatchField, GeoMesh>> mag
tmp<GeometricField<typename typeOfMag<Type>::type, PatchField, GeoMesh>>
mag
(
const GeometricField<Type, PatchField, GeoMesh>& gf
)
{
typedef typename typeOfMag<Type>::type magType;
auto tres =
tmp<GeometricField<scalar, PatchField, GeoMesh>>::New
tmp<GeometricField<magType, PatchField, GeoMesh>>::New
(
IOobject
(
@ -364,30 +353,13 @@ tmp<GeometricField<scalar, PatchField, GeoMesh>> mag
}
template<class Type, template<class> class PatchField, class GeoMesh>
tmp<GeometricField<scalar, PatchField, GeoMesh>> mag
tmp<GeometricField<typename typeOfMag<Type>::type, PatchField, GeoMesh>>
mag
(
const tmp<GeometricField<Type, PatchField, GeoMesh>>& tgf
)
{
const GeometricField<Type, PatchField, GeoMesh>& gf = tgf();
auto tres =
tmp<GeometricField<scalar, PatchField, GeoMesh>>::New
(
IOobject
(
"mag(" + gf.name() + ')',
gf.instance(),
gf.db(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
gf.mesh(),
gf.dimensions()
);
mag(tres.ref(), gf);
auto tres = mag(tgf.cref());
tgf.clear();
return tres;
@ -559,7 +531,7 @@ dimensioned<returnType> func \
UNARY_REDUCTION_FUNCTION(Type, sum, gSum)
UNARY_REDUCTION_FUNCTION(Type, average, gAverage)
UNARY_REDUCTION_FUNCTION(scalar, sumMag, gSumMag)
UNARY_REDUCTION_FUNCTION(typename typeOfMag<Type>::type, sumMag, gSumMag)
#undef UNARY_REDUCTION_FUNCTION

18
src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldFunctions.H
View File

@ -142,18 +142,20 @@ sqr(const tmp<GeometricField<Type, PatchField, GeoMesh>>& tgf);
template<class Type, template<class> class PatchField, class GeoMesh>
void magSqr
(
GeometricField<scalar, PatchField, GeoMesh>& gsf,
GeometricField<typename typeOfMag<Type>::type, PatchField, GeoMesh>& gsf,
const GeometricField<Type, PatchField, GeoMesh>& gf
);
template<class Type, template<class> class PatchField, class GeoMesh>
tmp<GeometricField<scalar, PatchField, GeoMesh>> magSqr
tmp<GeometricField<typename typeOfMag<Type>::type, PatchField, GeoMesh>>
magSqr
(
const GeometricField<Type, PatchField, GeoMesh>& gf
);
template<class Type, template<class> class PatchField, class GeoMesh>
tmp<GeometricField<scalar, PatchField, GeoMesh>> magSqr
tmp<GeometricField<typename typeOfMag<Type>::type, PatchField, GeoMesh>>
magSqr
(
const tmp<GeometricField<Type, PatchField, GeoMesh>>& tgf
);
@ -161,18 +163,20 @@ tmp<GeometricField<scalar, PatchField, GeoMesh>> magSqr
template<class Type, template<class> class PatchField, class GeoMesh>
void mag
(
GeometricField<scalar, PatchField, GeoMesh>& gsf,
GeometricField<typename typeOfMag<Type>::type, PatchField, GeoMesh>& gsf,
const GeometricField<Type, PatchField, GeoMesh>& gf
);
template<class Type, template<class> class PatchField, class GeoMesh>
tmp<GeometricField<scalar, PatchField, GeoMesh>> mag
tmp<GeometricField<typename typeOfMag<Type>::type, PatchField, GeoMesh>>
mag
(
const GeometricField<Type, PatchField, GeoMesh>& gf
);
template<class Type, template<class> class PatchField, class GeoMesh>
tmp<GeometricField<scalar, PatchField, GeoMesh>> mag
tmp<GeometricField<typename typeOfMag<Type>::type, PatchField, GeoMesh>>
mag
(
const tmp<GeometricField<Type, PatchField, GeoMesh>>& tgf
);
@ -252,7 +256,7 @@ dimensioned<returnType> func \
UNARY_REDUCTION_FUNCTION(Type, sum, gSum)
UNARY_REDUCTION_FUNCTION(Type, average, gAverage)
UNARY_REDUCTION_FUNCTION(scalar, sumMag, gSumMag)
UNARY_REDUCTION_FUNCTION(typename typeOfMag<Type>::type, sumMag, gSumMag)
#undef UNARY_REDUCTION_FUNCTION

4
src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterfaceFields/cyclicLduInterfaceField/cyclicLduInterfaceField.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2012 OpenFOAM Foundation
@ -46,7 +46,7 @@ Foam::cyclicLduInterfaceField::~cyclicLduInterfaceField()
void Foam::cyclicLduInterfaceField::transformCoupleField
(
scalarField& f,
solveScalarField& f,
const direction cmpt
) const
{

4
src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterfaceFields/cyclicLduInterfaceField/cyclicLduInterfaceField.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2012 OpenFOAM Foundation
@ -93,7 +93,7 @@ public:
//- Transform given patch component field
void transformCoupleField
(
scalarField& f,
solveScalarField& f,
const direction cmpt
) const;
};

6
src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterfaceFields/lduInterfaceField/lduInterfaceField.C
View File

@ -35,10 +35,4 @@ defineTypeNameAndDebug(lduInterfaceField, 0);
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::lduInterfaceField::~lduInterfaceField()
{}
// ************************************************************************* //

20
src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterfaceFields/lduInterfaceField/lduInterfaceField.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2013 OpenFOAM Foundation
@ -93,7 +93,7 @@ public:
//- Destructor
virtual ~lduInterfaceField();
virtual ~lduInterfaceField() = default;
// Member Functions
@ -137,11 +137,11 @@ public:
// or subtract coupled contributions to matrix
virtual void initInterfaceMatrixUpdate
(
scalarField&,
solveScalarField& result,
const bool add,
const scalarField&,
const scalarField&,
const direction,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
) const
{}
@ -150,11 +150,11 @@ public:
// or subtract coupled contributions to matrix
virtual void updateInterfaceMatrix
(
scalarField&,
solveScalarField& result,
const bool add,
const scalarField&,
const scalarField&,
const direction,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
) const = 0;

4
src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterfaceFields/processorLduInterfaceField/processorLduInterfaceField.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2012 OpenFOAM Foundation
@ -46,7 +46,7 @@ Foam::processorLduInterfaceField::~processorLduInterfaceField()
void Foam::processorLduInterfaceField::transformCoupleField
(
scalarField& f,
solveScalarField& f,
const direction cmpt
) const
{

4
src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterfaceFields/processorLduInterfaceField/processorLduInterfaceField.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2014 OpenFOAM Foundation
@ -99,7 +99,7 @@ public:
//- Transform given patch component field
void transformCoupleField
(
scalarField& f,
solveScalarField& f,
const direction cmpt
) const;
};

10
src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2017 OpenFOAM Foundation
@ -29,7 +29,7 @@ License
#include "IOstreams.H"
#include "Switch.H"
#include "objectRegistry.H"
#include "IOField.H"
#include "scalarIOField.H"
#include "Time.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -319,7 +319,7 @@ const Foam::scalarField& Foam::lduMatrix::upper() const
void Foam::lduMatrix::setResidualField
(
const Field<scalar>& residual,
const scalarField& residual,
const word& fieldName,
const bool initial
) const
@ -339,8 +339,8 @@ void Foam::lduMatrix::setResidualField
lookupName = word("residual:" + fieldName);
}
IOField<scalar>* residualPtr =
lduMesh_.thisDb().getObjectPtr<IOField<scalar>>(lookupName);
scalarIOField* residualPtr =
lduMesh_.thisDb().getObjectPtr<scalarIOField>(lookupName);
if (residualPtr)
{

59
src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2016-2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2017 OpenFOAM Foundation
@ -262,12 +262,12 @@ public:
//- Return the matrix norm used to normalise the residual for the
//- stopping criterion
scalar normFactor
solveScalarField::cmptType normFactor
(
const scalarField& psi,
const scalarField& source,
const scalarField& Apsi,
scalarField& tmpField
const solveScalarField& psi,
const solveScalarField& source,
const solveScalarField& Apsi,
solveScalarField& tmpField
) const;
};
@ -404,11 +404,20 @@ public:
//- Smooth the solution for a given number of sweeps
virtual void smooth
(
scalarField& psi,
solveScalarField& psi,
const scalarField& source,
const direction cmpt,
const label nSweeps
) const = 0;
//- Smooth the solution for a given number of sweeps
virtual void scalarSmooth
(
solveScalarField& psi,
const solveScalarField& source,
const direction cmpt,
const label nSweeps
) const = 0;
};
@ -495,8 +504,8 @@ public:
//- Return wA the preconditioned form of residual rA
virtual void precondition
(
scalarField& wA,
const scalarField& rA,
solveScalarField& wA,
const solveScalarField& rA,
const direction cmpt=0
) const = 0;
@ -505,8 +514,8 @@ public:
// This is only required for preconditioning asymmetric matrices.
virtual void preconditionT
(
scalarField& wT,
const scalarField& rT,
solveScalarField& wT,
const solveScalarField& rT,
const direction cmpt=0
) const
{
@ -624,8 +633,8 @@ public:
//- Matrix multiplication with updated interfaces.
void Amul
(
scalarField&,
const tmp<scalarField>&,
solveScalarField&,
const tmp<solveScalarField>&,
const FieldField<Field, scalar>&,
const lduInterfaceFieldPtrsList&,
const direction cmpt
@ -634,8 +643,8 @@ public:
//- Matrix transpose multiplication with updated interfaces.
void Tmul
(
scalarField&,
const tmp<scalarField>&,
solveScalarField&,
const tmp<solveScalarField>&,
const FieldField<Field, scalar>&,
const lduInterfaceFieldPtrsList&,
const direction cmpt
@ -645,7 +654,7 @@ public:
//- Sum the coefficients on each row of the matrix
void sumA
(
scalarField&,
solveScalarField&,
const FieldField<Field, scalar>&,
const lduInterfaceFieldPtrsList&
) const;
@ -653,17 +662,17 @@ public:
void residual
(
scalarField& rA,
const scalarField& psi,
solveScalarField& rA,
const solveScalarField& psi,
const scalarField& source,
const FieldField<Field, scalar>& interfaceBouCoeffs,
const lduInterfaceFieldPtrsList& interfaces,
const direction cmpt
) const;
tmp<scalarField> residual
tmp<solveScalarField> residual
(
const scalarField& psi,
const solveScalarField& psi,
const scalarField& source,
const FieldField<Field, scalar>& interfaceBouCoeffs,
const lduInterfaceFieldPtrsList& interfaces,
@ -678,8 +687,8 @@ public:
const bool add,
const FieldField<Field, scalar>& interfaceCoeffs,
const lduInterfaceFieldPtrsList& interfaces,
const scalarField& psiif,
scalarField& result,
const solveScalarField& psiif,
solveScalarField& result,
const direction cmpt
) const;
@ -689,8 +698,8 @@ public:
const bool add,
const FieldField<Field, scalar>& interfaceCoeffs,
const lduInterfaceFieldPtrsList& interfaces,
const scalarField& psiif,
scalarField& result,
const solveScalarField& psiif,
solveScalarField& result,
const direction cmpt
) const;
@ -698,7 +707,7 @@ public:
//- if it exists
void setResidualField
(
const Field<scalar>& residual,
const scalarField& residual,
const word& fieldName,
const bool initial
) const;

40
src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixATmul.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2017 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2017-2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
@ -35,17 +35,17 @@ Description
void Foam::lduMatrix::Amul
(
scalarField& Apsi,
const tmp<scalarField>& tpsi,
solveScalarField& Apsi,
const tmp<solveScalarField>& tpsi,
const FieldField<Field, scalar>& interfaceBouCoeffs,
const lduInterfaceFieldPtrsList& interfaces,
const direction cmpt
) const
{
scalar* __restrict__ ApsiPtr = Apsi.begin();
solveScalar* __restrict__ ApsiPtr = Apsi.begin();
const scalarField& psi = tpsi();
const scalar* const __restrict__ psiPtr = psi.begin();
const solveScalarField& psi = tpsi();
const solveScalar* const __restrict__ psiPtr = psi.begin();
const scalar* const __restrict__ diagPtr = diag().begin();
@ -98,17 +98,17 @@ void Foam::lduMatrix::Amul
void Foam::lduMatrix::Tmul
(
scalarField& Tpsi,
const tmp<scalarField>& tpsi,
solveScalarField& Tpsi,
const tmp<solveScalarField>& tpsi,
const FieldField<Field, scalar>& interfaceIntCoeffs,
const lduInterfaceFieldPtrsList& interfaces,
const direction cmpt
) const
{
scalar* __restrict__ TpsiPtr = Tpsi.begin();
solveScalar* __restrict__ TpsiPtr = Tpsi.begin();
const scalarField& psi = tpsi();
const scalar* const __restrict__ psiPtr = psi.begin();
const solveScalarField& psi = tpsi();
const solveScalar* const __restrict__ psiPtr = psi.begin();
const scalar* const __restrict__ diagPtr = diag().begin();
@ -159,12 +159,12 @@ void Foam::lduMatrix::Tmul
void Foam::lduMatrix::sumA
(
scalarField& sumA,
solveScalarField& sumA,
const FieldField<Field, scalar>& interfaceBouCoeffs,
const lduInterfaceFieldPtrsList& interfaces
) const
{
scalar* __restrict__ sumAPtr = sumA.begin();
solveScalar* __restrict__ sumAPtr = sumA.begin();
const scalar* __restrict__ diagPtr = diag().begin();
@ -208,17 +208,17 @@ void Foam::lduMatrix::sumA
void Foam::lduMatrix::residual
(
scalarField& rA,
const scalarField& psi,
solveScalarField& rA,
const solveScalarField& psi,
const scalarField& source,
const FieldField<Field, scalar>& interfaceBouCoeffs,
const lduInterfaceFieldPtrsList& interfaces,
const direction cmpt
) const
{
scalar* __restrict__ rAPtr = rA.begin();
solveScalar* __restrict__ rAPtr = rA.begin();
const scalar* const __restrict__ psiPtr = psi.begin();
const solveScalar* const __restrict__ psiPtr = psi.begin();
const scalar* const __restrict__ diagPtr = diag().begin();
const scalar* const __restrict__ sourcePtr = source.begin();
@ -278,16 +278,16 @@ void Foam::lduMatrix::residual
}
Foam::tmp<Foam::scalarField> Foam::lduMatrix::residual
Foam::tmp<Foam::Field<Foam::solveScalar>> Foam::lduMatrix::residual
(
const scalarField& psi,
const solveScalarField& psi,
const scalarField& source,
const FieldField<Field, scalar>& interfaceBouCoeffs,
const lduInterfaceFieldPtrsList& interfaces,
const direction cmpt
) const
{
tmp<scalarField> trA(new scalarField(psi.size()));
tmp<solveScalarField> trA(new solveScalarField(psi.size()));
residual(trA.ref(), psi, source, interfaceBouCoeffs, interfaces, cmpt);
return trA;
}

12
src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixSolver.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2016-2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2017 OpenFOAM Foundation
@ -170,12 +170,12 @@ void Foam::lduMatrix::solver::read(const dictionary& solverControls)
}
Foam::scalar Foam::lduMatrix::solver::normFactor
Foam::solveScalarField::cmptType Foam::lduMatrix::solver::normFactor
(
const scalarField& psi,
const scalarField& source,
const scalarField& Apsi,
scalarField& tmpField
const solveScalarField& psi,
const solveScalarField& source,
const solveScalarField& Apsi,
solveScalarField& tmpField
) const
{
// --- Calculate A dot reference value of psi

10
src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixUpdateMatrixInterfaces.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2017 OpenFOAM Foundation
@ -34,8 +34,8 @@ void Foam::lduMatrix::initMatrixInterfaces
const bool add,
const FieldField<Field, scalar>& coupleCoeffs,
const lduInterfaceFieldPtrsList& interfaces,
const scalarField& psiif,
scalarField& result,
const solveScalarField& psiif,
solveScalarField& result,
const direction cmpt
) const
{
@ -103,8 +103,8 @@ void Foam::lduMatrix::updateMatrixInterfaces
const bool add,
const FieldField<Field, scalar>& coupleCoeffs,
const lduInterfaceFieldPtrsList& interfaces,
const scalarField& psiif,
scalarField& result,
const solveScalarField& psiif,
solveScalarField& result,
const direction cmpt
) const
{

28
src/OpenFOAM/matrices/lduMatrix/preconditioners/DICPreconditioner/DICPreconditioner.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2015 OpenFOAM Foundation
@ -26,6 +26,7 @@ License
\*---------------------------------------------------------------------------*/
#include "DICPreconditioner.H"
#include <algorithm>
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -48,8 +49,11 @@ Foam::DICPreconditioner::DICPreconditioner
)
:
lduMatrix::preconditioner(sol),
rD_(sol.matrix().diag())
rD_(sol.matrix().diag().size())
{
const scalarField& diag = sol.matrix().diag();
std::copy(diag.begin(), diag.end(), rD_.begin());
calcReciprocalD(rD_, sol.matrix());
}
@ -58,11 +62,11 @@ Foam::DICPreconditioner::DICPreconditioner
void Foam::DICPreconditioner::calcReciprocalD
(
scalarField& rD,
solveScalarField& rD,
const lduMatrix& matrix
)
{
scalar* __restrict__ rDPtr = rD.begin();
solveScalar* __restrict__ rDPtr = rD.begin();
const label* const __restrict__ uPtr = matrix.lduAddr().upperAddr().begin();
const label* const __restrict__ lPtr = matrix.lduAddr().lowerAddr().begin();
@ -88,14 +92,14 @@ void Foam::DICPreconditioner::calcReciprocalD
void Foam::DICPreconditioner::precondition
(
scalarField& wA,
const scalarField& rA,
solveScalarField& wA,
const solveScalarField& rA,
const direction
) const
{
scalar* __restrict__ wAPtr = wA.begin();
const scalar* __restrict__ rAPtr = rA.begin();
const scalar* __restrict__ rDPtr = rD_.begin();
solveScalar* __restrict__ wAPtr = wA.begin();
const solveScalar* __restrict__ rAPtr = rA.begin();
const solveScalar* __restrict__ rDPtr = rD_.begin();
const label* const __restrict__ uPtr =
solver_.matrix().lduAddr().upperAddr().begin();
@ -104,9 +108,9 @@ void Foam::DICPreconditioner::precondition
const scalar* const __restrict__ upperPtr =
solver_.matrix().upper().begin();
label nCells = wA.size();
label nFaces = solver_.matrix().upper().size();
label nFacesM1 = nFaces - 1;
const label nCells = wA.size();
const label nFaces = solver_.matrix().upper().size();
const label nFacesM1 = nFaces - 1;
for (label cell=0; cell<nCells; cell++)
{

10
src/OpenFOAM/matrices/lduMatrix/preconditioners/DICPreconditioner/DICPreconditioner.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011 OpenFOAM Foundation
@ -60,7 +60,7 @@ class DICPreconditioner
// Private data
//- The reciprocal preconditioned diagonal
scalarField rD_;
solveScalarField rD_;
public:
@ -87,13 +87,13 @@ public:
// Member Functions
//- Calculate the reciprocal of the preconditioned diagonal
static void calcReciprocalD(scalarField& rD, const lduMatrix& matrix);
static void calcReciprocalD(solveScalarField&, const lduMatrix&);
//- Return wA the preconditioned form of residual rA
virtual void precondition
(
scalarField& wA,
const scalarField& rA,
solveScalarField& wA,
const solveScalarField& rA,
const direction cmpt=0
) const;
};

55
src/OpenFOAM/matrices/lduMatrix/preconditioners/DILUPreconditioner/DILUPreconditioner.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2015 OpenFOAM Foundation
@ -26,6 +26,7 @@ License
\*---------------------------------------------------------------------------*/
#include "DILUPreconditioner.H"
#include <algorithm>
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -48,8 +49,11 @@ Foam::DILUPreconditioner::DILUPreconditioner
)
:
lduMatrix::preconditioner(sol),
rD_(sol.matrix().diag())
rD_(sol.matrix().diag().size())
{
const scalarField& diag = sol.matrix().diag();
std::copy(diag.begin(), diag.end(), rD_.begin());
calcReciprocalD(rD_, sol.matrix());
}
@ -58,11 +62,11 @@ Foam::DILUPreconditioner::DILUPreconditioner
void Foam::DILUPreconditioner::calcReciprocalD
(
scalarField& rD,
solveScalarField& rD,
const lduMatrix& matrix
)
{
scalar* __restrict__ rDPtr = rD.begin();
solveScalar* __restrict__ rDPtr = rD.begin();
const label* const __restrict__ uPtr = matrix.lduAddr().upperAddr().begin();
const label* const __restrict__ lPtr = matrix.lduAddr().lowerAddr().begin();
@ -78,7 +82,7 @@ void Foam::DILUPreconditioner::calcReciprocalD
// Calculate the reciprocal of the preconditioned diagonal
label nCells = rD.size();
const label nCells = rD.size();
for (label cell=0; cell<nCells; cell++)
{
@ -89,14 +93,14 @@ void Foam::DILUPreconditioner::calcReciprocalD
void Foam::DILUPreconditioner::precondition
(
scalarField& wA,
const scalarField& rA,
solveScalarField& wA,
const solveScalarField& rA,
const direction
) const
{
scalar* __restrict__ wAPtr = wA.begin();
const scalar* __restrict__ rAPtr = rA.begin();
const scalar* __restrict__ rDPtr = rD_.begin();
solveScalar* __restrict__ wAPtr = wA.begin();
const solveScalar* __restrict__ rAPtr = rA.begin();
const solveScalar* __restrict__ rDPtr = rD_.begin();
const label* const __restrict__ uPtr =
solver_.matrix().lduAddr().upperAddr().begin();
@ -110,21 +114,18 @@ void Foam::DILUPreconditioner::precondition
const scalar* const __restrict__ lowerPtr =
solver_.matrix().lower().begin();
label nCells = wA.size();
label nFaces = solver_.matrix().upper().size();
label nFacesM1 = nFaces - 1;
const label nCells = wA.size();
const label nFaces = solver_.matrix().upper().size();
const label nFacesM1 = nFaces - 1;
for (label cell=0; cell<nCells; cell++)
{
wAPtr[cell] = rDPtr[cell]*rAPtr[cell];
}
label sface;
for (label face=0; face<nFaces; face++)
{
sface = losortPtr[face];
const label sface = losortPtr[face];
wAPtr[uPtr[sface]] -=
rDPtr[uPtr[sface]]*lowerPtr[sface]*wAPtr[lPtr[sface]];
}
@ -139,14 +140,14 @@ void Foam::DILUPreconditioner::precondition
void Foam::DILUPreconditioner::preconditionT
(
scalarField& wT,
const scalarField& rT,
solveScalarField& wT,
const solveScalarField& rT,
const direction
) const
{
scalar* __restrict__ wTPtr = wT.begin();
const scalar* __restrict__ rTPtr = rT.begin();
const scalar* __restrict__ rDPtr = rD_.begin();
solveScalar* __restrict__ wTPtr = wT.begin();
const solveScalar* __restrict__ rTPtr = rT.begin();
const solveScalar* __restrict__ rDPtr = rD_.begin();
const label* const __restrict__ uPtr =
solver_.matrix().lduAddr().upperAddr().begin();
@ -160,9 +161,9 @@ void Foam::DILUPreconditioner::preconditionT
const scalar* const __restrict__ lowerPtr =
solver_.matrix().lower().begin();
label nCells = wT.size();
label nFaces = solver_.matrix().upper().size();
label nFacesM1 = nFaces - 1;
const label nCells = wT.size();
const label nFaces = solver_.matrix().upper().size();
const label nFacesM1 = nFaces - 1;
for (label cell=0; cell<nCells; cell++)
{
@ -176,11 +177,9 @@ void Foam::DILUPreconditioner::preconditionT
}
label sface;
for (label face=nFacesM1; face>=0; face--)
{
sface = losortPtr[face];
const label sface = losortPtr[face];
wTPtr[lPtr[sface]] -=
rDPtr[lPtr[sface]]*lowerPtr[sface]*wTPtr[uPtr[sface]];
}

14
src/OpenFOAM/matrices/lduMatrix/preconditioners/DILUPreconditioner/DILUPreconditioner.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011 OpenFOAM Foundation
@ -60,7 +60,7 @@ class DILUPreconditioner
// Private data
//- The reciprocal preconditioned diagonal
scalarField rD_;
solveScalarField rD_;
public:
@ -87,21 +87,21 @@ public:
// Member Functions
//- Calculate the reciprocal of the preconditioned diagonal
static void calcReciprocalD(scalarField& rD, const lduMatrix& matrix);
static void calcReciprocalD(solveScalarField&, const lduMatrix&);
//- Return wA the preconditioned form of residual rA
virtual void precondition
(
scalarField& wA,
const scalarField& rA,
solveScalarField& wA,
const solveScalarField& rA,
const direction cmpt=0
) const;
//- Return wT the transpose-matrix preconditioned form of residual rT.
virtual void preconditionT
(
scalarField& wT,
const scalarField& rT,
solveScalarField& wT,
const solveScalarField& rT,
const direction cmpt=0
) const;
};

38
src/OpenFOAM/matrices/lduMatrix/preconditioners/FDICPreconditioner/FDICPreconditioner.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2015 OpenFOAM Foundation
@ -26,6 +26,7 @@ License
\*---------------------------------------------------------------------------*/
#include "FDICPreconditioner.H"
#include <algorithm>
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -48,13 +49,16 @@ Foam::FDICPreconditioner::FDICPreconditioner
)
:
lduMatrix::preconditioner(sol),
rD_(sol.matrix().diag()),
rD_(sol.matrix().diag().size()),
rDuUpper_(sol.matrix().upper().size()),
rDlUpper_(sol.matrix().upper().size())
{
scalar* __restrict__ rDPtr = rD_.begin();
scalar* __restrict__ rDuUpperPtr = rDuUpper_.begin();
scalar* __restrict__ rDlUpperPtr = rDlUpper_.begin();
const scalarField& diag = sol.matrix().diag();
std::copy(diag.begin(), diag.end(), rD_.begin());
solveScalar* __restrict__ rDPtr = rD_.begin();
solveScalar* __restrict__ rDuUpperPtr = rDuUpper_.begin();
solveScalar* __restrict__ rDlUpperPtr = rDlUpper_.begin();
const label* const __restrict__ uPtr =
solver_.matrix().lduAddr().upperAddr().begin();
@ -63,8 +67,8 @@ Foam::FDICPreconditioner::FDICPreconditioner
const scalar* const __restrict__ upperPtr =
solver_.matrix().upper().begin();
label nCells = rD_.size();
label nFaces = solver_.matrix().upper().size();
const label nCells = rD_.size();
const label nFaces = solver_.matrix().upper().size();
for (label face=0; face<nFaces; face++)
{
@ -89,26 +93,26 @@ Foam::FDICPreconditioner::FDICPreconditioner
void Foam::FDICPreconditioner::precondition
(
scalarField& wA,
const scalarField& rA,
solveScalarField& wA,
const solveScalarField& rA,
const direction
) const
{
scalar* __restrict__ wAPtr = wA.begin();
const scalar* __restrict__ rAPtr = rA.begin();
const scalar* __restrict__ rDPtr = rD_.begin();
solveScalar* __restrict__ wAPtr = wA.begin();
const solveScalar* __restrict__ rAPtr = rA.begin();
const solveScalar* __restrict__ rDPtr = rD_.begin();
const label* const __restrict__ uPtr =
solver_.matrix().lduAddr().upperAddr().begin();
const label* const __restrict__ lPtr =
solver_.matrix().lduAddr().lowerAddr().begin();
const scalar* const __restrict__ rDuUpperPtr = rDuUpper_.begin();
const scalar* const __restrict__ rDlUpperPtr = rDlUpper_.begin();
const solveScalar* const __restrict__ rDuUpperPtr = rDuUpper_.begin();
const solveScalar* const __restrict__ rDlUpperPtr = rDlUpper_.begin();
label nCells = wA.size();
label nFaces = solver_.matrix().upper().size();
label nFacesM1 = nFaces - 1;
const label nCells = wA.size();
const label nFaces = solver_.matrix().upper().size();
const label nFacesM1 = nFaces - 1;
for (label cell=0; cell<nCells; cell++)
{

12
src/OpenFOAM/matrices/lduMatrix/preconditioners/FDICPreconditioner/FDICPreconditioner.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011 OpenFOAM Foundation
@ -62,9 +62,9 @@ class FDICPreconditioner
// Private data
//- The reciprocal preconditioned diagonal
scalarField rD_;
scalarField rDuUpper_;
scalarField rDlUpper_;
solveScalarField rD_;
solveScalarField rDuUpper_;
solveScalarField rDlUpper_;
// Private Member Functions
@ -102,8 +102,8 @@ public:
//- Return wA the preconditioned form of residual rA
virtual void precondition
(
scalarField& wA,
const scalarField& rA,
solveScalarField& wA,
const solveScalarField& rA,
const direction cmpt=0
) const;
};

30
src/OpenFOAM/matrices/lduMatrix/preconditioners/GAMGPreconditioner/GAMGPreconditioner.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2013 OpenFOAM Foundation
@ -26,6 +26,7 @@ License
\*---------------------------------------------------------------------------*/
#include "GAMGPreconditioner.H"
#include "PrecisionAdaptor.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -82,29 +83,29 @@ void Foam::GAMGPreconditioner::readControls()
void Foam::GAMGPreconditioner::precondition
(
scalarField& wA,
const scalarField& rA,
solveScalarField& wA,
const solveScalarField& rA_ss,
const direction cmpt
) const
{
wA = 0.0;
scalarField AwA(wA.size());
scalarField finestCorrection(wA.size());
scalarField finestResidual(rA);
solveScalarField AwA(wA.size());
solveScalarField finestCorrection(wA.size());
solveScalarField finestResidual(rA_ss);
// Create coarse grid correction fields
PtrList<scalarField> coarseCorrFields;
PtrList<solveScalarField> coarseCorrFields;
// Create coarse grid sources
PtrList<scalarField> coarseSources;
PtrList<solveScalarField> coarseSources;
// Create the smoothers for all levels
PtrList<lduMatrix::smoother> smoothers;
// Scratch fields if processor-agglomerated coarse level meshes
// are bigger than original. Usually not needed
scalarField ApsiScratch;
scalarField finestCorrectionScratch;
solveScalarField ApsiScratch;
solveScalarField finestCorrectionScratch;
// Initialise the above data structures
initVcycle
@ -116,8 +117,15 @@ void Foam::GAMGPreconditioner::precondition
finestCorrectionScratch
);
// Storage area when solveScalar != scalar
scalarField rA_s;
for (label cycle=0; cycle<nVcycles_; cycle++)
{
const scalarField& rA =
ConstPrecisionAdaptor<scalar, solveScalar>::get(rA_ss, rA_s);
Vcycle
(
smoothers,
@ -143,7 +151,7 @@ void Foam::GAMGPreconditioner::precondition
{
// Calculate finest level residual field
matrix_.Amul(AwA, wA, interfaceBouCoeffs_, interfaces_, cmpt);
finestResidual = rA;
finestResidual = rA_ss;
finestResidual -= AwA;
}
}

6
src/OpenFOAM/matrices/lduMatrix/preconditioners/GAMGPreconditioner/GAMGPreconditioner.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
@ -93,8 +93,8 @@ public:
//- Return wA the preconditioned form of residual rA
virtual void precondition
(
scalarField& wA,
const scalarField& rA,
solveScalarField& wA,
const solveScalarField& rA,
const direction cmpt=0
) const;
};

18
src/OpenFOAM/matrices/lduMatrix/preconditioners/diagonalPreconditioner/diagonalPreconditioner.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2004-2010-2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2015 OpenFOAM Foundation
@ -54,10 +54,10 @@ Foam::diagonalPreconditioner::diagonalPreconditioner
lduMatrix::preconditioner(sol),
rD(sol.matrix().diag().size())
{
scalar* __restrict__ rDPtr = rD.begin();
solveScalar* __restrict__ rDPtr = rD.begin();
const scalar* __restrict__ DPtr = solver_.matrix().diag().begin();
label nCells = rD.size();
const label nCells = rD.size();
// Generate reciprocal diagonal
for (label cell=0; cell<nCells; cell++)
@ -71,16 +71,16 @@ Foam::diagonalPreconditioner::diagonalPreconditioner
void Foam::diagonalPreconditioner::precondition
(
scalarField& wA,
const scalarField& rA,
solveScalarField& wA,
const solveScalarField& rA,
const direction
) const
{
scalar* __restrict__ wAPtr = wA.begin();
const scalar* __restrict__ rAPtr = rA.begin();
const scalar* __restrict__ rDPtr = rD.begin();
solveScalar* __restrict__ wAPtr = wA.begin();
const solveScalar* __restrict__ rAPtr = rA.begin();
const solveScalar* __restrict__ rDPtr = rD.begin();
label nCells = wA.size();
const label nCells = wA.size();
for (label cell=0; cell<nCells; cell++)
{

12
src/OpenFOAM/matrices/lduMatrix/preconditioners/diagonalPreconditioner/diagonalPreconditioner.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011 OpenFOAM Foundation
@ -61,7 +61,7 @@ class diagonalPreconditioner
// Private data
//- The reciprocal diagonal
scalarField rD;
solveScalarField rD;
// Private Member Functions
@ -99,16 +99,16 @@ public:
//- Return wA the preconditioned form of residual rA
virtual void precondition
(
scalarField& wA,
const scalarField& rA,
solveScalarField& wA,
const solveScalarField& rA,
const direction cmpt=0
) const;
//- Return wT the transpose-matrix preconditioned form of residual rT.
virtual void preconditionT
(
scalarField& wT,
const scalarField& rT,
solveScalarField& wT,
const solveScalarField& rT,
const direction cmpt=0
) const
{

12
src/OpenFOAM/matrices/lduMatrix/preconditioners/noPreconditioner/noPreconditioner.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2015 OpenFOAM Foundation
@ -59,15 +59,15 @@ Foam::noPreconditioner::noPreconditioner
void Foam::noPreconditioner::precondition
(
scalarField& wA,
const scalarField& rA,
solveScalarField& wA,
const solveScalarField& rA,
const direction
) const
{
scalar* __restrict__ wAPtr = wA.begin();
const scalar* __restrict__ rAPtr = rA.begin();
solveScalar* __restrict__ wAPtr = wA.begin();
const solveScalar* __restrict__ rAPtr = rA.begin();
label nCells = wA.size();
const label nCells = wA.size();
for (label cell=0; cell<nCells; cell++)
{

10
src/OpenFOAM/matrices/lduMatrix/preconditioners/noPreconditioner/noPreconditioner.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011 OpenFOAM Foundation
@ -90,16 +90,16 @@ public:
//- Return wA the preconditioned form of residual rA
virtual void precondition
(
scalarField& wA,
const scalarField& rA,
solveScalarField& wA,
const solveScalarField& rA,
const direction cmpt=0
) const;
//- Return wT the transpose-matrix preconditioned form of residual rT.
virtual void preconditionT
(
scalarField& wT,
const scalarField& rT,
solveScalarField& wT,
const solveScalarField& rT,
const direction cmpt=0
) const
{

48
src/OpenFOAM/matrices/lduMatrix/smoothers/DIC/DICSmoother.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2015 OpenFOAM Foundation
@ -27,6 +27,8 @@ License
#include "DICSmoother.H"
#include "DICPreconditioner.H"
#include "PrecisionAdaptor.H"
#include <algorithm>
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -58,8 +60,11 @@ Foam::DICSmoother::DICSmoother
interfaceIntCoeffs,
interfaces
),
rD_(matrix_.diag())
rD_(matrix_.diag().size())
{
const scalarField& diag = matrix_.diag();
std::copy(diag.begin(), diag.end(), rD_.begin());
DICPreconditioner::calcReciprocalD(rD_, matrix_);
}
@ -68,13 +73,13 @@ Foam::DICSmoother::DICSmoother
void Foam::DICSmoother::smooth
(
scalarField& psi,
solveScalarField& psi,
const scalarField& source,
const direction cmpt,
const label nSweeps
) const
{
const scalar* const __restrict__ rDPtr = rD_.begin();
const solveScalar* const __restrict__ rDPtr = rD_.begin();
const scalar* const __restrict__ upperPtr = matrix_.upper().begin();
const label* const __restrict__ uPtr =
matrix_.lduAddr().upperAddr().begin();
@ -82,8 +87,8 @@ void Foam::DICSmoother::smooth
matrix_.lduAddr().lowerAddr().begin();
// Temporary storage for the residual
scalarField rA(rD_.size());
scalar* __restrict__ rAPtr = rA.begin();
solveScalarField rA(rD_.size());
solveScalar* __restrict__ rAPtr = rA.begin();
for (label sweep=0; sweep<nSweeps; sweep++)
{
@ -97,19 +102,22 @@ void Foam::DICSmoother::smooth
cmpt
);
rA *= rD_;
forAll(rA, i)
{
rA[i] *= rD_[i];
}
label nFaces = matrix_.upper().size();
const label nFaces = matrix_.upper().size();
for (label facei=0; facei<nFaces; facei++)
{
label u = uPtr[facei];
const label u = uPtr[facei];
rAPtr[u] -= rDPtr[u]*upperPtr[facei]*rAPtr[lPtr[facei]];
}
label nFacesM1 = nFaces - 1;
const label nFacesM1 = nFaces - 1;
for (label facei=nFacesM1; facei>=0; facei--)
{
label l = lPtr[facei];
const label l = lPtr[facei];
rAPtr[l] -= rDPtr[l]*upperPtr[facei]*rAPtr[uPtr[facei]];
}
@ -118,4 +126,22 @@ void Foam::DICSmoother::smooth
}
void Foam::DICSmoother::scalarSmooth
(
solveScalarField& psi,
const solveScalarField& source,
const direction cmpt,
const label nSweeps
) const
{
smooth
(
psi,
ConstPrecisionAdaptor<scalar, solveScalar>(source),
cmpt,
nSweeps
);
}
// ************************************************************************* //

15
src/OpenFOAM/matrices/lduMatrix/smoothers/DIC/DICSmoother.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011 OpenFOAM Foundation
@ -62,7 +62,7 @@ class DICSmoother
// Private data
//- The reciprocal preconditioned diagonal
scalarField rD_;
solveScalarField rD_;
public:
@ -89,11 +89,20 @@ public:
//- Smooth the solution for a given number of sweeps
void smooth
(
scalarField& psi,
solveScalarField& psi,
const scalarField& source,
const direction cmpt,
const label nSweeps
) const;
//- Smooth the solution for a given number of sweeps
void scalarSmooth
(
solveScalarField& psi,
const solveScalarField& source,
const direction cmpt,
const label nSweeps
) const;
};

17
src/OpenFOAM/matrices/lduMatrix/smoothers/DICGaussSeidel/DICGaussSeidelSmoother.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011 OpenFOAM Foundation
@ -80,7 +80,7 @@ Foam::DICGaussSeidelSmoother::DICGaussSeidelSmoother
void Foam::DICGaussSeidelSmoother::smooth
(
scalarField& psi,
solveScalarField& psi,
const scalarField& source,
const direction cmpt,
const label nSweeps
@ -91,4 +91,17 @@ void Foam::DICGaussSeidelSmoother::smooth
}
void Foam::DICGaussSeidelSmoother::scalarSmooth
(
solveScalarField& psi,
const solveScalarField& source,
const direction cmpt,
const label nSweeps
) const
{
dicSmoother_.scalarSmooth(psi, source, cmpt, nSweeps);
gsSmoother_.scalarSmooth(psi, source, cmpt, nSweeps);
}
// ************************************************************************* //

13
src/OpenFOAM/matrices/lduMatrix/smoothers/DICGaussSeidel/DICGaussSeidelSmoother.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011 OpenFOAM Foundation
@ -89,11 +89,20 @@ public:
//- Smooth the solution for a given number of sweeps
virtual void smooth
(
scalarField& psi,
solveScalarField& psi,
const scalarField& Source,
const direction cmpt,
const label nSweeps
) const;
//- Smooth the solution for a given number of sweeps
virtual void scalarSmooth
(
solveScalarField& psi,
const solveScalarField& Source,
const direction cmpt,
const label nSweeps
) const;
};

48
src/OpenFOAM/matrices/lduMatrix/smoothers/DILU/DILUSmoother.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2015 OpenFOAM Foundation
@ -27,6 +27,8 @@ License
#include "DILUSmoother.H"
#include "DILUPreconditioner.H"
#include "PrecisionAdaptor.H"
#include <algorithm>
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -58,8 +60,11 @@ Foam::DILUSmoother::DILUSmoother
interfaceIntCoeffs,
interfaces
),
rD_(matrix_.diag())
rD_(matrix_.diag().size())
{
const scalarField& diag = matrix_.diag();
std::copy(diag.begin(), diag.end(), rD_.begin());
DILUPreconditioner::calcReciprocalD(rD_, matrix_);
}
@ -68,13 +73,13 @@ Foam::DILUSmoother::DILUSmoother
void Foam::DILUSmoother::smooth
(
scalarField& psi,
solveScalarField& psi,
const scalarField& source,
const direction cmpt,
const label nSweeps
) const
{
const scalar* const __restrict__ rDPtr = rD_.begin();
const solveScalar* const __restrict__ rDPtr = rD_.begin();
const label* const __restrict__ uPtr =
matrix_.lduAddr().upperAddr().begin();
@ -85,8 +90,8 @@ void Foam::DILUSmoother::smooth
const scalar* const __restrict__ lowerPtr = matrix_.lower().begin();
// Temporary storage for the residual
scalarField rA(rD_.size());
scalar* __restrict__ rAPtr = rA.begin();
solveScalarField rA(rD_.size());
solveScalar* __restrict__ rAPtr = rA.begin();
for (label sweep=0; sweep<nSweeps; sweep++)
{
@ -100,19 +105,22 @@ void Foam::DILUSmoother::smooth
cmpt
);
rA *= rD_;
forAll(rA, i)
{
rA[i] *= rD_[i];
}
label nFaces = matrix_.upper().size();
const label nFaces = matrix_.upper().size();
for (label face=0; face<nFaces; face++)
{
label u = uPtr[face];
const label u = uPtr[face];
rAPtr[u] -= rDPtr[u]*lowerPtr[face]*rAPtr[lPtr[face]];
}
label nFacesM1 = nFaces - 1;
const label nFacesM1 = nFaces - 1;
for (label face=nFacesM1; face>=0; face--)
{
label l = lPtr[face];
const label l = lPtr[face];
rAPtr[l] -= rDPtr[l]*upperPtr[face]*rAPtr[uPtr[face]];
}
@ -121,4 +129,22 @@ void Foam::DILUSmoother::smooth
}
void Foam::DILUSmoother::scalarSmooth
(
solveScalarField& psi,
const solveScalarField& source,
const direction cmpt,
const label nSweeps
) const
{
smooth
(
psi,
ConstPrecisionAdaptor<scalar, solveScalar>(source),
cmpt,
nSweeps
);
}
// ************************************************************************* //

15
src/OpenFOAM/matrices/lduMatrix/smoothers/DILU/DILUSmoother.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2013 OpenFOAM Foundation
@ -58,7 +58,7 @@ class DILUSmoother
// Private data
//- The reciprocal preconditioned diagonal
scalarField rD_;
solveScalarField rD_;
public:
@ -85,11 +85,20 @@ public:
//- Smooth the solution for a given number of sweeps
void smooth
(
scalarField& psi,
solveScalarField& psi,
const scalarField& source,
const direction cmpt,
const label nSweeps
) const;
//- Smooth the solution for a given number of sweeps
void scalarSmooth
(
solveScalarField& psi,
const solveScalarField& source,
const direction cmpt,
const label nSweeps
) const;
};

27
src/OpenFOAM/matrices/lduMatrix/smoothers/DILUGaussSeidel/DILUGaussSeidelSmoother.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011 OpenFOAM Foundation
@ -26,6 +26,7 @@ License
\*---------------------------------------------------------------------------*/
#include "DILUGaussSeidelSmoother.H"
#include "PrecisionAdaptor.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -79,16 +80,34 @@ Foam::DILUGaussSeidelSmoother::DILUGaussSeidelSmoother
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::DILUGaussSeidelSmoother::scalarSmooth
(
solveScalarField& psi,
const solveScalarField& source,
const direction cmpt,
const label nSweeps
) const
{
diluSmoother_.scalarSmooth(psi, source, cmpt, nSweeps);
gsSmoother_.scalarSmooth(psi, source, cmpt, nSweeps);
}
void Foam::DILUGaussSeidelSmoother::smooth
(
scalarField& psi,
solveScalarField& psi,
const scalarField& source,
const direction cmpt,
const label nSweeps
) const
{
diluSmoother_.smooth(psi, source, cmpt, nSweeps);
gsSmoother_.smooth(psi, source, cmpt, nSweeps);
scalarSmooth
(
psi,
ConstPrecisionAdaptor<solveScalar, scalar>(source),
cmpt,
nSweeps
);
}

13
src/OpenFOAM/matrices/lduMatrix/smoothers/DILUGaussSeidel/DILUGaussSeidelSmoother.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011 OpenFOAM Foundation
@ -86,11 +86,20 @@ public:
//- Smooth the solution for a given number of sweeps
virtual void smooth
(
scalarField& psi,
solveScalarField& psi,
const scalarField& Source,
const direction cmpt,
const label nSweeps
) const;
//- Smooth the solution for a given number of sweeps
void scalarSmooth
(
solveScalarField& psi,
const solveScalarField& source,
const direction cmpt,
const label nSweeps
) const;
};

35
src/OpenFOAM/matrices/lduMatrix/smoothers/FDIC/FDICSmoother.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2015 OpenFOAM Foundation
@ -27,6 +27,7 @@ License
#include "FDICSmoother.H"
#include "FDICPreconditioner.H"
#include "PrecisionAdaptor.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -73,8 +74,8 @@ Foam::FDICSmoother::FDICSmoother
const scalar* const __restrict__ upperPtr =
matrix_.upper().begin();
label nCells = rD_.size();
label nFaces = matrix_.upper().size();
const label nCells = rD_.size();
const label nFaces = matrix_.upper().size();
for (label face=0; face<nFaces; face++)
{
@ -99,7 +100,7 @@ Foam::FDICSmoother::FDICSmoother
void Foam::FDICSmoother::smooth
(
scalarField& psi,
solveScalarField& psi,
const scalarField& source,
const direction cmpt,
const label nSweeps
@ -114,8 +115,8 @@ void Foam::FDICSmoother::smooth
matrix_.lduAddr().lowerAddr().begin();
// Temporary storage for the residual
scalarField rA(rD_.size());
scalar* __restrict__ rAPtr = rA.begin();
solveScalarField rA(rD_.size());
solveScalar* __restrict__ rAPtr = rA.begin();
for (label sweep=0; sweep<nSweeps; sweep++)
{
@ -131,13 +132,13 @@ void Foam::FDICSmoother::smooth
rA *= rD_;
label nFaces = matrix_.upper().size();
const label nFaces = matrix_.upper().size();
for (label face=0; face<nFaces; face++)
{
rAPtr[uPtr[face]] -= rDuUpperPtr[face]*rAPtr[lPtr[face]];
}
label nFacesM1 = nFaces - 1;
const label nFacesM1 = nFaces - 1;
for (label face=nFacesM1; face>=0; face--)
{
rAPtr[lPtr[face]] -= rDlUpperPtr[face]*rAPtr[uPtr[face]];
@ -148,4 +149,22 @@ void Foam::FDICSmoother::smooth
}
void Foam::FDICSmoother::scalarSmooth
(
solveScalarField& psi,
const solveScalarField& source,
const direction cmpt,
const label nSweeps
) const
{
smooth
(
psi,
ConstPrecisionAdaptor<scalar, solveScalar>(source),
cmpt,
nSweeps
);
}
// ************************************************************************* //

13
src/OpenFOAM/matrices/lduMatrix/smoothers/FDIC/FDICSmoother.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2012 OpenFOAM Foundation
@ -88,11 +88,20 @@ public:
//- Smooth the solution for a given number of sweeps
void smooth
(
scalarField& psi,
solveScalarField& psi,
const scalarField& source,
const direction cmpt,
const label nSweeps
) const;
//- Smooth the solution for a given number of sweeps
virtual void scalarSmooth
(
solveScalarField& psi,
const solveScalarField& source,
const direction cmpt,
const label nSweeps
) const;
};

39
src/OpenFOAM/matrices/lduMatrix/smoothers/GaussSeidel/GaussSeidelSmoother.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2017 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2017-2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2015 OpenFOAM Foundation
@ -26,6 +26,7 @@ License
\*---------------------------------------------------------------------------*/
#include "GaussSeidelSmoother.H"
#include "PrecisionAdaptor.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -68,21 +69,21 @@ Foam::GaussSeidelSmoother::GaussSeidelSmoother
void Foam::GaussSeidelSmoother::smooth
(
const word& fieldName_,
scalarField& psi,
solveScalarField& psi,
const lduMatrix& matrix_,
const scalarField& source,
const solveScalarField& source,
const FieldField<Field, scalar>& interfaceBouCoeffs_,
const lduInterfaceFieldPtrsList& interfaces_,
const direction cmpt,
const label nSweeps
)
{
scalar* __restrict__ psiPtr = psi.begin();
solveScalar* __restrict__ psiPtr = psi.begin();
const label nCells = psi.size();
scalarField bPrime(nCells);
scalar* __restrict__ bPrimePtr = bPrime.begin();
solveScalarField bPrime(nCells);
solveScalar* __restrict__ bPrimePtr = bPrime.begin();
const scalar* const __restrict__ diagPtr = matrix_.diag().begin();
const scalar* const __restrict__ upperPtr =
@ -133,7 +134,7 @@ void Foam::GaussSeidelSmoother::smooth
cmpt
);
scalar psii;
solveScalar psii;
label fStart;
label fEnd = ownStartPtr[0];
@ -169,11 +170,33 @@ void Foam::GaussSeidelSmoother::smooth
void Foam::GaussSeidelSmoother::smooth
(
scalarField& psi,
solveScalarField& psi,
const scalarField& source,
const direction cmpt,
const label nSweeps
) const
{
smooth
(
fieldName_,
psi,
matrix_,
ConstPrecisionAdaptor<solveScalar, scalar>(source),
interfaceBouCoeffs_,
interfaces_,
cmpt,
nSweeps
);
}
void Foam::GaussSeidelSmoother::scalarSmooth
(
solveScalarField& psi,
const solveScalarField& source,
const direction cmpt,
const label nSweeps
) const
{
smooth
(

17
src/OpenFOAM/matrices/lduMatrix/smoothers/GaussSeidel/GaussSeidelSmoother.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011 OpenFOAM Foundation
@ -81,9 +81,9 @@ public:
static void smooth
(
const word& fieldName,
scalarField& psi,
solveScalarField& psi,
const lduMatrix& matrix,
const scalarField& source,
const solveScalarField& source,
const FieldField<Field, scalar>& interfaceBouCoeffs,
const lduInterfaceFieldPtrsList& interfaces,
const direction cmpt,
@ -94,11 +94,20 @@ public:
//- Smooth the solution for a given number of sweeps
virtual void smooth
(
scalarField& psi,
solveScalarField& psi,
const scalarField& Source,
const direction cmpt,
const label nSweeps
) const;
//- Smooth the solution for a given number of sweeps
virtual void scalarSmooth
(
solveScalarField& psi,
const solveScalarField& Source,
const direction cmpt,
const label nSweeps
) const;
};

40
src/OpenFOAM/matrices/lduMatrix/smoothers/nonBlockingGaussSeidel/nonBlockingGaussSeidelSmoother.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2017 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2017-2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
@ -26,6 +26,7 @@ License
\*---------------------------------------------------------------------------*/
#include "nonBlockingGaussSeidelSmoother.H"
#include "PrecisionAdaptor.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -95,22 +96,22 @@ Foam::nonBlockingGaussSeidelSmoother::nonBlockingGaussSeidelSmoother
void Foam::nonBlockingGaussSeidelSmoother::smooth
(
const word& fieldName_,
scalarField& psi,
solveScalarField& psi,
const lduMatrix& matrix_,
const label blockStart,
const scalarField& source,
const solveScalarField& source,
const FieldField<Field, scalar>& interfaceBouCoeffs_,
const lduInterfaceFieldPtrsList& interfaces_,
const direction cmpt,
const label nSweeps
)
{
scalar* __restrict__ psiPtr = psi.begin();
solveScalar* __restrict__ psiPtr = psi.begin();
const label nCells = psi.size();
scalarField bPrime(nCells);
scalar* __restrict__ bPrimePtr = bPrime.begin();
solveScalarField bPrime(nCells);
solveScalar* __restrict__ bPrimePtr = bPrime.begin();
const scalar* const __restrict__ diagPtr = matrix_.diag().begin();
const scalar* const __restrict__ upperPtr =
@ -136,7 +137,6 @@ void Foam::nonBlockingGaussSeidelSmoother::smooth
// To compensate for this, it is necessary to turn the
// sign of the contribution.
for (label sweep=0; sweep<nSweeps; sweep++)
{
bPrime = source;
@ -151,7 +151,7 @@ void Foam::nonBlockingGaussSeidelSmoother::smooth
cmpt
);
scalar curPsi;
solveScalar curPsi;
label fStart;
label fEnd = ownStartPtr[0];
@ -223,10 +223,10 @@ void Foam::nonBlockingGaussSeidelSmoother::smooth
}
void Foam::nonBlockingGaussSeidelSmoother::smooth
void Foam::nonBlockingGaussSeidelSmoother::scalarSmooth
(
scalarField& psi,
const scalarField& source,
solveScalarField& psi,
const solveScalarField& source,
const direction cmpt,
const label nSweeps
) const
@ -246,4 +246,22 @@ void Foam::nonBlockingGaussSeidelSmoother::smooth
}
void Foam::nonBlockingGaussSeidelSmoother::smooth
(
solveScalarField& psi,
const scalarField& source,
const direction cmpt,
const label nSweeps
) const
{
scalarSmooth
(
psi,
ConstPrecisionAdaptor<solveScalar, scalar>(source),
cmpt,
nSweeps
);
}
// ************************************************************************* //

17
src/OpenFOAM/matrices/lduMatrix/smoothers/nonBlockingGaussSeidel/nonBlockingGaussSeidelSmoother.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011 OpenFOAM Foundation
@ -93,10 +93,10 @@ public:
static void smooth
(
const word& fieldName,
scalarField& psi,
solveScalarField& psi,
const lduMatrix& matrix,
const label blockStart,
const scalarField& source,
const solveScalarField& source,
const FieldField<Field, scalar>& interfaceBouCoeffs,
const lduInterfaceFieldPtrsList& interfaces,
const direction cmpt,
@ -106,11 +106,20 @@ public:
//- Smooth the solution for a given number of sweeps
virtual void smooth
(
scalarField& psi,
solveScalarField& psi,
const scalarField& Source,
const direction cmpt,
const label nSweeps
) const;
//- Smooth the solution for a given number of sweeps
virtual void scalarSmooth
(
solveScalarField& psi,
const solveScalarField& source,
const direction cmpt,
const label nSweeps
) const;
};

39
src/OpenFOAM/matrices/lduMatrix/smoothers/symGaussSeidel/symGaussSeidelSmoother.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2017 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2017-2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2012-2015 OpenFOAM Foundation
@ -26,6 +26,7 @@ License
\*---------------------------------------------------------------------------*/
#include "symGaussSeidelSmoother.H"
#include "PrecisionAdaptor.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -68,21 +69,21 @@ Foam::symGaussSeidelSmoother::symGaussSeidelSmoother
void Foam::symGaussSeidelSmoother::smooth
(
const word& fieldName_,
scalarField& psi,
solveScalarField& psi,
const lduMatrix& matrix_,
const scalarField& source,
const solveScalarField& source,
const FieldField<Field, scalar>& interfaceBouCoeffs_,
const lduInterfaceFieldPtrsList& interfaces_,
const direction cmpt,
const label nSweeps
)
{
scalar* __restrict__ psiPtr = psi.begin();
solveScalar* __restrict__ psiPtr = psi.begin();
const label nCells = psi.size();
scalarField bPrime(nCells);
scalar* __restrict__ bPrimePtr = bPrime.begin();
solveScalarField bPrime(nCells);
solveScalar* __restrict__ bPrimePtr = bPrime.begin();
const scalar* const __restrict__ diagPtr = matrix_.diag().begin();
const scalar* const __restrict__ upperPtr =
@ -133,7 +134,7 @@ void Foam::symGaussSeidelSmoother::smooth
cmpt
);
scalar psii;
solveScalar psii;
label fStart;
label fEnd = ownStartPtr[0];
@ -196,10 +197,10 @@ void Foam::symGaussSeidelSmoother::smooth
}
void Foam::symGaussSeidelSmoother::smooth
void Foam::symGaussSeidelSmoother::scalarSmooth
(
scalarField& psi,
const scalarField& source,
solveScalarField& psi,
const solveScalarField& source,
const direction cmpt,
const label nSweeps
) const
@ -218,4 +219,22 @@ void Foam::symGaussSeidelSmoother::smooth
}
void Foam::symGaussSeidelSmoother::smooth
(
solveScalarField& psi,
const scalarField& source,
const direction cmpt,
const label nSweeps
) const
{
scalarSmooth
(
psi,
ConstPrecisionAdaptor<solveScalar, scalar>(source),
cmpt,
nSweeps
);
}
// ************************************************************************* //

17
src/OpenFOAM/matrices/lduMatrix/smoothers/symGaussSeidel/symGaussSeidelSmoother.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2012 OpenFOAM Foundation
@ -81,9 +81,9 @@ public:
static void smooth
(
const word& fieldName,
scalarField& psi,
solveScalarField& psi,
const lduMatrix& matrix,
const scalarField& source,
const solveScalarField& source,
const FieldField<Field, scalar>& interfaceBouCoeffs,
const lduInterfaceFieldPtrsList& interfaces,
const direction cmpt,
@ -94,11 +94,20 @@ public:
//- Smooth the solution for a given number of sweeps
virtual void smooth
(
scalarField& psi,
solveScalarField& psi,
const scalarField& Source,
const direction cmpt,
const label nSweeps
) const;
//- Smooth the solution for a given number of sweeps
virtual void scalarSmooth
(
solveScalarField& psi,
const solveScalarField& source,
const direction cmpt,
const label nSweeps
) const;
};

46
src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
@ -230,8 +230,8 @@ class GAMGSolver
//- Interpolate the correction after injected prolongation
void interpolate
(
scalarField& psi,
scalarField& Apsi,
solveScalarField& psi,
solveScalarField& Apsi,
const lduMatrix& m,
const FieldField<Field, scalar>& interfaceBouCoeffs,
const lduInterfaceFieldPtrsList& interfaces,
@ -242,13 +242,13 @@ class GAMGSolver
// re-normalise
void interpolate
(
scalarField& psi,
scalarField& Apsi,
solveScalarField& psi,
solveScalarField& Apsi,
const lduMatrix& m,
const FieldField<Field, scalar>& interfaceBouCoeffs,
const lduInterfaceFieldPtrsList& interfaces,
const labelList& restrictAddressing,
const scalarField& psiC,
const solveScalarField& psiC,
const direction cmpt
) const;
@ -259,23 +259,23 @@ class GAMGSolver
// scaling factor.
void scale
(
scalarField& field,
scalarField& Acf,
solveScalarField& field,
solveScalarField& Acf,
const lduMatrix& A,
const FieldField<Field, scalar>& interfaceLevelBouCoeffs,
const lduInterfaceFieldPtrsList& interfaceLevel,
const scalarField& source,
const solveScalarField& source,
const direction cmpt
) const;
//- Initialise the data structures for the V-cycle
void initVcycle
(
PtrList<scalarField>& coarseCorrFields,
PtrList<scalarField>& coarseSources,
PtrList<solveScalarField>& coarseCorrFields,
PtrList<solveScalarField>& coarseSources,
PtrList<lduMatrix::smoother>& smoothers,
scalarField& scratch1,
scalarField& scratch2
solveScalarField& scratch1,
solveScalarField& scratch2
) const;
@ -283,17 +283,17 @@ class GAMGSolver
void Vcycle
(
const PtrList<lduMatrix::smoother>& smoothers,
scalarField& psi,
solveScalarField& psi,
const scalarField& source,
scalarField& Apsi,
scalarField& finestCorrection,
scalarField& finestResidual,
solveScalarField& Apsi,
solveScalarField& finestCorrection,
solveScalarField& finestResidual,
scalarField& scratch1,
scalarField& scratch2,
solveScalarField& scratch1,
solveScalarField& scratch2,
PtrList<scalarField>& coarseCorrFields,
PtrList<scalarField>& coarseSources,
PtrList<solveScalarField>& coarseCorrFields,
PtrList<solveScalarField>& coarseSources,
const direction cmpt=0
) const;
@ -316,8 +316,8 @@ class GAMGSolver
//- Solve the coarsest level with either an iterative or direct solver
void solveCoarsestLevel
(
scalarField& coarsestCorrField,
const scalarField& coarsestSource
solveScalarField& coarsestCorrField,
const solveScalarField& coarsestSource
) const;

28
src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverInterpolate.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2017 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2017-2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2013-2015 OpenFOAM Foundation
@ -31,15 +31,15 @@ License
void Foam::GAMGSolver::interpolate
(
scalarField& psi,
scalarField& Apsi,
solveScalarField& psi,
solveScalarField& Apsi,
const lduMatrix& m,
const FieldField<Field, scalar>& interfaceBouCoeffs,
const lduInterfaceFieldPtrsList& interfaces,
const direction cmpt
) const
{
scalar* __restrict__ psiPtr = psi.begin();
solveScalar* __restrict__ psiPtr = psi.begin();
const label* const __restrict__ uPtr = m.lduAddr().upperAddr().begin();
const label* const __restrict__ lPtr = m.lduAddr().lowerAddr().begin();
@ -49,7 +49,7 @@ void Foam::GAMGSolver::interpolate
const scalar* const __restrict__ lowerPtr = m.lower().begin();
Apsi = 0;
scalar* __restrict__ ApsiPtr = Apsi.begin();
solveScalar* __restrict__ ApsiPtr = Apsi.begin();
m.initMatrixInterfaces
(
@ -88,13 +88,13 @@ void Foam::GAMGSolver::interpolate
void Foam::GAMGSolver::interpolate
(
scalarField& psi,
scalarField& Apsi,
solveScalarField& psi,
solveScalarField& Apsi,
const lduMatrix& m,
const FieldField<Field, scalar>& interfaceBouCoeffs,
const lduInterfaceFieldPtrsList& interfaces,
const labelList& restrictAddressing,
const scalarField& psiC,
const solveScalarField& psiC,
const direction cmpt
) const
{
@ -109,17 +109,17 @@ void Foam::GAMGSolver::interpolate
);
const label nCells = m.diag().size();
scalar* __restrict__ psiPtr = psi.begin();
solveScalar* __restrict__ psiPtr = psi.begin();
const scalar* const __restrict__ diagPtr = m.diag().begin();
const scalar* const __restrict__ psiCPtr = psiC.begin();
const solveScalar* const __restrict__ psiCPtr = psiC.begin();
const label nCCells = psiC.size();
scalarField corrC(nCCells, 0);
scalar* __restrict__ corrCPtr = corrC.begin();
solveScalarField corrC(nCCells, 0);
solveScalar* __restrict__ corrCPtr = corrC.begin();
scalarField diagC(nCCells, 0);
scalar* __restrict__ diagCPtr = diagC.begin();
solveScalarField diagC(nCCells, 0);
solveScalar* __restrict__ diagCPtr = diagC.begin();
for (label celli=0; celli<nCells; celli++)
{

26
src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverScale.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2017 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2017-2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2017 OpenFOAM Foundation
@ -32,12 +32,12 @@ License
void Foam::GAMGSolver::scale
(
scalarField& field,
scalarField& Acf,
solveScalarField& field,
solveScalarField& Acf,
const lduMatrix& A,
const FieldField<Field, scalar>& interfaceLevelBouCoeffs,
const lduInterfaceFieldPtrsList& interfaceLevel,
const scalarField& source,
const solveScalarField& source,
const direction cmpt
) const
{
@ -52,13 +52,13 @@ void Foam::GAMGSolver::scale
const label nCells = field.size();
scalar* __restrict__ fieldPtr = field.begin();
const scalar* const __restrict__ sourcePtr = source.begin();
const scalar* const __restrict__ AcfPtr = Acf.begin();
solveScalar* __restrict__ fieldPtr = field.begin();
const solveScalar* const __restrict__ sourcePtr = source.begin();
const solveScalar* const __restrict__ AcfPtr = Acf.begin();
scalar scalingFactorNum = 0.0;
scalar scalingFactorDenom = 0.0;
solveScalar scalingFactorNum = 0.0;
solveScalar scalingFactorDenom = 0.0;
for (label i=0; i<nCells; i++)
{
@ -66,10 +66,12 @@ void Foam::GAMGSolver::scale
scalingFactorDenom += AcfPtr[i]*fieldPtr[i];
}
vector2D scalingVector(scalingFactorNum, scalingFactorDenom);
A.mesh().reduce(scalingVector, sumOp<vector2D>());
Vector2D<solveScalar> scalingVector(scalingFactorNum, scalingFactorDenom);
A.mesh().reduce(scalingVector, sumOp<Vector2D<solveScalar>>());
const scalar sf = scalingVector.x()/stabilise(scalingVector.y(), VSMALL);
const solveScalar sf =
scalingVector.x()
/stabilise(scalingVector.y(), pTraits<solveScalar>::vsmall);
if (debug >= 2)
{

123
src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2016-2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2017 OpenFOAM Foundation
@ -29,29 +29,36 @@ License
#include "PCG.H"
#include "PBiCGStab.H"
#include "SubField.H"
#include "PrecisionAdaptor.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::solverPerformance Foam::GAMGSolver::solve
(
scalarField& psi,
scalarField& psi_s,
const scalarField& source,
const direction cmpt
) const
{
PrecisionAdaptor<solveScalar, scalar> tpsi(psi_s);
solveScalarField& psi = tpsi.constCast();
ConstPrecisionAdaptor<solveScalar, scalar> tsource(source);
// Setup class containing solver performance data
solverPerformance solverPerf(typeName, fieldName_);
// Calculate A.psi used to calculate the initial residual
scalarField Apsi(psi.size());
solveScalarField Apsi(psi.size());
matrix_.Amul(Apsi, psi, interfaceBouCoeffs_, interfaces_, cmpt);
// Create the storage for the finestCorrection which may be used as a
// temporary in normFactor
scalarField finestCorrection(psi.size());
solveScalarField finestCorrection(psi.size());
// Calculate normalisation factor
scalar normFactor = this->normFactor(psi, source, Apsi, finestCorrection);
solveScalar normFactor =
this->normFactor(psi, tsource(), Apsi, finestCorrection);
if (debug >= 2)
{
@ -59,9 +66,14 @@ Foam::solverPerformance Foam::GAMGSolver::solve
}
// Calculate initial finest-grid residual field
scalarField finestResidual(source - Apsi);
solveScalarField finestResidual(tsource() - Apsi);
matrix().setResidualField(finestResidual, fieldName_, true);
matrix().setResidualField
(
ConstPrecisionAdaptor<scalar, solveScalar>(finestResidual)(),
fieldName_,
true
);
// Calculate normalised residual for convergence test
solverPerf.initialResidual() = gSumMag
@ -80,18 +92,18 @@ Foam::solverPerformance Foam::GAMGSolver::solve
)
{
// Create coarse grid correction fields
PtrList<scalarField> coarseCorrFields;
PtrList<solveScalarField> coarseCorrFields;
// Create coarse grid sources
PtrList<scalarField> coarseSources;
PtrList<solveScalarField> coarseSources;
// Create the smoothers for all levels
PtrList<lduMatrix::smoother> smoothers;
// Scratch fields if processor-agglomerated coarse level meshes
// are bigger than original. Usually not needed
scalarField scratch1;
scalarField scratch2;
solveScalarField scratch1;
solveScalarField scratch2;
// Initialise the above data structures
initVcycle
@ -124,7 +136,7 @@ Foam::solverPerformance Foam::GAMGSolver::solve
// Calculate finest level residual field
matrix_.Amul(Apsi, psi, interfaceBouCoeffs_, interfaces_, cmpt);
finestResidual = source;
finestResidual = tsource();
finestResidual -= Apsi;
solverPerf.finalResidual() = gSumMag
@ -147,7 +159,12 @@ Foam::solverPerformance Foam::GAMGSolver::solve
);
}
matrix().setResidualField(finestResidual, fieldName_, false);
matrix().setResidualField
(
ConstPrecisionAdaptor<scalar, solveScalar>(finestResidual)(),
fieldName_,
false
);
return solverPerf;
}
@ -156,17 +173,17 @@ Foam::solverPerformance Foam::GAMGSolver::solve
void Foam::GAMGSolver::Vcycle
(
const PtrList<lduMatrix::smoother>& smoothers,
scalarField& psi,
solveScalarField& psi,
const scalarField& source,
scalarField& Apsi,
scalarField& finestCorrection,
scalarField& finestResidual,
solveScalarField& Apsi,
solveScalarField& finestCorrection,
solveScalarField& finestResidual,
scalarField& scratch1,
scalarField& scratch2,
solveScalarField& scratch1,
solveScalarField& scratch2,
PtrList<scalarField>& coarseCorrFields,
PtrList<scalarField>& coarseSources,
PtrList<solveScalarField>& coarseCorrFields,
PtrList<solveScalarField>& coarseSources,
const direction cmpt
) const
{
@ -194,10 +211,10 @@ void Foam::GAMGSolver::Vcycle
{
coarseCorrFields[leveli] = 0.0;
smoothers[leveli + 1].smooth
smoothers[leveli + 1].scalarSmooth
(
coarseCorrFields[leveli],
coarseSources[leveli],
coarseSources[leveli], //coarseSource,
cmpt,
min
(
@ -206,7 +223,7 @@ void Foam::GAMGSolver::Vcycle
)
);
scalarField::subField ACf
solveScalarField::subField ACf
(
scratch1,
coarseCorrFields[leveli].size()
@ -219,9 +236,9 @@ void Foam::GAMGSolver::Vcycle
scale
(
coarseCorrFields[leveli],
const_cast<scalarField&>
const_cast<solveScalarField&>
(
ACf.operator const scalarField&()
ACf.operator const solveScalarField&()
),
matrixLevels_[leveli],
interfaceLevelsBouCoeffs_[leveli],
@ -234,9 +251,9 @@ void Foam::GAMGSolver::Vcycle
// Correct the residual with the new solution
matrixLevels_[leveli].Amul
(
const_cast<scalarField&>
const_cast<solveScalarField&>
(
ACf.operator const scalarField&()
ACf.operator const solveScalarField&()
),
coarseCorrFields[leveli],
interfaceLevelsBouCoeffs_[leveli],
@ -282,7 +299,7 @@ void Foam::GAMGSolver::Vcycle
// Smoothing and prolongation of the coarse correction fields
// (going to finer levels)
scalarField dummyField(0);
solveScalarField dummyField(0);
for (label leveli = coarsestLevel - 1; leveli >= 0; leveli--)
{
@ -291,7 +308,7 @@ void Foam::GAMGSolver::Vcycle
// Create a field for the pre-smoothed correction field
// as a sub-field of the finestCorrection which is not
// currently being used
scalarField::subField preSmoothedCoarseCorrField
solveScalarField::subField preSmoothedCoarseCorrField
(
scratch2,
coarseCorrFields[leveli].size()
@ -318,13 +335,16 @@ void Foam::GAMGSolver::Vcycle
// Create A.psi for this coarse level as a sub-field of Apsi
scalarField::subField ACf
solveScalarField::subField ACf
(
scratch1,
coarseCorrFields[leveli].size()
);
scalarField& ACfRef =
const_cast<scalarField&>(ACf.operator const scalarField&());
solveScalarField& ACfRef =
const_cast
<
solveScalarField&
>(ACf.operator const solveScalarField&());
if (interpolateCorrection_) //&& leveli < coarsestLevel - 2)
{
@ -383,10 +403,10 @@ void Foam::GAMGSolver::Vcycle
coarseCorrFields[leveli] += preSmoothedCoarseCorrField;
}
smoothers[leveli + 1].smooth
smoothers[leveli + 1].scalarSmooth
(
coarseCorrFields[leveli],
coarseSources[leveli],
coarseSources[leveli], //coarseSource,
cmpt,
min
(
@ -453,11 +473,11 @@ void Foam::GAMGSolver::Vcycle
void Foam::GAMGSolver::initVcycle
(
PtrList<scalarField>& coarseCorrFields,
PtrList<scalarField>& coarseSources,
PtrList<solveScalarField>& coarseCorrFields,
PtrList<solveScalarField>& coarseSources,
PtrList<lduMatrix::smoother>& smoothers,
scalarField& scratch1,
scalarField& scratch2
solveScalarField& scratch1,
solveScalarField& scratch2
) const
{
label maxSize = matrix_.diag().size();
@ -487,7 +507,7 @@ void Foam::GAMGSolver::initVcycle
{
label nCoarseCells = agglomeration_.nCells(leveli);
coarseSources.set(leveli, new scalarField(nCoarseCells));
coarseSources.set(leveli, new solveScalarField(nCoarseCells));
}
if (matrixLevels_.set(leveli))
@ -498,7 +518,7 @@ void Foam::GAMGSolver::initVcycle
maxSize = max(maxSize, nCoarseCells);
coarseCorrFields.set(leveli, new scalarField(nCoarseCells));
coarseCorrFields.set(leveli, new solveScalarField(nCoarseCells));
smoothers.set
(
@ -555,8 +575,8 @@ Foam::dictionary Foam::GAMGSolver::PBiCGStabSolverDict
void Foam::GAMGSolver::solveCoarsestLevel
(
scalarField& coarsestCorrField,
const scalarField& coarsestSource
solveScalarField& coarsestCorrField,
const solveScalarField& coarsestSource
) const
{
const label coarsestLevel = matrixLevels_.size() - 1;
@ -565,7 +585,18 @@ void Foam::GAMGSolver::solveCoarsestLevel
if (directSolveCoarsest_)
{
coarsestLUMatrixPtr_->solve(coarsestCorrField, coarsestSource);
PrecisionAdaptor<scalar, solveScalar> tcorrField
(
coarsestCorrField
);
coarsestLUMatrixPtr_->solve
(
tcorrField.constCast(),
ConstPrecisionAdaptor<scalar, solveScalar>
(
coarsestSource
)()
);
}
//else if
//(
@ -681,7 +712,7 @@ void Foam::GAMGSolver::solveCoarsestLevel
interfaceLevelsIntCoeffs_[coarsestLevel],
interfaceLevels_[coarsestLevel],
PBiCGStabSolverDict(tolerance_, relTol_)
).solve
).scalarSolve
(
coarsestCorrField,
coarsestSource
@ -697,7 +728,7 @@ void Foam::GAMGSolver::solveCoarsestLevel
interfaceLevelsIntCoeffs_[coarsestLevel],
interfaceLevels_[coarsestLevel],
PCGsolverDict(tolerance_, relTol_)
).solve
).scalarSolve
(
coarsestCorrField,
coarsestSource

8
src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/cyclicGAMGInterfaceField/cyclicGAMGInterfaceField.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2013 OpenFOAM Foundation
@ -110,16 +110,16 @@ Foam::cyclicGAMGInterfaceField::~cyclicGAMGInterfaceField()
void Foam::cyclicGAMGInterfaceField::updateInterfaceMatrix
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes
) const
{
// Get neighbouring field
scalarField pnf
solveScalarField pnf
(
cyclicInterface_.neighbPatch().interfaceInternalField(psiInternal)
);

6
src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/cyclicGAMGInterfaceField/cyclicGAMGInterfaceField.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2013 OpenFOAM Foundation
@ -147,9 +147,9 @@ public:
//- Update result field based on interface functionality
virtual void updateInterfaceMatrix
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType

18
src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/processorGAMGInterfaceField/processorGAMGInterfaceField.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2017 OpenFOAM Foundation
@ -94,9 +94,9 @@ Foam::processorGAMGInterfaceField::~processorGAMGInterfaceField()
void Foam::processorGAMGInterfaceField::initInterfaceMatrixUpdate
(
scalarField&,
solveScalarField&,
const bool,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField&,
const direction,
const Pstream::commsTypes commsType
@ -145,9 +145,9 @@ void Foam::processorGAMGInterfaceField::initInterfaceMatrixUpdate
void Foam::processorGAMGInterfaceField::updateInterfaceMatrix
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField&,
const solveScalarField&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
@ -187,9 +187,13 @@ void Foam::processorGAMGInterfaceField::updateInterfaceMatrix
}
else
{
scalarField pnf
solveScalarField pnf
(
procInterface_.compressedReceive<scalar>(commsType, coeffs.size())
procInterface_.compressedReceive<solveScalar>
(
commsType,
coeffs.size()
)
);
transformCoupleField(pnf, cmpt);

14
src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/processorGAMGInterfaceField/processorGAMGInterfaceField.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2014 OpenFOAM Foundation
@ -76,10 +76,10 @@ class processorGAMGInterfaceField
mutable label outstandingRecvRequest_;
//- Scalar send buffer
mutable Field<scalar> scalarSendBuf_;
mutable solveScalarField scalarSendBuf_;
//- Scalar receive buffer
mutable Field<scalar> scalarReceiveBuf_;
mutable solveScalarField scalarReceiveBuf_;
@ -139,9 +139,9 @@ public:
//- Initialise neighbour matrix update
virtual void initInterfaceMatrixUpdate
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
@ -150,9 +150,9 @@ public:
//- Update result field based on interface functionality
virtual void updateInterfaceMatrix
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType

61
src/OpenFOAM/matrices/lduMatrix/solvers/PBiCG/PBiCG.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2017 OpenFOAM Foundation
@ -26,6 +26,7 @@ License
\*---------------------------------------------------------------------------*/
#include "PBiCG.H"
#include "PrecisionAdaptor.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -66,11 +67,14 @@ Foam::PBiCG::PBiCG
Foam::solverPerformance Foam::PBiCG::solve
(
scalarField& psi,
scalarField& psi_s,
const scalarField& source,
const direction cmpt
) const
{
PrecisionAdaptor<solveScalar, scalar> tpsi(psi_s);
solveScalarField& psi = tpsi.constCast();
// --- Setup class containing solver performance data
solverPerformance solverPerf
(
@ -80,25 +84,31 @@ Foam::solverPerformance Foam::PBiCG::solve
const label nCells = psi.size();
scalar* __restrict__ psiPtr = psi.begin();
solveScalar* __restrict__ psiPtr = psi.begin();
scalarField pA(nCells);
scalar* __restrict__ pAPtr = pA.begin();
solveScalarField pA(nCells);
solveScalar* __restrict__ pAPtr = pA.begin();
scalarField wA(nCells);
scalar* __restrict__ wAPtr = wA.begin();
solveScalarField wA(nCells);
solveScalar* __restrict__ wAPtr = wA.begin();
// --- Calculate A.psi
matrix_.Amul(wA, psi, interfaceBouCoeffs_, interfaces_, cmpt);
// --- Calculate initial residual field
scalarField rA(source - wA);
scalar* __restrict__ rAPtr = rA.begin();
ConstPrecisionAdaptor<solveScalar, scalar> tsource(source);
solveScalarField rA(tsource() - wA);
solveScalar* __restrict__ rAPtr = rA.begin();
matrix().setResidualField(rA, fieldName_, true);
matrix().setResidualField
(
ConstPrecisionAdaptor<scalar, solveScalar>(rA)(),
fieldName_,
false
);
// --- Calculate normalisation factor
const scalar normFactor = this->normFactor(psi, source, wA, pA);
const solveScalar normFactor = this->normFactor(psi, tsource(), wA, pA);
if (lduMatrix::debug >= 2)
{
@ -118,21 +128,21 @@ Foam::solverPerformance Foam::PBiCG::solve
|| !solverPerf.checkConvergence(tolerance_, relTol_)
)
{
scalarField pT(nCells, 0);
scalar* __restrict__ pTPtr = pT.begin();
solveScalarField pT(nCells, 0);
solveScalar* __restrict__ pTPtr = pT.begin();
scalarField wT(nCells);
scalar* __restrict__ wTPtr = wT.begin();
solveScalarField wT(nCells);
solveScalar* __restrict__ wTPtr = wT.begin();
// --- Calculate T.psi
matrix_.Tmul(wT, psi, interfaceIntCoeffs_, interfaces_, cmpt);
// --- Calculate initial transpose residual field
scalarField rT(source - wT);
scalar* __restrict__ rTPtr = rT.begin();
solveScalarField rT(tsource() - wT);
solveScalar* __restrict__ rTPtr = rT.begin();
// --- Initial value not used
scalar wArT = 0;
solveScalar wArT = 0;
// --- Select and construct the preconditioner
autoPtr<lduMatrix::preconditioner> preconPtr =
@ -146,7 +156,7 @@ Foam::solverPerformance Foam::PBiCG::solve
do
{
// --- Store previous wArT
const scalar wArTold = wArT;
const solveScalar wArTold = wArT;
// --- Precondition residuals
preconPtr->precondition(wA, rA, cmpt);
@ -165,7 +175,7 @@ Foam::solverPerformance Foam::PBiCG::solve
}
else
{
const scalar beta = wArT/wArTold;
const solveScalar beta = wArT/wArTold;
for (label cell=0; cell<nCells; cell++)
{
@ -179,7 +189,7 @@ Foam::solverPerformance Foam::PBiCG::solve
matrix_.Amul(wA, pA, interfaceBouCoeffs_, interfaces_, cmpt);
matrix_.Tmul(wT, pT, interfaceIntCoeffs_, interfaces_, cmpt);
const scalar wApT = gSumProd(wA, pT, matrix().mesh().comm());
const solveScalar wApT = gSumProd(wA, pT, matrix().mesh().comm());
// --- Test for singularity
if (solverPerf.checkSingularity(mag(wApT)/normFactor))
@ -190,7 +200,7 @@ Foam::solverPerformance Foam::PBiCG::solve
// --- Update solution and residual:
const scalar alpha = wArT/wApT;
const solveScalar alpha = wArT/wApT;
for (label cell=0; cell<nCells; cell++)
{
@ -222,7 +232,12 @@ Foam::solverPerformance Foam::PBiCG::solve
<< exit(FatalError);
}
matrix().setResidualField(rA, fieldName_, false);
matrix().setResidualField
(
ConstPrecisionAdaptor<scalar, solveScalar>(rA)(),
fieldName_,
false
);
return solverPerf;
}

91
src/OpenFOAM/matrices/lduMatrix/solvers/PBiCGStab/PBiCGStab.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2016-2017 OpenFOAM Foundation
@ -26,6 +26,7 @@ License
\*---------------------------------------------------------------------------*/
#include "PBiCGStab.H"
#include "PrecisionAdaptor.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -67,10 +68,10 @@ Foam::PBiCGStab::PBiCGStab
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::solverPerformance Foam::PBiCGStab::solve
Foam::solverPerformance Foam::PBiCGStab::scalarSolve
(
scalarField& psi,
const scalarField& source,
solveScalarField& psi,
const solveScalarField& source,
const direction cmpt
) const
{
@ -83,25 +84,30 @@ Foam::solverPerformance Foam::PBiCGStab::solve
const label nCells = psi.size();
scalar* __restrict__ psiPtr = psi.begin();
solveScalar* __restrict__ psiPtr = psi.begin();
scalarField pA(nCells);
scalar* __restrict__ pAPtr = pA.begin();
solveScalarField pA(nCells);
solveScalar* __restrict__ pAPtr = pA.begin();
scalarField yA(nCells);
scalar* __restrict__ yAPtr = yA.begin();
solveScalarField yA(nCells);
solveScalar* __restrict__ yAPtr = yA.begin();
// --- Calculate A.psi
matrix_.Amul(yA, psi, interfaceBouCoeffs_, interfaces_, cmpt);
// --- Calculate initial residual field
scalarField rA(source - yA);
scalar* __restrict__ rAPtr = rA.begin();
solveScalarField rA(source - yA);
solveScalar* __restrict__ rAPtr = rA.begin();
matrix().setResidualField(rA, fieldName_, true);
matrix().setResidualField
(
ConstPrecisionAdaptor<scalar, solveScalar>(rA)(),
fieldName_,
true
);
// --- Calculate normalisation factor
const scalar normFactor = this->normFactor(psi, source, yA, pA);
const solveScalar normFactor = this->normFactor(psi, source, yA, pA);
if (lduMatrix::debug >= 2)
{
@ -121,25 +127,25 @@ Foam::solverPerformance Foam::PBiCGStab::solve
|| !solverPerf.checkConvergence(tolerance_, relTol_)
)
{
scalarField AyA(nCells);
scalar* __restrict__ AyAPtr = AyA.begin();
solveScalarField AyA(nCells);
solveScalar* __restrict__ AyAPtr = AyA.begin();
scalarField sA(nCells);
scalar* __restrict__ sAPtr = sA.begin();
solveScalarField sA(nCells);
solveScalar* __restrict__ sAPtr = sA.begin();
scalarField zA(nCells);
scalar* __restrict__ zAPtr = zA.begin();
solveScalarField zA(nCells);
solveScalar* __restrict__ zAPtr = zA.begin();
scalarField tA(nCells);
scalar* __restrict__ tAPtr = tA.begin();
solveScalarField tA(nCells);
solveScalar* __restrict__ tAPtr = tA.begin();
// --- Store initial residual
const scalarField rA0(rA);
const solveScalarField rA0(rA);
// --- Initial values not used
scalar rA0rA = 0;
scalar alpha = 0;
scalar omega = 0;
solveScalar rA0rA = 0;
solveScalar alpha = 0;
solveScalar omega = 0;
// --- Select and construct the preconditioner
autoPtr<lduMatrix::preconditioner> preconPtr =
@ -153,7 +159,7 @@ Foam::solverPerformance Foam::PBiCGStab::solve
do
{
// --- Store previous rA0rA
const scalar rA0rAold = rA0rA;
const solveScalar rA0rAold = rA0rA;
rA0rA = gSumProd(rA0, rA, matrix().mesh().comm());
@ -179,7 +185,7 @@ Foam::solverPerformance Foam::PBiCGStab::solve
break;
}
const scalar beta = (rA0rA/rA0rAold)*(alpha/omega);
const solveScalar beta = (rA0rA/rA0rAold)*(alpha/omega);
for (label cell=0; cell<nCells; cell++)
{
@ -194,7 +200,8 @@ Foam::solverPerformance Foam::PBiCGStab::solve
// --- Calculate AyA
matrix_.Amul(AyA, yA, interfaceBouCoeffs_, interfaces_, cmpt);
const scalar rA0AyA = gSumProd(rA0, AyA, matrix().mesh().comm());
const solveScalar rA0AyA =
gSumProd(rA0, AyA, matrix().mesh().comm());
alpha = rA0rA/rA0AyA;
@ -230,7 +237,7 @@ Foam::solverPerformance Foam::PBiCGStab::solve
// --- Calculate tA
matrix_.Amul(tA, zA, interfaceBouCoeffs_, interfaces_, cmpt);
const scalar tAtA = gSumSqr(tA, matrix().mesh().comm());
const solveScalar tAtA = gSumSqr(tA, matrix().mesh().comm());
// --- Calculate omega from tA and sA
// (cheaper than using zA with preconditioned tA)
@ -256,10 +263,34 @@ Foam::solverPerformance Foam::PBiCGStab::solve
);
}
matrix().setResidualField(rA, fieldName_, false);
matrix().setResidualField
(
ConstPrecisionAdaptor<scalar, solveScalar>(rA)(),
fieldName_,
false
);
return solverPerf;
}
Foam::solverPerformance Foam::PBiCGStab::solve
(
scalarField& psi_s,
const scalarField& source,
const direction cmpt
) const
{
PrecisionAdaptor<solveScalar, scalar> tpsi(psi_s);
solveScalarField& psi = tpsi.constCast();
return scalarSolve
(
psi,
ConstPrecisionAdaptor<solveScalar, scalar>(source)(),
cmpt
);
}
// ************************************************************************* //

8
src/OpenFOAM/matrices/lduMatrix/solvers/PBiCGStab/PBiCGStab.H
View File

@ -104,6 +104,14 @@ public:
// Member Functions
//- Solve the matrix with this solver
virtual solverPerformance scalarSolve
(
solveScalarField& psi,
const solveScalarField& source,
const direction cmpt = 0
) const;
//- Solve the matrix with this solver
virtual solverPerformance solve
(

68
src/OpenFOAM/matrices/lduMatrix/solvers/PCG/PCG.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2017 OpenFOAM Foundation
@ -26,6 +26,7 @@ License
\*---------------------------------------------------------------------------*/
#include "PCG.H"
#include "PrecisionAdaptor.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -64,10 +65,10 @@ Foam::PCG::PCG
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::solverPerformance Foam::PCG::solve
Foam::solverPerformance Foam::PCG::scalarSolve
(
scalarField& psi,
const scalarField& source,
solveScalarField& psi,
const solveScalarField& source,
const direction cmpt
) const
{
@ -80,28 +81,33 @@ Foam::solverPerformance Foam::PCG::solve
label nCells = psi.size();
scalar* __restrict__ psiPtr = psi.begin();
solveScalar* __restrict__ psiPtr = psi.begin();
scalarField pA(nCells);
scalar* __restrict__ pAPtr = pA.begin();
solveScalarField pA(nCells);
solveScalar* __restrict__ pAPtr = pA.begin();
scalarField wA(nCells);
scalar* __restrict__ wAPtr = wA.begin();
solveScalarField wA(nCells);
solveScalar* __restrict__ wAPtr = wA.begin();
scalar wArA = solverPerf.great_;
scalar wArAold = wArA;
solveScalar wArA = solverPerf.great_;
solveScalar wArAold = wArA;
// --- Calculate A.psi
matrix_.Amul(wA, psi, interfaceBouCoeffs_, interfaces_, cmpt);
// --- Calculate initial residual field
scalarField rA(source - wA);
scalar* __restrict__ rAPtr = rA.begin();
solveScalarField rA(source - wA);
solveScalar* __restrict__ rAPtr = rA.begin();
matrix().setResidualField(rA, fieldName_, true);
matrix().setResidualField
(
ConstPrecisionAdaptor<scalar, solveScalar>(rA)(),
fieldName_,
false
);
// --- Calculate normalisation factor
scalar normFactor = this->normFactor(psi, source, wA, pA);
solveScalar normFactor = this->normFactor(psi, source, wA, pA);
if (lduMatrix::debug >= 2)
{
@ -150,7 +156,7 @@ Foam::solverPerformance Foam::PCG::solve
}
else
{
scalar beta = wArA/wArAold;
solveScalar beta = wArA/wArAold;
for (label cell=0; cell<nCells; cell++)
{
@ -162,8 +168,7 @@ Foam::solverPerformance Foam::PCG::solve
// --- Update preconditioned residual
matrix_.Amul(wA, pA, interfaceBouCoeffs_, interfaces_, cmpt);
scalar wApA = gSumProd(wA, pA, matrix().mesh().comm());
solveScalar wApA = gSumProd(wA, pA, matrix().mesh().comm());
// --- Test for singularity
if (solverPerf.checkSingularity(mag(wApA)/normFactor)) break;
@ -171,7 +176,7 @@ Foam::solverPerformance Foam::PCG::solve
// --- Update solution and residual:
scalar alpha = wArA/wApA;
solveScalar alpha = wArA/wApA;
for (label cell=0; cell<nCells; cell++)
{
@ -193,10 +198,33 @@ Foam::solverPerformance Foam::PCG::solve
);
}
matrix().setResidualField(rA, fieldName_, false);
matrix().setResidualField
(
ConstPrecisionAdaptor<scalar, solveScalar>(rA)(),
fieldName_,
false
);
return solverPerf;
}
Foam::solverPerformance Foam::PCG::solve
(
scalarField& psi_s,
const scalarField& source,
const direction cmpt
) const
{
PrecisionAdaptor<solveScalar, scalar> tpsi(psi_s);
return scalarSolve
(
tpsi.constCast(),
ConstPrecisionAdaptor<solveScalar, scalar>(source)(),
cmpt
);
}
// ************************************************************************* //

8
src/OpenFOAM/matrices/lduMatrix/solvers/PCG/PCG.H
View File

@ -92,6 +92,14 @@ public:
// Member Functions
//- Solve the matrix with this solver
virtual solverPerformance scalarSolve
(
solveScalarField& psi,
const solveScalarField& source,
const direction cmpt=0
) const;
//- Solve the matrix with this solver
virtual solverPerformance solve
(

36
src/OpenFOAM/matrices/lduMatrix/solvers/smoothSolver/smoothSolver.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2016-2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2014 OpenFOAM Foundation
@ -27,6 +27,7 @@ License
#include "smoothSolver.H"
#include "profiling.H"
#include "PrecisionAdaptor.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -79,11 +80,14 @@ void Foam::smoothSolver::readControls()
Foam::solverPerformance Foam::smoothSolver::solve
(
scalarField& psi,
scalarField& psi_s,
const scalarField& source,
const direction cmpt
) const
{
PrecisionAdaptor<solveScalar, scalar> tpsi(psi_s);
solveScalarField& psi = tpsi.constCast();
// Setup class containing solver performance data
solverPerformance solverPerf(typeName, fieldName_);
@ -114,22 +118,29 @@ Foam::solverPerformance Foam::smoothSolver::solve
}
else
{
scalar normFactor = 0;
scalarField residual;
solveScalar normFactor = 0;
solveScalarField residual;
ConstPrecisionAdaptor<solveScalar, scalar> tsource(source);
{
scalarField Apsi(psi.size());
scalarField temp(psi.size());
solveScalarField Apsi(psi.size());
solveScalarField temp(psi.size());
// Calculate A.psi
matrix_.Amul(Apsi, psi, interfaceBouCoeffs_, interfaces_, cmpt);
// Calculate normalisation factor
normFactor = this->normFactor(psi, source, Apsi, temp);
normFactor = this->normFactor(psi, tsource(), Apsi, temp);
residual = source - Apsi;
residual = tsource() - Apsi;
matrix().setResidualField(residual, fieldName_, true);
matrix().setResidualField
(
ConstPrecisionAdaptor<scalar, solveScalar>(residual)(),
fieldName_,
false
);
// Calculate residual magnitude
solverPerf.initialResidual() =
@ -197,7 +208,12 @@ Foam::solverPerformance Foam::smoothSolver::solve
);
}
matrix().setResidualField(residual, fieldName_, false);
matrix().setResidualField
(
ConstPrecisionAdaptor<scalar, solveScalar>(residual)(),
fieldName_,
false
);
}
return solverPerf;

4
src/OpenFOAM/primitives/Scalar/Scalar.H
View File

@ -52,10 +52,14 @@ public:
//- Component type
typedef Scalar cmptType;
//- Magnitude type
typedef Scalar magType;
//- Equivalent type of labels used for valid component indexing
typedef label labelType;
// Member constants
//- Dimensionality of space

3
src/OpenFOAM/primitives/VectorSpace/VectorSpace.H
View File

@ -88,6 +88,9 @@ public:
//- Component type
typedef Cmpt cmptType;
//- Magnitude type
typedef Cmpt magType;
// Static Constants

26
src/OpenFOAM/primitives/VectorSpace/products.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
@ -44,13 +44,13 @@ namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//- Abstract template class to provide the form resulting from
// the inner-product of two forms
//- the inner-product of two forms
template<class Cmpt, class Form1, class Form2>
class typeOfInnerProduct
{};
//- Abstract template class to provide the form resulting from
// the outer-product of two forms
//- the outer-product of two forms
template<class Cmpt, class Form1, class Form2>
class typeOfOuterProduct
{};
@ -71,6 +71,16 @@ class symmTypeOfRank
{};
//- The magnitude type for given argument.
template<class arg1>
class typeOfMag
{
public:
typedef typename pTraits<typename pTraits<arg1>::cmptType>::magType type;
};
template<class arg1, class arg2>
class typeOfSum
{
@ -93,6 +103,16 @@ public:
};
//- Outer-product of identical types
template<class arg1>
class outerProduct1
:
public outerProduct<arg1, arg1>
{
public:
};
template<class arg1, class arg2>
class crossProduct
{

6
src/OpenFOAM/primitives/bools/bool/bool.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
@ -74,6 +74,10 @@ public:
//- Component type
typedef bool cmptType;
//- Magnitude type
typedef bool magType;
// Member constants
//- Dimensionality of space

3
src/OpenFOAM/primitives/complex/complex.H
View File

@ -252,6 +252,9 @@ public:
//- Component type
typedef complex cmptType;
//- Magnitude type
typedef scalar magType;
//- Equivalent type of labels used for valid component indexing
typedef label labelType;

3
src/OpenFOAM/primitives/ints/int32/int32.H
View File

@ -138,6 +138,9 @@ public:
//- Component type
typedef int32_t cmptType;
//- Magnitude type
typedef int32_t magType;
// Member constants

3
src/OpenFOAM/primitives/ints/int64/int64.H
View File

@ -139,6 +139,9 @@ public:
//- Component type
typedef int64_t cmptType;
//- Magnitude type
typedef int64_t magType;
// Member constants

3
src/OpenFOAM/primitives/quaternion/quaternion.H
View File

@ -96,6 +96,9 @@ public:
//- Component type
typedef scalar cmptType;
//- Magnitude type
typedef scalar magType;
//- Euler-angle rotation order
enum eulerOrder : unsigned char
{

1
src/OpenFOAM/primitives/ranges/MinMax/MinMax.H
View File

@ -132,6 +132,7 @@ public:
//- The value type the MinMax represents
typedef T value_type;
//- Component type
typedef typename pTraits<T>::cmptType cmptType;

47
src/finiteArea/faMatrices/faMatrix/faMatrixSolve.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2016-2017 Wikki Ltd
@ -28,6 +28,10 @@ Description
\*---------------------------------------------------------------------------*/
#include "PrecisionAdaptor.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
void Foam::faMatrix<Type>::setComponentReference
(
@ -96,25 +100,30 @@ Foam::SolverPerformance<Type> Foam::faMatrix<Type>::solve
// Use the initMatrixInterfaces and updateMatrixInterfaces to correct
// bouCoeffsCmpt for the explicit part of the coupled boundary
// conditions
initMatrixInterfaces
(
true,
bouCoeffsCmpt,
interfaces,
psiCmpt,
sourceCmpt,
cmpt
);
{
PrecisionAdaptor<solveScalar, scalar> sourceCmpt_ss(sourceCmpt);
ConstPrecisionAdaptor<solveScalar, scalar> psiCmpt_ss(psiCmpt);
updateMatrixInterfaces
(
true,
bouCoeffsCmpt,
interfaces,
psiCmpt,
sourceCmpt,
cmpt
);
initMatrixInterfaces
(
true,
bouCoeffsCmpt,
interfaces,
psiCmpt_ss(),
sourceCmpt_ss.ref(),
cmpt
);
updateMatrixInterfaces
(
true,
bouCoeffsCmpt,
interfaces,
psiCmpt_ss(),
sourceCmpt_ss.ref(),
cmpt
);
}
solverPerformance solverPerf;

19
src/finiteArea/faMatrices/faScalarMatrix/faScalarMatrix.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2016-2017 Wikki Ltd
@ -30,6 +30,7 @@ Description
#include "faScalarMatrix.H"
#include "zeroGradientFaPatchFields.H"
#include "PrecisionAdaptor.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -101,21 +102,27 @@ Foam::tmp<Foam::scalarField> Foam::faMatrix<Foam::scalar>::residual() const
scalarField boundaryDiag(psi_.size(), Zero);
addBoundaryDiag(boundaryDiag, 0);
tmp<scalarField> tres
const scalarField& psif = psi_.internalField();
ConstPrecisionAdaptor<solveScalar, scalar> tpsi(psif);
const solveScalarField& psi = tpsi();
tmp<solveScalarField> tres
(
lduMatrix::residual
(
psi_.internalField(),
source_ - boundaryDiag*psi_.internalField(),
psi,
source_ - boundaryDiag*psif,
boundaryCoeffs_,
psi_.boundaryField().scalarInterfaces(),
0
)
);
addBoundarySource(tres.ref());
ConstPrecisionAdaptor<scalar, solveScalar> tres_s(tres);
return tres;
addBoundarySource(tres_s.ref());
return tres_s;
}

8
src/finiteArea/fields/faPatchFields/basic/coupled/coupledFaPatchField.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2016-2017 Wikki Ltd
@ -183,16 +183,16 @@ public:
//- Transform given patch component field
virtual void transformCoupleField
(
scalarField& f,
solveScalarField& f,
const direction cmpt
) const = 0;
//- Update result field based on interface functionality
virtual void updateInterfaceMatrix
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction,
const Pstream::commsTypes commsType

8
src/finiteArea/fields/faPatchFields/constraint/cyclic/cyclicFaPatchField.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2016-2017 Wikki Ltd
@ -161,15 +161,15 @@ Foam::cyclicFaPatchField<Type>::patchNeighbourField() const
template<class Type>
void Foam::cyclicFaPatchField<Type>::updateInterfaceMatrix
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
) const
{
scalarField pnf(this->size());
solveScalarField pnf(this->size());
label sizeby2 = this->size()/2;
const labelUList& faceCells = cyclicPatch_.faceCells();

8
src/finiteArea/fields/faPatchFields/constraint/cyclic/cyclicFaPatchField.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2016-2017 Wikki Ltd
@ -166,7 +166,7 @@ public:
//- Transform neighbour field
virtual void transformCoupleField
(
scalarField& f,
solveScalarField& f,
const direction cmpt
) const
{
@ -176,9 +176,9 @@ public:
//- Update result field based on interface functionality
virtual void updateInterfaceMatrix
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType

14
src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchField.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2016-2017 Wikki Ltd
@ -189,9 +189,9 @@ Foam::tmp<Foam::Field<Type>> Foam::processorFaPatchField<Type>::snGrad() const
template<class Type>
void Foam::processorFaPatchField<Type>::initInterfaceMatrixUpdate
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction,
const Pstream::commsTypes commsType
@ -208,17 +208,17 @@ void Foam::processorFaPatchField<Type>::initInterfaceMatrixUpdate
template<class Type>
void Foam::processorFaPatchField<Type>::updateInterfaceMatrix
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField&,
const solveScalarField&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
) const
{
scalarField pnf
solveScalarField pnf
(
procPatch_.receive<scalar>(commsType, this->size())()
procPatch_.receive<solveScalar>(commsType, this->size())()
);
// Transform according to the transformation tensor

12
src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchField.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2016-2017 Wikki Ltd
@ -180,7 +180,7 @@ public:
//- Transform neighbour field
virtual void transformCoupleField
(
scalarField& f,
solveScalarField& f,
const direction cmpt
) const
{
@ -190,9 +190,9 @@ public:
//- Initialise neighbour matrix update
virtual void initInterfaceMatrixUpdate
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
@ -201,9 +201,9 @@ public:
//- Update result field based on interface functionality
virtual void updateInterfaceMatrix
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType

16
src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchScalarField.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2016-2017 Wikki Ltd
@ -32,7 +32,7 @@ License
template<>
void Foam::processorFaPatchField<Foam::scalar>::transformCoupleField
(
scalarField& f,
solveScalarField& f,
const direction cmpt
) const
{}
@ -41,9 +41,9 @@ void Foam::processorFaPatchField<Foam::scalar>::transformCoupleField
template<>
void Foam::processorFaPatchField<Foam::scalar>::initInterfaceMatrixUpdate
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction,
const Pstream::commsTypes commsType
@ -60,17 +60,17 @@ void Foam::processorFaPatchField<Foam::scalar>::initInterfaceMatrixUpdate
template<>
void Foam::processorFaPatchField<Foam::scalar>::updateInterfaceMatrix
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField&,
const solveScalarField&,
const scalarField& coeffs,
const direction,
const Pstream::commsTypes commsType
) const
{
scalarField pnf
solveScalarField pnf
(
procPatch_.receive<scalar>(commsType, this->size())()
procPatch_.receive<solveScalar>(commsType, this->size())()
);
const labelUList& edgeFaces = patch().edgeFaces();

12
src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchScalarField.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2016-2017 Wikki Ltd
@ -40,7 +40,7 @@ namespace Foam
template<>
void processorFaPatchField<scalar>::transformCoupleField
(
scalarField& f,
solveScalarField& f,
const direction cmpt
) const;
@ -48,9 +48,9 @@ void processorFaPatchField<scalar>::transformCoupleField
template<>
void processorFaPatchField<scalar>::initInterfaceMatrixUpdate
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField&,
const solveScalarField&,
const scalarField& coeffs,
const direction,
const Pstream::commsTypes commsType
@ -60,9 +60,9 @@ void processorFaPatchField<scalar>::initInterfaceMatrixUpdate
template<>
void processorFaPatchField<scalar>::updateInterfaceMatrix
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField&,
const solveScalarField&,
const scalarField& coeffs,
const direction,
const Pstream::commsTypes commsType

6
src/finiteVolume/fields/fvPatchFields/basic/coupled/coupledFvPatchField.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
@ -208,9 +208,9 @@ public:
//- Update result field based on interface functionality
virtual void updateInterfaceMatrix
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction,
const Pstream::commsTypes commsType

8
src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2017 OpenFOAM Foundation
@ -172,9 +172,9 @@ Foam::cyclicFvPatchField<Type>::neighbourPatchField() const
template<class Type>
void Foam::cyclicFvPatchField<Type>::updateInterfaceMatrix
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes
@ -183,7 +183,7 @@ void Foam::cyclicFvPatchField<Type>::updateInterfaceMatrix
const labelUList& nbrFaceCells =
cyclicPatch().cyclicPatch().neighbPatch().faceCells();
scalarField pnf(psiInternal, nbrFaceCells);
solveScalarField pnf(psiInternal, nbrFaceCells);
// Transform according to the transformation tensors
transformCoupleField(pnf, cmpt);

6
src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
@ -179,9 +179,9 @@ public:
//- Update result field based on interface functionality
virtual void updateInterfaceMatrix
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType

8
src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2013-2017 OpenFOAM Foundation
@ -215,9 +215,9 @@ Foam::cyclicACMIFvPatchField<Type>::nonOverlapPatchField() const
template<class Type>
void Foam::cyclicACMIFvPatchField<Type>::updateInterfaceMatrix
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes
@ -230,7 +230,7 @@ void Foam::cyclicACMIFvPatchField<Type>::updateInterfaceMatrix
const labelUList& nbrFaceCellsCoupled =
cpp.neighbPatch().faceCells();
scalarField pnf(psiInternal, nbrFaceCellsCoupled);
solveScalarField pnf(psiInternal, nbrFaceCellsCoupled);
// Transform according to the transformation tensors
transformCoupleField(pnf, cmpt);

6
src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2013-2016 OpenFOAM Foundation
@ -203,9 +203,9 @@ public:
//- Update result field based on interface functionality
virtual void updateInterfaceMatrix
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType

10
src/finiteVolume/fields/fvPatchFields/constraint/cyclicAMI/cyclicAMIFvPatchField.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2017 OpenFOAM Foundation
@ -178,9 +178,9 @@ Foam::cyclicAMIFvPatchField<Type>::neighbourPatchField() const
template<class Type>
void Foam::cyclicAMIFvPatchField<Type>::updateInterfaceMatrix
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes
@ -189,14 +189,14 @@ void Foam::cyclicAMIFvPatchField<Type>::updateInterfaceMatrix
const labelUList& nbrFaceCells =
cyclicAMIPatch_.cyclicAMIPatch().neighbPatch().faceCells();
scalarField pnf(psiInternal, nbrFaceCells);
solveScalarField pnf(psiInternal, nbrFaceCells);
// Transform according to the transformation tensors
transformCoupleField(pnf, cmpt);
if (cyclicAMIPatch_.applyLowWeightCorrection())
{
scalarField pif(psiInternal, cyclicAMIPatch_.faceCells());
solveScalarField pif(psiInternal, cyclicAMIPatch_.faceCells());
pnf = cyclicAMIPatch_.interpolate(pnf, pif);
}
else

6
src/finiteVolume/fields/fvPatchFields/constraint/cyclicAMI/cyclicAMIFvPatchField.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
@ -184,9 +184,9 @@ public:
//- Update result field based on interface functionality
virtual void updateInterfaceMatrix
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType

6
src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2017 OpenFOAM Foundation
@ -133,9 +133,9 @@ Foam::jumpCyclicFvPatchField<Type>::patchNeighbourField() const
template<class Type>
void Foam::jumpCyclicFvPatchField<Type>::updateInterfaceMatrix
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes

14
src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2018 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2018-2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
@ -129,9 +129,9 @@ public:
//- Update result field based on interface functionality
virtual void updateInterfaceMatrix
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
@ -152,9 +152,9 @@ public:
template<>
void jumpCyclicFvPatchField<scalar>::updateInterfaceMatrix
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
@ -164,9 +164,9 @@ void jumpCyclicFvPatchField<scalar>::updateInterfaceMatrix
template<>
void jumpCyclicFvPatchField<vector>::updateInterfaceMatrix
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType

26
src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchFields.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2018 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2018-2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
@ -41,21 +41,25 @@ namespace Foam
template<>
void Foam::jumpCyclicFvPatchField<Foam::scalar>::updateInterfaceMatrix
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes
) const
{
scalarField pnf(this->size());
solveScalarField pnf(this->size());
const labelUList& nbrFaceCells =
this->cyclicPatch().neighbFvPatch().faceCells();
// only apply jump to original field
if (&psiInternal == &this->primitiveField())
if
(
reinterpret_cast<const void*>(&psiInternal)
== reinterpret_cast<const void*>(&this->primitiveField())
)
{
Field<scalar> jf(this->jump());
@ -88,15 +92,15 @@ void Foam::jumpCyclicFvPatchField<Foam::scalar>::updateInterfaceMatrix
template<>
void Foam::jumpCyclicFvPatchField<Foam::vector>::updateInterfaceMatrix
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes
) const
{
scalarField pnf(this->size());
solveScalarField pnf(this->size());
const labelUList& nbrFaceCells =
this->cyclicPatch().neighbFvPatch().faceCells();
@ -104,7 +108,11 @@ void Foam::jumpCyclicFvPatchField<Foam::vector>::updateInterfaceMatrix
const Field<vector>& iField = this->primitiveField();
// only apply jump to original field
if (&psiInternal == &(iField.component(cmpt).ref()))
if
(
reinterpret_cast<const void*>(&psiInternal)
== reinterpret_cast<const void*>(&(iField.component(cmpt).ref()))
)
{
Field<vector> jf(this->jump());

6
src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchField.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2012-2017 OpenFOAM Foundation
@ -135,9 +135,9 @@ Foam::jumpCyclicAMIFvPatchField<Type>::patchNeighbourField() const
template<class Type>
void Foam::jumpCyclicAMIFvPatchField<Type>::updateInterfaceMatrix
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes

10
src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchField.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2012-2016 OpenFOAM Foundation
@ -132,9 +132,9 @@ public:
//- Update result field based on interface functionality
virtual void updateInterfaceMatrix
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
@ -156,9 +156,9 @@ public:
template<>
void jumpCyclicAMIFvPatchField<scalar>::updateInterfaceMatrix
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType

20
src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclicAMI/jumpCyclicAMIFvPatchFields.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2012-2016 OpenFOAM Foundation
@ -43,9 +43,9 @@ makePatchFieldTypeNames(jumpCyclicAMI);
template<>
void Foam::jumpCyclicAMIFvPatchField<scalar>::updateInterfaceMatrix
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes
@ -54,12 +54,16 @@ void Foam::jumpCyclicAMIFvPatchField<scalar>::updateInterfaceMatrix
const labelUList& nbrFaceCells =
this->cyclicAMIPatch().cyclicAMIPatch().neighbPatch().faceCells();
scalarField pnf(psiInternal, nbrFaceCells);
solveScalarField pnf(psiInternal, nbrFaceCells);
pnf = this->cyclicAMIPatch().interpolate(pnf);
// only apply jump to original field
if (&psiInternal == &this->primitiveField())
if
(
reinterpret_cast<const void*>(&psiInternal)
== reinterpret_cast<const void*>(&this->primitiveField())
)
{
Field<scalar> jf(this->jump());
@ -68,7 +72,11 @@ void Foam::jumpCyclicAMIFvPatchField<scalar>::updateInterfaceMatrix
jf *= -1.0;
}
pnf -= jf;
//pnf -= jf;
forAll(pnf, i)
{
pnf[i] -= jf[i];
}
}
// Transform according to the transformation tensors

16
src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchField.C
View File

@ -310,9 +310,9 @@ Foam::processorFvPatchField<Type>::snGrad
template<class Type>
void Foam::processorFvPatchField<Type>::initInterfaceMatrixUpdate
(
scalarField&,
solveScalarField&,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField&,
const direction,
const Pstream::commsTypes commsType
@ -371,9 +371,9 @@ void Foam::processorFvPatchField<Type>::initInterfaceMatrixUpdate
template<class Type>
void Foam::processorFvPatchField<Type>::updateInterfaceMatrix
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField&,
const solveScalarField&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
@ -415,9 +415,13 @@ void Foam::processorFvPatchField<Type>::updateInterfaceMatrix
}
else
{
scalarField pnf
solveScalarField pnf
(
procPatch_.compressedReceive<scalar>(commsType, this->size())()
procPatch_.compressedReceive<solveScalar>
(
commsType,
this->size()
)()
);
if (!std::is_arithmetic<Type>::value)

15
src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchField.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
@ -88,10 +88,11 @@ class processorFvPatchField
mutable label outstandingRecvRequest_;
//- Scalar send buffer
mutable Field<scalar> scalarSendBuf_;
mutable solveScalarField scalarSendBuf_;
//- Scalar receive buffer
mutable Field<scalar> scalarReceiveBuf_;
mutable solveScalarField scalarReceiveBuf_;
public:
@ -210,9 +211,9 @@ public:
//- Initialise neighbour matrix update
virtual void initInterfaceMatrixUpdate
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
@ -221,9 +222,9 @@ public:
//- Update result field based on interface functionality
virtual void updateInterfaceMatrix
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType

3
src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2017 OpenFOAM Foundation
@ -49,6 +49,7 @@ SourceFiles
#include "fvPatch.H"
#include "DimensionedField.H"
#include "fieldTypes.H"
#include "scalarField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

44
src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2016-2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2016 OpenFOAM Foundation
@ -28,6 +28,7 @@ License
#include "LduMatrix.H"
#include "diagTensorField.H"
#include "profiling.H"
#include "PrecisionAdaptor.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -173,25 +174,30 @@ Foam::SolverPerformance<Type> Foam::fvMatrix<Type>::solveSegregated
// Use the initMatrixInterfaces and updateMatrixInterfaces to correct
// bouCoeffsCmpt for the explicit part of the coupled boundary
// conditions
initMatrixInterfaces
(
true,
bouCoeffsCmpt,
interfaces,
psiCmpt,
sourceCmpt,
cmpt
);
{
PrecisionAdaptor<solveScalar, scalar> sourceCmpt_ss(sourceCmpt);
ConstPrecisionAdaptor<solveScalar, scalar> psiCmpt_ss(psiCmpt);
updateMatrixInterfaces
(
true,
bouCoeffsCmpt,
interfaces,
psiCmpt,
sourceCmpt,
cmpt
);
initMatrixInterfaces
(
true,
bouCoeffsCmpt,
interfaces,
psiCmpt_ss(),
sourceCmpt_ss.constCast(),
cmpt
);
updateMatrixInterfaces
(
true,
bouCoeffsCmpt,
interfaces,
psiCmpt_ss(),
sourceCmpt_ss.constCast(),
cmpt
);
}
solverPerformance solverPerf;

20
src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2016-2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2017 OpenFOAM Foundation
@ -28,6 +28,7 @@ License
#include "fvScalarMatrix.H"
#include "extrapolatedCalculatedFvPatchFields.H"
#include "profiling.H"
#include "PrecisionAdaptor.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -199,21 +200,24 @@ Foam::tmp<Foam::scalarField> Foam::fvMatrix<Foam::scalar>::residual() const
scalarField boundaryDiag(psi_.size(), Zero);
addBoundaryDiag(boundaryDiag, 0);
tmp<scalarField> tres
const scalarField& psif = psi_.primitiveField();
ConstPrecisionAdaptor<solveScalar, scalar> tpsi(psif);
const solveScalarField& psi = tpsi();
tmp<solveScalarField> tres
(
lduMatrix::residual
(
psi_.primitiveField(),
source_ - boundaryDiag*psi_.primitiveField(),
psi,
source_ - boundaryDiag*psif,
boundaryCoeffs_,
psi_.boundaryField().scalarInterfaces(),
0
)
);
addBoundarySource(tres.ref());
return tres;
ConstPrecisionAdaptor<scalar, solveScalar> tres_s(tres);
addBoundarySource(tres_s.constCast());
return tres_s;
}

8
src/meshTools/AMIInterpolation/GAMG/interfaceFields/cyclicACMIGAMGInterfaceField/cyclicACMIGAMGInterfaceField.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2013 OpenFOAM Foundation
@ -94,16 +94,16 @@ Foam::cyclicACMIGAMGInterfaceField::~cyclicACMIGAMGInterfaceField()
void Foam::cyclicACMIGAMGInterfaceField::updateInterfaceMatrix
(
scalarField& result,
solveScalarField& result,
const bool add,
const scalarField& psiInternal,
const solveScalarField& psiInternal,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes
) const
{
// Get neighbouring field
scalarField pnf
solveScalarField pnf
(
cyclicACMIInterface_.neighbPatch().interfaceInternalField(psiInternal)
);

Some files were not shown because too many files have changed in this diff Show More