Div mpdata dyn (do not merge)

.travis_scripts/sandbox.sh

@@ -10,4 +10,5 @@ VERBOSE=1 $make_j
 OMP_NUM_THREADS=4 ctest -R tgv_2d || cat Testing/Temporary/LastTest.log /
 OMP_NUM_THREADS=4 make -C convergence_2d_3d test || cat convergence_2d_3d/Testing/Temporary/LastTest.log /
 if [[ $TRAVIS_OS_NAME == 'linux' ]]; then OMP_NUM_THREADS=4 make -C convergence_spacetime test || cat convergence_spacetime/Testing/Temporary/LastTest.log /; fi
+OMP_NUM_THREADS=4 make -C bconds_div test || cat bconds_div/Testing/Temporary/LastTest.log /
 cd ../../..

libmpdata++/bcond/cyclic_1d.hpp

@@ -35,6 +35,11 @@ namespace libmpdataxx
       {
 	av[0](this->left_halo_vctr) = av[0](this->rght_intr_vctr);
       }
+
+      void fill_halos_vctr_alng_cyclic(arrvec_t<arr_t> &av, const bool ad = false)
+      {
+        fill_halos_vctr_alng(av, ad);
+      }
     };
 
     template <typename real_t, int halo, bcond_e knd, drctn_e dir, int n_dims, int dim>
@@ -62,6 +67,10 @@ namespace libmpdataxx
 	av[0](this->rght_halo_vctr) = av[0](this->left_intr_vctr);
       }
 
+      void fill_halos_vctr_alng_cyclic(arrvec_t<arr_t> &av, const bool ad = false)
+      {
+        fill_halos_vctr_alng(av, ad);
+      }
     };
   } // namespace bcond
 } // namespace libmpdataxx

libmpdata++/bcond/cyclic_2d.hpp

@@ -68,6 +68,16 @@ namespace libmpdataxx
       {
         fill_halos_sclr(a, j);
       }
+
+      void fill_halos_vctr_alng_cyclic(arrvec_t<arr_t> &av, const rng_t &j, const bool ad = false)
+      {
+        fill_halos_vctr_alng(av, j, ad);
+      }
+
+      void fill_halos_vctr_nrml_cyclic(arr_t &a, const rng_t &j)
+      {
+        fill_halos_vctr_nrml(a, j);
+      }
     };
 
     template <typename real_t, int halo, bcond_e knd, drctn_e dir, int n_dims, int d>
@@ -127,6 +137,16 @@ namespace libmpdataxx
       {
         fill_halos_sclr(a, j);
       }
+
+      void fill_halos_vctr_alng_cyclic(arrvec_t<arr_t> &av, const rng_t &j, const bool ad = false)
+      {
+        fill_halos_vctr_alng(av, j, ad);
+      }
+
+      void fill_halos_vctr_nrml_cyclic(arr_t &a, const rng_t &j)
+      {
+        fill_halos_vctr_nrml(a, j);
+      }
     };
   } // namespace bcond
 } // namespace libmpdataxx

libmpdata++/bcond/cyclic_3d.hpp

@@ -68,6 +68,16 @@ namespace libmpdataxx
       {
         fill_halos_sclr(a, j, k);
       }
+
+      void fill_halos_vctr_alng_cyclic(arrvec_t<arr_t> &av, const rng_t &j, const rng_t &k, const bool ad = false)
+      {
+        fill_halos_vctr_alng(av, j, k, ad);
+      }
+
+      void fill_halos_vctr_nrml_cyclic(arr_t &a, const rng_t &j, const rng_t &k)
+      {
+        fill_halos_vctr_nrml(a, j, k);
+      }
     };
 
     template <typename real_t, int halo, bcond_e knd, drctn_e dir, int n_dims, int d>
@@ -127,6 +137,16 @@ namespace libmpdataxx
       {
         fill_halos_sclr(a, j, k);
       }
+
+      void fill_halos_vctr_alng_cyclic(arrvec_t<arr_t> &av, const rng_t &j, const rng_t &k, const bool ad = false)
+      {
+        fill_halos_vctr_alng(av, j, k, ad);
+      }
+
+      void fill_halos_vctr_nrml_cyclic(arr_t &a, const rng_t &j, const rng_t &k)
+      {
+        fill_halos_vctr_nrml(a, j, k);
+      }
     };
   } // namespace bcond
 } // namespace libmpdataxx

libmpdata++/bcond/detail/bcond_common.hpp

