Direction Splitting: Self-Consistent J & EMF

src/picongpu/include/algorithms/LinearInterpolateWithUpper.hpp

@@ -0,0 +1,86 @@
+/**
+ * Copyright 2015 Heiko Burau, Rene Widera
+ *
+ * This file is part of PIConGPU.
+ *
+ * PIConGPU is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * PIConGPU is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with PIConGPU.
+ * If not, see <http://www.gnu.org/licenses/>.
+ */
+
+
+
+#pragma once
+
+#include "types.h"
+#include "math/Vector.hpp"
+
+namespace picongpu
+{
+
+/** Calculate linear interpolation to upper cell value
+ *
+ * @tparam T_Dim for how many dimensions does this operator interpolate
+ *
+ * If `GetDifference` is called for a direction greater or equal T_Dim,
+ * a zeroed value is returned (assumes symmetry in those directions).
+ */
+template<uint32_t T_Dim>
+struct LinearInterpolateWithUpper
+{
+    static const uint32_t dim = T_Dim;
+
+    typedef typename PMacc::math::CT::make_Int<dim, 0>::type OffsetOrigin;
+    typedef typename PMacc::math::CT::make_Int<dim, 1>::type OffsetEnd;
+
+    /** calculate the linear interpolation for a given direction
+     *
+     * @tparam T_direction direction for the interpolation operation
+     * @tparam T_isLesserThanDim not needed/ this is calculated by the compiler
+     */
+    template<uint32_t T_direction, bool T_isLesserThanDim = (T_direction < dim)>
+    struct GetInterpolatedValue
+    {
+        static const uint32_t direction = T_direction;
+
+        /** get interpolated value
+         * @return interpolated value
+         */
+        template<class Memory >
+        HDINLINE typename Memory::ValueType operator()(const Memory& mem) const
+        {
+            const DataSpace<dim> indexIdentity; /* defaults to (0, 0, 0) in 3D */
+            DataSpace<dim> indexUpper; /* e.g., (0, 1, 0) for direction y in 3D */
+            indexUpper[direction] = 1;
+
+            return ( mem(indexUpper) + mem(indexIdentity)) * Memory::ValueType::create(0.5);
+        }
+    };
+
+    /** special case for `direction >= simulation dimensions`*/
+    template<uint32_t T_direction>
+    struct GetInterpolatedValue<T_direction, false>
+    {
+
+        /** @return always identity
+         */
+        template<class Memory >
+        HDINLINE typename Memory::ValueType operator()(const Memory& mem) const
+        {
+            return *mem;
+        }
+    };
+
+};
+
+} //namespace picongpu

src/picongpu/include/fields/MaxwellSolver/DirSplitting/DirSplitting.hpp

@@ -150,7 +150,20 @@ class DirSplitting : private ConditionCheck<fieldSolver::FieldSolver>
     }
 
     void update_afterCurrent(uint32_t) const
-    {    }
+    {
+        DataConnector &dc = Environment<>::get().DataConnector();
+
+        FieldE& fieldE = dc.getData<FieldE > (FieldE::getName(), true);
+        FieldB& fieldB = dc.getData<FieldB > (FieldB::getName(), true);
+
+        EventTask eRfieldE = fieldE.asyncCommunication(__getTransactionEvent());
+        EventTask eRfieldB = fieldB.asyncCommunication(__getTransactionEvent());
+        __setTransactionEvent(eRfieldE);
+        __setTransactionEvent(eRfieldB);
+
+        dc.releaseData(FieldE::getName());
+        dc.releaseData(FieldB::getName());
+    }
 };
 
 } // dirSplitting

src/picongpu/include/fields/currentDeposition/ZigZag/ZigZag.hpp

@@ -195,7 +195,7 @@ struct ZigZag
              *   - run calculations in a shape optimized coordinate system
              *     with fixed interpolation points
              */
-            ShiftCoordinateSystem<Supports_direction>()(cursor, pos, fieldSolver::NumericalCellType::getEFieldPosition()[dir]);
+            ShiftCoordinateSystem<Supports_direction>()(cursor, pos, fieldSolver::NumericalCellType::getJFieldPosition()[dir]);
 
             /* define grid points where we evaluate the shape function*/
             typedef typename PMacc::math::CT::make_Vector<

