Merge pull request #424 from mwarusz/impl_tht_abs_fix

fix trapezoidal integration of theta absorber in implicit boussinesq solver

.travis.yml

@@ -34,8 +34,6 @@ addons:
       - ubuntu-toolchain-r-test
       - llvm-toolchain-trusty-4.0
 
-osx_image: xcode8
-
 before_install:
     # define a compiler variable for use in if statements
     - if [[ $TRAVIS_OS_NAME == 'linux' && $CXX == 'clang++' ]]; then export COMPILER=clang++; fi
@@ -48,6 +46,12 @@ before_install:
     - if [[ $TRAVIS_OS_NAME == 'linux' ]]; then sudo apt-get update; fi
     - if [[ $TRAVIS_OS_NAME == 'linux' ]]; then export apt_get_install="apt-get install --no-install-recommends -t trusty -y"; fi
     - if [[ $TRAVIS_OS_NAME == 'osx' ]]; then brew update; fi
+
+    # Workaround for upgrading gcc, see https://github.com/travis-ci/travis-ci/issues/8826
+    - if [[ $TRAVIS_OS_NAME == 'osx' ]]; then brew cask uninstall oclint; fi
+
+    # https://github.com/travis-ci/travis-ci/issues/10165
+    - if [[ $TRAVIS_OS_NAME == 'osx' ]]; then rvm get stable; fi
 
     # fixing broken OSX gzip tools
     - if [[ $TRAVIS_OS_NAME == 'osx' ]]; then export PATH=/usr/local/bin:$PATH; fi # so that brew-installed ones are first

.travis_scripts/UWLCM.sh

@@ -1,5 +1,5 @@
 #!/usr/bin/env sh
-set -e
+set -ex
 
 # install libmpata++
 cd libmpdata++/build

.travis_scripts/elliptic_drop.sh

@@ -1,5 +1,5 @@
 #!/usr/bin/env sh
-set -e
+set -ex
 # installing
 cd libmpdata++/build
 sudo make install

.travis_scripts/nair_jablonowski_2008.sh

@@ -1,5 +1,5 @@
 #!/usr/bin/env sh
-set -e
+set -ex
 cd tests/nair_jablonowski_2008
 mkdir build
 cd build

.travis_scripts/paper.sh

@@ -1,5 +1,5 @@
 #!/usr/bin/env sh
-set -e
+set -ex
 cd tests/paper_2015_GMD
 mkdir build
 cd build

.travis_scripts/pbl.sh

@@ -1,5 +1,5 @@
 #!/usr/bin/env sh
-set -e
+set -ex
 cd tests/sandbox
 mkdir build
 cd build

.travis_scripts/sandbox.sh

@@ -1,5 +1,5 @@
 #!/usr/bin/env sh
-set -e
+set -ex
 cd tests/sandbox
 mkdir build
 cd build

.travis_scripts/unit.sh

@@ -1,5 +1,5 @@
 #!/usr/bin/env sh
-set -e
+set -ex
 cd tests/unit
 mkdir build 
 cd build

libmpdata++/solvers/detail/boussinesq_impl.hpp

@@ -102,10 +102,10 @@ namespace libmpdataxx
             }
             case (1):
             {
-              rhs.at(ix::tht)(ijk) += (this->hflux_frc(ijk) - this->tht_abs(ijk) * tht(ijk))
-                                      / (1 + 0.5 * this->dt * this->tht_abs(ijk));
+              rhs.at(ix::tht)(ijk) += this->hflux_frc(ijk);
 
-              rhs.at(ix_w)(ijk) += this->g * (tht(ijk) + 0.5 * this->dt * rhs.at(ix::tht)(ijk)) / this->Tht_ref;
+              rhs.at(ix_w)(ijk) += this->g * (tht(ijk) + 0.5 * this->dt * rhs.at(ix::tht)(ijk))
+                                   / (this->Tht_ref * (1 + 0.5 * this->dt * this->tht_abs(ijk)));
               break;
             }
           }
@@ -116,8 +116,11 @@ namespace libmpdataxx
           parent_t::vip_rhs_impl_fnlz();
 
           const auto &w = this->vips()[ct_params_t::n_dims - 1];
-          this->state(ix::tht)(this->ijk) += - 0.5 * this->dt * w(this->ijk) * this->dtht_e(this->ijk);
-          this->rhs.at(ix::tht)(this->ijk) += -w(this->ijk) * this->dtht_e(this->ijk);
+          this->state(ix::tht)(this->ijk) = ( this->state(ix::tht)(this->ijk) 
+                                            - 0.5 * this->dt * w(this->ijk) * this->dtht_e(this->ijk))
+                                            / (1 + 0.5 * this->dt * this->tht_abs(this->ijk));
+          this->rhs.at(ix::tht)(this->ijk) += -w(this->ijk) * this->dtht_e(this->ijk)
+                                              -this->tht_abs(this->ijk) * this->state(ix::tht)(this->ijk);
         }
 
         public:

libmpdata++/solvers/detail/solver_3d.hpp

@@ -155,6 +155,16 @@ namespace libmpdataxx
           if (!cyclic)
           {
             for (auto &bc : this->bcs[1]) bc->fill_halos_vctr_nrml(arrvec[0], range_ijk[2]^ext^1, range_ijk[0]^ext^h);
+
+            // without this barrier, there is a race condition when some threads handle subdomains
+            // with one gridpoint width, the problem manifests itself, for example, in pbl test
+            // TODO: figure out the exact cause and try to avoid this barrier, what about 2D,
+            //       what about atypical boundary condition choices -- rigid/cyclic/rigid etc
+            if (parent_t::div3_mpdata)
+            {
+              this->mem->barrier();
+            }
+
             for (auto &bc : this->bcs[2]) bc->fill_halos_vctr_nrml(arrvec[0], range_ijk_0__ext_h, range_ijk[1]^ext^1);
 
             for (auto &bc : this->bcs[0]) bc->fill_halos_vctr_nrml(arrvec[1], range_ijk[1]^ext^h, range_ijk[2]^ext^1);

tests/sandbox/pbl/budget.py

@@ -98,7 +98,7 @@
 outname = dirname[4:]
 outfile = open('budget_' + outname + '.txt', 'w')
 for lev_d in zip(Z, hflux, tke_t, dssp, inb):
-    outfile.write('{:6.2f} {:10.3e} {:10.3e} {:10.3e} {:10.1e}\n'.format(*lev_d))
+    outfile.write('{:6.2f} {:10.1e} {:10.1e} {:10.1e} {:10.1e}\n'.format(*lev_d))
 
 fig, axarr = plt.subplots(1, 1, figsize= (10, 6))
 axarr.plot(hflux, Z, 'k-', lw = 1, label = 'B')