Dual mortar

[dewenyushu]

Reason

The dual basis is often utilized in the dual Lagrange approach. Here, the shape functions for the Lagrange multipliers (LMs) satisfy the biorthogonal condition across the interface. The biorthogonal condition results in a diagonalized local mass matrix, which allows to condense out the LMs and transform the saddle point problem in to a positive definite one.

Design

• [x] Add dual shape functions in libmesh (see issue 2550)

• [x] Proof of concept using diffusion with EqualValueConstraint

• [x] Enable dual mortar for mechanical contact

• [ ] Interface for condensing out the LMs

Impact

Improve the performance for mortar-based contact problems by producing a positive definite system and allows a broader range of preconditioners/solvers

[dewenyushu]

Mortar formulation results in a singular Jacobian for conforming mesh while using Dirichlet BCs. Two (nearly) zero singular values are associated with the Lagrange multipliers DOFs at the boundary. Refs https://github.com/idaholab/moose/pull/15216#issuecomment-630913039

[dewenyushu]

Here I am including some initial verifications of dual mortar implementation using MMS. The result is shown for a 2D diffusion problem in a unit square with analytical solution u(x,y) = sin(2pix)sin(2piy). The domain is composed of two subdomains (left and right), where solution continuity is enforced via using Lagrange multipliers (LMs) along the vertical interface. In this case the LM, \lambda=-2pisin(2pi*y).

The analytical solution is shown below:

The LM also converges with refined mesh:

The upper-right sub-block of Jacobian (which couples slave and LM) is simplified/diagonalized:

[dewenyushu]

@lindsayad @jwpeterson here are some initial verification results. I will be working on different element types and convergence studies as @lindsayad suggested.

[lindsayad]

@lindsayad @jwpeterson here are some initial verification results. I will be working on different element types and convergence studies as @lindsayad suggested.

Yea some L2/H1 error plots would be great. Very neat figures!

[jwpeterson]

@dewenyushu these results look very promising indeed, the sparsity patterns are especially interesting to see. I think the relative sizes of the u1 and u2 blocks indicate that you have chosen the LM/slave variable to be on the coarser mesh? From what I recall, this approach gives better results than the converse, but I'm also curious whether it makes any difference to your results. When the method is used in practice, sometimes it is not possible to pick a "coarse" side and obviously users can always set up problems incorrectly.

Also I was asking @lindsayad about it previously, but for the dual LM method, the dual basis functions are piecewise-discontinuous, so I would expect the LM solution to also be piecewise discontinuous. I can't tell from your figure if that's the case or whether the LM variable is actually continuous? Can you comment on this/help me understand this aspect better?

[lindsayad]

I think the relative sizes of the u1 and u2 blocks indicate that you have chosen the LM/slave variable to be on the coarser mesh? From what I recall, this approach gives better results than the converse, but I'm also curious whether it makes any difference to your results. When the method is used in practice, sometimes it is not possible to pick a "coarse" side and obviously users can always set up problems incorrectly.

I could have sworn that you advised picking the finer side as the slave side, but going back through your technical report, I actually don't see any statement regarding about a choice of side. In your test case 9.2 the slave side is about twice as fine as the master side.

[dewenyushu]

think the relative sizes of the u1 and u2 blocks indicate that you have chosen the LM/slave variable to be on the coarser mesh? From what I recall, this approach gives better results than the converse, but I'm also curious whether it makes any difference to your results. When the method is used in practice, sometimes it is not possible to pick a "coarse" side and obviously users can always set up problems incorrectly.

Yes, I am using a coarser mesh for the slave. Actually I did not choose it intentionally or know that would give me a better performance. So I agree it would be beneficial to try a finer mesh on the slave side. I will remember to report back what I see!

for the dual LM method, the dual basis functions are piecewise-discontinuous, so I would expect the LM solution to also be piecewise discontinuous. I can't tell from your figure if that's the case or whether the LM variable is actually continuous? Can you comment on this/help me understand this aspect better?

