Support arbitrary inflow at open boundary conditions (ft. psc_shock)

src/kg/include/kg/Vec3.h

@@ -38,6 +38,15 @@ struct Vec : gt::sarray<T, N>
     return res;
   }
 
+  // ----------------------------------------------------------------------
+  // unit vector ctors
+  KG_INLINE static Vec unit(int d)
+  {
+    Vec res;
+    res[d] = value_type(1);
+    return res;
+  }
+
   // ----------------------------------------------------------------------
   // converting to Vec of different type (e.g., float -> double)
 
@@ -277,6 +286,15 @@ struct Vec : gt::sarray<T, N>
     return res;
   }
 
+  KG_INLINE Vec reverse() const
+  {
+    Vec res;
+    for (int i = 0; i < N; i++) {
+      res[i] = (*this)[N - 1 - i];
+    }
+    return res;
+  }
+
   // conversion to pointer
 
   KG_INLINE operator const T*() const { return this->data(); }

src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx

@@ -3,6 +3,7 @@
 #include "kg/VecRange.hxx"
 #include "fields.hxx"
 #include "bnd_fields.hxx"
+#include "radiating_bnd.hxx"
 
 #include <mrc_bits.h>
 
@@ -84,6 +85,10 @@ struct BndFields_ : BndFieldsBase
     for (int p = 0; p < mflds.n_patches(); p++) {
       // lo
       for (int d = 0; d < 3; d++) {
+        if (grid.bc.fld_lo[d] == BND_FLD_OPEN && radiation) {
+          radiation->update_cache_lower(grid.time(), d);
+        }
+
         if (grid.atBoundaryLo(p, d)) {
           switch (grid.bc.fld_lo[d]) {
             case BND_FLD_PERIODIC: {
@@ -105,6 +110,10 @@ struct BndFields_ : BndFieldsBase
       }
       // hi
       for (int d = 0; d < 3; d++) {
+        if (grid.bc.fld_hi[d] == BND_FLD_OPEN && radiation) {
+          radiation->update_cache_upper(grid.time(), d);
+        }
+
         if (grid.atBoundaryHi(p, d)) {
           switch (grid.bc.fld_hi[d]) {
             case BND_FLD_PERIODIC: {
@@ -543,29 +552,33 @@ struct BndFields_ : BndFieldsBase
     stop[d0] = 0;
 
     for (Int3 i3 : VecRange(start, stop)) {
-      Int3 edge_idx = i3;
-      edge_idx[d0] = 0;
-
-      Int3 neg1_d1{0, 0, 0};
-      neg1_d1[d1] = -1;
-
-      Int3 neg1_d2{0, 0, 0};
-      neg1_d2[d2] = -1;
+      Int3 edge_idx = i3 + Int3::unit(d0);
+
+      real_t s = 0.0;
+      real_t p = 0.0;
+      if (radiation) {
+        Real3 x3_s = (Real3(edge_idx) + Real3::unit(d1) * real_t(0.5)) *
+                     Real3(grid.domain.dx);
+        Real3 x3_p = (Real3(edge_idx) + Real3::unit(d2) * real_t(0.5)) *
+                     Real3(grid.domain.dx);
+        s = radiation->pulse_s_lower(grid.time(), d0, p, x3_s);
+        p = radiation->pulse_p_lower(grid.time(), d0, p, x3_p);
+      }
 
-      F(H2, i3) = (/* + 4.f * C_s_pulse_y1(x,y,z+0.5*dz,t), where d0=y */
-                   -2.f * (F(E1, edge_idx) - background_e[d1]) -
-                   dtdx[d2] * (F(H0, edge_idx) - F(H0, edge_idx + neg1_d2)) -
-                   (1.f - dtdx[d0]) * (F(H2, edge_idx) - background_h[d2]) +
-                   dt * F(J1, edge_idx)) /
-                    (1.f + dtdx[d0]) +
-                  background_h[d2];
-      F(H1, i3) = (/* + 4.f * C_p_pulse_y1(x+.5*dx,y,z,t), where d0=y */
-                   +2.f * (F(E2, edge_idx) - background_e[d2]) -
-                   dtdx[d1] * (F(H0, edge_idx) - F(H0, edge_idx + neg1_d1)) -
-                   (1.f - dtdx[d0]) * (F(H1, edge_idx) - background_h[d1]) +
-                   dt * F(J2, edge_idx)) /
-                    (1.f + dtdx[d0]) +
-                  background_h[d1];
+      F(H2, i3) =
+        (4.f * s - 2.f * (F(E1, edge_idx) - background_e[d1]) -
+         dtdx[d2] * (F(H0, edge_idx) - F(H0, edge_idx - Int3::unit(d2))) -
+         (1.f - dtdx[d0]) * (F(H2, edge_idx) - background_h[d2]) +
+         dt * F(J1, edge_idx)) /
+          (1.f + dtdx[d0]) +
+        background_h[d2];
+      F(H1, i3) =
+        (-4.f * p + 2.f * (F(E2, edge_idx) - background_e[d2]) -
+         dtdx[d1] * (F(H0, edge_idx) - F(H0, edge_idx - Int3::unit(d1))) -
+         (1.f - dtdx[d0]) * (F(H1, edge_idx) - background_h[d1]) +
+         dt * F(J2, edge_idx)) /
+          (1.f + dtdx[d0]) +
+        background_h[d1];
     }
   }
 
@@ -590,25 +603,27 @@ struct BndFields_ : BndFieldsBase
     stop[d0] = grid.ldims[d0] + 1;
 
     for (Int3 i3 : VecRange(start, stop)) {
-      Int3 edge_idx = i3;
-      edge_idx[d0] -= 1;
-
-      Int3 neg1_d1{0, 0, 0};
-      neg1_d1[d1] = -1;
-
-      Int3 neg1_d2{0, 0, 0};
-      neg1_d2[d2] = -1;
+      Int3 edge_idx = i3 - Int3::unit(d0);
+
+      real_t s = 0.0;
+      real_t p = 0.0;
+      if (radiation) {
+        Real3 x3_s = (Real3(edge_idx) + Real3::unit(d1) * real_t(0.5)) *
+                     Real3(grid.domain.dx);
+        Real3 x3_p = (Real3(edge_idx) + Real3::unit(d2) * real_t(0.5)) *
+                     Real3(grid.domain.dx);
+        s = radiation->pulse_s_upper(grid.time(), d0, p, x3_s);
+        p = radiation->pulse_p_upper(grid.time(), d0, p, x3_p);
+      }
 
-      F(H2, i3) = (/* + 4.f * C_s_pulse_y2(x,y,z+0.5*dz,t), where d0=y */
-                   +2.f * (F(E1, i3) - background_e[d1]) +
-                   dtdx[d2] * (F(H0, i3) - F(H0, i3 + neg1_d2)) -
+      F(H2, i3) = (-4.f * s + 2.f * (F(E1, i3) - background_e[d1]) +
+                   dtdx[d2] * (F(H0, i3) - F(H0, i3 - Int3::unit(d2))) -
                    (1.f - dtdx[d0]) * (F(H2, edge_idx) - background_h[d2]) -
                    dt * F(J1, i3)) /
                     (1.f + dtdx[d0]) +
                   background_h[d2];
-      F(H1, i3) = (/* + 4.f * C_p_pulse_y2(x+.5*dx,y,z,t), where d0=y */
-                   -2.f * (F(E2, i3) - background_e[d2]) +
-                   dtdx[d1] * (F(H0, i3) - F(H0, i3 + neg1_d1)) -
+      F(H1, i3) = (4.f * p - 2.f * (F(E2, i3) - background_e[d2]) +
+                   dtdx[d1] * (F(H0, i3) - F(H0, i3 - Int3::unit(d1))) -
                    (1.f - dtdx[d0]) * (F(H1, edge_idx) - background_h[d1]) -
                    dt * F(J2, i3)) /
                     (1.f + dtdx[d0]) +
@@ -618,6 +633,8 @@ struct BndFields_ : BndFieldsBase
 
   Vec3<real_t> background_e = {0.0, 0.0, 0.0};
   Vec3<real_t> background_h = {0.0, 0.0, 0.0};
+
+  RadiatingBoundary<real_t>* radiation = nullptr;
 };
 
 // ======================================================================

src/libpsc/psc_bnd_fields/radiating_bnd.hxx

@@ -0,0 +1,17 @@
+#include "kg/Vec3.h"
+
+template <typename real_t>
+struct RadiatingBoundary
+{
+  using Real3 = Vec3<real_t>;
+
+  virtual real_t pulse_s_lower(double t, int d, int p, Real3 x3) = 0;
+  virtual real_t pulse_p_lower(double t, int d, int p, Real3 x3) = 0;
+
+  virtual real_t pulse_s_upper(double t, int d, int p, Real3 x3) = 0;
+  virtual real_t pulse_p_upper(double t, int d, int p, Real3 x3) = 0;
+
+  // FIXME these are a hack for a specific subclass to work
+  virtual void update_cache_lower(double t, int d) {}
+  virtual void update_cache_upper(double t, int d) {}
+};

src/psc_shock.cxx

@@ -32,6 +32,8 @@ using Collision = PscConfig::Collision;
 using Checks = PscConfig::Checks;
 using Marder = PscConfig::Marder;
 using OutputParticles = PscConfig::OutputParticles;
+using real_t = PscConfig::Mfields::real_t;
+using Real3 = Vec3<real_t>;
 
 // ======================================================================
 // Global parameters
@@ -54,6 +56,9 @@ double b_x;
 double b_y;
 double b_z;
 
+Real3 background_e;
+Real3 background_h;
+
 int nx;
 int ny;
 int nz;
@@ -112,6 +117,12 @@ void setupParameters(int argc, char** argv)
   b_y = b_mag * cos(b_angle_y_to_x_rad);
   b_z = 0.0;
 
+  Real3 v_upstream = {v_upstream_x, v_upstream_y, v_upstream_z};
+  double gamma = 1 / sqrt(1 - v_upstream.mag2());
+  background_e = -gamma * v_upstream.cross({b_x, b_y, b_z});
+  // note: this only holds for vx=vz=0
+  background_h = {b_x * gamma, b_y, b_z * gamma};
+
   nx = inputParams.get<int>("nx");
   ny = inputParams.get<int>("ny");
   nz = inputParams.get<int>("nz");
@@ -662,6 +673,145 @@ void initializeFields(MfieldsState& mflds)
   boost_fields(mflds, -v_upstream_y);
 }
 
+struct AdvectedPeriodicFields : RadiatingBoundary<real_t>
+{
+  static const int DIM_Y = 1;
+
+  AdvectedPeriodicFields(MfieldsState& mflds, real_t v_advect,
+                         Real3 background_e, Real3 background_h)
+    : v_advect(v_advect), grid(mflds.grid())
+  {
+    // FIXME would be better to exclude J, but the interpolator uses EX, etc.
+    auto&& e_b_fields = mflds.storage().view(_all, _all, _all, _all, _all);
+    // FIXME this probably isn't the best way to copy a gtensor array
+    cycled_fields = gt::zeros_like(e_b_fields);
+    cycled_fields.view(_all, _all, _all, _all, _all) = e_b_fields;
+
+    for (int d = 0; d < 3; d++) {
+      cycled_fields.view(_all, _all, _all, EX + d, _all) =
+        cycled_fields.view(_all, _all, _all, EX + d, _all) - background_e[d];
+      cycled_fields.view(_all, _all, _all, HX + d, _all) =
+        cycled_fields.view(_all, _all, _all, HX + d, _all) - background_h[d];
+    }
+  }
+
+  Real3 advect_x3(Real3 x3, double t)
+  {
+    return x3 - Real3::unit(DIM_Y) * real_t(t * v_advect);
+  }
+
+  int shift_to_patch_local(Real3& x3_advected)
+  {
+    real_t patch_size = grid.domain.length[DIM_Y] / grid.domain.np[DIM_Y];
+    int n_patches_to_the_left = 0;
+    while (x3_advected[DIM_Y] < grid.domain.corner[DIM_Y]) {
+      x3_advected[DIM_Y] += patch_size;
+      n_patches_to_the_left += 1;
+    }
+    return n_patches_to_the_left;
+  }
+
+  void calc_e_h(double t, int p, Real3 x3, int d_e, real_t& e, int d_h,
+                real_t& h)
+  {
+    Real3 x3_advected = advect_x3(x3, t);
+    int n_patches_to_the_left = shift_to_patch_local(x3_advected);
+    assert(n_patches_to_the_left == n_patch_cycles);
+
+    ip.set_coeffs(x3_advected * grid.domain.dx_inv);
+    auto em = decltype(ip)::fields_t(
+      cycled_fields.view(_all, _all, _all, _all, p), -grid.ibn);
+
+    switch (d_e) {
+      case 0: e = ip.ex(em); break;
+      case 1: e = ip.ey(em); break;
+      case 2: e = ip.ez(em); break;
+    }
+    switch (d_h) {
+      case 0: h = ip.hx(em); break;
+      case 1: h = ip.hy(em); break;
+      case 2: h = ip.hz(em); break;
+    }
+  }
+
+  real_t pulse_s_lower(double t, int d, int p, Real3 x3) override
+  {
+    int d1 = (d + 1) % 3;
+    int d2 = (d + 2) % 3;
+
+    real_t e, h;
+    calc_e_h(t, p, x3, d1, e, d2, h);
+
+    return (e + h) / 2.0;
+  }
+
+  real_t pulse_p_lower(double t, int d, int p, Real3 x3) override
+  {
+    int d1 = (d + 1) % 3;
+    int d2 = (d + 2) % 3;
+
+    real_t e, h;
+    calc_e_h(t, p, x3, d2, e, d1, h);
+
+    return (e - h) / 2.0;
+  }
+
+  real_t pulse_s_upper(double t, int d, int p, Real3 x3) override
+  {
+    return 0.0;
+  }
+
+  real_t pulse_p_upper(double t, int d, int p, Real3 x3) override
+  {
+    return 0.0;
+  }
+
+  void update_cache_lower(double t, int d) override
+  {
+    if (d != DIM_Y) {
+      return;
+    }
+
+    // TODO make this work with >1 patch per process?
+    assert(grid.n_patches() == 1);
+
+    Real3 x3_advected = advect_x3({0.0, 0.0, 0.0}, t);
+    int n_patches_to_the_left = shift_to_patch_local(x3_advected);
+
+    if (n_patches_to_the_left == n_patch_cycles) {
+      // "cache hit"; no need to cycle
+      return;
+    }
+
+    LOG_INFO("cycling turbulence...\n");
+
+    // hack: guess the rank based on how mrc does it for simple domains
+    // (can't use mrc, because it wouldn't apply periodicity)
+    Int3 np = grid.domain.np;
+    Int3 proc = grid.localPatchInfo(0).idx3;
+    Int3 dest_proc = (proc + Int3::unit(DIM_Y)) % np;
+    int dest_rank = flatten_index(dest_proc.reverse(), np.reverse());
+    Int3 source_proc = (proc - Int3::unit(DIM_Y) + np) % np;
+    int source_rank = flatten_index(source_proc.reverse(), np.reverse());
+
+    MPI_Status status;
+    int err = MPI_Sendrecv_replace(cycled_fields.data(), cycled_fields.size(),
+                                   MpiDtypeTraits<real_t>::value(), dest_rank,
+                                   0, source_rank, 0, grid.comm(), &status);
+
+    n_patch_cycles += 1;
+  }
+
+  const Grid_t& grid;
+  gt::gtensor<real_t, 5> cycled_fields;
+  int n_patch_cycles = 0;
+  real_t v_advect;
+  InterpolateEM<
+    Fields3d<decltype(cycled_fields.view(_all, _all, _all, _all, 0)), Dim>,
+    opt_ip_1st_ec, Dim>
+    ip;
+};
+
 // ======================================================================
 // run
 
@@ -760,10 +910,10 @@ static void run(int argc, char** argv)
 
   psc.add_gauss_corrector(&marder);
 
-  double gamma = 1 / sqrt(1 - sqr(v_upstream_y));
-  psc.bndf.background_e =
-    -gamma * Double3{0, v_upstream_y, 0}.cross({b_x, b_y, b_z});
-  psc.bndf.background_h = {b_x * gamma, b_y, b_z * gamma};
+  psc.bndf.background_e = background_e;
+  psc.bndf.background_h = background_h;
+  psc.bndf.radiation =
+    new AdvectedPeriodicFields{mflds, 1.0, background_e, background_h};
 
   psc.add_diagnostic(&outf);
   psc.add_diagnostic(&outp);