Merge pull request #454 from babrodtk/openpm_experiment

Parallel assembly (partial)

CMakeLists.txt

@@ -18,6 +18,7 @@ cmake_minimum_required (VERSION 2.8)
 
 set( OPM_COMMON_ROOT "" CACHE PATH "Root directory containing OPM related cmake modules")
 option(SIBLING_SEARCH "Search for other modules in sibling directories?" ON)
+set( USE_OPENMP_DEFAULT OFF ) # Use of OpenMP is considered experimental
 
 if(NOT OPM_COMMON_ROOT)
    find_package(opm-common QUIET)

examples/flow.cpp

@@ -95,6 +95,10 @@
 #include <boost/filesystem.hpp>
 #include <boost/algorithm/string.hpp>
 
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
 #include <memory>
 #include <algorithm>
 #include <cstdlib>
@@ -150,6 +154,20 @@ try
         std::cout << "**********************************************************************\n\n";
     }
 
+#ifdef _OPENMP
+    if (!getenv("OMP_NUM_THREADS")) {
+        //Default to max(4,#cores) threads,
+        //not number of cores (unless ENV(OMP_NUM_THREADS) is defined)
+        int num_cores = omp_get_num_procs();
+        int num_threads = std::min(4, num_cores);
+        omp_set_num_threads(num_threads);
+    }
+#pragma omp parallel
+    if (omp_get_thread_num() == 0){
+        std::cout << "OpenMP using " << omp_get_num_threads() << " threads." << std::endl;
+    }
+#endif
+
     // Read parameters, see if a deck was specified on the command line.
     if ( output_cout )
     {

opm/autodiff/AutoDiffBlock.hpp

@@ -223,6 +223,7 @@ namespace Opm
                 assert (value().size() == rhs.value().size());
 
                 const int num_blocks = numBlocks();
+#pragma omp parallel for schedule(static)
                 for (int block = 0; block < num_blocks; ++block) {
                     assert(jac_[block].rows() == rhs.jac_[block].rows());
                     assert(jac_[block].cols() == rhs.jac_[block].cols());
@@ -241,6 +242,7 @@ namespace Opm
             if (jac_.empty()) {
                 const int num_blocks = rhs.numBlocks();
                 jac_.resize(num_blocks);
+#pragma omp parallel for schedule(static)
                 for (int block = 0; block < num_blocks; ++block) {
                     jac_[block] = rhs.jac_[block] * (-1.0);
                 }
@@ -249,6 +251,7 @@ namespace Opm
                 assert (value().size() == rhs.value().size());
 
                 const int num_blocks = numBlocks();
+#pragma omp parallel for schedule(static)
                 for (int block = 0; block < num_blocks; ++block) {
                     assert(jac_[block].rows() == rhs.jac_[block].rows());
                     assert(jac_[block].cols() == rhs.jac_[block].cols());
@@ -276,6 +279,7 @@ namespace Opm
             std::vector<M> jac = jac_;
             assert(numBlocks() == rhs.numBlocks());
             int num_blocks = numBlocks();
+#pragma omp parallel for schedule(static)
             for (int block = 0; block < num_blocks; ++block) {
                 assert(jac[block].rows() == rhs.jac_[block].rows());
                 assert(jac[block].cols() == rhs.jac_[block].cols());
@@ -299,6 +303,7 @@ namespace Opm
             std::vector<M> jac = jac_;
             assert(numBlocks() == rhs.numBlocks());
             int num_blocks = numBlocks();
+#pragma omp parallel for schedule(static)
             for (int block = 0; block < num_blocks; ++block) {
                 assert(jac[block].rows() == rhs.jac_[block].rows());
                 assert(jac[block].cols() == rhs.jac_[block].cols());
@@ -324,6 +329,7 @@ namespace Opm
             assert(numBlocks() == rhs.numBlocks());
             M D1(val_.matrix().asDiagonal());
             M D2(rhs.val_.matrix().asDiagonal());
+#pragma omp parallel for schedule(dynamic)
             for (int block = 0; block < num_blocks; ++block) {
                 assert(jac_[block].rows() == rhs.jac_[block].rows());
                 assert(jac_[block].cols() == rhs.jac_[block].cols());
@@ -360,6 +366,7 @@ namespace Opm
             M D1(val_.matrix().asDiagonal());
             M D2(rhs.val_.matrix().asDiagonal());
             M D3((1.0/(rhs.val_*rhs.val_)).matrix().asDiagonal());
+#pragma omp parallel for schedule(dynamic)
             for (int block = 0; block < num_blocks; ++block) {
                 assert(jac_[block].rows() == rhs.jac_[block].rows());
                 assert(jac_[block].cols() == rhs.jac_[block].cols());
@@ -484,6 +491,7 @@ namespace Opm
         int num_blocks = rhs.numBlocks();
         std::vector<typename AutoDiffBlock<Scalar>::M> jac(num_blocks);
         assert(lhs.cols() == rhs.value().rows());
+#pragma omp parallel for schedule(dynamic)
         for (int block = 0; block < num_blocks; ++block) {
             fastSparseProduct(lhs, rhs.derivative()[block], jac[block]);
         }

opm/autodiff/BlackoilModelBase_impl.hpp

@@ -914,6 +914,7 @@ namespace detail {
         trans_all << transi, trans_nnc;
 
         const std::vector<ADB> kr = asImpl().computeRelPerm(state);
+#pragma omp parallel for schedule(static)
         for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
             asImpl().computeMassFlux(phaseIdx, trans_all, kr[canph_[phaseIdx]], state.canonical_phase_pressures[canph_[phaseIdx]], state);
 
@@ -1312,6 +1313,7 @@ namespace detail {
             const int nw = wells.number_of_wells;
             ADB::V retval = ADB::V::Zero(nw);
 
+#pragma omp parallel for schedule(static)
             for (int i=0; i<nw; ++i) {
                 retval[i] = computeHydrostaticCorrection(wells, i, vfp_ref_depth[i], well_perforation_densities, gravity);
             }
@@ -1365,6 +1367,7 @@ namespace detail {
         const int np = wells().number_of_phases;
         const int nw = wells().number_of_wells;
         const Opm::PhaseUsage& pu = fluid_.phaseUsage();
+#pragma omp parallel for schedule(dynamic)
         for (int w = 0; w < nw; ++w) {
             const WellControls* wc = wells().ctrls[w];
             // The current control in the well state overrides
@@ -2075,6 +2078,7 @@ namespace detail {
             // Thp update
             const Opm::PhaseUsage& pu = fluid_.phaseUsage();
             //Loop over all wells
+#pragma omp parallel for schedule(static)
             for (int w=0; w<nw; ++w) {
                 const WellControls* wc = wells().ctrls[w];
                 const int nwc = well_controls_get_num(wc);
@@ -2747,6 +2751,7 @@ namespace detail {
         if (rock_comp_props_ && rock_comp_props_->isActive()) {
             V pm(n);
             V dpm(n);
+#pragma omp parallel for schedule(static)
             for (int i = 0; i < n; ++i) {
                 pm[i] = rock_comp_props_->poroMult(p.value()[i]);
                 dpm[i] = rock_comp_props_->poroMultDeriv(p.value()[i]);
@@ -2754,6 +2759,7 @@ namespace detail {
             ADB::M dpm_diag(dpm.matrix().asDiagonal());
             const int num_blocks = p.numBlocks();
             std::vector<ADB::M> jacs(num_blocks);
+#pragma omp parallel for schedule(dynamic)
             for (int block = 0; block < num_blocks; ++block) {
                 fastSparseProduct(dpm_diag, p.derivative()[block], jacs[block]);
             }
@@ -2775,6 +2781,7 @@ namespace detail {
         if (rock_comp_props_ && rock_comp_props_->isActive()) {
             V tm(n);
             V dtm(n);
+#pragma omp parallel for schedule(static)
             for (int i = 0; i < n; ++i) {
                 tm[i] = rock_comp_props_->transMult(p.value()[i]);
                 dtm[i] = rock_comp_props_->transMultDeriv(p.value()[i]);
@@ -2782,6 +2789,7 @@ namespace detail {
             ADB::M dtm_diag(dtm.matrix().asDiagonal());
             const int num_blocks = p.numBlocks();
             std::vector<ADB::M> jacs(num_blocks);
+#pragma omp parallel for schedule(dynamic)
             for (int block = 0; block < num_blocks; ++block) {
                 fastSparseProduct(dtm_diag, p.derivative()[block], jacs[block]);
             }