Merge pull request #45 from WENO-OF/development

Merge development branch for update to version 4.1

.gitignore

@@ -67,7 +67,11 @@ tests/Cases/performanceTest/Case/0/
 tests/Cases/performanceTest/Case/constant/
 tests/Cases/performanceTest/Case/performanceRun.log
 tests/Cases/performanceTest/Case/plotPerformance.dat
+tests/Cases/2DMesh-cyclic/PLOT/results.dat
 
+# ignore test case executable
+tests/src/WENO_TEST
 
 # ignore WENO directories
 tests/**/WENOBase*
+

.gitmodules

@@ -0,0 +1,3 @@
+[submodule "tests/Catch2"]
+	path = tests/Catch2
+	url = https://github.com/catchorg/Catch2.git

Allwmake

@@ -121,7 +121,7 @@ if [[ ${FORCE} != "YES" ]]; then
     checkRepoStatus
 fi
 
-cmake .. ${POSITIONAL}
+cmake .. ${POSITIONAL[@]}
 if [[ ${PARALLEL} == "YES" ]]; then
     make -j ${NUM_CORES}
 else

CI/Dockerfile.OF5

@@ -12,6 +12,9 @@ WORKDIR /home/gitlab/build
 # Bundle lib source
 COPY --chown=gitlab:gitlab . .
 
+# Update git submodules
+RUN git submodule init && git submodule update
+
 # First compile the code
 RUN . $foamDotFile && ./Allwclean && ./Allwmake -f -j
 

CI/Dockerfile.OF7

@@ -12,6 +12,9 @@ WORKDIR /home/gitlab/build
 # Bundle lib source
 COPY --chown=gitlab:gitlab . .
 
+# Update git submodules
+RUN git submodule init && git submodule update
+
 # First compile the code
 RUN . $foamDotFile && ./Allwclean && ./Allwmake -f -j
 

CI/Dockerfile.OF8

@@ -12,6 +12,9 @@ WORKDIR /home/gitlab/build
 # Bundle lib source
 COPY --chown=gitlab:gitlab . .
 
+# Update git submodules
+RUN git submodule init && git submodule update
+
 # First compile the code
 RUN . $foamDotFile && ./Allwclean && ./Allwmake -f -j
 

CI/Dockerfile.OFv1912

@@ -12,6 +12,9 @@ WORKDIR /home/gitlab/build
 # Bundle lib source
 COPY --chown=gitlab:gitlab . .
 
+# Update git submodules
+RUN git submodule init && git submodule update
+
 # First compile the code
 RUN . $foamDotFile && ./Allwclean && ./Allwmake -f -j
 

CI/Dockerfile.OFv2006

@@ -12,6 +12,9 @@ WORKDIR /home/gitlab/build
 # Bundle lib source
 COPY --chown=gitlab:gitlab . .
 
+# Update git submodules
+RUN git submodule init && git submodule update
+
 # First compile the code
 RUN . $foamDotFile && ./Allwclean && ./Allwmake -f -j
 

CI/Dockerfile.OFv2012

@@ -12,6 +12,9 @@ WORKDIR /home/gitlab/build
 # Bundle lib source
 COPY --chown=gitlab:gitlab . .
 
+# Update git submodules
+RUN git submodule init && git submodule update
+
 # First compile the code
 RUN . $foamDotFile && ./Allwclean && ./Allwmake -f -j
 

README.md

@@ -52,24 +52,27 @@ When using this work please cite:
 
 Parallel compilation can be activated with:
 `./Allwmake -j <# of cores>`
-if the number of cores is omitted then the number of cores is determined
-automatically.
+if the number of cores is omitted, it is determined automatically.
 
 On the master branch only clean git commits can be compiled. 
 If a dirty git state shall be compiled the `-f|--force` option has to be used
 on the master branch. 
 
-By default the `march=native` compiler flag is activated. For some cases it is 
-necessary to deactive this flag. To deactivate this flag use:
-`./Allwmake -DMARCH_NATIVE=OFF`
+### CMake Options
 
