From 70ed6ba9677091df89a0b68b66088f1aa4a27b3d Mon Sep 17 00:00:00 2001
From: Henrik Rusche <h.rusche@wikki-gmbh.de>
Date: Tue, 19 Feb 2019 11:06:50 +0100
Subject: [PATCH 1/5] BUGFIX: View factor smoothing now takes into account
 reciprocity condition

---
 .../radiationModel/viewFactor/viewFactor.C    | 247 +++++++++++++++++-
 1 file changed, 239 insertions(+), 8 deletions(-)

diff --git a/src/thermophysicalModels/radiation/radiationModel/viewFactor/viewFactor.C b/src/thermophysicalModels/radiation/radiationModel/viewFactor/viewFactor.C
index 60a8ba652..3c81fc7c1 100644
--- a/src/thermophysicalModels/radiation/radiationModel/viewFactor/viewFactor.C
+++ b/src/thermophysicalModels/radiation/radiationModel/viewFactor/viewFactor.C
@@ -64,7 +64,7 @@ void Foam::radiation::viewFactor::initialise()
         }
     }
 
-    selectedPatches_.resize(count--);
+    selectedPatches_.resize(count);
 
     if (debug)
     {
@@ -173,6 +173,82 @@ void Foam::radiation::viewFactor::initialise()
 
     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())
     {
         Fmatrix_.reset
@@ -200,18 +276,173 @@ void Foam::radiation::viewFactor::initialise()
         {
             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;
-                for (label j=0; j<totalNCoarseFaces_; j++)
+                scalar maxErr = -GREAT;
+                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;
             }
         }
 

From a7f9e38d5cff40d77ecac0f89ea1050975887aa8 Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Wed, 1 May 2019 09:14:36 +0100
Subject: [PATCH 2/5] Change reconstruct zero option name

---
 tutorials/mesh/cfMesh/cartesianMesh/sawOctree/Allrun | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tutorials/mesh/cfMesh/cartesianMesh/sawOctree/Allrun b/tutorials/mesh/cfMesh/cartesianMesh/sawOctree/Allrun
index 11b1fcaf8..e41175eac 100755
--- a/tutorials/mesh/cfMesh/cartesianMesh/sawOctree/Allrun
+++ b/tutorials/mesh/cfMesh/cartesianMesh/sawOctree/Allrun
@@ -9,5 +9,5 @@ if [ "$WM_PROJECT" = "OpenFOAM" ]
 then
 runApplication reconstructParMesh -constant -fullMatch
 else
-runApplication reconstructParMesh -zeroTime
+runApplication reconstructParMesh -withZero
 fi

From 076919ae5575906db9cc39a4f5a8b58a74bfc538 Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Wed, 1 May 2019 09:40:29 +0100
Subject: [PATCH 3/5] Bugfix: call to completeAssembly in wrong place.  Henrik
 Rusche

---
 .../coupledFvMatrices/coupledFvMatrix/coupledFvMatrix.C       | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

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

From b71851619f653d31d63326e46fdd1747d4eed14b Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Wed, 1 May 2019 13:55:18 +0200
Subject: [PATCH 4/5] Porting: Clean up of initialisation of static data,
 casting and comparison with enum, Mac OSX

---
 src/OSspecific/POSIX/signals/sigFpe.C                     | 4 ++--
 src/foam/containers/NamedEnum/NamedEnum.C                 | 8 +++++++-
 .../chemistryReaders/chemkinReader/chemkinReader.C        | 2 +-
 3 files changed, 10 insertions(+), 4 deletions(-)

