Condensed representation for Dirichlet BCs

[jtostie]

We would like to be able to construct the linear system without having to insert values of 1.0 in the system associated with the essential BCs.

[ikalash]

I believe @jfike did this for MOR a few weeks ago -- see commit b18b881cca750c2f237b46e1b0bb522d1054215f . It may be useful to look at this for this task.

[lxmota]

This may also be related to (#53). Still investigating.

[bartgol]

Is this desire related to loss of symmetry in the matrix? Would it help to have a feature in Albany that maintains symmetry of the matrix (diagonalizing the columns of Dirichlet dofs as well)?

[ibaned]

Although I'm not an expert on this, I think another concern may be matrix condition number. These 1.0 values may not be in the same order of magnitude as the other entries which are influenced by physical constants. It would be nice to remove this as a source of concerns when debugging convergence issues.

[ikalash]

There are two issues for LCM:

1.) The 1s on the diagonal can ruin the condition number. I did a study of this a few months ago and the condition number can be O(10^11) with the 1s on the diagonal and O(10) w/ the 1s removed (!!!!). 2.) The 1s on the diagonal are ruining convergence rates for one of the implementations of the Schwarz coupling in Albany, leading to slow, non-monotonic convergence of the Newton-Schwarz loop.

The loss of symmetry would be a third reason to change the BCs.

[ibaned]

@ikalash just curious, how did you compute/estimate the condition number ? I'd really like to be able to do that in the future.

[ikalash]

I dumped the Jacobians to matrix market file (there is an option you can add to your input file to do that) and computed it in MATLAB -- so not a very elegant approach. It may be good to add the (optional) computation of condition number in the code as a debug option. It should be easy to do if Tpetra/Epetra have a routine for the condition number. I could add this option, if it is of interest.

Of course, you would not be able to look at the relative condition number with the 1s removed from the diagonal, until that implementation exists in Albany.

[maxrpi]

At RPI, we would find Albany having an option of reporting condition numbers to be of enormous interest.

On 2/6/2017 11:57 AM, Irina K. Tezaur wrote:

I dumped the Jacobians to matrix market file (there is an option you can add to your input file to do that) and computed it in MATLAB -- so not a very elegant approach. It may be good to add the (optional) computation of condition number in the code as a debug option. It should be easy to do if Tpetra/Epetra have a routine for the condition number. I could add this option, if it is of interest.

Of course, you would not be able to look at the relative condition number with the 1s removed from the diagonal, until that implementation exists in Albany.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/gahansen/Albany/issues/13#issuecomment-277743309, or mute the thread https://github.com/notifications/unsubscribe-auth/AGw83OzC82jObvkuB5X_A4eR7a5zZEB2ks5rZ1D2gaJpZM4JoA-P.

[lxmota]

This is of great interest to us for at least three reasons:

• Improve the condition number of the systems to solve.
• It appears that the way Dirichlet BCs are implemented now affects negatively the convergence of the Schwarz method.
• Application of Dirichlet BCs for explicit dynamics.

[csiefer2]

As a note, the condition number of the diagonally scaled matrix should be what actually matters w.r.t. solving, not the condition number of the matrix proper.

Also, AztecOO / Belos have condition number estimators built into variants of their CG routines that can be used.

---

From: Alejandro Mota notifications@github.com Sent: Monday, February 6, 2017 10:45 AM To: gahansen/Albany Cc: Subscribed Subject: [EXTERNAL] Re: [gahansen/Albany] Condensed representation for Dirichlet BCs (#13)

This is of great interest to us for at least three reasons:

• Improve the condition number of the systems to solve.
• It appears that the way Dirichlet BCs are implemented now affects negatively the convergence of the Schwarz method.
• Application of Dirichlet BCs for explicit dynamics.

- You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHubhttps://github.com/gahansen/Albany/issues/13#issuecomment-277757202, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AH63V8dGrSn4zWNhGwhCUQXZFcIBB93fks5rZ1w6gaJpZM4JoA-P.

[bartgol]

I may be over-optimistic, but I would argue that the condition number would not deteriorate if one also diagonalizes the column.

Also, a couple of remarks:

1) we can (easily?) add the possibility of choosing a different scaling coefficient, perhaps user defined. If not provided, we could use the average of the row elements instead of 1.0. That should keep the condition number close to the original one. 2) condition number for GMRES (or Krylov methods in general) does not necessarily mean bad convergence. A nxn system with only two distinct e-values={1e-5, 1e5} has condition number 1e10, but would converge (in exact arithmetic) in 2 iterations. What actually matters is the clustering of the eigenvalues. So unless the diagonalization causes the e-values to spread out more, GMRES should have similar convergence properties. Which is why I believe that a full diagonalization (also the column) of the matrix should not deteriorate convergence.

[ikalash]

In response to @bartgol and @csiefer2 's comments: I just checked in some input files that demonstrate how to specify different scalings in Albany. I put in the option to do various scalings to try to improve Jacobian condition numbers and hence performance into Albany some months ago. I think this feature is largely unknown outside the LCM group, so I wanted to bring it to people's attention.

The scaling feature can be combined with the new in-code condition number estimation utility (see commit https://github.com/gahansen/Albany/commit/9b2a4b8ce3a05672efb156bd74be437b4ad409e5) to study the condition number as a function of the scaling. Here are some results for the 3 available scalings vs. no scaling for an LCM test case, the Rubiks Cube problem: . You can see the condition number goes from O(10^5) with no scaling to O(10)-O(100) with scaling. Diagonal scaling gives the lowest condition number, but, somewhat oddly, it leads to the most Newton iterations.

I encourage for those interested to play with the scalings to see if it improves the issues caused by the 1s on the diagonal in the Jacobians in the current BC implementation, or by differences in scales for, e.g., a multiphysics problem, for their application. I will look at this for Schwarz coupling once we resolve a couple other recently-found bugs in the Schwarz implementation (see Issue SNLComputation/Albany#67).

[djlittl]

Thanks for this, Irina.

Does the scaling affect the value of the residual? I am wondering if scaling effectively tightens the convergence tolerance in cases when a relative tolerance is used. This could perhaps explain the increase in the number of Newton iterations you're seeing with diagonal scaling.

[lxmota]

Created new issue (#68) to focus the discussion here about Dirichlet boundary conditions.

[bgranzow]

I wrote a short set of notes related to this issue that may be of some use.

[lxmota]

@bgranzow Your notes are clear and excellent. This is exactly what we want to do.

I agree that zeroing out the row and column of the Jacobian, keeping the diagonal entry intact, and using a zero residual entry for that row would achieve the imposition of the Dirichlet BC. All this without introducing the horrible scaling problems that are hampering our progress in so many fronts, despite endless reassurances that preconditioners would take care of this.

I also agree that this seems to be much easier than to try to form the reduced system; and that perhaps the best way to achieve this is to use NOX PrePost operators. I'm using those operators right now for Schwarz coupling. Within the operator you have access to the solution vector. There are methods called getX() and setX() in the solution group class that would be useful for this, and as you mention this would happen only once before NOX solves the system.

This issue has now become of the highest priority for the LCM project. @ikalash and I have seen that even implicit dynamics has severe problems when using the current implementation of the BCs. We are supposed to deliver a dynamic Schwarz capability by the end of the FY. Without resolving this issue, we are stuck.

Bringing @jwfoulk and @calleman21 into the discussion so that they are aware of developments.

[jtostie]

I agree that Brian's notes are good, and a reasonable path forward. It seems like we should be able to write a Dirichlet BC that modifies the residual and solution vector accordingly. Can we do that @lxmota?

I still strongly prefer creating an alternate map that only considers the active DOFs and only assembling the system of active DOFs.

[lxmota]

@jtostie Yes, my strong preference would be to assemble the reduced system from the start.

But apparently this would require too much infrastructure change.

@bgranzow's proposed method can be implemented in two parts: A new set of DBCs that would set their corresponding residual entries to zero, and their Jacobian to zero out their rows and columns, except for the diagonal entries.

The other part is setting the solution to the desired value. Because the ModelEvaluator does not allow the modifications of the solution, the proposal is to do it by means of a NOX PrePost operator, where the solution vector becomes available for modification. This is jumping through hoops, but this is a constraint imposed by the ModelEvaluator.

Maybe it is time to think about following in CTM's path and create a main LCM executable that by-passes the ModelEvaluator altogether.

But since we are in crunch time due to our deliverables, I propose we implement what @bgranzow outlines in his notes for the time being.

[bgranzow]

I'd like to make a few notes:

• I proposed the approach in my notes solely because it seemed like an easy first step.
• The approach in my notes is essentially the same thing that is done in this deal ii step, so someone has used it successfully.
• I've attempted to implement the approach in my notes in a prototyping code and have not yet been able to get it to work. Prof. Mark Shephard here at RPI brought up the point that pre-applying the Dirichlet boundary conditions to the initial guess x^(0) may actually be a poor initial guess for Newton's method, and may cause it to diverge. I'm in the process of trying to understand what is going wrong, and @ibaned has helped out with that process.

[lxmota]

@bgranzow Wrt to your 3rd bullet, that applying DBCs as initial guess may not be good. This is not unique to what you are proposing. If the DBCs impose large deformations wrt to the reference configuration, then yes, Newton's method will not converge.

That is why it is important to apply them gradually, specially in inelasticity.

Even in simulations where one expects abrupt initial displacements, they are applied gradually by using small time steps. Otherwise you have serious problems with element distortion.

I don't see why you would have trouble in your prototype code. Maybe your DBC displacements are too large?

I have a Matlab code where I can relatively quickly prototype this as well. I'll give it a try.

[lxmota]

@bgranzow It works! I just implemented what you wrote in your notes in our Matlab code, and it works even for Schwarz coupling!

Of course, since you form the whole system, it is slower than the implementation that uses the reduced system. That is why in the end we would like to have an implementation that uses the reduced system. There are many problems of interest that have many DOFs constrained as DBCs. In those problems, the reduced system would perform significantly faster.

In any case, the full and reduced system implementation give the same solution up to 10 to 12 significant digits at the end of Schwarz iterations where the residuals are very small. It appears that the whole system introduces some insignificant numerical noise. I would call that as equivalent solutions.

[lxmota]

I'm working on the implementation outlined above.

Since this is not strictly a condensed representation, I will open a new issue to track it.

This issue remains open until we implement the true condensed representation. Hopefully sooner rather than later. Maybe at the same time we remove the ModelEvaluator.