Intrepid2
Intrepid2_TestUtils.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Intrepid2 Package
5// Copyright (2007) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Kyungjoo Kim (kyukim@sandia.gov),
38// Mauro Perego (mperego@sandia.gov), or
39// Nate Roberts (nvrober@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
49#ifndef Intrepid2_TestUtils_h
50#define Intrepid2_TestUtils_h
51
52#include "Kokkos_Core.hpp"
53#include "Kokkos_DynRankView.hpp"
54
55#include "Intrepid2_Basis.hpp"
60#include "Intrepid2_Sacado.hpp" // Sacado includes, guarded by the appropriate preprocessor variable
61#include "Intrepid2_Utils.hpp"
62
63#include "Teuchos_UnitTestHarness.hpp"
64
65namespace Intrepid2
66{
69
71#if defined(KOKKOS_ENABLE_CUDA)
72 using DefaultTestDeviceType = Kokkos::Device<Kokkos::Cuda,Kokkos::CudaSpace>;
73#elif defined(KOKKOS_ENABLE_HIP)
74 using DefaultTestDeviceType = Kokkos::Device<Kokkos::Experimental::HIP,Kokkos::Experimental::HIPSpace>;
75#else
76 using DefaultTestDeviceType = typename Kokkos::DefaultExecutionSpace::device_type;
77#endif
78
80 template <class Scalar>
81 const typename Teuchos::ScalarTraits<Scalar>::magnitudeType
83 {
84 using ST = Teuchos::ScalarTraits<Scalar>;
85 return ST::magnitude(Teuchos::RelErrSmallNumber<ST::hasMachineParameters,Scalar>::smallNumber());
86 }
87
89 template <class Scalar1, class Scalar2>
90 KOKKOS_INLINE_FUNCTION
91 bool
92 relErrMeetsTol( const Scalar1 &s1, const Scalar2 &s2, const typename Teuchos::ScalarTraits< typename std::common_type<Scalar1,Scalar2>::type >::magnitudeType &smallNumber, const double &tol )
93 {
94 using Scalar = typename std::common_type<Scalar1,Scalar2>::type;
95 const Scalar s1Abs = fabs(s1);
96 const Scalar s2Abs = fabs(s2);
97 const Scalar maxAbs = (s1Abs > s2Abs) ? s1Abs : s2Abs;
98 const Scalar relErr = fabs( s1 - s2 ) / ( smallNumber + maxAbs );
99 return relErr < tol;
100 }
101
102 template <class Scalar1, class Scalar2>
103 KOKKOS_INLINE_FUNCTION
104 bool
105 errMeetsAbsAndRelTol( const Scalar1 &s1, const Scalar2 &s2, const double &relTol, const double &absTol )
106 {
107 return fabs( s1 - s2 ) <= absTol + fabs(s1) * relTol;
108 }
109
110 static const double TEST_TOLERANCE_TIGHT = 1.e2 * std::numeric_limits<double>::epsilon();
111
112 // we use DynRankView for both input points and values
113 template<typename ScalarType, typename DeviceType>
114 using ViewType = Kokkos::DynRankView<ScalarType,DeviceType>;
115
116 template<typename ScalarType, typename DeviceType>
117 using FixedRankViewType = Kokkos::View<ScalarType,DeviceType>;
118
119 template<typename ScalarType>
120 KOKKOS_INLINE_FUNCTION bool valuesAreSmall(const ScalarType &a, const ScalarType &b, const double &epsilon)
121 {
122 using std::abs;
123 return (abs(a) < epsilon) && (abs(b) < epsilon);
124 }
125
126 inline bool approximatelyEqual(double a, double b, double epsilon)
127 {
128 const double larger_magnitude = (std::abs(a) < std::abs(b) ? std::abs(b) : std::abs(a));
129 return std::abs(a - b) <= larger_magnitude * epsilon;
130 }
131
132 inline bool essentiallyEqual(double a, double b, double epsilon)
133 {
134 const double smaller_magnitude = (std::abs(a) > std::abs(b) ? std::abs(b) : std::abs(a));
135 return std::abs(a - b) <= smaller_magnitude * epsilon;
136 }
137
138 // conversion from the ref element [0,1] (used by ESEAS) to the ref element [-1,1] (used by Intrepid2)
139 KOKKOS_INLINE_FUNCTION double fromZeroOne(double x_zero_one)
140 {
141 return x_zero_one * 2.0 - 1.0;
142 }
143
144 // conversion from the ref element [-1,1] (used by Intrepid2) to the ref element [0,1] (used by ESEAS)
145 KOKKOS_INLINE_FUNCTION double toZeroOne(double x_minus_one_one)
146 {
147 return (x_minus_one_one + 1.0) / 2.0;
148 }
149
150 // conversion from the ref element [0,1] (used by ESEAS) to the ref element [-1,1] (used by Intrepid2)
151 KOKKOS_INLINE_FUNCTION double fromZeroOne_dx(double dx_zero_one)
152 {
153 return dx_zero_one / 2.0;
154 }
155
156 // conversion from the ref element [-1,1] (used by Intrepid2) to the ref element [0,1] (used by ESEAS)
157 KOKKOS_INLINE_FUNCTION double toZeroOne_dx(double dx_minus_one_one)
158 {
159 return dx_minus_one_one * 2.0;
160 }
161
162 template<class DeviceViewType>
163 typename DeviceViewType::HostMirror getHostCopy(const DeviceViewType &deviceView)
164 {
165 typename DeviceViewType::HostMirror hostView = Kokkos::create_mirror(deviceView);
166 Kokkos::deep_copy(hostView, deviceView);
167 return hostView;
168 }
169
170 template<class BasisFamily>
171 inline Teuchos::RCP< Intrepid2::Basis<DefaultTestDeviceType,double,double> > getBasisUsingFamily(shards::CellTopology cellTopo, Intrepid2::EFunctionSpace fs,
172 int polyOrder_x, int polyOrder_y=-1, int polyOrder_z = -1)
173 {
174 using BasisPtr = typename BasisFamily::BasisPtr;
175
176 BasisPtr basis;
177 using namespace Intrepid2;
178
179 if (cellTopo.getBaseKey() == shards::Line<>::key)
180 {
181 basis = getLineBasis<BasisFamily>(fs,polyOrder_x);
182 }
183 else if (cellTopo.getBaseKey() == shards::Quadrilateral<>::key)
184 {
185 INTREPID2_TEST_FOR_EXCEPTION(polyOrder_y < 0, std::invalid_argument, "polyOrder_y must be specified");
186 basis = getQuadrilateralBasis<BasisFamily>(fs,polyOrder_x,polyOrder_y);
187 }
188 else if (cellTopo.getBaseKey() == shards::Triangle<>::key)
189 {
190 basis = getTriangleBasis<BasisFamily>(fs,polyOrder_x);
191 }
192 else if (cellTopo.getBaseKey() == shards::Hexahedron<>::key)
193 {
194 INTREPID2_TEST_FOR_EXCEPTION(polyOrder_y < 0, std::invalid_argument, "polyOrder_y must be specified");
195 INTREPID2_TEST_FOR_EXCEPTION(polyOrder_z < 0, std::invalid_argument, "polyOrder_z must be specified");
196 basis = getHexahedronBasis<BasisFamily>(fs,polyOrder_x,polyOrder_y,polyOrder_z);
197 }
198 else if (cellTopo.getBaseKey() == shards::Tetrahedron<>::key)
199 {
200 basis = getTetrahedronBasis<BasisFamily>(fs, polyOrder_x);
201 }
202 else if (cellTopo.getBaseKey() == shards::Wedge<>::key)
203 {
204 INTREPID2_TEST_FOR_EXCEPTION(polyOrder_y < 0, std::invalid_argument, "polyOrder_y must be specified");
205 basis = getWedgeBasis<BasisFamily>(fs,polyOrder_x,polyOrder_y);
206 }
207 else
208 {
209 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported cell topology");
210 }
211 return basis;
212 }
213
214 template<bool defineVertexFunctions>
215 inline Teuchos::RCP< Intrepid2::Basis<DefaultTestDeviceType,double,double> > getHierarchicalBasis(shards::CellTopology cellTopo, Intrepid2::EFunctionSpace fs,
216 int polyOrder_x, int polyOrder_y=-1, int polyOrder_z = -1)
217 {
218 using DeviceType = DefaultTestDeviceType;
219 using Scalar = double;
220 using namespace Intrepid2;
221
226
228
229 return getBasisUsingFamily<BasisFamily>(cellTopo, fs, polyOrder_x, polyOrder_y, polyOrder_z);
230 }
231
232 template<typename ValueType, typename DeviceType, class ... DimArgs>
233 inline ViewType<ValueType,DeviceType> getView(const std::string &label, DimArgs... dims)
234 {
235 const bool allocateFadStorage = !std::is_pod<ValueType>::value;
236 if (!allocateFadStorage)
237 {
238 return ViewType<ValueType,DeviceType>(label,dims...);
239 }
240 else
241 {
242 return ViewType<ValueType,DeviceType>(label,dims...,MAX_FAD_DERIVATIVES_FOR_TESTS+1);
243 }
244 }
245
246 // this method is to allow us to switch tests over incrementally; should collapse with ViewType once everything has been switched
247 template<typename ValueType, class ... DimArgs>
248 inline FixedRankViewType< typename RankExpander<ValueType, sizeof...(DimArgs) >::value_type, DefaultTestDeviceType > getFixedRankView(const std::string &label, DimArgs... dims)
249 {
250 const bool allocateFadStorage = !std::is_pod<ValueType>::value;
251 using value_type = typename RankExpander<ValueType, sizeof...(dims) >::value_type;
252 if (!allocateFadStorage)
253 {
254 return FixedRankViewType<value_type,DefaultTestDeviceType>(label,dims...);
255 }
256 else
257 {
258 return FixedRankViewType<value_type,DefaultTestDeviceType>(label,dims...,MAX_FAD_DERIVATIVES_FOR_TESTS+1);
259 }
260 }
261
268 template <typename PointValueType, typename DeviceType>
269 inline ViewType<PointValueType,DeviceType> getInputPointsView(shards::CellTopology &cellTopo, int numPoints_1D)
270 {
271 if (cellTopo.getBaseKey() == shards::Wedge<>::key)
272 {
273 shards::CellTopology lineTopo = shards::CellTopology(shards::getCellTopologyData<shards::Line<> >() );
274 shards::CellTopology triTopo = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<> >() );
275
276 const ordinal_type order = numPoints_1D - 1;
277 ordinal_type numPoints_tri = PointTools::getLatticeSize(triTopo, order);
278 ordinal_type numPoints_line = PointTools::getLatticeSize(lineTopo, order);
279 ordinal_type numPoints = numPoints_tri * numPoints_line;
280 ordinal_type spaceDim = cellTopo.getDimension();
281
282 ViewType<PointValueType,DeviceType> inputPointsTri = getView<PointValueType,DeviceType>("input points",numPoints_tri, 2);
283 ViewType<PointValueType,DeviceType> inputPointsLine = getView<PointValueType,DeviceType>("input points",numPoints_line,1);
284 PointTools::getLattice(inputPointsTri, triTopo, order, 0, POINTTYPE_EQUISPACED );
285 PointTools::getLattice(inputPointsLine, lineTopo, order, 0, POINTTYPE_EQUISPACED );
286
287 ViewType<PointValueType,DeviceType> inputPoints = getView<PointValueType,DeviceType>("input points",numPoints,spaceDim);
288
289 using ExecutionSpace = typename ViewType<PointValueType,DeviceType>::execution_space;
290
291 Kokkos::RangePolicy < ExecutionSpace > policy(0,numPoints_tri);
292 Kokkos::parallel_for( policy,
293 KOKKOS_LAMBDA (const ordinal_type &triPointOrdinal )
294 {
295 ordinal_type pointOrdinal = triPointOrdinal * numPoints_line;
296 for (ordinal_type linePointOrdinal=0; linePointOrdinal<numPoints_line; linePointOrdinal++)
297 {
298 inputPoints(pointOrdinal,0) = inputPointsTri( triPointOrdinal,0);
299 inputPoints(pointOrdinal,1) = inputPointsTri( triPointOrdinal,1);
300 inputPoints(pointOrdinal,2) = inputPointsLine(linePointOrdinal,0);
301 pointOrdinal++;
302 }
303 }
304 );
305
306 return inputPoints;
307 }
308 else
309 {
310 const ordinal_type order = numPoints_1D - 1;
311 ordinal_type numPoints = PointTools::getLatticeSize(cellTopo, order);
312 ordinal_type spaceDim = cellTopo.getDimension();
313
314 ViewType<PointValueType,DeviceType> inputPoints = getView<PointValueType,DeviceType>("input points",numPoints,spaceDim);
315 PointTools::getLattice(inputPoints, cellTopo, order, 0, POINTTYPE_EQUISPACED );
316
317 return inputPoints;
318 }
319 }
320
321 template<typename OutputValueType, typename DeviceType>
322 inline ViewType<OutputValueType,DeviceType> getOutputView(Intrepid2::EFunctionSpace fs, Intrepid2::EOperator op, int basisCardinality, int numPoints, int spaceDim)
323 {
324 switch (fs) {
325 case Intrepid2::FUNCTION_SPACE_HGRAD:
326 switch (op) {
327 case Intrepid2::OPERATOR_VALUE:
328 return getView<OutputValueType,DeviceType>("H^1 value output",basisCardinality,numPoints);
329 case Intrepid2::OPERATOR_GRAD:
330 return getView<OutputValueType,DeviceType>("H^1 derivative output",basisCardinality,numPoints,spaceDim);
331 case Intrepid2::OPERATOR_D1:
332 case Intrepid2::OPERATOR_D2:
333 case Intrepid2::OPERATOR_D3:
334 case Intrepid2::OPERATOR_D4:
335 case Intrepid2::OPERATOR_D5:
336 case Intrepid2::OPERATOR_D6:
337 case Intrepid2::OPERATOR_D7:
338 case Intrepid2::OPERATOR_D8:
339 case Intrepid2::OPERATOR_D9:
340 case Intrepid2::OPERATOR_D10:
341 {
342 const auto dkcard = getDkCardinality(op, spaceDim);
343 return getView<OutputValueType,DeviceType>("H^1 derivative output",basisCardinality,numPoints,dkcard);
344 }
345 default:
346 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported op/fs combination");
347 }
348 case Intrepid2::FUNCTION_SPACE_HCURL:
349 switch (op) {
350 case Intrepid2::OPERATOR_VALUE:
351 return getView<OutputValueType,DeviceType>("H(curl) value output",basisCardinality,numPoints,spaceDim);
352 case Intrepid2::OPERATOR_CURL:
353 if (spaceDim == 2)
354 return getView<OutputValueType,DeviceType>("H(curl) derivative output",basisCardinality,numPoints);
355 else if (spaceDim == 3)
356 return getView<OutputValueType,DeviceType>("H(curl) derivative output",basisCardinality,numPoints,spaceDim);
357 default:
358 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported op/fs combination");
359 }
360 case Intrepid2::FUNCTION_SPACE_HDIV:
361 switch (op) {
362 case Intrepid2::OPERATOR_VALUE:
363 return getView<OutputValueType,DeviceType>("H(div) value output",basisCardinality,numPoints,spaceDim);
364 case Intrepid2::OPERATOR_DIV:
365 return getView<OutputValueType,DeviceType>("H(div) derivative output",basisCardinality,numPoints);
366 default:
367 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported op/fs combination");
368 }
369
370 case Intrepid2::FUNCTION_SPACE_HVOL:
371 switch (op) {
372 case Intrepid2::OPERATOR_VALUE:
373 return getView<OutputValueType,DeviceType>("H(vol) value output",basisCardinality,numPoints);
374 case Intrepid2::OPERATOR_D1:
375 case Intrepid2::OPERATOR_D2:
376 case Intrepid2::OPERATOR_D3:
377 case Intrepid2::OPERATOR_D4:
378 case Intrepid2::OPERATOR_D5:
379 case Intrepid2::OPERATOR_D6:
380 case Intrepid2::OPERATOR_D7:
381 case Intrepid2::OPERATOR_D8:
382 case Intrepid2::OPERATOR_D9:
383 case Intrepid2::OPERATOR_D10:
384 {
385 const auto dkcard = getDkCardinality(op, spaceDim);
386 return getView<OutputValueType,DeviceType>("H(vol) derivative output",basisCardinality,numPoints,dkcard);
387 }
388 default:
389 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported op/fs combination");
390 }
391 default:
392 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported op/fs combination");
393 }
394 }
395
396 // ! This returns a vector whose entries are vector<int>s containing 1-3 polynomial orders from 1 up to and including those specified
397 // ! Intended for testing bases that support anisotropic polynomial degree, such as the hierarchical bases
398 inline std::vector< std::vector<int> > getBasisTestCasesUpToDegree(int spaceDim, int minDegree, int polyOrder_x, int polyOrder_y=-1, int polyOrder_z=-1)
399 {
400 std::vector<int> degrees(spaceDim);
401 degrees[0] = polyOrder_x;
402 if (spaceDim > 1) degrees[1] = polyOrder_y;
403 if (spaceDim > 2) degrees[2] = polyOrder_z;
404
405 int numCases = degrees[0];
406 for (unsigned d=1; d<degrees.size(); d++)
407 {
408 INTREPID2_TEST_FOR_EXCEPTION(degrees[d] < minDegree, std::invalid_argument, "Unsupported degree/minDegree combination");
409 numCases = numCases * (degrees[d] + 1 - minDegree);
410 }
411 std::vector< std::vector<int> > subBasisDegreeTestCases(numCases);
412 for (int caseOrdinal=0; caseOrdinal<numCases; caseOrdinal++)
413 {
414 std::vector<int> subBasisDegrees(degrees.size());
415 int caseRemainder = caseOrdinal;
416 for (int d=degrees.size()-1; d>=0; d--)
417 {
418 int subBasisDegree = caseRemainder % (degrees[d] + 1 - minDegree);
419 caseRemainder = caseRemainder / (degrees[d] + 1 - minDegree);
420 subBasisDegrees[d] = subBasisDegree + minDegree;
421 }
422 subBasisDegreeTestCases[caseOrdinal] = subBasisDegrees;
423 }
424 return subBasisDegreeTestCases;
425 }
426
428 template<class Functor, class Scalar, int rank>
429 typename ViewType<Scalar,DefaultTestDeviceType>::HostMirror copyFunctorToHostView(const Functor &deviceFunctor)
430 {
431 INTREPID2_TEST_FOR_EXCEPTION(rank != getFunctorRank(deviceFunctor), std::invalid_argument, "functor rank must match the template argument");
432
433 using DeviceType = DefaultTestDeviceType;
434 ViewType<Scalar,DeviceType> view;
435 const std::string label = "functor copy";
436 const auto &f = deviceFunctor;
437 switch (rank)
438 {
439 case 0:
440 view = getView<Scalar,DeviceType>(label);
441 break;
442 case 1:
443 view = getView<Scalar,DeviceType>(label, f.extent_int(0));
444 break;
445 case 2:
446 view = getView<Scalar,DeviceType>(label, f.extent_int(0), f.extent_int(1));
447 break;
448 case 3:
449 view = getView<Scalar,DeviceType>(label, f.extent_int(0), f.extent_int(1), f.extent_int(2));
450 break;
451 case 4:
452 view = getView<Scalar,DeviceType>(label, f.extent_int(0), f.extent_int(1), f.extent_int(2), f.extent_int(3));
453 break;
454 case 5:
455 view = getView<Scalar,DeviceType>(label, f.extent_int(0), f.extent_int(1), f.extent_int(2), f.extent_int(3), f.extent_int(4));
456 break;
457 case 6:
458 view = getView<Scalar,DeviceType>(label, f.extent_int(0), f.extent_int(1), f.extent_int(2), f.extent_int(3), f.extent_int(4), f.extent_int(5));
459 break;
460 case 7:
461 view = getView<Scalar,DeviceType>(label, f.extent_int(0), f.extent_int(1), f.extent_int(2), f.extent_int(3), f.extent_int(4), f.extent_int(5), f.extent_int(6));
462 break;
463 default:
464 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported functor rank");
465 }
466
467 int entryCount = view.size();
468
469 using ExecutionSpace = typename ViewType<Scalar,DeviceType>::execution_space;
470
471 using ViewIteratorScalar = Intrepid2::ViewIterator<ViewType<Scalar,DeviceType>, Scalar>;
472 using FunctorIteratorScalar = FunctorIterator<Functor, Scalar, rank>;
473
474 Kokkos::RangePolicy < ExecutionSpace > policy(0,entryCount);
475 Kokkos::parallel_for( policy,
476 KOKKOS_LAMBDA (const int &enumerationIndex )
477 {
478 ViewIteratorScalar vi(view);
479 FunctorIteratorScalar fi(f);
480
481 vi.setEnumerationIndex(enumerationIndex);
482 fi.setEnumerationIndex(enumerationIndex);
483
484 vi.set(fi.get());
485 }
486 );
487
488 auto hostView = Kokkos::create_mirror_view(view);
489 Kokkos::deep_copy(hostView, view);
490 return hostView;
491 }
492
493 template<class FunctorType, typename Scalar, int rank>
494 void printFunctor(const FunctorType &functor, std::ostream &out, const std::string &functorName = "")
495 {
496 auto functorHostCopy = copyFunctorToHostView<FunctorType, Scalar, rank>(functor);
497
498 std::string name = (functorName == "") ? "Functor" : functorName;
499
500 out << "\n******** " << name << " (rank " << rank << ") ********\n";
501 out << "dimensions: (";
502 for (int r=0; r<rank; r++)
503 {
504 out << functor.extent_int(r);
505 if (r<rank-1) out << ",";
506 }
507 out << ")\n";
508
509 ViewIterator<decltype(functorHostCopy),Scalar> vi(functorHostCopy);
510
511 bool moreEntries = true;
512 while (moreEntries)
513 {
514 Scalar value = vi.get();
515
516 auto location = vi.getLocation();
517 out << functorName << "(";
518 for (ordinal_type i=0; i<rank; i++)
519 {
520 out << location[i];
521 if (i<rank-1)
522 {
523 out << ",";
524 }
525 }
526 out << ") " << value << std::endl;
527
528 moreEntries = (vi.increment() != -1);
529 }
530 out << "\n";
531 }
532
533 template<class FunctorType>
534 void printFunctor1(const FunctorType &functor, std::ostream &out, const std::string &functorName = "")
535 {
536 using Scalar = typename std::remove_reference<decltype(functor(0))>::type;
537 printFunctor<FunctorType, Scalar, 1>(functor, out, functorName);
538 }
539
540 template<class FunctorType>
541 void printFunctor2(const FunctorType &functor, std::ostream &out, const std::string &functorName = "")
542 {
543 using Scalar = typename std::remove_const<typename std::remove_reference<decltype(functor(0,0))>::type>::type;
544 printFunctor<FunctorType, Scalar, 2>(functor, out, functorName);
545 }
546
547 template<class FunctorType>
548 void printFunctor3(const FunctorType &functor, std::ostream &out, const std::string &functorName = "")
549 {
550 using Scalar = typename std::remove_const<typename std::remove_reference<decltype(functor(0,0,0))>::type>::type;
551 printFunctor<FunctorType, Scalar, 3>(functor, out, functorName);
552 }
553
554 template<class FunctorType>
555 void printFunctor4(const FunctorType &functor, std::ostream &out, const std::string &functorName = "")
556 {
557 using Scalar = typename std::remove_const<typename std::remove_reference<decltype(functor(0,0,0,0))>::type>::type;
558 printFunctor<FunctorType, Scalar, 4>(functor, out, functorName);
559 }
560
561 template<class FunctorType>
562 void printFunctor5(const FunctorType &functor, std::ostream &out, const std::string &functorName = "")
563 {
564 using Scalar = typename std::remove_const<typename std::remove_reference<decltype(functor(0,0,0,0,0))>::type>::type;
565 printFunctor<FunctorType, Scalar, 5>(functor, out, functorName);
566 }
567
568 template<class FunctorType>
569 void printFunctor6(const FunctorType &functor, std::ostream &out, const std::string &functorName = "")
570 {
571 using Scalar = typename std::remove_const<typename std::remove_reference<decltype(functor(0,0,0,0,0,0))>::type>::type;
572 printFunctor<FunctorType, Scalar, 6>(functor, out, functorName);
573 }
574
575 template<class FunctorType>
576 void printFunctor7(const FunctorType &functor, std::ostream &out, const std::string &functorName = "")
577 {
578 using Scalar = typename std::remove_const<typename std::remove_reference<decltype(functor(0,0,0,0,0,0,0))>::type>::type;
579 printFunctor<FunctorType, Scalar, 7>(functor, out, functorName);
580 }
581
582 template<class View>
583 void printView(const View &view, std::ostream &out, const std::string &viewName = "")
584 {
585 using Scalar = typename View::value_type;
586 using HostView = typename View::HostMirror;
587 using HostViewIteratorScalar = Intrepid2::ViewIterator<HostView, Scalar>;
588
589 auto hostView = getHostCopy(view);
590
591 HostViewIteratorScalar vi(hostView);
592
593 bool moreEntries = (vi.nextIncrementRank() != -1);
594 while (moreEntries)
595 {
596 Scalar value = vi.get();
597
598 auto location = vi.getLocation();
599 out << viewName << "(";
600 for (unsigned i=0; i<getFunctorRank(view); i++)
601 {
602 out << location[i];
603 if (i<getFunctorRank(view)-1)
604 {
605 out << ",";
606 }
607 }
608 out << ") " << value << std::endl;
609
610 moreEntries = (vi.increment() != -1);
611 }
612 }
613
614 template <class FunctorType1, class FunctorType2, int rank, typename Scalar=typename FunctorType1::value_type, class ExecutionSpace = typename FunctorType1::execution_space>
616 testFloatingEquality(const FunctorType1 &functor1, const FunctorType2 &functor2, double relTol, double absTol, Teuchos::FancyOStream &out, bool &success,
617 std::string functor1Name = "Functor 1", std::string functor2Name = "Functor 2")
618 {
619 INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "testFloatingEquality() called on FunctorType1 or FunctorType2 that does not support the specified rank.\n");
620 }
621
623 template <class FunctorType1, class FunctorType2, int rank, typename Scalar=typename FunctorType1::value_type, class ExecutionSpace = typename FunctorType1::execution_space>
625 testFloatingEquality(const FunctorType1 &functor1, const FunctorType2 &functor2, double relTol, double absTol, Teuchos::FancyOStream &out, bool &success,
626 std::string functor1Name = "Functor 1", std::string functor2Name = "Functor 2")
627 {
628 static_assert( supports_rank<FunctorType1,rank>::value, "Functor1 must support the specified rank through operator().");
629 static_assert( supports_rank<FunctorType2,rank>::value, "Functor2 must support the specified rank through operator().");
630
631 using Functor1IteratorScalar = FunctorIterator<FunctorType1, Scalar, rank>;
632 using Functor2IteratorScalar = FunctorIterator<FunctorType2, Scalar, rank>;
633
634 // check that rank/size match
635 TEUCHOS_TEST_FOR_EXCEPTION(getFunctorRank(functor1) != rank, std::invalid_argument, "functor1 and functor2 must agree in rank"); // Kokkos::View does not provide rank() method; getFunctorRank() provides common interface
636 TEUCHOS_TEST_FOR_EXCEPTION(getFunctorRank(functor2) != rank, std::invalid_argument, "functor1 and functor2 must agree in rank"); // Kokkos::View does not provide rank() method; getFunctorRank() provides common interface
637
638 int entryCount = 1;
639 for (unsigned i=0; i<getFunctorRank(functor1); i++)
640 {
641 TEUCHOS_TEST_FOR_EXCEPTION(functor1.extent_int(i) != functor2.extent_int(i), std::invalid_argument,
642 "functor1 and functor2 must agree in size in each dimension; functor1 has extent "
643 + std::to_string(functor1.extent_int(i)) + " in dimension " + std::to_string(i)
644 + "; functor2 has extent " + std::to_string(functor2.extent_int(i)) );
645 entryCount *= functor1.extent_int(i);
646 }
647 if (entryCount == 0) return; // nothing to test
648
649 ViewType<bool,ExecutionSpace> valuesMatch = getView<bool,ExecutionSpace>("valuesMatch", entryCount);
650
651 Kokkos::RangePolicy < ExecutionSpace > policy(0,entryCount);
652 Kokkos::parallel_for( policy,
653 KOKKOS_LAMBDA (const int &enumerationIndex )
654 {
655 Functor1IteratorScalar vi1(functor1);
656 Functor2IteratorScalar vi2(functor2);
657
658 vi1.setEnumerationIndex(enumerationIndex);
659 vi2.setEnumerationIndex(enumerationIndex);
660
661 const Scalar & value1 = vi1.get();
662 const Scalar & value2 = vi2.get();
663
664 bool errMeetsTol = errMeetsAbsAndRelTol(value1, value2, relTol, absTol);
665 valuesMatch(enumerationIndex) = errMeetsTol;
666 }
667 );
668
669 int failureCount = 0;
670 Kokkos::RangePolicy<ExecutionSpace > reducePolicy(0, entryCount);
671 Kokkos::parallel_reduce( reducePolicy,
672 KOKKOS_LAMBDA( const int &enumerationIndex, int &reducedValue )
673 {
674 reducedValue += valuesMatch(enumerationIndex) ? 0 : 1;
675 }, failureCount);
676
677 const bool allValuesMatch = (failureCount == 0);
678 success = success && allValuesMatch;
679
680 if (!allValuesMatch)
681 {
682 // copy functors to host views
683 auto functor1CopyHost = copyFunctorToHostView<FunctorType1,Scalar,rank>(functor1);
684 auto functor2CopyHost = copyFunctorToHostView<FunctorType2,Scalar,rank>(functor2);
685
686 auto valuesMatchHost = getHostCopy(valuesMatch);
687
688 FunctorIterator<decltype(functor1CopyHost),Scalar,rank> vi1(functor1CopyHost);
689 FunctorIterator<decltype(functor2CopyHost),Scalar,rank> vi2(functor2CopyHost);
690 Intrepid2::ViewIterator<decltype(valuesMatchHost), bool> viMatch(valuesMatchHost);
691
692 bool moreEntries = true;
693 while (moreEntries)
694 {
695 bool errMeetsTol = viMatch.get();
696
697 if (!errMeetsTol)
698 {
699 const Scalar value1 = vi1.get();
700 const Scalar value2 = vi2.get();
701
702 success = false;
703 auto location = vi1.getLocation();
704 out << "At location (";
705 for (unsigned i=0; i<getFunctorRank(functor1); i++)
706 {
707 out << location[i];
708 if (i<getFunctorRank(functor1)-1)
709 {
710 out << ",";
711 }
712 }
713 out << ") " << functor1Name << " value != " << functor2Name << " value (";
714 out << value1 << " != " << value2 << "); diff is " << std::abs(value1-value2) << std::endl;
715 }
716
717 moreEntries = (vi1.increment() != -1);
718 moreEntries = moreEntries && (vi2.increment() != -1);
719 moreEntries = moreEntries && (viMatch.increment() != -1);
720 }
721 }
722 }
723
724 template <class ViewType, class FunctorType>
725 void testFloatingEquality1(const ViewType &view, const FunctorType &functor, double relTol, double absTol, Teuchos::FancyOStream &out, bool &success,
726 std::string view1Name = "View", std::string view2Name = "Functor")
727 {
728 testFloatingEquality<ViewType, FunctorType, 1>(view, functor, relTol, absTol, out, success, view1Name, view2Name);
729 }
730
731 template <class ViewType, class FunctorType>
732 void testFloatingEquality2(const ViewType &view, const FunctorType &functor, double relTol, double absTol, Teuchos::FancyOStream &out, bool &success,
733 std::string view1Name = "View", std::string view2Name = "Functor")
734 {
735 testFloatingEquality<ViewType, FunctorType, 2>(view, functor, relTol, absTol, out, success, view1Name, view2Name);
736 }
737
738 template <class ViewType, class FunctorType>
739 void testFloatingEquality3(const ViewType &view, const FunctorType &functor, double relTol, double absTol, Teuchos::FancyOStream &out, bool &success,
740 std::string view1Name = "View", std::string view2Name = "Functor")
741 {
742 testFloatingEquality<ViewType, FunctorType, 3>(view, functor, relTol, absTol, out, success, view1Name, view2Name);
743 }
744
745 template <class ViewType, class FunctorType>
746 void testFloatingEquality4(const ViewType &view, const FunctorType &functor, double relTol, double absTol, Teuchos::FancyOStream &out, bool &success,
747 std::string view1Name = "View", std::string view2Name = "Functor")
748 {
749 testFloatingEquality<ViewType, FunctorType, 4>(view, functor, relTol, absTol, out, success, view1Name, view2Name);
750 }
751
752 template <class ViewType, class FunctorType>
753 void testFloatingEquality5(const ViewType &view, const FunctorType &functor, double relTol, double absTol, Teuchos::FancyOStream &out, bool &success,
754 std::string view1Name = "View", std::string view2Name = "Functor")
755 {
756 testFloatingEquality<ViewType, FunctorType, 5>(view, functor, relTol, absTol, out, success, view1Name, view2Name);
757 }
758
759 template <class ViewType, class FunctorType>
760 void testFloatingEquality6(const ViewType &view, const FunctorType &functor, double relTol, double absTol, Teuchos::FancyOStream &out, bool &success,
761 std::string view1Name = "View", std::string view2Name = "Functor")
762 {
763 testFloatingEquality<ViewType, FunctorType, 6>(view, functor, relTol, absTol, out, success, view1Name, view2Name);
764 }
765
766 template <class ViewType, class FunctorType>
767 void testFloatingEquality7(const ViewType &view, const FunctorType &functor, double relTol, double absTol, Teuchos::FancyOStream &out, bool &success,
768 std::string view1Name = "View", std::string view2Name = "Functor")
769 {
770 testFloatingEquality<ViewType, FunctorType, 7>(view, functor, relTol, absTol, out, success, view1Name, view2Name);
771 }
772
773 template <class ViewType1, class ViewType2>
774 void testViewFloatingEquality(const ViewType1 &view1, const ViewType2 &view2, double relTol, double absTol, Teuchos::FancyOStream &out, bool &success,
775 std::string view1Name = "View 1", std::string view2Name = "View 2")
776 {
777 // check that rank/size match
778 TEUCHOS_TEST_FOR_EXCEPTION(view1.rank() != view2.rank(), std::invalid_argument, "views must agree in rank");
779 for (unsigned i=0; i<view1.rank(); i++)
780 {
781 TEUCHOS_TEST_FOR_EXCEPTION(view1.extent_int(i) != view2.extent_int(i), std::invalid_argument, "views must agree in size in each dimension");
782 }
783
784 if (view1.size() == 0) return; // nothing to test
785
786 const int rank = view1.rank();
787 switch (rank)
788 {
789 case 0: testFloatingEquality<ViewType1, ViewType2, 0>(view1, view2, relTol, absTol, out, success, view1Name, view2Name); break;
790 case 1: testFloatingEquality<ViewType1, ViewType2, 1>(view1, view2, relTol, absTol, out, success, view1Name, view2Name); break;
791 case 2: testFloatingEquality<ViewType1, ViewType2, 2>(view1, view2, relTol, absTol, out, success, view1Name, view2Name); break;
792 case 3: testFloatingEquality<ViewType1, ViewType2, 3>(view1, view2, relTol, absTol, out, success, view1Name, view2Name); break;
793 case 4: testFloatingEquality<ViewType1, ViewType2, 4>(view1, view2, relTol, absTol, out, success, view1Name, view2Name); break;
794 case 5: testFloatingEquality<ViewType1, ViewType2, 5>(view1, view2, relTol, absTol, out, success, view1Name, view2Name); break;
795 case 6: testFloatingEquality<ViewType1, ViewType2, 6>(view1, view2, relTol, absTol, out, success, view1Name, view2Name); break;
796 case 7: testFloatingEquality<ViewType1, ViewType2, 7>(view1, view2, relTol, absTol, out, success, view1Name, view2Name); break;
797 default: INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument, "Unsupported rank");
798 }
799 }
800
801} // namespace Intrepid2
802
803#ifdef HAVE_INTREPID2_SACADO
804using Sacado_Fad_DFadType = Sacado::Fad::DFad<double>;
805#define INTREPID2_POINTSCALAR_TEST_INSTANT(GROUP_NAME, TEST_NAME) \
806\
807TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( GROUP_NAME, TEST_NAME, double ) \
808\
809TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( GROUP_NAME, TEST_NAME, Sacado_Fad_DFadType ) \
810
811#define INTREPID2_OUTPUTSCALAR_POINTSCALAR_TEST_INSTANT(GROUP_NAME, TEST_NAME) \
812\
813TEUCHOS_UNIT_TEST_TEMPLATE_2_INSTANT( GROUP_NAME, TEST_NAME, double, double ) \
814\
815TEUCHOS_UNIT_TEST_TEMPLATE_2_INSTANT( GROUP_NAME, TEST_NAME, Sacado_Fad_DFadType, Sacado_Fad_DFadType ) \
816
817#else
818#define INTREPID2_POINTSCALAR_TEST_INSTANT(GROUP_NAME, TEST_NAME) \
819\
820TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( GROUP_NAME, TEST_NAME, double ) \
821
822#define INTREPID2_OUTPUTSCALAR_POINTSCALAR_TEST_INSTANT(GROUP_NAME, TEST_NAME) \
823\
824TEUCHOS_UNIT_TEST_TEMPLATE_2_INSTANT( GROUP_NAME, TEST_NAME, double, double ) \
825
826#endif
827
828#endif /* Intrepid2_TestUtils_h */
Header file for the abstract base class Intrepid2::Basis.
Teuchos::RCP< Basis< DeviceType, OutputType, PointType > > BasisPtr
Basis Pointer.
KOKKOS_INLINE_FUNCTION ordinal_type getDkCardinality(const EOperator operatorType, const ordinal_type spaceDim)
Returns cardinality of Dk, i.e., the number of all derivatives of order k.
Stateless class representing a family of basis functions, templated on H(vol) and H(grad) on the line...
Defines the Intrepid2::FunctorIterator class, as well as the Intrepid2::functor_returns_ref SFINAE he...
Stateless classes that act as factories for two families of hierarchical bases. HierarchicalBasisFami...
Header file for Intrepid2::PointTools class to provide utilities for barycentric coordinates,...
Header file to include all Sacado headers that are required if using Intrepid2 with Sacado types.
ViewType< Scalar, DefaultTestDeviceType >::HostMirror copyFunctorToHostView(const Functor &deviceFunctor)
Copy the values for the specified functor.
ViewType< PointValueType, DeviceType > getInputPointsView(shards::CellTopology &cellTopo, int numPoints_1D)
Returns a DynRankView containing regularly-spaced points on the specified cell topology.
typename Kokkos::DefaultExecutionSpace::device_type DefaultTestDeviceType
Default Kokkos::Device to use for tests; depends on platform.
constexpr int MAX_FAD_DERIVATIVES_FOR_TESTS
Maximum number of derivatives to track for Fad types in tests.
const Teuchos::ScalarTraits< Scalar >::magnitudeType smallNumber()
Use Teuchos small number determination on host; pass this to Intrepid2::relErr() on device.
KOKKOS_INLINE_FUNCTION bool relErrMeetsTol(const Scalar1 &s1, const Scalar2 &s2, const typename Teuchos::ScalarTraits< typename std::common_type< Scalar1, Scalar2 >::type >::magnitudeType &smallNumber, const double &tol)
Adapted from Teuchos::relErr(); for use in code that may be executed on device.
Header function for Intrepid2::Util class and other utility functions.
enable_if_t< has_rank_method< Functor >::value, unsigned > KOKKOS_INLINE_FUNCTION getFunctorRank(const Functor &functor)
A family of basis functions, constructed from H(vol) and H(grad) bases on the line.
essentially, a read-only variant of ViewIterator, for a general functor (extent_int() and rank() supp...
Basis defining integrated Legendre basis on the line, a polynomial subspace of H(grad) on the line.
Basis defining Legendre basis on the line, a polynomial subspace of L^2 (a.k.a. H(vol)) on the line.
static ordinal_type getLatticeSize(const shards::CellTopology cellType, const ordinal_type order, const ordinal_type offset=0)
Computes the number of points in a lattice of a given order on a simplex (currently disabled for othe...
static void getLattice(Kokkos::DynRankView< pointValueType, pointProperties... > points, const shards::CellTopology cellType, const ordinal_type order, const ordinal_type offset=0, const EPointType pointType=POINTTYPE_EQUISPACED)
Computes a lattice of points of a given order on a reference simplex, quadrilateral or hexahedron (cu...
A helper class that allows iteration over some part of a Kokkos View, while allowing the calling code...
SFINAE helper to detect whether a type supports a rank-integral-argument operator().
Helper to get Scalar[*+] where the number of *'s matches the given rank.