+Several options can be set over CMake commands with,
+```bash
+./Allwmake -D<CMAKE_OPTION>
+```
+
+|Cmake Command |Options| Description |
+|:-----------------------|:---------|:----------------|
+|MARCH_NATIVE    |ON/OFF | Activates `march=native` flag. Default ON </br> Use this flag if you get an "illegal instruction error" during execution.|
+|USE_LAPACK        |ON/OFF| Use LAPACK library for matrix operations such as eigen values</br>If switched on, check with the [WENO-PerformanceTests](https://github.com/WENO-OF/WENO-PerformanceTests) if the performance improves or decreases.|
+|USE_FMA|ON/OFF|Use std::fma for WENO math functions. Default ON|
+|CMAKE_BUILD_TYPE|Release/Debug/None|When the debug option is selected the OpenFOAM FULLDEBUG flag is activated|
 
-To compute the eigenvalues of a matrix and the Jacobi inverse matrices own
-functions defined in mathFunctionsWENO are used. It is possible to switch on 
-the usage of LAPACK libraries by setting 
-`./Allwmake -DUSE_LAPACK=ON`
-If switched on, check with the [WENO-PerformanceTests](https://github.com/WENO-OF/WENO-PerformanceTests)
-if the performance improves or decreases.
+Commands not listed in the table are forwarded to cmake, allowing to use all standard CMake commands.
 
 ### Note to GNU compiler:
 

libWENOEXT/WENOBase/geometryWENO/mathFunctionsWENO.H

@@ -35,6 +35,9 @@ Author
 #define mathFunctions_H
 
 #include "geometryWENO.H"
+#include "cubicEqn.H"
+#include "mathematicalConstants.H"
+#include "realEigenValues.H"
 #include <math.h>
 
 #ifdef USE_FMA
@@ -56,29 +59,30 @@ namespace Foam
 namespace mathFunctionsWENO
 {
     //- Calculate eigen values of 3x3 matrix
-    //  see: https://en.wikipedia.org/wiki/Eigenvalue_algorithm#3%C3%973_matrices
-    blaze::DynamicVector<double,blaze::columnVector> eigen
+    //  Based on Numerical Recipes 3rd Edition Chapter 11.6
+    blaze::DynamicVector<blaze::complex<double>,blaze::columnVector> eigen
     (
         const geometryWENO::scalarSquareMatrix& A
     );
 
     // Calculate the 3x3 determinant
     double det(const blaze::DynamicMatrix<double>& A);
 
-    // Calculate the 2x2 determinant
-    // solves the system A = a b
-    //                       c d 
-    // with det2 = a*d - c*b
+    //- Calculate the 2x2 determinant
+    //  solves the system A = a b
+    //                        c d 
+    //  with det2 = a*d - c*b
     double det2 (
         const double a, 
         const double d,
         const double c,
         const double b
     );
 
-    // Stablize a 3x3 matrix by calculating the pivoted matrix A' = P*A
+    //- Stablize a 3x3 matrix by calculating the pivoted matrix A' = P*A
     void pivot(geometryWENO::scalarSquareMatrix& A);
 
+    //- Invert a 3x3 matrix
     geometryWENO::scalarSquareMatrix inv(const geometryWENO::scalarSquareMatrix& A);
 
 } // End namespace mathFunctionsWENO 
@@ -90,48 +94,12 @@ namespace mathFunctionsWENO
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-blaze::DynamicVector<double,blaze::columnVector> 
+blaze::DynamicVector<blaze::complex<double>,blaze::columnVector> 
 Foam::mathFunctionsWENO::eigen(const geometryWENO::scalarSquareMatrix& A)
 {
-    // Create blaze vector with three entries
-    blaze::DynamicVector<double,blaze::columnVector> eig(3);
+    realEigenValues eigenVal(A);
 
-    // Create identity matrix
-    geometryWENO::scalarSquareMatrix I;
-    I = blaze::declid(I);
-
-    const double p1 = A(0,1)*A(0,1) + A(0,2)*A(0,2) + A(1,2)*A(1,2);
-    if (p1 == 0)
-    {
-        eig[0] = A(0,0);
-        eig[1] = A(1,1);
-        eig[2] = A(2,2);
-    }
-    else
-    {
-        const double q  = blaze::trace(A)/3.0;
-        const double p2 =   std::pow((A(0,0) - q),2)
-                          + std::pow((A(1,1) - q),2)
-                          + std::pow((A(2,2) - q),2)
-                          + 2.*p1;
-        const double p  = sqrt(p2/6.0);
-        const auto   B  = (A - (q*I))/p;
-        const double r  = 0.5*det(B);
-        // In exact arithmetic for a symmetric matrix  -1 <= r <= 1
-        // but computation error can leave it slightly outside this range
-        double phi;
-        if (r <= -1)
-            phi = M_PI/3.0;
-        else if (r >= 1)
-            phi = 0;
-        else
-            phi = std::acos(r)/3.0;
-
-        eig[0] = q + 2.0*p*std::cos(phi);
-        eig[1] = q + 2.0*p*std::cos(phi + (2.0*M_PI/3.0));
-        eig[2] = 3*q - eig[0] - eig[1];
-    }
-    return eig;
+    return eigenVal.eig();
 }
 
 
@@ -244,21 +212,30 @@ void Foam::mathFunctionsWENO::pivot(geometryWENO::scalarSquareMatrix& A)
 Foam::geometryWENO::scalarSquareMatrix 
 Foam::mathFunctionsWENO::inv(const geometryWENO::scalarSquareMatrix& A)
 {
-    const double invDet = 1.0/det(A);
-
+    // Based on blaze::invertGeneral3x3
     geometryWENO::scalarSquareMatrix AInv;
+
+    AInv(0, 0) = det2(A(1, 1),A(2, 2),A(2, 1),A(1, 2));
+    AInv(1, 0) = det2(A(1, 2),A(2, 0),A(1, 0),A(2, 2));
+    AInv(2, 0) = det2(A(1, 0),A(2, 1),A(2, 0),A(1, 1));
+
+    const double Det = A(0,0)*AInv(0,0) + A(0,1)*AInv(1,0) + A(0,2)*AInv(2,0);
+
+    #ifdef FULLDEBUG
+        if (Det == 0)
+            FatalErrorInFunction 
+                << "Determinant is zero! Division by zero occurs"
+                << exit(FatalError);
+    #endif
+
+    AInv(0, 1) = det2(A(0, 2),A(2, 1),A(0, 1),A(2, 2));
+    AInv(0, 2) = det2(A(0, 1),A(1, 2),A(0, 2),A(1, 1));
+    AInv(1, 1) = det2(A(0, 0),A(2, 2),A(0, 2),A(2, 0));
+    AInv(1, 2) = det2(A(1, 0),A(0, 2),A(0, 0),A(1, 2));
+    AInv(2, 1) = det2(A(2, 0),A(0, 1),A(0, 0),A(2, 1));
+    AInv(2, 2) = det2(A(0, 0),A(1, 1),A(1, 0),A(0, 1));
 
-    AInv(0, 0) = det2(A(1, 1),A(2, 2),A(2, 1),A(1, 2)) * invDet;
-    AInv(0, 1) = det2(A(0, 2),A(2, 1),A(0, 1),A(2, 2)) * invDet;
-    AInv(0, 2) = det2(A(0, 1),A(1, 2),A(0, 2),A(1, 1)) * invDet;
-    AInv(1, 0) = det2(A(1, 2),A(2, 0),A(1, 0),A(2, 2)) * invDet;
-    AInv(1, 1) = det2(A(0, 0),A(2, 2),A(0, 2),A(2, 0)) * invDet;
-    AInv(1, 2) = det2(A(1, 0),A(0, 2),A(0, 0),A(1, 2)) * invDet;
-    AInv(2, 0) = det2(A(1, 0),A(2, 1),A(2, 0),A(1, 1)) * invDet;
-    AInv(2, 1) = det2(A(2, 0),A(0, 1),A(0, 0),A(2, 1)) * invDet;
-    AInv(2, 2) = det2(A(0, 0),A(1, 1),A(1, 0),A(0, 1)) * invDet;
-
-    return AInv;
+    return AInv/Det;
 }