Add CONDUCTOR boundary

cmake/kokkosConfig.cmake

@@ -37,41 +37,6 @@ set(Kokkos_ENABLE_OPENMP
     ${default_KOKKOS_ENABLE_OPENMP}
     CACHE BOOL "Enable OpenMP")
 
-# set memory space
-if(${Kokkos_ENABLE_CUDA})
-  add_compile_definitions(CUDA_ENABLED)
-  set(ACC_MEM_SPACE Kokkos::CudaSpace)
-elseif(${Kokkos_ENABLE_HIP})
-  add_compile_definitions(HIP_ENABLED)
-  set(ACC_MEM_SPACE Kokkos::HIPSpace)
-else()
-  set(ACC_MEM_SPACE Kokkos::HostSpace)
-endif()
-
-set(HOST_MEM_SPACE Kokkos::HostSpace)
-
-# set execution space
-if(${Kokkos_ENABLE_CUDA})
-  set(ACC_EXE_SPACE Kokkos::Cuda)
-elseif(${Kokkos_ENABLE_HIP})
-  set(ACC_EXE_SPACE Kokkos::HIP)
-elseif(${Kokkos_ENABLE_OPENMP})
-  set(ACC_EXE_SPACE Kokkos::OpenMP)
-else()
-  set(ACC_EXE_SPACE Kokkos::Serial)
-endif()
-
-if(${Kokkos_ENABLE_OPENMP})
-  set(HOST_EXE_SPACE Kokkos::OpenMP)
-else()
-  set(HOST_EXE_SPACE Kokkos::Serial)
-endif()
-
-add_compile_options("-D AccelExeSpace=${ACC_EXE_SPACE}")
-add_compile_options("-D AccelMemSpace=${ACC_MEM_SPACE}")
-add_compile_options("-D HostExeSpace=${HOST_EXE_SPACE}")
-add_compile_options("-D HostMemSpace=${HOST_MEM_SPACE}")
-
 if(${BUILD_TESTING} STREQUAL "OFF")
   set(Kokkos_ENABLE_TESTS
       OFF

input.example.toml

@@ -90,7 +90,7 @@
     # Boundary conditions for fields:
     #   @required
     #   @type: 1/2/3-size array of string tuples, each of size 1 or 2
-    #   @valid: "PERIODIC", "MATCH", "FIXED", "ATMOSPHERE", "CUSTOM", "HORIZON"
+    #   @valid: "PERIODIC", "MATCH", "FIXED", "ATMOSPHERE", "CUSTOM", "HORIZON", "CONDUCTOR"
     #   @example: [["CUSTOM", "MATCH"]] (for 2D spherical [[rmin, rmax]])
     #   @note: When periodic in any of the directions, you should only set one value: [..., ["PERIODIC"], ...]
     #   @note: In spherical, bondaries in theta/phi are set automatically (only specify bc @ [rmin, rmax]): [["ATMOSPHERE", "MATCH"]]

legacy/tests/kernels-sr.cpp

@@ -1,17 +1,17 @@
-#include "wrapper.h"
-
 #include <Kokkos_Core.hpp>
 
 #include <iostream>
 #include <stdexcept>
 
+#include "wrapper.h"
+
 #include METRIC_HEADER
 #include PGEN_HEADER
 
-#include "particle_macros.h"
-
 #include "kernels/particle_pusher_sr.hpp"
 
+#include "particle_macros.h"
+
 template <typename T>
 void put_value(ntt::array_t<T*>& arr, T value, int i) {
   auto arr_h = Kokkos::create_mirror_view(arr);
@@ -221,4 +221,4 @@ auto main(int argc, char* argv[]) -> int {
   ntt::GlobalFinalize();
 
   return 0;
-}
+}

setups/srpic/shock/pgen.hpp

@@ -79,9 +79,13 @@ namespace user {
     using arch::ProblemGenerator<S, M>::C;
     using arch::ProblemGenerator<S, M>::params;
 
-    const real_t drift_ux, temperature;
-
-    const real_t  Btheta, Bphi, Bmag;
+    // gas properties
+    const real_t drift_ux, temperature, filling_fraction;
+    // injector properties
+    const real_t injector_velocity, injection_start, dt;
+    const int   injection_frequency;
+    // magnetic field properties
+    real_t        Btheta, Bphi, Bmag;
     InitFields<D> init_flds;
 
     inline PGen(const SimulationParams& p, const Metadomain<S, M>& m)
@@ -91,40 +95,249 @@ namespace user {
       , Bmag { p.template get<real_t>("setup.Bmag", ZERO) }
       , Btheta { p.template get<real_t>("setup.Btheta", ZERO) }
       , Bphi { p.template get<real_t>("setup.Bphi", ZERO) }
-      , init_flds { Bmag, Btheta, Bphi, drift_ux } {}
+      , init_flds { Bmag, Btheta, Bphi, drift_ux }
+      , filling_fraction { p.template get<real_t>("setup.filling_fraction", 1.0) }
+      , injector_velocity { p.template get<real_t>("setup.injector_velocity", 1.0) }
+      , injection_start { p.template get<real_t>("setup.injection_start", 0.0) }
+      , injection_frequency { p.template get<int>("setup.injection_frequency", 100) }
+      , dt { p.template get<real_t>("algorithms.timestep.dt") } {}
 
     inline PGen() {}
 
-    auto FixFieldsConst(const bc_in&, const em& comp) const
-      -> std::pair<real_t, bool> {
-      if (comp == em::ex2) {
-        return { init_flds.ex2({ ZERO }), true };
+    auto MatchFields(real_t time) const -> InitFields<D> {
+      return init_flds;
+    }
+
+    auto ResetFields(const em& comp) const -> real_t {
+      if (comp == em::ex1) {
+        return init_flds.ex1({ ZERO });
+      } else if (comp == em::ex2) {
+        return init_flds.ex2({ ZERO });
       } else if (comp == em::ex3) {
-        return { init_flds.ex3({ ZERO }), true };
+        return init_flds.ex3({ ZERO });
+      } else if (comp == em::bx1) {
+        return init_flds.bx1({ ZERO });
+      } else if (comp == em::bx2) {
+        return init_flds.bx2({ ZERO });
+      } else if (comp == em::bx3) {
+        return init_flds.bx3({ ZERO });
       } else {
-        return { ZERO, false };
+        raise::Error("Invalid component", HERE);
+        return ZERO;
       }
     }
 
-    auto MatchFields(real_t time) const -> InitFields<D> {
-      return init_flds;
-    }
-
     inline void InitPrtls(Domain<S, M>& local_domain) {
+
+      // minimum and maximum position of particles
+      real_t xg_min = local_domain.mesh.extent(in::x1).first;
+      real_t xg_max = local_domain.mesh.extent(in::x1).first +
+                      filling_fraction * (local_domain.mesh.extent(in::x1).second -
+                                          local_domain.mesh.extent(in::x1).first);
+
+      // define box to inject into
+      boundaries_t<real_t> box;
+      // loop over all dimensions
+      for (unsigned short d { 0 }; d < static_cast<unsigned short>(M::Dim); ++d) {
+        // compute the range for the x-direction
+        if (d == static_cast<unsigned short>(in::x1)) {
+          box.push_back({ xg_min, xg_max });
+        } else {
+          // inject into full range in other directions
+          box.push_back(Range::All);
+        }
+      }
+
+      // spatial distribution of the particles
+      // -> hack to use the uniform distribution in NonUniformInjector
+      const auto spatial_dist = arch::Piston<S, M>(local_domain.mesh.metric,
+                                                   xg_min,
+                                                   xg_max,
+                                                   in::x1);
+
+      // energy distribution of the particles
       const auto energy_dist = arch::Maxwellian<S, M>(local_domain.mesh.metric,
                                                       local_domain.random_pool,
                                                       temperature,
                                                       -drift_ux,
                                                       in::x1);
 
-      const auto injector = arch::UniformInjector<S, M, arch::Maxwellian>(
+      const auto injector = arch::NonUniformInjector<S, M, arch::Maxwellian, arch::Piston>(
         energy_dist,
+        spatial_dist,
         { 1, 2 });
-      arch::InjectUniform<S, M, arch::UniformInjector<S, M, arch::Maxwellian>>(
+
+      arch::InjectNonUniform<S, M, arch::NonUniformInjector<S, M, arch::Maxwellian, arch::Piston>>(
         params,
         local_domain,
         injector,
-        1.0);
+        1.0,   // target density
+        false, // no weights
+        box);
+    }
+
+    void CustomPostStep(std::size_t step, long double time, Domain<S, M>& domain) {
+
+      // check if the injector should be active
+      if (step % injection_frequency != 0) {
+        return;
+      }
+
+      // initial position of injector
+      const auto x_init = domain.mesh.extent(in::x1).first +
+                          filling_fraction * (domain.mesh.extent(in::x1).second -
+                                              domain.mesh.extent(in::x1).first);
+
+      // check if injector is supposed to start moving already
+      const auto dt_inj = time - injection_start > ZERO ? 
+                          time - injection_start : ZERO;
+
+      // define box to inject into
+      boundaries_t<real_t> box;
+
+      // loop over all dimension
+      for (auto d = 0u; d < M::Dim; ++d) {
+        if (d == 0) {
+          box.push_back({ x_init + injector_velocity * dt_inj - 
+                          drift_ux / math::sqrt(1 + SQR(drift_ux)) * dt -
+                          1.5 * injection_frequency * dt,
+                          x_init + injector_velocity * (dt_inj + dt) });
+        } else {
+          box.push_back(Range::All);
+        }
+      }
+
+      // define indice range to reset fields
+      boundaries_t<bool> incl_ghosts;
+      for (auto d = 0; d < M::Dim; ++d) {
+        incl_ghosts.push_back({ true, true });
+      }
+      const auto extent = domain.mesh.ExtentToRange(box, incl_ghosts);
+      tuple_t<std::size_t, M::Dim> x_min { 0 }, x_max { 0 };
+      for (auto d = 0; d < M::Dim; ++d) {
+        x_min[d] = extent[d].first;
+        x_max[d] = extent[d].second;
+      }
+
+      // reset fields
+      std::vector<unsigned short> comps = { em::ex1, em::ex2, em::ex3,
+                                            em::bx1, em::bx2, em::bx3 };
+
+      // loop over all components
+      for (const auto& comp : comps) {
+        auto value = ResetFields((em)comp);
+
+        if constexpr (M::Dim == Dim::_1D) {
+          Kokkos::deep_copy(Kokkos::subview(domain.fields.em,
+                                            std::make_pair(x_min[0], x_max[0]),
+                                            comp),
+                            value);
+        } else if constexpr (M::Dim == Dim::_2D) {
+          Kokkos::deep_copy(Kokkos::subview(domain.fields.em,
+                                            std::make_pair(x_min[0], x_max[0]),
+                                            std::make_pair(x_min[1], x_max[1]),
+                                            comp),
+                            value);
+        } else if constexpr (M::Dim == Dim::_3D) {
+          Kokkos::deep_copy(Kokkos::subview(domain.fields.em,
+                                            std::make_pair(x_min[0], x_max[0]),
+                                            std::make_pair(x_min[1], x_max[1]),
+                                            std::make_pair(x_min[2], x_max[2]),
+                                            comp),
+                            value);
+        } else {
+          raise::Error("Invalid dimension", HERE);
+        }
+      }
+
+      /* 
+        tag particles inside the injection zone as dead 
+      */
+
+      // loop over particle species
+      for (std::size_t s { 0 }; s < 2; ++s) {
+
+        // get particle properties
+        auto& species = domain.species[s];
+        auto i1 = species.i1;
+        auto tag = species.tag;
+
+
+        // tag all particles with x > box[0].first as dead
+        Kokkos::parallel_for(
+          "RemoveParticles",
+          species.rangeActiveParticles(),
+          Lambda(index_t p) {
+            // check if the particle is already dead
+            if (tag(p) == ParticleTag::dead) {
+              return;
+            }
+            // select the x-coordinate index
+            auto x_i1 = i1(p);
+            // check if the particle is inside the box of new plasma
+            if (x_i1 > x_min[0]) {
+              tag(p) = ParticleTag::dead;
+            }
+          }
+        );
+      }
+
+      /* 
+          Inject piston of fresh plasma
+      */
+
+      // same maxwell distribution as above
+      const auto energy_dist  = arch::Maxwellian<S, M>(domain.mesh.metric,
+                                                      domain.random_pool,
+                                                      temperature,
+                                                      -drift_ux,
+                                                      in::x1);
+      // spatial distribution of the particles
+      // -> hack to use the uniform distribution in NonUniformInjector
+      const auto spatial_dist = arch::Piston<S, M>(domain.mesh.metric,
+                                                   box[0].first,
+                                                   box[0].second,
+                                                   in::x1);
+
+      // inject piston of fresh plasma
+      const auto injector = arch::NonUniformInjector<S, M, arch::Maxwellian, arch::Piston>(
+        energy_dist,
+        spatial_dist,
+        { 1, 2 });
+
+      // inject non-uniformly within the defined box
+      arch::InjectNonUniform<S, M, decltype(injector)>(params,
+                                                       domain,
+                                                       injector,
+                                                       ONE,
+                                                       false,
+                                                       box);
+
+      /*
+          I thought this option would be better, but I can't get it to work
+      */
+
+      //   const auto spatial_dist = arch::Replenish<S, M,
+      //   decltype(TargetDensityProfile)>(domain.mesh.metric,
+      //                                                   domain.fields.bckp,
+      //                                                   box,
+      //                                                   TargetDensity,
+      //                                                   1.0);
+
+      //   const auto injector = arch::NonUniformInjector<S, M, arch::Maxwellian, arch::Replenish>(
+      //       energy_dist,
+      //       spatial_dist,
+      //       {1, 2});
+
+      // const auto injector = arch::MovingInjector<S, M, in::x1> {
+      //   domain.mesh.metric,
+      //   domain.fields.bckp,
+      //   energy_dist,
+      //   box[0].first,
+      //   box[0].second,
+      //   1.0,
+      //   { 1, 2 }
+      // };
     }
   };