Automatically force full matrix when using AD with global indexing #16395
Labels
C: Automatic Differentiation
Tickets pertaining the MetaPhysicL based forward mode AD system
C: Framework
P: normal
A defect affecting operation with a low possibility of significantly affects.
T: task
An enhancement to the software.
Reason
AD with global indexing adds to the Jacobian for any non-zero dependence. It would be difficult to prevent additions to the matrix if the user has not specified coupling between variables which exists in practice; we would have to have query two maps: a map from the i index to variable number, and a map from the j index to the variable number. This could induce undesirable expense in most cases. Consequently, we will allow AD with global indexing to add any non-zero coupling that exists. With that in mind we want to, by default, assume a full coupling sparsity pattern whenever we have AD objects and we are using global indexing. By assuming a full coupling sparsity pattern, we will prevent new nonzero allocations in PETSc that can greatly slow down simulations.
Design
Automatic triggering of full coupling has to happen early in the simulation, ideally before
EquationSystems::init
or else we will have to reinit theEquationSystems
later, which is typically a fairly expensive operation. Unfortunately AD objects are typically not added until well afterEquationSystems::init
. To circumvent this, I will leverage ourRelationshipManager/GhostingFunctor
system sinceGhostingFunctors
are designed to be used to as coupling functors for determination of the sparsity pattern/coupling.Additionally, I think it makes sense to allow the user to turn off the full matrix if they trust the sparsity pattern they're providing in their input file. This could be very important if there are a lot of variables, and the coupling between these variables is sparse. In that case, preallocating a fully sparsity pattern would cause an unnecessary initial memory spike (although later assembly of the matrix would shrink out the unnecessary memory).
Impact
In NavierStokes for example, we have erroring on new nonzero allocations turned on by default. This resulted in @smharper encountering a PETSc message that he didn't understand (see #15644 (comment) and #15644 (comment)); almost all of our users would be equally if not more confused. We want to avoid that. So just like we did in #13411 we'll automatically construct a full matrix (by default) behind the scenes when using AD with global indexing whether we are doing NEWTON or PJFNK.
The text was updated successfully, but these errors were encountered: