44 #ifndef KOKKOS_EXAMPLE_FENLFUNCTORS_HPP 45 #define KOKKOS_EXAMPLE_FENLFUNCTORS_HPP 56 #include <Kokkos_Pair.hpp> 57 #include <Kokkos_UnorderedMap.hpp> 59 #include <impl/Kokkos_Timer.hpp> 80 Dim3(
const size_t x_,
const size_t y_ = 1,
const size_t z_ = 1) :
81 x(x_),
y(y_),
z(z_) {}
89 const size_t threads_per_block_x_ = 0,
90 const size_t threads_per_block_y_ = 0,
91 const size_t threads_per_block_z_ = 1) :
92 block_dim(threads_per_block_x_,threads_per_block_y_,threads_per_block_z_),
98 template<
typename ValueType ,
class Space >
110 ,
values(
"crs_matrix_values" , arg_graph.entries.dimension_0() )
116 template <
typename ViewType,
typename Enabled =
void>
126 KOKKOS_INLINE_FUNCTION
128 const unsigned local_rank)
132 #if defined( KOKKOS_HAVE_CUDA ) 134 template <
typename ViewType>
135 struct LocalViewTraits<
137 typename std::enable_if< std::is_same<typename ViewType::execution_space,
138 Kokkos::Cuda>::value &&
139 Kokkos::is_view_mp_vector<ViewType>::value
146 KOKKOS_INLINE_FUNCTION
148 const unsigned local_rank)
150 return Kokkos::partition<1>(v, local_rank);
157 template <
typename ScalarType>
167 template <
typename StorageType>
172 static const unsigned VectorSize = StorageType::static_size;
173 #if defined( KOKKOS_HAVE_CUDA ) 174 enum { is_cuda = Kokkos::Impl::is_same< execution_space, Kokkos::Cuda >::value };
176 enum { is_cuda =
false };
189 template<
class ElemNodeIdView ,
class CrsGraphType ,
unsigned ElemNode >
196 typedef Kokkos::UnorderedMap< key_type, void , execution_space >
SetType ;
197 typedef typename CrsGraphType::row_map_type::non_const_type
RowMapType ;
201 typedef Kokkos::View< unsigned*[ElemNode][ElemNode] , execution_space >
236 const unsigned arg_node_count,
252 Kokkos::Impl::Timer wall_clock ;
258 size_t set_capacity = (((28ull *
node_count) / 2ull)*4ull)/3ull;
271 Kokkos::parallel_for(
elem_node_id.dimension_0() , *this );
274 execution_space::fence();
290 unsigned graph_entry_count = 0 ;
296 graph.entries =
typename CrsGraphType::entries_type(
"graph_entries" , graph_entry_count );
301 execution_space::fence();
308 execution_space::fence();
325 execution_space::fence();
333 Kokkos::parallel_for(
elem_node_id.dimension_0() , *this );
335 execution_space::fence();
342 KOKKOS_INLINE_FUNCTION
346 for (
unsigned row_local_node = 0 ; row_local_node <
elem_node_id.dimension_1() ; ++row_local_node ) {
348 const unsigned row_node =
elem_node_id( ielem , row_local_node );
350 for (
unsigned col_local_node = row_local_node ; col_local_node <
elem_node_id.dimension_1() ; ++col_local_node ) {
352 const unsigned col_node =
elem_node_id( ielem , col_local_node );
358 const key_type key = (row_node < col_node) ? make_pair( row_node, col_node ) : make_pair( col_node, row_node ) ;
360 const typename SetType::insert_result result =
node_node_set.insert( key );
362 if ( result.success() ) {
363 if ( row_node <
row_count.dimension_0() ) { atomic_fetch_add( &
row_count( row_node ) , 1 ); }
364 if ( col_node <
row_count.dimension_0() && col_node != row_node ) { atomic_fetch_add( &
row_count( col_node ) , 1 ); }
371 KOKKOS_INLINE_FUNCTION
376 const unsigned row_node = key.first ;
377 const unsigned col_node = key.second ;
379 if ( row_node <
row_count.dimension_0() ) {
380 const unsigned offset =
graph.row_map( row_node ) + atomic_fetch_add( &
row_count( row_node ) , 1 );
381 graph.entries( offset ) = col_node ;
384 if ( col_node <
row_count.dimension_0() && col_node != row_node ) {
385 const unsigned offset =
graph.row_map( col_node ) + atomic_fetch_add( &
row_count( col_node ) , 1 );
386 graph.entries( offset ) = row_node ;
391 KOKKOS_INLINE_FUNCTION
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);
400 for ( ; row_beg <
j && col <
graph.entries(
j-1) ; --
j ) {
407 KOKKOS_INLINE_FUNCTION
410 typedef typename CrsGraphType::data_type entry_type;
411 for (
unsigned row_local_node = 0 ; row_local_node <
elem_node_id.dimension_1() ; ++row_local_node ) {
413 const unsigned row_node =
elem_node_id( ielem , row_local_node );
415 for (
unsigned col_local_node = 0 ; col_local_node <
elem_node_id.dimension_1() ; ++col_local_node ) {
417 const unsigned col_node =
elem_node_id( ielem , col_local_node );
419 entry_type entry = 0 ;
421 if ( row_node + 1 <
graph.row_map.dimension_0() ) {
423 const entry_type entry_end =
static_cast<entry_type
> (
graph.row_map( row_node + 1 ));
425 entry =
graph.row_map( row_node );
427 for ( ; entry < entry_end &&
graph.entries(entry) !=
static_cast<entry_type
> (col_node) ; ++entry );
429 if ( entry == entry_end ) entry = ~0u ;
432 elem_graph( ielem , row_local_node , col_local_node ) = entry ;
437 KOKKOS_INLINE_FUNCTION
459 KOKKOS_INLINE_FUNCTION
468 if ( irow + 1 ==
row_count.dimension_0() ) {
475 KOKKOS_INLINE_FUNCTION
478 KOKKOS_INLINE_FUNCTION
479 void join(
volatile unsigned &
update ,
const volatile unsigned & input )
const {
update += input ; }
500 KOKKOS_INLINE_FUNCTION
501 double operator()(
double pt[],
unsigned ensemble_rank)
const 512 template <
typename Scalar,
typename MeshScalar,
typename Device >
517 template <
typename T1,
typename T2 = MeshScalar,
typename T3 = Device>
524 typedef typename RandomVariableView::size_type
size_type;
541 const MeshScalar mean ,
542 const MeshScalar variance ,
543 const MeshScalar correlation_length ,
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));
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);
573 KOKKOS_INLINE_FUNCTION
576 KOKKOS_INLINE_FUNCTION
579 KOKKOS_INLINE_FUNCTION
593 class CoordinateMap ,
typename ScalarType >
608 typedef Kokkos::View< scalar_type* , Kokkos::LayoutLeft, execution_space >
vector_type ;
622 typedef Kokkos::View< scalar_type*[FunctionCount][FunctionCount] , execution_space >
elem_matrices_type ;
673 KOKKOS_INLINE_FUNCTION
681 double dpsidz[] )
const 683 enum { j11 = 0 , j12 = 1 , j13 = 2 ,
684 j21 = 3 , j22 = 4 , j23 = 5 ,
685 j31 = 6 , j32 = 7 , j33 = 8 };
689 double J[
TensorDim ] = { 0, 0, 0, 0, 0, 0, 0, 0, 0 };
692 const double x1 =
x[i] ;
693 const double x2 =
y[i] ;
694 const double x3 = z[i] ;
696 const double g1 = grad[0][i] ;
697 const double g2 = grad[1][i] ;
698 const double g3 = grad[2][i] ;
716 static_cast<double>( J[j22] * J[j33] - J[j23] * J[j32] ) ,
717 static_cast<double>( J[j13] * J[j32] - J[j12] * J[j33] ) ,
718 static_cast<double>( J[j12] * J[j23] - J[j13] * J[j22] ) ,
720 static_cast<double>( J[j23] * J[j31] - J[j21] * J[j33] ) ,
721 static_cast<double>( J[j11] * J[j33] - J[j13] * J[j31] ) ,
722 static_cast<double>( J[j13] * J[j21] - J[j11] * J[j23] ) ,
724 static_cast<double>( J[j21] * J[j32] - J[j22] * J[j31] ) ,
725 static_cast<double>( J[j12] * J[j31] - J[j11] * J[j32] ) ,
726 static_cast<double>( J[j11] * J[j22] - J[j12] * J[j21] ) };
728 const double detJ = J[j11] * invJ[j11] +
732 const double detJinv = 1.0 / detJ ;
734 for (
unsigned i = 0 ; i <
TensorDim ; ++i ) { invJ[i] *= detJinv ; }
739 const double g0 = grad[0][i];
740 const double g1 = grad[1][i];
741 const double g2 = grad[2][i];
743 dpsidx[i] = g0 * invJ[j11] + g1 * invJ[j12] + g2 * invJ[j13];
744 dpsidy[i] = g0 * invJ[j21] + g1 * invJ[j22] + g2 * invJ[j23];
745 dpsidz[i] = g0 * invJ[j31] + g1 * invJ[j32] + g2 * invJ[j33];
760 template<
class FiniteElementMeshType ,
761 class SparseMatrixType ,
763 class CoeffFunctionType = ElementComputationConstantCoefficient
768 class CoordinateMap ,
typename ScalarType ,
class CoeffFunctionType >
771 CrsMatrix< ScalarType , ExecutionSpace > ,
783 static const unsigned FunctionCount = base_type::FunctionCount;
784 static const unsigned IntegrationCount = base_type::IntegrationCount;
785 static const unsigned ElemNodeCount = base_type::ElemNodeCount;
792 coeff_function( rhs.coeff_function ) ,
793 dev_config( rhs.dev_config ) {}
797 const CoeffFunctionType & arg_coeff_function ,
803 base_type(arg_mesh, arg_solution, arg_elem_graph, arg_jacobian,
805 coeff_function( arg_coeff_function ) ,
806 dev_config( arg_dev_config ) {}
812 const size_t nelem = this->elem_node_ids.dimension_0();
813 parallel_for( nelem , *
this );
816 KOKKOS_INLINE_FUNCTION
819 unsigned node_index[],
820 double x[],
double y[],
double z[],
824 for (
unsigned i = 0 ; i < ElemNodeCount ; ++i ) {
825 const unsigned ni = this->elem_node_ids( ielem , i );
829 x[i] = this->node_coords( ni , 0 );
830 y[i] = this->node_coords( ni , 1 );
831 z[i] = this->node_coords( ni , 2 );
833 val[i] = this->solution( ni ) ;
836 for(
unsigned j = 0;
j < FunctionCount ;
j++){
842 KOKKOS_INLINE_FUNCTION
844 const unsigned node_index[],
848 for(
unsigned i = 0 ; i < FunctionCount ; i++ ) {
849 const unsigned row = node_index[i] ;
850 if ( row < this->residual.dimension_0() ) {
851 atomic_add( & this->residual( row ) , res[i] );
853 for(
unsigned j = 0 ;
j < FunctionCount ;
j++ ) {
854 const unsigned entry = this->elem_graph( ielem , i ,
j );
855 if ( entry != ~0u ) {
856 atomic_add( & this->jacobian.values( entry ) , mat[i][
j] );
863 KOKKOS_INLINE_FUNCTION
873 double coeff_src = 1.234;
874 double advection[] = { 1.1, 1.2, 1.3 };
875 double dpsidx[ FunctionCount ] ;
876 double dpsidy[ FunctionCount ] ;
877 double dpsidz[ FunctionCount ] ;
878 double pt[] = {0.0, 0.0, 0.0};
879 for (
unsigned i = 0 ; i < IntegrationCount ; ++i ) {
883 if ( ! coeff_function.is_constant ) {
884 for (
unsigned j = 0 ;
j < FunctionCount ; ++
j ) {
885 pt[0] +=
x[
j] * this->elem_data.values[i][
j] ;
886 pt[1] +=
y[
j] * this->elem_data.values[i][
j] ;
887 pt[2] += z[
j] * this->elem_data.values[i][
j] ;
890 coeff_k = coeff_function(pt, 0);
892 const double integ_weight = this->elem_data.weights[i];
893 const double* bases_vals = this->elem_data.values[i];
895 this->transform_gradients( this->elem_data.gradients[i] ,
897 dpsidx , dpsidy , dpsidz );
898 const double detJ_weight = detJ * integ_weight;
899 const scalar_type detJ_weight_coeff_k = detJ_weight * coeff_k;
905 for (
unsigned m = 0 ; m < FunctionCount ; m++ ) {
906 value_at_pt += dof_values[m] * bases_vals[m] ;
907 gradx_at_pt += dof_values[m] * dpsidx[m] ;
908 grady_at_pt += dof_values[m] * dpsidy[m] ;
909 gradz_at_pt += dof_values[m] * dpsidz[m] ;
913 coeff_src * value_at_pt * value_at_pt ;
915 2.0 * coeff_src * value_at_pt ;
922 advection_x*gradx_at_pt +
923 advection_y*grady_at_pt +
924 advection_z*gradz_at_pt ;
926 for (
unsigned m = 0; m < FunctionCount; ++m) {
928 const double bases_val_m = bases_vals[m] * detJ_weight ;
929 const double dpsidx_m = dpsidx[m] ;
930 const double dpsidy_m = dpsidy[m] ;
931 const double dpsidz_m = dpsidz[m] ;
934 detJ_weight_coeff_k * ( dpsidx_m * gradx_at_pt +
935 dpsidy_m * grady_at_pt +
936 dpsidz_m * gradz_at_pt ) +
937 bases_val_m * ( advection_term + source_term ) ;
939 for(
unsigned n = 0; n < FunctionCount; n++) {
940 const double dpsidx_n = dpsidx[n] ;
941 const double dpsidy_n = dpsidy[n] ;
942 const double dpsidz_n = dpsidz[n] ;
944 detJ_weight_coeff_k * ( dpsidx_m * dpsidx_n +
945 dpsidy_m * dpsidy_n +
946 dpsidz_m * dpsidz_n ) +
947 bases_val_m * ( advection_x * dpsidx_n +
948 advection_y * dpsidy_n +
949 advection_z * dpsidz_n +
950 source_deriv * bases_vals[n] ) ;
956 KOKKOS_INLINE_FUNCTION
959 double x[ FunctionCount ] ;
960 double y[ FunctionCount ] ;
961 double z[ FunctionCount ] ;
962 unsigned node_index[ ElemNodeCount ];
966 scalar_type elem_mat[ FunctionCount ][ FunctionCount ] ;
969 gatherSolution(ielem,
val, node_index,
x,
y, z, elem_res, elem_mat);
972 computeElementResidualJacobian(
val,
x,
y, z, elem_res , elem_mat );
975 scatterResidual( ielem, node_index, elem_res, elem_mat );
980 class CoordinateMap ,
typename ScalarType ,
class CoeffFunctionType >
983 CrsMatrix< ScalarType , ExecutionSpace > ,
995 static const unsigned FunctionCount = base_type::FunctionCount;
996 static const unsigned IntegrationCount = base_type::IntegrationCount;
997 static const unsigned ElemNodeCount = base_type::ElemNodeCount;
1007 coeff_function( rhs.coeff_function ) ,
1008 dev_config( rhs.dev_config ) {}
1012 const CoeffFunctionType & arg_coeff_function ,
1018 base_type(arg_mesh, arg_solution, arg_elem_graph,
1019 arg_jacobian, arg_residual),
1020 coeff_function( arg_coeff_function ) ,
1021 dev_config( arg_dev_config ) {}
1027 const size_t nelem = this->elem_node_ids.dimension_0();
1028 parallel_for( nelem , *
this );
1031 KOKKOS_INLINE_FUNCTION
1034 unsigned node_index[],
1035 double x[],
double y[],
double z[],
1038 for (
unsigned i = 0 ; i < ElemNodeCount ; ++i ) {
1039 const unsigned ni = this->elem_node_ids( ielem , i );
1041 node_index[i] = ni ;
1043 x[i] = this->node_coords( ni , 0 );
1044 y[i] = this->node_coords( ni , 1 );
1045 z[i] = this->node_coords( ni , 2 );
1047 val[i].val() = this->solution( ni );
1048 val[i].diff( i, FunctionCount );
1052 KOKKOS_INLINE_FUNCTION
1054 const unsigned node_index[],
1057 for(
unsigned i = 0 ; i < FunctionCount ; i++ ) {
1058 const unsigned row = node_index[i] ;
1059 if ( row < this->residual.dimension_0() ) {
1062 for(
unsigned j = 0 ;
j < FunctionCount ;
j++ ) {
1063 const unsigned entry = this->elem_graph( ielem , i ,
j );
1064 if ( entry != ~0u ) {
1065 atomic_add( & this->jacobian.values( entry ) ,
1066 res[i].fastAccessDx(
j) );
1073 template <
typename local_scalar_type>
1074 KOKKOS_INLINE_FUNCTION
1079 local_scalar_type elem_res[] )
const 1082 double coeff_src = 1.234;
1083 double advection[] = { 1.1, 1.2, 1.3 };
1084 double dpsidx[ FunctionCount ] ;
1085 double dpsidy[ FunctionCount ] ;
1086 double dpsidz[ FunctionCount ] ;
1087 double pt[] = {0.0, 0.0, 0.0};
1088 for (
unsigned i = 0 ; i < IntegrationCount ; ++i ) {
1092 if ( ! coeff_function.is_constant ) {
1093 for (
unsigned j = 0 ;
j < FunctionCount ; ++
j ) {
1094 pt[0] +=
x[
j] * this->elem_data.values[i][
j] ;
1095 pt[1] +=
y[
j] * this->elem_data.values[i][
j] ;
1096 pt[2] += z[
j] * this->elem_data.values[i][
j] ;
1099 coeff_k = coeff_function(pt, 0);
1101 const double integ_weight = this->elem_data.weights[i];
1102 const double* bases_vals = this->elem_data.values[i];
1104 this->transform_gradients( this->elem_data.gradients[i] ,
1106 dpsidx , dpsidy , dpsidz );
1107 const double detJ_weight = detJ * integ_weight;
1108 const scalar_type detJ_weight_coeff_k = detJ_weight * coeff_k;
1110 local_scalar_type value_at_pt = 0 ;
1111 local_scalar_type gradx_at_pt = 0 ;
1112 local_scalar_type grady_at_pt = 0 ;
1113 local_scalar_type gradz_at_pt = 0 ;
1114 for (
unsigned m = 0 ; m < FunctionCount ; m++ ) {
1115 value_at_pt += dof_values[m] * bases_vals[m] ;
1116 gradx_at_pt += dof_values[m] * dpsidx[m] ;
1117 grady_at_pt += dof_values[m] * dpsidy[m] ;
1118 gradz_at_pt += dof_values[m] * dpsidz[m] ;
1121 const local_scalar_type source_term =
1122 coeff_src * value_at_pt * value_at_pt ;
1124 const local_scalar_type advection_term =
1125 advection[0]*gradx_at_pt +
1126 advection[1]*grady_at_pt +
1127 advection[2]*gradz_at_pt;
1129 for (
unsigned m = 0; m < FunctionCount; ++m) {
1130 const double bases_val_m = bases_vals[m] * detJ_weight ;
1131 const double dpsidx_m = dpsidx[m] ;
1132 const double dpsidy_m = dpsidy[m] ;
1133 const double dpsidz_m = dpsidz[m] ;
1136 detJ_weight_coeff_k * ( dpsidx_m * gradx_at_pt +
1137 dpsidy_m * grady_at_pt +
1138 dpsidz_m * gradz_at_pt ) +
1139 bases_val_m * ( advection_term + source_term ) ;
1144 KOKKOS_INLINE_FUNCTION
1147 double x[ FunctionCount ] ;
1148 double y[ FunctionCount ] ;
1149 double z[ FunctionCount ] ;
1150 unsigned node_index[ ElemNodeCount ];
1156 gatherSolution( ielem,
val, node_index,
x,
y, z, elem_res );
1159 computeElementResidual(
val,
x,
y, z, elem_res );
1162 scatterResidual( ielem, node_index, elem_res );
1167 class CoordinateMap ,
typename ScalarType ,
class CoeffFunctionType >
1170 CrsMatrix< ScalarType , ExecutionSpace > ,
1172 public ElementComputation< Kokkos::Example::BoxElemFixture< ExecutionSpace , Order , CoordinateMap > ,
1173 CrsMatrix< ScalarType , ExecutionSpace > ,
1174 FadElement, CoeffFunctionType > {
1184 static const unsigned FunctionCount = base_type::FunctionCount;
1185 static const unsigned IntegrationCount = base_type::IntegrationCount;
1186 static const unsigned ElemNodeCount = base_type::ElemNodeCount;
1195 const CoeffFunctionType & arg_coeff_function ,
1201 base_type(arg_mesh, arg_coeff_function, arg_solution, arg_elem_graph,
1202 arg_jacobian, arg_residual, arg_dev_config) {}
1208 const size_t nelem = this->elem_node_ids.dimension_0();
1209 parallel_for( nelem , *
this );
1212 KOKKOS_INLINE_FUNCTION
1215 unsigned node_index[],
1216 double x[],
double y[],
double z[],
1219 for (
unsigned i = 0 ; i < ElemNodeCount ; ++i ) {
1220 const unsigned ni = this->elem_node_ids( ielem , i );
1222 node_index[i] = ni ;
1224 x[i] = this->node_coords( ni , 0 );
1225 y[i] = this->node_coords( ni , 1 );
1226 z[i] = this->node_coords( ni , 2 );
1228 val[i] = this->solution( ni );
1232 template <
typename local_scalar_type>
1233 KOKKOS_INLINE_FUNCTION
1238 local_scalar_type elem_res[] )
const 1241 double coeff_src = 1.234;
1242 double advection[] = { 1.1, 1.2, 1.3 };
1243 double dpsidx[ FunctionCount ] ;
1244 double dpsidy[ FunctionCount ] ;
1245 double dpsidz[ FunctionCount ] ;
1246 double pt[] = {0.0, 0.0, 0.0};
1247 for (
unsigned i = 0 ; i < IntegrationCount ; ++i ) {
1251 if ( ! this->coeff_function.is_constant ) {
1252 for (
unsigned j = 0 ;
j < FunctionCount ; ++
j ) {
1253 pt[0] +=
x[
j] * this->elem_data.values[i][
j] ;
1254 pt[1] +=
y[
j] * this->elem_data.values[i][
j] ;
1255 pt[2] += z[
j] * this->elem_data.values[i][
j] ;
1258 coeff_k = this->coeff_function(pt, 0);
1260 const double integ_weight = this->elem_data.weights[i];
1261 const double* bases_vals = this->elem_data.values[i];
1263 this->transform_gradients( this->elem_data.gradients[i] ,
1265 dpsidx , dpsidy , dpsidz );
1266 const double detJ_weight = detJ * integ_weight;
1267 const scalar_type detJ_weight_coeff_k = detJ_weight * coeff_k;
1269 local_scalar_type value_at_pt(FunctionCount,
scalar_type(0.0), Sacado::NoInitDerivArray) ;
1270 local_scalar_type gradx_at_pt(FunctionCount,
scalar_type(0.0), Sacado::NoInitDerivArray) ;
1271 local_scalar_type grady_at_pt(FunctionCount,
scalar_type(0.0), Sacado::NoInitDerivArray) ;
1272 local_scalar_type gradz_at_pt(FunctionCount,
scalar_type(0.0), Sacado::NoInitDerivArray) ;
1273 for (
unsigned m = 0 ; m < FunctionCount ; m++ ) {
1274 value_at_pt.val() += dof_values[m] * bases_vals[m] ;
1275 value_at_pt.fastAccessDx(m) = bases_vals[m] ;
1277 gradx_at_pt.val() += dof_values[m] * dpsidx[m] ;
1278 gradx_at_pt.fastAccessDx(m) = dpsidx[m] ;
1280 grady_at_pt.val() += dof_values[m] * dpsidy[m] ;
1281 grady_at_pt.fastAccessDx(m) = dpsidy[m] ;
1283 gradz_at_pt.val() += dof_values[m] * dpsidz[m] ;
1284 gradz_at_pt.fastAccessDx(m) = dpsidz[m] ;
1287 const local_scalar_type source_term =
1288 coeff_src * value_at_pt * value_at_pt ;
1290 const local_scalar_type advection_term =
1291 advection[0]*gradx_at_pt +
1292 advection[1]*grady_at_pt +
1293 advection[2]*gradz_at_pt;
1295 for (
unsigned m = 0; m < FunctionCount; ++m) {
1296 const double bases_val_m = bases_vals[m] * detJ_weight ;
1297 const double dpsidx_m = dpsidx[m] ;
1298 const double dpsidy_m = dpsidy[m] ;
1299 const double dpsidz_m = dpsidz[m] ;
1302 detJ_weight_coeff_k * ( dpsidx_m * gradx_at_pt +
1303 dpsidy_m * grady_at_pt +
1304 dpsidz_m * gradz_at_pt ) +
1305 bases_val_m * ( advection_term + source_term ) ;
1310 KOKKOS_INLINE_FUNCTION
1313 double x[ FunctionCount ] ;
1314 double y[ FunctionCount ] ;
1315 double z[ FunctionCount ] ;
1316 unsigned node_index[ ElemNodeCount ];
1322 gatherSolution( ielem,
val, node_index,
x,
y, z, elem_res );
1325 computeElementResidual(
val,
x,
y, z, elem_res );
1328 this->scatterResidual( ielem, node_index, elem_res );
1333 class CoordinateMap ,
typename ScalarType ,
class CoeffFunctionType >
1336 CrsMatrix< ScalarType , ExecutionSpace > ,
1338 public ElementComputation< Kokkos::Example::BoxElemFixture< ExecutionSpace , Order , CoordinateMap > ,
1339 CrsMatrix< ScalarType , ExecutionSpace > ,
1340 Analytic , CoeffFunctionType > {
1350 static const unsigned FunctionCount = base_type::FunctionCount;
1351 static const unsigned IntegrationCount = base_type::IntegrationCount;
1352 static const unsigned ElemNodeCount = base_type::ElemNodeCount;
1361 const CoeffFunctionType & arg_coeff_function ,
1367 base_type(arg_mesh, arg_coeff_function, arg_solution, arg_elem_graph,
1368 arg_jacobian, arg_residual, arg_dev_config) {}
1374 const size_t nelem = this->elem_node_ids.dimension_0();
1375 parallel_for( nelem , *
this );
1378 KOKKOS_INLINE_FUNCTION
1388 double coeff_src = 1.234;
1389 double advection[] = { 1.1, 1.2, 1.3 };
1390 double dpsidx[ FunctionCount ] ;
1391 double dpsidy[ FunctionCount ] ;
1392 double dpsidz[ FunctionCount ] ;
1393 double pt[] = {0.0, 0.0, 0.0};
1399 for (
unsigned i = 0 ; i < IntegrationCount ; ++i ) {
1403 if ( ! this->coeff_function.is_constant ) {
1404 for (
unsigned j = 0 ;
j < FunctionCount ; ++
j ) {
1405 pt[0] +=
x[
j] * this->elem_data.values[i][
j] ;
1406 pt[1] +=
y[
j] * this->elem_data.values[i][
j] ;
1407 pt[2] += z[
j] * this->elem_data.values[i][
j] ;
1410 coeff_k = this->coeff_function(pt, 0);
1412 const double integ_weight = this->elem_data.weights[i];
1413 const double* bases_vals = this->elem_data.values[i];
1415 this->transform_gradients( this->elem_data.gradients[i] ,
1417 dpsidx , dpsidy , dpsidz );
1418 const double detJ_weight = detJ * integ_weight;
1419 const scalar_type detJ_weight_coeff_k = detJ_weight * coeff_k;
1421 value_at_pt.val() = 0.0 ;
1422 gradx_at_pt.val() = 0.0 ;
1423 grady_at_pt.val() = 0.0 ;
1424 gradz_at_pt.val() = 0.0 ;
1425 for (
unsigned m = 0 ; m < FunctionCount ; m++ ) {
1426 value_at_pt.val() += dof_values[m] * bases_vals[m] ;
1427 gradx_at_pt.val() += dof_values[m] * dpsidx[m] ;
1428 grady_at_pt.val() += dof_values[m] * dpsidy[m] ;
1429 gradz_at_pt.val() += dof_values[m] * dpsidz[m] ;
1433 coeff_src * value_at_pt * value_at_pt ;
1436 advection[0]*gradx_at_pt +
1437 advection[1]*grady_at_pt +
1438 advection[2]*gradz_at_pt;
1440 for (
unsigned m = 0; m < FunctionCount; ++m) {
1441 const double bases_val_m = bases_vals[m] * detJ_weight ;
1443 detJ_weight_coeff_k * ( dpsidx[m] * gradx_at_pt +
1444 dpsidy[m] * grady_at_pt +
1445 dpsidz[m] * gradz_at_pt ) +
1446 bases_val_m * ( advection_term + source_term ) ;
1448 elem_res[m] += res.val();
1451 for(
unsigned n = 0; n < FunctionCount; n++) {
1452 mat[n] += res.fastAccessDx(0) * bases_vals[n] +
1453 res.fastAccessDx(1) * dpsidx[n] +
1454 res.fastAccessDx(2) * dpsidy[n] +
1455 res.fastAccessDx(3) * dpsidz[n];
1461 KOKKOS_INLINE_FUNCTION
1464 double x[ FunctionCount ] ;
1465 double y[ FunctionCount ] ;
1466 double z[ FunctionCount ] ;
1467 unsigned node_index[ ElemNodeCount ];
1471 scalar_type elem_mat[ FunctionCount ][ FunctionCount ] ;
1474 this->gatherSolution( ielem,
val, node_index,
x,
y, z, elem_res, elem_mat );
1477 computeElementResidualJacobian(
val,
x,
y, z, elem_res, elem_mat );
1480 this->scatterResidual( ielem, node_index, elem_res, elem_mat );
1485 template<
class FiniteElementMeshType ,
class SparseMatrixType
1486 ,
class CoeffFunctionType = ElementComputationConstantCoefficient
1488 class ElementComputation ;
1492 typename ScalarType ,
class CoeffFunctionType >
1493 class ElementComputation
1496 , CoeffFunctionType >
1509 typedef typename sparse_matrix_type::StaticCrsGraphType sparse_graph_type ;
1510 typedef typename sparse_matrix_type::values_type matrix_values_type ;
1511 typedef Kokkos::View< scalar_type* , Kokkos::LayoutLeft, execution_space > vector_type ;
1515 typedef LocalViewTraits< vector_type > local_vector_view_traits;
1516 typedef LocalViewTraits< matrix_values_type> local_matrix_view_traits;
1517 typedef typename local_vector_view_traits::local_view_type local_vector_type;
1518 typedef typename local_matrix_view_traits::local_view_type local_matrix_type;
1519 typedef typename local_vector_view_traits::local_value_type local_scalar_type;
1520 static const bool use_team = local_vector_view_traits::use_team;
1522 static const unsigned SpatialDim = element_data_type::spatial_dimension ;
1523 static const unsigned TensorDim = SpatialDim * SpatialDim ;
1524 static const unsigned ElemNodeCount = element_data_type::element_node_count ;
1525 static const unsigned FunctionCount = element_data_type::function_count ;
1526 static const unsigned IntegrationCount = element_data_type::integration_count ;
1530 typedef typename mesh_type::node_coord_type node_coord_type ;
1531 typedef typename mesh_type::elem_node_type elem_node_type ;
1532 typedef Kokkos::View< scalar_type*[FunctionCount][FunctionCount] , execution_space > elem_matrices_type ;
1533 typedef Kokkos::View< scalar_type*[FunctionCount] , execution_space > elem_vectors_type ;
1535 typedef LocalViewTraits< elem_matrices_type > local_elem_matrices_traits;
1536 typedef LocalViewTraits< elem_vectors_type > local_elem_vectors_traits;
1537 typedef typename local_elem_matrices_traits::local_view_type local_elem_matrices_type;
1538 typedef typename local_elem_vectors_traits::local_view_type local_elem_vectors_type;
1548 const element_data_type elem_data ;
1549 const elem_node_type elem_node_ids ;
1550 const node_coord_type node_coords ;
1551 const elem_graph_type elem_graph ;
1552 const elem_matrices_type elem_jacobians ;
1553 const elem_vectors_type elem_residuals ;
1554 const vector_type solution ;
1555 const vector_type residual ;
1556 const sparse_matrix_type jacobian ;
1557 const CoeffFunctionType coeff_function ;
1560 ElementComputation(
const ElementComputation & rhs )
1562 , elem_node_ids( rhs.elem_node_ids )
1563 , node_coords( rhs.node_coords )
1564 , elem_graph( rhs.elem_graph )
1565 , elem_jacobians( rhs.elem_jacobians )
1566 , elem_residuals( rhs.elem_residuals )
1567 , solution( rhs.solution )
1568 , residual( rhs.residual )
1569 , jacobian( rhs.jacobian )
1570 , coeff_function( rhs.coeff_function )
1571 , dev_config( rhs.dev_config )
1576 ElementComputation(
const mesh_type & arg_mesh ,
1577 const CoeffFunctionType & arg_coeff_function ,
1578 const vector_type & arg_solution ,
1579 const elem_graph_type & arg_elem_graph ,
1580 const sparse_matrix_type & arg_jacobian ,
1581 const vector_type & arg_residual ,
1584 , elem_node_ids( arg_mesh.elem_node() )
1585 , node_coords( arg_mesh.node_coord() )
1586 , elem_graph( arg_elem_graph )
1589 , solution( arg_solution )
1590 , residual( arg_residual )
1591 , jacobian( arg_jacobian )
1592 , coeff_function( arg_coeff_function )
1593 , dev_config( arg_dev_config )
1600 const size_t nelem = elem_node_ids.dimension_0();
1603 const size_t league_size =
1605 Kokkos::TeamPolicy< execution_space > config( league_size, team_size );
1606 parallel_for( config , *
this );
1609 parallel_for( nelem , *
this );
1615 static const unsigned FLOPS_transform_gradients =
1616 FunctionCount * TensorDim * 2 +
1618 FunctionCount * 15 ;
1620 KOKKOS_INLINE_FUNCTION
1621 double transform_gradients(
1622 const double grad[][ FunctionCount ] ,
1628 double dpsidz[] )
const 1630 enum { j11 = 0 , j12 = 1 , j13 = 2 ,
1631 j21 = 3 , j22 = 4 , j23 = 5 ,
1632 j31 = 6 , j32 = 7 , j33 = 8 };
1636 double J[ TensorDim ] = { 0, 0, 0, 0, 0, 0, 0, 0, 0 };
1638 for(
unsigned i = 0; i < FunctionCount ; ++i ) {
1639 const double x1 =
x[i] ;
1640 const double x2 =
y[i] ;
1641 const double x3 = z[i] ;
1643 const double g1 = grad[0][i] ;
1644 const double g2 = grad[1][i] ;
1645 const double g3 = grad[2][i] ;
1662 double invJ[ TensorDim ] = {
1663 static_cast<double>( J[j22] * J[j33] - J[j23] * J[j32] ) ,
1664 static_cast<double>( J[j13] * J[j32] - J[j12] * J[j33] ) ,
1665 static_cast<double>( J[j12] * J[j23] - J[j13] * J[j22] ) ,
1667 static_cast<double>( J[j23] * J[j31] - J[j21] * J[j33] ) ,
1668 static_cast<double>( J[j11] * J[j33] - J[j13] * J[j31] ) ,
1669 static_cast<double>( J[j13] * J[j21] - J[j11] * J[j23] ) ,
1671 static_cast<double>( J[j21] * J[j32] - J[j22] * J[j31] ) ,
1672 static_cast<double>( J[j12] * J[j31] - J[j11] * J[j32] ) ,
1673 static_cast<double>( J[j11] * J[j22] - J[j12] * J[j21] ) };
1675 const double detJ = J[j11] * invJ[j11] +
1676 J[j21] * invJ[j12] +
1677 J[j31] * invJ[j13] ;
1679 const double detJinv = 1.0 / detJ ;
1681 for (
unsigned i = 0 ; i < TensorDim ; ++i ) { invJ[i] *= detJinv ; }
1685 for(
unsigned i = 0; i < FunctionCount ; ++i ) {
1686 const double g0 = grad[0][i];
1687 const double g1 = grad[1][i];
1688 const double g2 = grad[2][i];
1690 dpsidx[i] = g0 * invJ[j11] + g1 * invJ[j12] + g2 * invJ[j13];
1691 dpsidy[i] = g0 * invJ[j21] + g1 * invJ[j22] + g2 * invJ[j23];
1692 dpsidz[i] = g0 * invJ[j31] + g1 * invJ[j32] + g2 * invJ[j33];
1698 KOKKOS_INLINE_FUNCTION
1699 void contributeResidualJacobian(
1700 const local_scalar_type dof_values[] ,
1701 const double dpsidx[] ,
1702 const double dpsidy[] ,
1703 const double dpsidz[] ,
1705 const local_scalar_type coeff_k ,
1706 const double integ_weight ,
1707 const double bases_vals[] ,
1708 local_scalar_type elem_res[] ,
1709 local_scalar_type elem_mat[][ FunctionCount ] )
const 1711 local_scalar_type value_at_pt = 0 ;
1712 local_scalar_type gradx_at_pt = 0 ;
1713 local_scalar_type grady_at_pt = 0 ;
1714 local_scalar_type gradz_at_pt = 0 ;
1716 for (
unsigned m = 0 ; m < FunctionCount ; m++ ) {
1717 value_at_pt += dof_values[m] * bases_vals[m] ;
1718 gradx_at_pt += dof_values[m] * dpsidx[m] ;
1719 grady_at_pt += dof_values[m] * dpsidy[m] ;
1720 gradz_at_pt += dof_values[m] * dpsidz[m] ;
1723 const local_scalar_type k_detJ_weight = coeff_k * detJ * integ_weight ;
1724 const local_scalar_type res_val = value_at_pt * value_at_pt * detJ * integ_weight ;
1725 const local_scalar_type mat_val = 2.0 * value_at_pt * detJ * integ_weight ;
1730 for (
unsigned m = 0; m < FunctionCount; ++m) {
1731 local_scalar_type *
const mat = elem_mat[m] ;
1732 const double bases_val_m = bases_vals[m];
1733 const double dpsidx_m = dpsidx[m] ;
1734 const double dpsidy_m = dpsidy[m] ;
1735 const double dpsidz_m = dpsidz[m] ;
1737 elem_res[m] += k_detJ_weight * ( dpsidx_m * gradx_at_pt +
1738 dpsidy_m * grady_at_pt +
1739 dpsidz_m * gradz_at_pt ) +
1740 res_val * bases_val_m ;
1742 for(
unsigned n = 0; n < FunctionCount; n++) {
1744 mat[n] += k_detJ_weight * ( dpsidx_m * dpsidx[n] +
1745 dpsidy_m * dpsidy[n] +
1746 dpsidz_m * dpsidz[n] ) +
1747 mat_val * bases_val_m * bases_vals[n];
1752 KOKKOS_INLINE_FUNCTION
1753 void operator()(
const typename TeamPolicy< execution_space >::member_type & dev )
const 1756 const unsigned num_ensemble_threads = dev_config.
block_dim.
x ;
1757 const unsigned num_element_threads = dev_config.
block_dim.
y ;
1758 const unsigned element_rank = dev.team_rank() / num_ensemble_threads ;
1759 const unsigned ensemble_rank = dev.team_rank() % num_ensemble_threads ;
1761 const unsigned ielem =
1762 dev.league_rank() * num_element_threads + element_rank;
1764 if (ielem >= elem_node_ids.dimension_0())
1767 (*this)( ielem, ensemble_rank );
1771 KOKKOS_INLINE_FUNCTION
1772 void operator()(
const unsigned ielem ,
1773 const unsigned ensemble_rank = 0 )
const 1775 local_vector_type local_solution =
1776 local_vector_view_traits::create_local_view(solution,
1778 local_vector_type local_residual =
1779 local_vector_view_traits::create_local_view(residual,
1781 local_matrix_type local_jacobian_values =
1782 local_matrix_view_traits::create_local_view(jacobian.values,
1787 double x[ FunctionCount ] ;
1788 double y[ FunctionCount ] ;
1789 double z[ FunctionCount ] ;
1790 local_scalar_type
val[ FunctionCount ] ;
1791 unsigned node_index[ ElemNodeCount ];
1793 for (
unsigned i = 0 ; i < ElemNodeCount ; ++i ) {
1794 const unsigned ni = elem_node_ids( ielem , i );
1796 node_index[i] = ni ;
1798 x[i] = node_coords( ni , 0 );
1799 y[i] = node_coords( ni , 1 );
1800 z[i] = node_coords( ni , 2 );
1802 val[i] = local_solution( ni );
1806 local_scalar_type elem_vec[ FunctionCount ] ;
1807 local_scalar_type elem_mat[ FunctionCount ][ FunctionCount ] ;
1809 for(
unsigned i = 0; i < FunctionCount ; i++ ) {
1811 for(
unsigned j = 0;
j < FunctionCount ;
j++){
1812 elem_mat[i][
j] = 0 ;
1817 for (
unsigned i = 0 ; i < IntegrationCount ; ++i ) {
1818 double dpsidx[ FunctionCount ] ;
1819 double dpsidy[ FunctionCount ] ;
1820 double dpsidz[ FunctionCount ] ;
1822 local_scalar_type coeff_k = 0 ;
1825 double pt[] = {0.0, 0.0, 0.0};
1829 if ( ! coeff_function.is_constant ) {
1830 for (
unsigned j = 0 ;
j < FunctionCount ; ++
j ) {
1831 pt[0] +=
x[
j] * elem_data.values[i][
j] ;
1832 pt[1] +=
y[
j] * elem_data.values[i][
j] ;
1833 pt[2] += z[
j] * elem_data.values[i][
j] ;
1838 coeff_k = coeff_function(pt, ensemble_rank);
1842 transform_gradients( elem_data.gradients[i] ,
x ,
y , z ,
1843 dpsidx , dpsidy , dpsidz );
1845 contributeResidualJacobian(
val , dpsidx , dpsidy , dpsidz ,
1847 elem_data.weights[i] ,
1848 elem_data.values[i] ,
1849 elem_vec , elem_mat );
1852 for(
unsigned i = 0 ; i < FunctionCount ; i++ ) {
1853 const unsigned row = node_index[i] ;
1854 if ( row < residual.dimension_0() ) {
1855 atomic_add( & local_residual( row ) , elem_vec[i] );
1857 for(
unsigned j = 0 ;
j < FunctionCount ;
j++ ) {
1858 const unsigned entry = elem_graph( ielem , i ,
j );
1859 if ( entry != ~0u ) {
1860 atomic_add( & local_jacobian_values( entry ) , elem_mat[i][
j] );
1871 template<
class FixtureType ,
class SparseMatrixType >
1875 typename ScalarType >
1900 static const bool use_team = local_vector_view_traits::use_team;
1925 const unsigned arg_bc_plane ,
1929 : node_coords( arg_mesh.node_coord() )
1930 , solution( arg_solution )
1931 , jacobian( arg_jacobian )
1932 , residual( arg_residual )
1933 , bc_lower_value( arg_bc_lower_value )
1934 , bc_upper_value( arg_bc_upper_value )
1937 , bc_plane( arg_bc_plane )
1938 , node_count( arg_mesh.node_count_owned() )
1940 , dev_config( arg_dev_config )
1942 parallel_for( node_count , *
this );
1950 const size_t league_size =
1952 Kokkos::TeamPolicy< execution_space > config( league_size, team_size );
1953 parallel_for( config , *
this );
1956 parallel_for( node_count , *
this );
1961 KOKKOS_INLINE_FUNCTION
1962 void operator()(
const typename TeamPolicy< execution_space >::member_type & dev )
const 1965 const unsigned num_ensemble_threads = dev_config.
block_dim.
x ;
1966 const unsigned num_node_threads = dev_config.
block_dim.
y ;
1967 const unsigned node_rank = dev.team_rank() / num_ensemble_threads ;
1968 const unsigned ensemble_rank = dev.team_rank() % num_ensemble_threads ;
1970 const unsigned inode = dev.league_rank() * num_node_threads + node_rank;
1972 if (inode >= node_count)
1975 (*this)( inode, ensemble_rank );
1979 KOKKOS_INLINE_FUNCTION
1981 const unsigned ensemble_rank = 0)
const 1984 local_vector_view_traits::create_local_view(residual,
1987 local_matrix_view_traits::create_local_view(jacobian.
values,
1995 const unsigned iBeg = jacobian.
graph.row_map[inode];
1996 const unsigned iEnd = jacobian.
graph.row_map[inode+1];
1999 const bool bc_lower = c <= bc_lower_limit ;
2000 const bool bc_upper = bc_upper_limit <= c ;
2003 solution(inode) = bc_lower ? bc_lower_value : (
2004 bc_upper ? bc_upper_value : 0 );
2007 if ( bc_lower || bc_upper ) {
2009 local_residual(inode) = 0 ;
2014 for(
unsigned i = iBeg ; i < iEnd ; ++i ) {
2015 local_jacobian_values(i) =
int(inode) ==
int(jacobian.
graph.entries(i)) ? 1 : 0 ;
2023 for(
unsigned i = iBeg ; i < iEnd ; ++i ) {
2024 const unsigned cnode = jacobian.
graph.entries(i) ;
2027 if ( ( cc <= bc_lower_limit ) || ( bc_upper_limit <= cc ) ) {
2028 local_jacobian_values(i) = 0 ;
2043 #if defined( __CUDACC__ )
static KOKKOS_INLINE_FUNCTION local_view_type create_local_view(const view_type &v, const unsigned local_rank)
KOKKOS_INLINE_FUNCTION void join(volatile unsigned &update, const volatile unsigned &input) const
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION void operator()(const unsigned iwork) const
KOKKOS_INLINE_FUNCTION void computeElementResidualJacobian(const scalar_type dof_values[], const double x[], const double y[], const double z[], scalar_type elem_res[], scalar_type elem_mat[][FunctionCount]) const
LocalViewTraits< vector_type > local_vector_view_traits
ExponentialKLCoefficient(const ExponentialKLCoefficient &rhs)
const scalar_coord_type bc_upper_limit
static const unsigned SpatialDim
Kokkos::View< scalar_type *, execution_space > vector_type
const Kokkos::Example::FENL::DeviceConfig dev_config
KOKKOS_INLINE_FUNCTION void setRandomVariables(const RandomVariableView &rv)
LocalViewTraits< matrix_values_type > local_matrix_view_traits
Sacado::Fad::SLFad< scalar_type, FunctionCount > fad_scalar_type
const unsigned node_count
ElementComputation(const ElementComputation &rhs)
KOKKOS_INLINE_FUNCTION void operator()(const unsigned ielem) const
LocalViewTraits< RandomVariableView > local_rv_view_traits
const elem_vectors_type elem_residuals
const vector_type solution
static const bool use_team
KOKKOS_INLINE_FUNCTION RandomVariableView getRandomVariables() const
Kokkos::View< scalar_type *[FunctionCount], execution_space > elem_vectors_type
static const unsigned element_node_count
CrsMatrix< ScalarType, ExecutionSpace > sparse_matrix_type
const elem_matrices_type elem_jacobians
view_type::value_type local_value_type
static const unsigned function_count
static void eval(Kokkos::Example::FENL::DeviceConfig &dev_config_elem, Kokkos::Example::FENL::DeviceConfig &dev_config_bc)
ElementComputation(const typename base_type::mesh_type &arg_mesh, const CoeffFunctionType &arg_coeff_function, const typename base_type::vector_type &arg_solution, const typename base_type::elem_graph_type &arg_elem_graph, const typename base_type::sparse_matrix_type &arg_jacobian, const typename base_type::vector_type &arg_residual, const Kokkos::Example::FENL::DeviceConfig arg_dev_config)
base_type::execution_space execution_space
ElementComputationBase< ExecutionSpace, Order, CoordinateMap, ScalarType > base_type
Kokkos::View< const double *[SpaceDim], Device > node_coord_type
Kokkos::DefaultExecutionSpace execution_space
Kokkos::View< const unsigned *[ElemNode], Device > elem_node_type
static const unsigned ElemNodeCount
Sacado::Fad::SLFad< scalar_type, 4 > fad_scalar_type
static void eval(Kokkos::Example::FENL::DeviceConfig &dev_config_elem, Kokkos::Example::FENL::DeviceConfig &dev_config_bc)
KOKKOS_INLINE_FUNCTION double operator()(double pt[], unsigned ensemble_rank) const
const Kokkos::Example::FENL::DeviceConfig dev_config
KOKKOS_INLINE_FUNCTION void gatherSolution(const unsigned ielem, fad_scalar_type val[], unsigned node_index[], double x[], double y[], double z[], fad_scalar_type res[]) const
const CoeffFunctionType coeff_function
ExecutionSpace execution_space
ElementComputation(const ElementComputation &rhs)
Kokkos::View< scalar_type *, Kokkos::LayoutLeft, execution_space > vector_type
KOKKOS_INLINE_FUNCTION void atomic_add(volatile Sacado::UQ::PCE< Storage > *const dest, const Sacado::UQ::PCE< Storage > &src)
const vector_type residual
const bc_scalar_type bc_lower_value
ElementComputation(const typename base_type::mesh_type &arg_mesh, const CoeffFunctionType &arg_coeff_function, const typename base_type::vector_type &arg_solution, const typename base_type::elem_graph_type &arg_elem_graph, const typename base_type::sparse_matrix_type &arg_jacobian, const typename base_type::vector_type &arg_residual, const Kokkos::Example::FENL::DeviceConfig arg_dev_config)
Kokkos::View< Scalar *, Kokkos::LayoutLeft, Device > RandomVariableView
KOKKOS_INLINE_FUNCTION local_scalar_type operator()(const MeshScalar point[], const size_type ensemble_rank) const
const elem_node_type elem_node_ids
const vector_type residual
const unsigned node_count
const bc_scalar_type bc_upper_value
KOKKOS_INLINE_FUNCTION void fill_graph_entries(const unsigned iset) const
base_type::scalar_type scalar_type
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
base_type::execution_space execution_space
NodeNodeGraph< elem_node_type, sparse_graph_type, ElemNodeCount >::ElemGraphType elem_graph_type
const vector_type solution
KOKKOS_INLINE_FUNCTION void operator()(const typename TeamPolicy< execution_space >::member_type &dev) const
const MeshScalar m_variance
local_rv_view_traits::local_view_type local_rv_view_type
ElementComputation< Kokkos::Example::BoxElemFixture< ExecutionSpace, Order, CoordinateMap >, CrsMatrix< ScalarType, ExecutionSpace >, Analytic, CoeffFunctionType > base_type
mesh_type::node_coord_type node_coord_type
ElemNodeIdView::execution_space execution_space
node_coord_type::value_type scalar_coord_type
static const unsigned TensorDim
KOKKOS_INLINE_FUNCTION void computeElementResidual(const scalar_type dof_values[], const double x[], const double y[], const double z[], local_scalar_type elem_res[]) const
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
ElementComputationConstantCoefficient(const ElementComputationConstantCoefficient &rhs)
RandomVariableView::size_type size_type
KOKKOS_INLINE_FUNCTION void gatherSolution(const unsigned ielem, scalar_type val[], unsigned node_index[], double x[], double y[], double z[], fad_scalar_type res[]) const
ExecutionSpace execution_space
double sort_graph_entries
size_t num_threads_per_block
double fill_graph_entries
ElementComputationBase(const ElementComputationBase &rhs)
KOKKOS_INLINE_FUNCTION void operator()(const unsigned ielem) const
KOKKOS_INLINE_FUNCTION void operator()(const unsigned irow, unsigned &update, const bool final) const
Dim3(const size_t x_, const size_t y_=1, const size_t z_=1)
KOKKOS_INLINE_FUNCTION void scatterResidual(const unsigned ielem, const unsigned node_index[], fad_scalar_type res[]) const
Kokkos::Example::HexElement_Data< mesh_type::ElemNode > element_data_type
CrsGraphType::row_map_type::non_const_type RowMapType
base_type::scalar_type scalar_type
const element_data_type elem_data
const MeshScalar m_corr_len
KOKKOS_INLINE_FUNCTION void sort_graph_entries(const unsigned irow) const
ElementComputationBase(const mesh_type &arg_mesh, const vector_type &arg_solution, const elem_graph_type &arg_elem_graph, const sparse_matrix_type &arg_jacobian, const vector_type &arg_residual)
ElementComputation(const ElementComputation &rhs)
base_type::execution_space execution_space
CrsMatrix(const StaticCrsGraphType &arg_graph)
const sparse_matrix_type jacobian
mesh_type::elem_node_type elem_node_type
Kokkos::Example::FENL::CrsMatrix< ScalarType, ExecutionSpace > sparse_matrix_type
void deep_copy(const Stokhos::CrsMatrix< ValueType, DstDevice, Layout > &dst, const Stokhos::CrsMatrix< ValueType, SrcDevice, Layout > &src)
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_set(const unsigned ielem) const
Kokkos::StaticCrsGraph< unsigned, Space, void, unsigned > StaticCrsGraphType
KOKKOS_INLINE_FUNCTION void gatherSolution(const unsigned ielem, scalar_type val[], unsigned node_index[], double x[], double y[], double z[], scalar_type res[], scalar_type mat[][FunctionCount]) const
base_type::execution_space execution_space
base_type::scalar_type scalar_type
const CoeffFunctionType coeff_function
const scalar_coord_type bc_lower_limit
KOKKOS_INLINE_FUNCTION void computeElementResidualJacobian(const scalar_type dof_values[], const double x[], const double y[], const double z[], scalar_type elem_res[], scalar_type elem_mat[][FunctionCount]) const
const ElemNodeIdView elem_node_id
static const unsigned integration_count
base_type::scalar_type scalar_type
Kokkos::Example::BoxElemFixture< ExecutionSpace, Order, CoordinateMap > mesh_type
mesh_type::node_coord_type node_coord_type
const elem_graph_type elem_graph
KOKKOS_INLINE_FUNCTION void operator()(const unsigned ielem) const
sparse_matrix_type::StaticCrsGraphType sparse_graph_type
const Kokkos::Example::FENL::DeviceConfig dev_config
Kokkos::Example::BoxElemFixture< ExecutionSpace, Order, CoordinateMap > mesh_type
const view_type & local_view_type
NodeNodeGraph(const ElemNodeIdView &arg_elem_node_id, const unsigned arg_node_count, Times &results)
const sparse_matrix_type jacobian
const node_coord_type node_coords
static const unsigned FunctionCount
ElementComputation< Kokkos::Example::BoxElemFixture< ExecutionSpace, Order, CoordinateMap >, CrsMatrix< ScalarType, ExecutionSpace >, FadElement, CoeffFunctionType > base_type
KOKKOS_INLINE_FUNCTION void fill_elem_graph_map(const unsigned ielem) const
sparse_matrix_type::StaticCrsGraphType sparse_graph_type
Kokkos::View< unsigned, execution_space > UnsignedValue
const node_coord_type node_coords
KOKKOS_INLINE_FUNCTION void scatterResidual(const unsigned ielem, const unsigned node_index[], const scalar_type res[], const scalar_type mat[][FunctionCount]) const
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)
ExponentialKLCoefficient< T1, T2, T3 > type
ElementComputation(const typename base_type::mesh_type &arg_mesh, const CoeffFunctionType &arg_coeff_function, const typename base_type::vector_type &arg_solution, const typename base_type::elem_graph_type &arg_elem_graph, const typename base_type::sparse_matrix_type &arg_jacobian, const typename base_type::vector_type &arg_residual, const Kokkos::Example::FENL::DeviceConfig arg_dev_config)
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)
ElementComputationBase< ExecutionSpace, Order, CoordinateMap, ScalarType > base_type
KOKKOS_INLINE_FUNCTION PCE< Storage > sin(const PCE< Storage > &a)
ElementComputation(const typename base_type::mesh_type &arg_mesh, const CoeffFunctionType &arg_coeff_function, const typename base_type::vector_type &arg_solution, const typename base_type::elem_graph_type &arg_elem_graph, const typename base_type::sparse_matrix_type &arg_jacobian, const typename base_type::vector_type &arg_residual, const Kokkos::Example::FENL::DeviceConfig arg_dev_config)
View< ValueType *, Space > values_type
KOKKOS_INLINE_FUNCTION void operator()(const unsigned ielem) const
StorageType::execution_space execution_space
local_matrix_view_traits::local_view_type local_matrix_type
KOKKOS_INLINE_FUNCTION void computeElementResidual(const local_scalar_type dof_values[], const double x[], const double y[], const double z[], local_scalar_type elem_res[]) const
double fill_element_graph
Kokkos::View< scalar_type *[FunctionCount][FunctionCount], execution_space > elem_matrices_type
const unsigned VectorSize
KOKKOS_INLINE_FUNCTION void operator()(const unsigned inode, const unsigned ensemble_rank=0) const
ElementComputationConstantCoefficient(const double val)
Kokkos::View< unsigned *[ElemNode][ElemNode], execution_space > ElemGraphType
Sacado::Fad::SLFad< scalar_type, FunctionCount > fad_scalar_type
Generate a distributed unstructured finite element mesh from a partitioned NX*NY*NZ box of elements...
Stokhos::KL::ExponentialRandomField< MeshScalar, Device > rf_type
sparse_matrix_type::values_type matrix_values_type
ExponentialKLCoefficient(const MeshScalar mean, const MeshScalar variance, const MeshScalar correlation_length, const size_type num_rv)
local_vector_view_traits::local_view_type local_vector_type
KOKKOS_INLINE_FUNCTION void init(unsigned &update) const
ElementComputation(const ElementComputation &rhs)
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType ValueType * y
void update(const ValueType &alpha, VectorType &x, const ValueType &beta, const VectorType &y)
Kokkos::UnorderedMap< key_type, void, execution_space > SetType
pair< unsigned, unsigned > key_type
static const unsigned IntegrationCount
static const unsigned spatial_dimension