src/picongpu/include/fields/currentInterpolation/CurrentInterpolation.def

@@ -21,3 +21,4 @@
 
 #include "fields/currentInterpolation/None/None.def"
 #include "fields/currentInterpolation/Binomial/Binomial.def"
+#include "fields/currentInterpolation/NoneDS/NoneDS.def"

src/picongpu/include/fields/currentInterpolation/CurrentInterpolation.hpp

@@ -21,3 +21,4 @@
 
 #include "fields/currentInterpolation/None/None.hpp"
 #include "fields/currentInterpolation/Binomial/Binomial.hpp"
+#include "fields/currentInterpolation/NoneDS/NoneDS.hpp"

src/picongpu/include/fields/currentInterpolation/NoneDS/NoneDS.def

@@ -0,0 +1,58 @@
+/**
+ * Copyright 2015 Axel Huebl
+ *
+ * This file is part of PIConGPU.
+ *
+ * PIConGPU is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * PIConGPU is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with PIConGPU.
+ * If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+namespace picongpu
+{
+namespace currentInterpolation
+{
+
+/* The standard interpolation for Directional Splitting
+ *
+ * In the Directional Splitting scheme (see MaxwellSolver and
+ * numericalCellTypes
+ */
+template<uint32_t T_simDim>
+struct NoneDS;
+
+} /* namespace currentInterpolation */
+
+namespace traits
+{
+
+/* Get margin of the current interpolation
+ *
+ * This class defines a LowerMargin and an UpperMargin.
+ */
+template<uint32_t T_dim>
+struct GetMargin<picongpu::currentInterpolation::NoneDS<T_dim > >
+{
+private:
+    typedef picongpu::currentInterpolation::NoneDS<T_dim> MyInterpolation;
+
+public:
+    typedef typename MyInterpolation::LowerMargin LowerMargin;
+    typedef typename MyInterpolation::UpperMargin UpperMargin;
+};
+
+} /* namespace traits */
+
+} /* namespace picongpu */

src/picongpu/include/fields/currentInterpolation/NoneDS/NoneDS.hpp

