search

brendanjmeade / celeri

Next generation earthquake cycle kinematics
BSD 3-Clause "New" or "Revised" License
25 stars 7 forks source link

Fixed four errors! Now, the sparse matvec is working correctly. #90

Closed tbenthompson closed 2 years ago

tbenthompson commented 2 years ago
  1. H-matrix construction was failing in rare situations where an approximate block turned out to have full rank.
  2. The TDE smoothing and TDE slip rate constraint rows were wrong. This was affecting non-sparse solving too (the smoothing was wrong!!).
  3. Use smoothing_matrix_simple instead of smoothing_matrix?
  4. Use the negative of the tde velocity matrix?

For these last two issues, I don't understand the underlying reasons, but I'm just copying the behavior in assemble_and_solve_dense into the sparse matvec and H-matrix construction functions.

See the "#TBTERROR" comments in the notebook. I am a terror.

review-notebook-app[bot] commented 2 years ago

Check out this pull request on  ReviewNB

See visual diffs & provide feedback on Jupyter Notebooks.

Powered by ReviewNB

jploveless commented 2 years ago

This is strange: I thought I had fixed the indexing for TDEs last week, as well as the sign of the tde velocity matrix (it should indeed be negative). But now I see that I was thinking about the assembly incorrectly: we are actually doing smoothing then constraints as a pair for each mesh, rather than smoothing for all meshes, then constraints for all meshes. This could explain the fact that the edge constraints seemed to be ignored in the dense solve! Thanks Ben!

brendanjmeade commented 2 years ago

@jploveless I had thought that we were doing all of the smoothing matrices and then all of the TDE constraints after that as well. As long as we pick a standard and keep with it, we should be good. This is also a good reminder to me to create a notebook that shows the structure of the linear system for the multiple meshes, case so that we have that for reference and documentaiton. See: https://github.com/brendanjmeade/celeri/issues/97

jploveless commented 2 years ago

I think either method can make sense conceptually, but I do feel that Ben's new indexing (smoothing, constraint pairs) is more consistent overall with how meshes are treated within celeri: they are referenced separately, unlike in Blocks where they are treated as part of a single structure.

tbenthompson commented 2 years ago

Just wanted to weigh in and say that I like Jack's argument for (smoothing1, constraint1, smoothing2, constraint2)! But, I also didn't realize that I was modifying the intended indexing when I changed that. The previous implementation was just wrong since it was overlapping the indices and everything was mushy bananas.

brendanjmeade commented 2 years ago

Added what I think is clean LaTex description of how we structure this in the notebook reference_linear_operator_structure.ipynb with https://github.com/brendanjmeade/celeri/commit/b12e6769130abded0430597c831753e3675446c0.

I haven't looked at the code for the indexingere. @tbenthompson and @jploveless are you in agreement on this structure?

jploveless commented 2 years ago

I know this isn't helpful, but I really could go either way!

brendanjmeade commented 2 years ago

@jploveless

  • The new LaTeX description looks fantastic and clearly lays out the structure. However, it would be cleaner in terms of LaTeX matrix structure to do all smoothing, then all slip rate constraints (no need for m_1 through m_n references)

I hear you that it would look cleaner but I think it shows the current assembly accurately?

  • The construction of indices via for loop through the meshes (as is currently done) is more consistent with using smoothing, constraint pairs for each mesh sequentially.

Are you thinking that it's worthwhile to change this to all of the smoothing and then all of the constraints?

jploveless commented 2 years ago

Sorry — I was trying to point out some pros and cons in my previous bullet points. Right now, Ben's indexing works well, and I think that conceptually it is more consistent with our more individualized treatment of meshes. Your LaTeX summary is accurate as is; I was just thinking that writing the assembly would be more compact if we did all smoothing, then all constraints. But I don't think it's necessary to change to that scheme, so I say we leave things as they are now.

brendanjmeade commented 2 years ago

I understand now! Let's leave it the way it is. Thanks!