@@ -52,6 +52,9 @@ namespace libmpdataxx
 	{ 
 	  assert(false && "bcond::fill_halos_vctr() called!"); 
 	};
+
+        virtual void fill_halos_vctr_alng_cyclic(arrvec_t<blitz::Array<real_t, 1>> &, const bool ad = false)
+	{};
 
 	// 2D
 	virtual void fill_halos_sclr(arr_2d_t &, const rng_t &, const bool deriv = false) 
@@ -98,6 +101,15 @@ namespace libmpdataxx
 	{
 	  assert(false && "bcond::fill_halos_vctr_nrml() called!");
 	};
+
+        virtual void fill_halos_vctr_alng_cyclic(arrvec_t<blitz::Array<real_t, 2>> &, const rng_t &, const bool ad = false) 
+	{};
+
+	virtual void fill_halos_vctr_nrml_cyclic(blitz::Array<real_t, 2> &, const rng_t &) 
+	{};
+
+        virtual void fill_halos_flux(arrvec_t<blitz::Array<real_t, 2>> &, const rng_t &) 
+	{};
 
 	// 3D
 	virtual void fill_halos_sclr(arr_3d_t &, const rng_t &, const rng_t &, const bool deriv = false) 
@@ -153,6 +165,15 @@ namespace libmpdataxx
 	{
 	  assert(false && "bcond::fill_halos_vctr_nrml() called!");
 	};
+
+        virtual void fill_halos_vctr_alng_cyclic(arrvec_t<blitz::Array<real_t, 3>> &, const rng_t &, const rng_t &, const bool ad = false) 
+	{};
+
+	virtual void fill_halos_vctr_nrml_cyclic(blitz::Array<real_t, 3> &, const rng_t &, const rng_t &) 
+	{};
+
+        virtual void fill_halos_flux(arrvec_t<blitz::Array<real_t, 3>> &, const rng_t &, const rng_t &) 
+	{};
 
 	protected:
 	  // sclr

libmpdata++/bcond/open_1d.hpp

@@ -39,7 +39,7 @@ namespace libmpdataxx
 
       void fill_halos_vctr_alng(arrvec_t<arr_t> &av, const bool ad = false)
       {
-        for (int i = this->left_halo_vctr.first(); i <= this->left_halo_vctr.last(); ++i)
+        for (int i = this->left_halo_vctr.first(); i <= this->left_halo_vctr.last() - (ad ? 1 : 0); ++i)
 	  av[0](rng_t(i, i)) = av[0](this->left_intr_vctr.first());
       }
     };
@@ -72,7 +72,7 @@ namespace libmpdataxx
 
       void fill_halos_vctr_alng(arrvec_t<arr_t> &av, const bool ad = false)
       {
-        for (int i = this->rght_halo_vctr.first(); i <= this->rght_halo_vctr.last(); ++i)
+        for (int i = this->rght_halo_vctr.first() + (ad ? 1 : 0); i <= this->rght_halo_vctr.last(); ++i)
 	  av[0](rng_t(i, i)) = av[0](this->rght_intr_vctr.first());
       }
     };

libmpdata++/bcond/open_2d.hpp

@@ -62,6 +62,12 @@ namespace libmpdataxx
       {
         using namespace idxperm;
         a(pi<d>(this->left_edge_sclr, j)) = sign * edge_velocity(pi<d>(0, j));
+
+        if (halo > 1)
+        {
+          a(pi<d>(this->left_halo_sclr.last() - 1, j)) =   3 * a(pi<d>(this->left_edge_sclr,     j))
+                                                         - 2 * a(pi<d>(this->left_edge_sclr + 1, j));
+        }
       }
 
       void fill_halos_vctr_alng(arrvec_t<arr_t> &av, const rng_t &j, const bool ad = false)
@@ -77,7 +83,7 @@ namespace libmpdataxx
         }
 
 	// zero-divergence condition
-        for (int ii = this->left_halo_vctr.first(); ii <= this->left_halo_vctr.last(); ++ii)
+        for (int ii = this->left_halo_vctr.first(); ii <= this->left_halo_vctr.last() - (ad ? 1 : 0); ++ii)
         {
 	  av[d](pi<d>(ii, j)) = 
 	    av[d](pi<d>(i+h, j)) 
@@ -133,6 +139,12 @@ namespace libmpdataxx
         // equivalent to one-sided derivatives at the boundary
         a(pi<d>(this->rght_halo_sclr.first(), j)) = 2 * a(pi<d>(this->rght_edge_sclr,     j))
                                                       - a(pi<d>(this->rght_edge_sclr - 1, j));
+        if (halo > 1)
+        {
+          a(pi<d>(this->rght_halo_sclr.first() + 1, j)) =   3 * a(pi<d>(this->rght_edge_sclr,     j))
+                                                          - 2 * a(pi<d>(this->rght_edge_sclr - 1, j));
+        }
+
       }
 
       void save_edge_vel(const arr_t &a, const rng_t &j)