@@ -0,0 +1,177 @@
+/**
+ * Copyright 2015 Axel Huebl
+ *
+ * This file is part of PIConGPU.
+ *
+ * PIConGPU is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * PIConGPU is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with PIConGPU.
+ * If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include "simulation_defines.hpp"
+#include "types.h"
+#include "basicOperations.hpp"
+
+#include "fields/currentInterpolation/None/None.def"
+#include "algorithms/DifferenceToUpper.hpp"
+#include "algorithms/LinearInterpolateWithUpper.hpp"
+#include "traits/GetComponentsType.hpp"
+
+#include "fields/MaxwellSolver/Yee/Curl.hpp"
+
+
+namespace picongpu
+{
+namespace currentInterpolation
+{
+using namespace PMacc;
+
+namespace detail
+{
+    template<uint32_t T_simDim, uint32_t T_plane>
+    struct LinearInterpolateComponentPlaneUpper
+    {
+        static const uint32_t dim = T_simDim;
+
+        /* UpperMargin is actually 0 in direction of T_plane */
+        typedef typename PMacc::math::CT::make_Int<dim, 0>::type LowerMargin;
+        typedef typename PMacc::math::CT::make_Int<dim, 1>::type UpperMargin;
+
+        template<typename DataBox>
+        HDINLINE float_X operator()(DataBox field) const
+        {
+            const DataSpace<dim> self;
+            DataSpace<dim> up;
+            up[(T_plane + 1) % dim] = 1;
+
+            typedef LinearInterpolateWithUpper<dim> Avg;
+
+            const typename Avg::template GetInterpolatedValue<(T_plane + 2) % dim> avg;
+
+            return float_X(0.5) * ( avg(field)[T_plane] + avg(field.shift(up))[T_plane] );
+        }
+    };
+
+    /* shift a databox along a specific direction
+     *
+     * returns the identity (assume periodic symmetry) if direction is not
+     * available, such as in a 2D simulation
+     *
+     * \todo accept a full CT::Vector and shift if possible
+     * \todo call with CT::Vector of correct dimensionality that was created
+     *       with AssignIfInRange...
+     *
+     * \tparam T_simDim maximum dimensionality of the mesh
+     * \tparam T_direction (0)X (1)Y or (2)Z for the direction one wants to
+     *                     shift to
+     * \tparam isShiftAble auto-filled value that decides if this direction
+     *                     is actually non-existent == periodic
+     */
+    template<uint32_t T_simDim, uint32_t T_direction, bool isShiftAble=(T_direction<T_simDim) >
+    struct ShiftMeIfYouCan
+    {
+        static const uint32_t dim = T_simDim;
+        static const uint32_t dir = T_direction;
+
+        template<class T_DataBox >
+        HDINLINE T_DataBox operator()(const T_DataBox& dataBox) const
+        {
+            DataSpace<dim> shift;
+            shift[dir] = 1;
+            return dataBox.shift(shift);
+        }
+    };
+
+    template<uint32_t T_simDim, uint32_t T_direction>
+    struct ShiftMeIfYouCan<T_simDim, T_direction, false>
+    {
+        template<class T_DataBox >
+        HDINLINE T_DataBox operator()(const T_DataBox& dataBox) const
+        {
+            return dataBox;
+        }
+    };
+
+    /* that is not a "real" yee curl, but it looks a bit like it */
+    template<class Difference>
+    struct ShiftCurl
+    {
+        typedef typename Difference::OffsetOrigin LowerMargin;
+        typedef typename Difference::OffsetEnd UpperMargin;
+
+        template<class DataBox >
+        HDINLINE typename DataBox::ValueType operator()(const DataBox& mem) const
+        {
+            const typename Difference::template GetDifference<0> Dx;
+            const typename Difference::template GetDifference<1> Dy;
+            const typename Difference::template GetDifference<2> Dz;
+
+            const ShiftMeIfYouCan<simDim, 0> sx;
+            const ShiftMeIfYouCan<simDim, 1> sy;
+            const ShiftMeIfYouCan<simDim, 2> sz;
+
+            return float3_X(Dy(sx(mem)).z() - Dz(sx(mem)).y(),
+                            Dz(sy(mem)).x() - Dx(sy(mem)).z(),
+                            Dx(sz(mem)).y() - Dy(sz(mem)).x());
+        }
+    };
+} /* namespace detail */
+
+template<uint32_t T_simDim>
+struct NoneDS
+{
+    static const uint32_t dim = T_simDim;
+
+    typedef typename PMacc::math::CT::make_Int<dim, 0>::type LowerMargin;
+    typedef typename PMacc::math::CT::make_Int<dim, 1>::type UpperMargin;
+
+    template<typename DataBoxE, typename DataBoxB, typename DataBoxJ>
+    HDINLINE void operator()(DataBoxE fieldE,
+                             DataBoxB fieldB,
+                             DataBoxJ fieldJ )
+    {
+        typedef typename DataBoxJ::ValueType TypeJ;
+        typedef typename GetComponentsType<TypeJ>::type ComponentJ;
+
+        const DataSpace<dim> self;
+
+        const ComponentJ deltaT = DELTA_T;
+        const ComponentJ constE = (float_X(1.0)  / EPS0) * deltaT;
+        const ComponentJ constB = (float_X(0.25) / EPS0) * deltaT * deltaT;
+
+        const detail::LinearInterpolateComponentPlaneUpper<dim, 0> avgX;
+        const ComponentJ jXavg = avgX(fieldJ);
+        const detail::LinearInterpolateComponentPlaneUpper<dim, 1> avgY;
+        const ComponentJ jYavg = avgY(fieldJ);
+        const detail::LinearInterpolateComponentPlaneUpper<dim, 2> avgZ;
+        const ComponentJ jZavg = avgZ(fieldJ);
+
+        const TypeJ jAvgE = TypeJ(jXavg, jYavg, jZavg);
+        fieldE(self) -= jAvgE * constE;
+
+        typedef yeeSolver::Curl<DifferenceToUpper<dim> > CurlRight;
+        typedef detail::ShiftCurl<DifferenceToUpper<dim> > ShiftCurlRight;
+        CurlRight curl;
+        ShiftCurlRight shiftCurl;
+
+        const TypeJ jAvgB = curl(fieldJ) + shiftCurl(fieldJ);
+        fieldB(self) += jAvgB * constB;
+    }
+
+};
+
+} /* namespace currentInterpolation */
+
+} /* namespace picongpu */