diff --git a/src/OSspecific/POSIX/signals/sigFpe.C b/src/OSspecific/POSIX/signals/sigFpe.C
index 6b48544ec..d849cb74c 100644
--- a/src/OSspecific/POSIX/signals/sigFpe.C
+++ b/src/OSspecific/POSIX/signals/sigFpe.C
@@ -326,7 +326,7 @@ void Foam::sigFpe::set(const bool verbose)
         {
             vm_protect(
                 mach_task_self(),
-                (uintptr_t)zone,
+                reinterpret_cast<uintptr_t>(zone),
                 sizeof(malloc_zone_t),
                 0,
                 VM_PROT_READ | VM_PROT_WRITE
@@ -338,7 +338,7 @@ void Foam::sigFpe::set(const bool verbose)
         {
             vm_protect(
                 mach_task_self(),
-                (uintptr_t)zone,
+                reinterpret_cast<uintptr_t>(zone),
                 sizeof(malloc_zone_t),
                 0,
                 VM_PROT_READ
diff --git a/src/foam/containers/NamedEnum/NamedEnum.C b/src/foam/containers/NamedEnum/NamedEnum.C
index d5421cea2..776162ef9 100644
--- a/src/foam/containers/NamedEnum/NamedEnum.C
+++ b/src/foam/containers/NamedEnum/NamedEnum.C
@@ -23,9 +23,15 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "NamedEnum.H"
 #include "stringList.H"
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+
+template<class Enum, int nEnum>
+const char* Foam::NamedEnum<Enum, nEnum>::names[nEnum];
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class Enum, int nEnum>
diff --git a/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.C b/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.C
index b69fe5aaa..bdaa6471d 100644
--- a/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.C
+++ b/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.C
@@ -392,7 +392,7 @@ void Foam::chemkinReader::addPressureDependentReaction
 
         default:
         {
-            if (fofType < 4)
+            if (fofType < unknownFallOffFunctionType)
             {
                 FatalErrorIn("chemkinReader::addPressureDependentReaction")
                     << "Fall-off function type "

From 99a5e3ecc80fe3d318914f8367b02ff9ad3043d0 Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Wed, 1 May 2019 13:55:56 +0200
Subject: [PATCH 5/5] Darwin compilation settings; optimiser bug requires lower
 compiler opt level

---
 wmake/rules/darwinIntel64Clang/c++    | 4 ++--
 wmake/rules/darwinIntel64Clang/c++Opt | 5 +----
 wmake/rules/darwinIntel64Gcc/c++      | 4 ++--
 wmake/rules/darwinIntel64Gcc/c++Debug | 2 +-
 wmake/rules/darwinIntel64Gcc/c++Opt   | 4 +---
 wmake/rules/darwinIntel64Gcc/cDebug   | 2 +-
 6 files changed, 8 insertions(+), 13 deletions(-)

diff --git a/wmake/rules/darwinIntel64Clang/c++ b/wmake/rules/darwinIntel64Clang/c++
index 0e55fe6a1..afd204ae3 100644
--- a/wmake/rules/darwinIntel64Clang/c++
+++ b/wmake/rules/darwinIntel64Clang/c++
@@ -1,8 +1,8 @@
 .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
 
 include $(RULES)/c++$(WM_COMPILE_OPTION)
diff --git a/wmake/rules/darwinIntel64Clang/c++Opt b/wmake/rules/darwinIntel64Clang/c++Opt
index 120b1006d..5319335a0 100644
--- a/wmake/rules/darwinIntel64Clang/c++Opt
+++ b/wmake/rules/darwinIntel64Clang/c++Opt
@@ -1,5 +1,2 @@
 c++DBUG     =
-# c++OPT      = -O4
-c++OPT      = -O3
-#c++OPT      = -march=nocona -O3
-# -ftree-vectorize -ftree-vectorizer-verbose=3
+c++OPT      = -O1 -ftree-vectorize
diff --git a/wmake/rules/darwinIntel64Gcc/c++ b/wmake/rules/darwinIntel64Gcc/c++
index f4ad3fbfa..021e26b60 100644
--- a/wmake/rules/darwinIntel64Gcc/c++
+++ b/wmake/rules/darwinIntel64Gcc/c++
@@ -1,8 +1,8 @@
 .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)
 
diff --git a/wmake/rules/darwinIntel64Gcc/c++Debug b/wmake/rules/darwinIntel64Gcc/c++Debug
index 684957e11..83d3e4224 100644
--- a/wmake/rules/darwinIntel64Gcc/c++Debug
+++ b/wmake/rules/darwinIntel64Gcc/c++Debug
@@ -1,2 +1,2 @@
 c++DBUG    = -ggdb2 -DFULLDEBUG
-c++OPT      = -O0 -fdefault-inline
+c++OPT      = -O0
diff --git a/wmake/rules/darwinIntel64Gcc/c++Opt b/wmake/rules/darwinIntel64Gcc/c++Opt
index 3446f7f58..a54f7458b 100644
--- a/wmake/rules/darwinIntel64Gcc/c++Opt
+++ b/wmake/rules/darwinIntel64Gcc/c++Opt
@@ -1,4 +1,2 @@
 c++DBUG     =
-c++OPT      = -O3
-#c++OPT      = -march=nocona -O3
-# -ftree-vectorize -ftree-vectorizer-verbose=3
+c++OPT      = -O1
diff --git a/wmake/rules/darwinIntel64Gcc/cDebug b/wmake/rules/darwinIntel64Gcc/cDebug
index 72b638f45..6c71ed93e 100644
--- a/wmake/rules/darwinIntel64Gcc/cDebug
+++ b/wmake/rules/darwinIntel64Gcc/cDebug
@@ -1,2 +1,2 @@
 cDBUG       = -ggdb -DFULLDEBUG
-cOPT        = -O1 -fdefault-inline -finline-functions
+cOPT        = -O1 -finline-functions