@@ -163,7 +175,7 @@ namespace libmpdataxx
         }
 
 	// zero-divergence condition
-        for (int ii = this->rght_halo_vctr.first(); ii <= this->rght_halo_vctr.last(); ++ii)
+        for (int ii = this->rght_halo_vctr.first() + (ad ? 1 : 0); ii <= this->rght_halo_vctr.last(); ++ii)
         {
 	  av[d](pi<d>(ii, j)) = 
 	    av[d](pi<d>(i-h, j)) + (

libmpdata++/bcond/open_3d.hpp

@@ -47,6 +47,11 @@ namespace libmpdataxx
         // equivalent to one-sided derivatives at the boundary
         a(pi<d>(this->left_halo_sclr.last(), j, k)) = 2 * a(pi<d>(this->left_edge_sclr,     j, k))
                                                         - a(pi<d>(this->left_edge_sclr + 1, j, k));
+        if (halo > 1)
+        {
+          a(pi<d>(this->left_halo_sclr.last() - 1, j, k)) =   3 * a(pi<d>(this->left_edge_sclr,     j, k))
+                                                            - 2 * a(pi<d>(this->left_edge_sclr + 1, j, k));
+        }
       }
 
       void save_edge_vel(const arr_t &a, const rng_t &j, const rng_t &k)
@@ -93,7 +98,7 @@ namespace libmpdataxx
 	assert(std::isfinite(sum(av[d+2](pi<d>(i, j, k+h)))));
 
         // zero-divergence condition
-        for (int ii = this->left_halo_vctr.first(); ii <= this->left_halo_vctr.last(); ++ii)
+        for (int ii = this->left_halo_vctr.first(); ii <= this->left_halo_vctr.last() - (ad ? 1 : 0); ++ii)
         {
           av[d](pi<d>(ii, j, k)) = 
             av[d](pi<d>(i+h, j, k)) 
@@ -153,6 +158,12 @@ namespace libmpdataxx
         // equivalent to one-sided derivatives at the boundary
         a(pi<d>(this->rght_halo_sclr.first(), j, k)) = 2 * a(pi<d>(this->rght_edge_sclr,     j, k))
                                                          - a(pi<d>(this->rght_edge_sclr - 1, j, k));
+
+        if (halo > 1)
+        {
+          a(pi<d>(this->rght_halo_sclr.first() + 1, j, k)) =   3 * a(pi<d>(this->rght_edge_sclr,     j, k))
+                                                             - 2 * a(pi<d>(this->rght_edge_sclr - 1, j, k));
+        }
       }
 
       void save_edge_vel(const arr_t &a, const rng_t &j, const rng_t &k)
@@ -197,7 +208,7 @@ namespace libmpdataxx
 	assert(std::isfinite(sum(av[d+2](pi<d>(i, j, k-h)))));
 	assert(std::isfinite(sum(av[d+2](pi<d>(i, j, k+h)))));
 
-        for (int ii = this->rght_halo_vctr.first(); ii <= this->rght_halo_vctr.last(); ++ii)
+        for (int ii = this->rght_halo_vctr.first() + (ad ? 1 : 0); ii <= this->rght_halo_vctr.last(); ++ii)
         {
           av[d](pi<d>(ii, j, k)) = 
             av[d](pi<d>(i-h, j, k)) 

libmpdata++/bcond/rigid_2d.hpp

@@ -35,14 +35,17 @@ namespace libmpdataxx
           a(pi<d>(i, j)) = a(pi<d>(this->left_edge_sclr + n, j));
         }
       }
-      
+
       void fill_halos_pres(arr_t &a, const rng_t &j)
       {
         using namespace idxperm;
         // equivalent to one-sided derivatives at the boundary
-        a(pi<d>(this->left_halo_sclr.last(), j)) = 2 * a(pi<d>(this->left_edge_sclr,     j))
-                                                     - a(pi<d>(this->left_edge_sclr + 1, j));
-      }
+        for (int i = this->left_halo_sclr.first(), n = halo; i <= this->left_halo_sclr.last(); ++i, --n)
+        {
+          a(pi<d>(i, j)) = 2 * a(pi<d>(this->left_edge_sclr,     j))
+                             - a(pi<d>(this->left_edge_sclr + n, j));
+        }
+      }     
 
       void save_edge_vel(const arr_t &, const rng_t &) {}
 
