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