Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
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 : // App includes
11 : #include "AugmentSparsityOnInterface.h"
12 : #include "Executioner.h"
13 : #include "FEProblemBase.h"
14 : #include "MooseApp.h"
15 :
16 : // libMesh includes
17 : #include "libmesh/elem.h"
18 : #include "libmesh/mesh_base.h"
19 : #include "libmesh/boundary_info.h"
20 :
21 : registerMooseObject("MooseApp", AugmentSparsityOnInterface);
22 :
23 : using namespace libMesh;
24 :
25 : InputParameters
26 32604 : AugmentSparsityOnInterface::validParams()
27 : {
28 32604 : InputParameters params = RelationshipManager::validParams();
29 32604 : params.addRequiredParam<BoundaryName>("primary_boundary",
30 : "The name of the primary boundary sideset.");
31 32604 : params.addRequiredParam<BoundaryName>("secondary_boundary",
32 : "The name of the secondary boundary sideset.");
33 32604 : params.addRequiredParam<SubdomainName>("primary_subdomain",
34 : "The name of the primary lower dimensional subdomain.");
35 32604 : params.addRequiredParam<SubdomainName>("secondary_subdomain",
36 : "The name of the secondary lower dimensional subdomain.");
37 97812 : params.addParam<bool>(
38 : "ghost_point_neighbors",
39 65208 : false,
40 : "Whether we should ghost point neighbors of secondary lower-dimensional elements and "
41 : "also their mortar interface couples for applications such as mortar nodal auxiliary "
42 : "kernels.");
43 97812 : params.addParam<bool>(
44 : "ghost_higher_d_neighbors",
45 65208 : false,
46 : "Whether we should ghost higher-dimensional neighbors. This is necessary when we are doing "
47 : "second order mortar with finite volume primal variables, because in order for the method to "
48 : "be second order we must use cell gradients, which couples in the neighbor cells.");
49 :
50 : // We want to wait until our mortar mesh has been built before trying to delete remote elements.
51 : // And our mortar mesh cannot be built until the entire mesh has been generated. By setting this
52 : // parameter to false we will make sure that any prepare_for_use calls during the mesh generation
53 : // phase will not delete remote elements *and* we will set a flag on the moose mesh saying that we
54 : // need to delete remote elements after the addition of late geometric ghosting functors
55 : // (including this ghosting functor)
56 32604 : params.set<bool>("attach_geometric_early") = false;
57 32604 : return params;
58 0 : }
59 :
60 7736 : AugmentSparsityOnInterface::AugmentSparsityOnInterface(const InputParameters & params)
61 : : RelationshipManager(params),
62 7736 : _primary_boundary_name(getParam<BoundaryName>("primary_boundary")),
63 7736 : _secondary_boundary_name(getParam<BoundaryName>("secondary_boundary")),
64 7736 : _primary_subdomain_name(getParam<SubdomainName>("primary_subdomain")),
65 7736 : _secondary_subdomain_name(getParam<SubdomainName>("secondary_subdomain")),
66 7736 : _is_coupling_functor(isType(Moose::RelationshipManagerType::COUPLING)),
67 7736 : _ghost_point_neighbors(getParam<bool>("ghost_point_neighbors")),
68 15472 : _ghost_higher_d_neighbors(getParam<bool>("ghost_higher_d_neighbors"))
69 : {
70 7736 : }
71 :
72 2867 : AugmentSparsityOnInterface::AugmentSparsityOnInterface(const AugmentSparsityOnInterface & other)
73 : : RelationshipManager(other),
74 2867 : _primary_boundary_name(other._primary_boundary_name),
75 2867 : _secondary_boundary_name(other._secondary_boundary_name),
76 2867 : _primary_subdomain_name(other._primary_subdomain_name),
77 2867 : _secondary_subdomain_name(other._secondary_subdomain_name),
78 2867 : _is_coupling_functor(other._is_coupling_functor),
79 2867 : _ghost_point_neighbors(other._ghost_point_neighbors),
80 2867 : _ghost_higher_d_neighbors(other._ghost_higher_d_neighbors)
81 : {
82 2867 : }
83 :
84 : void
85 2867 : AugmentSparsityOnInterface::internalInitWithMesh(const MeshBase &)
86 : {
87 2867 : }
88 :
89 : std::string
90 0 : AugmentSparsityOnInterface::getInfo() const
91 : {
92 0 : std::ostringstream oss;
93 0 : oss << "AugmentSparsityOnInterface";
94 0 : return oss.str();
95 0 : }
96 :
97 : void
98 348391 : AugmentSparsityOnInterface::ghostMortarInterfaceCouplings(
99 : const processor_id_type p,
100 : const Elem * const elem,
101 : map_type & coupled_elements,
102 : const AutomaticMortarGeneration & amg) const
103 : {
104 : // Look up elem in the mortar_interface_coupling data structure.
105 348391 : const auto & mic = amg.mortarInterfaceCoupling();
106 348391 : auto find_it = mic.find(elem->id());
107 348391 : if (find_it == mic.end())
108 302848 : return;
109 :
110 45543 : const auto & coupled_set = find_it->second;
111 :
112 227513 : for (const auto coupled_elem_id : coupled_set)
113 : {
114 181970 : const Elem * coupled_elem = _mesh->elem_ptr(coupled_elem_id);
115 : mooseAssert(coupled_elem,
116 : "The coupled element with id " << coupled_elem_id << " doesn't exist!");
117 :
118 181970 : if (coupled_elem->processor_id() != p)
119 144589 : coupled_elements.emplace(coupled_elem, _null_mat);
120 : }
121 : }
122 :
123 : void
124 350 : AugmentSparsityOnInterface::ghostLowerDSecondaryElemPointNeighbors(
125 : const processor_id_type p,
126 : const Elem * const query_elem,
127 : map_type & coupled_elements,
128 : const BoundaryID secondary_boundary_id,
129 : const SubdomainID secondary_subdomain_id,
130 : const AutomaticMortarGeneration & amg) const
131 : {
132 : // I hypothesize that node processor ids are tied to higher dimensional element processor
133 : // ids over lower dimensional element processor ids based on debugging experience.
134 : // Consequently we need to be checking for higher dimensional elements and whether they are
135 : // along our secondary boundary. From there we will query the AMG object's
136 : // mortar-interface-coupling container to get the secondary lower-dimensional element, and
137 : // then we will ghost it's point neighbors and their mortar interface couples
138 :
139 : // It's possible that one higher-dimensional element could have multiple lower-dimensional
140 : // elements from multiple sides. Morever, even if there is only one lower-d element per
141 : // higher-d element, the unordered_multimap that holds the coupling information can have
142 : // duplicate key-value pairs if there are multiple mortar segments per secondary face. So to
143 : // prevent attempting to insert into the coupled elements map multiple times with the same
144 : // element, we'll keep track of the elements we've handled. We're going to use a tree-based
145 : // set here since the number of lower-d elements handled should never exceed the number of
146 : // element sides (which is small)
147 350 : std::set<dof_id_type> secondary_lower_elems_handled;
148 350 : const BoundaryInfo & binfo = _mesh->get_boundary_info();
149 980 : for (auto side : query_elem->side_index_range())
150 : {
151 735 : if (!binfo.has_boundary_id(query_elem, side, secondary_boundary_id))
152 : // We're not a higher-dimensional element along the secondary face, or at least this
153 : // side isn't
154 630 : continue;
155 :
156 105 : const auto & mic = amg.mortarInterfaceCoupling();
157 105 : auto find_it = mic.find(query_elem->id());
158 105 : if (find_it == mic.end())
159 0 : continue;
160 :
161 105 : const auto & coupled_set = find_it->second;
162 350 : for (const auto coupled_elem_id : coupled_set)
163 : {
164 245 : auto * const coupled_elem = _mesh->elem_ptr(coupled_elem_id);
165 :
166 245 : if (coupled_elem->subdomain_id() != secondary_subdomain_id)
167 : {
168 : // We support higher-d-secondary to higher-d-primary coupling now, e.g.
169 : // if we get here, coupled_elem is not actually a secondary lower elem; it's a
170 : // primary higher-d elem
171 : mooseAssert(coupled_elem->dim() == query_elem->dim(), "These should be matching dim");
172 140 : continue;
173 : }
174 :
175 105 : auto insert_pr = secondary_lower_elems_handled.insert(coupled_elem_id);
176 :
177 : // If insertion didn't happen, then we've already handled this element
178 105 : if (!insert_pr.second)
179 0 : continue;
180 :
181 : // We've already ghosted the secondary lower-d element itself if it needed to be
182 : // outside of the _ghost_point_neighbors logic. But now we must make sure to ghost the
183 : // point neighbors of the secondary lower-d element and their mortar interface
184 : // couplings
185 105 : std::set<const Elem *> secondary_lower_elem_point_neighbors;
186 105 : coupled_elem->find_point_neighbors(secondary_lower_elem_point_neighbors);
187 :
188 350 : for (const Elem * const neigh : secondary_lower_elem_point_neighbors)
189 : {
190 245 : if (neigh->processor_id() != p)
191 0 : coupled_elements.emplace(neigh, _null_mat);
192 :
193 245 : ghostMortarInterfaceCouplings(p, neigh, coupled_elements, amg);
194 : }
195 105 : } // end iteration over mortar interface couplings
196 :
197 : // We actually should have added all the lower-dimensional elements associated with the
198 : // higher-dimensional element, so we can stop iterating over sides
199 105 : return;
200 :
201 : } // end for side_index_range
202 350 : }
203 :
204 : void
205 78204 : AugmentSparsityOnInterface::ghostHigherDNeighbors(const processor_id_type p,
206 : const Elem * const query_elem,
207 : map_type & coupled_elements,
208 : const BoundaryID secondary_boundary_id,
209 : const SubdomainID secondary_subdomain_id,
210 : const AutomaticMortarGeneration & amg) const
211 : {
212 : // The coupling is this for second order FV: secondary lower dimensional elem dofs will depend on
213 : // higher-d secondary elem neighbor primal dofs and higher-d primary elem neighbor primal dofs.
214 : // Higher dimensional primal dofs only depend on the secondary lower dimensional LM dofs for the
215 : // current FV gap heat transfer use cases
216 :
217 78204 : if (query_elem->subdomain_id() != secondary_subdomain_id)
218 77322 : return;
219 :
220 882 : const BoundaryInfo & binfo = _mesh->get_boundary_info();
221 882 : const auto which_side = query_elem->interior_parent()->which_side_am_i(query_elem);
222 882 : if (!binfo.has_boundary_id(query_elem->interior_parent(), which_side, secondary_boundary_id))
223 0 : return;
224 :
225 882 : const auto & mic = amg.mortarInterfaceCoupling();
226 882 : auto find_it = mic.find(query_elem->id());
227 882 : if (find_it == mic.end())
228 : // Perhaps no projection onto primary
229 0 : return;
230 :
231 882 : const auto & lower_d_coupled_set = find_it->second;
232 3528 : for (const auto coupled_elem_id : lower_d_coupled_set)
233 : {
234 2646 : auto * const coupled_elem = _mesh->elem_ptr(coupled_elem_id);
235 2646 : if (coupled_elem->dim() == query_elem->dim())
236 : // lower-d-elem
237 882 : continue;
238 :
239 1764 : std::vector<const Elem *> active_neighbors;
240 :
241 8820 : for (auto s : coupled_elem->side_index_range())
242 : {
243 7056 : const Elem * const neigh = coupled_elem->neighbor_ptr(s);
244 :
245 7056 : if (!neigh || neigh == remote_elem)
246 1932 : continue;
247 :
248 : // With any kind of neighbor, we need to couple to all the
249 : // active descendants on our side.
250 5124 : neigh->active_family_tree_by_neighbor(active_neighbors, coupled_elem);
251 :
252 10248 : for (const auto & neighbor : active_neighbors)
253 5124 : if (neighbor->processor_id() != p)
254 : {
255 : mooseAssert(
256 : neighbor->subdomain_id() != secondary_subdomain_id,
257 : "Ensure that we aren't missing potential erasures from the secondary-to-msms map");
258 5124 : coupled_elements.emplace(neighbor, _null_mat);
259 : }
260 : }
261 1764 : }
262 : }
263 :
264 : void
265 241457 : AugmentSparsityOnInterface::operator()(const MeshBase::const_element_iterator & range_begin,
266 : const MeshBase::const_element_iterator & range_end,
267 : const processor_id_type p,
268 : map_type & coupled_elements)
269 : {
270 : // We ask the user to pass boundary names instead of ids to our constraint object. However, We
271 : // are unable to get the boundary ids from boundary names until we've attached the MeshBase object
272 : // to the MooseMesh
273 241457 : const bool generating_mesh = !_moose_mesh->getMeshPtr();
274 : const auto primary_boundary_id = generating_mesh
275 241457 : ? Moose::INVALID_BOUNDARY_ID
276 241457 : : _moose_mesh->getBoundaryID(_primary_boundary_name);
277 : const auto secondary_boundary_id = generating_mesh
278 241457 : ? Moose::INVALID_BOUNDARY_ID
279 241457 : : _moose_mesh->getBoundaryID(_secondary_boundary_name);
280 : const auto primary_subdomain_id = generating_mesh
281 241457 : ? Moose::INVALID_BLOCK_ID
282 241457 : : _moose_mesh->getSubdomainID(_primary_subdomain_name);
283 : const auto secondary_subdomain_id = generating_mesh
284 241457 : ? Moose::INVALID_BLOCK_ID
285 241457 : : _moose_mesh->getSubdomainID(_secondary_subdomain_name);
286 :
287 : const AutomaticMortarGeneration * const amg =
288 482914 : _app.getExecutioner() ? &_app.getExecutioner()->feProblem().getMortarInterface(
289 241457 : std::make_pair(primary_boundary_id, secondary_boundary_id),
290 0 : std::make_pair(primary_subdomain_id, secondary_subdomain_id),
291 241457 : _use_displaced_mesh)
292 241457 : : nullptr;
293 :
294 : // If we're on a dynamic mesh or we have not yet constructed the mortar mesh, we need to ghost the
295 : // entire interface because we don't know a priori what elements will project onto what. We *do
296 : // not* add the whole interface if we are a coupling functor because it is very expensive. This is
297 : // because when building the sparsity pattern, we call through to the ghosting functors with one
298 : // element at a time (and then below we do a loop over all the mesh's active elements). It's
299 : // perhaps faster in this case to deal with mallocs coming out of MatSetValues, especially if the
300 : // mesh displacements are relatively small
301 241457 : if ((!amg || _use_displaced_mesh) && !_is_coupling_functor)
302 : {
303 104972 : for (const Elem * const elem : _mesh->active_element_ptr_range())
304 : {
305 52288 : if (generating_mesh)
306 : {
307 : // We are still generating the mesh, so it's possible we don't even have the right boundary
308 : // ids created yet! So we actually ghost all boundary elements and all lower dimensional
309 : // elements who have parents on a boundary
310 0 : if (elem->on_boundary())
311 0 : coupled_elements.insert(std::make_pair(elem, _null_mat));
312 0 : else if (const Elem * const ip = elem->interior_parent())
313 : {
314 0 : if (ip->on_boundary())
315 0 : coupled_elements.insert(std::make_pair(elem, _null_mat));
316 : }
317 : }
318 : else
319 : {
320 : // We've finished generating our mesh so we can be selective and only ghost elements lying
321 : // in our lower-dimensional subdomains and their interior parents
322 :
323 : mooseAssert(primary_boundary_id != Moose::INVALID_BOUNDARY_ID,
324 : "Primary boundary id should exist by now.");
325 : mooseAssert(secondary_boundary_id != Moose::INVALID_BOUNDARY_ID,
326 : "Secondary boundary id should exist by now.");
327 : mooseAssert(primary_subdomain_id != Moose::INVALID_BLOCK_ID,
328 : "Primary subdomain id should exist by now.");
329 : mooseAssert(secondary_subdomain_id != Moose::INVALID_BLOCK_ID,
330 : "Secondary subdomain id should exist by now.");
331 :
332 : // Higher-dimensional boundary elements
333 52288 : const BoundaryInfo & binfo = _mesh->get_boundary_info();
334 :
335 212180 : for (auto side : elem->side_index_range())
336 239824 : if ((elem->processor_id() != p) &&
337 79932 : (binfo.has_boundary_id(elem, side, primary_boundary_id) ||
338 78588 : binfo.has_boundary_id(elem, side, secondary_boundary_id)))
339 2808 : coupled_elements.insert(std::make_pair(elem, _null_mat));
340 :
341 : // Lower dimensional subdomain elements
342 77152 : if ((elem->processor_id() != p) && (elem->subdomain_id() == primary_subdomain_id ||
343 24864 : elem->subdomain_id() == secondary_subdomain_id))
344 : {
345 2808 : coupled_elements.insert(std::make_pair(elem, _null_mat));
346 :
347 : #ifndef NDEBUG
348 : // let's do some safety checks
349 : const Elem * const ip = elem->interior_parent();
350 : mooseAssert(ip,
351 : "We should have set interior parents for all of our lower-dimensional mortar "
352 : "subdomains");
353 : auto side = ip->which_side_am_i(elem);
354 : auto bnd_id = elem->subdomain_id() == primary_subdomain_id ? primary_boundary_id
355 : : secondary_boundary_id;
356 : mooseAssert(_mesh->get_boundary_info().has_boundary_id(ip, side, bnd_id),
357 : "The interior parent for the lower-dimensional element does not lie on the "
358 : "boundary");
359 : #endif
360 : }
361 : }
362 396 : }
363 396 : }
364 : // For a static mesh (or for determining a sparsity pattern approximation on a displaced mesh) we
365 : // can just ghost the coupled elements determined during mortar mesh generation
366 241061 : else if (amg)
367 : {
368 937353 : for (const Elem * const elem : as_range(range_begin, range_end))
369 : {
370 348146 : ghostMortarInterfaceCouplings(p, elem, coupled_elements, *amg);
371 :
372 348146 : if (_ghost_point_neighbors)
373 350 : ghostLowerDSecondaryElemPointNeighbors(
374 : p, elem, coupled_elements, secondary_boundary_id, secondary_subdomain_id, *amg);
375 348146 : if (_ghost_higher_d_neighbors)
376 78204 : ghostHigherDNeighbors(
377 : p, elem, coupled_elements, secondary_boundary_id, secondary_subdomain_id, *amg);
378 241061 : } // end for loop over input range
379 : } // end if amg
380 241457 : }
381 :
382 : bool
383 41294 : AugmentSparsityOnInterface::operator>=(const RelationshipManager & other) const
384 : {
385 41294 : if (auto asoi = dynamic_cast<const AugmentSparsityOnInterface *>(&other))
386 : {
387 34366 : if (_primary_boundary_name == asoi->_primary_boundary_name &&
388 26580 : _secondary_boundary_name == asoi->_secondary_boundary_name &&
389 26580 : _primary_subdomain_name == asoi->_primary_subdomain_name &&
390 13290 : _secondary_subdomain_name == asoi->_secondary_subdomain_name &&
391 13290 : _ghost_point_neighbors >= asoi->_ghost_point_neighbors &&
392 34366 : _ghost_higher_d_neighbors >= asoi->_ghost_higher_d_neighbors && baseGreaterEqual(*asoi))
393 5890 : return true;
394 : }
395 35404 : return false;
396 : }
397 :
398 : std::unique_ptr<GhostingFunctor>
399 2867 : AugmentSparsityOnInterface::clone() const
400 : {
401 2867 : return _app.getFactory().copyConstruct(*this);
402 : }
|
