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 "KokkosFESystem.h"
11 : #include "KokkosNodalBCBase.h"
12 :
13 : #include "MooseMesh.h"
14 : #include "Assembly.h"
15 : #include "NonlinearSystemBase.h"
16 : #include "LinearSystem.h"
17 : #include "FEProblemBase.h"
18 :
19 : #include "libmesh/system.h"
20 : #include "libmesh/fe_interface.h"
21 :
22 : namespace Moose::Kokkos
23 : {
24 :
25 3150 : FESystem::FESystem(SystemBase & system)
26 1714 : : System(system), AssemblyHolder(system.feProblem().kokkosAssembly())
27 : {
28 3150 : setupVariables();
29 3150 : setupDofs();
30 3150 : setupCoupling();
31 :
32 3150 : if (dynamic_cast<NonlinearSystemBase *>(&_system))
33 2213 : setupNodalBCDofs();
34 :
35 3150 : _qp_solutions.create(MAX_TAG);
36 3150 : _qp_solutions_grad.create(MAX_TAG);
37 3150 : _qp_vector_solutions.create(MAX_TAG);
38 3150 : _qp_vector_solutions_grad.create(MAX_TAG);
39 3150 : _qp_vector_solutions_curl.create(MAX_TAG);
40 3150 : }
41 :
42 49 : FESystem::FESystem(System & base, SystemBase & system)
43 26 : : System(base), AssemblyHolder(system.feProblem().kokkosAssembly())
44 : {
45 49 : setupVariables();
46 49 : setupDofs();
47 :
48 49 : if (dynamic_cast<NonlinearSystemBase *>(&_system))
49 2 : setupNodalBCDofs();
50 :
51 49 : _qp_solutions.create(MAX_TAG);
52 49 : _qp_solutions_grad.create(MAX_TAG);
53 49 : _qp_vector_solutions.create(MAX_TAG);
54 49 : _qp_vector_solutions_grad.create(MAX_TAG);
55 49 : _qp_vector_solutions_curl.create(MAX_TAG);
56 49 : }
57 :
58 : void
59 3199 : FESystem::setupVariables()
60 : {
61 3199 : auto & sys = _system.system();
62 :
63 3199 : _var_fe_types.create(_num_vars);
64 3199 : _var_is_vector.create(_num_vars);
65 :
66 7523 : for (unsigned int var = 0; var < _num_vars; ++var)
67 : {
68 4324 : const auto fe_type = sys.variable_type(var);
69 :
70 4324 : _var_fe_types[var] = kokkosAssembly().getFETypeID(fe_type);
71 4324 : _var_is_vector[var] = FEInterface::field_type(fe_type) == libMesh::TYPE_VECTOR;
72 : }
73 :
74 3199 : _var_fe_types.copyToDevice();
75 3199 : _var_is_vector.copyToDevice();
76 3199 : }
77 :
78 : void
79 3150 : FESystem::setupCoupling()
80 : {
81 3150 : if (auto * const nl_system = dynamic_cast<NonlinearSystemBase *>(&_system))
82 : {
83 2213 : _coupling.create(_num_vars);
84 :
85 2213 : std::map<unsigned int, std::vector<unsigned int>> coupling;
86 :
87 2213 : auto & ce = _system.feProblem().couplingEntries(0, nl_system->number());
88 :
89 6099 : for (const auto & [ivar, jvar] : ce)
90 3886 : if (ivar->number() != jvar->number())
91 1116 : coupling[ivar->number()].push_back(jvar->number());
92 :
93 4987 : for (const auto var : make_range(_num_vars))
94 2774 : _coupling[var] = coupling[var];
95 :
96 2213 : _coupling.copyToDevice();
97 2213 : }
98 3150 : }
99 :
100 : void
101 3199 : FESystem::setupDofs()
102 : {
103 3199 : auto & sys = _system.system();
104 :
105 3199 : auto num_nodes = kokkosMesh().getNumLocalNodes();
106 :
107 3199 : _local_node_dof_index.create(_num_vars);
108 :
109 1459 : auto * const solution =
110 1740 : libMesh::cast_ptr<PetscVector<Number> *>(sys.current_local_solution.get());
111 :
112 : #ifdef MOOSE_ENABLE_KOKKOS_GPU
113 : // Kokkos array thinks OpenMP clause is device code when using OpenMP backend
114 1459 : #pragma omp parallel for
115 : #endif
116 4077 : for (unsigned int var = 0; var < _num_vars; ++var)
117 : {
118 2337 : std::vector<dof_id_type> dof_indices;
119 :
120 2337 : _local_node_dof_index[var].create(num_nodes);
121 2337 : _local_node_dof_index[var] = libMesh::DofObject::invalid_id;
122 :
123 914767 : for (const auto node : kokkosMesh().getLocalNodes())
124 912430 : if (node->processor_id() == _comm.rank())
125 : {
126 904183 : const auto id = kokkosMesh().getContiguousNodeID(node);
127 :
128 904183 : _dof_map.dof_indices(node, dof_indices, var);
129 :
130 904183 : if (dof_indices.size())
131 : {
132 873981 : for (unsigned int i = 1; i < dof_indices.size(); ++i)
133 14165 : if (dof_indices[i] != dof_indices[i - 1] + 1)
134 0 : mooseError("Kokkos system error: a variable has multiple DOFs on a node, but the DOF "
135 : "indices are discontiguous. This is not supported.");
136 :
137 859816 : _local_node_dof_index[var][id] = solution->map_global_to_local_index(dof_indices[0]);
138 : }
139 : }
140 2337 : }
141 :
142 3199 : _local_node_dof_index.copyToDeviceNested();
143 3199 : }
144 :
145 : void
146 2215 : FESystem::setupNodalBCDofs()
147 : {
148 2215 : auto & nl_system = static_cast<NonlinearSystemBase &>(_system);
149 :
150 2215 : _nbc_matrix_tag_dof.create(MAX_TAG);
151 :
152 5952 : for (auto bc : nl_system.getKokkosNodalBCWarehouse().getActiveObjects())
153 : {
154 3737 : auto nbc = static_cast<NodalBCBase *>(bc.get());
155 :
156 3737 : auto matrix_tags = nbc->getMatrixTags({});
157 :
158 11398 : for (auto tag : matrix_tags)
159 7661 : getNodalBCDofs(nbc, _nbc_matrix_tag_dof[tag]);
160 3737 : }
161 :
162 2215 : _nbc_matrix_tag_dof.copyToDeviceNested();
163 2215 : }
164 :
165 : void
166 7661 : FESystem::getNodalBCDofs(const NodalBCBase * nbc, Array<bool> & dofs)
167 : {
168 7661 : auto var = nbc->variable().number();
169 7661 : auto num_dofs = _var_is_vector[var] ? kokkosAssembly().getDimension() : 1;
170 :
171 7661 : if (!dofs.isAlloc())
172 : {
173 4012 : dofs.create(_num_local_dofs + _num_ghost_dofs);
174 4012 : dofs = false;
175 : }
176 :
177 58087 : for (auto node : nbc->getContiguousNodes())
178 : {
179 50426 : auto dof = _local_node_dof_index[var][node];
180 :
181 50426 : if (dof != libMesh::DofObject::invalid_id)
182 102204 : for (const auto i : make_range(num_dofs))
183 53858 : dofs[dof + i] = true;
184 7661 : }
185 :
186 : // Let remote processes know about ghost DOFs associated with nodal BCs not to contribute on them
187 :
188 7661 : auto num_procs = _comm.size();
189 :
190 22983 : std::vector<std::vector<char>> send(num_procs), recv(num_procs);
191 :
192 19074 : for (processor_id_type proc = 0; proc < num_procs; proc++)
193 36421 : for (auto dof : _local_comm_list[proc])
194 25008 : send[proc].push_back(dofs[dof]);
195 :
196 19074 : for (processor_id_type proc = 0; proc < num_procs; proc++)
197 11413 : _comm.scatter(send, recv[proc], proc);
198 :
199 19074 : for (processor_id_type proc = 0; proc < num_procs; proc++)
200 36421 : for (dof_id_type i = 0; i < _ghost_comm_list[proc].size(); ++i)
201 25008 : dofs[_ghost_comm_list[proc][i]] = recv[proc][i];
202 7661 : }
203 :
204 : void
205 144746 : FESystem::reinit()
206 : {
207 494167 : for (auto tag : _active_solution_tags)
208 : {
209 349421 : if (!_qp_solutions[tag].isAlloc())
210 4936 : _qp_solutions[tag].create(kokkosMesh().getNumSubdomains(), _num_vars);
211 :
212 349421 : if (!_qp_solutions_grad[tag].isAlloc())
213 4936 : _qp_solutions_grad[tag].create(kokkosMesh().getNumSubdomains(), _num_vars);
214 :
215 349421 : if (!_qp_vector_solutions[tag].isAlloc())
216 4936 : _qp_vector_solutions[tag].create(kokkosMesh().getNumSubdomains(), _num_vars);
217 :
218 349421 : if (!_qp_vector_solutions_grad[tag].isAlloc())
219 4936 : _qp_vector_solutions_grad[tag].create(kokkosMesh().getNumSubdomains(), _num_vars);
220 :
221 349421 : if (!_qp_vector_solutions_curl[tag].isAlloc())
222 4936 : _qp_vector_solutions_curl[tag].create(kokkosMesh().getNumSubdomains(), _num_vars);
223 :
224 730233 : for (auto subdomain : _mesh.meshSubdomains())
225 : {
226 380812 : auto sid = kokkosMesh().getContiguousSubdomainID(subdomain);
227 :
228 797223 : for (auto var : _active_variables)
229 : {
230 416411 : if (!_var_subdomain_active(var, sid))
231 626 : continue;
232 :
233 415785 : if (_var_is_vector[var])
234 : {
235 37930 : if (!_qp_vector_solutions[tag](sid, var).isAlloc())
236 515 : _qp_vector_solutions[tag](sid, var).createDevice(kokkosAssembly().getNumQps(sid));
237 :
238 37930 : if (!_qp_vector_solutions_grad[tag](sid, var).isAlloc())
239 515 : _qp_vector_solutions_grad[tag](sid, var).createDevice(kokkosAssembly().getNumQps(sid));
240 :
241 37930 : if (!_qp_vector_solutions_curl[tag](sid, var).isAlloc())
242 515 : _qp_vector_solutions_curl[tag](sid, var).createDevice(kokkosAssembly().getNumQps(sid));
243 : }
244 : else
245 : {
246 377855 : if (!_qp_solutions[tag](sid, var).isAlloc())
247 6720 : _qp_solutions[tag](sid, var).createDevice(kokkosAssembly().getNumQps(sid));
248 :
249 377855 : if (!_qp_solutions_grad[tag](sid, var).isAlloc())
250 6720 : _qp_solutions_grad[tag](sid, var).createDevice(kokkosAssembly().getNumQps(sid));
251 : }
252 : }
253 : }
254 :
255 349421 : _qp_solutions[tag].copyToDevice();
256 349421 : _qp_solutions_grad[tag].copyToDevice();
257 349421 : _qp_vector_solutions[tag].copyToDevice();
258 349421 : _qp_vector_solutions_grad[tag].copyToDevice();
259 349421 : _qp_vector_solutions_curl[tag].copyToDevice();
260 : }
261 :
262 144746 : _qp_solutions.copyToDevice();
263 144746 : _qp_solutions_grad.copyToDevice();
264 144746 : _qp_vector_solutions.copyToDevice();
265 144746 : _qp_vector_solutions_grad.copyToDevice();
266 144746 : _qp_vector_solutions_curl.copyToDevice();
267 :
268 144746 : dof_id_type num_elems = kokkosMesh().getNumLocalElements();
269 :
270 144746 : _thread.resize(kokkosAssembly().getMaxQpsPerElem(),
271 : num_elems,
272 : _active_variables.size(),
273 : _active_solution_tags.size());
274 :
275 144746 : ::Kokkos::RangePolicy<ExecSpace, ::Kokkos::IndexType<ThreadID>> policy(0, _thread.size());
276 144746 : ::Kokkos::parallel_for(policy, *this);
277 144746 : ::Kokkos::fence();
278 144746 : }
279 :
280 : KOKKOS_FUNCTION void
281 70666403 : FESystem::operator()(const ThreadID tid) const
282 : {
283 70666403 : auto qp = _thread(tid, 0);
284 70666403 : auto elem = _thread(tid, 1);
285 70666403 : auto var = _active_variables(_thread(tid, 2));
286 70666403 : auto tag = _active_solution_tags(_thread(tid, 3));
287 :
288 70666403 : auto info = kokkosMesh().getElementInfo(elem);
289 70666403 : auto sid = info.subdomain;
290 70666403 : auto elem_type = info.type;
291 :
292 70666403 : if (!_var_subdomain_active(var, sid))
293 3872 : return;
294 :
295 70662531 : auto fe_type = _var_fe_types[var];
296 70662531 : auto num_dofs = kokkosAssembly().getNumDofs(elem_type, fe_type);
297 70662531 : auto num_qps = kokkosAssembly().getNumQps(info);
298 70662531 : auto qp_offset = kokkosAssembly().getQpOffset(info);
299 :
300 70662531 : if (qp >= num_qps)
301 0 : return;
302 :
303 70662531 : if (_var_is_vector[var])
304 : {
305 2516965 : auto & phi = kokkosAssembly().getVectorPhi(sid, elem_type, fe_type);
306 2516965 : auto & grad_phi = kokkosAssembly().getVectorGradPhi(sid, elem_type, fe_type);
307 2516965 : auto jacobian = kokkosAssembly().getJacobian(info, qp).transpose();
308 :
309 2516965 : Real3 value = 0;
310 2516965 : Real33 grad = 0;
311 :
312 30885525 : for (unsigned int i = 0; i < num_dofs; ++i)
313 : {
314 28368560 : auto vector = getVectorDofValue(getElemLocalDofIndex(elem, i, var), tag);
315 :
316 28368560 : value += vector * phi(i, qp);
317 28368560 : grad += vector * grad_phi(i, qp);
318 : }
319 :
320 2516965 : grad = grad * jacobian;
321 :
322 2516965 : getVectorQpVectorValue(info, qp_offset + qp, var, tag) = value;
323 2516965 : getVectorQpVectorGrad(info, qp_offset + qp, var, tag) = grad;
324 2516965 : getVectorQpVectorCurl(info, qp_offset + qp, var, tag) =
325 5033930 : curlFromVectorGradient(grad, kokkosAssembly().getDimension());
326 : }
327 : else
328 : {
329 68145566 : auto & phi = kokkosAssembly().getPhi(sid, elem_type, fe_type);
330 68145566 : auto & grad_phi = kokkosAssembly().getGradPhi(sid, elem_type, fe_type);
331 68145566 : auto jacobian = kokkosAssembly().getJacobian(info, qp);
332 :
333 68145566 : Real value = 0;
334 68145566 : Real3 grad = 0;
335 :
336 372949230 : for (unsigned int i = 0; i < num_dofs; ++i)
337 : {
338 304803664 : auto vector = getVectorDofValue(getElemLocalDofIndex(elem, i, var), tag);
339 :
340 304803664 : value += vector * phi(i, qp);
341 304803664 : grad += vector * grad_phi(i, qp);
342 : }
343 :
344 68145566 : grad = jacobian * grad;
345 :
346 68145566 : getVectorQpValue(info, qp_offset + qp, var, tag) = value;
347 68145566 : getVectorQpGrad(info, qp_offset + qp, var, tag) = grad;
348 : }
349 : }
350 :
351 : } // namespace Moose::Kokkos
|
