44 #ifndef KOKKOS_EXAMPLE_FENLFUNCTORS_HPP 45 #define KOKKOS_EXAMPLE_FENLFUNCTORS_HPP 56 #include <Kokkos_Core.hpp> 57 #include <Kokkos_Pair.hpp> 58 #include <Kokkos_UnorderedMap.hpp> 59 #include <Kokkos_StaticCrsGraph.hpp> 61 #include <impl/Kokkos_Timer.hpp> 75 template<
typename ValueType ,
class Space >
87 ,
coeff(
"crs_matrix_coeff" , arg_graph.entries.dimension_0() )
91 template<
class ElemNodeIdView ,
class CrsGraphType ,
unsigned ElemNode >
98 typedef Kokkos::UnorderedMap< key_type, void , execution_space >
SetType ;
99 typedef typename CrsGraphType::row_map_type::non_const_type
RowMapType ;
103 typedef Kokkos::View< unsigned*[ElemNode][ElemNode] , execution_space >
107 struct TagScanNodeCount {};
108 struct TagFillGraphEntries {};
109 struct TagSortGraphEntries {};
110 struct TagFillElementGraph {};
144 const unsigned arg_node_count,
160 Kokkos::Impl::Timer wall_clock ;
166 size_t set_capacity = (28ull *
node_count) / 2;
167 unsigned failed_insert_count = 0 ;
178 Kokkos::parallel_reduce( Kokkos::RangePolicy<execution_space,TagFillNodeSet>(0,
elem_node_id.dimension_0())
180 , failed_insert_count );
182 }
while ( failed_insert_count );
184 execution_space::fence();
186 results.fill_node_set = wall_clock.seconds();
200 unsigned graph_entry_count = 0 ;
202 Kokkos::deep_copy( graph_entry_count ,
row_total );
206 graph.entries =
typename CrsGraphType::entries_type(
"graph_entries" , graph_entry_count );
211 execution_space::fence();
212 results.scan_node_count = wall_clock.seconds();
218 execution_space::fence();
219 results.fill_graph_entries = wall_clock.seconds();
235 execution_space::fence();
236 results.sort_graph_entries = wall_clock.seconds();
243 Kokkos::parallel_for(
elem_node_id.dimension_0() , *this );
245 execution_space::fence();
246 results.fill_element_graph = wall_clock.seconds();
256 for (
unsigned row_local_node = 0 ; row_local_node <
elem_node_id.dimension_1() ; ++row_local_node ) {
258 const unsigned row_node =
elem_node_id( ielem , row_local_node );
260 for (
unsigned col_local_node = row_local_node ; col_local_node <
elem_node_id.dimension_1() ; ++col_local_node ) {
262 const unsigned col_node =
elem_node_id( ielem , col_local_node );
268 const key_type key = (row_node < col_node) ? make_pair( row_node, col_node ) : make_pair( col_node, row_node ) ;
270 const typename SetType::insert_result result =
node_node_set.insert( key );
273 if ( result.success() ) {
276 if ( row_node <
row_count.dimension_0() ) { atomic_fetch_add( &
row_count( row_node ) , 1 ); }
279 if ( col_node <
row_count.dimension_0() && col_node != row_node ) { atomic_fetch_add( &
row_count( col_node ) , 1 ); }
281 else if ( result.failed() ) {
296 const unsigned row_node = key.first ;
297 const unsigned col_node = key.second ;
299 if ( row_node <
row_count.dimension_0() ) {
300 const unsigned offset =
graph.row_map( row_node ) + atomic_fetch_add( &
row_count( row_node ) , 1 );
301 graph.entries( offset ) = col_node ;
304 if ( col_node <
row_count.dimension_0() && col_node != row_node ) {
305 const unsigned offset =
graph.row_map( col_node ) + atomic_fetch_add( &
row_count( col_node ) , 1 );
306 graph.entries( offset ) = row_node ;
314 const unsigned row_beg =
graph.row_map( irow );
315 const unsigned row_end =
graph.row_map( irow + 1 );
316 for (
unsigned i = row_beg + 1 ; i < row_end ; ++i ) {
317 const unsigned col =
graph.entries(i);
319 for ( ; row_beg < j && col <
graph.entries(j-1) ; --j ) {
322 graph.entries(j) = col ;
329 for (
unsigned row_local_node = 0 ; row_local_node <
elem_node_id.dimension_1() ; ++row_local_node ) {
331 const unsigned row_node =
elem_node_id( ielem , row_local_node );
333 for (
unsigned col_local_node = 0 ; col_local_node <
elem_node_id.dimension_1() ; ++col_local_node ) {
335 const unsigned col_node =
elem_node_id( ielem , col_local_node );
337 unsigned entry = ~0u ;
339 if ( row_node + 1 <
graph.row_map.dimension_0() ) {
341 const unsigned entry_end =
graph.row_map( row_node + 1 );
343 entry =
graph.row_map( row_node );
345 for ( ; entry < entry_end &&
graph.entries(entry) != col_node ; ++entry );
347 if ( entry == entry_end ) entry = ~0u ;
350 elem_graph( ielem , row_local_node , col_local_node ) = entry ;
380 void operator()(
const unsigned irow ,
unsigned & update ,
const bool final )
const 383 if (
final ) {
row_map( irow ) = update ; }
388 if ( irow + 1 ==
row_count.dimension_0() ) {
401 ,
volatile unsigned & update
402 ,
volatile const unsigned & input )
const { update += input ; }
406 void init(
unsigned & update )
const { update = 0 ; }
409 void join(
volatile unsigned & update
410 ,
volatile const unsigned & input )
const { update += input ; }
427 class CoordinateMap ,
typename ScalarType >
428 class ElementComputationBase
442 typedef Kokkos::View< scalar_type* , Kokkos::LayoutLeft, execution_space >
vector_type ;
456 typedef Kokkos::View< scalar_type*[FunctionCount][FunctionCount] , execution_space >
elem_matrices_type ;
515 double dpsidz[] )
const 517 enum { j11 = 0 , j12 = 1 , j13 = 2 ,
518 j21 = 3 , j22 = 4 , j23 = 5 ,
519 j31 = 6 , j32 = 7 , j33 = 8 };
523 double J[
TensorDim ] = { 0, 0, 0, 0, 0, 0, 0, 0, 0 };
526 const double x1 = x[i] ;
527 const double x2 = y[i] ;
528 const double x3 = z[i] ;
530 const double g1 = grad[0][i] ;
531 const double g2 = grad[1][i] ;
532 const double g3 = grad[2][i] ;
550 static_cast<double>( J[j22] * J[j33] - J[j23] * J[j32] ) ,
551 static_cast<double>( J[j13] * J[j32] - J[j12] * J[j33] ) ,
552 static_cast<double>( J[j12] * J[j23] - J[j13] * J[j22] ) ,
554 static_cast<double>( J[j23] * J[j31] - J[j21] * J[j33] ) ,
555 static_cast<double>( J[j11] * J[j33] - J[j13] * J[j31] ) ,
556 static_cast<double>( J[j13] * J[j21] - J[j11] * J[j23] ) ,
558 static_cast<double>( J[j21] * J[j32] - J[j22] * J[j31] ) ,
559 static_cast<double>( J[j12] * J[j31] - J[j11] * J[j32] ) ,
560 static_cast<double>( J[j11] * J[j22] - J[j12] * J[j21] ) };
562 const double detJ = J[j11] * invJ[j11] +
566 const double detJinv = 1.0 / detJ ;
568 for (
unsigned i = 0 ; i <
TensorDim ; ++i ) { invJ[i] *= detJinv ; }
573 const double g0 = grad[0][i];
574 const double g1 = grad[1][i];
575 const double g2 = grad[2][i];
577 dpsidx[i] = g0 * invJ[j11] + g1 * invJ[j12] + g2 * invJ[j13];
578 dpsidy[i] = g0 * invJ[j21] + g1 * invJ[j22] + g2 * invJ[j23];
579 dpsidz[i] = g0 * invJ[j31] + g1 * invJ[j32] + g2 * invJ[j33];
594 template<
class FiniteElementMeshType ,
595 class SparseMatrixType ,
598 class ElementComputation ;
601 class CoordinateMap ,
typename ScalarType >
602 class ElementComputation
604 CrsMatrix< ScalarType , ExecutionSpace > ,
606 public ElementComputationBase<ExecutionSpace, Order, CoordinateMap,
610 typedef ElementComputationBase<ExecutionSpace, Order, CoordinateMap,
616 static const unsigned FunctionCount = base_type::FunctionCount;
617 static const unsigned IntegrationCount = base_type::IntegrationCount;
618 static const unsigned ElemNodeCount = base_type::ElemNodeCount;
620 typedef Kokkos::View<scalar_type[FunctionCount],Kokkos::LayoutRight,execution_space,Kokkos::MemoryUnmanaged>
elem_vec_type;
621 typedef Kokkos::View<scalar_type[FunctionCount][FunctionCount],Kokkos::LayoutRight,execution_space,Kokkos::MemoryUnmanaged>
elem_mat_type;
631 base_type(arg_mesh, arg_solution, arg_elem_graph,
632 arg_jacobian, arg_residual) {}
638 const size_t nelem = this->elem_node_ids.dimension_0();
639 parallel_for( nelem , *
this );
645 unsigned node_index[],
646 double x[],
double y[],
double z[],
650 for (
unsigned i = 0 ; i < ElemNodeCount ; ++i ) {
651 const unsigned ni = this->elem_node_ids( ielem , i );
655 x[i] = this->node_coords( ni , 0 );
656 y[i] = this->node_coords( ni , 1 );
657 z[i] = this->node_coords( ni , 2 );
659 val(i) = this->solution( ni ) ;
662 for(
unsigned j = 0; j < FunctionCount ; j++){
670 const unsigned node_index[],
674 for(
unsigned i = 0 ; i < FunctionCount ; i++ ) {
675 const unsigned row = node_index[i] ;
676 if ( row < this->residual.dimension_0() ) {
677 atomic_add( & this->residual( row ) , res(i) );
679 for(
unsigned j = 0 ; j < FunctionCount ; j++ ) {
680 const unsigned entry = this->elem_graph( ielem , i , j );
681 if ( entry != ~0u ) {
682 atomic_add( & this->jacobian.coeff( entry ) , mat(i,j) );
698 double coeff_k = 3.456;
699 double coeff_src = 1.234;
700 double advection[] = { 1.1, 1.2, 1.3 };
701 double dpsidx[ FunctionCount ] ;
702 double dpsidy[ FunctionCount ] ;
703 double dpsidz[ FunctionCount ] ;
704 for (
unsigned i = 0 ; i < IntegrationCount ; ++i ) {
706 const double integ_weight = this->elem_data.weights[i];
707 const double* bases_vals = this->elem_data.values[i];
709 this->transform_gradients( this->elem_data.gradients[i] ,
711 dpsidx , dpsidy , dpsidz );
712 const double detJ_weight = detJ * integ_weight;
713 const double detJ_weight_coeff_k = detJ_weight * coeff_k;
719 for (
unsigned m = 0 ; m < FunctionCount ; m++ ) {
720 value_at_pt += dof_values(m) * bases_vals[m] ;
721 gradx_at_pt += dof_values(m) * dpsidx[m] ;
722 grady_at_pt += dof_values(m) * dpsidy[m] ;
723 gradz_at_pt += dof_values(m) * dpsidz[m] ;
727 coeff_src * value_at_pt * value_at_pt ;
729 2.0 * coeff_src * value_at_pt ;
736 advection_x*gradx_at_pt +
737 advection_y*grady_at_pt +
738 advection_z*gradz_at_pt ;
740 for (
unsigned m = 0; m < FunctionCount; ++m) {
741 const double bases_val_m = bases_vals[m] * detJ_weight ;
742 const double dpsidx_m = dpsidx[m] ;
743 const double dpsidy_m = dpsidy[m] ;
744 const double dpsidz_m = dpsidz[m] ;
747 detJ_weight_coeff_k * ( dpsidx_m * gradx_at_pt +
748 dpsidy_m * grady_at_pt +
749 dpsidz_m * gradz_at_pt ) +
750 bases_val_m * ( advection_term + source_term ) ;
752 for(
unsigned n = 0;
n < FunctionCount;
n++) {
753 const double dpsidx_n = dpsidx[
n] ;
754 const double dpsidy_n = dpsidy[
n] ;
755 const double dpsidz_n = dpsidz[
n] ;
757 detJ_weight_coeff_k * ( dpsidx_m * dpsidx_n +
758 dpsidy_m * dpsidy_n +
759 dpsidz_m * dpsidz_n ) +
760 bases_val_m * ( advection_x * dpsidx_n +
761 advection_y * dpsidy_n +
762 advection_z * dpsidz_n +
763 source_deriv * bases_vals[
n] ) ;
772 double x[ FunctionCount ] ;
773 double y[ FunctionCount ] ;
774 double z[ FunctionCount ] ;
775 unsigned node_index[ ElemNodeCount ];
777 scalar_type local_val[FunctionCount], local_res[FunctionCount], local_mat[FunctionCount*FunctionCount];
780 elem_mat_type elem_mat(local_mat, FunctionCount, FunctionCount ) ;
783 gatherSolution(ielem,
val, node_index, x, y, z, elem_res, elem_mat);
786 computeElementResidualJacobian(
val, x, y, z, elem_res , elem_mat );
789 scatterResidual( ielem, node_index, elem_res, elem_mat );
794 class CoordinateMap ,
typename ScalarType >
795 class ElementComputation
797 CrsMatrix< ScalarType , ExecutionSpace > ,
798 FadElement > :
public ElementComputationBase<ExecutionSpace, Order, CoordinateMap,
802 typedef ElementComputationBase<ExecutionSpace, Order, CoordinateMap,
808 static const unsigned FunctionCount = base_type::FunctionCount;
809 static const unsigned IntegrationCount = base_type::IntegrationCount;
810 static const unsigned ElemNodeCount = base_type::ElemNodeCount;
815 typedef Kokkos::View<fad_scalar_type*,Kokkos::LayoutRight,execution_space,Kokkos::MemoryUnmanaged>
elem_vec_type;
825 base_type(arg_mesh, arg_solution, arg_elem_graph,
826 arg_jacobian, arg_residual) {}
832 const size_t nelem = this->elem_node_ids.dimension_0();
833 parallel_for( nelem , *
this );
839 unsigned node_index[],
840 double x[],
double y[],
double z[],
843 for (
unsigned i = 0 ; i < ElemNodeCount ; ++i ) {
844 const unsigned ni = this->elem_node_ids( ielem , i );
848 x[i] = this->node_coords( ni , 0 );
849 y[i] = this->node_coords( ni , 1 );
850 z[i] = this->node_coords( ni , 2 );
852 val(i).val() = this->solution( ni );
853 val(i).diff( i, FunctionCount );
860 const unsigned node_index[],
863 for(
unsigned i = 0 ; i < FunctionCount ; i++ ) {
864 const unsigned row = node_index[i] ;
865 if ( row < this->residual.dimension_0() ) {
866 atomic_add( & this->residual( row ) , res(i).
val() );
868 for(
unsigned j = 0 ; j < FunctionCount ; j++ ) {
869 const unsigned entry = this->elem_graph( ielem , i , j );
870 if ( entry != ~0u ) {
871 atomic_add( & this->jacobian.coeff( entry ) ,
872 res(i).fastAccessDx(j) );
886 double coeff_k = 3.456;
887 double coeff_src = 1.234;
888 double advection[] = { 1.1, 1.2, 1.3 };
889 double dpsidx[ FunctionCount ] ;
890 double dpsidy[ FunctionCount ] ;
891 double dpsidz[ FunctionCount ] ;
892 for (
unsigned i = 0 ; i < IntegrationCount ; ++i ) {
894 const double integ_weight = this->elem_data.weights[i];
895 const double* bases_vals = this->elem_data.values[i];
897 this->transform_gradients( this->elem_data.gradients[i] ,
899 dpsidx , dpsidy , dpsidz );
900 const double detJ_weight = detJ * integ_weight;
901 const double detJ_weight_coeff_k = detJ_weight * coeff_k;
907 for (
unsigned m = 0 ; m < FunctionCount ; m++ ) {
908 value_at_pt += dof_values(m) * bases_vals[m] ;
909 gradx_at_pt += dof_values(m) * dpsidx[m] ;
910 grady_at_pt += dof_values(m) * dpsidy[m] ;
911 gradz_at_pt += dof_values(m) * dpsidz[m] ;
915 coeff_src * value_at_pt * value_at_pt ;
918 advection[0]*gradx_at_pt +
919 advection[1]*grady_at_pt +
920 advection[2]*gradz_at_pt;
922 for (
unsigned m = 0; m < FunctionCount; ++m) {
923 const double bases_val_m = bases_vals[m] * detJ_weight ;
924 const double dpsidx_m = dpsidx[m] ;
925 const double dpsidy_m = dpsidy[m] ;
926 const double dpsidz_m = dpsidz[m] ;
929 detJ_weight_coeff_k * ( dpsidx_m * gradx_at_pt +
930 dpsidy_m * grady_at_pt +
931 dpsidz_m * gradz_at_pt ) +
932 bases_val_m * ( advection_term + source_term ) ;
940 double x[ FunctionCount ] ;
941 double y[ FunctionCount ] ;
942 double z[ FunctionCount ] ;
943 unsigned node_index[ ElemNodeCount ];
945 scalar_type local_val[FunctionCount*(FunctionCount+1)], local_res[FunctionCount*(FunctionCount+1)];
947 elem_vec_type elem_res(local_res, FunctionCount, FunctionCount+1) ;
950 gatherSolution( ielem,
val, node_index, x, y, z, elem_res );
953 computeElementResidual(
val, x, y, z, elem_res );
956 scatterResidual( ielem, node_index, elem_res );
961 class CoordinateMap ,
typename ScalarType >
962 class ElementComputation
964 CrsMatrix< ScalarType , ExecutionSpace > ,
966 public ElementComputation< Kokkos::Example::BoxElemFixture< ExecutionSpace , Order , CoordinateMap > ,
967 CrsMatrix< ScalarType , ExecutionSpace > ,
971 typedef ElementComputation< Kokkos::Example::BoxElemFixture< ExecutionSpace , Order , CoordinateMap > ,
972 CrsMatrix< ScalarType , ExecutionSpace > ,
978 static const unsigned FunctionCount = base_type::FunctionCount;
979 static const unsigned IntegrationCount = base_type::IntegrationCount;
980 static const unsigned ElemNodeCount = base_type::ElemNodeCount;
984 typedef Kokkos::View<scalar_type[FunctionCount],Kokkos::LayoutRight,execution_space,Kokkos::MemoryUnmanaged>
scalar_elem_vec_type;
995 base_type(arg_mesh, arg_solution, arg_elem_graph,
996 arg_jacobian, arg_residual) {}
1002 const size_t nelem = this->elem_node_ids.dimension_0();
1003 parallel_for( nelem , *
this );
1009 unsigned node_index[],
1010 double x[],
double y[],
double z[],
1013 for (
unsigned i = 0 ; i < ElemNodeCount ; ++i ) {
1014 const unsigned ni = this->elem_node_ids( ielem , i );
1016 node_index[i] = ni ;
1018 x[i] = this->node_coords( ni , 0 );
1019 y[i] = this->node_coords( ni , 1 );
1020 z[i] = this->node_coords( ni , 2 );
1022 val(i) = this->solution( ni );
1034 double coeff_k = 3.456;
1035 double coeff_src = 1.234;
1036 double advection[] = { 1.1, 1.2, 1.3 };
1037 double dpsidx[ FunctionCount ] ;
1038 double dpsidy[ FunctionCount ] ;
1039 double dpsidz[ FunctionCount ] ;
1040 for (
unsigned i = 0 ; i < IntegrationCount ; ++i ) {
1042 const double integ_weight = this->elem_data.weights[i];
1043 const double* bases_vals = this->elem_data.values[i];
1045 this->transform_gradients( this->elem_data.gradients[i] ,
1047 dpsidx , dpsidy , dpsidz );
1048 const double detJ_weight = detJ * integ_weight;
1049 const double detJ_weight_coeff_k = detJ_weight * coeff_k;
1055 for (
unsigned m = 0 ; m < FunctionCount ; m++ ) {
1056 value_at_pt.val() += dof_values(m) * bases_vals[m] ;
1057 value_at_pt.fastAccessDx(m) = bases_vals[m] ;
1059 gradx_at_pt.val() += dof_values(m) * dpsidx[m] ;
1060 gradx_at_pt.fastAccessDx(m) = dpsidx[m] ;
1062 grady_at_pt.val() += dof_values(m) * dpsidy[m] ;
1063 grady_at_pt.fastAccessDx(m) = dpsidy[m] ;
1065 gradz_at_pt.val() += dof_values(m) * dpsidz[m] ;
1066 gradz_at_pt.fastAccessDx(m) = dpsidz[m] ;
1070 coeff_src * value_at_pt * value_at_pt ;
1073 advection[0]*gradx_at_pt +
1074 advection[1]*grady_at_pt +
1075 advection[2]*gradz_at_pt;
1077 for (
unsigned m = 0; m < FunctionCount; ++m) {
1078 const double bases_val_m = bases_vals[m] * detJ_weight ;
1079 const double dpsidx_m = dpsidx[m] ;
1080 const double dpsidy_m = dpsidy[m] ;
1081 const double dpsidz_m = dpsidz[m] ;
1084 detJ_weight_coeff_k * ( dpsidx_m * gradx_at_pt +
1085 dpsidy_m * grady_at_pt +
1086 dpsidz_m * gradz_at_pt ) +
1087 bases_val_m * ( advection_term + source_term ) ;
1095 double x[ FunctionCount ] ;
1096 double y[ FunctionCount ] ;
1097 double z[ FunctionCount ] ;
1098 unsigned node_index[ ElemNodeCount ];
1100 scalar_type local_val[FunctionCount], local_res[FunctionCount*(FunctionCount+1)];
1105 gatherSolution( ielem,
val, node_index, x, y, z, elem_res );
1108 computeElementResidual(
val, x, y, z, elem_res );
1111 this->scatterResidual( ielem, node_index, elem_res );
1116 class CoordinateMap ,
typename ScalarType >
1117 class ElementComputation
1119 CrsMatrix< ScalarType , ExecutionSpace > ,
1121 public ElementComputation< Kokkos::Example::BoxElemFixture< ExecutionSpace , Order , CoordinateMap > ,
1122 CrsMatrix< ScalarType , ExecutionSpace > ,
1126 typedef ElementComputation< Kokkos::Example::BoxElemFixture< ExecutionSpace , Order , CoordinateMap > ,
1127 CrsMatrix< ScalarType , ExecutionSpace > ,
1133 static const unsigned FunctionCount = base_type::FunctionCount;
1134 static const unsigned IntegrationCount = base_type::IntegrationCount;
1135 static const unsigned ElemNodeCount = base_type::ElemNodeCount;
1150 base_type(arg_mesh, arg_solution, arg_elem_graph,
1151 arg_jacobian, arg_residual) {}
1157 const size_t nelem = this->elem_node_ids.dimension_0();
1158 parallel_for( nelem , *
this );
1170 double coeff_k = 3.456;
1171 double coeff_src = 1.234;
1172 double advection[] = { 1.1, 1.2, 1.3 };
1173 double dpsidx[ FunctionCount ] ;
1174 double dpsidy[ FunctionCount ] ;
1175 double dpsidz[ FunctionCount ] ;
1181 for (
unsigned i = 0 ; i < IntegrationCount ; ++i ) {
1183 const double integ_weight = this->elem_data.weights[i];
1184 const double* bases_vals = this->elem_data.values[i];
1186 this->transform_gradients( this->elem_data.gradients[i] ,
1188 dpsidx , dpsidy , dpsidz );
1189 const double detJ_weight = detJ * integ_weight;
1190 const double detJ_weight_coeff_k = detJ_weight * coeff_k;
1192 value_at_pt.val() = 0.0 ;
1193 gradx_at_pt.val() = 0.0 ;
1194 grady_at_pt.val() = 0.0 ;
1195 gradz_at_pt.val() = 0.0 ;
1196 for (
unsigned m = 0 ; m < FunctionCount ; m++ ) {
1197 value_at_pt.val() += dof_values(m) * bases_vals[m] ;
1198 gradx_at_pt.val() += dof_values(m) * dpsidx[m] ;
1199 grady_at_pt.val() += dof_values(m) * dpsidy[m] ;
1200 gradz_at_pt.val() += dof_values(m) * dpsidz[m] ;
1204 coeff_src * value_at_pt * value_at_pt ;
1207 advection[0]*gradx_at_pt +
1208 advection[1]*grady_at_pt +
1209 advection[2]*gradz_at_pt;
1211 for (
unsigned m = 0; m < FunctionCount; ++m) {
1212 const double bases_val_m = bases_vals[m] * detJ_weight ;
1214 detJ_weight_coeff_k * ( dpsidx[m] * gradx_at_pt +
1215 dpsidy[m] * grady_at_pt +
1216 dpsidz[m] * gradz_at_pt ) +
1217 bases_val_m * ( advection_term + source_term ) ;
1219 elem_res(m) += res.val();
1221 for(
unsigned n = 0;
n < FunctionCount;
n++) {
1222 elem_mat(m,
n) += res.fastAccessDx(0) * bases_vals[
n] +
1223 res.fastAccessDx(1) * dpsidx[
n] +
1224 res.fastAccessDx(2) * dpsidy[
n] +
1225 res.fastAccessDx(3) * dpsidz[
n];
1234 double x[ FunctionCount ] ;
1235 double y[ FunctionCount ] ;
1236 double z[ FunctionCount ] ;
1237 unsigned node_index[ ElemNodeCount ];
1239 scalar_type local_val[FunctionCount], local_res[FunctionCount], local_mat[FunctionCount*FunctionCount];
1242 elem_mat_type elem_mat(local_mat, FunctionCount, FunctionCount) ;
1245 this->gatherSolution( ielem,
val, node_index, x, y, z, elem_res, elem_mat );
1248 computeElementResidualJacobian(
val, x, y, z, elem_res, elem_mat );
1251 this->scatterResidual( ielem, node_index, elem_res, elem_mat );
KOKKOS_INLINE_FUNCTION void operator()(const unsigned iwork) const
KOKKOS_INLINE_FUNCTION void operator()(const unsigned ielem) const
ElementComputation(const typename base_type::mesh_type &arg_mesh, 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)
KOKKOS_INLINE_FUNCTION void join(volatile unsigned &update, volatile const unsigned &input) const
static const unsigned TensorDim
KOKKOS_INLINE_FUNCTION void gatherSolution(const unsigned ielem, const elem_vec_type &val, unsigned node_index[], double x[], double y[], double z[], const elem_vec_type &res) const
base_type::scalar_type scalar_type
base_type::execution_space execution_space
const elem_vectors_type elem_residuals
static const unsigned SpatialDim
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
static const unsigned function_count
KOKKOS_INLINE_FUNCTION void operator()(const unsigned ielem) const
ElementComputation< Kokkos::Example::BoxElemFixture< ExecutionSpace, Order, CoordinateMap >, CrsMatrix< ScalarType, ExecutionSpace >, Analytic > base_type
Kokkos::View< const double *[SpaceDim], Device > node_coord_type
Kokkos::View< const unsigned *[ElemNode], Device > elem_node_type
ElementComputationBase< ExecutionSpace, Order, CoordinateMap, ScalarType > base_type
ExecutionSpace execution_space
Kokkos::View< scalar_type[FunctionCount], Kokkos::LayoutRight, execution_space, Kokkos::MemoryUnmanaged > scalar_elem_vec_type
Kokkos::View< scalar_type *, Kokkos::LayoutLeft, execution_space > vector_type
base_type::execution_space execution_space
const vector_type residual
Sacado::Fad::SFad< scalar_type, 4 > fad_scalar_type
static const unsigned IntegrationCount
const elem_node_type elem_node_ids
KOKKOS_INLINE_FUNCTION void init(const TagFillNodeSet &, unsigned &update) const
const unsigned node_count
KOKKOS_INLINE_FUNCTION void join(const TagFillNodeSet &, volatile unsigned &update, volatile const unsigned &input) const
KOKKOS_INLINE_FUNCTION void fill_graph_entries(const unsigned iset) const
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::elem_vec_type fad_elem_vec_type
NodeNodeGraph< elem_node_type, sparse_graph_type, ElemNodeCount >::ElemGraphType elem_graph_type
const vector_type solution
base_type::execution_space execution_space
#define KOKKOS_INLINE_FUNCTION
ElemNodeIdView::execution_space execution_space
KOKKOS_INLINE_FUNCTION void operator()(const unsigned ielem) const
ElementComputationBase< ExecutionSpace, Order, CoordinateMap, ScalarType > base_type
Kokkos::View< scalar_type[FunctionCount], Kokkos::LayoutRight, execution_space, Kokkos::MemoryUnmanaged > elem_vec_type
double sort_graph_entries
double fill_graph_entries
KOKKOS_INLINE_FUNCTION void operator()(const TagFillNodeSet &, unsigned ielem, unsigned &count) const
ElementComputationBase(const ElementComputationBase &rhs)
KOKKOS_INLINE_FUNCTION void operator()(const unsigned irow, unsigned &update, const bool final) const
Kokkos::Example::HexElement_Data< mesh_type::ElemNode > element_data_type
CrsGraphType::row_map_type::non_const_type RowMapType
KOKKOS_INLINE_FUNCTION void gatherSolution(const unsigned ielem, const elem_vec_type &val, unsigned node_index[], double x[], double y[], double z[], const elem_vec_type &res, const elem_mat_type &mat) const
const element_data_type elem_data
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)
CrsMatrix(const StaticCrsGraphType &arg_graph)
mesh_type::elem_node_type elem_node_type
Do not initialize the derivative array.
Kokkos::View< fad_scalar_type *, Kokkos::LayoutRight, execution_space, Kokkos::MemoryUnmanaged > elem_vec_type
Kokkos::StaticCrsGraph< unsigned, Space, void, unsigned > StaticCrsGraphType
static const unsigned ElemNodeCount
const ElemNodeIdView elem_node_id
static const unsigned integration_count
mesh_type::node_coord_type node_coord_type
const elem_graph_type elem_graph
ElementComputation(const ElementComputation &rhs)
base_type::elem_vec_type elem_vec_type
KOKKOS_INLINE_FUNCTION void operator()(const unsigned ielem) const
Kokkos::Example::BoxElemFixture< ExecutionSpace, Order, CoordinateMap > mesh_type
Kokkos::View< scalar_type[FunctionCount][FunctionCount], Kokkos::LayoutRight, execution_space, Kokkos::MemoryUnmanaged > elem_mat_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
KOKKOS_INLINE_FUNCTION void computeElementResidual(const elem_vec_type &dof_values, const double x[], const double y[], const double z[], const elem_vec_type &elem_res) const
KOKKOS_INLINE_FUNCTION void computeElementResidualJacobian(const elem_vec_type &dof_values, const double x[], const double y[], const double z[], const elem_vec_type &elem_res, const elem_mat_type &elem_mat) const
base_type::execution_space execution_space
Sacado::Fad::SFad< scalar_type, FunctionCount > fad_scalar_type
KOKKOS_INLINE_FUNCTION void fill_elem_graph_map(const unsigned ielem) const
View< ValueType *, Space > coeff_type
sparse_matrix_type::StaticCrsGraphType sparse_graph_type
Kokkos::View< unsigned, execution_space > UnsignedValue
ElementComputation(const typename base_type::mesh_type &arg_mesh, 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)
KOKKOS_INLINE_FUNCTION void scatterResidual(const unsigned ielem, const unsigned node_index[], const elem_vec_type &res) const
ElementComputation< Kokkos::Example::BoxElemFixture< ExecutionSpace, Order, CoordinateMap >, CrsMatrix< ScalarType, ExecutionSpace >, FadElement > base_type
ElementComputation(const typename base_type::mesh_type &arg_mesh, 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)
ElementComputation(const ElementComputation &rhs)
base_type::fad_scalar_type fad_scalar_type
ElementComputation(const ElementComputation &rhs)
ElementComputation(const typename base_type::mesh_type &arg_mesh, 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)
KOKKOS_INLINE_FUNCTION void computeElementResidualJacobian(const elem_vec_type &dof_values, const double x[], const double y[], const double z[], const elem_vec_type &elem_res, const elem_mat_type &elem_mat) const
base_type::scalar_type scalar_type
base_type::scalar_type scalar_type
double fill_element_graph
Kokkos::View< scalar_type *[FunctionCount][FunctionCount], execution_space > elem_matrices_type
base_type::elem_mat_type elem_mat_type
KOKKOS_INLINE_FUNCTION void gatherSolution(const unsigned ielem, const scalar_elem_vec_type &val, unsigned node_index[], double x[], double y[], double z[], const fad_elem_vec_type &res) const
Kokkos::View< unsigned *[ElemNode][ElemNode], execution_space > ElemGraphType
static const unsigned FunctionCount
Generate a distributed unstructured finite element mesh from a partitioned NX*NY*NZ box of elements...
KOKKOS_INLINE_FUNCTION void scatterResidual(const unsigned ielem, const unsigned node_index[], const elem_vec_type &res, const elem_mat_type &mat) const
KOKKOS_INLINE_FUNCTION void computeElementResidual(const scalar_elem_vec_type dof_values, const double x[], const double y[], const double z[], const fad_elem_vec_type &elem_res) const
KOKKOS_INLINE_FUNCTION void init(unsigned &update) const
base_type::scalar_type scalar_type
Kokkos::UnorderedMap< key_type, void, execution_space > SetType
pair< unsigned, unsigned > key_type
static const unsigned spatial_dimension