Merge pull request #362 from mwarusz/lauritzen_et_al_2012

Lauritzen et al 2012 test (do not merge yet)

Author: trontrytel

tests/lauritzen_et_al_2012_GMD/CMakeLists.txt

@@ -0,0 +1,51 @@
+if(APPLE)
+  # needed for the XCode clang to be identified as AppleClang and not Clang
+  cmake_minimum_required(VERSION 3.0) 
+else()
+  # needed for the OpenMP test to work in C++-only project 
+  # (see http://public.kitware.com/Bug/view.php?id=11910)
+  cmake_minimum_required(VERSION 2.8.8) 
+endif()
+
+project(libmpdata++-tests-lauritzen_et_al_2012 CXX)
+
+include(${CMAKE_SOURCE_DIR}/../../libmpdata++-config.cmake)
+if(NOT libmpdataxx_FOUND) 
+  message(FATAL_ERROR "local libmpdata++-config.cmake not found!")
+endif()
+
+if(NOT CMAKE_BUILD_TYPE STREQUAL "Debug")
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${libmpdataxx_CXX_FLAGS_RELEASE}")
+  set(CMAKE_CXX_FLAGS_RELEASE "")
+else()
+  set(CMAKE_CXX_FLAGS_DEBUG ${libmpdataxx_CXX_FLAGS_DEBUG})
+endif()
+
+# to make <libmpdata++/...> work
+set(CMAKE_CXX_FLAGS "-I${CMAKE_CURRENT_SOURCE_DIR}/../.. ${CMAKE_CXX_FLAGS}")
+
+# macro to be used in the subdirectories
+function(libmpdataxx_add_test test)
+  add_executable(${test} ${test}.cpp)
+  target_link_libraries(${test} ${libmpdataxx_LIBRARIES})
+  target_include_directories(${test} PUBLIC ${libmpdataxx_INCLUDE_DIRS})
+  add_test(${test} ${test})
+endfunction()
+
+enable_testing()
+
+libmpdataxx_add_test(gmd_2012)
+
+add_test(gmd_2012_stats_and_plots bash -c "
+  for i in nug_i2_120 nug_i2_240 nug_iga_fct_i2_120 nug_iga_fct_i2_240; do 
+    python  ${CMAKE_CURRENT_SOURCE_DIR}/stats_and_plots.py $i;
+  done
+") 
+
+add_test(gmd_2012_stats_diff bash -c "
+  for i in nug_i2_120 nug_i2_240 nug_iga_fct_i2_120 nug_iga_fct_i2_240; do 
+    echo  ${CMAKE_CURRENT_SOURCE_DIR}/stats_$i.txt.gz;
+    echo  ${CMAKE_CURRENT_BINARY_DIR}/stats_$i.txt;
+    zdiff ${CMAKE_CURRENT_SOURCE_DIR}/refdata/stats_$i.txt.gz ${CMAKE_CURRENT_BINARY_DIR}/stats_$i.txt || exit 1; 
+  done
+") 

tests/lauritzen_et_al_2012_GMD/filaments.py

@@ -0,0 +1,30 @@
+import numpy as np
+import matplotlib
+matplotlib.use('Pdf')
+import matplotlib.pyplot as plt
+
+def calc_filament_diags(g, field_data):
+    taus = np.linspace(0.10, 1.0, 19)
+    lfs = []
+    eps = 1e-12
+    for tau in taus:
+        area_0 = np.sum(g[field_data['cb']['initial'] >= tau - eps])
+        area_h = np.sum(g[field_data['cb']['halfway'] >= tau - eps])
+        if area_0 < eps:
+            lf = 0
+        else:
+            lf = 100 * area_h / area_0
+        lfs.append(lf)
+    return list(zip(taus, lfs))
+
+def plot_filament_diags(dirname, filament_diags, deg):
+    plt.clf()
+    taus, lfs = zip(*filament_diags)
+    plt.plot(taus, lfs, 'ks-', lw = 2, ms = 8)
+    fs = 22
+    plt.xlabel(r'$\tau$', fontsize = fs)
+    plt.ylabel(r'$l_f(\tau, T / 2)$', fontsize = fs)
+    plt.ylim([0, 140])
+    plt.title('filament preservation diagnostics at ${}\\degree$'.format(deg), fontsize = fs)
+    plt.grid()
+    plt.savefig('filaments_' + dirname + '.pdf')

tests/lauritzen_et_al_2012_GMD/gmd_2012.cpp

@@ -0,0 +1,177 @@
+/* 
+ * @file
+ * @copyright University of Warsaw
+ * @section LICENSE
+ * GPLv3+ (see the COPYING file or http://www.gnu.org/licenses/)
+ */
+
+#include <cmath>
+#include <boost/math/constants/constants.hpp>
+
+#include <libmpdata++/solvers/mpdata.hpp>
+#include <libmpdata++/concurr/threads.hpp>
+#include <libmpdata++/output/hdf5_xdmf.hpp>
+
+#include "gmd_2012.hpp"
+
+using namespace libmpdataxx;
+
+struct ct_test_params_t
+{
+  // choosing non-dimensional sphere radius and evolution period
+  static constexpr ::T
+    R = 1,
+    T = 5;
+};
+// to make refering to compile time test params easier
+using tp = ct_test_params_t;
+
+// gaussian hills initial condition
+struct gauss_t
+{
+  T x0, y0;
+  static constexpr T b = 5, hmax = 0.95;
+  T xc(const T x, const T y) const { return tp::R * std::cos(y) * std::cos(x); }
+  T yc(const T x, const T y) const { return tp::R * std::cos(y) * std::sin(x); }
+  T zc(const T x, const T y) const { return tp::R * std::sin(y); }
+  
+  T operator()(const T x, const T y) const
+  {
+    return hmax * std::exp(-b * ( std::pow(xc(x, y) - xc(x0, y0), 2)
+                                + std::pow(yc(x, y) - yc(x0, y0), 2)
+                                + std::pow(zc(x, y) - zc(x0, y0), 2)));
+  }
+  BZ_DECLARE_FUNCTOR2(gauss_t);
+};
+
+// cosine bells initial condition
+struct bells_t
+{
+  T x0, y0;
+  static constexpr T r = tp::R / 2, hmax = 1.0, c = 0.9;
+  T operator()(const T x, const T y) const
+  {
+    T ri = tp::R * std::acos(std::sin(y0) * std::sin(y) + std::cos(y0) * std::cos(y) * std::cos(x - x0));
+    return ri < r ? c * hmax / 2 * (1 + std::cos(pi * ri / r)) : 0;
+  }
+  BZ_DECLARE_FUNCTOR2(bells_t);
+};
+
+// slotted cylinders initial condition
+struct scyls_t
+{
+  T x0, y0;
+  bool flip;
+  static constexpr T r = tp::R / 2, c = 1.0;
+  T operator()(const T x, const T y) const
+  {
+    T ri = tp::R * std::acos(std::sin(y0) * std::sin(y) + std::cos(y0) * std::cos(y) * std::cos(x - x0));
+    bool cond1 = ri < r && std::abs(x - x0) >= r / (6 * tp::R);
+    bool cond2 = flip ? (y - y0) < -5. / 12 * r / tp::R : (y - y0) > 5. / 12 * r / tp::R;
+    bool cond3 = ri < r && (std::abs(x - x0) <  r / (6 * tp::R) && cond2);
+
+    return cond1 || cond3 ? c : 0;
+  }
+  BZ_DECLARE_FUNCTOR2(scyls_t);
+};
+
+template <bool var_dt_arg, int opts_arg, int opts_iters>
+void test(const std::string &base_name, const int ny, const T max_cfl)
+{
+  auto dir_name = base_name + "_" + std::to_string(ny);
+  std::cout << "Calculating: " << dir_name << std::endl;
+
+  struct ct_params_t : ct_params_default_t
+  {
+    using real_t = double;
+    enum { var_dt = var_dt_arg };
+    enum { n_dims = 2 };
+    enum { n_eqns = 4 };
+    enum { opts = opts_arg };
+  };
+  
+
+  const int nx = 2 * ny + 1;
+  const T dx = 2 * pi / (nx - 1), dy = pi / ny;
+  const T time = tp::T;
+  const T dt = dx / (48 * 2 * pi);
+  const int nt = time / dt;
+
+  using slv_out_t = 
+      output::hdf5_xdmf<
+        gmd_2012<ct_params_t, ct_test_params_t>
+    >;
+  typename slv_out_t::rt_params_t p;
+
+  p.n_iters = opts_iters; 
+  p.grid_size = {nx, ny};
+  p.dt = dt;
+  p.di = dx;
+  p.dj = dy;
+  p.max_courant = max_cfl;
+
+  p.outfreq = var_dt_arg ? tp::T / 2 : nt / 2; 
+  p.outvars[0].name = "gh";
+  p.outvars[1].name = "cb";
+  p.outvars[2].name = "sc";
+  p.outvars[3].name = "ccb";
+  p.outdir = dir_name;
+  
+  concurr::threads<
+    slv_out_t, 
+    bcond::cyclic, bcond::cyclic,
+    bcond::polar, bcond::polar
+  > run(p); 
+
+  blitz::firstIndex i;
+  blitz::secondIndex j;
+
+  // coordinates
+  decltype(run.advectee()) X(run.advectee().extent()), Y(run.advectee().extent());
+  X = i * dx;
+  Y = (j + 0.5) * dy - pi / 2;
+
+  std::array<T, 2> x0s = {5 * pi / 6, 7 * pi / 6};
+  std::array<T, 2> y0s = {0, 0};
+
+  // initial conditions
+  const T b = 0.1;
+  run.advectee(0) = 0;
+  run.advectee(1) = b;
+  run.advectee(2) = b;
+
+  for (int m = 0; m < 2; ++m)
+  {
+    run.advectee(0) += gauss_t{x0s[m], y0s[m]}(X, Y);
+    run.advectee(1) += bells_t{x0s[m], y0s[m]}(X, Y);
+    run.advectee(2) += scyls_t{x0s[m], y0s[m], m > 0}(X, Y);
+  }
+  run.advectee(3) = -0.8 * blitz::pow(run.advectee(1), 2) + 0.9;
+  run.g_factor() = tp::R * blitz::cos(Y) * dx * dy;
+  
+  // to avoid divergence check
+  run.advector(0) = 0;
+  run.advector(1) = 0;
+
+  // integration
+  run.advance(var_dt_arg ? time : nt);
+}
+
+int main()
+{
+  const bool var_dt = true;
+  const T max_cfl = 0.90;
+
+  std::vector<int> nys = {120, 240};
+  {
+    enum { opts = opts::nug};
+    const int opts_iters = 2;
+    for (const auto ny : nys) test<var_dt, opts, opts_iters>("nug_i2", ny, max_cfl);
+  }
+  
+  {
+    enum { opts = opts::nug | opts::iga | opts::fct};
+    const int opts_iters = 2;
+    for (const auto ny : nys) test<var_dt, opts, opts_iters>("nug_iga_fct_i2", ny, max_cfl);
+  }
+}

tests/lauritzen_et_al_2012_GMD/gmd_2012.hpp

@@ -0,0 +1,101 @@
+/** 
+ * @file
+ * @copyright University of Warsaw
+ * @section LICENSE
+ * GPLv3+ (see the COPYING file or http://www.gnu.org/licenses/)
+ */
+
+#pragma once
+#include <libmpdata++/solvers/mpdata.hpp>
+
+using T = double;
+constexpr T pi = boost::math::constants::pi<T>();
+
+template <class ct_params_t, class ct_test_params_t>
+class gmd_2012 : public libmpdataxx::solvers::mpdata<ct_params_t>
+{
+  using parent_t = libmpdataxx::solvers::mpdata<ct_params_t>;
+  using parent_t::parent_t;
+  
+  public:
+
+  using real_t = typename ct_params_t::real_t;
+
+  protected:
+
+  using tp = ct_test_params_t;
+
+  libmpdataxx::arrvec_t<typename parent_t::arr_t> &gc_node;
+
+  std::pair<real_t, real_t> gc_exact(const real_t t, const real_t x, const real_t y)
+  {
+    const auto xp = x - 2 * pi * t / tp::T;
+
+    return {
+             (   10 * tp::R / tp::T * std::pow(std::sin(xp), 2) * std::sin(2 * y) * std::cos(pi * t / tp::T)
+               + 2 * pi * tp::R / tp::T * std::cos(y)
+             ) * this->dj * this->dt
+           ,
+             (   10 * tp::R / tp::T * std::sin(2 * xp) * std::pow(std::cos(y), 2) * std::cos(pi * t / tp::T)
+             ) * this->di * this->dt
+           };
+  }
+
+  bool calc_gc() final
+  {
+    using namespace libmpdataxx::arakawa_c;
+    const auto t = this->time + 0.5 * this->dt;
+
+    // calculate exact advectors at grid points
+    for (int i = this->i.first(); i <= this->i.last(); ++i)
+    {
+      for (int j = this->j.first(); j <= this->j.last(); ++j)
+      {
+        const auto x = i * this->di;
+        const auto y = (j+ 0.5) * this->dj - pi / 2;
+
+        auto gc_ij = gc_exact(t, x, y);
+
+        gc_node[0](i, j) = gc_ij.first;
+        gc_node[1](i, j) = gc_ij.second;
+      }
+    }
+
+    this->xchng_sclr(gc_node[0], this->i, this->j);
+    this->xchng_sclr(gc_node[1], this->i, this->j);
+
+    // calculate staggered advectors
+    for (int i = this->i.first()-1; i <= this->i.last(); ++i)
+    {
+      for (int j = this->j.first()-1; j <= this->j.last(); ++j)
+      {
+        this->mem->GC[0](i+h, j) = 0.5 * (gc_node[0](i, j) + gc_node[0](i+1, j));
+        this->mem->GC[1](i, j+h) = 0.5 * (gc_node[1](i, j) + gc_node[1](i, j+1));
+      }
+    }
+
+    auto ex = this->halo - 1;
+    this->xchng_vctr_nrml(this->mem->GC, this->i^ex, this->j^ex);
+
+    return true;
+  }
+
+  public:
+
+  // ctor
+  gmd_2012( 
+    typename parent_t::ctor_args_t args, 
+    const typename parent_t::rt_params_t &p
+  ) :
+    parent_t(args, p),
+    gc_node(args.mem->tmp[__FILE__][0])
+  {}
+
+  static void alloc(
+    typename parent_t::mem_t *mem, 
+    const int &n_iters
+  ) {
+    parent_t::alloc(mem, n_iters);
+    parent_t::alloc_tmp_sclr(mem, __FILE__, 2); // gc_node
+  }
+};