FEAT 3
Finite Element Analysis Toolbox
Loading...
Searching...
No Matches
burgers_assembler.cpp
1// FEAT3: Finite Element Analysis Toolbox, Version 3
2// Copyright (C) 2010 by Stefan Turek & the FEAT group
3// FEAT3 is released under the GNU General Public License version 3,
4// see the file 'copyright.txt' in the top level directory for details.
5
6#include <kernel/voxel_assembly/burgers_assembler.hpp>
7#include <kernel/lafem/matrix_gather_scatter_helper.hpp>
8#include <kernel/lafem/vector_gather_scatter_helper.hpp>
9#include <kernel/util/math.hpp>
10#include <kernel/lafem/vector_mirror.hpp>
11#include <kernel/global/vector.hpp>
12
13namespace FEAT
14{
15 namespace VoxelAssembly
16 {
17 namespace Kernel
18 {
19
20 template<typename Space_, typename DT_, typename IT_, FEAT::Intern::MatrixGatherScatterPolicy pol_ = FEAT::Intern::MatrixGatherScatterPolicy::useLocalOps>
21 void full_burgers_assembler_matrix1_bcsr_host(DT_* matrix_data, const DT_* conv_data,
22 const IT_* matrix_row_ptr, const IT_* matrix_col_idx, Index matrix_num_rows, Index matrix_num_cols,
24 const DT_* cub_wg, int num_cubs, DT_ alpha,
25 const IT_* cell_to_dof, [[maybe_unused]] Index cell_num,
26 const Tiny::Vector<DT_, Space_::world_dim>* nodes, [[maybe_unused]] Index node_size,
27 const int* coloring_map, Index coloring_size, const IT_* cell_to_dof_sorter,
28 const VoxelAssembly::AssemblyBurgersData<DT_>& burgers_params)
29 {
30 //define types
31 typedef Space_ SpaceType;
32 typedef DT_ DataType;
33 typedef IT_ IndexType;
34
35 const DataType& beta{burgers_params.beta};
36 const DataType& sd_delta{burgers_params.sd_delta};
37 const DataType& sd_v_norm{burgers_params.sd_v_norm};
38
39 const DataType tol_eps = Math::sqrt(Math::eps<DataType>());
40 // const DataType tol_eps = CudaMath::cuda_sqrt(DBL_EPSILON);
41 // const DataType tol_eps = DBL_EPSILON;
42 const bool need_streamline = (Math::abs(sd_delta) > DataType(0)) && (sd_v_norm > tol_eps);
43 const bool need_convection = Math::abs(beta) > DataType(0);
44
46
47 constexpr int dim = SpaceHelp::dim;
48
49 //define local sizes
50 constexpr int num_loc_dofs = SpaceType::DofMappingType::dof_count;
51 //get number of nodes per element
52 constexpr int num_loc_verts = SpaceType::MeshType::template IndexSet<dim, 0>::Type::num_indices;
53
54 //our local datatypes
55 typedef Tiny::Vector<DataType, dim> VecValueType;
56 typedef Tiny::Matrix<DataType, dim, dim> MatValueType;
57 typedef Tiny::Vector<VecValueType, num_loc_dofs> LocalVectorType;
59
60 FEAT_PRAGMA_OMP(parallel for)
61 for(Index idx = 0; idx < coloring_size; ++idx)
62 {
63 // define local coefficients
65 // and local matrix
66 LocalMatrixType loc_mat;
67 LocalVectorType local_conv_dofs(DataType(0));
68 // now do work for this cell
69 IndexType cell = IndexType(coloring_map[idx]);
70 // std::cout << "Starting with cell " << cell << "\n";
71 const IndexSetWrapper<IndexType> local_dofs_w{cell_to_dof, IndexType(num_loc_dofs)};
72 const IndexType* local_dofs = cell_to_dof + cell*num_loc_dofs;
73 const IndexType* local_dof_sorter = cell_to_dof_sorter + cell*num_loc_dofs;
74
75 SpaceHelp::set_coefficients(local_coeffs, local_dofs_w, nodes, cell);
76 //if we need to, gather local convection vector
77 if(need_convection || need_streamline) //need stream diff or convection?
78 {
79 LAFEM::template VectorGatherScatterHelper<SpaceType, DataType, IndexType>::gather_vector_dense(local_conv_dofs,
80 (const VecValueType*)conv_data, IndexType(matrix_num_rows), local_dofs,DataType(1));
81 }
82
83 VoxelAssembly::Kernel::burgers_mat_assembly_kernel<SpaceHelp>(loc_mat, local_conv_dofs, local_coeffs, cub_pt, cub_wg, num_cubs, burgers_params,
84 need_streamline, need_convection, tol_eps);
85
86 //scatter
87 LAFEM::template MatrixGatherScatterHelper<SpaceType, DataType, IndexType, pol_>::scatter_matrix_csr(loc_mat, (MatValueType*)matrix_data, local_dofs, local_dofs, IndexType(matrix_num_rows), IndexType(matrix_num_cols), matrix_row_ptr, matrix_col_idx, alpha, local_dof_sorter);
88
89 }
90 }
91
92 template<typename Space_, typename DT_, typename IT_>
93 void full_burgers_assembler_vector_bd_host(DT_* vector_data,
94 const DT_* conv_data, const DT_* primal_data, Index vec_size,
96 const DT_* cub_wg, int num_cubs, DT_ alpha,
97 const IT_* cell_to_dof, [[maybe_unused]] Index cell_num,
98 const Tiny::Vector<DT_, Space_::world_dim>* nodes, [[maybe_unused]] Index node_size,
99 const int* coloring_map, Index coloring_size,
100 const VoxelAssembly::AssemblyBurgersData<DT_>& burgers_params)
101 {
102 //define types
103 typedef Space_ SpaceType;
104 typedef DT_ DataType;
105 typedef IT_ IndexType;
106
107 const DataType& beta{burgers_params.beta};
108
109 const bool need_convection = Math::abs(beta) > DataType(0);
110
112
113 constexpr int dim = SpaceHelp::dim;
114
115 //define local sizes
116 constexpr int num_loc_dofs = SpaceType::DofMappingType::dof_count;
117 //get number of nodes per element
118 constexpr int num_loc_verts = SpaceType::MeshType::template IndexSet<dim, 0>::Type::num_indices;
119
120 //our local datatypes
121 typedef Tiny::Vector<DataType, dim> VecValueType;
122 // typedef Tiny::Matrix<DataType, dim, dim> MatValueType;
123 typedef Tiny::Vector<VecValueType, num_loc_dofs> LocalVectorType;
124 // typedef Tiny::Matrix<MatValueType, num_loc_dofs, num_loc_dofs> LocalMatrixType;
125
126 FEAT_PRAGMA_OMP(parallel for)
127 for(Index idx = 0; idx < coloring_size; ++idx)
128 {
129 //define local array
130 LocalVectorType loc_vec(DataType(0));
131 LocalVectorType local_conv_dofs(DataType(0));
132 LocalVectorType local_prim_dofs(DataType(0));
133 // typename NewSpaceHelp::ImagePointType img_point;
135
136 //now do work for this cell
137 IndexType cell = IndexType(coloring_map[idx]);
138 const IndexSetWrapper<IndexType> local_dofs_w{cell_to_dof, IndexType(num_loc_dofs)};
139 const IndexType* local_dofs = cell_to_dof + cell*num_loc_dofs;
140 SpaceHelp::set_coefficients(local_coeffs, local_dofs_w, nodes, cell);
141 LAFEM::template VectorGatherScatterHelper<SpaceType, DataType, IndexType>::gather_vector_dense(local_prim_dofs,
142 (const VecValueType*)primal_data, IndexType(vec_size), local_dofs, DataType(1));
143
144 //if we need to, gather local convection vector
145 if(need_convection) //need stream diff or convection?
146 {
147 LAFEM::template VectorGatherScatterHelper<SpaceType, DataType, IndexType>::gather_vector_dense(local_conv_dofs,
148 (const VecValueType*)conv_data, IndexType(vec_size), local_dofs, DataType(1));
149 }
150
151 VoxelAssembly::Kernel::burgers_defect_assembly_kernel<SpaceHelp>(loc_vec, local_prim_dofs, local_conv_dofs, local_coeffs, cub_pt, cub_wg, num_cubs,
152 burgers_params, need_convection);
153 //scatter
154 LAFEM::template VectorGatherScatterHelper<SpaceType, DataType, IndexType>::scatter_vector_dense(loc_vec,
155 (VecValueType*)vector_data, IndexType(vec_size), local_dofs, alpha);
156
157 }
158 }
159
173 template<typename DT_, int dim_>
174 void set_sd_v_norm_host(const Tiny::Vector<DT_, dim_>* convect, DT_* result, Index vec_size)
175 {
176 DT_ max_val(DT_(0));
177 // simply reduce over max_val
178 FEAT_PRAGMA_OMP(parallel for reduction(max : max_val))
179 for(int i = 0; i < int(vec_size); ++i)
180 {
181 //synchronize all threads in block, to guarentee that all values are written
182 max_val = Math::max(convect[i].norm_euclid(), max_val);
183 }
184 //and overwrite
185 *result = max_val;
186
187 }
188
189 }
190
191 namespace Arch
192 {
193 template<typename Space_, typename DT_, typename IT_>
194 void assemble_burgers_csr_host([[maybe_unused]]const Space_& space,
195 const CSRMatrixData<DT_, IT_>& matrix_data,
196 const DT_* conv_data,
197 const AssemblyCubatureData<DT_>& cubature,
198 const AssemblyMappingData<DT_, IT_>& dof_mapping,
199 const std::vector<int*>& coloring_maps,
200 const std::vector<Index>& coloring_map_sizes,
201 DT_ alpha, const AssemblyBurgersData<DT_>& burgers_params)
202 {
203 for(Index col = 0; col < Index(coloring_maps.size()); ++col)
204 {
205 VoxelAssembly::Kernel::template full_burgers_assembler_matrix1_bcsr_host<Space_, DT_, IT_, FEAT::Intern::MatrixGatherScatterPolicy::useLocalSortHelper>(
206 matrix_data.data, conv_data, matrix_data.row_ptr, matrix_data.col_idx, matrix_data.num_rows, matrix_data.num_cols,
207 (const typename Tiny::Vector<DT_, Space_::world_dim>*) cubature.cub_pt,
208 cubature.cub_wg, cubature.num_cubs, alpha,
209 dof_mapping.cell_to_dof, dof_mapping.cell_num,
210 (const typename Tiny::Vector<DT_, Space_::world_dim>*) dof_mapping.nodes, dof_mapping.node_size,
211 (const int*) coloring_maps[col], coloring_map_sizes[col], dof_mapping.cell_to_dof_sorter,
212 burgers_params
213 );
214 }
215
216 }
217
218 template<typename Space_, typename DT_, typename IT_>
219 void assemble_burgers_defect_host([[maybe_unused]] const Space_& space,
220 DT_* vector_data,
221 const DT_* conv_data,
222 const DT_* primal_data,
223 const AssemblyCubatureData<DT_>& cubature,
224 const AssemblyMappingData<DT_, IT_>& dof_mapping,
225 const std::vector<int*>& coloring_maps,
226 const std::vector<Index>& coloring_map_sizes,
227 DT_ alpha, const AssemblyBurgersData<DT_>& burgers_params)
228 {
229 for(Index col = 0; col < Index(coloring_maps.size()); ++col)
230 {
231 VoxelAssembly::Kernel::template full_burgers_assembler_vector_bd_host<Space_, DT_, IT_>(vector_data, conv_data, primal_data,
232 space.get_num_dofs(), (const typename Tiny::Vector<DT_, Space_::world_dim>*) cubature.cub_pt,
233 cubature.cub_wg, cubature.num_cubs, alpha,
234 dof_mapping.cell_to_dof, dof_mapping.cell_num,
235 (const typename Tiny::Vector<DT_, Space_::world_dim>*) dof_mapping.nodes, dof_mapping.node_size,
236 (const int*) coloring_maps[col], coloring_map_sizes[col], burgers_params
237 );
238 }
239 }
240
241 template<typename DT_, typename IT_, int dim_>
243 {
244 //extract pointer
245 const Tiny::Vector<DT_, dim_>* vec_data = (const Tiny::Vector<DT_, dim_>*)convect.template elements<LAFEM::Perspective::native>();
246
247 DT_ glob_res = DT_(0);
248
249 VoxelAssembly::Kernel::template set_sd_v_norm_host<DT_, dim_>(vec_data, &glob_res, convect.template size<LAFEM::Perspective::native>());
250
251 return glob_res;
252 }
253
254 template<typename DT_, typename IT_, int dim_>
256 {
257 auto local_norm = get_sd_v_norm_host(convect.local());
258 const auto* gate = convect.get_gate();
259 if(gate != nullptr)
260 {
261 local_norm = gate->max(local_norm);
262 }
263 return local_norm;
264 }
265 }
266 }
267}
268
269//instantiate the templates
270using namespace FEAT;
271using namespace FEAT::VoxelAssembly;
272
273
274/*******************************************************2D implementations***************************************************/
276 const std::vector<int*>&, const std::vector<Index>&, double, const AssemblyBurgersData<double>&);
278 const std::vector<int*>&, const std::vector<Index>&, float, const AssemblyBurgersData<float>&);
280 const std::vector<int*>&, const std::vector<Index>&, double, const AssemblyBurgersData<double>&);
282 const std::vector<int*>&, const std::vector<Index>&, float, const AssemblyBurgersData<float>&);
283#ifdef FEAT_HAVE_HALFMATH
285 const std::vector<int*>&, const std::vector<Index>&, Half, const AssemblyBurgersData<Half>&);
287 const std::vector<int*>&, const std::vector<Index>&, Half, const AssemblyBurgersData<Half>&);
288#endif
289
290template void Arch::assemble_burgers_defect_host(const Q2StandardQuad&, double*, const double*, const double*, const AssemblyCubatureData<double>&, const AssemblyMappingData<double, std::uint32_t>&,
291 const std::vector<int*>&, const std::vector<Index>&, double, const AssemblyBurgersData<double>&);
292template void Arch::assemble_burgers_defect_host(const Q2StandardQuad&, float*, const float*, const float*, const AssemblyCubatureData<float>&, const AssemblyMappingData<float, std::uint32_t>&,
293 const std::vector<int*>&, const std::vector<Index>&, float, const AssemblyBurgersData<float>&);
294template void Arch::assemble_burgers_defect_host(const Q2StandardQuad&, double*, const double*, const double*, const AssemblyCubatureData<double>&, const AssemblyMappingData<double, std::uint64_t>&,
295 const std::vector<int*>&, const std::vector<Index>&, double, const AssemblyBurgersData<double>&);
296template void Arch::assemble_burgers_defect_host(const Q2StandardQuad&, float*, const float*, const float*, const AssemblyCubatureData<float>&, const AssemblyMappingData<float, std::uint64_t>&,
297 const std::vector<int*>&, const std::vector<Index>&, float, const AssemblyBurgersData<float>&);
298#ifdef FEAT_HAVE_HALFMATH
300 const std::vector<int*>&, const std::vector<Index>&, Half, const AssemblyBurgersData<Half>&);
302 const std::vector<int*>&, const std::vector<Index>&, Half, const AssemblyBurgersData<Half>&);
303#endif
304
309#ifdef FEAT_HAVE_HALFMATH
312#endif
313
318#ifdef FEAT_HAVE_HALFMATH
321#endif
322
323/*********************************************************3D implementations**************************************************************************************/
325 const std::vector<int*>&, const std::vector<Index>&, double, const AssemblyBurgersData<double>&);
327 const std::vector<int*>&, const std::vector<Index>&, float, const AssemblyBurgersData<float>&);
329 const std::vector<int*>&, const std::vector<Index>&, double, const AssemblyBurgersData<double>&);
331 const std::vector<int*>&, const std::vector<Index>&, float, const AssemblyBurgersData<float>&);
332#ifdef FEAT_HAVE_HALFMATH
334 const std::vector<int*>&, const std::vector<Index>&, Half, const AssemblyBurgersData<Half>&);
336 const std::vector<int*>&, const std::vector<Index>&, Half, const AssemblyBurgersData<Half>&);
337#endif
338
339template void Arch::assemble_burgers_defect_host(const Q2StandardHexa&, double*, const double*, const double*, const AssemblyCubatureData<double>&, const AssemblyMappingData<double, std::uint32_t>&,
340 const std::vector<int*>&, const std::vector<Index>&, double, const AssemblyBurgersData<double>&);
341template void Arch::assemble_burgers_defect_host(const Q2StandardHexa&, float*, const float*, const float*, const AssemblyCubatureData<float>&, const AssemblyMappingData<float, std::uint32_t>&,
342 const std::vector<int*>&, const std::vector<Index>&, float, const AssemblyBurgersData<float>&);
343template void Arch::assemble_burgers_defect_host(const Q2StandardHexa&, double*, const double*, const double*, const AssemblyCubatureData<double>&, const AssemblyMappingData<double, std::uint64_t>&,
344 const std::vector<int*>&, const std::vector<Index>&, double, const AssemblyBurgersData<double>&);
345template void Arch::assemble_burgers_defect_host(const Q2StandardHexa&, float*, const float*, const float*, const AssemblyCubatureData<float>&, const AssemblyMappingData<float, std::uint64_t>&,
346 const std::vector<int*>&, const std::vector<Index>&, float, const AssemblyBurgersData<float>&);
347#ifdef FEAT_HAVE_HALFMATH
349 const std::vector<int*>&, const std::vector<Index>&, Half, const AssemblyBurgersData<Half>&);
351 const std::vector<int*>&, const std::vector<Index>&, Half, const AssemblyBurgersData<Half>&);
352#endif
353
358#ifdef FEAT_HAVE_HALFMATH
361#endif
362
367#ifdef FEAT_HAVE_HALFMATH
370#endif
Global vector wrapper class template.
Definition: vector.hpp:68
Blocked Dense data vector class template.
Handles vector prolongation, restriction and serialization.
Standard Lagrange-2 Finite-Element space class template.
Definition: element.hpp:39
Tiny Matrix class template.
Tiny Vector class template.
T_ sqrt(T_ x)
Returns the square-root of a value.
Definition: math.hpp:300
T_ abs(T_ x)
Returns the absolute value.
Definition: math.hpp:275
T_ max(T_ a, T_ b)
Returns the maximum of two values.
Definition: math.hpp:137
void assemble_burgers_defect_host(const Space_ &space, DT_ *vector_data, const DT_ *conv_data, const DT_ *primal_data, const AssemblyCubatureData< DT_ > &cubature, const AssemblyMappingData< DT_, IT_ > &dof_mapping, const std::vector< int * > &coloring_maps, const std::vector< Index > &coloring_map_sizes, DT_ alpha, const AssemblyBurgersData< DT_ > &burgers_params)
Host kernel wrapper for the defect burgers assembler.
void assemble_burgers_csr_host(const Space_ &space, const CSRMatrixData< DT_, IT_ > &matrix_data, const DT_ *conv_data, const AssemblyCubatureData< DT_ > &cubature, const AssemblyMappingData< DT_, IT_ > &dof_mapping, const std::vector< int * > &coloring_maps, const std::vector< Index > &coloring_map_sizes, DT_ alpha, const AssemblyBurgersData< DT_ > &burgers_params)
Host kernel wrapper for the full matrix burgers assembler.
DT_ get_sd_v_norm_host(const LAFEM::DenseVectorBlocked< DT_, IT_, dim_ > &convect)
Host kernel wrapper for the local sd_v_norm calculation.
void set_sd_v_norm_host(const Tiny::Vector< DT_, dim_ > *convect, DT_ *result, Index vec_size)
Reduces the max local vector norm of a convection vector.
Namespace for different voxel based assembly methods.
FEAT namespace.
Definition: adjactor.hpp:12
__half Half
Half data type.
Definition: half.hpp:25
std::uint64_t Index
Index data type.
const void * cub_pt
The cubature point data array.
A data field for all necessary values that define the dof mapping for assembly.
const IT_ * cell_to_dof
The cell to dof, where cell_to_dof[i],..., cell_to_dof[i+cell_dofs-1] are the dofs of one cell.
const IT_ * cell_to_dof_sorter
Array of sortingindices of cell_to_dof.
const void * nodes
An array of the nodes fitting to the cell_to_dof mapping.