-
Notifications
You must be signed in to change notification settings - Fork 285
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Static condensation #3883
base: devel
Are you sure you want to change the base?
Static condensation #3883
Conversation
bdb4d4e
to
21e8140
Compare
aa8246a
to
1468055
Compare
Use the same mutex for both `operator<<` overloads. Refs CIVET failure https://civet.inl.gov/job/2295949/ on libMesh/libmesh#3883
3a98105
to
11f7769
Compare
Job Test mac on 11f7769 : invalidated by @lindsayad |
Job Coverage on 31388ef wanted to post the following: Coverage
Warnings
This comment will be updated on new commits. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
There may be other issues I've missed; this is a pretty massive PR, and I'm running out of pre-vacation time, and I've already hit you with a ton of complaints. ;-)
There are a few obvious utility things we could split into smaller PRs if you wanted to get those in ASAP.
* Fills the vector \p di with the global degree of freedom indices | ||
* for the element. For one variable, and potentially for a | ||
* non-default element p refinement level | ||
*/ |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
We've got to have some documentation here on how those functors work.
I think what I will do is split out the utilities into their own PRs. And then I'm going to add more work that's coming from using this in the HDG example. I'm pretty sure in that example I will not be allowed to remove the pressure dofs (even though they are all "internal") from the global system, so some more work will still need to be done to handle this in the library. Given that, I'm returning this to draft for now |
this will go along with @roystgnr's desire to use "condensed"/"uncondensed" verbiage over internal/trace |
9491d1f
to
6818bf3
Compare
f4f4afe
to
da85e83
Compare
32b4e29
to
e60b60c
Compare
Job Min gcc on e60b60c : invalidated by @lindsayad |
Job Test 32bit on e60b60c : invalidated by @lindsayad |
Job Test debug on e60b60c : invalidated by @lindsayad |
finally ready for review again 😄 |
include/base/dof_map.h
Outdated
@@ -1616,7 +1634,7 @@ class DofMap : public ReferenceCountedObject<DofMap>, | |||
* constraints or sufficiently many user constraints. | |||
*/ | |||
std::unique_ptr<SparsityPattern::Build> build_sparsity(const MeshBase & mesh, | |||
bool calculate_constrained = false) const; | |||
const bool calculate_constrained = false) const; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
At some point we've got to make a global decision about when pass-by-value parameters should be const
. On the one hand it's helpful for protection the implementation, but on the other hand it just muddies up the declaration with an implementation-specific detail that doesn't actually promise anything.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yeah, since we don't pass const
non-reference parameters like this consistently elsewhere in the library, I agree it looks a bit out of place to make this change. I could be talked into inconsistently having const
parameters in the .C file based on the author's preference, but I'm not convinced that really helps catch lots of bugs either.
/// If there are "spider" nodes in the mesh (i.e. a single node which | ||
/// is connected to many 1D elements) and Constraints, we can end up | ||
/// sorting the same set of DOFs multiple times in handle_vi_vj(), | ||
/// only to find that it has no net effect on the final sparsity. In | ||
/// such cases it is much faster to keep track of (element_dofs_i, | ||
/// element_dofs_j) pairs which have already been handled and not | ||
/// repeat the computation. We use this data structure to keep track | ||
/// of hashes of sets of dofs we have already seen, thus avoiding | ||
/// unnecessary calculations. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Surely this isn't necessary?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I didn't realize that Doxygen ignored C++ comments. Could we just switch to the ASCII-art-C-style that we use e.g. for clear_full_sparsity()
above?
include/numerics/petsc_matrix.h
Outdated
@@ -83,7 +83,7 @@ class PetscMatrix final : public PetscMatrixBase<T> | |||
*/ | |||
explicit | |||
PetscMatrix (Mat m, | |||
const Parallel::Communicator & comm_in); | |||
const Parallel::Communicator & comm_in); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Clean this up
static void uncondensed_dofs_from_scalar_dofs(std::vector<dof_id_type> & uncondensed_dofs, | ||
const std::vector<dof_id_type> & scalar_dofs); | ||
|
||
static void total_dofs_from_field_dof(std::vector<dof_id_type> & dofs, |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yeah, they're private, but they all need documentation. Total dofs per what? What is a field dof?
* have the same sparsity pattern as the matrix used during | ||
* solution. When not \p System but the user wants to | ||
* initialize the mayor matrix, then all the additional matrices, | ||
* if existent, have to be initialized by the user, too. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
What's the value of adding a matrix to the system that the system isn't going to be maintaining it during reinit()
etc? Just throwing the memory management hot potato to System
to catch? That doesn't seem worth inconsistent behavior (an uninitialized matrix after reinit()
if we skip it) or broken behavior (we aren't skipping anything in _matrices
in System::reinit()
, right? what happens when we compute_sparsity()
for one you didn't want us to mess with?).
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I just copy-pasted the doxygen from the method above. So I don't know what to do here. It seems like you signed off on this back when it was added in #2638
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The "mayor" ("main") in the comment looks like it was your own version of the brown M&Ms trick, and I missed it. But I didn't read the comment carefully enough this time either, just freaked out at imagining trying to deal with a mix of System-initialized and user/subclass-initialized sparsity patterns.
So I guess my remaining concern is just - if we have all user/subclass-initialized matrices, aren't we still doing a compute_sparsity()
anyway, and are we sure that's never going to be redundant?
src/numerics/petsc_matrix.C
Outdated
const numeric_index_type n_in, | ||
const numeric_index_type m_l, | ||
const numeric_index_type n_l, | ||
const std::vector<numeric_index_type> & n_nz, | ||
const std::vector<numeric_index_type> & n_oz, | ||
const numeric_index_type blocksize_in) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Okay, my new policy is that if I see brown M&Ms backstage I'm not putting equipment or people at risk onstage.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I can understand the active objection to labeling const
in the headers, but I don't understand the active complaint about marking variables that don't change as const
in the implementation. If it's an objection for the method parameters, then it seems like it should be a general objection
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Does Github show you something different than it shows me here? I didn't see any change having anything to do with const
, just a bunch of added whitespace that broke the formatting in the file.
@@ -257,6 +258,16 @@ class PetscLinearSolver : public LinearSolver<T> | |||
*/ | |||
virtual LinearConvergenceReason get_converged_reason() const override; | |||
|
|||
protected: | |||
#if defined(LIBMESH_ENABLE_AMR) && defined(LIBMESH_HAVE_METAPHYSICL) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Do we need these guards? I thought the wrapper was good for field splits even on an unrefined mesh, and I know MetaPhysicL is only for the refinement+GMG case.
include/systems/implicit_system.h
Outdated
@@ -335,13 +336,37 @@ class ImplicitSystem : public ExplicitSystem | |||
*/ | |||
mutable std::unique_ptr<LinearSolver<Number>> linear_solver; | |||
|
|||
/** | |||
* Retrieve the static condensation class. Errors if the user has not | |||
* requested it on the command line, so calls to this should be guarded by |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Surely we should have both a CLI and a C++ setter interface to that?
@@ -398,7 +398,7 @@ extern "C" | |||
sys.update(); | |||
X.swap(X_sys); | |||
|
|||
// Let's also wrap J and Jpre in PetscMatrix objects for convenience | |||
// Let's also wrap J and Jpre in PetscMatrixBase objects for convenience |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Not sure if I missed these last time or if you got overzealous search-and-replacing comments to fix a converse complaint of mine last time?
{ | ||
if (libMesh::on_command_line("--" + name_in + "-static-condensation")) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Can we move this block to a separate method so static condensation can be enabled from code? Or would it have to be a constructor argument? It would seem to me that turning SC on any time before System::init()
ought to be feasible, at least.
Use the same mutex for both `operator<<` overloads. Refs CIVET failure https://civet.inl.gov/job/2295949/ on libMesh/libmesh#3883
I added this at a time when I thought I would add a test with a field split of a statically condensed system, but I never did. So the need is no longer there and Roy hates it so why not remove
0b22214
to
31388ef
Compare
MOOSE failure is unrelated and already codified in idaholab/moose#27555 |
This should be ready for another round of review |
Adds the ability to condense out element internal degrees of freedom. This can be used to eliminate internal degrees of freedom in high order CFEM or to eliminate internal degrees of freedom in HDG methods