Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://www.mooseframework.org
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "KokkosAssembly.h"
11 :
12 : #include "MooseMesh.h"
13 : #include "FEProblemBase.h"
14 : #include "NonlinearSystemBase.h"
15 : #include "LinearSystem.h"
16 : #include "AuxiliarySystem.h"
17 : #include "Assembly.h"
18 : #include "BoundaryRestrictable.h"
19 :
20 : #include "libmesh/fe_interface.h"
21 : #include "libmesh/reference_elem.h"
22 :
23 : namespace Moose::Kokkos
24 : {
25 :
26 47924 : Assembly::Assembly(FEProblemBase & problem)
27 45530 : : MeshHolder(*problem.mesh().getKokkosMesh()),
28 45530 : _problem(problem),
29 45530 : _mesh(problem.mesh()),
30 136590 : _dimension(_mesh.dimension())
31 : {
32 47924 : }
33 :
34 : void
35 2622 : Assembly::init()
36 : {
37 : // Cache mesh information
38 :
39 2622 : const auto num_subdomains = kokkosMesh().getNumSubdomains();
40 :
41 2622 : _coord_type.create(num_subdomains);
42 :
43 5974 : for (auto subdomain : _mesh.meshSubdomains())
44 3352 : _coord_type[kokkosMesh().getContiguousSubdomainID(subdomain)] = _mesh.getCoordSystem(subdomain);
45 :
46 2622 : _coord_type.copyToDevice();
47 :
48 2622 : if (_mesh.usingGeneralAxisymmetricCoordAxes())
49 : {
50 0 : _rz_axis.create(num_subdomains);
51 :
52 0 : for (auto subdomain : _mesh.meshSubdomains())
53 0 : _rz_axis[kokkosMesh().getContiguousSubdomainID(subdomain)] =
54 0 : _mesh.getGeneralAxisymmetricCoordAxis(subdomain);
55 :
56 0 : _rz_axis.copyToDevice();
57 : }
58 : else
59 2622 : _rz_radial_coord = _mesh.getAxisymmetricRadialCoord();
60 :
61 : // Initialize quadrature and shape data
62 :
63 2622 : initQuadrature();
64 2622 : initShape();
65 2622 : cachePhysicalMap();
66 2622 : }
67 :
68 : void
69 2622 : Assembly::initQuadrature()
70 : {
71 2622 : const auto num_subdomains = kokkosMesh().getNumSubdomains();
72 2622 : const auto num_elem_types = kokkosMesh().getNumLocalElementTypes();
73 :
74 2622 : _q_points.create(num_subdomains, num_elem_types);
75 2622 : _q_points_face.create(num_subdomains, num_elem_types);
76 2622 : _weights.create(num_subdomains, num_elem_types);
77 2622 : _weights_face.create(num_subdomains, num_elem_types);
78 :
79 : // Find boundaries where material properties should be computed
80 :
81 1208 : auto & boundary_objects =
82 1414 : _problem.getKokkosMaterialPropertyStorageConsumers(Moose::BOUNDARY_MATERIAL_DATA);
83 :
84 3476 : for (auto object : boundary_objects)
85 : {
86 854 : auto boundary_restriction = dynamic_cast<const BoundaryRestrictable *>(object);
87 :
88 854 : if (boundary_restriction)
89 1810 : for (auto boundary : boundary_restriction->boundaryIDs())
90 956 : _material_boundaries.insert(boundary);
91 : else
92 0 : mooseError("Kokkos assembly error: ", object->name(), " is not boundary-restricted.");
93 : }
94 :
95 : // Cache quadrature data
96 :
97 2622 : std::map<SubdomainID, std::map<ElemType, unsigned int>> n_qps;
98 2622 : std::map<SubdomainID, std::map<ElemType, std::vector<unsigned int>>> n_qps_face;
99 :
100 2622 : _max_qps_per_elem = 0;
101 :
102 5974 : for (auto subdomain : _mesh.meshSubdomains())
103 : {
104 3352 : auto sid = kokkosMesh().getContiguousSubdomainID(subdomain);
105 :
106 3352 : auto & assembly = _problem.assembly(0, 0);
107 3352 : auto qrule = assembly.writeableQRule(_dimension, subdomain, {});
108 3352 : auto qrule_face = assembly.writeableQRuleFace(_dimension, subdomain, {});
109 :
110 6696 : for (auto & [elem_type, elem_type_id] : kokkosMesh().getElementTypeMap())
111 : {
112 3344 : auto elem = &libMesh::ReferenceElem::get(elem_type);
113 :
114 3344 : _q_points_face(sid, elem_type_id).create(elem->n_sides());
115 3344 : _weights_face(sid, elem_type_id).create(elem->n_sides());
116 :
117 : // Cache volume quadrature of each reference element
118 :
119 3344 : qrule->init(*elem, /* p-level */ 0);
120 3344 : n_qps[subdomain][elem_type] = qrule->n_points();
121 :
122 3344 : _q_points(sid, elem_type_id).create(qrule->n_points());
123 3344 : _weights(sid, elem_type_id).create(qrule->n_points());
124 :
125 15778 : for (const auto qp : make_range(qrule->n_points()))
126 : {
127 12434 : _q_points(sid, elem_type_id)[qp] = qrule->qp(qp);
128 12434 : _weights(sid, elem_type_id)[qp] = qrule->w(qp);
129 : }
130 :
131 : // Cache face quadrature of each reference element
132 :
133 16088 : for (const auto side : elem->side_index_range())
134 : {
135 12744 : qrule_face->init(*elem->side_ptr(side), /* p-level */ 0);
136 12744 : n_qps_face[subdomain][elem_type].push_back(qrule_face->n_points());
137 :
138 12744 : _q_points_face(sid, elem_type_id)[side].create(qrule_face->n_points());
139 12744 : _weights_face(sid, elem_type_id)[side].create(qrule_face->n_points());
140 :
141 37434 : for (const auto qp : make_range(qrule_face->n_points()))
142 : {
143 24690 : _q_points_face(sid, elem_type_id)[side][qp] = qrule_face->qp(qp);
144 24690 : _weights_face(sid, elem_type_id)[side][qp] = qrule_face->w(qp);
145 : }
146 : }
147 :
148 3344 : _max_qps_per_elem = std::max(_max_qps_per_elem, n_qps[subdomain][elem_type]);
149 : }
150 : }
151 :
152 2622 : const auto num_elems = _mesh.nActiveLocalElem();
153 :
154 2622 : _n_subdomain_qps.create(num_subdomains);
155 2622 : _n_subdomain_qps_face.create(num_subdomains);
156 2622 : _n_subdomain_qps = 0;
157 2622 : _n_subdomain_qps_face = 0;
158 :
159 2622 : _n_qps.create(num_elems);
160 2622 : _n_qps_face.create(_mesh.getMaxSidesPerElem(), num_elems);
161 2622 : _n_qps_face = 0;
162 :
163 2622 : _qp_offset.create(num_elems);
164 2622 : _qp_offset_face.create(_mesh.getMaxSidesPerElem(), num_elems);
165 2622 : _qp_offset_face = libMesh::DofObject::invalid_id;
166 :
167 2622 : _elem_face_property_idx.create(_mesh.getMaxSidesPerElem(), num_elems);
168 2622 : _elem_face_property_idx = libMesh::DofObject::invalid_id;
169 :
170 2622 : _n_elem_face_properties.create(num_subdomains);
171 2622 : _n_elem_face_properties = 0;
172 :
173 681222 : for (auto elem : *_mesh.getActiveLocalElementRange())
174 : {
175 678600 : auto eid = kokkosMesh().getContiguousElementID(elem);
176 678600 : auto sid = kokkosMesh().getContiguousSubdomainID(elem->subdomain_id());
177 :
178 678600 : _n_qps[eid] = n_qps[elem->subdomain_id()][elem->type()];
179 678600 : _qp_offset[eid] = _n_subdomain_qps[sid];
180 678600 : _n_subdomain_qps[sid] += _n_qps[eid];
181 :
182 3383492 : for (const auto side : elem->side_index_range())
183 2704892 : _n_qps_face(side, eid) = n_qps_face[elem->subdomain_id()][elem->type()][side];
184 : }
185 :
186 3476 : for (auto boundary : _material_boundaries)
187 6769 : for (auto elem_id : _mesh.getBoundaryActiveSemiLocalElemIds(boundary))
188 : {
189 5915 : auto elem = _mesh.elemPtr(elem_id);
190 :
191 5915 : if (elem->processor_id() == _problem.processor_id())
192 : {
193 5037 : auto sid = kokkosMesh().getContiguousSubdomainID(elem->subdomain_id());
194 5037 : auto eid = kokkosMesh().getContiguousElementID(elem);
195 5037 : auto side = _mesh.sideWithBoundaryID(elem, boundary);
196 :
197 5037 : _qp_offset_face(side, eid) = _n_subdomain_qps_face[sid];
198 5037 : _n_subdomain_qps_face[sid] += _n_qps_face(side, eid);
199 :
200 5037 : _elem_face_property_idx(side, eid) = _n_elem_face_properties[sid];
201 5037 : ++_n_elem_face_properties[sid];
202 : }
203 854 : }
204 :
205 2622 : _q_points.copyToDeviceNested();
206 2622 : _q_points_face.copyToDeviceNested();
207 2622 : _weights.copyToDeviceNested();
208 2622 : _weights_face.copyToDeviceNested();
209 :
210 2622 : _n_qps.copyToDevice();
211 2622 : _n_qps_face.copyToDevice();
212 2622 : _n_subdomain_qps.copyToDevice();
213 2622 : _n_subdomain_qps_face.copyToDevice();
214 2622 : _qp_offset.copyToDevice();
215 2622 : _qp_offset_face.copyToDevice();
216 :
217 2622 : _elem_face_property_idx.copyToDevice();
218 2622 : _n_elem_face_properties.copyToDevice();
219 2622 : }
220 :
221 : void
222 2622 : Assembly::initShape()
223 : {
224 : // Generate the list of unique FE types
225 :
226 2622 : std::set<FEType> fe_types;
227 :
228 5273 : auto getFETypes = [&](::System & system)
229 : {
230 9774 : for (const auto var : make_range(system.n_vars()))
231 : {
232 4501 : const auto fe_type = system.variable_type(var);
233 4501 : const auto field_type = FEInterface::field_type(fe_type);
234 :
235 4501 : if (field_type == libMesh::TYPE_VECTOR)
236 : {
237 172 : if (fe_type.family != libMesh::LAGRANGE_VEC)
238 0 : mooseError("Kokkos currently only supports LAGRANGE_VEC for vector FE families.");
239 : }
240 : else
241 : {
242 4329 : if (fe_type.family != libMesh::LAGRANGE && fe_type.family != libMesh::L2_LAGRANGE &&
243 338 : fe_type.family != libMesh::MONOMIAL && fe_type.family != libMesh::SCALAR)
244 0 : mooseError("Kokkos currently only supports LAGRANGE, L2_LAGRNAGE, MONOMIAL, SCALAR for "
245 : "scalar FE families.");
246 : }
247 :
248 4501 : fe_types.insert(fe_type);
249 : }
250 5273 : };
251 :
252 5078 : for (const auto nl : make_range(_problem.numNonlinearSystems()))
253 2456 : getFETypes(_problem.getNonlinearSystemBase(nl).system());
254 :
255 2817 : for (const auto linear : make_range(_problem.numLinearSystems()))
256 195 : getFETypes(_problem.getLinearSystem(linear).system());
257 :
258 2622 : getFETypes(_problem.getAuxiliarySystem().system());
259 :
260 2622 : _fe_type_map.clear();
261 :
262 5653 : for (auto & fet : fe_types)
263 3031 : _fe_type_map[fet] = _fe_type_map.size();
264 :
265 : // Cache reference shape data
266 :
267 2622 : const auto num_subdomains = kokkosMesh().getNumSubdomains();
268 2622 : const auto num_elem_types = kokkosMesh().getNumLocalElementTypes();
269 :
270 2622 : _phi.create(num_subdomains, num_elem_types, _fe_type_map.size());
271 2622 : _phi_face.create(num_subdomains, num_elem_types, _fe_type_map.size());
272 2622 : _grad_phi.create(num_subdomains, num_elem_types, _fe_type_map.size());
273 2622 : _grad_phi_face.create(num_subdomains, num_elem_types, _fe_type_map.size());
274 2622 : _vector_phi.create(num_subdomains, num_elem_types, _fe_type_map.size());
275 2622 : _vector_phi_face.create(num_subdomains, num_elem_types, _fe_type_map.size());
276 2622 : _vector_grad_phi.create(num_subdomains, num_elem_types, _fe_type_map.size());
277 2622 : _vector_grad_phi_face.create(num_subdomains, num_elem_types, _fe_type_map.size());
278 2622 : _is_vector_fe_type.create(_fe_type_map.size());
279 :
280 2622 : _map_phi.create(num_subdomains, num_elem_types);
281 2622 : _map_phi_face.create(num_subdomains, num_elem_types);
282 2622 : _map_psi_face.create(num_subdomains, num_elem_types);
283 2622 : _map_grad_phi.create(num_subdomains, num_elem_types);
284 2622 : _map_grad_phi_face.create(num_subdomains, num_elem_types);
285 2622 : _map_grad_psi_face.create(num_subdomains, num_elem_types);
286 :
287 2622 : _normal_dx_dxi.create(num_subdomains, num_elem_types);
288 2622 : _normal_dx_deta.create(num_subdomains, num_elem_types);
289 :
290 2622 : _n_dofs.create(num_elem_types, _fe_type_map.size());
291 2622 : _n_dofs = 0;
292 :
293 5653 : for (auto & [fe_type, fe_type_id] : _fe_type_map)
294 3031 : _is_vector_fe_type[fe_type_id] = FEInterface::field_type(fe_type) == libMesh::TYPE_VECTOR;
295 :
296 5974 : for (auto subdomain : _mesh.meshSubdomains())
297 : {
298 3352 : auto sid = kokkosMesh().getContiguousSubdomainID(subdomain);
299 :
300 3352 : auto & assembly = _problem.assembly(0, 0);
301 3352 : auto qrule = assembly.writeableQRule(_dimension, subdomain, {});
302 3352 : auto qrule_face = assembly.writeableQRuleFace(_dimension, subdomain, {});
303 :
304 7317 : for (auto & [fe_type, fe_type_id] : _fe_type_map)
305 : {
306 3965 : const bool is_vector_fe = FEInterface::field_type(fe_type) == libMesh::TYPE_VECTOR;
307 :
308 3965 : std::unique_ptr<FEBase> fe;
309 3965 : std::unique_ptr<FEBase> fe_face;
310 3965 : std::unique_ptr<FEVectorBase> vector_fe;
311 3965 : std::unique_ptr<FEVectorBase> vector_fe_face;
312 :
313 3965 : if (is_vector_fe)
314 : {
315 150 : vector_fe = FEVectorBase::build(_dimension, fe_type);
316 150 : vector_fe_face = FEVectorBase::build(_dimension, fe_type);
317 :
318 150 : vector_fe->attach_quadrature_rule(qrule);
319 150 : vector_fe_face->attach_quadrature_rule(qrule_face);
320 : }
321 : else
322 : {
323 3815 : fe = FEBase::build(_dimension, fe_type);
324 3815 : fe_face = FEBase::build(_dimension, fe_type);
325 :
326 3815 : fe->attach_quadrature_rule(qrule);
327 3815 : fe_face->attach_quadrature_rule(qrule_face);
328 : }
329 :
330 7918 : for (auto & [elem_type, elem_type_id] : kokkosMesh().getElementTypeMap())
331 : {
332 3953 : auto elem = &libMesh::ReferenceElem::get(elem_type);
333 :
334 3953 : if (is_vector_fe)
335 : {
336 146 : auto & phi = vector_fe->get_phi();
337 146 : auto & grad_phi = vector_fe->get_dphi();
338 :
339 146 : vector_fe->reinit(elem);
340 :
341 146 : _n_dofs(elem_type_id, fe_type_id) = phi.size();
342 :
343 146 : _vector_phi(sid, elem_type_id, fe_type_id).create(phi.size(), qrule->n_points());
344 146 : _vector_grad_phi(sid, elem_type_id, fe_type_id).create(phi.size(), qrule->n_points());
345 :
346 1339 : for (const auto i : index_range(phi))
347 7408 : for (const auto qp : make_range(qrule->n_points()))
348 : {
349 6215 : _vector_phi(sid, elem_type_id, fe_type_id)(i, qp) = phi[i][qp];
350 6215 : _vector_grad_phi(sid, elem_type_id, fe_type_id)(i, qp) = grad_phi[i][qp];
351 : }
352 :
353 146 : _vector_phi_face(sid, elem_type_id, fe_type_id).create(elem->n_sides());
354 146 : _vector_grad_phi_face(sid, elem_type_id, fe_type_id).create(elem->n_sides());
355 :
356 672 : for (const auto side : elem->side_index_range())
357 : {
358 526 : auto & phi = vector_fe_face->get_phi();
359 526 : auto & grad_phi = vector_fe_face->get_dphi();
360 :
361 526 : vector_fe_face->reinit(elem, side);
362 :
363 526 : _vector_phi_face(sid, elem_type_id, fe_type_id)(side).create(phi.size(),
364 : qrule_face->n_points());
365 526 : _vector_grad_phi_face(sid, elem_type_id, fe_type_id)(side).create(
366 : phi.size(), qrule_face->n_points());
367 :
368 5124 : for (const auto i : index_range(phi))
369 14844 : for (const auto qp : make_range(qrule_face->n_points()))
370 : {
371 10246 : _vector_phi_face(sid, elem_type_id, fe_type_id)(side)(i, qp) = phi[i][qp];
372 10246 : _vector_grad_phi_face(sid, elem_type_id, fe_type_id)(side)(i, qp) = grad_phi[i][qp];
373 : }
374 : }
375 : }
376 : else
377 : {
378 3807 : auto & phi = fe->get_phi();
379 3807 : auto & grad_phi = fe->get_dphi();
380 :
381 3807 : fe->reinit(elem);
382 :
383 3807 : _n_dofs(elem_type_id, fe_type_id) = phi.size();
384 :
385 3807 : _phi(sid, elem_type_id, fe_type_id).create(phi.size(), qrule->n_points());
386 3807 : _grad_phi(sid, elem_type_id, fe_type_id).create(phi.size(), qrule->n_points());
387 :
388 17447 : for (const auto i : index_range(phi))
389 74100 : for (const auto qp : make_range(qrule->n_points()))
390 : {
391 60460 : _phi(sid, elem_type_id, fe_type_id)(i, qp) = phi[i][qp];
392 60460 : _grad_phi(sid, elem_type_id, fe_type_id)(i, qp) = grad_phi[i][qp];
393 : }
394 :
395 3807 : _phi_face(sid, elem_type_id, fe_type_id).create(elem->n_sides());
396 3807 : _grad_phi_face(sid, elem_type_id, fe_type_id).create(elem->n_sides());
397 :
398 18539 : for (const auto side : elem->side_index_range())
399 : {
400 14732 : auto & phi = fe_face->get_phi();
401 14732 : auto & grad_phi = fe_face->get_dphi();
402 :
403 14732 : fe_face->reinit(elem, side);
404 :
405 14732 : _phi_face(sid, elem_type_id, fe_type_id)(side).create(phi.size(),
406 : qrule_face->n_points());
407 14732 : _grad_phi_face(sid, elem_type_id, fe_type_id)(side).create(phi.size(),
408 : qrule_face->n_points());
409 :
410 69526 : for (const auto i : index_range(phi))
411 173916 : for (const auto qp : make_range(qrule_face->n_points()))
412 : {
413 119122 : _phi_face(sid, elem_type_id, fe_type_id)(side)(i, qp) = phi[i][qp];
414 119122 : _grad_phi_face(sid, elem_type_id, fe_type_id)(side)(i, qp) = grad_phi[i][qp];
415 : }
416 : }
417 : }
418 : }
419 3965 : }
420 :
421 6696 : for (auto & [elem_type, elem_type_id] : kokkosMesh().getElementTypeMap())
422 : {
423 3344 : auto elem = &libMesh::ReferenceElem::get(elem_type);
424 3344 : auto fe_type = FEType(elem->default_order(), LAGRANGE);
425 3344 : auto shape_deriv_ptr = libMesh::FEInterface::shape_deriv_function(fe_type, elem);
426 :
427 3344 : std::unique_ptr<FEBase> fe(FEBase::build(_dimension, fe_type));
428 3344 : std::unique_ptr<FEBase> fe_face(FEBase::build(_dimension, fe_type));
429 :
430 3344 : fe->attach_quadrature_rule(qrule);
431 3344 : fe_face->attach_quadrature_rule(qrule_face);
432 :
433 3344 : auto & phi = fe->get_phi();
434 3344 : auto & grad_phi = fe->get_dphi();
435 :
436 3344 : fe->reinit(elem);
437 :
438 3344 : _map_phi(sid, elem_type_id).create(phi.size(), qrule->n_points());
439 3344 : _map_grad_phi(sid, elem_type_id).create(phi.size(), qrule->n_points());
440 :
441 16922 : for (const auto i : index_range(phi))
442 69731 : for (const auto qp : make_range(qrule->n_points()))
443 : {
444 56153 : _map_phi(sid, elem_type_id)(i, qp) = phi[i][qp];
445 56153 : _map_grad_phi(sid, elem_type_id)(i, qp) = grad_phi[i][qp];
446 : }
447 :
448 3344 : _map_phi_face(sid, elem_type_id).create(elem->n_sides());
449 3344 : _map_grad_phi_face(sid, elem_type_id).create(elem->n_sides());
450 3344 : _map_psi_face(sid, elem_type_id).create(elem->n_sides());
451 3344 : _map_grad_psi_face(sid, elem_type_id).create(elem->n_sides());
452 :
453 3344 : _normal_dx_dxi(sid, elem_type_id).create(elem->n_sides());
454 3344 : _normal_dx_deta(sid, elem_type_id).create(elem->n_sides());
455 :
456 16088 : for (const auto side : elem->side_index_range())
457 : {
458 12744 : auto side_ptr = elem->build_side_ptr(side);
459 :
460 12744 : auto & phi = fe_face->get_phi();
461 12744 : auto & grad_phi = fe_face->get_dphi();
462 :
463 12744 : fe_face->reinit(elem, side);
464 :
465 12744 : _map_phi_face(sid, elem_type_id)(side).create(phi.size(), qrule_face->n_points());
466 12744 : _map_grad_phi_face(sid, elem_type_id)(side).create(phi.size(), qrule_face->n_points());
467 :
468 66754 : for (const auto i : index_range(phi))
469 166504 : for (const auto qp : make_range(qrule_face->n_points()))
470 : {
471 112494 : _map_phi_face(sid, elem_type_id)(side)(i, qp) = phi[i][qp];
472 112494 : _map_grad_phi_face(sid, elem_type_id)(side)(i, qp) = grad_phi[i][qp];
473 : }
474 :
475 12744 : auto & psi = fe_face->get_fe_map().get_psi();
476 12744 : auto & dpsidxi = fe_face->get_fe_map().get_dpsidxi();
477 12744 : auto & dpsideta = fe_face->get_fe_map().get_dpsideta();
478 :
479 12744 : _map_psi_face(sid, elem_type_id)(side).create(psi.size(), qrule_face->n_points());
480 12744 : _map_grad_psi_face(sid, elem_type_id)(side).create(psi.size(), qrule_face->n_points());
481 :
482 38958 : for (const auto i : index_range(psi))
483 80296 : for (const auto qp : make_range(qrule_face->n_points()))
484 : {
485 54082 : _map_psi_face(sid, elem_type_id)(side)(i, qp) = psi[i][qp];
486 54082 : if (elem->dim() > 1)
487 53280 : _map_grad_psi_face(sid, elem_type_id)(side)(i, qp)(0) = dpsidxi[i][qp];
488 54082 : if (elem->dim() > 2)
489 8160 : _map_grad_psi_face(sid, elem_type_id)(side)(i, qp)(1) = dpsideta[i][qp];
490 : }
491 :
492 12744 : _normal_dx_dxi(sid, elem_type_id)(side).create(phi.size(), qrule_face->n_points());
493 12744 : _normal_dx_deta(sid, elem_type_id)(side).create(phi.size(), qrule_face->n_points());
494 :
495 37434 : for (const auto qp : make_range(qrule_face->n_points()))
496 : {
497 24690 : Point reference_point;
498 :
499 24690 : if (_dimension == 1)
500 802 : reference_point = side ? Point(1) : Point(-1);
501 23888 : else if (_dimension == 2)
502 66968 : for (const auto i : index_range(psi))
503 45120 : reference_point.add_scaled(side_ptr->point(i), psi[i][qp]);
504 :
505 137184 : for (const auto i : index_range(phi))
506 : {
507 112494 : if (_dimension < 3)
508 96174 : _normal_dx_dxi(sid, elem_type_id)(side)(i, qp) =
509 51922 : shape_deriv_ptr(fe_type, elem, i, 0, reference_point, false);
510 112494 : if (_dimension == 2)
511 94512 : _normal_dx_deta(sid, elem_type_id)(side)(i, qp) =
512 51056 : shape_deriv_ptr(fe_type, elem, i, 1, reference_point, false);
513 : }
514 : }
515 12744 : }
516 3344 : }
517 : }
518 :
519 2622 : _phi.copyToDeviceNested();
520 2622 : _phi_face.copyToDeviceNested();
521 2622 : _grad_phi.copyToDeviceNested();
522 2622 : _grad_phi_face.copyToDeviceNested();
523 2622 : _vector_phi.copyToDeviceNested();
524 2622 : _vector_phi_face.copyToDeviceNested();
525 2622 : _vector_grad_phi.copyToDeviceNested();
526 2622 : _vector_grad_phi_face.copyToDeviceNested();
527 2622 : _is_vector_fe_type.copyToDevice();
528 :
529 2622 : _map_phi.copyToDeviceNested();
530 2622 : _map_phi_face.copyToDeviceNested();
531 2622 : _map_psi_face.copyToDeviceNested();
532 2622 : _map_grad_phi.copyToDeviceNested();
533 2622 : _map_grad_phi_face.copyToDeviceNested();
534 2622 : _map_grad_psi_face.copyToDeviceNested();
535 :
536 2622 : _normal_dx_dxi.copyToDeviceNested();
537 2622 : _normal_dx_deta.copyToDeviceNested();
538 :
539 2622 : _n_dofs.copyToDevice();
540 2622 : }
541 :
542 : void
543 2622 : Assembly::cachePhysicalMap()
544 : {
545 2622 : const auto num_subdomains = kokkosMesh().getNumSubdomains();
546 2622 : const auto num_elems = kokkosMesh().getNumLocalElements();
547 :
548 2622 : _jacobian.create(num_subdomains);
549 2622 : _jxw.create(num_subdomains);
550 2622 : _xyz.create(num_subdomains);
551 :
552 5974 : for (auto subdomain : _mesh.meshSubdomains())
553 : {
554 3352 : auto sid = kokkosMesh().getContiguousSubdomainID(subdomain);
555 :
556 3352 : _jacobian[sid].createDevice(_n_subdomain_qps[sid]);
557 3352 : _jxw[sid].createDevice(_n_subdomain_qps[sid]);
558 3352 : _xyz[sid].createDevice(_n_subdomain_qps[sid]);
559 : }
560 :
561 2622 : _jacobian.copyToDeviceNested();
562 2622 : _jxw.copyToDeviceNested();
563 2622 : _xyz.copyToDeviceNested();
564 :
565 2622 : ::Kokkos::RangePolicy<ExecSpace, ::Kokkos::IndexType<ThreadID>> policy(0, num_elems);
566 2622 : ::Kokkos::parallel_for(policy, *this);
567 2622 : ::Kokkos::fence();
568 2622 : }
569 :
570 : KOKKOS_FUNCTION void
571 403307 : Assembly::operator()(const ThreadID tid) const
572 : {
573 403307 : auto info = kokkosMesh().getElementInfo(tid);
574 403307 : auto offset = getQpOffset(info);
575 :
576 403307 : auto jacobian = &_jacobian[info.subdomain][offset];
577 403307 : auto jxw = &_jxw[info.subdomain][offset];
578 403307 : auto xyz = &_xyz[info.subdomain][offset];
579 :
580 1069721 : for (unsigned int qp = 0; qp < getNumQps(info); ++qp)
581 666414 : computePhysicalMap(info, qp, &jacobian[qp], &jxw[qp], &xyz[qp]);
582 403307 : }
583 :
584 : } // namespace Moose::Kokkos
|
