Merge pull request ComputationalRadiationPhysics#880 from psychocoderHPC/fix-fieldSolverAssert

move field solver condition check to field solver

Author: ax3l

examples/KelvinHelmholtz/cmakeFlags

@@ -36,6 +36,7 @@ flags[3]="-DCUDA_ARCH=sm_35 -DPARAM_OVERWRITES:LIST=-DPARAM_CURRENTSOLVER=ZigZag
 flags[4]="-DCUDA_ARCH=sm_35 -DPARAM_OVERWRITES:LIST=-DPARAM_CURRENTSOLVER=ZigZag<UsedParticleShape>;-DPARAM_PARTICLESHAPE=PCS"
 flags[5]="-DCUDA_ARCH=sm_35 -DPARAM_OVERWRITES:LIST=-DPARAM_CURRENTSOLVER=ZigZag<UsedParticleShape>;-DPARAM_PARTICLESHAPE=TSC"
 flags[6]="-DCUDA_ARCH=sm_35 -DPARAM_OVERWRITES:LIST=-DPARAM_CURRENTSOLVER=ZigZag<UsedParticleShape>;-DPARAM_PARTICLESHAPE=TSC;-DPARAM_DIMENSION=DIM2"
+flags[7]="-DCUDA_ARCH=sm_20 -DPARAM_OVERWRITES:LIST=-DPARAM_FIELDSOLVER=fieldSolverDirSplitting"
 
 
 ################################################################################

examples/KelvinHelmholtz/include/simulation_defines/param/componentsConfig.param

@@ -50,7 +50,13 @@ namespace picongpu
      *  - fieldSolverYeeNative : generic version of fieldSolverYee
      *    (need more shared memory per GPU and is slow)
      */
-    namespace fieldSolver = fieldSolverYee;
+#ifndef PARAM_FIELDSOLVER
+/* WARNING: if you change field solver by hand please update your CELL_WIDTH_SI
+ * in `gridConfig.param` to fulfill the convergence condition (CFL)
+ */
+#define PARAM_FIELDSOLVER fieldSolverYee
+#endif
+namespace fieldSolver = PARAM_FIELDSOLVER;
 
     /*enable (1) or disable (0) current calculation*/
 #define ENABLE_CURRENT 1

examples/KelvinHelmholtz/include/simulation_defines/param/gridConfig.param

@@ -1,23 +1,23 @@
 /**
  * Copyright 2013-2015 Axel Huebl, Rene Widera, Richard Pausch
  *
- * 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/>. 
- */ 
- 
+ * 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
@@ -32,7 +32,27 @@ namespace picongpu
         const double DELTA_T_SI = 1.79e-16;
         /** equals X
          *  unit: meter */
+#define fieldSolverDirSplitting 1
+#if (PARAM_FIELDSOLVER == 1)
+        /* THIS CODE PATH IS ONLY USED IF `PARAM_FIELDSOLVER` IS CHANGED IN
+         * `cmakeFlags` and the field solver there is set to fieldSolverDirSplitting
+         *
+         * Directional Splitting requires a fixed ratio between dt and dx
+         * and in addition cubic cells.
+         * conditions: dX == dt * c
+         *             dX == dY
+         *             dX == dZ
+         */
+        const double CELL_WIDTH_SI = DELTA_T_SI*SPEED_OF_LIGHT_SI;
+#else
+        /* cell size for Yee solver (must fulfill CFL)
+         * WARNING: if you change the field solver in `componentsConfig` you
+         * have to change the CELL_SIZE in this code path
+         */
         const double CELL_WIDTH_SI = 9.34635e-8;
+#endif
+#undef fieldSolverDirSplitting
+
         /** equals Y - the propagation direction
          *  unit: meter */
         const double CELL_HEIGHT_SI = CELL_WIDTH_SI;
@@ -47,22 +67,22 @@ namespace picongpu
         {32, 32},  /*y direction [negative,positive]*/
         {32, 32}   /*z direction [negative,positive]*/
     }; //unit: number of cells
-    
+
     //! Define the strange of the absober for any direction
     const float_X ABSORBER_STRENGTH[3][2] = {
         {1.0e-3, 1.0e-3}, /*x direction [negative,positive]*/
         {1.0e-3, 1.0e-3}, /*y direction [negative,positive]*/
         {1.0e-3, 1.0e-3}  /*z direction [negative,positive]*/
-    }; //unit: none 
+    }; //unit: none
 
     const uint32_t ABSORBER_FADE_IN_STEPS = 16;
-    
+
     /** When to move the co-moving window.
      *  An initial pseudo particle, flying with the speed of light,
      *  is fired at the begin of the simulation.
      *  When it reaches slide_point % of the absolute(*) simulation area,
      *  the co-moving window starts to move with the speed of light.
-     * 
+     *
      *  (*) Note: beware, that there is one "hidden" row of gpus at the y-front,
      *            when you use the co-moving window
      *  0.75 means only 75% of simulation area is used for real simulation

examples/LaserWakefield/cmakeFlags

@@ -31,17 +31,15 @@
 
 flags[0]="-DCUDA_ARCH=sm_20"
 flags[1]="-DCUDA_ARCH=sm_20 -DPARAM_OVERWRITES:LIST=-DPARAM_FIELDSOLVER=fieldSolverLehe;-DPARAM_PARTICLEPUSHER=Vay"
-flags[2]="-DCUDA_ARCH=sm_20 -DPARAM_OVERWRITES:LIST=-DPARAM_FIELDSOLVER=fieldSolverDirSplitting"
-flags[3]="-DCUDA_ARCH=sm_20 -DPARAM_OVERWRITES:LIST=-DPARAM_CURRENTSOLVER=Esirkepov;-DPARAM_PARTICLESHAPE=CIC"
-flags[4]="-DCUDA_ARCH=sm_20 -DPARAM_OVERWRITES:LIST=-DPARAM_CURRENTSOLVER=VillaBune;-DPARAM_PARTICLESHAPE=CIC"
-flags[5]="-DCUDA_ARCH=sm_20 -DPARAM_OVERWRITES:LIST=-DPARAM_PRECISION=precision64Bit"
-flags[6]="-DCUDA_ARCH=sm_20 -DPARAM_OVERWRITES:LIST=-DPARAM_FIELDSOLVER=fieldSolverDirSplitting;-DPARAM_PRECISION=precision64Bit"
-flags[7]="-DCUDA_ARCH=sm_35 -DPARAM_OVERWRITES:LIST=-DENABLE_CURRENT=0"
-flags[8]="-DCUDA_ARCH=sm_20 -DPARAM_OVERWRITES:LIST=-DPARAM_DIMENSION=DIM2"
-flags[9]="-DCUDA_ARCH=sm_20 -DPARAM_OVERWRITES:LIST=-DPARAM_CURRENTSOLVER=Esirkepov;-DPARAM_PARTICLESHAPE=CIC;-DPARAM_DIMENSION=DIM2"
-flags[10]="-DCUDA_ARCH=sm_20 -DPARAM_OVERWRITES:LIST=-DPARAM_PRECISION=precision64Bit;-DPARAM_DIMENSION=DIM2"
-flags[11]="-DCUDA_ARCH=sm_35 -DPARAM_OVERWRITES:LIST=-DENABLE_CURRENT=0;-DPARAM_DIMENSION=DIM2"
-flags[12]="-DCUDA_ARCH=sm_35 -DPARAM_OVERWRITES:LIST=-DPARAM_IONS=1;-DPARAM_IONIZATION=1"
+flags[2]="-DCUDA_ARCH=sm_20 -DPARAM_OVERWRITES:LIST=-DPARAM_CURRENTSOLVER=Esirkepov;-DPARAM_PARTICLESHAPE=CIC"
+flags[3]="-DCUDA_ARCH=sm_20 -DPARAM_OVERWRITES:LIST=-DPARAM_CURRENTSOLVER=VillaBune;-DPARAM_PARTICLESHAPE=CIC"
+flags[4]="-DCUDA_ARCH=sm_20 -DPARAM_OVERWRITES:LIST=-DPARAM_PRECISION=precision64Bit"
+flags[5]="-DCUDA_ARCH=sm_35 -DPARAM_OVERWRITES:LIST=-DENABLE_CURRENT=0"
+flags[6]="-DCUDA_ARCH=sm_20 -DPARAM_OVERWRITES:LIST=-DPARAM_DIMENSION=DIM2"
+flags[7]="-DCUDA_ARCH=sm_20 -DPARAM_OVERWRITES:LIST=-DPARAM_CURRENTSOLVER=Esirkepov;-DPARAM_PARTICLESHAPE=CIC;-DPARAM_DIMENSION=DIM2"
+flags[8]="-DCUDA_ARCH=sm_20 -DPARAM_OVERWRITES:LIST=-DPARAM_PRECISION=precision64Bit;-DPARAM_DIMENSION=DIM2"
+flags[9]="-DCUDA_ARCH=sm_35 -DPARAM_OVERWRITES:LIST=-DENABLE_CURRENT=0;-DPARAM_DIMENSION=DIM2"
+flags[10]="-DCUDA_ARCH=sm_35 -DPARAM_OVERWRITES:LIST=-DPARAM_IONS=1;-DPARAM_IONIZATION=1"
 ################################################################################
 # execution
 

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

@@ -40,7 +40,34 @@ namespace dirSplitting
 {
 using namespace PMacc;
 
-class DirSplitting
+/** Check Directional Splitting grid and time conditions
+ *
+ * This is a workaround that the condition check is only
+ * triggered if the current used solver is `DirSplitting`
+ */
+template<typename T_UsedSolver, typename T_Dummy=void>
+struct ConditionCheck
+{
+};
+
+template<typename T_Dummy>
+struct ConditionCheck<DirSplitting, T_Dummy>
+{
+    /* Directional Splitting conditions:
+     *
+     * using SI units to avoid round off errors
+     */
+    PMACC_CASSERT_MSG(DirectionSplitting_Set_dX_equal_dt_times_c____check_your_gridConfig_param_file,
+                      (SI::SPEED_OF_LIGHT_SI * SI::DELTA_T_SI) == SI::CELL_WIDTH_SI);
+    PMACC_CASSERT_MSG(DirectionSplitting_use_cubic_cells____check_your_gridConfig_param_file,
+                      SI::CELL_HEIGHT_SI == SI::CELL_WIDTH_SI);
+#if (SIMDIM == DIM3)
+    PMACC_CASSERT_MSG(DirectionSplitting_use_cubic_cells____check_your_gridConfig_param_file,
+                      SI::CELL_DEPTH_SI == SI::CELL_WIDTH_SI);
+#endif
+};
+
+class DirSplitting : private ConditionCheck<fieldSolver::FieldSolver>
 {
 private:
     template<typename CursorE, typename CursorB, typename GridSize>

src/picongpu/include/fields/MaxwellSolver/None/NoSolver.hpp

@@ -34,26 +34,31 @@ namespace picongpu
     {
         using namespace PMacc;
 
+        /** Check Yee grid and time conditions
+         *
+         * This is a workaround that the condition check is only
+         * triggered if the current used solver is `NoSolver`
+         */
+        template<typename T_UsedSolver, typename T_Dummy=void>
+        struct ConditionCheck
+        {
+        };
+
+        template<typename T_Dummy>
+        struct ConditionCheck<NoSolver, T_Dummy>
+        {
+            /* Courant-Friedrichs-Levy-Condition for Yee Field Solver: */
+            PMACC_CASSERT_MSG(Courant_Friedrichs_Levy_condition_failure____check_your_gridConfig_param_file,
+                (SPEED_OF_LIGHT*SPEED_OF_LIGHT*DELTA_T*DELTA_T*INV_CELL2_SUM)<=1.0);
+        };
 
-        class NoSolver
+        class NoSolver : private ConditionCheck<fieldSolver::FieldSolver>
         {
         private:
             typedef MappingDesc::SuperCellSize SuperCellSize;
 
             MappingDesc cellDescription;
 
-            template<uint32_t AREA>
-            void updateE()
-            {
-                return;
-            }
-
-            template<uint32_t AREA>
-            void updateBHalf()
-            {
-                return;
-            }
-
         public:
 
             NoSolver(MappingDesc cellDescription) : cellDescription(cellDescription)

src/picongpu/include/fields/MaxwellSolver/Yee/YeeSolver.hpp

@@ -57,13 +57,18 @@ class YeeSolver
 private:
     typedef MappingDesc::SuperCellSize SuperCellSize;
 
+
     FieldE* fieldE;
     FieldB* fieldB;
     MappingDesc cellDescription;
 
     template<uint32_t AREA>
     void updateE()
     {
+        /* Courant-Friedrichs-Levy-Condition for Yee Field Solver: */
+        PMACC_CASSERT_MSG(Courant_Friedrichs_Levy_condition_failure____check_your_gridConfig_param_file,
+            (SPEED_OF_LIGHT*SPEED_OF_LIGHT*DELTA_T*DELTA_T*INV_CELL2_SUM)<=1.0);
+        
         typedef SuperCellDescription<
                 SuperCellSize,
                 typename CurlB::LowerMargin,

src/picongpu/include/simulation_defines/unitless/gridConfig.unitless

@@ -47,8 +47,4 @@ namespace picongpu
     const float_X INV_CELL2_SUM = 1.0 / ( CELL_WIDTH  * CELL_WIDTH  );
 #endif
 
-    /* Courant-Friedrichs-Levy-Condition for Yee Field Solver: */
-    PMACC_CASSERT_MSG(Courant_Friedrichs_Levy_condition_failure____check_your_gridConfig_param_file,
-                     (SPEED_OF_LIGHT*SPEED_OF_LIGHT*DELTA_T*DELTA_T*INV_CELL2_SUM)<=1.0);
-
 }