Current Interpolation: During add to EMF

ax3l · ax3l · commit c2709369b6bb · 2015-05-26T13:40:28.000+02:00
This commit adds the general infrastructure to call functors
before adding `FieldJ` to the electro-magnetic fields (in
Yee: before manipulating the `E` field).

The previous method `None` without interpolation is added as a
default in `fieldSolver.param`. Additionally, a simple second order
Binomial filter for 2D and 3D is added to test the new functionality.
diff --git a/src/picongpu/include/fields/FieldJ.hpp b/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();
 
@@ -126,6 +126,7 @@ class FieldJ : public SimulationFieldHelper<MappingDesc>, public ISimulationData
     GridBuffer<ValueType, simDim> fieldJ;
 
     FieldE *fieldE;
+    FieldB *fieldB;
 };
 
 template<typename T_SpeciesName, typename T_Area>
diff --git a/src/picongpu/include/fields/FieldJ.kernel b/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>
diff --git a/src/picongpu/include/fields/FieldJ.tpp b/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,7 +53,7 @@ using namespace PMacc;
 
 FieldJ::FieldJ( MappingDesc cellDescription ) :
 SimulationFieldHelper<MappingDesc>( cellDescription ),
-fieldJ( cellDescription.getGridLayout( ) ), fieldE( NULL )
+fieldJ( cellDescription.getGridLayout( ) ), fieldE( NULL ), fieldB( NULL )
 {
     const DataSpace<simDim> coreBorderSize = cellDescription.getGridLayout( ).getDataSpaceWithoutGuarding( );
 
@@ -166,9 +166,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 +251,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 );
 }
 
 }
diff --git a/src/picongpu/include/fields/currentInterpolation/Binomial/Binomial.def b/src/picongpu/include/fields/currentInterpolation/Binomial/Binomial.def
@@ -0,0 +1,34 @@
+/**
+ * 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 picongpu */
diff --git a/src/picongpu/include/fields/currentInterpolation/Binomial/Binomial.hpp b/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 */
diff --git a/src/picongpu/include/fields/currentInterpolation/CurrentInterpolation.def b/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"
diff --git a/src/picongpu/include/fields/currentInterpolation/CurrentInterpolation.hpp b/src/picongpu/include/fields/currentInterpolation/CurrentInterpolation.hpp
@@ -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.hpp"
+#include "fields/currentInterpolation/Binomial/Binomial.hpp"
diff --git a/src/picongpu/include/fields/currentInterpolation/None/None.def b/src/picongpu/include/fields/currentInterpolation/None/None.def
@@ -0,0 +1,33 @@
+/**
+ * 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
+{
+
+template<uint32_t T_dim>
+struct None;
+
+} /* namespace currentInterpolation */
+
+} /* namespace picongpu */
diff --git a/src/picongpu/include/fields/currentInterpolation/None/None.hpp b/src/picongpu/include/fields/currentInterpolation/None/None.hpp
diff --git a/src/picongpu/include/simulationControl/MySimulation.hpp b/src/picongpu/include/simulationControl/MySimulation.hpp
diff --git a/src/picongpu/include/simulation_defines/param/fieldSolver.param b/src/picongpu/include/simulation_defines/param/fieldSolver.param

Original file line number	Diff line number	Diff line change
`@@ -1,5 +1,5 @@`
`1`	`1`	`/**`
`2`		`- * Copyright 2013-2014 Axel Huebl, Heiko Burau, Rene Widera, Felix Schmitt,`
	`2`	`+ * Copyright 2013-2015 Axel Huebl, Heiko Burau, Rene Widera, Felix Schmitt,`
`3`	`3`	`* Richard Pausch`
`4`	`4`	`*`
`5`	`5`	`* This file is part of PIConGPU.`
`@@ -53,7 +53,7 @@ using namespace PMacc;`
`53`	`53`
`54`	`54`	`FieldJ::FieldJ( MappingDesc cellDescription ) :`
`55`	`55`	`SimulationFieldHelper<MappingDesc>( cellDescription ),`
`56`		`-fieldJ( cellDescription.getGridLayout( ) ), fieldE( NULL )`
	`56`	`+fieldJ( cellDescription.getGridLayout( ) ), fieldE( NULL ), fieldB( NULL )`
`57`	`57`	`{`
`58`	`58`	`const DataSpace<simDim> coreBorderSize = cellDescription.getGridLayout( ).getDataSpaceWithoutGuarding( );`
`59`	`59`
`@@ -166,9 +166,10 @@ void FieldJ::insertField( uint32_t exchangeType )`
`166`	`166`	`direction, mapper );`
`167`	`167`	`}`
`168`	`168`
`169`		`-void FieldJ::init( FieldE &fieldE )`
	`169`	`+void FieldJ::init( FieldE &fieldE, FieldB &fieldB )`
`170`	`170`	`{`
`171`	`171`	`this->fieldE = &fieldE;`
	`172`	`+ this->fieldB = &fieldB;`
`172`	`173`
`173`	`174`	`Environment<>::get( ).DataConnector( ).registerData( *this );`
`174`	`175`	`}`
`@@ -250,15 +251,17 @@ void FieldJ::computeCurrent( ParticlesClass &parClass, uint32_t ) throw (std::in`
`250`	`251`	`__setTransactionEvent( __endTransaction( ) );`
`251`	`252`	`}`
`252`	`253`
`253`		`-template<uint32_t AREA>`
`254`		`-void FieldJ::addCurrentToE( )`
	`254`	`+template<uint32_t AREA, class T_CurrentInterpolation>`
	`255`	`+void FieldJ::addCurrentToEMF( T_CurrentInterpolation& myCurrentInterpolation )`
`255`	`256`	`{`
`256`		`- __picKernelArea( ( kernelAddCurrentToE ),`
	`257`	`+ __picKernelArea( ( kernelAddCurrentToEMF ),`
`257`	`258`	`cellDescription,`
`258`	`259`	`AREA )`
`259`	`260`	`( MappingDesc::SuperCellSize::toRT( ).toDim3( ) )`
`260`	`261`	`( this->fieldE->getDeviceDataBox( ),`
`261`		`- this->fieldJ.getDeviceBuffer( ).getDataBox( ) );`
	`262`	`+ this->fieldB->getDeviceDataBox( ),`
	`263`	`+ this->fieldJ.getDeviceBuffer( ).getDataBox( ),`
	`264`	`+ myCurrentInterpolation );`
`262`	`265`	`}`
`263`	`266`
`264`	`267`	`}`