StrongDBC

[bartgol]

From what I understand, the StrongDBC evaluator in LCM diagonalizes also of the columns of the Jacobian corresponding to Dirichlet DOFs, is that correct? I have two questions:

1) Where does the rhs get modified accordingly? 2) Can we promote this evaluator to PHAL?

I would really like to have the ability to diagonalize the Dirichlet dofs columns (symmetry preservation, and conditioning being the biggest concerns). Assuming I understood the evaluator, I think it would be a big addition to the core of Albany.

[ikalash]

Re 1.) The residual (rhs) gets set to 0 in this evaluator at the dirichlet DOFs, and also the solution (x) gets set to the dirichlet value - see around like 56 of StrongDBC_Def.hpp.

Re 2.) I think so. I believe the reason Alejandro put it in LCM is to not be forced to support the BC for other non-LCM physics, which can take time.

[bartgol]

I mean, the rhs at the other dofs (the non-dirichlet ones) must be modified too. In particular, for each Dirichlet column j that is diagonalized, one needs to set every other non-dirichelt dof i to

b_i = b_i - J_ij * x_j.

I was wondering where this modification happens, since it's clearly not in the StrongDBC evaluator.

[ikalash]

@bartgol: this is only if we remove the constrained DOFs from the RHS. This is not done at the present time - the Jacobian retains the full system size. Ultimately for efficiency one would want to remove the constrained DOFs and do what you suggest but it hasn't been done yet.

[bgranzow]

@bartgol, I think the idea is that Albany pre-imposes the Dirichlet values in the solution vector, leading to a formulation where the solution increment, and thus for the residual, is always zero for nodes corresponding to DBC values. In your notation, this leads to x_j being zero, meaning those values don't need to get modified.

I wrote some notes about this here.

[bartgol]

@bgranzow Ah, I see what you mean. Since the solution is imposed, the linear system you need to solve is J*delta=-f, and delta_j=0. Therefore it makes no difference. You are right. I forgot that this was no general linear system, but it was indeed a Newton step system.

Also, I didn't notice that in the residual you did something different from PHAL_Dirichlet, namely, you actually set x_j=datum and f_j=0, rather than f=x_j-datum. Since the residual is always evaluated at least once before the Jacobian, then we are indeed sure that x_j=datum, so delta_j=0.

Ok, it makes sense now.

As I said, I think this should be moved to PHAL. I don't see any LCM-specific code in this bc, and I can only foresee an improvement of PHAL by its addition.

[mperego]

I agree that this should go into PHAL. Thanks @bgranzow for sharing your writeup.

[lxmota]

No objections here in moving it to PHAL. As @ikalash says, it would be better to reduce the size of the system by removing those DOFs altogether. Also, @ibaned introduced a much more efficient algorithm for this for MPI.

Feel free to improve it.

[mperego]

I'm not so sure that removing the DOFs to reduce the size of the system will save you that much in general, and it is not trivial. Are you thinking of assembling directly the problem using only the active DOFs and accounting for the terms associated to the lifting (as done in the FE theory)?

[lxmota]

For some solid mechanics problems, it is fairly common to have a lot of DBCs (entire surfaces or regions may be constrained), and thus the savings can be significant if the system could be reduced just as it is done in FE theory.

But yes, I agree that this wouldn't be trivial. So we went ahead and implemented @bgranzow's solution as described in his notes as a compromise.

[bartgol]

We could have a separate evaluator for strong dbc with dof removal. The user can then choose which dbc to use: DBC, SDBC, SDBCDR (terrible acronym probably, not a proposal).

[ibaned]

@bartgol we don't currently remove DOFs, that is already what the existing SDBC is. The trouble is its difficult to remove DOFs, and no one yet needs it enough to write the code for it.

[bartgol]

