Merge pull request #120 from kaschau/copyRotate

Copy rotate

examples/couetteFlow.py

@@ -41,41 +41,24 @@ def simulate():
         mbDims=[1, 1, 1],
         dimsPerBlock=[2, ny, 2],
         lengths=[0.01, h, 0.01],
+        periodic=[True, False, False],
     )
 
     mb.initSolverArrays(config)
 
     blk = mb[0]
 
     # face 1
-    blk.getFace(1).bcFam = None
-    blk.getFace(1).bcType = "b1"
-    blk.getFace(1).neighbor = 0
-    blk.getFace(1).orientation = "123"
     blk.getFace(1).commRank = 0
     # face 2
-    blk.getFace(2).bcFam = None
-    blk.getFace(2).bcType = "b1"
-    blk.getFace(2).neighbor = 0
-    blk.getFace(2).orientation = "123"
     blk.getFace(2).commRank = 0
     # face 3
-    blk.getFace(3).bcFam = None
     blk.getFace(3).bcType = "adiabaticNoSlipWall"
-    blk.getFace(3).neighbor = None
-    blk.getFace(3).orientation = None
-    blk.getFace(3).commRank = 0
     # face 4 isoT moving wall
-    blk.getFace(4).bcFam = "whoosh"
     blk.getFace(4).bcType = "adiabaticMovingWall"
-    blk.getFace(4).neighbor = None
-    blk.getFace(4).orientation = None
-    blk.getFace(4).commRank = 0
 
     for face in [5, 6]:
-        blk.getFace(face).bcFam = None
         blk.getFace(face).bcType = "adiabaticSlipWall"
-        blk.getFace(face).commRank = 0
 
     valueDict = {"u": wallSpeed, "v": 0.0, "w": 0.0}
     face4 = blk.getFace(4)

examples/oneDAdvection.py

@@ -28,21 +28,15 @@ def simulate():
         mbDims=[1, 1, 1],
         dimsPerBlock=[41, 2, 2],
         lengths=[1, 0.01, 0.01],
+        periodic=[True, False, False],
     )
     mb.initSolverArrays(config)
 
     blk = mb[0]
-    for face in blk.faces:
+    for face in blk.faces[2::]:
         face.bcType = "adiabaticSlipWall"
 
-    blk.getFace(1).bcType = "b1"
-    blk.getFace(1).neighbor = 0
-    blk.getFace(1).orientation = "123"
     blk.getFace(1).commRank = 0
-
-    blk.getFace(2).bcType = "b1"
-    blk.getFace(2).neighbor = 0
-    blk.getFace(2).orientation = "123"
     blk.getFace(2).commRank = 0
 
     mb.setBlockCommunication()

examples/threeDTaylorGreen.py

@@ -33,6 +33,7 @@ def simulate():
         mbDims=[1, 1, 1],
         dimsPerBlock=[NE, NN, NX],
         lengths=[2 * np.pi for _ in range(3)],
+        periodic=[True, True, True],
     )
 
     mb.initSolverArrays(config)
@@ -41,13 +42,7 @@ def simulate():
     ng = blk.ng
 
     for face in blk.faces:
-        face.bcType = "b1"
-        face.neighbor = 0
-        face.orientation = "123"
         face.commRank = 0
-    pg.writers.writeGrid(mb)
-    pg.writers.writeConnectivity(mb)
-    mb.generateHalo()
 
     mb.setBlockCommunication()
 

examples/threeDTaylorGreenSkewed.py

@@ -39,6 +39,7 @@ def simulate():
         mbDims=[1, 1, 1],
         dimsPerBlock=[NE, NN, NX],
         lengths=[2 * np.pi for _ in range(3)],
+        periodic=[True, True, True],
     )
 
     mb.initSolverArrays(config)
@@ -85,10 +86,6 @@ def simulate():
                 )
 
     for face in blk.faces:
-        face.bcType = "b1"
-        face.bcFam = None
-        face.neighbor = 0
-        face.orientation = "123"
         face.commRank = 0
 
     mb.setBlockCommunication()

examples/twoDEulerVortex.py

@@ -31,7 +31,8 @@ def simulate():
         mb,
         mbDims=[1, 1, 1],
         dimsPerBlock=[NE, NN, 2],
-        lengths=[1, 1, 0.01],
+        lengths=[12, 12, 0.01],
+        periodic=[True, True, False],
     )
 
     blk = mb[0]
@@ -59,19 +60,11 @@ def simulate():
                 + Ay * np.sin(2 * np.pi * kappa) * np.sin(lamY * np.pi * E * delX / Lx)
             )
 
-    for face in blk.faces:
-        face.bcType = "b1"
-        face.bcFam = None
-        face.neighbor = 0
-        face.orientation = "123"
+    for face in blk.faces[0:4]:
         face.commRank = 0
     for f in [5, 6]:
         face = blk.getFace(f)
         face.bcType = "adiabaticSlipWall"
-        face.bcFam = None
-        face.neighbor = None
-        face.orientation = "000"
-        face.commRank = None
 
     mb.initSolverArrays(config)
     mb.setBlockCommunication()

src/bindings.cpp

@@ -169,6 +169,8 @@ PYBIND11_MODULE(compute, m) {
     .def_readwrite("recvBuffer4", &face_::recvBuffer4)
     .def_readwrite("tempRecvBuffer3", &face_::tempRecvBuffer3)
     .def_readwrite("tempRecvBuffer4", &face_::tempRecvBuffer4)
+    .def_readwrite("periodicRotMatrixUp", &face_::periodicRotMatrixUp)
+    .def_readwrite("periodicRotMatrixDown", &face_::periodicRotMatrixDown)
 
     .def_readwrite("intervalDt", &face_::intervalDt)
     .def_readwrite("cubicSplineAlphas", &face_::cubicSplineAlphas)

src/compute/boundaryConditions/bindingsBoundaryConditions.cpp

@@ -101,4 +101,21 @@ void bindBoundaryConditions(py::module_ &m) {
             py::arg("thtrdat_ object"),
             py::arg("terms"),
             py::arg("tme"));
+
+  //  |----> periodics.cpp
+  py::module periodics = bcs.def_submodule("periodics", "periodics boundary conditions module");
+  periodics.def("periodicRotHigh", &periodicRotHigh, "High rotational periodic face",
+            py::arg("block_ object"),
+            py::arg("face_ object"),
+            py::arg("eos pointer"),
+            py::arg("thtrdat_ object"),
+            py::arg("terms"),
+            py::arg("tme"));
+  periodics.def("periodicRotLow", &periodicRotLow, "Low rotational periodic face",
+            py::arg("block_ object"),
+            py::arg("face_ object"),
+            py::arg("eos pointer"),
+            py::arg("thtrdat_ object"),
+            py::arg("terms"),
+            py::arg("tme"));
 }

src/compute/boundaryConditions/exits.cpp

@@ -91,7 +91,7 @@ void constantPressureSubsonicExit(block_ b,
       threeDsubview dqdz0 = getHaloSlice(b.dqdz, face._nface, s0);
 
       Kokkos::parallel_for(
-          "Constant pressure subsonic exit euler terms", range_face,
+          "Constant pressure subsonic exit viscous terms", range_face,
           KOKKOS_LAMBDA(const int i, const int j, const int l) {
             // neumann all gradients
             dqdx0(i, j, l) = dqdx1(i, j, l);

src/compute/boundaryConditions/periodics.cpp

@@ -0,0 +1,165 @@
+#include "Kokkos_Core.hpp"
+#include "block_.hpp"
+#include "compute.hpp"
+#include "face_.hpp"
+#include "kokkos_types.hpp"
+#include "thtrdat_.hpp"
+
+void periodicRotHigh(block_ b,
+                     face_ face,
+                     const std::function<void(block_, thtrdat_, int, std::string)> &/*eos*/,
+                     const thtrdat_ /*th*/,
+                     const std::string terms,
+                     const double /*tme*/) {
+  //-------------------------------------------------------------------------------------------|
+  // Apply BC to face, slice by slice.
+  //-------------------------------------------------------------------------------------------|
+  const int ng = b.ng;
+  int s0, s1, s2, plus;
+  setHaloSlices(s0, s1, s2, plus, b.ni, b.nj, b.nk, ng, face._nface);
+
+  if (terms.compare("euler") == 0) {
+
+    threeDsubview q1 = getHaloSlice(b.q, face._nface, s1);
+    MDRange2 range_face = MDRange2({0, 0}, {q1.extent(0), q1.extent(1)});
+    for (int g = 0; g < b.ng; g++) {
+      s0 -= plus * g;
+
+      threeDsubview q0 = getHaloSlice(b.q, face._nface, s0);
+      threeDsubview Q0 = getHaloSlice(b.Q, face._nface, s0);
+
+      Kokkos::parallel_for(
+          "Rotate periodic euler terms", range_face,
+          KOKKOS_LAMBDA(const int i, const int j) {
+            // rotate velocity vector up to high face
+            double tempU, tempV, tempW;
+            double u = q0(i, j, 1);
+            double v = q0(i, j, 2);
+            double w = q0(i, j, 3);
+            tempU = face.periodicRotMatrixUp(0,0)*u + face.periodicRotMatrixUp(0,1)*v + face.periodicRotMatrixUp(0,2)*w;
+            tempV = face.periodicRotMatrixUp(1,0)*u + face.periodicRotMatrixUp(1,1)*v + face.periodicRotMatrixUp(1,2)*w;
+            tempW = face.periodicRotMatrixUp(2,0)*u + face.periodicRotMatrixUp(2,1)*v + face.periodicRotMatrixUp(2,2)*w;
+
+            // Update velocity
+            q0(i, j, 1) = tempU;
+            q0(i, j, 2) = tempV;
+            q0(i, j, 3) = tempW;
+
+            // Update momentum
+            Q0(i, j, 1) = tempU*Q0(i,j,0);
+            Q0(i, j, 2) = tempV*Q0(i,j,0);
+            Q0(i, j, 3) = tempW*Q0(i,j,0);
+          });
+    }
+
+  } else if (terms.compare("viscous") == 0) {
+
+    threeDsubview dqdx1 = getHaloSlice(b.dqdx, face._nface, s1);
+    MDRange3 range_face =
+        MDRange3({0, 0, 0}, {static_cast<long>(dqdx1.extent(0)),
+                             static_cast<long>(dqdx1.extent(1)), b.ne});
+    for (int g = 0; g < b.ng; g++) {
+      s0 -= plus * g;
+      threeDsubview dqdx0 = getHaloSlice(b.dqdx, face._nface, s0);
+      threeDsubview dqdy0 = getHaloSlice(b.dqdy, face._nface, s0);
+      threeDsubview dqdz0 = getHaloSlice(b.dqdz, face._nface, s0);
+
+      Kokkos::parallel_for(
+          "Periodic viscous terms", range_face,
+          KOKKOS_LAMBDA(const int i, const int j, const int l) {
+            // rotate gradients vector up to high face
+            double tempdx, tempdy, tempdz;
+            double dx = dqdx0(i, j, l);
+            double dy = dqdy0(i, j, l);
+            double dz = dqdz0(i, j, l);
+            tempdx = face.periodicRotMatrixUp(0,0)*dx + face.periodicRotMatrixUp(0,1)*dy + face.periodicRotMatrixUp(0,2)*dz;
+            tempdy = face.periodicRotMatrixUp(1,0)*dx + face.periodicRotMatrixUp(1,1)*dy + face.periodicRotMatrixUp(1,2)*dz;
+            tempdz = face.periodicRotMatrixUp(2,0)*dx + face.periodicRotMatrixUp(2,1)*dy + face.periodicRotMatrixUp(2,2)*dz;
+
+            dqdx0(i, j, l) = tempdx;
+            dqdy0(i, j, l) = tempdy;
+            dqdz0(i, j, l) = tempdz;
+          });
+    }
+  }
+}
+
+
+void periodicRotLow(block_ b,
+                     face_ face,
+                     const std::function<void(block_, thtrdat_, int, std::string)> &/*eos*/,
+                     const thtrdat_ /*th*/,
+                     const std::string terms,
+                     const double /*tme*/) {
+  //-------------------------------------------------------------------------------------------|
+  // Apply BC to face, slice by slice.
+  //-------------------------------------------------------------------------------------------|
+  const int ng = b.ng;
+  int s0, s1, s2, plus;
+  setHaloSlices(s0, s1, s2, plus, b.ni, b.nj, b.nk, ng, face._nface);
+
+  if (terms.compare("euler") == 0) {
+
+    threeDsubview q1 = getHaloSlice(b.q, face._nface, s1);
+    MDRange2 range_face = MDRange2({0, 0}, {q1.extent(0), q1.extent(1)});
+    for (int g = 0; g < b.ng; g++) {
+      s0 -= plus * g;
+
+      threeDsubview q0 = getHaloSlice(b.q, face._nface, s0);
+      threeDsubview Q0 = getHaloSlice(b.Q, face._nface, s0);
+
+      Kokkos::parallel_for(
+          "Rotate periodic euler terms", range_face,
+          KOKKOS_LAMBDA(const int i, const int j) {
+            // rotate velocity vector up to high face
+            double tempU, tempV, tempW;
+            double u = q0(i, j, 1);
+            double v = q0(i, j, 2);
+            double w = q0(i, j, 3);
+            tempU = face.periodicRotMatrixDown(0,0)*u + face.periodicRotMatrixDown(0,1)*v + face.periodicRotMatrixDown(0,2)*w;
+            tempV = face.periodicRotMatrixDown(1,0)*u + face.periodicRotMatrixDown(1,1)*v + face.periodicRotMatrixDown(1,2)*w;
+            tempW = face.periodicRotMatrixDown(2,0)*u + face.periodicRotMatrixDown(2,1)*v + face.periodicRotMatrixDown(2,2)*w;
+
+            // Update velocity
+            q0(i, j, 1) = tempU;
+            q0(i, j, 2) = tempV;
+            q0(i, j, 3) = tempW;
+
+            // Update momentum
+            Q0(i, j, 1) = tempU*Q0(i,j,0);
+            Q0(i, j, 2) = tempV*Q0(i,j,0);
+            Q0(i, j, 3) = tempW*Q0(i,j,0);
+          });
+    }
+
+  } else if (terms.compare("viscous") == 0) {
+
+    threeDsubview dqdx1 = getHaloSlice(b.dqdx, face._nface, s1);
+    MDRange3 range_face =
+        MDRange3({0, 0, 0}, {static_cast<long>(dqdx1.extent(0)),
+                             static_cast<long>(dqdx1.extent(1)), b.ne});
+    for (int g = 0; g < b.ng; g++) {
+      s0 -= plus * g;
+      threeDsubview dqdx0 = getHaloSlice(b.dqdx, face._nface, s0);
+      threeDsubview dqdy0 = getHaloSlice(b.dqdy, face._nface, s0);
+      threeDsubview dqdz0 = getHaloSlice(b.dqdz, face._nface, s0);
+
+      Kokkos::parallel_for(
+          "Periodic viscous terms", range_face,
+          KOKKOS_LAMBDA(const int i, const int j, const int l) {
+            // rotate gradients vector up to high face
+            double tempdx, tempdy, tempdz;
+            double dx = dqdx0(i, j, l);
+            double dy = dqdy0(i, j, l);
+            double dz = dqdz0(i, j, l);
+            tempdx = face.periodicRotMatrixDown(0,0)*dx + face.periodicRotMatrixDown(0,1)*dy + face.periodicRotMatrixDown(0,2)*dz;
+            tempdy = face.periodicRotMatrixDown(1,0)*dx + face.periodicRotMatrixDown(1,1)*dy + face.periodicRotMatrixDown(1,2)*dz;
+            tempdz = face.periodicRotMatrixDown(2,0)*dx + face.periodicRotMatrixDown(2,1)*dy + face.periodicRotMatrixDown(2,2)*dz;
+
+            dqdx0(i, j, l) = tempdx;
+            dqdy0(i, j, l) = tempdy;
+            dqdz0(i, j, l) = tempdz;
+          });
+    }
+  }
+}

src/compute/compute.hpp

@@ -89,6 +89,15 @@ void supersonicExit(
     block_ b, face_ face,
     const std::function<void(block_, thtrdat_, int, std::string)> &eos,
     const thtrdat_ th, const std::string terms, const double tme);
+//    |------> periodics
+void periodicRotHigh(
+    block_ b, face_ face,
+    const std::function<void(block_, thtrdat_, int, std::string)> &eos,
+    const thtrdat_ th, const std::string terms, const double tme);
+void periodicRotLow(
+    block_ b, face_ face,
+    const std::function<void(block_, thtrdat_, int, std::string)> &eos,
+    const thtrdat_ th, const std::string terms, const double tme);
 
 // ./chemistry
 //    |------> CH4_O2_Stanford_Skeletal

src/compute/face_.hpp

@@ -25,6 +25,9 @@ struct face_ {
   double intervalDt;
   int currentInterval = -1;
 
+  // Periodic rotation matricies
+  twoDview periodicRotMatrixUp, periodicRotMatrixDown;
+
 };
 
 #endif

src/peregrinepy/bcs/bcfam_templates/walls/periodicRot.yaml

@@ -0,0 +1,7 @@
+---
+
+periodic:
+  bcType: periodicRot
+  bcVals:
+    periodicSpan: 45  # degrees
+    periodicAxis: [1, 0, 0]

src/peregrinepy/bcs/bcfam_templates/walls/periodicTrans.yaml

@@ -0,0 +1,7 @@
+---
+
+periodic:
+  bcType: periodicTrans
+  bcVals:
+    periodicSpan: 0.1
+    periodicAxis: [1, 0, 0]

src/peregrinepy/consistify.py

@@ -27,7 +27,6 @@ def consistify(mb):
         # Update spatial derivatives
         mb.dqdxyz(blk)
 
-        # TODO: can we get rid of this if check?
         if mb.config["RHS"]["diffusion"]:
             # Apply viscous boundary conditions
             for face in blk.faces: