Make use of the DeviceStreamManager
[alexxy/gromacs.git] / src / gromacs / ewald / tests / pmegathertest.cpp
1 /*
2  * This file is part of the GROMACS molecular simulation package.
3  *
4  * Copyright (c) 2016,2017,2018,2019,2020, by the GROMACS development team, led by
5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6  * and including many others, as listed in the AUTHORS file in the
7  * top-level source directory and at http://www.gromacs.org.
8  *
9  * GROMACS is free software; you can redistribute it and/or
10  * modify it under the terms of the GNU Lesser General Public License
11  * as published by the Free Software Foundation; either version 2.1
12  * of the License, or (at your option) any later version.
13  *
14  * GROMACS is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17  * Lesser General Public License for more details.
18  *
19  * You should have received a copy of the GNU Lesser General Public
20  * License along with GROMACS; if not, see
21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
23  *
24  * If you want to redistribute modifications to GROMACS, please
25  * consider that scientific software is very special. Version
26  * control is crucial - bugs must be traceable. We will be happy to
27  * consider code for inclusion in the official distribution, but
28  * derived work must not be called official GROMACS. Details are found
29  * in the README & COPYING files - if they are missing, get the
30  * official version at http://www.gromacs.org.
31  *
32  * To help us fund GROMACS development, we humbly ask that you cite
33  * the research papers on the package. Check out http://www.gromacs.org.
34  */
35 /*! \internal \file
36  * \brief
37  * Implements PME force gathering tests.
38  *
39  * \author Aleksei Iupinov <a.yupinov@gmail.com>
40  * \ingroup module_ewald
41  */
42
43 #include "gmxpre.h"
44
45 #include <string>
46
47 #include <gmock/gmock.h>
48
49 #include "gromacs/mdtypes/inputrec.h"
50 #include "gromacs/utility/stringutil.h"
51
52 #include "testutils/refdata.h"
53 #include "testutils/testasserts.h"
54
55 #include "pmetestcommon.h"
56 #include "testhardwarecontexts.h"
57
58 namespace gmx
59 {
60 namespace test
61 {
62 namespace
63 {
64
65 /* Valid input instances */
66
67 //! A couple of valid inputs for boxes.
68 std::vector<Matrix3x3> const c_sampleBoxes{
69     // normal box
70     Matrix3x3{ { 8.0F, 0.0F, 0.0F, 0.0F, 3.4F, 0.0F, 0.0F, 0.0F, 2.0F } },
71     // triclinic box
72     Matrix3x3{ { 7.0F, 0.0F, 0.0F, 0.0F, 4.1F, 0.0F, 3.5F, 2.0F, 12.2F } },
73 };
74
75 //! A couple of valid inputs for grid sizes
76 std::vector<IVec> const c_sampleGridSizes{ IVec{ 16, 12, 14 }, IVec{ 13, 15, 11 } };
77 //! Random charges
78 std::vector<real> const c_sampleChargesFull{ 4.95F, 3.11F, 3.97F, 1.08F, 2.09F, 1.1F,
79                                              4.13F, 3.31F, 2.8F,  5.83F, 5.09F, 6.1F,
80                                              2.86F, 0.24F, 5.76F, 5.19F, 0.72F };
81
82 //! All the input atom gridline indices
83 std::vector<IVec> const c_sampleGridLineIndicesFull{
84     IVec{ 4, 2, 6 }, IVec{ 1, 4, 10 }, IVec{ 0, 6, 6 }, IVec{ 0, 1, 4 }, IVec{ 6, 3, 0 },
85     IVec{ 7, 2, 2 }, IVec{ 8, 3, 1 },  IVec{ 4, 0, 3 }, IVec{ 0, 0, 0 }, IVec{ 8, 5, 8 },
86     IVec{ 4, 4, 2 }, IVec{ 7, 1, 7 },  IVec{ 8, 5, 5 }, IVec{ 2, 6, 5 }, IVec{ 1, 6, 2 },
87     IVec{ 7, 1, 8 }, IVec{ 3, 5, 1 },
88 };
89
90 // Spline values/derivatives below are also generated randomly, so they are bogus,
91 // but that should not affect the reproducibility, which we're after
92
93 //! A lot of bogus input spline values - should have at list (max PME order = 5) * (DIM = 3) * (total unique atom number in all test cases = 16) values
94 std::vector<real> const c_sampleSplineValuesFull{
95     0.12F, 0.81F, 0.29F, 0.22F, 0.13F, 0.19F, 0.12F, 0.8F,  0.44F, 0.38F, 0.32F, 0.36F, 0.27F,
96     0.11F, 0.17F, 0.94F, 0.07F, 0.9F,  0.98F, 0.96F, 0.07F, 0.94F, 0.77F, 0.24F, 0.84F, 0.16F,
97     0.77F, 0.57F, 0.52F, 0.27F, 0.39F, 0.45F, 0.6F,  0.59F, 0.44F, 0.91F, 0.97F, 0.43F, 0.24F,
98     0.52F, 0.73F, 0.55F, 0.99F, 0.39F, 0.97F, 0.35F, 0.1F,  0.68F, 0.19F, 0.1F,  0.77F, 0.2F,
99     0.43F, 0.69F, 0.76F, 0.32F, 0.31F, 0.94F, 0.53F, 0.6F,  0.93F, 0.57F, 0.94F, 0.88F, 0.75F,
100     0.77F, 0.91F, 0.72F, 0.07F, 0.78F, 0.09F, 0.02F, 0.48F, 0.97F, 0.89F, 0.39F, 0.48F, 0.19F,
101     0.02F, 0.92F, 0.8F,  0.41F, 0.53F, 0.32F, 0.38F, 0.58F, 0.36F, 0.46F, 0.92F, 0.91F, 0.01F,
102     0.86F, 0.54F, 0.86F, 0.94F, 0.37F, 0.35F, 0.81F, 0.89F, 0.48F, 0.34F, 0.18F, 0.11F, 0.02F,
103     0.87F, 0.95F, 0.66F, 0.67F, 0.38F, 0.45F, 0.04F, 0.94F, 0.54F, 0.76F, 0.58F, 0.83F, 0.31F,
104     0.73F, 0.71F, 0.06F, 0.35F, 0.32F, 0.35F, 0.61F, 0.27F, 0.98F, 0.83F, 0.11F, 0.3F,  0.42F,
105     0.95F, 0.69F, 0.58F, 0.29F, 0.1F,  0.68F, 0.94F, 0.62F, 0.51F, 0.47F, 0.04F, 0.47F, 0.34F,
106     0.71F, 0.52F, 0.19F, 0.69F, 0.5F,  0.59F, 0.05F, 0.74F, 0.11F, 0.4F,  0.81F, 0.24F, 0.53F,
107     0.71F, 0.07F, 0.17F, 0.41F, 0.23F, 0.78F, 0.27F, 0.1F,  0.71F, 0.36F, 0.67F, 0.6F,  0.94F,
108     0.69F, 0.19F, 0.58F, 0.68F, 0.5F,  0.62F, 0.38F, 0.29F, 0.44F, 0.04F, 0.89F, 0.0F,  0.76F,
109     0.22F, 0.16F, 0.08F, 0.62F, 0.51F, 0.62F, 0.83F, 0.72F, 0.96F, 0.99F, 0.4F,  0.79F, 0.83F,
110     0.21F, 0.43F, 0.32F, 0.44F, 0.72F, 0.21F, 0.4F,  0.93F, 0.07F, 0.11F, 0.41F, 0.24F, 0.04F,
111     0.36F, 0.15F, 0.92F, 0.08F, 0.99F, 0.35F, 0.42F, 0.7F,  0.17F, 0.39F, 0.69F, 0.0F,  0.86F,
112     0.89F, 0.59F, 0.81F, 0.77F, 0.15F, 0.89F, 0.17F, 0.76F, 0.67F, 0.58F, 0.78F, 0.26F, 0.19F,
113     0.69F, 0.18F, 0.46F, 0.6F,  0.69F, 0.23F, 0.34F, 0.3F,  0.64F, 0.34F, 0.6F,  0.99F, 0.69F,
114     0.57F, 0.75F, 0.07F, 0.36F, 0.75F, 0.81F, 0.8F,  0.42F, 0.09F, 0.94F, 0.66F, 0.35F, 0.67F,
115     0.34F, 0.66F, 0.02F, 0.47F, 0.78F, 0.21F, 0.02F, 0.18F, 0.42F, 0.2F,  0.46F, 0.34F, 0.4F,
116     0.46F, 0.96F, 0.86F, 0.25F, 0.25F, 0.22F, 0.37F, 0.59F, 0.19F, 0.45F, 0.61F, 0.04F, 0.71F,
117     0.77F, 0.51F, 0.77F, 0.15F, 0.78F, 0.36F, 0.62F, 0.24F, 0.86F, 0.2F,  0.77F, 0.08F, 0.09F,
118     0.3F,  0.0F,  0.6F,  0.99F, 0.69F,
119 };
120
121 //! A lot of bogus input spline derivatives - should have at list (max PME order = 5) * (DIM = 3) * (total unique atom number in all test cases = 16) values
122 std::vector<real> const c_sampleSplineDerivativesFull{
123     0.82F, 0.88F, 0.83F, 0.11F, 0.93F, 0.32F, 0.71F, 0.37F, 0.69F, 0.88F, 0.11F, 0.38F, 0.25F,
124     0.5F,  0.36F, 0.81F, 0.78F, 0.31F, 0.66F, 0.32F, 0.27F, 0.35F, 0.53F, 0.83F, 0.08F, 0.08F,
125     0.94F, 0.71F, 0.65F, 0.24F, 0.13F, 0.01F, 0.33F, 0.65F, 0.24F, 0.53F, 0.45F, 0.84F, 0.33F,
126     0.97F, 0.31F, 0.7F,  0.03F, 0.31F, 0.41F, 0.76F, 0.12F, 0.3F,  0.57F, 0.65F, 0.87F, 0.99F,
127     0.42F, 0.97F, 0.32F, 0.39F, 0.73F, 0.23F, 0.03F, 0.67F, 0.97F, 0.57F, 0.42F, 0.38F, 0.54F,
128     0.17F, 0.53F, 0.54F, 0.18F, 0.8F,  0.76F, 0.13F, 0.29F, 0.83F, 0.77F, 0.56F, 0.4F,  0.87F,
129     0.36F, 0.18F, 0.59F, 0.04F, 0.05F, 0.61F, 0.26F, 0.91F, 0.62F, 0.16F, 0.89F, 0.23F, 0.26F,
130     0.59F, 0.33F, 0.2F,  0.49F, 0.41F, 0.25F, 0.4F,  0.16F, 0.83F, 0.44F, 0.82F, 0.21F, 0.95F,
131     0.14F, 0.8F,  0.37F, 0.31F, 0.41F, 0.53F, 0.15F, 0.85F, 0.78F, 0.17F, 0.92F, 0.03F, 0.13F,
132     0.2F,  0.03F, 0.33F, 0.87F, 0.38F, 0,     0.08F, 0.79F, 0.36F, 0.53F, 0.05F, 0.07F, 0.94F,
133     0.23F, 0.85F, 0.13F, 0.27F, 0.23F, 0.22F, 0.26F, 0.38F, 0.15F, 0.48F, 0.18F, 0.33F, 0.23F,
134     0.62F, 0.1F,  0.36F, 0.99F, 0.07F, 0.02F, 0.04F, 0.09F, 0.29F, 0.52F, 0.29F, 0.83F, 0.97F,
135     0.61F, 0.81F, 0.49F, 0.56F, 0.08F, 0.09F, 0.03F, 0.65F, 0.46F, 0.1F,  0.06F, 0.06F, 0.39F,
136     0.29F, 0.04F, 0.03F, 0.1F,  0.83F, 0.94F, 0.59F, 0.97F, 0.82F, 0.2F,  0.66F, 0.23F, 0.11F,
137     0.03F, 0.16F, 0.27F, 0.53F, 0.94F, 0.46F, 0.43F, 0.29F, 0.97F, 0.64F, 0.46F, 0.37F, 0.43F,
138     0.48F, 0.37F, 0.93F, 0.5F,  0.2F,  0.92F, 0.09F, 0.74F, 0.55F, 0.44F, 0.05F, 0.13F, 0.17F,
139     0.79F, 0.44F, 0.11F, 0.6F,  0.64F, 0.05F, 0.96F, 0.3F,  0.45F, 0.47F, 0.42F, 0.74F, 0.91F,
140     0.06F, 0.89F, 0.24F, 0.26F, 0.68F, 0.4F,  0.88F, 0.5F,  0.65F, 0.48F, 0.15F, 0.0F,  0.41F,
141     0.67F, 0.4F,  0.31F, 0.73F, 0.77F, 0.36F, 0.26F, 0.74F, 0.46F, 0.56F, 0.78F, 0.92F, 0.32F,
142     0.9F,  0.06F, 0.55F, 0.6F,  0.13F, 0.38F, 0.93F, 0.5F,  0.92F, 0.96F, 0.82F, 0.0F,  0.04F,
143     0.9F,  0.55F, 0.97F, 1.0F,  0.23F, 0.46F, 0.52F, 0.49F, 0.0F,  0.32F, 0.16F, 0.4F,  0.62F,
144     0.36F, 0.03F, 0.63F, 0.16F, 0.58F, 0.97F, 0.03F, 0.44F, 0.07F, 0.22F, 0.75F, 0.32F, 0.61F,
145     0.94F, 0.33F, 0.7F,  0.57F, 0.5F,  0.84F, 0.7F,  0.47F, 0.18F, 0.09F, 0.25F, 0.77F, 0.94F,
146     0.85F, 0.09F, 0.83F, 0.02F, 0.91F,
147 };
148
149 //! 2 c_sample grids - only non-zero values have to be listed
150 std::vector<SparseRealGridValuesInput> const c_sampleGrids{ SparseRealGridValuesInput{
151                                                                     { IVec{ 0, 0, 0 }, 3.5F },
152                                                                     { IVec{ 7, 0, 0 }, -2.5F },
153                                                                     { IVec{ 3, 5, 7 }, -0.006F },
154                                                                     { IVec{ 1, 6, 7 }, -2.76F },
155                                                                     { IVec{ 3, 1, 2 }, 0.6F },
156                                                                     { IVec{ 6, 2, 4 }, 7.1F },
157                                                                     { IVec{ 4, 5, 6 }, 4.1F },
158                                                                     { IVec{ 4, 4, 6 }, -3.7F },
159                                                             },
160                                                             SparseRealGridValuesInput{
161                                                                     { IVec{ 0, 4, 0 }, 6.F },
162                                                                     { IVec{ 4, 2, 7 }, 13.76F },
163                                                                     { IVec{ 0, 6, 7 }, 3.6F },
164                                                                     { IVec{ 3, 2, 8 }, 0.61F },
165                                                                     { IVec{ 5, 4, 3 }, 2.1F },
166                                                                     { IVec{ 2, 5, 10 }, 3.6F },
167                                                                     { IVec{ 5, 3, 6 }, 2.1F },
168                                                                     { IVec{ 6, 6, 6 }, 6.1F },
169                                                             } };
170
171 //! Input forces for reduction
172 std::vector<RVec> const c_sampleForcesFull{
173     RVec{ 0.02F, 0.87F, 0.95F }, RVec{ 0.66F, 0.67F, 0.38F }, RVec{ 0.45F, 0.04F, 0.94F },
174     RVec{ 0.54F, 0.76F, 0.58F }, RVec{ 0.83F, 0.31F, 0.73F }, RVec{ 0.71F, 0.06F, 0.35F },
175     RVec{ 0.32F, 0.35F, 0.61F }, RVec{ 0.27F, 0.98F, 0.83F }, RVec{ 0.11F, 0.3F, 0.42F },
176     RVec{ 0.95F, 0.69F, 0.58F }, RVec{ 0.29F, 0.1F, 0.68F },  RVec{ 0.94F, 0.62F, 0.51F },
177     RVec{ 0.47F, 0.04F, 0.47F }, RVec{ 0.34F, 0.71F, 0.52F }
178 };
179
180 //! PME orders to test
181 std::vector<int> const pmeOrders{ 3, 4, 5 };
182 //! Atom counts to test
183 std::vector<size_t> const atomCounts{ 1, 2, 13 };
184
185 /* Helper structures for test input */
186
187 //! A structure for all the spline data which depends in size both on the PME order and atom count
188 struct AtomAndPmeOrderSizedData
189 {
190     //! Spline values
191     SplineParamsVector splineValues;
192     //! Spline derivatives
193     SplineParamsVector splineDerivatives;
194 };
195
196 //! A structure for all the input atom data which depends in size on atom count - including a range of spline data for different PME orders
197 struct AtomSizedData
198 {
199     //! Gridline indices
200     GridLineIndicesVector gridLineIndices;
201     //! Charges
202     ChargesVector charges;
203     //! Coordinates
204     CoordinatesVector coordinates;
205     //! Spline data for different orders
206     std::map<int, AtomAndPmeOrderSizedData> splineDataByPmeOrder;
207 };
208
209 //! A range of the test input data sets, uniquely identified by the atom count
210 typedef std::map<size_t, AtomSizedData> InputDataByAtomCount;
211
212 /*! \brief Convenience typedef of the test input parameters - unit cell box, PME interpolation
213  * order, grid dimensions, grid values, overwriting/reducing the input forces, atom count.
214  *
215  * The rest of the atom-related input data - gridline indices, spline theta values, spline dtheta
216  * values, atom charges - is looked up in the inputAtomDataSets_ test fixture variable.
217  */
218 typedef std::tuple<Matrix3x3, int, IVec, SparseRealGridValuesInput, size_t> GatherInputParameters;
219
220 //! Test fixture
221 class PmeGatherTest : public ::testing::TestWithParam<GatherInputParameters>
222 {
223 private:
224     //! Storage of all the input atom datasets
225     static InputDataByAtomCount s_inputAtomDataSets_;
226
227 public:
228     PmeGatherTest() = default;
229     //! Sets the input atom data references once
230     static void SetUpTestCase()
231     {
232         size_t start = 0;
233         for (auto atomCount : atomCounts)
234         {
235             AtomSizedData atomData;
236             atomData.gridLineIndices =
237                     GridLineIndicesVector(c_sampleGridLineIndicesFull).subArray(start, atomCount);
238             atomData.charges = ChargesVector(c_sampleChargesFull).subArray(start, atomCount);
239             start += atomCount;
240             atomData.coordinates.resize(atomCount, RVec{ 1e6, 1e7, -1e8 });
241             /* The coordinates are intentionally bogus in this test - only the size matters; the gridline indices are fed directly as inputs */
242             for (auto pmeOrder : pmeOrders)
243             {
244                 AtomAndPmeOrderSizedData splineData;
245                 const size_t             dimSize = atomCount * pmeOrder;
246                 for (int dimIndex = 0; dimIndex < DIM; dimIndex++)
247                 {
248                     splineData.splineValues[dimIndex] =
249                             SplineParamsDimVector(c_sampleSplineValuesFull).subArray(dimIndex * dimSize, dimSize);
250                     splineData.splineDerivatives[dimIndex] =
251                             SplineParamsDimVector(c_sampleSplineDerivativesFull).subArray(dimIndex * dimSize, dimSize);
252                 }
253                 atomData.splineDataByPmeOrder[pmeOrder] = splineData;
254             }
255             s_inputAtomDataSets_[atomCount] = atomData;
256         }
257     }
258     //! The test
259     void runTest()
260     {
261         /* Getting the input */
262         Matrix3x3                 box;
263         int                       pmeOrder;
264         IVec                      gridSize;
265         size_t                    atomCount;
266         SparseRealGridValuesInput nonZeroGridValues;
267         std::tie(box, pmeOrder, gridSize, nonZeroGridValues, atomCount) = GetParam();
268         auto inputAtomData       = s_inputAtomDataSets_[atomCount];
269         auto inputAtomSplineData = inputAtomData.splineDataByPmeOrder[pmeOrder];
270
271         /* Storing the input where it's needed, running the test */
272         t_inputrec inputRec;
273         inputRec.nkx         = gridSize[XX];
274         inputRec.nky         = gridSize[YY];
275         inputRec.nkz         = gridSize[ZZ];
276         inputRec.pme_order   = pmeOrder;
277         inputRec.coulombtype = eelPME;
278         inputRec.epsilon_r   = 1.0;
279
280         TestReferenceData refData;
281         for (const auto& context : getPmeTestEnv()->getHardwareContexts())
282         {
283             CodePath   codePath = context->codePath();
284             const bool supportedInput =
285                     pmeSupportsInputForMode(*getPmeTestEnv()->hwinfo(), &inputRec, codePath);
286             if (!supportedInput)
287             {
288                 /* Testing the failure for the unsupported input */
289                 EXPECT_THROW_GMX(pmeInitWrapper(&inputRec, codePath, nullptr, nullptr, nullptr, box),
290                                  NotImplementedError);
291                 continue;
292             }
293
294             /* Describing the test uniquely */
295             SCOPED_TRACE(
296                     formatString("Testing force gathering with %s %sfor PME grid size %d %d %d"
297                                  ", order %d, %zu atoms",
298                                  codePathToString(codePath), context->description().c_str(),
299                                  gridSize[XX], gridSize[YY], gridSize[ZZ], pmeOrder, atomCount));
300
301             PmeSafePointer pmeSafe =
302                     pmeInitWrapper(&inputRec, codePath, context->deviceContext(),
303                                    context->deviceStream(), context->pmeGpuProgram(), box);
304             std::unique_ptr<StatePropagatorDataGpu> stateGpu =
305                     (codePath == CodePath::GPU)
306                             ? makeStatePropagatorDataGpu(*pmeSafe.get(), context->deviceContext(),
307                                                          context->deviceStream())
308                             : nullptr;
309
310             pmeInitAtoms(pmeSafe.get(), stateGpu.get(), codePath, inputAtomData.coordinates,
311                          inputAtomData.charges);
312
313             /* Setting some more inputs */
314             pmeSetRealGrid(pmeSafe.get(), codePath, nonZeroGridValues);
315             pmeSetGridLineIndices(pmeSafe.get(), codePath, inputAtomData.gridLineIndices);
316             for (int dimIndex = 0; dimIndex < DIM; dimIndex++)
317             {
318                 pmeSetSplineData(pmeSafe.get(), codePath, inputAtomSplineData.splineValues[dimIndex],
319                                  PmeSplineDataType::Values, dimIndex);
320                 pmeSetSplineData(pmeSafe.get(), codePath, inputAtomSplineData.splineDerivatives[dimIndex],
321                                  PmeSplineDataType::Derivatives, dimIndex);
322             }
323
324             /* Explicitly copying the sample forces to be able to modify them */
325             auto inputForcesFull(c_sampleForcesFull);
326             GMX_RELEASE_ASSERT(inputForcesFull.size() >= atomCount, "Bad input forces size");
327             auto forces = ForcesVector(inputForcesFull).subArray(0, atomCount);
328
329             /* Running the force gathering itself */
330             pmePerformGather(pmeSafe.get(), codePath, forces);
331             pmeFinalizeTest(pmeSafe.get(), codePath);
332
333             /* Check the output forces correctness */
334             TestReferenceChecker forceChecker(refData.rootChecker());
335             const auto           ulpTolerance = 3 * pmeOrder;
336             forceChecker.setDefaultTolerance(relativeToleranceAsUlp(1.0, ulpTolerance));
337             forceChecker.checkSequence(forces.begin(), forces.end(), "Forces");
338         }
339     }
340 };
341
342 // An instance of static atom data
343 InputDataByAtomCount PmeGatherTest::s_inputAtomDataSets_;
344
345 //! Test for PME force gathering
346 TEST_P(PmeGatherTest, ReproducesOutputs)
347 {
348     EXPECT_NO_THROW_GMX(runTest());
349 }
350
351 //! Instantiation of the PME gathering test
352 INSTANTIATE_TEST_CASE_P(SaneInput,
353                         PmeGatherTest,
354                         ::testing::Combine(::testing::ValuesIn(c_sampleBoxes),
355                                            ::testing::ValuesIn(pmeOrders),
356                                            ::testing::ValuesIn(c_sampleGridSizes),
357                                            ::testing::ValuesIn(c_sampleGrids),
358                                            ::testing::ValuesIn(atomCounts)));
359
360 } // namespace
361 } // namespace test
362 } // namespace gmx