@@ -55,10 +58,20 @@ namespace libmpdataxx
       void fill_halos_vctr_alng(arrvec_t<arr_t> &av, const rng_t &j, const bool ad = false)
       {
 	using namespace idxperm;
-        // zero velocity condition
-        for (int i = this->left_halo_vctr.first(), n = halo; i <= this->left_halo_vctr.last(); ++i, --n)
+        int i = this->left_halo_vctr.last();
+
+        if (!ad)
         {
-	  av[d](pi<d>(i, j)) = -av[d](pi<d>(this->left_edge_sclr + n - h, j));
+          // zero velocity condition
+          av[d](pi<d>(i, j)) = -av[d](pi<d>(this->left_edge_sclr + h, j));
+        }
+        if (halo > 1)
+        {
+          // extrapolation
+          av[d](pi<d>(i - 1, j)) = av[d](pi<d>(this->left_edge_sclr + 1 + h, j)) + 3 * (
+                                    av[d](pi<d>(this->left_edge_sclr - h, j))
+                                  - av[d](pi<d>(this->left_edge_sclr + h, j))
+                                  );
         }
       }
 
@@ -67,6 +80,12 @@ namespace libmpdataxx
         // note intentional sclr
         fill_halos_sclr(a, j);
       }
+
+      void fill_halos_flux(arrvec_t<arr_t> &av, const rng_t &j)
+      {
+        using namespace idxperm;
+	av[d](pi<d>(this->left_halo_vctr.last(), j)) = -av[d](pi<d>(this->left_edge_sclr + h, j));
+      }
     };
 
     template <typename real_t, int halo, bcond_e knd, drctn_e dir, int n_dims, int d>
@@ -98,8 +117,11 @@ namespace libmpdataxx
       {
         using namespace idxperm;
         // equivalent to one-sided derivatives at the boundary
-        a(pi<d>(this->rght_halo_sclr.first(), j)) = 2 * a(pi<d>(this->rght_edge_sclr,     j))
-                                                      - a(pi<d>(this->rght_edge_sclr - 1, j));
+        for (int i = this->rght_halo_sclr.first(), n = 1; i <= this->rght_halo_sclr.last(); ++i, ++n)
+        {
+          a(pi<d>(i, j)) = 2 * a(pi<d>(this->rght_edge_sclr,     j))
+                             - a(pi<d>(this->rght_edge_sclr - n, j));
+        }
       }
 
       void save_edge_vel(const arr_t &, const rng_t &) {}
@@ -113,10 +135,20 @@ namespace libmpdataxx
       void fill_halos_vctr_alng(arrvec_t<arr_t> &av, const rng_t &j, const bool ad = false)
       {
 	using namespace idxperm;
-        // zero velocity condition
-        for (int i = this->rght_halo_vctr.first(), n = 1; i <= this->rght_halo_vctr.last(); ++i, ++n)
+
+        int i = this->rght_halo_vctr.first();
+        if (!ad)
+        {
+          // zero velocity condition
+	  av[d](pi<d>(i, j)) = -av[d](pi<d>(this->rght_edge_sclr - h, j));
+        }
+        if (halo > 1)
         {
-	  av[d](pi<d>(i, j)) = -av[d](pi<d>(this->rght_edge_sclr - n + h, j));
+          // extrapolation
+          av[d](pi<d>(i + 1, j)) = av[d](pi<d>(this->rght_edge_sclr - 1 - h, j)) + 3 * (
+                                    av[d](pi<d>(this->rght_edge_sclr + h, j))
+                                  - av[d](pi<d>(this->rght_edge_sclr - h, j))
+                                  );
         }
       }
 
@@ -125,6 +157,13 @@ namespace libmpdataxx
         // note intentional sclr
         fill_halos_sclr(a, j);
       }
+
+      void fill_halos_flux(arrvec_t<arr_t> &av, const rng_t &j)
+      {
+	using namespace idxperm;
+        // zero flux condition
+	av[d](pi<d>(this->rght_halo_vctr.first(), j)) = -av[d](pi<d>(this->rght_edge_sclr - h, j));
+      }
     };
   } // namespace bcond
 } // namespace libmpdataxx