Merge pull request ComputationalRadiationPhysics#888 from ax3l/topic-currentInterpolation

Current Interpolation: During add to EMF

Author: psychocoderHPC

src/picongpu/include/fields/FieldJ.hpp

@@ -1,5 +1,5 @@
 /**
- * Copyright 2013-2014 Axel Huebl, Heiko Burau, Rene Widera, Richard Pausch
+ * Copyright 2013-2015 Axel Huebl, Heiko Burau, Rene Widera, Richard Pausch
  *
  * This file is part of PIConGPU.
  *
@@ -70,7 +70,7 @@ class FieldJ : public SimulationFieldHelper<MappingDesc>, public ISimulationData
 
     virtual EventTask asyncCommunication(EventTask serialEvent);
 
-    void init(FieldE &fieldE);
+    void init(FieldE &fieldE, FieldB &fieldB);
 
     GridLayout<simDim> getGridLayout();
 
@@ -87,8 +87,8 @@ class FieldJ : public SimulationFieldHelper<MappingDesc>, public ISimulationData
     template<uint32_t AREA, class ParticlesClass>
     void computeCurrent(ParticlesClass &parClass, uint32_t currentStep) throw (std::invalid_argument);
 
-    template<uint32_t AREA>
-    void addCurrentToE();
+    template<uint32_t AREA, class T_CurrentInterpolation>
+    void addCurrentToEMF( T_CurrentInterpolation& myCurrentInterpolation );
 
     SimulationDataId getUniqueId();
 
@@ -124,8 +124,10 @@ class FieldJ : public SimulationFieldHelper<MappingDesc>, public ISimulationData
 private:
 
     GridBuffer<ValueType, simDim> fieldJ;
+    GridBuffer<ValueType, simDim>* fieldJrecv;
 
     FieldE *fieldE;
+    FieldB *fieldB;
 };
 
 template<typename T_SpeciesName, typename T_Area>

src/picongpu/include/fields/FieldJ.kernel

@@ -1,5 +1,5 @@
 /**
- * Copyright 2013-2014 Axel Huebl, Heiko Burau, Rene Widera
+ * Copyright 2013-2015 Axel Huebl, Heiko Burau, Rene Widera
  *
  * This file is part of PIConGPU.
  *
@@ -160,24 +160,45 @@ private:
     const PMACC_ALIGN(deltaTime, float);
 };
 
-template<class Mapping>
-__global__ void kernelAddCurrentToE(typename FieldE::DataBoxType fieldE,
-                                    J_DataBox fieldJ,
-                                    Mapping mapper)
+template<class T_CurrentInterpolation, class T_Mapping>
+__global__ void kernelAddCurrentToEMF(typename FieldE::DataBoxType fieldE,
+                                      typename FieldB::DataBoxType fieldB,
+                                      J_DataBox fieldJ,
+                                      T_CurrentInterpolation currentInterpolation,
+                                      T_Mapping mapper)
 {
+    /* Caching of fieldJ */
+    typedef SuperCellDescription<
+                SuperCellSize,
+                typename T_CurrentInterpolation::LowerMargin,
+                typename T_CurrentInterpolation::UpperMargin
+                > BlockArea;
 
-    const DataSpace<simDim> blockCell(
-                                      mapper.getSuperCellIndex(DataSpace<simDim > (blockIdx))
-                                      * Mapping::SuperCellSize::toRT()
-                                      );
-    const DataSpace<Mapping::Dim> cell(blockCell + DataSpace<simDim > (threadIdx));
+    PMACC_AUTO(cachedJ, CachedBox::create < 0, typename J_DataBox::ValueType > (BlockArea()));
+
+    nvidia::functors::Assign assign;
+    const DataSpace<simDim> block(mapper.getSuperCellIndex(DataSpace<simDim > (blockIdx)));
+    const DataSpace<simDim> blockCell = block * MappingDesc::SuperCellSize::toRT();
+
+    const DataSpace<simDim > threadIndex(threadIdx);
+    PMACC_AUTO(fieldJBlock, fieldJ.shift(blockCell));
+
+    ThreadCollective<BlockArea> collectiv(threadIndex);
+    collectiv(
+              assign,
+              cachedJ,
+              fieldJBlock
+              );
+
+    __syncthreads();
+
+    const DataSpace<T_Mapping::Dim> cell(blockCell + threadIndex);
 
     // Amperes Law:
     //   Change of the dE = - j / EPS0 * dt
     //                        j = current density (= current per area)
     //                          = fieldJ