Great question! The dual basis I am using is actually piecewise continuous for the first order Lagrange. According to my understanding, the dual basis must satisfy the biorthogonal condition: where \phi is the primal basis and \psi is the dual basis. Therefore, theoretically there are many choices of the dual basis. Here, I am referring to Alexander Popp's paper and assuming the dual basis to be linear combinations of the primal basis, i.e.,

This is why my dual basis is actually continuous, and so is the LM field. I think this choice should have some benefit in terms of convergence, but I did not find it explicitly stated in literature. But yea, if you choose the primal basis to be discontinuous, then you will get discontinuous dual basis as well.

[dewenyushu]

Here is the 1st order and 2nd order LAGRANGE primal-dual basis comparison:

[jwpeterson]

In your biorthogonality equation, is the domain "\Gamma_e" a single edge in the discretization (Popp's Eqn (18)) or is it the entire contact interface (Popp's Eqn (17))?

I think you are referring to Popp's Eqn (17), in which case the dual basis will be continuous, but you have to solve a "global" (i.e. on the entire contact interface) L2-projection problem to get the continuous dual basis function coefficients. From what I can tell, you have actually implemented Popp's (18) in libmesh, i.e. local L2-projections on an element-by-element basis, which should in general lead to a discontinuous dual basis.

The dual Lagrange basis functions you posted above are discontinuous because, unlike the primal Lagrange basis functions, they do not go to zero at the element boundary on either side of the element... This is difficult to explain in words, but you can also look at e.g. Fig. 3.2 in (Wohmuth 2000) to try and get some idea of what I am talking about.

All that being said, your method seems to work so there is probably more that I am just not understanding about the approach...

[dewenyushu]

In your biorthogonality equation, is the domain "\Gamma_e" a single edge in the discretization (Popp's Eqn (18)) or is it the entire contact interface (Popp's Eqn (17))?

Here the \Gamma_e refers to a single edge in the discretization (i.e., Eqn (18)).

Yes, I agree with you. What I meant "continuous" previously was in terms of the local continuity in one element. Sorry about the confusion. Globally this could lead to discontinuous basis.

All that being said, your method seems to work so there is probably more that I am just not understanding about the approach...

I do not think I know more than you do haha. Thanks! All are really insightful questions and I will try to resolve them one by one through continuous research

[jwpeterson]

I could have sworn that you advised picking the finer side as the slave side, but going back through your technical report, I actually don't see any statement regarding about a choice of side. In your test case 9.2 the slave side is about twice as fine as the master side.

Yeah, Section 9.2 in the report was meant to test the robustness of the method when you do the "wrong" thing, i.e. make the slave side have a finer discretization than the master side. I guess I should have been more clear about that... anyway the results show that you lose a half power of h in the flux approximation error when doing this. I'm curious if there will be something similar for the dual method.

[lindsayad]

Yeah, Section 9.2 in the report was meant to test the robustness of the method when you do the "wrong" thing, i.e. make the slave side have a finer discretization than the master side.

I'm going to have to dig in an old branch and see whether the convergence studies that are shown in #13080 (which show super-convergence with the additional half power in most cases) had the finer side as the slave or master side. My inclination is that the finer side was the slave side, but I'm going to have to double check. You don't have a plot in that technical report that I missed that shows the mismatch switched, do you? E.g. twice as fine on the master side? Basically I'm looking for confirmation that the half-order reduction occurs because of making the "wrong" choice for the slave-master side, or whether the reduction simply comes about because of the large discrepancy in h between the sides (regardless of which one is slave or master).

[jwpeterson]

Globally this could lead to discontinuous basis.

OK, thanks for confirming this. My issue is therefore the following: if the dual basis functions "\Phi_j" are piecewise discontinuous, then in general the LM variable must also be piecewise discontinuous, do you agree? In other words, we could not use a (LINEAR, LAGRANGE) variable for the LM...

I suppose it's possible that you are introducing some kind of approximation by doing this, one which is very convenient for the implementation, and therefore if everything still converges at the expected rates, then everything is "fine". Anyway, I am curious about this for my understanding but the results are what actually matters of course.

[jwpeterson]

You don't have a plot in that technical report that I missed that shows the mismatch switched, do you?

No, I just have the results showing O(h^1.5) when the grids are non-conforming but approximately equal-sized elements, and the non-smooth interface case which shows a O(h^0.5) for P1-P1 formulation. I'll check through some of my other results to see whether I explicitly tested the mismatched case with slave side coarser.

[lindsayad]

Ok, looking back I did have the slave on the finer side, but the h mismatch was quite small

[lindsayad]

No, I just have the results showing an extra half power when the grids are non-conforming but approximately equal-sized elements, and the non-smooth interface case which shows a potential full power of h lost (for P1-P1 formulation). I'll check through some of my other results to see whether I explicitly tested the mismatched case with slave side coarser.

If results really are better when the slave side is coarser, then that would be awesome to know because I believe everywhere I do mortar (most importantly in mechanical contact, but also in thermal contact as well) I have the slave side as the finer side. Boy what a potential dumbo I am

[jwpeterson]

I'll check through some of my other results to see whether I explicitly tested the mismatched case with slave side coarser.

Definitely don't have any convergence results for the 2:1 mismatch with slave side coarser sitting around, but would definitely be interested to see what you get if you try it out.

[jwpeterson]

@lindsayad I managed to dig up some old results CSV files for the 2:1 mismatch case and regenerated the PDFs for them. In summary, I get flux error convergence at ~O(h^1.5) when the slave side is coarser and O(h) when the slave side is finer. There does not seem to be much effect on the primal variable's internal L2 or H1 error based on the choice of slave vs. master, only the convergence rate of the flux changes. I'd be curious to see if you get the same result for the dual LM method.

coarse_slave.pdf fine_slave.pdf

[dewenyushu]

@lindsayad @jwpeterson ok, here are some convergence performances for the dual mortar approach:

• conforming mesh:

• non-conforming mesh (finer master):

• non-conforming mesh (finer slave):

[dewenyushu]

The convergence looks pretty encouraging for the conforming case, where the both primal and LM have the same order of convergence for the 1st order Lagrange element. However, the 2nd order convergence for LM is converging much slower for some reason...

For the non-conforming mesh cases, it does not seem to show significant difference between choosing a finer slave side or a finer master side. However, it seems the LM convergences for non matching, second order cases are problematic.

[jwpeterson]

Thanks for making these plots. I have a few comments/questions:

1. For the dual method, it is required that the primal and Lagrange Multiplier variables are approximated with the same order, correct? That makes sense to me, but I think it means @lindsayad and I should focus on equal-order approximations for the non-dual method when attempting to make comparisons...
2. I don't have many results for the "conforming" mesh case in my test code, although it should be pretty easy to make some. It's interesting to me that you would get full 2nd-order accuracy for the flux in this case, but perhaps this is the natural extension of the behavior we normally see for linear elements, which ranges from somewhere between O(h) for "2:1 mismatch", and O(h^1.5) for "slight mismatch".
3. It's hard to tell but it looks like the error of the "finer slave" case is slightly lower than "finer master" for linear elements. For example the "finer master" starts off above 1, whereas the "finer slave" starts off below 1... plotting them on the same graph would help make it a bit more clear I guess.
4. Something seems not quite right for the Second-Order results, since we should not observe less than O(h^3) for the primal variable in my experience. Also for this case I would expect O(h^2) for the flux variable, so this may account for why the primal error does not converge at the correct rate.

[dewenyushu]

@jwpeterson Thanks for looking closely into these initial results!

1. For the dual method, it is required that the primal and Lagrange Multiplier variables are approximated with the same order, correct?

Yes. The primal and LMs are approximated with same order.

3. It's hard to tell but it looks like the error of the "finer slave" case is slightly lower than "finer master" for linear elements.

Yes, the error of LM for the 'finer slave' case is lower than the 'finer master' case with linear element. The values are around 1-2 orders of magnitude lower in this case. Yeah I agree putting them in the same plot would clearly show this. Will try to do that in future plots.

4. Something seems not quite right for the Second-Order results,

I agree. There must be something wrong with the second order cases. I am actually expecting equal-order convergence for the primal and the LMs while using a dual basis. So I am not sure if any of the 2nd order conforming cases is problem-free. I will need to look into that

[dewenyushu]

Here are comparisons between dual mortar and standard mortar:

• Conforming mesh:

• Non-conforming mesh (finer master):

[dewenyushu]

For the 2nd order nonconforming case, it seems the dual and non-dual produces different values of u, which is confusing to me.. is it because the lm does not converge in this case? Also, does this have to do with the mesh issue you pointed out in #15338? @lindsayad

[lindsayad]

What mismatch are you using in the non-conforming case? And for all these plots, are you using equal order discretization for primal and LM even for standard mortar?

[dewenyushu]

What mismatch are you using in the non-conforming case?

In terms of element size, the mismatch is h_slave:h_master=5:3

And for all these plots, are you using equal order discretization for primal and LM even for standard mortar?

Yes, I am using same order of approximation for both variables in all cases. Not sure if this makes sense for standard mortar in terms of stability. But I did this in order to compare with dual mortar.

[lindsayad]

As @jwpeterson said technically even mixed-order discretizations are unstable, so nothing really wrong with doing the equal-order.

It is very odd that you have that 2.273 slope for the standard LM in the non-conforming first order case. But given results in #15338 I believe there's some bugs in the grass that I need to snuff out

[dewenyushu]

It is very odd that you have that 2.273 slope for the standard LM in the non-conforming first order case.

Good point. I agree this lm super convergence also indicates something problematic.

But given results in #15338 I believe there's some bugs in the grass that I need to snuff out

lol like your metaphor! I think I will need to look into the 2nd order issue with dual mortar. But I guess your fix would also help with the performance on both sides. Very exciting. Let me know if I could be of some help!

[jwpeterson]

@dewenyushu I have one comment on your curve fits: in case you don't already do this, it's best if you compute the slope based only on the last 3-4 data points that are in the "asymptotic regime". In your Non-conforming, First-Order results, the red and green lines look almost parallel to me, so I think they should have basically exactly the same slope.

[dewenyushu]

@dewenyushu I have one comment on your curve fits: in case you don't already do this, it's best if you compute the slope based only on the last 3-4 data points that are in the "asymptotic regime". In your Non-conforming, First-Order results, the red and green lines look almost parallel to me, so I think they should have basically exactly the same slope.

Yes, that is a good point. I agree. It does not make sense to include points outside the asymptotic range for convergence study. Besides the one you just pointed out, the slope of the green curve in the conforming second order case should be based on the last few points, too. I will correct that

[lindsayad]

You could also rebase onto #15338. I've added some handy arguments to Andrew Slaughter's mms module that allow you to do things like limit the number of fitted points

On Wed, May 27, 2020 at 8:12 AM dewenyushu notifications@github.com wrote:

@dewenyushu https://github.com/dewenyushu I have one comment on your curve fits: in case you don't already do this, it's best if you compute the slope based only on the last 3-4 data points that are in the "asymptotic regime". In your Non-conforming, First-Order results, the red and green lines look almost parallel to me, so I think they should have basically exactly the same slope.

Yes, that is a good point. I agree. It does not make sense to include points outside the asymptotic range for convergence study. Besides the one you just pointed out, the slope of the green curve in the conforming second order case should be based on the last few points, too. I will correct that

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/idaholab/moose/issues/15215#issuecomment-634729023, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACOGA4HHPJNIFY2KBKZGBF3RTUUVDANCNFSM4MZCLYPA .

[dewenyushu]

@jwpeterson @lindsayad correct me if I am wrong: while integrating on the slave/master, we use the integration points (qp) from all mortar segments, not the slave/master element's standard qp, right?

If so, then we should expect to see 'super convergence' effect for the non-conforming case, as we are actually using more integration points (i.e., smaller size elements) than integrating directly on the slave or the master elements, which should imply a better convergence, right?