Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : #include "libmesh/petsc_macro.h"
21 :
22 : #ifdef LIBMESH_HAVE_PETSC
23 :
24 : #include "libmesh/ignore_warnings.h"
25 :
26 : // PETSc includes
27 : #if !PETSC_VERSION_LESS_THAN(3,6,0)
28 : # include <petsc/private/dmimpl.h>
29 : #else
30 : # include <petsc-private/dmimpl.h>
31 : #endif
32 :
33 : #include "libmesh/restore_warnings.h"
34 :
35 : // Local Includes
36 : #include "libmesh/libmesh_common.h"
37 : #include "libmesh/nonlinear_implicit_system.h"
38 : #include "libmesh/petsc_nonlinear_solver.h"
39 : #include "libmesh/petsc_linear_solver.h"
40 : #include "libmesh/petsc_vector.h"
41 : #include "libmesh/petsc_matrix_base.h"
42 : #include "libmesh/petscdmlibmesh.h"
43 : #include "libmesh/dof_map.h"
44 : #include "libmesh/constraint_elimination_dof_map.h"
45 : #include "libmesh/preconditioner.h"
46 : #include "libmesh/elem.h"
47 : #include "libmesh/parallel.h"
48 :
49 :
50 : using namespace libMesh;
51 :
52 :
53 : #define DMLIBMESH_NO_DECOMPOSITION 0
54 : #define DMLIBMESH_FIELD_DECOMPOSITION 1
55 : #define DMLIBMESH_DOMAIN_DECOMPOSITION 2
56 :
57 : #define DMLIBMESH_NO_EMBEDDING 0
58 : #define DMLIBMESH_FIELD_EMBEDDING 1
59 : #define DMLIBMESH_DOMAIN_EMBEDDING 2
60 :
61 : struct DM_libMesh
62 : {
63 : NonlinearImplicitSystem * sys;
64 : std::map<std::string, unsigned int> * varids;
65 : std::map<unsigned int, std::string> * varnames;
66 : std::map<std::string, unsigned int> * blockids;
67 : std::map<unsigned int, std::string> * blocknames;
68 : unsigned int decomposition_type;
69 : std::vector<std::set<unsigned int>> * decomposition;
70 : unsigned int embedding_type;
71 : IS embedding;
72 : unsigned int vec_count;
73 : };
74 :
75 : struct DMVec_libMesh {
76 : std::string label;
77 : };
78 :
79 : #undef __FUNCT__
80 : #define __FUNCT__ "DMlibMeshGetVec_Private"
81 0 : PetscErrorCode DMlibMeshGetVec_Private(DM /*dm*/, const char * /*name*/, Vec *)
82 : {
83 : PetscFunctionBegin;
84 :
85 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
86 : }
87 :
88 :
89 :
90 : PetscErrorCode DMlibMeshSetUpName_Private(DM dm);
91 :
92 : #undef __FUNCT__
93 : #define __FUNCT__ "DMlibMeshSetSystem_libMesh"
94 2574 : PetscErrorCode DMlibMeshSetSystem_libMesh(DM dm, NonlinearImplicitSystem & sys)
95 : {
96 1716 : const Parallel::Communicator & comm = sys.comm();
97 :
98 : PetscFunctionBegin;
99 : PetscValidHeaderSpecific(dm,DM_CLASSID,1);
100 : PetscBool islibmesh;
101 2574 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH,&islibmesh));
102 2574 : if (!islibmesh) LIBMESH_SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Got DM of type %s, not of type %s", ((PetscObject)dm)->type_name, DMLIBMESH);
103 :
104 2574 : if (dm->setupcalled) SETERRQ(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONGSTATE, "Cannot reset the libMesh system after DM has been set up.");
105 2574 : DM_libMesh * dlm = (DM_libMesh *)(dm->data);
106 2574 : dlm->sys =&sys;
107 : /* Initially populate the sets of active blockids and varids using all of the
108 : existing blocks/variables (only variables are supported at the moment). */
109 858 : DofMap & dofmap = dlm->sys->get_dof_map();
110 2574 : dlm->varids->clear();
111 2574 : dlm->varnames->clear();
112 5526 : for (auto v : make_range(dofmap.n_variables())) {
113 2952 : std::string vname = dofmap.variable(v).name();
114 2952 : dlm->varids->insert(std::pair<std::string,unsigned int>(vname,v));
115 4920 : dlm->varnames->insert(std::pair<unsigned int,std::string>(v,vname));
116 : }
117 2574 : const MeshBase & mesh = dlm->sys->get_mesh();
118 2574 : dlm->blockids->clear();
119 2574 : dlm->blocknames->clear();
120 1716 : std::set<subdomain_id_type> blocks;
121 : /* The following effectively is a verbatim copy of MeshBase::n_subdomains(). */
122 : // This requires an inspection on every processor
123 858 : libmesh_parallel_only(mesh.comm());
124 3590898 : for (const auto & elem : mesh.active_element_ptr_range())
125 2690814 : blocks.insert(elem->subdomain_id());
126 : // Some subdomains may only live on other processors
127 2574 : comm.set_union(blocks);
128 :
129 858 : std::set<subdomain_id_type>::iterator bit = blocks.begin();
130 858 : std::set<subdomain_id_type>::iterator bend = blocks.end();
131 2574 : if (bit == bend) SETERRQ(((PetscObject)dm)->comm, PETSC_ERR_PLIB, "No mesh blocks found.");
132 :
133 5196 : for (; bit != bend; ++bit) {
134 2622 : subdomain_id_type bid = *bit;
135 2622 : std::string bname = mesh.subdomain_name(bid);
136 2622 : if (!bname.length()) {
137 : /* Block names are currently implemented for Exodus II meshes
138 : only, so we might have to make up our own block names and
139 : maintain our own mapping of block ids to names.
140 : */
141 2622 : std::ostringstream ss;
142 1748 : ss << "dm" << bid;
143 3496 : bname = ss.str();
144 874 : }
145 2622 : dlm->blockids->insert(std::pair<std::string,unsigned int>(bname,bid));
146 4370 : dlm->blocknames->insert(std::pair<unsigned int,std::string>(bid,bname));
147 : }
148 2574 : LibmeshPetscCallQ(DMlibMeshSetUpName_Private(dm));
149 858 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
150 : }
151 :
152 : #undef __FUNCT__
153 : #define __FUNCT__ "DMlibMeshGetSystem_libMesh"
154 1224 : PetscErrorCode DMlibMeshGetSystem_libMesh(DM dm, NonlinearImplicitSystem *& sys)
155 : {
156 : PetscFunctionBegin;
157 : PetscValidHeaderSpecific(dm,DM_CLASSID,1);
158 : PetscBool islibmesh;
159 1224 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH,&islibmesh));
160 1224 : if (!islibmesh) LIBMESH_SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Got DM of type %s, not of type %s", ((PetscObject)dm)->type_name, DMLIBMESH);
161 1224 : DM_libMesh * dlm = (DM_libMesh *)(dm->data);
162 1224 : sys = dlm->sys;
163 1224 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
164 : }
165 :
166 :
167 : #undef __FUNCT__
168 : #define __FUNCT__ "DMlibMeshGetBlocks"
169 0 : PetscErrorCode DMlibMeshGetBlocks(DM dm, PetscInt * n, char *** blocknames)
170 : {
171 : PetscInt i;
172 : PetscFunctionBegin;
173 : PetscValidHeaderSpecific(dm,DM_CLASSID,1);
174 : PetscBool islibmesh;
175 0 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH,&islibmesh));
176 0 : if (!islibmesh) LIBMESH_SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Got DM of type %s, not of type %s", ((PetscObject)dm)->type_name, DMLIBMESH);
177 0 : DM_libMesh * dlm = (DM_libMesh *)(dm->data);
178 : #if PETSC_RELEASE_GREATER_EQUALS(3,20,0)
179 : PetscAssertPointer(n,2);
180 : #else
181 : PetscValidPointer(n,2);
182 : #endif
183 0 : *n = cast_int<unsigned int>(dlm->blockids->size());
184 0 : if (!blocknames) PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
185 0 : LibmeshPetscCallQ(PetscMalloc(*n*sizeof(char *), blocknames));
186 0 : i = 0;
187 0 : for (const auto & pr : *(dlm->blockids))
188 : {
189 0 : LibmeshPetscCallQ(PetscStrallocpy(pr.first.c_str(), *blocknames+i));
190 0 : ++i;
191 : }
192 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
193 : }
194 :
195 : #undef __FUNCT__
196 : #define __FUNCT__ "DMlibMeshGetVariables"
197 0 : PetscErrorCode DMlibMeshGetVariables(DM dm, PetscInt * n, char *** varnames)
198 : {
199 : PetscFunctionBegin;
200 : PetscValidHeaderSpecific(dm,DM_CLASSID,1);
201 : PetscBool islibmesh;
202 : PetscInt i;
203 0 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH,&islibmesh));
204 0 : if (!islibmesh) LIBMESH_SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Got DM of type %s, not of type %s", ((PetscObject)dm)->type_name, DMLIBMESH);
205 0 : DM_libMesh * dlm = (DM_libMesh *)(dm->data);
206 : #if PETSC_RELEASE_GREATER_EQUALS(3,20,0)
207 : PetscAssertPointer(n,2);
208 : #else
209 : PetscValidPointer(n,2);
210 : #endif
211 0 : *n = cast_int<unsigned int>(dlm->varids->size());
212 0 : if (!varnames) PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
213 0 : LibmeshPetscCallQ(PetscMalloc(*n*sizeof(char *), varnames));
214 0 : i = 0;
215 0 : for (const auto & pr : *(dlm->varids))
216 : {
217 0 : LibmeshPetscCallQ(PetscStrallocpy(pr.first.c_str(), *varnames+i));
218 0 : ++i;
219 : }
220 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
221 : }
222 :
223 : #undef __FUNCT__
224 : #define __FUNCT__ "DMlibMeshSetUpName_Private"
225 2580 : PetscErrorCode DMlibMeshSetUpName_Private(DM dm)
226 : {
227 2580 : DM_libMesh * dlm = (DM_libMesh *)dm->data;
228 : PetscFunctionBegin;
229 2580 : std::string name = dlm->sys->name();
230 860 : std::map<unsigned int, std::string> * dnames = LIBMESH_PETSC_NULLPTR,
231 860 : * enames = LIBMESH_PETSC_NULLPTR;
232 2580 : if (dlm->decomposition_type == DMLIBMESH_FIELD_DECOMPOSITION) {
233 2 : name += ":dec:var:";
234 6 : dnames = dlm->varnames;
235 : }
236 2580 : if (dlm->decomposition_type == DMLIBMESH_DOMAIN_DECOMPOSITION) {
237 0 : name += ":dec:block:";
238 0 : dnames = dlm->blocknames;
239 : }
240 2580 : if (dnames) {
241 18 : for (auto decomp : *dlm->decomposition) {
242 12 : for (std::set<unsigned int>::iterator dit_begin = decomp.begin(),
243 4 : dit = dit_begin,
244 4 : dit_end = decomp.end();
245 28 : dit != dit_end; ++dit) {
246 12 : unsigned int id = *dit;
247 12 : if (dit != dit_begin)
248 0 : name += ",";
249 12 : name += (*dnames)[id];
250 : }
251 4 : name += ";";
252 : }
253 : }
254 2580 : if (dlm->embedding_type == DMLIBMESH_FIELD_EMBEDDING) {
255 0 : name += ":emb:var:";
256 0 : enames = dlm->varnames;
257 : }
258 2580 : if (dlm->embedding_type == DMLIBMESH_DOMAIN_EMBEDDING) {
259 0 : name += ":emb:block:";
260 0 : enames = dlm->blocknames;
261 : }
262 2580 : if (enames) {
263 0 : for (auto eit = enames->begin(),
264 0 : eit_end = enames->end(); eit != eit_end; ++eit) {
265 0 : std::string & ename = eit->second;
266 0 : if (eit != enames->begin())
267 0 : name += ",";
268 0 : name += ename;
269 : }
270 0 : name += ";";
271 : }
272 2580 : LibmeshPetscCallQ(PetscObjectSetName((PetscObject)dm, name.c_str()));
273 3440 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
274 : }
275 :
276 :
277 : /**
278 : * Inserts into \p indices the local indices, in the numbering of the vectors SNES solves on, of
279 : * variable \p var on \p elem: reduced indices when \p sys eliminates constrained degrees of
280 : * freedom, and full ones otherwise
281 : */
282 114 : static void DMlibMeshInsertLocalDofs_Private(const NonlinearImplicitSystem & sys,
283 : const Elem * elem,
284 : unsigned int var,
285 : std::set<numeric_index_type> & indices)
286 : {
287 38 : const DofMap & dof_map = sys.get_dof_map();
288 76 : std::vector<dof_id_type> evindices;
289 114 : dof_map.dof_indices(elem, evindices, var);
290 :
291 : const ConstraintEliminationDofMap * reduced =
292 76 : dof_map.has_constraint_elimination() ? &dof_map.get_constraint_elimination() : nullptr;
293 :
294 855 : for (const auto dof : evindices)
295 494 : if (dof_map.local_index(dof))
296 : {
297 660 : if (!reduced)
298 0 : indices.insert(dof);
299 660 : else if (const auto r = reduced->find_reduced_dof(dof); r != DofObject::invalid_id)
300 370 : indices.insert(r);
301 : }
302 114 : }
303 :
304 : #undef __FUNCT__
305 : #define __FUNCT__ "DMCreateFieldDecomposition_libMesh"
306 6 : static PetscErrorCode DMCreateFieldDecomposition_libMesh(DM dm, PetscInt * len, char *** namelist, IS ** islist, DM ** dmlist)
307 : {
308 : PetscFunctionBegin;
309 6 : DM_libMesh * dlm = (DM_libMesh *)(dm->data);
310 6 : NonlinearImplicitSystem * sys = dlm->sys;
311 : IS emb;
312 6 : if (dlm->decomposition_type != DMLIBMESH_FIELD_DECOMPOSITION) PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
313 :
314 6 : *len = cast_int<unsigned int>(dlm->decomposition->size());
315 6 : if (namelist) {LibmeshPetscCallQ(PetscMalloc(*len*sizeof(char *), namelist));}
316 6 : if (islist) {LibmeshPetscCallQ(PetscMalloc(*len*sizeof(IS), islist));}
317 6 : if (dmlist) {LibmeshPetscCallQ(PetscMalloc(*len*sizeof(DM), dmlist));}
318 18 : for (auto d : index_range(*dlm->decomposition)) {
319 4 : std::set<numeric_index_type> dindices;
320 4 : std::string dname;
321 4 : std::map<std::string, unsigned int> dvarids;
322 4 : std::map<unsigned int, std::string> dvarnames;
323 4 : unsigned int dvcount = 0;
324 24 : for (const auto & v : (*dlm->decomposition)[d]) {
325 12 : std::string vname = (*dlm->varnames)[v];
326 16 : dvarids.insert(std::pair<std::string, unsigned int>(vname,v));
327 12 : dvarnames.insert(std::pair<unsigned int,std::string>(v,vname));
328 12 : if (!dvcount) dname = vname;
329 0 : else dname += "_" + vname;
330 12 : ++dvcount;
331 12 : if (!islist) continue;
332 : // Iterate only over this DM's blocks.
333 24 : for (const auto & pr : *(dlm->blockids)) {
334 12 : const subdomain_id_type sbd_id = cast_int<subdomain_id_type>(pr.second);
335 38 : for (const auto & elem :
336 20 : as_range(sys->get_mesh().active_local_subdomain_elements_begin(sbd_id),
337 230 : sys->get_mesh().active_local_subdomain_elements_end(sbd_id))) {
338 114 : DMlibMeshInsertLocalDofs_Private(*sys, elem, v, dindices);
339 4 : }
340 : }
341 : }
342 12 : if (namelist) {
343 0 : LibmeshPetscCallQ(PetscStrallocpy(dname.c_str(),(*namelist)+d));
344 : }
345 12 : if (islist) {
346 : IS dis;
347 : PetscInt * darray;
348 16 : LibmeshPetscCallQ(PetscMalloc(sizeof(PetscInt)*dindices.size(), &darray));
349 4 : numeric_index_type i = 0;
350 315 : for (const auto & id : dindices) {
351 303 : darray[i] = id;
352 303 : ++i;
353 : }
354 16 : LibmeshPetscCallQ(ISCreateGeneral(((PetscObject)dm)->comm,
355 : cast_int<PetscInt>(dindices.size()),
356 : darray, PETSC_OWN_POINTER, &dis));
357 12 : if (dlm->embedding) {
358 : /* Create a relative embedding into the parent's index space. */
359 0 : LibmeshPetscCallQ(ISEmbed(dis,dlm->embedding, PETSC_TRUE, &emb));
360 : PetscInt elen, dlen;
361 0 : LibmeshPetscCallQ(ISGetLocalSize(emb, &elen));
362 0 : LibmeshPetscCallQ(ISGetLocalSize(dis, &dlen));
363 0 : if (elen != dlen) LIBMESH_SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_PLIB, "Failed to embed subdomain %zu", d);
364 0 : LibmeshPetscCallQ(ISDestroy(&dis));
365 0 : dis = emb;
366 : }
367 : else {
368 12 : emb = dis;
369 : }
370 12 : (*islist)[d] = dis;
371 : }
372 12 : if (dmlist) {
373 : DM ddm;
374 0 : LibmeshPetscCallQ(DMCreate(((PetscObject)dm)->comm, &ddm));
375 0 : LibmeshPetscCallQ(DMSetType(ddm, DMLIBMESH));
376 0 : DM_libMesh * ddlm = (DM_libMesh *)(ddm->data);
377 0 : ddlm->sys = dlm->sys;
378 : /* copy over the block ids and names */
379 0 : *ddlm->blockids = *dlm->blockids;
380 0 : *ddlm->blocknames = *dlm->blocknames;
381 : /* set the vars from the d-th part of the decomposition. */
382 0 : *ddlm->varids = dvarids;
383 0 : *ddlm->varnames = dvarnames;
384 0 : LibmeshPetscCallQ(PetscObjectReference((PetscObject)emb));
385 0 : ddlm->embedding = emb;
386 0 : ddlm->embedding_type = DMLIBMESH_FIELD_EMBEDDING;
387 :
388 0 : LibmeshPetscCallQ(DMlibMeshSetUpName_Private(ddm));
389 0 : LibmeshPetscCallQ(DMSetFromOptions(ddm));
390 0 : (*dmlist)[d] = ddm;
391 : }
392 : }
393 2 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
394 : }
395 :
396 : #undef __FUNCT__
397 : #define __FUNCT__ "DMCreateDomainDecomposition_libMesh"
398 0 : static PetscErrorCode DMCreateDomainDecomposition_libMesh(DM dm, PetscInt * len, char *** namelist, IS ** innerislist, IS ** outerislist, DM ** dmlist)
399 : {
400 : PetscFunctionBegin;
401 0 : DM_libMesh * dlm = (DM_libMesh *)(dm->data);
402 0 : NonlinearImplicitSystem * sys = dlm->sys;
403 : IS emb;
404 0 : if (dlm->decomposition_type != DMLIBMESH_DOMAIN_DECOMPOSITION) PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
405 0 : *len = cast_int<unsigned int>(dlm->decomposition->size());
406 0 : if (namelist) {LibmeshPetscCallQ(PetscMalloc(*len*sizeof(char *), namelist));}
407 0 : if (innerislist) {LibmeshPetscCallQ(PetscMalloc(*len*sizeof(IS), innerislist));}
408 0 : if (outerislist) *outerislist = LIBMESH_PETSC_NULLPTR; /* FIX: allow mesh-based overlap. */
409 0 : if (dmlist) {LibmeshPetscCallQ(PetscMalloc(*len*sizeof(DM), dmlist));}
410 0 : for (auto d : index_range(*dlm->decomposition)) {
411 0 : std::set<numeric_index_type> dindices;
412 0 : std::string dname;
413 0 : std::map<std::string, unsigned int> dblockids;
414 0 : std::map<unsigned int,std::string> dblocknames;
415 0 : unsigned int dbcount = 0;
416 0 : for (const auto & b : (*dlm->decomposition)[d]) {
417 0 : std::string bname = (*dlm->blocknames)[b];
418 0 : dblockids.insert(std::pair<std::string, unsigned int>(bname,b));
419 0 : dblocknames.insert(std::pair<unsigned int,std::string>(b,bname));
420 0 : if (!dbcount) dname = bname;
421 0 : else dname += "_" + bname;
422 0 : ++dbcount;
423 0 : if (!innerislist) continue;
424 0 : const subdomain_id_type b_sbd_id = cast_int<subdomain_id_type>(b);
425 0 : MeshBase::const_element_iterator el = sys->get_mesh().active_local_subdomain_elements_begin(b_sbd_id);
426 0 : const MeshBase::const_element_iterator end_el = sys->get_mesh().active_local_subdomain_elements_end(b_sbd_id);
427 0 : for ( ; el != end_el; ++el) {
428 0 : const Elem * elem = *el;
429 : // Iterate only over this DM's variables.
430 0 : for (const auto & pr : *(dlm->varids))
431 0 : DMlibMeshInsertLocalDofs_Private(*sys, elem, pr.second, dindices);
432 : }
433 : }
434 0 : if (namelist) {
435 0 : LibmeshPetscCallQ(PetscStrallocpy(dname.c_str(),(*namelist)+d));
436 : }
437 0 : if (innerislist) {
438 : PetscInt * darray;
439 : IS dis;
440 0 : LibmeshPetscCallQ(PetscMalloc(sizeof(PetscInt)*dindices.size(), &darray));
441 0 : numeric_index_type i = 0;
442 0 : for (const auto & id : dindices) {
443 0 : darray[i] = id;
444 0 : ++i;
445 : }
446 0 : LibmeshPetscCallQ(ISCreateGeneral(((PetscObject)dm)->comm,
447 : cast_int<PetscInt>(dindices.size()),
448 : darray, PETSC_OWN_POINTER, &dis));
449 0 : if (dlm->embedding) {
450 : /* Create a relative embedding into the parent's index space. */
451 0 : LibmeshPetscCallQ(ISEmbed(dis,dlm->embedding, PETSC_TRUE, &emb));
452 : PetscInt elen, dlen;
453 0 : LibmeshPetscCallQ(ISGetLocalSize(emb, &elen));
454 0 : LibmeshPetscCallQ(ISGetLocalSize(dis, &dlen));
455 0 : if (elen != dlen) LIBMESH_SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_PLIB, "Failed to embed field %zu" , d);
456 0 : LibmeshPetscCallQ(ISDestroy(&dis));
457 0 : dis = emb;
458 : }
459 : else {
460 0 : emb = dis;
461 : }
462 0 : if (innerislist) {
463 0 : LibmeshPetscCallQ(PetscObjectReference((PetscObject)dis));
464 0 : (*innerislist)[d] = dis;
465 : }
466 0 : LibmeshPetscCallQ(ISDestroy(&dis));
467 : }
468 0 : if (dmlist) {
469 : DM ddm;
470 0 : LibmeshPetscCallQ(DMCreate(((PetscObject)dm)->comm, &ddm));
471 0 : LibmeshPetscCallQ(DMSetType(ddm, DMLIBMESH));
472 0 : DM_libMesh * ddlm = (DM_libMesh *)(ddm->data);
473 0 : ddlm->sys = dlm->sys;
474 : /* copy over the varids and varnames */
475 0 : *ddlm->varids = *dlm->varids;
476 0 : *ddlm->varnames = *dlm->varnames;
477 : /* set the blocks from the d-th part of the decomposition. */
478 0 : *ddlm->blockids = dblockids;
479 0 : *ddlm->blocknames = dblocknames;
480 0 : LibmeshPetscCallQ(PetscObjectReference((PetscObject)emb));
481 0 : ddlm->embedding = emb;
482 0 : ddlm->embedding_type = DMLIBMESH_DOMAIN_EMBEDDING;
483 :
484 0 : LibmeshPetscCallQ(DMlibMeshSetUpName_Private(ddm));
485 0 : LibmeshPetscCallQ(DMSetFromOptions(ddm));
486 0 : (*dmlist)[d] = ddm;
487 : }
488 : }
489 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
490 : }
491 :
492 :
493 : #undef __FUNCT__
494 : #define __FUNCT__ "DMlibMeshCreateFieldDecompositionDM"
495 6 : PetscErrorCode DMlibMeshCreateFieldDecompositionDM(DM dm, PetscInt dnumber, PetscInt * dsizes, char *** dvarlists, DM * ddm)
496 : {
497 : PetscBool islibmesh;
498 : PetscFunctionBegin;
499 : PetscValidHeaderSpecific(dm,DM_CLASSID,1);
500 6 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH,&islibmesh));
501 6 : if (!islibmesh) LIBMESH_SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Got DM of type %s, not of type %s", ((PetscObject)dm)->type_name, DMLIBMESH);
502 6 : if (dnumber < 0) LIBMESH_SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Negative number %" LIBMESH_PETSCINT_FMT " of decomposition parts", dnumber);
503 : #if PETSC_RELEASE_GREATER_EQUALS(3,20,0)
504 : PetscAssertPointer(ddm,5);
505 : #else
506 : PetscValidPointer(ddm,5);
507 : #endif
508 6 : DM_libMesh * dlm = (DM_libMesh *)(dm->data);
509 6 : LibmeshPetscCallQ(DMCreate(((PetscObject)dm)->comm, ddm));
510 6 : LibmeshPetscCallQ(DMSetType(*ddm, DMLIBMESH));
511 6 : DM_libMesh * ddlm = (DM_libMesh *)((*ddm)->data);
512 6 : ddlm->sys = dlm->sys;
513 : // Each DM owns, and DMDestroy_libMesh() deletes, its own variable and block maps
514 6 : *ddlm->varids = *dlm->varids;
515 6 : *ddlm->varnames = *dlm->varnames;
516 6 : *ddlm->blockids = *dlm->blockids;
517 6 : *ddlm->blocknames = *dlm->blocknames;
518 6 : ddlm->decomposition = new(std::vector<std::set<unsigned int>>);
519 6 : ddlm->decomposition_type = DMLIBMESH_FIELD_DECOMPOSITION;
520 6 : if (dnumber) {
521 0 : for (PetscInt d = 0; d < dnumber; ++d) {
522 0 : if (dsizes[d] < 0) LIBMESH_SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Negative size %" LIBMESH_PETSCINT_FMT " of decomposition part %" LIBMESH_PETSCINT_FMT, dsizes[d],d);
523 0 : ddlm->decomposition->push_back(std::set<unsigned int>());
524 0 : for (PetscInt v = 0; v < dsizes[d]; ++v) {
525 0 : std::string vname(dvarlists[d][v]);
526 0 : std::map<std::string, unsigned int>::const_iterator vit = dlm->varids->find(vname);
527 0 : if (vit == dlm->varids->end())
528 0 : LIBMESH_SETERRQ3(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Variable %" LIBMESH_PETSCINT_FMT " on the %" LIBMESH_PETSCINT_FMT "-th list with name %s is not owned by this DM", v, d, dvarlists[d][v]);
529 0 : unsigned int vid = vit->second;
530 0 : (*ddlm->decomposition)[d].insert(vid);
531 : }
532 : }
533 : }
534 : else { // Empty splits indicate default: split all variables with one per split.
535 2 : PetscInt d = 0;
536 18 : for (const auto & pr : (*ddlm->varids)) {
537 12 : ddlm->decomposition->push_back(std::set<unsigned int>());
538 12 : unsigned int vid = pr.second;
539 12 : (*ddlm->decomposition)[d].insert(vid);
540 12 : ++d;
541 : }
542 : }
543 6 : LibmeshPetscCallQ(DMlibMeshSetUpName_Private(*ddm));
544 6 : LibmeshPetscCallQ(DMSetFromOptions(*ddm));
545 6 : LibmeshPetscCallQ(DMSetUp(*ddm));
546 2 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
547 : }
548 :
549 : #undef __FUNCT__
550 : #define __FUNCT__ "DMlibMeshCreateDomainDecompositionDM"
551 0 : PetscErrorCode DMlibMeshCreateDomainDecompositionDM(DM dm, PetscInt dnumber, PetscInt * dsizes, char *** dblocklists, DM * ddm)
552 : {
553 : PetscBool islibmesh;
554 : PetscFunctionBegin;
555 : PetscValidHeaderSpecific(dm,DM_CLASSID,1);
556 0 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH,&islibmesh));
557 0 : if (!islibmesh) LIBMESH_SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Got DM of type %s, not of type %s", ((PetscObject)dm)->type_name, DMLIBMESH);
558 0 : if (dnumber < 0) LIBMESH_SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Negative number %" LIBMESH_PETSCINT_FMT " of decomposition parts", dnumber);
559 : #if PETSC_RELEASE_GREATER_EQUALS(3,20,0)
560 : PetscAssertPointer(ddm,5);
561 : #else
562 : PetscValidPointer(ddm,5);
563 : #endif
564 0 : DM_libMesh * dlm = (DM_libMesh *)(dm->data);
565 0 : LibmeshPetscCallQ(DMCreate(((PetscObject)dm)->comm, ddm));
566 0 : LibmeshPetscCallQ(DMSetType(*ddm, DMLIBMESH));
567 0 : DM_libMesh * ddlm = (DM_libMesh *)((*ddm)->data);
568 0 : ddlm->sys = dlm->sys;
569 : // Each DM owns, and DMDestroy_libMesh() deletes, its own variable and block maps
570 0 : *ddlm->varids = *dlm->varids;
571 0 : *ddlm->varnames = *dlm->varnames;
572 0 : *ddlm->blockids = *dlm->blockids;
573 0 : *ddlm->blocknames = *dlm->blocknames;
574 0 : ddlm->decomposition = new(std::vector<std::set<unsigned int>>);
575 0 : ddlm->decomposition_type = DMLIBMESH_DOMAIN_DECOMPOSITION;
576 0 : if (dnumber) {
577 0 : for (PetscInt d = 0; d < dnumber; ++d) {
578 0 : if (dsizes[d] < 0) LIBMESH_SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Negative size %" LIBMESH_PETSCINT_FMT " of decomposition part %" LIBMESH_PETSCINT_FMT, dsizes[d],d);
579 0 : ddlm->decomposition->push_back(std::set<unsigned int>());
580 0 : for (PetscInt b = 0; b < dsizes[d]; ++b) {
581 0 : std::string bname(dblocklists[d][b]);
582 0 : std::map<std::string, unsigned int>::const_iterator bit = dlm->blockids->find(bname);
583 0 : if (bit == dlm->blockids->end())
584 0 : LIBMESH_SETERRQ3(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "Block %" LIBMESH_PETSCINT_FMT " on the %" LIBMESH_PETSCINT_FMT "-th list with name %s is not owned by this DM", b, d, dblocklists[d][b]);
585 0 : unsigned int bid = bit->second;
586 0 : (*ddlm->decomposition)[d].insert(bid);
587 : }
588 : }
589 : }
590 : else { // Empty splits indicate default: split all blocks with one per split.
591 0 : PetscInt d = 0;
592 0 : for (const auto & pr : (*ddlm->blockids)) {
593 0 : ddlm->decomposition->push_back(std::set<unsigned int>());
594 0 : unsigned int bid = pr.second;
595 0 : (*ddlm->decomposition)[d].insert(bid);
596 0 : ++d;
597 : }
598 : }
599 0 : LibmeshPetscCallQ(DMlibMeshSetUpName_Private(*ddm));
600 0 : LibmeshPetscCallQ(DMSetFromOptions(*ddm));
601 0 : LibmeshPetscCallQ(DMSetUp(*ddm));
602 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
603 : }
604 :
605 : struct token {
606 : const char * s;
607 : struct token * next;
608 : };
609 :
610 :
611 :
612 : #undef __FUNCT__
613 : #define __FUNCT__ "DMlibMeshFunction"
614 0 : static PetscErrorCode DMlibMeshFunction(DM dm, Vec x, Vec r)
615 : {
616 : PetscFunctionBegin;
617 0 : libmesh_assert(x);
618 0 : libmesh_assert(r);
619 :
620 : NonlinearImplicitSystem * _sys;
621 0 : LibmeshPetscCallQ(DMlibMeshGetSystem(dm, _sys));
622 0 : NonlinearImplicitSystem & sys = *_sys;
623 0 : if (sys.has_constraint_elimination())
624 0 : SETERRQ(((PetscObject)dm)->comm, PETSC_ERR_SUP, "DMlibMesh does not evaluate residuals of systems that eliminate constraints; PetscNonlinearSolver does");
625 0 : PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(sys.solution.get());
626 0 : PetscVector<Number> & R_sys = *cast_ptr<PetscVector<Number> *>(sys.rhs);
627 0 : PetscVector<Number> X_global(x, _sys->comm()), R(r, _sys->comm());
628 :
629 : // Use the systems update() to get a good local version of the parallel solution
630 0 : X_global.swap(X_sys);
631 0 : R.swap(R_sys);
632 :
633 0 : _sys->get_dof_map().enforce_constraints_exactly(*_sys);
634 0 : _sys->update();
635 :
636 : // Swap back
637 0 : X_global.swap(X_sys);
638 0 : R.swap(R_sys);
639 0 : R.zero();
640 :
641 : // if the user has provided both function pointers and objects only the pointer
642 : // will be used, so catch that as an error
643 0 : libmesh_error_msg_if(_sys->nonlinear_solver->residual && _sys->nonlinear_solver->residual_object,
644 : "ERROR: cannot specify both a function and object to compute the Residual!");
645 :
646 0 : libmesh_error_msg_if(_sys->nonlinear_solver->matvec && _sys->nonlinear_solver->residual_and_jacobian_object,
647 : "ERROR: cannot specify both a function and object to compute the combined Residual & Jacobian!");
648 :
649 0 : if (_sys->nonlinear_solver->residual != nullptr)
650 0 : _sys->nonlinear_solver->residual(*(_sys->current_local_solution.get()), R, *_sys);
651 :
652 0 : else if (_sys->nonlinear_solver->residual_object != nullptr)
653 0 : _sys->nonlinear_solver->residual_object->residual(*(_sys->current_local_solution.get()), R, *_sys);
654 :
655 0 : else if (_sys->nonlinear_solver->matvec != nullptr)
656 0 : _sys->nonlinear_solver->matvec(*(_sys->current_local_solution.get()), &R, nullptr, *_sys);
657 :
658 0 : else if (_sys->nonlinear_solver->residual_and_jacobian_object != nullptr)
659 0 : _sys->nonlinear_solver->residual_and_jacobian_object->residual_and_jacobian(*(_sys->current_local_solution.get()), &R, nullptr, *_sys);
660 :
661 : else
662 0 : libmesh_error_msg("Error! Unable to compute residual and/or Jacobian!");
663 :
664 0 : R.close();
665 0 : X_global.close();
666 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
667 0 : }
668 :
669 : #undef __FUNCT__
670 : #define __FUNCT__ "SNESFunction_DMlibMesh"
671 0 : static PetscErrorCode SNESFunction_DMlibMesh(SNES, Vec x, Vec r, void * ctx)
672 : {
673 0 : DM dm = (DM)ctx;
674 : PetscFunctionBegin;
675 0 : LibmeshPetscCallQ(DMlibMeshFunction(dm,x,r));
676 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
677 : }
678 :
679 :
680 : #undef __FUNCT__
681 : #define __FUNCT__ "DMlibMeshJacobian"
682 0 : static PetscErrorCode DMlibMeshJacobian(DM dm, Vec x, Mat jac, Mat pc)
683 : {
684 : PetscFunctionBegin;
685 : NonlinearImplicitSystem * _sys;
686 0 : LibmeshPetscCallQ(DMlibMeshGetSystem(dm, _sys));
687 0 : NonlinearImplicitSystem & sys = *_sys;
688 0 : if (sys.has_constraint_elimination())
689 0 : SETERRQ(((PetscObject)dm)->comm, PETSC_ERR_SUP, "DMlibMesh does not evaluate Jacobians of systems that eliminate constraints; PetscNonlinearSolver does");
690 :
691 0 : libmesh_assert(pc);
692 0 : libmesh_assert(jac);
693 0 : PetscMatrixBase<Number> & the_pc = *PetscMatrixBase<Number>::get_context(pc, sys.comm());
694 0 : PetscMatrixBase<Number> & Jac = *PetscMatrixBase<Number>::get_context(jac, sys.comm());
695 0 : PetscVector<Number> & X_sys = *cast_ptr<PetscVector<Number> *>(sys.solution.get());
696 0 : PetscMatrixBase<Number> & Jac_sys = *cast_ptr<PetscMatrixBase<Number> *>(sys.matrix);
697 0 : PetscVector<Number> X_global(x, sys.comm());
698 :
699 : // Set the dof maps
700 0 : the_pc.attach_dof_map(sys.get_dof_map());
701 0 : Jac.attach_dof_map(sys.get_dof_map());
702 :
703 : // Use the systems update() to get a good local version of the parallel solution
704 0 : X_global.swap(X_sys);
705 0 : Jac.swap(Jac_sys);
706 :
707 0 : sys.get_dof_map().enforce_constraints_exactly(sys);
708 0 : sys.update();
709 :
710 0 : X_global.swap(X_sys);
711 0 : Jac.swap(Jac_sys);
712 :
713 0 : the_pc.zero();
714 :
715 : // if the user has provided both function pointers and objects only the pointer
716 : // will be used, so catch that as an error
717 0 : libmesh_error_msg_if(sys.nonlinear_solver->jacobian && sys.nonlinear_solver->jacobian_object,
718 : "ERROR: cannot specify both a function and object to compute the Jacobian!");
719 :
720 0 : libmesh_error_msg_if(sys.nonlinear_solver->matvec && sys.nonlinear_solver->residual_and_jacobian_object,
721 : "ERROR: cannot specify both a function and object to compute the combined Residual & Jacobian!");
722 :
723 0 : if (sys.nonlinear_solver->jacobian != nullptr)
724 0 : sys.nonlinear_solver->jacobian(*(sys.current_local_solution.get()), the_pc, sys);
725 :
726 0 : else if (sys.nonlinear_solver->jacobian_object != nullptr)
727 0 : sys.nonlinear_solver->jacobian_object->jacobian(*(sys.current_local_solution.get()), the_pc, sys);
728 :
729 0 : else if (sys.nonlinear_solver->matvec != nullptr)
730 0 : sys.nonlinear_solver->matvec(*(sys.current_local_solution.get()), nullptr, &the_pc, sys);
731 :
732 0 : else if (sys.nonlinear_solver->residual_and_jacobian_object != nullptr)
733 0 : sys.nonlinear_solver->residual_and_jacobian_object->residual_and_jacobian(*(sys.current_local_solution.get()), nullptr, &the_pc, sys);
734 :
735 : else
736 0 : libmesh_error_msg("Error! Unable to compute residual and/or Jacobian!");
737 :
738 0 : the_pc.close();
739 0 : Jac.close();
740 0 : X_global.close();
741 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
742 0 : }
743 :
744 : #undef __FUNCT__
745 : #define __FUNCT__ "SNESJacobian_DMlibMesh"
746 0 : static PetscErrorCode SNESJacobian_DMlibMesh(SNES, Vec x, Mat jac, Mat pc, void * ctx)
747 : {
748 0 : DM dm = (DM)ctx;
749 : PetscFunctionBegin;
750 0 : LibmeshPetscCallQ(DMlibMeshJacobian(dm,x,jac,pc));
751 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
752 : }
753 :
754 : #undef __FUNCT__
755 : #define __FUNCT__ "DMVariableBounds_libMesh"
756 1224 : static PetscErrorCode DMVariableBounds_libMesh(DM dm, Vec xl, Vec xu)
757 : {
758 : NonlinearImplicitSystem * _sys;
759 1224 : LibmeshPetscCallQ(DMlibMeshGetSystem(dm, _sys));
760 1224 : NonlinearImplicitSystem & sys = *_sys;
761 : PetscFunctionBegin;
762 : // Workaround for nonstandard Q suffix warning with quad precision
763 408 : const PetscReal petsc_inf = std::numeric_limits<PetscReal>::max() / 4;
764 :
765 : // Users compute bounds on full vectors. When eliminating constraints, SNES's xl and xu are
766 : // reduced, and they take the bounds of the unconstrained degrees of freedom
767 1224 : const bool reduced = sys.has_constraint_elimination();
768 816 : std::unique_ptr<NumericVector<Number>> full_l, full_u;
769 1224 : if (reduced)
770 : {
771 8 : full_l = sys.solution->zero_clone();
772 8 : full_u = sys.solution->zero_clone();
773 : }
774 2040 : PetscVector<Number> XL(xl, sys.comm()), XU(xu, sys.comm());
775 1224 : NumericVector<Number> & L = reduced ? *full_l : static_cast<NumericVector<Number> &>(XL);
776 816 : NumericVector<Number> & U = reduced ? *full_u : static_cast<NumericVector<Number> &>(XU);
777 1224 : L.add(-petsc_inf);
778 1224 : U.add(petsc_inf);
779 1224 : L.close();
780 1224 : U.close();
781 :
782 1224 : if (sys.nonlinear_solver->bounds != nullptr)
783 0 : sys.nonlinear_solver->bounds(L,U,sys);
784 1224 : else if (sys.nonlinear_solver->bounds_object != nullptr)
785 1224 : sys.nonlinear_solver->bounds_object->bounds(L,U, sys);
786 : else
787 0 : SETERRQ(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "No bounds calculation in this libMesh object");
788 :
789 1224 : if (reduced)
790 : {
791 2 : auto & solver = cast_ref<PetscNonlinearSolver<Number> &>(*sys.nonlinear_solver);
792 6 : L.close();
793 6 : U.close();
794 6 : solver.restrict_vector(L, xl);
795 6 : solver.restrict_vector(U, xu);
796 : }
797 :
798 408 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
799 408 : }
800 :
801 :
802 : #undef __FUNCT__
803 : #define __FUNCT__ "DMCreateGlobalVector_libMesh"
804 1224 : static PetscErrorCode DMCreateGlobalVector_libMesh(DM dm, Vec *x)
805 : {
806 : PetscFunctionBegin;
807 1224 : DM_libMesh * dlm = (DM_libMesh *)(dm->data);
808 : PetscBool eq;
809 :
810 1224 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH, &eq));
811 :
812 1224 : if (!eq)
813 0 : LIBMESH_SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "DM of type %s, not of type %s", ((PetscObject)dm)->type_name, DMLIBMESH);
814 :
815 1224 : if (!dlm->sys)
816 0 : SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "No libMesh system set for DM_libMesh");
817 :
818 408 : NumericVector<Number> * nv = (dlm->sys->solution).get();
819 1224 : PetscVector<Number> * pv = dynamic_cast<PetscVector<Number> *>(nv);
820 816 : Vec v = pv->vec();
821 408 : const DofMap & dof_map = dlm->sys->get_dof_map();
822 : /* Unfortunately, currently this does not produce a ghosted vector, so nonlinear subproblem solves aren't going to be easily available.
823 : Should work fine for getting vectors out for linear subproblem solvers. */
824 1224 : if (dlm->embedding || dof_map.has_constraint_elimination()) {
825 : /* The size of an embedded subproblem, or of the reduced space SNES solves on */
826 : PetscInt n;
827 6 : LibmeshPetscCallQ(VecCreate(((PetscObject)v)->comm, x));
828 6 : if (dlm->embedding)
829 0 : LibmeshPetscCallQ(ISGetLocalSize(dlm->embedding, &n));
830 : else
831 6 : n = cast_int<PetscInt>(dof_map.get_constraint_elimination().n_local_dofs());
832 6 : LibmeshPetscCallQ(VecSetSizes(*x,n,PETSC_DETERMINE));
833 6 : LibmeshPetscCallQ(VecSetType(*x,((PetscObject)v)->type_name));
834 6 : LibmeshPetscCallQ(VecSetFromOptions(*x));
835 6 : LibmeshPetscCallQ(VecSetUp(*x));
836 : }
837 : else {
838 1218 : LibmeshPetscCallQ(VecDuplicate(v,x));
839 : }
840 :
841 : #if PETSC_VERSION_LESS_THAN(3,13,0)
842 816 : LibmeshPetscCallQ(PetscObjectCompose((PetscObject)*x,"DM",(PetscObject)dm));
843 : #else
844 408 : LibmeshPetscCallQ(VecSetDM(*x, dm));
845 : #endif
846 :
847 408 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
848 : }
849 :
850 :
851 :
852 :
853 : #undef __FUNCT__
854 : #define __FUNCT__ "DMCreateMatrix_libMesh"
855 0 : static PetscErrorCode DMCreateMatrix_libMesh(DM dm, Mat * A)
856 : {
857 : PetscFunctionBegin;
858 0 : DM_libMesh * dlm = (DM_libMesh *)(dm->data);
859 : PetscBool eq;
860 :
861 0 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH, &eq));
862 :
863 0 : if (!eq)
864 0 : LIBMESH_SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "DM of type %s, not of type %s", ((PetscObject)dm)->type_name, DMLIBMESH);
865 :
866 0 : if (!dlm->sys)
867 0 : SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "No libMesh system set for DM_libMesh");
868 :
869 0 : *A = (dynamic_cast<PetscMatrixBase<Number> *>(dlm->sys->matrix))->mat();
870 0 : LibmeshPetscCallQ(PetscObjectReference((PetscObject)(*A)));
871 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
872 : }
873 :
874 :
875 : #undef __FUNCT__
876 : #define __FUNCT__ "DMView_libMesh"
877 0 : static PetscErrorCode DMView_libMesh(DM dm, PetscViewer viewer)
878 : {
879 : PetscBool isascii;
880 : const char * name, * prefix;
881 0 : DM_libMesh * dlm = (DM_libMesh *)dm->data;
882 : PetscFunctionBegin;
883 0 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
884 0 : if (isascii) {
885 0 : LibmeshPetscCallQ(PetscObjectGetName((PetscObject)dm, &name));
886 0 : LibmeshPetscCallQ(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
887 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "DM libMesh with name %s and prefix %s\n", name, prefix));
888 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "blocks:"));
889 0 : std::map<std::string,unsigned int>::iterator bit = dlm->blockids->begin();
890 0 : std::map<std::string,unsigned int>::const_iterator bend = dlm->blockids->end();
891 0 : for (; bit != bend; ++bit) {
892 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "(%s,%d) ", bit->first.c_str(), bit->second));
893 : }
894 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "\n"));
895 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "variables:"));
896 0 : std::map<std::string,unsigned int>::iterator vit = dlm->varids->begin();
897 0 : std::map<std::string,unsigned int>::const_iterator vend = dlm->varids->end();
898 0 : for (; vit != vend; ++vit) {
899 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "(%s,%d) ", vit->first.c_str(), vit->second));
900 : }
901 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "\n"));
902 0 : if (dlm->decomposition_type == DMLIBMESH_NO_DECOMPOSITION) {
903 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "No decomposition\n"));
904 : }
905 : else {
906 0 : if (dlm->decomposition_type == DMLIBMESH_FIELD_DECOMPOSITION) {
907 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "Field decomposition by variable: "));
908 : }
909 0 : else if (dlm->decomposition_type == DMLIBMESH_DOMAIN_DECOMPOSITION) {
910 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "Domain decomposition by block: "));
911 : }
912 0 : else LIBMESH_SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_PLIB, "Unexpected decomposition type: %d", dlm->decomposition_type);
913 : /* FIX: decompositions might have different sizes and components on different ranks. */
914 0 : for (auto d : index_range(*dlm->decomposition)) {
915 0 : std::set<unsigned int>::iterator dbegin = (*dlm->decomposition)[d].begin();
916 0 : std::set<unsigned int>::iterator dit = (*dlm->decomposition)[d].begin();
917 0 : std::set<unsigned int>::iterator dend = (*dlm->decomposition)[d].end();
918 0 : for (; dit != dend; ++dit) {
919 0 : if (dit != dbegin) {
920 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, ","));
921 : }
922 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "%u", *dit));
923 : }
924 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, ";"));
925 : }
926 0 : LibmeshPetscCallQ(PetscViewerASCIIPrintf(viewer, "\n"));
927 : }
928 : }
929 :
930 0 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
931 : }
932 :
933 : #undef __FUNCT__
934 : #define __FUNCT__ "DMSetUp_libMesh"
935 2580 : static PetscErrorCode DMSetUp_libMesh(DM dm)
936 : {
937 : PetscFunctionBegin;
938 2580 : DM_libMesh * dlm = (DM_libMesh *)(dm->data);
939 : PetscBool eq;
940 :
941 2580 : LibmeshPetscCallQ(PetscObjectTypeCompare((PetscObject)dm, DMLIBMESH, &eq));
942 :
943 2580 : if (!eq)
944 0 : LIBMESH_SETERRQ2(((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "DM of type %s, not of type %s", ((PetscObject)dm)->type_name, DMLIBMESH);
945 :
946 2580 : if (!dlm->sys)
947 0 : SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "No libMesh system set for DM_libMesh");
948 : /*
949 : Do not evaluate function, Jacobian or bounds for an embedded DM -- the subproblem might not have enough information for that.
950 : */
951 2580 : if (!dlm->embedding) {
952 2580 : LibmeshPetscCallQ(DMSNESSetFunction(dm, SNESFunction_DMlibMesh, (void *)dm));
953 2580 : LibmeshPetscCallQ(DMSNESSetJacobian(dm, SNESJacobian_DMlibMesh, (void *)dm));
954 3440 : if (dlm->sys->nonlinear_solver->bounds || dlm->sys->nonlinear_solver->bounds_object)
955 2436 : LibmeshPetscCallQ(DMSetVariableBounds(dm, DMVariableBounds_libMesh));
956 : }
957 : else {
958 : /*
959 : Fow now we don't implement even these, although a linear "Dirichlet" subproblem is well-defined.
960 : Creating the submatrix, however, might require extracting the submatrix preallocation from an unassembled matrix.
961 : */
962 0 : dm->ops->createglobalvector = 0;
963 0 : dm->ops->creatematrix = 0;
964 : }
965 860 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
966 : }
967 :
968 :
969 :
970 :
971 : #undef __FUNCT__
972 : #define __FUNCT__ "DMDestroy_libMesh"
973 2580 : static PetscErrorCode DMDestroy_libMesh(DM dm)
974 : {
975 2580 : DM_libMesh * dlm = (DM_libMesh *)(dm->data);
976 : PetscFunctionBegin;
977 4300 : delete dlm->varids;
978 4300 : delete dlm->varnames;
979 4300 : delete dlm->blockids;
980 4300 : delete dlm->blocknames;
981 2580 : delete dlm->decomposition;
982 2580 : LibmeshPetscCallQ(ISDestroy(&dlm->embedding));
983 2580 : LibmeshPetscCallQ(PetscFree(dm->data));
984 :
985 2580 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
986 : }
987 :
988 : EXTERN_C_BEGIN
989 : #undef __FUNCT__
990 : #define __FUNCT__ "DMCreate_libMesh"
991 2580 : PetscErrorCode DMCreate_libMesh(DM dm)
992 : {
993 : DM_libMesh * dlm;
994 :
995 : PetscFunctionBegin;
996 : PetscValidHeaderSpecific(dm,DM_CLASSID,1);
997 2580 : LibmeshPetscCallQ(PetscNew(&dlm));
998 2580 : dm->data = dlm;
999 :
1000 : /* DMlibMesh impl */
1001 2580 : dlm->varids = new(std::map<std::string, unsigned int>);
1002 2580 : dlm->blockids = new(std::map<std::string, unsigned int>);
1003 2580 : dlm->varnames = new(std::map<unsigned int, std::string>);
1004 2580 : dlm->blocknames = new(std::map<unsigned int, std::string>);
1005 2580 : dlm->decomposition = LIBMESH_PETSC_NULLPTR;
1006 2580 : dlm->decomposition_type = DMLIBMESH_NO_DECOMPOSITION;
1007 :
1008 : /* DM API */
1009 2580 : dm->ops->createglobalvector = DMCreateGlobalVector_libMesh;
1010 2580 : dm->ops->createlocalvector = 0; // DMCreateLocalVector_libMesh;
1011 2580 : dm->ops->getcoloring = 0; // DMGetColoring_libMesh;
1012 2580 : dm->ops->creatematrix = DMCreateMatrix_libMesh;
1013 2580 : dm->ops->createinterpolation= 0; // DMCreateInterpolation_libMesh;
1014 :
1015 2580 : dm->ops->refine = 0; // DMRefine_libMesh;
1016 2580 : dm->ops->coarsen = 0; // DMCoarsen_libMesh;
1017 :
1018 : // * dm->ops->getinjection was renamed to dm->ops->createinjection in PETSc 5a84ad338 (5 Jul 2019)
1019 : // * dm->ops-getaggregates was removed in PETSc 97779f9a (5 Jul 2019)
1020 : // * Both changes were merged into PETSc master in 94aad3ce (7 Jul 2019).
1021 : #if PETSC_VERSION_LESS_THAN(3,12,0)
1022 1720 : dm->ops->getinjection = 0; // DMGetInjection_libMesh;
1023 1720 : dm->ops->getaggregates = 0; // DMGetAggregates_libMesh;
1024 : #else
1025 860 : dm->ops->createinjection = 0;
1026 : #endif
1027 :
1028 :
1029 2580 : dm->ops->createfielddecomposition = DMCreateFieldDecomposition_libMesh;
1030 2580 : dm->ops->createdomaindecomposition = DMCreateDomainDecomposition_libMesh;
1031 :
1032 2580 : dm->ops->destroy = DMDestroy_libMesh;
1033 2580 : dm->ops->view = DMView_libMesh;
1034 2580 : dm->ops->setfromoptions = 0; // DMSetFromOptions_libMesh;
1035 2580 : dm->ops->setup = DMSetUp_libMesh;
1036 :
1037 : /* DMlibMesh API */
1038 2580 : LibmeshPetscCallQ(PetscObjectComposeFunction((PetscObject)dm,"DMlibMeshSetSystem_C",DMlibMeshSetSystem_libMesh));
1039 2580 : LibmeshPetscCallQ(PetscObjectComposeFunction((PetscObject)dm,"DMlibMeshGetSystem_C",DMlibMeshGetSystem_libMesh));
1040 :
1041 2580 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
1042 : }
1043 : EXTERN_C_END
1044 :
1045 : #endif // LIBMESH_HAVE_PETSC
|
