_DEUTLICH_ mehr Schwuppdizität!!

Interface zwischen Python und C++ verbessert, so dass weniger kopiert und
transponiert wird. Scheint ne Menge zu bringen :-)

fresnel/solver.py

@@ -136,7 +136,9 @@ class Solver3d():
         u[:, 0] = boundary[2]
         u[:,-1] = boundary[3]
 
-        self.u = _solver.step_ACF(up, rxp, ryp, rcp, u, rx, ry, rc)
+        # the C++ functions use Eigen3 internally, which requires C-Major data storage
+        u_T = _solver.step_ACF(up.T, rxp.T, ryp.T, rcp.T, u.T, rx.T, ry.T, rc.T)
+        self.u = u_T.T
 
         return self.u
 

source/finite_difference.h

@@ -110,14 +110,14 @@ namespace finite_differences{
     
 
   // Two-step alternating-direction Peaceman-Rachford scheme (J. W. homas: Numerical Partial Differential Equations)
-  const array_2D step_ACF(const Eigen::Ref<const array_2D> up, 
-                          const Eigen::Ref<const array_2D> rxp,
-                          const Eigen::Ref<const array_2D> ryp, 
-                          const Eigen::Ref<const array_2D> rcp, 
-                          const Eigen::Ref<const array_2D> u, 
-                          const Eigen::Ref<const array_2D> rx,
-                          const Eigen::Ref<const array_2D> ry, 
-                          const Eigen::Ref<const array_2D> rc) {
+  const array_2D step_ACF(const Eigen::Ref<const array_2D> &up, 
+                          const Eigen::Ref<const array_2D> &rxp,
+                          const Eigen::Ref<const array_2D> &ryp, 
+                          const Eigen::Ref<const array_2D> &rcp, 
+                          const Eigen::Ref<const array_2D> &u, 
+                          const Eigen::Ref<const array_2D> &rx,
+                          const Eigen::Ref<const array_2D> &ry, 
+                          const Eigen::Ref<const array_2D> &rc) {
                           
     // rx and ry differ by a factor of 1/2 from the definitions in Fuhse et al. (2005)
  
@@ -125,23 +125,19 @@ namespace finite_differences{
     // .  :  i+1 
     // TODO: interpolate for half-step?
       
-    // TODO: ausrichtung überprüfen!
-      
     // coordinate system:
     // numpy (z,y,x), eigen (rows, cols): x <-> cols, y <-> rows
     
-    size_t nx = u.cols() - 2;
-    size_t ny = u.rows() - 2;
+    size_t nx = u.rows() - 2;
+    size_t ny = u.cols() - 2;
       
-    //std::cout << "nx " << u.cols() << ", ny " << u.rows() << std::endl; 
+    std::cout << "nx " << u.cols() << ", ny " << u.rows() << std::endl; 
     
     // first halfstep
     // x at i + 1/2
     // y at i
     
-    // note: internal storage is column-major such that column-wise access would be much faster
-      
-    array_2D uhalf_t = array_2D::Zero(u.cols(), u.rows()); // transposed shape
+    array_2D uhalf = array_2D::Zero(u.rows(), u.cols());
     
     #pragma omp parallel for
     for (size_t iy = 1; iy <= ny; ++iy) {
@@ -150,26 +146,27 @@ namespace finite_differences{
       array_1D r = array_1D::Zero(nx,1);
         
       for (size_t ix = 1; ix <= nx; ++ ix) {
-        m(0,ix-1) = - rx(iy, ix - 1); 
-        m(1,ix-1) = 1. + 2. * rx(iy, ix) - rc(iy, ix);
-        m(2,ix-1) = - rx(iy, ix); 
+        m(0,ix-1) = - rx(ix - 1,iy); 
+        m(1,ix-1) = 1. + 2. * rx(ix, iy) - rc(ix, iy);
+        m(2,ix-1) = - rx(ix, iy); 
         
-        r(ix-1) = (1. - 2. * ryp(iy, ix) + rcp(iy, ix)) * up(iy, ix)
-                + ryp(iy, ix) * (up(iy - 1, ix) + up(iy + 1, ix));
+        r(ix-1) = (1. - 2. * ryp(ix, iy) + rcp(ix, iy)) * up(ix, iy)
+                + ryp(ix, iy) * (up(ix, iy - 1) + up(ix, iy + 1));
       }
     
       // boundary conditions
-      r(0) += rx(iy, 0) * u(iy, 0);
-      r(nx-1) += rx(iy, nx+1) * u(iy, nx+1);
+      // TODO: use interpolated values?
+      r(0) += rx(0, iy) * u(0, iy);
+      r(nx-1) += rx(nx+1, iy) * u(nx+1, iy);
         
-      uhalf_t.col(iy).segment(1,nx) = algebra::tridiagonal(m, r);
+      uhalf.col(iy).segment(1,nx) = algebra::tridiagonal(m, r);
     }
       
-    array_2D uhalf {uhalf_t.transpose()};
     
     // second halfstep
       
-    array_2D u_new { u };
+    // note: internal storage is column-major such that column-wise access would be much faster
+    array_2D u_new_trans = array_2D::Zero(u.cols(), u.rows()); //transposed shape
 
     #pragma omp parallel for
     for (size_t ix = 1; ix <= nx; ++ix) {
@@ -178,20 +175,22 @@ namespace finite_differences{
       array_1D r = array_1D::Zero(ny,1);
         
       for (size_t iy = 1; iy <= ny; ++ iy) {
-         m(0,iy-1) = -ry(iy-1, ix);
-         m(1,iy-1) = 1. + 2. * ry(iy, ix) - rc(iy, ix);
-         m(2,iy-1) = -ry(iy, ix);
+         m(0,iy-1) = -ry(ix, iy-1);
+         m(1,iy-1) = 1. + 2. * ry(ix, iy) - rc(ix, iy);
+         m(2,iy-1) = -ry(ix, iy);
           
-         r(iy-1) = ( 1. - 2. * rxp(iy, ix) + rcp(iy, ix) ) * uhalf(iy, ix)
-                 + rxp(iy, ix) * ( uhalf(iy, ix-1) + uhalf(iy, ix+1) ); 
+         r(iy-1) = ( 1. - 2. * rxp(ix, iy) + rcp(ix, iy) ) * uhalf(ix, iy)
+                 + rxp(ix, iy) * ( uhalf(ix-1,iy) + uhalf(ix+1, iy) ); 
       }
         
-      r(0) += ry(0, ix) * u(0, ix);
-      r(ny-1) += ry(ny+1, ix) * u(ny+1, ix);
+      r(0) += ry(ix, 0) * u(ix, 0);
+      r(ny-1) += ry(ix, ny+1) * u(ix, ny+1);
 
-      u_new.col(ix).segment(1,ny) = algebra::tridiagonal(m, r);
+      u_new_trans.col(ix).segment(1,ny) = algebra::tridiagonal(m, r);
     }
     
+    array_2D u_new {u_new_trans.transpose()};
+      
     return u_new;
   }
 

source/python.cpp

@@ -11,11 +11,25 @@ using namespace finite_differences;
 PYBIND11_MODULE(_solver, m){
   
   m.def("step_AF",
-        &finite_differences::step_AF
+        &finite_differences::step_AF,
+        py::arg("up").noconvert(),
+        py::arg("rxp").noconvert(),
+        py::arg("rcp").noconvert(),
+        py::arg("u").noconvert(),
+        py::arg("rx").noconvert(),
+        py::arg("rc").noconvert()
   );
   
   m.def("step_ACF",
-        &finite_differences::step_ACF
+        &finite_differences::step_ACF,
+        py::arg("up").noconvert(),
+        py::arg("rxp").noconvert(),
+        py::arg("ryp").noconvert(),
+        py::arg("rcp").noconvert(),
+        py::arg("u").noconvert(),
+        py::arg("rx").noconvert(),
+        py::arg("ry").noconvert(),
+        py::arg("rc").noconvert()
   );
  
 }