@ibaned Well, SDBC already does something more than DBC, namely diagonalizes the column, which has some advantages (simmetry and conditioning). I would even argue about making SDBC (as it is now, without dof removal) the new "standard" DBC (I mean the general x_i = g_i bc from the PDE textbooks). That is, unless someone has a strong reason not to do it, and keep the current DBC (i.e., diagonalize only the row).

[ibaned]

Sounds good to me, other people can comment on possible issues

[mperego]

I guess the only issue with the new implementation is the cost of "diagonalizing" the matrix. In theory, if the initial guess satisfies the dirichlet conditions, both the diagonalized and nondiagonalized jacobian should preform in the exact same way when using a (non preconditioned) Krylov method. But I guess preconditioning might change this. So, I think I'm also for making the new implementation the default one.

[bartgol]

Yes, the only "objection" would be the extra cost of the column diagonalization. However, I would argue that that cost is negligible in a "standard" application.

[bgranzow]

My concern would be the possibility of many regression test values getting slightly bumped, so that existing tests no longer pass.

[bartgol]

That's a good point. I guess we can always change the expected values, but it may take some time.

[lxmota]

I strongly vote to make SDBC the default.

[ikalash]

Per our discussion at the ProSPect meeting this morning, I volunteer to move the SDBC to the Albany/src/evaluator class from the LCM directory to make it accessible to all Albany physics.

Question, also per the discussion at our meeting this morning: what would be a good name for the evaluator, since StrongDBC is somewhat of a misnomer? @bartgol suggested SymmetricDBC. Thoughts?

[mperego]

Thanks @ikalash. Would it be hard to make the symmetrization (or diagonalization) optional so that it could be set in the input file? If so we can get rid of the current implementation.

[ikalash]

@mperego : oh, yes, that was the plan - not to overwrite PHAL_Dirichlet* but to have also the option to use the SDBC. I know we had the discussion about which should be the default above. I think the first step would be to move the evaluator; we can talk about defaults later.

[mperego]

@ikalash I was actually proposing to override PHAL_Dirichlet*, if possible.. to avoid code duplications and having to maintain both classes.

[bartgol]

I'm also ok making it the default DBC behavior. But we need a broad consensus for that. As @ikalash said, I'm for 'Symmetric', which allows us to still use the SDBC acronym, but I think is more verbose on what it does. Otherwise SymmetryPreserving, but it gets a little long. In the end, I'm fine also with 'Strong', so long as we get it in PHAL.

[ikalash]

Ah, I see, so you mean to have one class with both implementations, with each one possible to select at run-time? Yes, this is probably the best option and I think it can be done. In this case, the only name to decide on is what to put in the parameterlist to define which form of BC to use. Currently it's "DBC" for regular Dirichlet BC, and "SDBC". I think this is probably OK, since S could stand for "Symmetric" :).

[mperego]

@ikalash That would be cool.

[agsalin]

How does SDBC work? Isn't the solution vector a const in the model evaluator? how do you set its values if it is a const?

My main concern is that sensitivities, such as tangent predictor, might not be correct if this isn't done consistently. Have you tested sensitivities w.r.t. DBC values for SDBC?

[ikalash]

Very good question. @lxmota did a cast from const to non-const to make the modification in x:

Teuchos::RCP x = Teuchos::rcpFromRef(const_cast<Tpetra_Vector &>(*dirichlet_workset.xT));

This is dangerous, I agree, but I am not sure how to get around it (@lxmota, feel free to comment). We did not try sensitivities, no, which are not of interest in LCM (I believe). This was one argument for having the capability in just LCM, I guess.

[lxmota]

Hi all,

Yes, SDBC is totally a misnomer, but one that is short and sweet.

We did not test for sensitivities. Irina is right that I'm doing a forceful and ugly casting to get this to work.

For LCM, we do use sensitivities. So for LCM this is ok.

Alejandro

On Nov 9, 2017, at 13:05, Irina K. Tezaur notifications@github.com<mailto:notifications@github.com> wrote:

