Merge branch 'nextRelease' of ssh://git.code.sf.net/p/foam-extend/foam-extend-4.0 into nextRelease

This commit is contained in:
Hrvoje Jasak 2019-05-01 13:00:40 +01:00
commit 19b5b8cfec
11 changed files with 260 additions and 26 deletions

View file

@ -326,7 +326,7 @@ void Foam::sigFpe::set(const bool verbose)
{ {
vm_protect( vm_protect(
mach_task_self(), mach_task_self(),
(uintptr_t)zone, reinterpret_cast<uintptr_t>(zone),
sizeof(malloc_zone_t), sizeof(malloc_zone_t),
0, 0,
VM_PROT_READ | VM_PROT_WRITE VM_PROT_READ | VM_PROT_WRITE
@ -338,7 +338,7 @@ void Foam::sigFpe::set(const bool verbose)
{ {
vm_protect( vm_protect(
mach_task_self(), mach_task_self(),
(uintptr_t)zone, reinterpret_cast<uintptr_t>(zone),
sizeof(malloc_zone_t), sizeof(malloc_zone_t),
0, 0,
VM_PROT_READ VM_PROT_READ

View file

@ -100,12 +100,14 @@ Foam::coupledFvMatrix<Type>::solve(const dictionary& solverControls)
fvMatrix<Type>& curMatrix = fvMatrix<Type>& curMatrix =
static_cast<fvMatrix<Type>& >(matrices[rowI]); static_cast<fvMatrix<Type>& >(matrices[rowI]);
// HR 12.03.19: Complete assembly before making copies.
curMatrix.completeAssembly();
saveDiag.set(rowI, new scalarField(curMatrix.diag())); saveDiag.set(rowI, new scalarField(curMatrix.diag()));
psiCmpt.set(rowI, new scalarField(curMatrix.psi().size())); psiCmpt.set(rowI, new scalarField(curMatrix.psi().size()));
source.set(rowI, new Field<Type>(curMatrix.source())); source.set(rowI, new Field<Type>(curMatrix.source()));
sourceCmpt.set(rowI, new scalarField(curMatrix.psi().size())); sourceCmpt.set(rowI, new scalarField(curMatrix.psi().size()));
curMatrix.completeAssembly();
curMatrix.addBoundarySource(source[rowI]); curMatrix.addBoundarySource(source[rowI]);
interfaces[rowI] = curMatrix.psi().boundaryField().interfaces(); interfaces[rowI] = curMatrix.psi().boundaryField().interfaces();

View file

@ -23,9 +23,15 @@ License
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#include "NamedEnum.H"
#include "stringList.H" #include "stringList.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
template<class Enum, int nEnum>
const char* Foam::NamedEnum<Enum, nEnum>::names[nEnum];
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Enum, int nEnum> template<class Enum, int nEnum>

View file

@ -64,7 +64,7 @@ void Foam::radiation::viewFactor::initialise()
} }
} }
selectedPatches_.resize(count--); selectedPatches_.resize(count);
if (debug) if (debug)
{ {
@ -173,6 +173,82 @@ void Foam::radiation::viewFactor::initialise()
globalIndex globalNumbering(nLocalCoarseFaces_); globalIndex globalNumbering(nLocalCoarseFaces_);
// Calculate global face areas
scalarField compactCoarseMagSf(map_->constructSize(), 0.0);
DynamicList<scalar> localCoarseMagSf(nLocalCoarseFaces_);
forAll (selectedPatches_, i)
{
label patchID = selectedPatches_[i];
const scalarField& sf = mesh_.magSf().boundaryField()[patchID];
const polyPatch& pp = coarseMesh_.boundaryMesh()[patchID];
const labelList& coarsePatchFace = coarseMesh_.patchFaceMap()[patchID];
scalarList magSf(pp.size(), 0.0);
if (pp.size() > 0)
{
const labelList& agglom = finalAgglom_[patchID];
label nAgglom = max(agglom) + 1;
labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
forAll (coarseToFine, coarseI)
{
const label coarseFaceID = coarsePatchFace[coarseI];
const labelList& fineFaces = coarseToFine[coarseFaceID];
UIndirectList<scalar> fineSf
(
sf,
fineFaces
);
magSf[coarseI] = sum(fineSf());
}
}
localCoarseMagSf.append(magSf);
}
// Fill the local values to distribute
SubList<scalar>(compactCoarseMagSf, nLocalCoarseFaces_).assign(localCoarseMagSf);
// Distribute data
map_->distribute(compactCoarseMagSf);
// Distribute local global ID
labelList compactGlobalIds(map_->constructSize(), 0.0);
labelList localGlobalIds(nLocalCoarseFaces_);
for(label k = 0; k < nLocalCoarseFaces_; k++)
{
localGlobalIds[k] = globalNumbering.toGlobal(Pstream::myProcNo(), k);
}
SubList<label>
(
compactGlobalIds,
nLocalCoarseFaces_
).assign(localGlobalIds);
map_->distribute(compactGlobalIds);
// Create global size vectors
scalarField coarseMagSf(totalNCoarseFaces_, 0.0);
// Fill lists from compact to global indexes.
forAll (compactCoarseMagSf, i)
{
coarseMagSf[compactGlobalIds[i]] = compactCoarseMagSf[i];
}
Pstream::listCombineGather(coarseMagSf, maxEqOp<scalar>());
Pstream::listCombineScatter(coarseMagSf);
if (Pstream::master()) if (Pstream::master())
{ {
Fmatrix_.reset Fmatrix_.reset
@ -200,18 +276,173 @@ void Foam::radiation::viewFactor::initialise()
{ {
Info<< "Smoothing the matrix..." << endl; Info<< "Smoothing the matrix..." << endl;
for (label i=0; i<totalNCoarseFaces_; i++) // HR 19.02.19: Old smoothing messed up the reciprocity condition:
// A_i*F_ij = A_j*F_ji. The new algorithm distributes the error
// in two passes moving downwards and upwards through the rows.
// On the way down, only the upper triangular is updated from the
// closedness condition: sum_j(F_ij) = 1; the lower triangular is
// updated from reciprocity. The reverse happen on the way up.
if(debug)
{ {
scalar sumF = 0.0; scalar maxErr = -GREAT;
for (label j=0; j<totalNCoarseFaces_; j++) scalar minSumF = GREAT;
scalar maxSumF = -GREAT;
for (label i=0; i<totalNCoarseFaces_; i++)
{ {
sumF += Fmatrix_()[i][j]; scalar sumF = 0.0;
for (label j=0; j<totalNCoarseFaces_; j++)
{
sumF += Fmatrix_()[i][j];
}
maxErr = max(maxErr, mag(sumF - 1.0));
minSumF = min(minSumF, sumF);
maxSumF = max(maxSumF, sumF);
} }
scalar delta = 1.0 - sumF;
for (label j=0; j<totalNCoarseFaces_; j++) Info << "max err = " << maxErr
<< " min sumF = " << minSumF
<< " max sumF = " << maxSumF << endl;
}
label iter = 0;
label maxIter = 20;
while(true)
{
iter++;
for (label i=0; i<totalNCoarseFaces_; i++)
{ {
Fmatrix_()[i][j] *= (1.0 - delta/(sumF + 0.001)); scalar sumF1 = 0.0;
for (label j=0; j<i; j++)
{
sumF1 += Fmatrix_()[i][j];
}
if(sumF1 > 1.0)
{
continue;
}
scalar sumF2 = 0.0;
for (label j=i+1; j<totalNCoarseFaces_; j++)
{
sumF2 += Fmatrix_()[i][j];
}
if(sumF2 < SMALL)
{
continue;
}
const scalar corr = (1.0 - sumF1)/sumF2;
const scalar magSfi = coarseMagSf[i];
for (label j=i+1; j<totalNCoarseFaces_; j++)
{
Fmatrix_()[i][j] *= corr;
Fmatrix_()[j][i] = Fmatrix_()[i][j]*magSfi/coarseMagSf[j];
}
} }
{
scalar maxErr = -GREAT;
scalar minSumF = GREAT;
scalar maxSumF = -GREAT;
for (label i=0; i<totalNCoarseFaces_; i++)
{
scalar sumF = 0.0;
for (label j=0; j<totalNCoarseFaces_; j++)
{
sumF += Fmatrix_()[i][j];
}
maxErr = max(maxErr, mag(sumF - 1.0));
minSumF = min(minSumF, sumF);
maxSumF = max(maxSumF, sumF);
}
if(debug)
{
Info << "max err = " << maxErr
<< " min sumF = " << minSumF
<< " max sumF = " << maxSumF << endl;
}
if(maxErr < 1e-6 || iter == maxIter)
{
break;
}
}
iter++;
for (label i=totalNCoarseFaces_-1; i>=0; i--)
{
scalar sumF1 = 0.0;
for (label j=0; j<i; j++)
{
sumF1 += Fmatrix_()[i][j];
}
if(sumF1 < SMALL)
{
continue;
}
scalar sumF2 = 0.0;
for (label j=i+1; j<totalNCoarseFaces_; j++)
{
sumF2 += Fmatrix_()[i][j];
}
if(sumF2 > 1.0)
{
continue;
}
const scalar corr = (1.0 - sumF2)/sumF1;
const scalar magSfi = coarseMagSf[i];
for (label j=0; j<i; j++)
{
Fmatrix_()[i][j] *= corr;
Fmatrix_()[j][i] = Fmatrix_()[i][j]*magSfi/coarseMagSf[j];
}
}
{
scalar maxErr = -GREAT;
scalar minSumF = GREAT;
scalar maxSumF = -GREAT;
for (label i=0; i<totalNCoarseFaces_; i++)
{
scalar sumF = 0.0;
for (label j=0; j<totalNCoarseFaces_; j++)
{
sumF += Fmatrix_()[i][j];
}
maxErr = max(maxErr, mag(sumF - 1.0));
minSumF = min(minSumF, sumF);
maxSumF = max(maxSumF, sumF);
}
if(debug)
{
Info << "max err = " << maxErr
<< " min sumF = " << minSumF
<< " max sumF = " << maxSumF << endl;
}
if(maxErr < 1e-6 || iter == maxIter)
{
break;
}
}
}
if(iter != maxIter)
{
Info << "converged in " << iter << " iterations" << endl;
}
else
{
Info << "not converged in " << iter << " iterations" << endl;
} }
} }

View file

@ -9,5 +9,5 @@ if [ "$WM_PROJECT" = "OpenFOAM" ]
then then
runApplication reconstructParMesh -constant -fullMatch runApplication reconstructParMesh -constant -fullMatch
else else
runApplication reconstructParMesh -zeroTime runApplication reconstructParMesh -withZero
fi fi

View file

@ -1,8 +1,8 @@
.SUFFIXES: .C .cxx .cc .cpp .SUFFIXES: .C .cxx .cc .cpp
c++WARN = -Wall -Wextra -Wno-unused-parameter -Wold-style-cast -Wnon-virtual-dtor -Wno-overloaded-virtual -Wno-unused-comparison c++WARN = -Wall -Wextra -Wno-unused-parameter -Wold-style-cast -Wno-overloaded-virtual -Wno-unused-comparison -Wno-undefined-var-template -Wno-null-dereference
CC = $(WM_CXX) -m64 -ftrapping-math CC = $(WM_CXX) -m64 -std=c++11 -ftrapping-math
# -fsignaling-nans # -fsignaling-nans
include $(RULES)/c++$(WM_COMPILE_OPTION) include $(RULES)/c++$(WM_COMPILE_OPTION)

View file

@ -1,5 +1,2 @@
c++DBUG = c++DBUG =
# c++OPT = -O4 c++OPT = -O1 -ftree-vectorize
c++OPT = -O3
#c++OPT = -march=nocona -O3
# -ftree-vectorize -ftree-vectorizer-verbose=3

View file

@ -1,8 +1,8 @@
.SUFFIXES: .C .cxx .cc .cpp .SUFFIXES: .C .cxx .cc .cpp
c++WARN = -Wall -Wextra -Wno-unused-parameter -Wold-style-cast -Wno-overloaded-virtual -Wno-unsequenced -Wno-c++11-extensions -Wno-unused-comparison -Wno-undefined-var-template -Wno-null-dereference c++WARN = -Wall -Wextra -Wno-unused-parameter -Wold-style-cast -Wno-overloaded-virtual -Wno-unsequenced -Wno-unused-comparison -Wno-undefined-var-template -Wno-null-dereference
CC = $(WM_CXX) $(WM_CXXFLAGS) -ftrapping-math CC = $(WM_CXX) -g -std=c++11 $(WM_CXXFLAGS) -ftrapping-math
include $(RULES)/c++$(WM_COMPILE_OPTION) include $(RULES)/c++$(WM_COMPILE_OPTION)

View file

@ -1,2 +1,2 @@
c++DBUG = -ggdb2 -DFULLDEBUG c++DBUG = -ggdb2 -DFULLDEBUG
c++OPT = -O0 -fdefault-inline c++OPT = -O0

View file

@ -1,4 +1,2 @@
c++DBUG = c++DBUG =
c++OPT = -O3 c++OPT = -O1
#c++OPT = -march=nocona -O3
# -ftree-vectorize -ftree-vectorizer-verbose=3

View file

@ -1,2 +1,2 @@
cDBUG = -ggdb -DFULLDEBUG cDBUG = -ggdb -DFULLDEBUG
cOPT = -O1 -fdefault-inline -finline-functions cOPT = -O1 -finline-functions