-    const float_X deltaT = DELTA_T;
-    fieldE(cell) -= fieldJ(cell) * (float_X(1.0) / EPS0) * deltaT;
+    currentInterpolation( fieldE.shift(cell), fieldB.shift(cell), cachedJ.shift(threadIndex) );
 }
 
 template<class Mapping>

src/picongpu/include/fields/FieldJ.tpp

@@ -1,5 +1,5 @@
 /**
- * Copyright 2013-2014 Axel Huebl, Heiko Burau, Rene Widera, Felix Schmitt,
+ * Copyright 2013-2015 Axel Huebl, Heiko Burau, Rene Widera, Felix Schmitt,
  *                     Richard Pausch
  *
  * This file is part of PIConGPU.
@@ -53,21 +53,35 @@ using namespace PMacc;
 
 FieldJ::FieldJ( MappingDesc cellDescription ) :
 SimulationFieldHelper<MappingDesc>( cellDescription ),
-fieldJ( cellDescription.getGridLayout( ) ), fieldE( NULL )
+fieldJ( cellDescription.getGridLayout( ) ), fieldE( NULL ), fieldB( NULL ), fieldJrecv( NULL )
 {
     const DataSpace<simDim> coreBorderSize = cellDescription.getGridLayout( ).getDataSpaceWithoutGuarding( );
 
-
+    /* cell margins the current might spread to due to particle shapes */
     typedef typename bmpl::accumulate<
         VectorAllSpecies,
         typename PMacc::math::CT::make_Int<simDim, 0>::type,
         PMacc::math::CT::max<bmpl::_1, GetLowerMargin< GetCurrentSolver<bmpl::_2> > >
-        >::type LowerMargin;
+        >::type LowerMarginShapes;
 
     typedef typename bmpl::accumulate<
         VectorAllSpecies,
         typename PMacc::math::CT::make_Int<simDim, 0>::type,
         PMacc::math::CT::max<bmpl::_1, GetUpperMargin< GetCurrentSolver<bmpl::_2> > >
+        >::type UpperMarginShapes;
+
+    /* margins are always positive, also for lower margins
+     * additional current interpolations and current filters on FieldJ might
+     * spread the dependencies on neighboring cells
+     *   -> use max(shape,filter) */
+    typedef typename PMacc::math::CT::max<
+        LowerMarginShapes,
+        GetMargin<fieldSolver::CurrentInterpolation>::LowerMargin
+        >::type LowerMargin;
+
+    typedef typename PMacc::math::CT::max<
+        UpperMarginShapes,
+        GetMargin<fieldSolver::CurrentInterpolation>::UpperMargin
         >::type UpperMargin;
 
     const DataSpace<simDim> originGuard( LowerMargin( ).toRT( ) );
@@ -101,10 +115,34 @@ fieldJ( cellDescription.getGridLayout( ) ), fieldE( NULL )
         // std::cout << "ex " << i << " x=" << guardingCells[0] << " y=" << guardingCells[1] << " z=" << guardingCells[2] << std::endl;
         fieldJ.addExchangeBuffer( i, guardingCells, FIELD_J );
     }
+
+    /* Receive border values in own guard for "receive" communication pattern - necessary for current interpolation/filter */
+    const DataSpace<simDim> originRecvGuard( GetMargin<fieldSolver::CurrentInterpolation>::LowerMargin( ).toRT( ) );
+    const DataSpace<simDim> endRecvGuard( GetMargin<fieldSolver::CurrentInterpolation>::UpperMargin( ).toRT( ) );
+    if( originRecvGuard != DataSpace<simDim>::create(0) ||
+        endRecvGuard != DataSpace<simDim>::create(0) )
+    {
+        fieldJrecv = new GridBuffer<ValueType, simDim > ( fieldJ.getDeviceBuffer(), cellDescription.getGridLayout( ) );
+
+        /*go over all directions*/
+        for ( uint32_t i = 1; i < NumberOfExchanges<simDim>::value; ++i )
+        {
+            DataSpace<simDim> relativMask = Mask::getRelativeDirections<simDim > ( i );
+            /* guarding cells depend on direction
+             * for negative direction use originGuard else endGuard (relative direction ZERO is ignored)
+             * don't switch end and origin because this is a read buffer and no send buffer
+             */
+            DataSpace<simDim> guardingCells;
+            for ( uint32_t d = 0; d < simDim; ++d )
+                guardingCells[d] = ( relativMask[d] == -1 ? originRecvGuard[d] : endRecvGuard[d] );
+            fieldJrecv->addExchange( GUARD, i, guardingCells, FIELD_JRECV );
+        }
+    }
 }
 
 FieldJ::~FieldJ( )
 {
+    __delete(fieldJrecv);
 }
 
 SimulationDataId FieldJ::getUniqueId( )
@@ -132,7 +170,14 @@ EventTask FieldJ::asyncCommunication( EventTask serialEvent )
     __startTransaction( serialEvent );
     FieldFactory::getInstance( ).createTaskFieldSend( *this );
     ret += __endTransaction( );
-    return ret;
+
+    if( fieldJrecv != NULL )
+    {
+        EventTask eJ = fieldJrecv->asyncCommunication( ret );
+        return eJ;
+    }
+    else
+        return ret;
 }
 
 void FieldJ::bashField( uint32_t exchangeType )
@@ -166,9 +211,10 @@ void FieldJ::insertField( uint32_t exchangeType )
           direction, mapper );
 }
 
-void FieldJ::init( FieldE &fieldE )
+void FieldJ::init( FieldE &fieldE, FieldB &fieldB )
 {
     this->fieldE = &fieldE;
+    this->fieldB = &fieldB;
 
     Environment<>::get( ).DataConnector( ).registerData( *this );
 }
@@ -250,15 +296,17 @@ void FieldJ::computeCurrent( ParticlesClass &parClass, uint32_t ) throw (std::in
     __setTransactionEvent( __endTransaction( ) );
 }
 
-template<uint32_t AREA>
-void FieldJ::addCurrentToE( )
+template<uint32_t AREA, class T_CurrentInterpolation>
+void FieldJ::addCurrentToEMF( T_CurrentInterpolation& myCurrentInterpolation )
 {
-    __picKernelArea( ( kernelAddCurrentToE ),
+    __picKernelArea( ( kernelAddCurrentToEMF ),
                      cellDescription,
                      AREA )
         ( MappingDesc::SuperCellSize::toRT( ).toDim3( ) )
         ( this->fieldE->getDeviceDataBox( ),
-          this->fieldJ.getDeviceBuffer( ).getDataBox( ) );
+          this->fieldB->getDeviceDataBox( ),
+          this->fieldJ.getDeviceBuffer( ).getDataBox( ),
+          myCurrentInterpolation );
 }
 
 }

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

@@ -0,0 +1,54 @@
+/**
+ * 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
+{
+
+/* 2nd order Binomial filter */
+template<uint32_t T_dim>
+struct Binomial;
+
+} /* 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::Binomial<T_dim > >
+{
+private:
+    typedef picongpu::currentInterpolation::Binomial<T_dim> MyInterpolation;
+
+public:
+    typedef typename MyInterpolation::LowerMargin LowerMargin;
+    typedef typename MyInterpolation::UpperMargin UpperMargin;
+};
+
+} /* namespace traits */
+
+} /* namespace picongpu */

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

@@ -0,0 +1,79 @@
+/**
+ * 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"
+
+namespace picongpu
+{
+namespace currentInterpolation
+{
+using namespace PMacc;
+
+template<uint32_t T_dim>
+struct Binomial
+{
+    static const uint32_t dim = T_dim;
+
+    typedef typename PMacc::math::CT::make_Int<dim, 1>::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,
+                             DataBoxJ fieldJ )
+    {
+        const DataSpace<dim> self;
+        typedef typename DataBoxJ::ValueType TypeJ;
+
+        /* 1 2 1 weighting for "left"(1x) "center"(2x) "right"(1x),
+         * see Pascal's triangle level N=2 */
+        TypeJ dirSum( TypeJ::create(0.0) );
+        for( uint32_t d = 0; d < dim; ++d )
+        {
+            DataSpace<dim> dw;
+            dw[d] = -1;
+            DataSpace<dim> up;
+            up[d] =  1;
+            const TypeJ dirDw = fieldJ(dw) + fieldJ(self);
+            const TypeJ dirUp = fieldJ(up) + fieldJ(self);
+
+            /* each fieldJ component is added individually */
+            dirSum += dirDw + dirUp;
+        }
+
+        /* component-wise division by sum of all weightings,
+         * in the second order binomial filter these are 4 values per direction
+         * (1D: 4 values; 2D: 8 values; 3D: 12 values) */
+        const TypeJ filteredJ = dirSum / TypeJ::create(4.0 * dim);
+
+        const float_X deltaT = DELTA_T;
+        fieldE(self) -= filteredJ * (float_X(1.0) / EPS0) * deltaT;
+    }
+};
+
+} /* namespace currentInterpolation */
+
+} /* namespace picongpu */

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

@@ -0,0 +1,23 @@
+/**
+ * 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/>.
+ */
+
+
+#include "fields/currentInterpolation/None/None.def"
+#include "fields/currentInterpolation/Binomial/Binomial.def"