Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MPAssembly/fenl_functors.hpp
Go to the documentation of this file.
1 /*
2 //@HEADER
3 // ************************************************************************
4 //
5 // Kokkos: Manycore Performance-Portable Multidimensional Arrays
6 // Copyright (2012) Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact H. Carter Edwards (hcedwar@sandia.gov)
39 //
40 // ************************************************************************
41 //@HEADER
42 */
43 
44 #ifndef KOKKOS_EXAMPLE_FENLFUNCTORS_HPP
45 #define KOKKOS_EXAMPLE_FENLFUNCTORS_HPP
46 
47 #include <stdio.h>
48 
49 #include <iostream>
50 #include <fstream>
51 #include <iomanip>
52 #include <cstdlib>
53 #include <cmath>
54 #include <limits>
55 
56 #include <Kokkos_Pair.hpp>
57 #include <Kokkos_UnorderedMap.hpp>
58 
59 #include <impl/Kokkos_Timer.hpp>
60 
61 #include <BoxElemFixture.hpp>
62 #include <HexElement.hpp>
63 
65 
66 //----------------------------------------------------------------------------
67 //----------------------------------------------------------------------------
68 
69 namespace Kokkos {
70 namespace Example {
71 namespace FENL {
72 
73 struct DeviceConfig {
74  struct Dim3 {
75  size_t x, y, z;
76  Dim3(const size_t x_, const size_t y_ = 1, const size_t z_ = 1) :
77  x(x_), y(y_), z(z_) {}
78  };
79 
80  Dim3 block_dim;
81  size_t num_blocks;
82  size_t num_threads_per_block;
83 
84  DeviceConfig(const size_t num_blocks_ = 0,
85  const size_t threads_per_block_x_ = 0,
86  const size_t threads_per_block_y_ = 0,
87  const size_t threads_per_block_z_ = 1) :
88  block_dim(threads_per_block_x_,threads_per_block_y_,threads_per_block_z_),
89  num_blocks(num_blocks_),
91  {}
92 };
93 
94 template< typename ValueType , class Space >
95 struct CrsMatrix {
96 #ifdef KOKKOS_ENABLE_DEPRECATED_CODE // Don't remove this until Kokkos has removed the deprecated code path probably around September 2018
97  typedef Kokkos::StaticCrsGraph< unsigned , Space , void , unsigned > StaticCrsGraphType ;
98 #else
99  typedef Kokkos::StaticCrsGraph< unsigned , Space , void , void , unsigned > StaticCrsGraphType ;
100 #endif
101  typedef View< ValueType * , Space > values_type ;
102 
105 
106  CrsMatrix() : graph(), values() {}
107 
108  CrsMatrix( const StaticCrsGraphType & arg_graph )
109  : graph( arg_graph )
110  , values( "crs_matrix_values" , arg_graph.entries.extent(0) )
111  {}
112 };
113 
114 // Traits class for creating strided local views for embedded ensemble-based,
115 // specialized for ensemble UQ scalar type
116 template <typename ViewType, typename Enabled = void>
117 struct LocalViewTraits {
118  typedef ViewType view_type;
119  // typedef Kokkos::View<typename view_type::data_type,
120  // typename view_type::array_layout,
121  // typename view_type::execution_space,
122  // Kokkos::MemoryUnmanaged> local_view_type;
123  typedef const view_type& local_view_type;
125  static const bool use_team = false;
126  KOKKOS_INLINE_FUNCTION
128  const unsigned local_rank)
129  { return v; }
130 };
131 
132 #if defined( KOKKOS_ENABLE_CUDA )
133 
134 template <typename ViewType>
135 struct LocalViewTraits<
136  ViewType,
137  typename std::enable_if< std::is_same<typename ViewType::execution_space,
138  Kokkos::Cuda>::value &&
139  Kokkos::is_view_mp_vector<ViewType>::value
140  >::type > {
141  typedef ViewType view_type;
144  static const bool use_team = true;
145 
146  KOKKOS_INLINE_FUNCTION
148  const unsigned local_rank)
149  {
150  return Kokkos::partition<1>(v, local_rank);
151  }
152 };
153 
154 #endif /* #if defined( KOKKOS_ENABLE_CUDA ) */
155 
156 // Compute DeviceConfig struct's based on scalar type
157 template <typename ScalarType>
158 struct CreateDeviceConfigs {
159  static void eval( Kokkos::Example::FENL::DeviceConfig& dev_config_elem,
160  Kokkos::Example::FENL::DeviceConfig& dev_config_bc ) {
161  dev_config_elem = Kokkos::Example::FENL::DeviceConfig( 0 , 1 , 1 );
162  dev_config_bc = Kokkos::Example::FENL::DeviceConfig( 0 , 1 , 1 );
163  }
164 };
165 
166 // Compute DeviceConfig struct's based on scalar type
167 template <typename StorageType>
168 struct CreateDeviceConfigs< Sacado::MP::Vector<StorageType> > {
170  static void eval( Kokkos::Example::FENL::DeviceConfig& dev_config_elem,
171  Kokkos::Example::FENL::DeviceConfig& dev_config_bc ) {
172  static const unsigned VectorSize = StorageType::static_size;
173 #if defined( KOKKOS_ENABLE_CUDA )
174  enum { is_cuda = Kokkos::Impl::is_same< execution_space, Kokkos::Cuda >::value };
175 #else
176  enum { is_cuda = false };
177 #endif /* #if defined( KOKKOS_ENABLE_CUDA ) */
178  if ( is_cuda ) {
179  dev_config_elem = Kokkos::Example::FENL::DeviceConfig( 0 , VectorSize , 64/VectorSize );
180  dev_config_bc = Kokkos::Example::FENL::DeviceConfig( 0 , VectorSize , 256/VectorSize );
181  }
182  else {
183  dev_config_elem = Kokkos::Example::FENL::DeviceConfig( 0 , 1 , 1 );
184  dev_config_bc = Kokkos::Example::FENL::DeviceConfig( 0 , 1 , 1 );
185  }
186  }
187 };
188 
189 template< class ElemNodeIdView , class CrsGraphType , unsigned ElemNode >
190 class NodeNodeGraph {
191 public:
192 
194  typedef pair<unsigned,unsigned> key_type ;
195 
196  typedef Kokkos::UnorderedMap< key_type, void , execution_space > SetType ;
197  typedef typename CrsGraphType::row_map_type::non_const_type RowMapType ;
198  typedef Kokkos::View< unsigned , execution_space > UnsignedValue ;
199 
200  // Static dimensions of 0 generate compiler warnings or errors.
201  typedef Kokkos::View< unsigned*[ElemNode][ElemNode] , execution_space >
203 
204 private:
205 
211 
212  const unsigned node_count ;
213  const ElemNodeIdView elem_node_id ;
218  PhaseType phase ;
219 
220 public:
221 
222  CrsGraphType graph ;
224 
225  struct Times
226  {
227  double ratio;
228  double fill_node_set;
229  double scan_node_count;
230  double fill_graph_entries;
231  double sort_graph_entries;
232  double fill_element_graph;
233  };
234 
235  NodeNodeGraph( const ElemNodeIdView & arg_elem_node_id ,
236  const unsigned arg_node_count,
237  Times & results
238  )
239  : node_count(arg_node_count)
240  , elem_node_id( arg_elem_node_id )
241  , row_total( "row_total" )
242  , row_count(Kokkos::ViewAllocateWithoutInitializing("row_count") , node_count ) // will deep_copy to 0 inside loop
243  , row_map( "graph_row_map" , node_count + 1 )
244  , node_node_set()
245  , phase( FILL_NODE_SET )
246  , graph()
247  , elem_graph()
248  {
249  //--------------------------------
250  // Guess at capacity required for the map:
251 
252  Kokkos::Impl::Timer wall_clock ;
253 
254  wall_clock.reset();
255  phase = FILL_NODE_SET ;
256 
257  // upper bound on the capacity
258  size_t set_capacity = (((28ull * node_count) / 2ull)*4ull)/3ull;
259 
260 
261  // Increase capacity until the (node,node) map is successfully filled.
262  {
263  // Zero the row count to restart the fill
265 
266  node_node_set = SetType( set_capacity );
267 
268  // May be larger that requested:
269  set_capacity = node_node_set.capacity();
270 
271  Kokkos::parallel_for( elem_node_id.extent(0) , *this );
272  }
273 
274  execution_space::fence();
275  results.ratio = (double)node_node_set.size() / (double)node_node_set.capacity();
276  results.fill_node_set = wall_clock.seconds();
277  //--------------------------------
278 
279  wall_clock.reset();
281 
282  // Exclusive scan of row_count into row_map
283  // including the final total in the 'node_count + 1' position.
284  // Zero the 'row_count' values.
285  Kokkos::parallel_scan( node_count , *this );
286 
287  // Zero the row count for the fill:
289 
290  unsigned graph_entry_count = 0 ;
291 
292  Kokkos::deep_copy( graph_entry_count , row_total );
293 
294  // Assign graph's row_map and allocate graph's entries
295  graph.row_map = row_map ;
296  graph.entries = typename CrsGraphType::entries_type( "graph_entries" , graph_entry_count );
297 
298  //--------------------------------
299  // Fill graph's entries from the (node,node) set.
300 
301  execution_space::fence();
302  results.scan_node_count = wall_clock.seconds();
303 
304  wall_clock.reset();
306  Kokkos::parallel_for( node_node_set.capacity() , *this );
307 
308  execution_space::fence();
309  results.fill_graph_entries = wall_clock.seconds();
310 
311  //--------------------------------
312  // Done with the temporary sets and arrays
313  wall_clock.reset();
315 
317  row_count = RowMapType();
318  row_map = RowMapType();
319  node_node_set.clear();
320 
321  //--------------------------------
322 
323  Kokkos::parallel_for( node_count , *this );
324 
325  execution_space::fence();
326  results.sort_graph_entries = wall_clock.seconds();
327 
328  //--------------------------------
329  // Element-to-graph mapping:
330  wall_clock.reset();
332  elem_graph = ElemGraphType("elem_graph", elem_node_id.extent(0) );
333  Kokkos::parallel_for( elem_node_id.extent(0) , *this );
334 
335  execution_space::fence();
336  results.fill_element_graph = wall_clock.seconds();
337  }
338 
339  //------------------------------------
340  // parallel_for: create map and count row length
341 
342  KOKKOS_INLINE_FUNCTION
343  void fill_set( const unsigned ielem ) const
344  {
345  // Loop over element's (row_local_node,col_local_node) pairs:
346  for ( unsigned row_local_node = 0 ; row_local_node < elem_node_id.extent(1) ; ++row_local_node ) {
347 
348  const unsigned row_node = elem_node_id( ielem , row_local_node );
349 
350  for ( unsigned col_local_node = row_local_node ; col_local_node < elem_node_id.extent(1) ; ++col_local_node ) {
351 
352  const unsigned col_node = elem_node_id( ielem , col_local_node );
353 
354  // If either node is locally owned then insert the pair into the unordered map:
355 
356  if ( row_node < row_count.extent(0) || col_node < row_count.extent(0) ) {
357 
358  const key_type key = (row_node < col_node) ? make_pair( row_node, col_node ) : make_pair( col_node, row_node ) ;
359 
360  const typename SetType::insert_result result = node_node_set.insert( key );
361 
362  if ( result.success() ) {
363  if ( row_node < row_count.extent(0) ) { atomic_fetch_add( & row_count( row_node ) , (typename RowMapType::value_type)1 ); }
364  if ( col_node < row_count.extent(0) && col_node != row_node ) { atomic_fetch_add( & row_count( col_node ) , (typename RowMapType::value_type)1 ); }
365  }
366  }
367  }
368  }
369  }
370 
371  KOKKOS_INLINE_FUNCTION
372  void fill_graph_entries( const unsigned iset ) const
373  {
374  if ( node_node_set.valid_at(iset) ) {
375  const key_type key = node_node_set.key_at(iset) ;
376  const unsigned row_node = key.first ;
377  const unsigned col_node = key.second ;
378 
379  if ( row_node < row_count.extent(0) ) {
380  const unsigned offset = graph.row_map( row_node ) + atomic_fetch_add( & row_count( row_node ) , (typename RowMapType::value_type)1 );
381  graph.entries( offset ) = col_node ;
382  }
383 
384  if ( col_node < row_count.extent(0) && col_node != row_node ) {
385  const unsigned offset = graph.row_map( col_node ) + atomic_fetch_add( & row_count( col_node ) , (typename RowMapType::value_type)1 );
386  graph.entries( offset ) = row_node ;
387  }
388  }
389  }
390 
391  KOKKOS_INLINE_FUNCTION
392  void sort_graph_entries( const unsigned irow ) const
393  {
394  typedef typename CrsGraphType::size_type size_type;
395  const size_type row_beg = graph.row_map( irow );
396  const size_type row_end = graph.row_map( irow + 1 );
397  for ( size_type i = row_beg + 1 ; i < row_end ; ++i ) {
398  const typename CrsGraphType::data_type col = graph.entries(i);
399  size_type j = i ;
400  for ( ; row_beg < j && col < graph.entries(j-1) ; --j ) {
401  graph.entries(j) = graph.entries(j-1);
402  }
403  graph.entries(j) = col ;
404  }
405  }
406 
407  KOKKOS_INLINE_FUNCTION
408  void fill_elem_graph_map( const unsigned ielem ) const
409  {
410  typedef typename CrsGraphType::data_type entry_type;
411  for ( unsigned row_local_node = 0 ; row_local_node < elem_node_id.extent(1) ; ++row_local_node ) {
412 
413  const unsigned row_node = elem_node_id( ielem , row_local_node );
414 
415  for ( unsigned col_local_node = 0 ; col_local_node < elem_node_id.extent(1) ; ++col_local_node ) {
416 
417  const unsigned col_node = elem_node_id( ielem , col_local_node );
418 
419  entry_type entry = 0 ;
420 
421  if ( row_node + 1 < graph.row_map.extent(0) ) {
422 
423  const entry_type entry_end = static_cast<entry_type> (graph.row_map( row_node + 1 ));
424 
425  entry = graph.row_map( row_node );
426 
427  for ( ; entry < entry_end && graph.entries(entry) != static_cast<entry_type> (col_node) ; ++entry );
428 
429  if ( entry == entry_end ) entry = ~0u ;
430  }
431 
432  elem_graph( ielem , row_local_node , col_local_node ) = entry ;
433  }
434  }
435  }
436 
437  KOKKOS_INLINE_FUNCTION
438  void operator()( const unsigned iwork ) const
439  {
440  if ( phase == FILL_NODE_SET ) {
441  fill_set( iwork );
442  }
443  else if ( phase == FILL_GRAPH_ENTRIES ) {
444  fill_graph_entries( iwork );
445  }
446  else if ( phase == SORT_GRAPH_ENTRIES ) {
447  sort_graph_entries( iwork );
448  }
449  else if ( phase == FILL_ELEMENT_GRAPH ) {
450  fill_elem_graph_map( iwork );
451  }
452  }
453 
454  //------------------------------------
455  // parallel_scan: row offsets
456 
457  typedef unsigned value_type ;
458 
459  KOKKOS_INLINE_FUNCTION
460  void operator()( const unsigned irow , unsigned & update , const bool final ) const
461  {
462  // exclusive scan
463  if ( final ) { row_map( irow ) = update ; }
464 
465  update += row_count( irow );
466 
467  if ( final ) {
468  if ( irow + 1 == row_count.extent(0) ) {
469  row_map( irow + 1 ) = update ;
470  row_total() = update ;
471  }
472  }
473  }
474 
475  KOKKOS_INLINE_FUNCTION
476  void init( unsigned & update ) const { update = 0 ; }
477 
478  KOKKOS_INLINE_FUNCTION
479  void join( volatile unsigned & update , const volatile unsigned & input ) const { update += input ; }
480 
481  //------------------------------------
482 };
483 
484 } /* namespace FENL */
485 } /* namespace Example */
486 } /* namespace Kokkos */
487 
488 //----------------------------------------------------------------------------
489 //----------------------------------------------------------------------------
490 
491 namespace Kokkos {
492 namespace Example {
493 namespace FENL {
494 
495 struct ElementComputationConstantCoefficient {
496  enum { is_constant = false };
497 
498  const double coeff_k ;
499 
500  KOKKOS_INLINE_FUNCTION
501  double operator()( double pt[], unsigned ensemble_rank) const
502  { return coeff_k * std::sin(pt[0]) * std::sin(pt[1]) * std::sin(pt[2]); }
503 
505  : coeff_k( val ) {}
506 
508  : coeff_k( rhs.coeff_k ) {}
509 };
510 
511 // Exponential KL from Stokhos
512 template < typename Scalar, typename MeshScalar, typename Device >
513 class ExponentialKLCoefficient {
514 public:
515 
516  // Turn into a meta-function class usable with Sacado::mpl::apply
517  template <typename T1, typename T2 = MeshScalar, typename T3 = Device>
518  struct apply {
520  };
521 
522  enum { is_constant = false };
523  typedef Kokkos::View<Scalar*, Kokkos::LayoutLeft, Device> RandomVariableView;
524  typedef typename RandomVariableView::size_type size_type;
525 
530 
531  rf_type m_rf; // Exponential random field
532  const MeshScalar m_mean; // Mean of random field
533  const MeshScalar m_variance; // Variance of random field
534  const MeshScalar m_corr_len; // Correlation length of random field
535  const size_type m_num_rv; // Number of random variables
536  RandomVariableView m_rv; // KL random variables
537 
538 public:
539 
541  const MeshScalar mean ,
542  const MeshScalar variance ,
543  const MeshScalar correlation_length ,
544  const size_type num_rv ) :
545  m_mean( mean ),
546  m_variance( variance ),
547  m_corr_len( correlation_length ),
548  m_num_rv( num_rv ),
549  m_rv( "KL Random Variables", m_num_rv )
550  {
551  Teuchos::ParameterList solverParams;
552  solverParams.set("Number of KL Terms", int(num_rv));
553  solverParams.set("Mean", mean);
554  solverParams.set("Standard Deviation", std::sqrt(variance));
555  int ndim = 3;
556  Teuchos::Array<double> domain_upper(ndim, 1.0), domain_lower(ndim, 0.0),
557  correlation_lengths(ndim, correlation_length);
558  solverParams.set("Domain Upper Bounds", domain_upper);
559  solverParams.set("Domain Lower Bounds", domain_lower);
560  solverParams.set("Correlation Lengths", correlation_lengths);
561 
562  m_rf = rf_type(solverParams);
563  }
564 
566  m_rf( rhs.m_rf ) ,
567  m_mean( rhs.m_mean ) ,
568  m_variance( rhs.m_variance ) ,
569  m_corr_len( rhs.m_corr_len ) ,
570  m_num_rv( rhs.m_num_rv ) ,
571  m_rv( rhs.m_rv ) {}
572 
573  KOKKOS_INLINE_FUNCTION
574  void setRandomVariables( const RandomVariableView& rv) { m_rv = rv; }
575 
576  KOKKOS_INLINE_FUNCTION
578 
579  KOKKOS_INLINE_FUNCTION
580  local_scalar_type operator() ( const MeshScalar point[],
581  const size_type ensemble_rank ) const
582  {
583  local_rv_view_type local_rv =
585 
586  local_scalar_type val = m_rf.evaluate(point, local_rv);
587 
588  return val;
589  }
590 };
591 
592 template< class FiniteElementMeshType , class SparseMatrixType
593  , class CoeffFunctionType = ElementComputationConstantCoefficient
594  >
595 class ElementComputation ;
596 
597 
598 template< class ExecutionSpace , BoxElemPart::ElemOrder Order , class CoordinateMap ,
599  typename ScalarType , class CoeffFunctionType >
601  < Kokkos::Example::BoxElemFixture< ExecutionSpace , Order , CoordinateMap >
602  , Kokkos::Example::FENL::CrsMatrix< ScalarType , ExecutionSpace >
603  , CoeffFunctionType >
604 {
605 public:
606 
609 
610  //------------------------------------
611 
612  typedef ExecutionSpace execution_space ;
613  typedef ScalarType scalar_type ;
614 
618  typedef Kokkos::View< scalar_type* , Kokkos::LayoutLeft, execution_space > vector_type ;
619 
620  //------------------------------------
621 
627  static const bool use_team = local_vector_view_traits::use_team;
628 
629  static const unsigned SpatialDim = element_data_type::spatial_dimension ;
630  static const unsigned TensorDim = SpatialDim * SpatialDim ;
631  static const unsigned ElemNodeCount = element_data_type::element_node_count ;
632  static const unsigned FunctionCount = element_data_type::function_count ;
633  static const unsigned IntegrationCount = element_data_type::integration_count ;
634 
635  //------------------------------------
636 
639  typedef Kokkos::View< scalar_type*[FunctionCount][FunctionCount] , execution_space > elem_matrices_type ;
640  typedef Kokkos::View< scalar_type*[FunctionCount] , execution_space > elem_vectors_type ;
641 
646 
648 
649  //------------------------------------
650 
651 
652  //------------------------------------
653  // Computational data:
654 
664  const CoeffFunctionType coeff_function ;
666 
668  : elem_data()
669  , elem_node_ids( rhs.elem_node_ids )
670  , node_coords( rhs.node_coords )
671  , elem_graph( rhs.elem_graph )
672  , elem_jacobians( rhs.elem_jacobians )
673  , elem_residuals( rhs.elem_residuals )
674  , solution( rhs.solution )
675  , residual( rhs.residual )
676  , jacobian( rhs.jacobian )
677  , coeff_function( rhs.coeff_function )
678  , dev_config( rhs.dev_config )
679  {}
680 
681  // If the element->sparse_matrix graph is provided then perform atomic updates
682  // Otherwise fill per-element contributions for subequent gather-add into a residual and jacobian.
683  ElementComputation( const mesh_type & arg_mesh ,
684  const CoeffFunctionType & arg_coeff_function ,
685  const vector_type & arg_solution ,
686  const elem_graph_type & arg_elem_graph ,
687  const sparse_matrix_type & arg_jacobian ,
688  const vector_type & arg_residual ,
689  const Kokkos::Example::FENL::DeviceConfig arg_dev_config )
690  : elem_data()
691  , elem_node_ids( arg_mesh.elem_node() )
692  , node_coords( arg_mesh.node_coord() )
693  , elem_graph( arg_elem_graph )
694  , elem_jacobians()
695  , elem_residuals()
696  , solution( arg_solution )
697  , residual( arg_residual )
698  , jacobian( arg_jacobian )
699  , coeff_function( arg_coeff_function )
700  , dev_config( arg_dev_config )
701  {}
702 
703  //------------------------------------
704 
705  void apply() const
706  {
707  const size_t nelem = elem_node_ids.extent(0);
708  if ( use_team ) {
709  const size_t team_size = dev_config.block_dim.x * dev_config.block_dim.y;
710  const size_t league_size =
711  (nelem + dev_config.block_dim.y-1) / dev_config.block_dim.y;
712  Kokkos::TeamPolicy< execution_space > config( league_size, team_size );
713  parallel_for( config , *this );
714  }
715  else {
716  parallel_for( nelem , *this );
717  }
718  }
719 
720  //------------------------------------
721 
722  static const unsigned FLOPS_transform_gradients =
723  /* Jacobian */ FunctionCount * TensorDim * 2 +
724  /* Inverse jacobian */ TensorDim * 6 + 6 +
725  /* Gradient transform */ FunctionCount * 15 ;
726 
727  KOKKOS_INLINE_FUNCTION
729  const double grad[][ FunctionCount ] , // Gradient of bases master element
730  const double x[] ,
731  const double y[] ,
732  const double z[] ,
733  double dpsidx[] ,
734  double dpsidy[] ,
735  double dpsidz[] ) const
736  {
737  enum { j11 = 0 , j12 = 1 , j13 = 2 ,
738  j21 = 3 , j22 = 4 , j23 = 5 ,
739  j31 = 6 , j32 = 7 , j33 = 8 };
740 
741  // Jacobian accumulation:
742 
743  double J[ TensorDim ] = { 0, 0, 0, 0, 0, 0, 0, 0, 0 };
744 
745  for( unsigned i = 0; i < FunctionCount ; ++i ) {
746  const double x1 = x[i] ;
747  const double x2 = y[i] ;
748  const double x3 = z[i] ;
749 
750  const double g1 = grad[0][i] ;
751  const double g2 = grad[1][i] ;
752  const double g3 = grad[2][i] ;
753 
754  J[j11] += g1 * x1 ;
755  J[j12] += g1 * x2 ;
756  J[j13] += g1 * x3 ;
757 
758  J[j21] += g2 * x1 ;
759  J[j22] += g2 * x2 ;
760  J[j23] += g2 * x3 ;
761 
762  J[j31] += g3 * x1 ;
763  J[j32] += g3 * x2 ;
764  J[j33] += g3 * x3 ;
765  }
766 
767  // Inverse jacobian:
768 
769  double invJ[ TensorDim ] = {
770  static_cast<double>( J[j22] * J[j33] - J[j23] * J[j32] ) ,
771  static_cast<double>( J[j13] * J[j32] - J[j12] * J[j33] ) ,
772  static_cast<double>( J[j12] * J[j23] - J[j13] * J[j22] ) ,
773 
774  static_cast<double>( J[j23] * J[j31] - J[j21] * J[j33] ) ,
775  static_cast<double>( J[j11] * J[j33] - J[j13] * J[j31] ) ,
776  static_cast<double>( J[j13] * J[j21] - J[j11] * J[j23] ) ,
777 
778  static_cast<double>( J[j21] * J[j32] - J[j22] * J[j31] ) ,
779  static_cast<double>( J[j12] * J[j31] - J[j11] * J[j32] ) ,
780  static_cast<double>( J[j11] * J[j22] - J[j12] * J[j21] ) };
781 
782  const double detJ = J[j11] * invJ[j11] +
783  J[j21] * invJ[j12] +
784  J[j31] * invJ[j13] ;
785 
786  const double detJinv = 1.0 / detJ ;
787 
788  for ( unsigned i = 0 ; i < TensorDim ; ++i ) { invJ[i] *= detJinv ; }
789 
790  // Transform gradients:
791 
792  for( unsigned i = 0; i < FunctionCount ; ++i ) {
793  const double g0 = grad[0][i];
794  const double g1 = grad[1][i];
795  const double g2 = grad[2][i];
796 
797  dpsidx[i] = g0 * invJ[j11] + g1 * invJ[j12] + g2 * invJ[j13];
798  dpsidy[i] = g0 * invJ[j21] + g1 * invJ[j22] + g2 * invJ[j23];
799  dpsidz[i] = g0 * invJ[j31] + g1 * invJ[j32] + g2 * invJ[j33];
800  }
801 
802  return detJ ;
803  }
804 
805  KOKKOS_INLINE_FUNCTION
807  const local_scalar_type dof_values[] ,
808  const double dpsidx[] ,
809  const double dpsidy[] ,
810  const double dpsidz[] ,
811  const double detJ ,
812  const local_scalar_type coeff_k ,
813  const double integ_weight ,
814  const double bases_vals[] ,
815  local_scalar_type elem_res[] ,
816  local_scalar_type elem_mat[][ FunctionCount ] ) const
817  {
818  local_scalar_type value_at_pt = 0 ;
819  local_scalar_type gradx_at_pt = 0 ;
820  local_scalar_type grady_at_pt = 0 ;
821  local_scalar_type gradz_at_pt = 0 ;
822 
823  for ( unsigned m = 0 ; m < FunctionCount ; m++ ) {
824  value_at_pt += dof_values[m] * bases_vals[m] ;
825  gradx_at_pt += dof_values[m] * dpsidx[m] ;
826  grady_at_pt += dof_values[m] * dpsidy[m] ;
827  gradz_at_pt += dof_values[m] * dpsidz[m] ;
828  }
829 
830  const local_scalar_type k_detJ_weight = coeff_k * detJ * integ_weight ;
831  const local_scalar_type res_val = value_at_pt * value_at_pt * detJ * integ_weight ;
832  const local_scalar_type mat_val = 2.0 * value_at_pt * detJ * integ_weight ;
833 
834  // $$ R_i = \int_{\Omega} \nabla \phi_i \cdot (k \nabla T) + \phi_i T^2 d \Omega $$
835  // $$ J_{i,j} = \frac{\partial R_i}{\partial T_j} = \int_{\Omega} k \nabla \phi_i \cdot \nabla \phi_j + 2 \phi_i \phi_j T d \Omega $$
836 
837  for ( unsigned m = 0; m < FunctionCount; ++m) {
838  local_scalar_type * const mat = elem_mat[m] ;
839  const double bases_val_m = bases_vals[m];
840  const double dpsidx_m = dpsidx[m] ;
841  const double dpsidy_m = dpsidy[m] ;
842  const double dpsidz_m = dpsidz[m] ;
843 
844  elem_res[m] += k_detJ_weight * ( dpsidx_m * gradx_at_pt +
845  dpsidy_m * grady_at_pt +
846  dpsidz_m * gradz_at_pt ) +
847  res_val * bases_val_m ;
848 
849  for( unsigned n = 0; n < FunctionCount; n++) {
850 
851  mat[n] += k_detJ_weight * ( dpsidx_m * dpsidx[n] +
852  dpsidy_m * dpsidy[n] +
853  dpsidz_m * dpsidz[n] ) +
854  mat_val * bases_val_m * bases_vals[n];
855  }
856  }
857  }
858 
859  KOKKOS_INLINE_FUNCTION
860  void operator()( const typename TeamPolicy< execution_space >::member_type & dev ) const
861  {
862 
863  const unsigned num_ensemble_threads = dev_config.block_dim.x ;
864  const unsigned num_element_threads = dev_config.block_dim.y ;
865  const unsigned element_rank = dev.team_rank() / num_ensemble_threads ;
866  const unsigned ensemble_rank = dev.team_rank() % num_ensemble_threads ;
867 
868  const unsigned ielem =
869  dev.league_rank() * num_element_threads + element_rank;
870 
871  if (ielem >= elem_node_ids.extent(0))
872  return;
873 
874  (*this)( ielem, ensemble_rank );
875 
876  }
877 
878  KOKKOS_INLINE_FUNCTION
879  void operator()( const unsigned ielem ,
880  const unsigned ensemble_rank = 0 ) const
881  {
882  local_vector_type local_solution =
883  local_vector_view_traits::create_local_view(solution,
884  ensemble_rank);
885  local_vector_type local_residual =
886  local_vector_view_traits::create_local_view(residual,
887  ensemble_rank);
888  local_matrix_type local_jacobian_values =
889  local_matrix_view_traits::create_local_view(jacobian.values,
890  ensemble_rank);
891 
892  // Gather nodal coordinates and solution vector:
893 
894  double x[ FunctionCount ] ;
895  double y[ FunctionCount ] ;
896  double z[ FunctionCount ] ;
897  local_scalar_type val[ FunctionCount ] ;
898  unsigned node_index[ ElemNodeCount ];
899 
900  for ( unsigned i = 0 ; i < ElemNodeCount ; ++i ) {
901  const unsigned ni = elem_node_ids( ielem , i );
902 
903  node_index[i] = ni ;
904 
905  x[i] = node_coords( ni , 0 );
906  y[i] = node_coords( ni , 1 );
907  z[i] = node_coords( ni , 2 );
908 
909  val[i] = local_solution( ni );
910  }
911 
912 
913  local_scalar_type elem_vec[ FunctionCount ] ;
914  local_scalar_type elem_mat[ FunctionCount ][ FunctionCount ] ;
915 
916  for( unsigned i = 0; i < FunctionCount ; i++ ) {
917  elem_vec[i] = 0 ;
918  for( unsigned j = 0; j < FunctionCount ; j++){
919  elem_mat[i][j] = 0 ;
920  }
921  }
922 
923 
924  for ( unsigned i = 0 ; i < IntegrationCount ; ++i ) {
925  double dpsidx[ FunctionCount ] ;
926  double dpsidy[ FunctionCount ] ;
927  double dpsidz[ FunctionCount ] ;
928 
929  local_scalar_type coeff_k = 0 ;
930 
931  {
932  double pt[] = {0.0, 0.0, 0.0};
933 
934  // If function is not constant
935  // then compute physical coordinates of integration point
936  if ( ! coeff_function.is_constant ) {
937  for ( unsigned j = 0 ; j < FunctionCount ; ++j ) {
938  pt[0] += x[j] * elem_data.values[i][j] ;
939  pt[1] += y[j] * elem_data.values[i][j] ;
940  pt[2] += z[j] * elem_data.values[i][j] ;
941  }
942  }
943 
944  // Need to fix this for local_scalar_type!!!!!!
945  coeff_k = coeff_function(pt, ensemble_rank);
946  }
947 
948  const double detJ =
949  transform_gradients( elem_data.gradients[i] , x , y , z ,
950  dpsidx , dpsidy , dpsidz );
951 
952  contributeResidualJacobian( val , dpsidx , dpsidy , dpsidz ,
953  detJ , coeff_k ,
954  elem_data.weights[i] ,
955  elem_data.values[i] ,
956  elem_vec , elem_mat );
957  }
958 
959  for( unsigned i = 0 ; i < FunctionCount ; i++ ) {
960  const unsigned row = node_index[i] ;
961  if ( row < residual.extent(0) ) {
962  atomic_add( & local_residual( row ) , elem_vec[i] );
963 
964  for( unsigned j = 0 ; j < FunctionCount ; j++ ) {
965  const unsigned entry = elem_graph( ielem , i , j );
966  if ( entry != ~0u ) {
967  atomic_add( & local_jacobian_values( entry ) , elem_mat[i][j] );
968  }
969  }
970  }
971  }
972  }
973 }; /* ElementComputation */
974 
975 //----------------------------------------------------------------------------
976 
977 template< class FixtureType , class SparseMatrixType >
978 class DirichletComputation ;
979 
980 template< class ExecutionSpace , BoxElemPart::ElemOrder Order , class CoordinateMap ,
981  typename ScalarType >
982 class DirichletComputation<
983  Kokkos::Example::BoxElemFixture< ExecutionSpace , Order , CoordinateMap > ,
984  Kokkos::Example::FENL::CrsMatrix< ScalarType , ExecutionSpace > >
985 {
986 public:
987 
991 
992  typedef ExecutionSpace execution_space ;
993  typedef ScalarType scalar_type ;
994 
998  typedef Kokkos::View< scalar_type* , execution_space > vector_type ;
999 
1000  //------------------------------------
1001 
1006  static const bool use_team = local_vector_view_traits::use_team;
1007 
1008  typedef double bc_scalar_type ;
1009 
1010  //------------------------------------
1011  // Computational data:
1012 
1013  const node_coord_type node_coords ;
1014  const vector_type solution ;
1015  const sparse_matrix_type jacobian ;
1016  const vector_type residual ;
1017  const bc_scalar_type bc_lower_value ;
1018  const bc_scalar_type bc_upper_value ;
1019  const scalar_coord_type bc_lower_limit ;
1020  const scalar_coord_type bc_upper_limit ;
1021  const unsigned bc_plane ;
1022  const unsigned node_count ;
1023  bool init ;
1024  const Kokkos::Example::FENL::DeviceConfig dev_config ;
1025 
1026 
1027  DirichletComputation( const mesh_type & arg_mesh ,
1028  const vector_type & arg_solution ,
1029  const sparse_matrix_type & arg_jacobian ,
1030  const vector_type & arg_residual ,
1031  const unsigned arg_bc_plane ,
1032  const bc_scalar_type arg_bc_lower_value ,
1033  const bc_scalar_type arg_bc_upper_value ,
1034  const Kokkos::Example::FENL::DeviceConfig arg_dev_config )
1035  : node_coords( arg_mesh.node_coord() )
1036  , solution( arg_solution )
1037  , jacobian( arg_jacobian )
1038  , residual( arg_residual )
1039  , bc_lower_value( arg_bc_lower_value )
1040  , bc_upper_value( arg_bc_upper_value )
1041  , bc_lower_limit( std::numeric_limits<scalar_coord_type>::epsilon() )
1042  , bc_upper_limit( scalar_coord_type(1) - std::numeric_limits<scalar_coord_type>::epsilon() )
1043  , bc_plane( arg_bc_plane )
1044  , node_count( arg_mesh.node_count_owned() )
1045  , init( false )
1046  , dev_config( arg_dev_config )
1047  {
1048  parallel_for( node_count , *this );
1049  init = true ;
1050  }
1051 
1052  void apply() const
1053  {
1054  if ( use_team ) {
1055  const size_t team_size = dev_config.block_dim.x * dev_config.block_dim.y;
1056  const size_t league_size =
1057  (node_count + dev_config.block_dim.y-1) / dev_config.block_dim.y;
1058  Kokkos::TeamPolicy< execution_space > config( league_size, team_size );
1059  parallel_for( config , *this );
1060  }
1061  else
1062  parallel_for( node_count , *this );
1063  }
1064 
1065  //------------------------------------
1066 
1067  KOKKOS_INLINE_FUNCTION
1068  void operator()( const typename TeamPolicy< execution_space >::member_type & dev ) const
1069  {
1070 
1071  const unsigned num_ensemble_threads = dev_config.block_dim.x ;
1072  const unsigned num_node_threads = dev_config.block_dim.y ;
1073  const unsigned node_rank = dev.team_rank() / num_ensemble_threads ;
1074  const unsigned ensemble_rank = dev.team_rank() % num_ensemble_threads ;
1075 
1076  const unsigned inode = dev.league_rank() * num_node_threads + node_rank;
1077 
1078  if (inode >= node_count)
1079  return;
1080 
1081  (*this)( inode, ensemble_rank );
1082 
1083  }
1084 
1085  KOKKOS_INLINE_FUNCTION
1086  void operator()( const unsigned inode ,
1087  const unsigned ensemble_rank = 0) const
1088  {
1089  local_vector_type local_residual =
1090  local_vector_view_traits::create_local_view(residual,
1091  ensemble_rank);
1092  local_matrix_type local_jacobian_values =
1093  local_matrix_view_traits::create_local_view(jacobian.values,
1094  ensemble_rank);
1095 
1096  // Apply dirichlet boundary condition on the Solution and Residual vectors.
1097  // To maintain the symmetry of the original global stiffness matrix,
1098  // zero out the columns that correspond to boundary conditions, and
1099  // update the residual vector accordingly
1100 
1101  const unsigned iBeg = jacobian.graph.row_map[inode];
1102  const unsigned iEnd = jacobian.graph.row_map[inode+1];
1103 
1104  const scalar_coord_type c = node_coords(inode,bc_plane);
1105  const bool bc_lower = c <= bc_lower_limit ;
1106  const bool bc_upper = bc_upper_limit <= c ;
1107 
1108  if ( ! init ) {
1109  solution(inode) = bc_lower ? bc_lower_value : (
1110  bc_upper ? bc_upper_value : 0 );
1111  }
1112  else {
1113  if ( bc_lower || bc_upper ) {
1114 
1115  local_residual(inode) = 0 ;
1116 
1117  // zero each value on the row, and leave a one
1118  // on the diagonal
1119 
1120  for( unsigned i = iBeg ; i < iEnd ; ++i ) {
1121  local_jacobian_values(i) = int(inode) == int(jacobian.graph.entries(i)) ? 1 : 0 ;
1122  }
1123  }
1124  else {
1125 
1126  // Find any columns that are boundary conditions.
1127  // Clear them and adjust the residual vector
1128 
1129  for( unsigned i = iBeg ; i < iEnd ; ++i ) {
1130  const unsigned cnode = jacobian.graph.entries(i) ;
1131  const scalar_coord_type cc = node_coords(cnode,bc_plane);
1132 
1133  if ( ( cc <= bc_lower_limit ) || ( bc_upper_limit <= cc ) ) {
1134  local_jacobian_values(i) = 0 ;
1135  }
1136  }
1137  }
1138  }
1139  }
1140 };
1141 
1142 template< typename FixtureType , typename VectorType >
1144 {
1145 public:
1146 
1147  typedef FixtureType fixture_type ;
1148  typedef VectorType vector_type ;
1151 
1153  static const unsigned SpatialDim = element_data_type::spatial_dimension ;
1154  static const unsigned TensorDim = SpatialDim * SpatialDim ;
1155  static const unsigned ElemNodeCount = element_data_type::element_node_count ;
1156  static const unsigned IntegrationCount = element_data_type::integration_count ;
1157 
1158  //------------------------------------
1159  // Computational data:
1160 
1164 
1166  : elem_data()
1167  , fixture( rhs.fixture )
1168  , solution( rhs.solution )
1169  {}
1170 
1171  ResponseComputation( const fixture_type& arg_fixture ,
1172  const vector_type & arg_solution )
1173  : elem_data()
1174  , fixture( arg_fixture )
1175  , solution( arg_solution )
1176  {}
1177 
1178  //------------------------------------
1179 
1181  {
1182  value_type response = 0;
1183  //Kokkos::parallel_reduce( fixture.elem_count() , *this , response );
1184  Kokkos::parallel_reduce( solution.extent(0) , *this , response );
1185  return response;
1186  }
1187 
1188  //------------------------------------
1189 
1190  KOKKOS_INLINE_FUNCTION
1192  const double grad[][ ElemNodeCount ] , // Gradient of bases master element
1193  const double x[] ,
1194  const double y[] ,
1195  const double z[] ) const
1196  {
1197  enum { j11 = 0 , j12 = 1 , j13 = 2 ,
1198  j21 = 3 , j22 = 4 , j23 = 5 ,
1199  j31 = 6 , j32 = 7 , j33 = 8 };
1200 
1201  // Jacobian accumulation:
1202 
1203  double J[ TensorDim ] = { 0, 0, 0, 0, 0, 0, 0, 0, 0 };
1204 
1205  for( unsigned i = 0; i < ElemNodeCount ; ++i ) {
1206  const double x1 = x[i] ;
1207  const double x2 = y[i] ;
1208  const double x3 = z[i] ;
1209 
1210  const double g1 = grad[0][i] ;
1211  const double g2 = grad[1][i] ;
1212  const double g3 = grad[2][i] ;
1213 
1214  J[j11] += g1 * x1 ;
1215  J[j12] += g1 * x2 ;
1216  J[j13] += g1 * x3 ;
1217 
1218  J[j21] += g2 * x1 ;
1219  J[j22] += g2 * x2 ;
1220  J[j23] += g2 * x3 ;
1221 
1222  J[j31] += g3 * x1 ;
1223  J[j32] += g3 * x2 ;
1224  J[j33] += g3 * x3 ;
1225  }
1226 
1227  // Inverse jacobian:
1228 
1229  double invJ[ TensorDim ] = {
1230  static_cast<double>( J[j22] * J[j33] - J[j23] * J[j32] ) ,
1231  static_cast<double>( J[j13] * J[j32] - J[j12] * J[j33] ) ,
1232  static_cast<double>( J[j12] * J[j23] - J[j13] * J[j22] ) ,
1233 
1234  static_cast<double>( J[j23] * J[j31] - J[j21] * J[j33] ) ,
1235  static_cast<double>( J[j11] * J[j33] - J[j13] * J[j31] ) ,
1236  static_cast<double>( J[j13] * J[j21] - J[j11] * J[j23] ) ,
1237 
1238  static_cast<double>( J[j21] * J[j32] - J[j22] * J[j31] ) ,
1239  static_cast<double>( J[j12] * J[j31] - J[j11] * J[j32] ) ,
1240  static_cast<double>( J[j11] * J[j22] - J[j12] * J[j21] ) };
1241 
1242  const double detJ = J[j11] * invJ[j11] +
1243  J[j21] * invJ[j12] +
1244  J[j31] * invJ[j13] ;
1245 
1246  return detJ ;
1247  }
1248 
1249  KOKKOS_INLINE_FUNCTION
1251  const value_type dof_values[] ,
1252  const double detJ ,
1253  const double integ_weight ,
1254  const double bases_vals[] ) const
1255  {
1256  // $$ g_i = \int_{\Omega} T^2 d \Omega $$
1257 
1258  value_type value_at_pt = 0 ;
1259  for ( unsigned m = 0 ; m < ElemNodeCount ; m++ ) {
1260  value_at_pt += dof_values[m] * bases_vals[m] ;
1261  }
1262 
1263  value_type elem_response =
1264  value_at_pt * value_at_pt * detJ * integ_weight ;
1265 
1266  return elem_response;
1267  }
1268 
1269  /*
1270  KOKKOS_INLINE_FUNCTION
1271  void operator()( const unsigned ielem , value_type& response ) const
1272  {
1273  // Gather nodal coordinates and solution vector:
1274 
1275  double x[ ElemNodeCount ] ;
1276  double y[ ElemNodeCount ] ;
1277  double z[ ElemNodeCount ] ;
1278  value_type val[ ElemNodeCount ] ;
1279 
1280  for ( unsigned i = 0 ; i < ElemNodeCount ; ++i ) {
1281  const unsigned ni = fixture.elem_node( ielem , i );
1282 
1283  x[i] = fixture.node_coord( ni , 0 );
1284  y[i] = fixture.node_coord( ni , 1 );
1285  z[i] = fixture.node_coord( ni , 2 );
1286 
1287  val[i] = solution( ni );
1288  }
1289 
1290  for ( unsigned i = 0 ; i < IntegrationCount ; ++i ) {
1291 
1292  const double detJ = compute_detJ( elem_data.gradients[i] , x , y , z );
1293 
1294  response += contributeResponse( val , detJ , elem_data.weights[i] ,
1295  elem_data.values[i] );
1296  }
1297  }
1298  */
1299 
1300  KOKKOS_INLINE_FUNCTION
1301  void operator()( const unsigned i , value_type& response ) const
1302  {
1303  value_type u = solution(i);
1304  response += (u * u) / fixture.node_count_global();
1305  }
1306 
1307  KOKKOS_INLINE_FUNCTION
1308  void init( value_type & response ) const
1309  { response = 0 ; }
1310 
1311  KOKKOS_INLINE_FUNCTION
1312  void join( volatile value_type & response ,
1313  volatile const value_type & input ) const
1314  { response += input ; }
1315 
1316 }; /* ResponseComputation */
1317 
1318 } /* namespace FENL */
1319 } /* namespace Example */
1320 } /* namespace Kokkos */
1321 
1322 //----------------------------------------------------------------------------
1323 
1324 /* A Cuda-specific specialization for the element computation functor. */
1325 #if defined( __CUDACC__ )
1326 // #include <NonlinearElement_Cuda.hpp>
1327 #endif
1328 
1329 //----------------------------------------------------------------------------
1330 
1331 #endif /* #ifndef KOKKOS_EXAMPLE_FENLFUNCTORS_HPP */
static KOKKOS_INLINE_FUNCTION local_view_type create_local_view(const view_type &v, const unsigned local_rank)
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION double compute_detJ(const double grad[][ElemNodeCount], const double x[], const double y[], const double z[]) const
KOKKOS_INLINE_FUNCTION void operator()(const unsigned i, value_type &response) const
ExponentialKLCoefficient(const ExponentialKLCoefficient &rhs)
KOKKOS_INLINE_FUNCTION void join(volatile value_type &response, volatile const value_type &input) const
KOKKOS_INLINE_FUNCTION void setRandomVariables(const RandomVariableView &rv)
LocalViewTraits< RandomVariableView > local_rv_view_traits
KOKKOS_INLINE_FUNCTION Teuchos::PromotionTraits< typename rv_type::value_type, value_type >::promote evaluate(const point_type &point, const rv_type &random_variables) const
Evaluate random field at a point.
KOKKOS_INLINE_FUNCTION void fill_graph_entries(const unsigned iset) const
static void eval(Kokkos::Example::FENL::DeviceConfig &dev_config_elem, Kokkos::Example::FENL::DeviceConfig &dev_config_bc)
KOKKOS_INLINE_FUNCTION value_type contributeResponse(const value_type dof_values[], const double detJ, const double integ_weight, const double bases_vals[]) const
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Kokkos::View< const double *[SpaceDim], Device > node_coord_type
Kokkos::DefaultExecutionSpace execution_space
Kokkos::View< const unsigned *[ElemNode], Device > elem_node_type
KOKKOS_INLINE_FUNCTION RandomVariableView getRandomVariables() const
static void eval(Kokkos::Example::FENL::DeviceConfig &dev_config_elem, Kokkos::Example::FENL::DeviceConfig &dev_config_bc)
KOKKOS_INLINE_FUNCTION void atomic_add(volatile Sacado::UQ::PCE< Storage > *const dest, const Sacado::UQ::PCE< Storage > &src)
Kokkos::View< Scalar *, Kokkos::LayoutLeft, Device > RandomVariableView
KOKKOS_INLINE_FUNCTION local_scalar_type operator()(const MeshScalar point[], const size_type ensemble_rank) const
local_rv_view_traits::local_view_type local_rv_view_type
ElemNodeIdView::execution_space execution_space
ResponseComputation(const fixture_type &arg_fixture, const vector_type &arg_solution)
KOKKOS_INLINE_FUNCTION void operator()(const unsigned irow, unsigned &update, const bool final) const
ElementComputationConstantCoefficient(const ElementComputationConstantCoefficient &rhs)
KOKKOS_INLINE_FUNCTION void init(unsigned &update) const
KOKKOS_INLINE_FUNCTION void fill_elem_graph_map(const unsigned ielem) const
Dim3(const size_t x_, const size_t y_=1, const size_t z_=1)
CrsGraphType::row_map_type::non_const_type RowMapType
KOKKOS_INLINE_FUNCTION void join(volatile unsigned &update, const volatile unsigned &input) const
KOKKOS_INLINE_FUNCTION void init(value_type &response) const
CrsMatrix(const StaticCrsGraphType &arg_graph)
Kokkos::StaticCrsGraph< unsigned, Space, void, void, unsigned > StaticCrsGraphType
void deep_copy(const Stokhos::CrsMatrix< ValueType, DstDevice, Layout > &dst, const Stokhos::CrsMatrix< ValueType, SrcDevice, Layout > &src)
KOKKOS_INLINE_FUNCTION void operator()(const unsigned iwork) const
KOKKOS_INLINE_FUNCTION void contributeResidualJacobian(const local_scalar_type dof_values[], const double dpsidx[], const double dpsidy[], const double dpsidz[], const double detJ, const local_scalar_type coeff_k, const double integ_weight, const double bases_vals[], local_scalar_type elem_res[], local_scalar_type elem_mat[][FunctionCount]) const
KOKKOS_INLINE_FUNCTION double operator()(double pt[], unsigned ensemble_rank) const
NodeNodeGraph(const ElemNodeIdView &arg_elem_node_id, const unsigned arg_node_count, Times &results)
expr val()
ElementComputation(const mesh_type &arg_mesh, const CoeffFunctionType &arg_coeff_function, const vector_type &arg_solution, const elem_graph_type &arg_elem_graph, const sparse_matrix_type &arg_jacobian, const vector_type &arg_residual, const Kokkos::Example::FENL::DeviceConfig arg_dev_config)
KOKKOS_INLINE_FUNCTION double transform_gradients(const double grad[][FunctionCount], const double x[], const double y[], const double z[], double dpsidx[], double dpsidy[], double dpsidz[]) const
Kokkos::View< unsigned, execution_space > UnsignedValue
DeviceConfig(const size_t num_blocks_=0, const size_t threads_per_block_x_=0, const size_t threads_per_block_y_=0, const size_t threads_per_block_z_=1)
Kokkos::Example::HexElement_Data< fixture_type::ElemNode > element_data_type
local_rv_view_traits::local_value_type local_scalar_type
DirichletComputation(const mesh_type &arg_mesh, const vector_type &arg_solution, const sparse_matrix_type &arg_jacobian, const vector_type &arg_residual, const unsigned arg_bc_plane, const bc_scalar_type arg_bc_lower_value, const bc_scalar_type arg_bc_upper_value, const Kokkos::Example::FENL::DeviceConfig arg_dev_config)
KOKKOS_INLINE_FUNCTION PCE< Storage > sin(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION void sort_graph_entries(const unsigned irow) const
KOKKOS_INLINE_FUNCTION void fill_set(const unsigned ielem) const
Kokkos::View< unsigned *[ElemNode][ElemNode], execution_space > ElemGraphType
int n
Generate a distributed unstructured finite element mesh from a partitioned NX*NY*NZ box of elements...
Stokhos::KL::ExponentialRandomField< MeshScalar, Device > rf_type
ExponentialKLCoefficient(const MeshScalar mean, const MeshScalar variance, const MeshScalar correlation_length, const size_type num_rv)
void update(const ValueType &alpha, VectorType &x, const ValueType &beta, const VectorType &y)
Kokkos::UnorderedMap< key_type, void, execution_space > SetType