Very good question. @lxmotahttps://github.com/lxmota did a cast from const to non-const to make the modification in x:

Teuchos::RCP x = Teuchos::rcpFromRef(const_cast<Tpetra_Vector &>(*dirichlet_workset.xT));

This is dangerous, I agree, but I am not sure how to get around it (@lxmotahttps://github.com/lxmota, feel free to comment). We did not try sensitivities, no, which are not of interest in LCM (I believe). This was one argument for having the capability in just LCM, I guess.

- You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/gahansen/Albany/issues/208#issuecomment-343291181, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AGykgTuPzOfZ2ddM9cp2jP0XBWH0u2kmks5s02kigaJpZM4QPz73.

[ikalash]

I can test what happens with sensitivities once I move the BC to the main part of Albany. I would say if it doesn't work with sensitivities, we should not make it the default.

[bartgol]

Uhm, I didn't realize there's a const cast. I'm not saying that's a no-no, cause there's always a valid use for each feature of the language, but perhaps we can think whether there is another way.

I'm going to throw an alternative to the const cast: leave the behavior of DBC on the Dirichlet nodes (i.e., f_i = x_i - datum_i), and for all the other rows, do the update

f_i -= J_ij * (x_j - datum_j),

for all dirichlet dof's j and all non-dirichlet dof's i. This is basically doing f = f - J*f_D, where f_D is the restriction of f to the Dirichlet nodes. I believe this should be done outside the evaluator chain, right before the call to the solver, I think NOX has some callback functionalities for this. I think @ibaned and @ikalash were discussing with Roger this when they needed to scale the residual. Perhaps we can use something like that here too. This would be a single mat-vec product, so not costly, avoiding the need of const-casting x.

But I'm wondering whether this makes it all more complicated, and we should just be fine with a const cast.

[lxmota]

The above is fine with me for the DBC that would be alternate to the original Albany DBC and that would go in PHAL.

But then, let's leave alone the current SDBC within LCM. It is much more important to us that it simulates what we see as the correct application of the DBC than many other things. In lieu of the reduced system, this is the next best thing. Many of our problems now actually work, including dynamics, and those that worked with the original DBC see better numerical behavior due to better conditioning.

[mperego]

@lxmota, I would bet that Luca's approach would give the same result (almost to machine precision) of the SDBC approach, if modified such that on the diagonal you have the same value you have in the current approach (something like f_i = J_ii (x_i - datum_i) ). Anyway, I'm not necessarily opting for that approach.

[ibaned]

I believe @bgranzow and I considered both approaches and had good reasons for selecting what we did, though I can't recall the details at the moment.

[mperego]

@ibaned for sure it is simpler because you do not have to perform this f_i -= J_ij * (x_j - datum_j).

[ikalash]

I suppose if @bartgol's approach is implemented in phal dirichlet we could do numerical studies comparing its performance to the lcm one.

[ikalash]

I spent some time this weekend looking at the SDBC, and would propose for now:

• Renaming it and moving it to the src/evaluators class (as is). I think having it as a separate evaluator will be cleaner than integrating it into PHAL::Dirichlet.
• Removing ALBANY_LCM ifdef guards around SDBC-specific code in Albany::Application so the SDBC runs correctly for all physics.
• Adding asserts / throwing of exceptions to the new evaluator in the case one tries to run it in an unsupported regime (e.g., where Tangent and DistParam are required, which are currently not implemented in the SDBC).

We will need to check/double check if the following work with the SDBC:

• Sensitivities (probably not since Tangent specialization is not implemented).
• If works with Epetra executable, Albany.
• If works with "Single Step" nonlinear solver.

[ikalash]

I have extended SDBC to be callable from all Albany problems regardless of whether the code has been build w/ or w/o LCM. If any problems are encountered using the BC, please open a new Albany issue.

[bartgol]

Since SDBC are now in the main Albany folders, I think we can close this for now. Please, reopen if there's more to discuss about this.