Enable usage of dual shape functions for the Lagrange multipliers #15216

dewenyushu · 2020-05-04T21:29:40Z

Address part of the issues 2550 and #15215

Note: To be merged after this PR in libmesh.

dewenyushu · 2020-05-04T22:25:25Z

Hi @bwspenc, would you mind helping me review this PR?

dewenyushu · 2020-05-11T14:36:02Z

framework/src/base/Assembly.C

+    FEShapeData * fesd = _fe_shape_data_dual_lower[fe_type];
+
+    fe_lower->set_calculate_dual(true);
+    fe_lower->reinit(elem, pts, weights);


@lindsayad here is where the dualLowerD elems are re-initialized with specific pts. dual_coeff is not computed in this case thus causes a seg fault...

See my comment libMesh/libmesh#2552 (comment). Is it valid to calculate the dual coefficients with points that don't correspond to a quadrature rule for that element? Is the correct solution to call fe_lower->reinit(elem) first before calling fe_lower->reinit(elem, pts, weights)?

The reinitLowerDElemRef and reinitLowerDElemDualRef are supposed to compute the quadrature rule from the customized integration points for the slave side. Apologies for the earlier confusion.. I also just realized that the dual_coeff should be computed using the same integration pts as the primal. So it is the "full" element for general cases, and becomes complicated for the slave elements when we are integrating using projected points (custom_xi1_pts). And I guess the local matrices totally depend on the 'correctness' of the customized integration points in order to be non-singular..

So what are you proposing? One can easily imagine (just thinking about the linear Lagrange case for simplicity) that you could have a mortar elem which occupies a very tiny fraction of the total slave element "volume", e.g. where the custom_xi1_pts are clustered on the farthest LHS or RHS of the slave element (or in the middle, or whatever).

My question to you is: is the dual method valid for any set of "integration" points, putting aside "practical" concerns like whether the biorthogonality condition yields a non-singular matrix?

Theoretically it should be valid for the set of integration points that are the same of the primal element. I suppose this is a better way to put it..

Is "primal element" language that people use? What is a "primal element"?

Good question. The "primal and dual" concept is something people often use. The "dual basis" is relative to its "primal basis". Here I think it is more appropriate to say the "element that uses the primal basis" and the "element that uses the dual basis"

dewenyushu · 2020-05-13T18:34:47Z

Hi @lindsayad, it seems libmesh is not updated for the MOOSE checks here.. do you know when this will be updated, and who would be the best person to talk to about this?

lindsayad · 2020-05-13T18:48:27Z

I am the best person to talk to. See #15182

dewenyushu · 2020-05-13T19:51:12Z

I am the best person to talk to. See #15182

Great! So how frequent will the checks be based on an updated libmesh? (I assume that is when I will be able to pass the checks here, right?)

lindsayad · 2020-05-13T19:54:06Z

you'll have to wait until the libmesh update is merged and until it passes testing on next and is "merged" into devel. Alternatively as soon as the update is merged into next, you could rebase your PR onto next

…

On Wed, May 13, 2020 at 12:51 PM dewenyushu ***@***.***> wrote: I am the best person to talk to. See #15182 <#15182> Great! So how frequent will the checks be based on an updated libmesh? (I assume that is when I will be able to pass the checks here, right?) — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#15216 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACOGA4GFVDQP2WTGZPMHQVDRRL235ANCNFSM4MZDFBEA> .

moosebuild · 2020-05-18T20:02:07Z

Job Documentation on be92a2c wanted to post the following:

View the site here

This comment will be updated on new commits.

dewenyushu · 2020-05-18T20:55:57Z

Hi @lindsayad, would you mind taking a look at the failing parallel test? I am not sure why the solver fails using 2cores... It does not help by changing to other solvers, like -pc_factor_mat_solver_type superlu_dist so it should not be the same problem as with libMesh/libmesh#2552 (comment)

lindsayad · 2020-05-19T01:32:03Z

Lol I can't get this to fail locally. Are you using conda libmesh? Oh I was in the wrong directory lol. non-conforming one looks great!

lindsayad · 2020-05-19T02:06:07Z

Well if I turn on -pc_type svd -pc_svd_monitor, then it appears that there are two singular values. Moreover, if I do Variables/lm/use_dual=false, there are still two singular values. These singularities don't go away if you use -snes_fd. For fun I tried the equalgradient.i test and that one also has a singular value. I checked the matrix using -ksp_view_pmat and there are no zero rows, so the singularities must be due to linear combinations. I can't spot any of them off the top of my head, but @fdkong is the matrix reading expert. There are no singularities in the non-conforming mesh tests.

@jwpeterson off the top of your head do you see any reason why mortar with a perfect conforming mesh would produce a singular matrix? These tests are just doing solution continuity.

@fdkong I don't suppose you can immediately see why this matrix is singular can you?

row 0: (0, 2.66667)  (1, -0.333333)  (2, -0.333333)  (3, -0.333333)  (4, -0.333333)  (5, -0.333333)  (12, -0.333333)  (13, -0.333333)  (14, -0.333333)  (30, 0.)  (31, 0.)  (34, 0.) 
row 1: (0, 0.)  (1, 1.)  (2, 0.)  (3, 0.)  (12, 0.)  (14, 0.) 
row 2: (0, 0.)  (1, 0.)  (2, 1.)  (3, 0.) 
row 3: (0, 0.)  (1, 0.)  (2, 0.)  (3, 1.)  (4, 0.)  (5, 0.)  (30, 0.)  (31, 0.) 
row 4: (0, 0.)  (3, 0.)  (4, 1.)  (5, 0.)  (30, 0.)  (31, 0.) 
row 5: (0, -0.333333)  (3, -0.333333)  (4, -0.166667)  (5, 1.33333)  (12, -0.333333)  (13, -0.166667)  (30, 0.)  (31, -0.25)  (34, -1.38778e-17) 
row 6: (6, 2.66667)  (7, -0.333333)  (8, -0.333333)  (9, -0.333333)  (10, -0.333333)  (11, -0.333333)  (12, -0.333333)  (13, -0.333333)  (14, -0.333333)  (32, 0.)  (33, 0.)  (34, 0.) 
row 7: (6, 0.)  (7, 1.)  (8, 0.)  (9, 0.)  (10, 0.)  (11, 0.)  (32, 0.)  (33, 0.) 
row 8: (6, 0.)  (7, 0.)  (8, 1.)  (9, 0.) 
row 9: (6, 0.)  (7, 0.)  (8, 0.)  (9, 1.)  (12, 0.)  (14, 0.) 
row 10: (6, -0.333333)  (7, -0.333333)  (10, 1.33333)  (11, -0.166667)  (12, -0.333333)  (13, -0.166667)  (32, -0.25)  (33, -1.38778e-17)  (34, 0.) 
row 11: (6, 0.)  (7, 0.)  (10, 0.)  (11, 1.)  (32, 0.)  (33, 0.) 
row 12: (0, -0.333333)  (1, -0.333333)  (5, -0.333333)  (6, -0.333333)  (9, -0.333333)  (10, -0.333333)  (12, 2.66667)  (13, -0.333333)  (14, -0.333333)  (31, 0.)  (32, 0.)  (34, 0.) 
row 13: (0, -0.333333)  (5, -0.166667)  (6, -0.333333)  (10, -0.166667)  (12, -0.333333)  (13, 1.33333)  (31, 0.)  (32, -1.38778e-17)  (34, -0.25) 
row 14: (0, 0.)  (1, 0.)  (6, 0.)  (9, 0.)  (12, 0.)  (14, 1.) 
row 15: (15, 1.)  (16, 0.)  (17, 0.)  (18, 0.) 
row 16: (15, 0.)  (16, 1.)  (17, 0.)  (18, 0.)  (28, 0.)  (29, 0.) 
row 17: (15, -0.333333)  (16, -0.333333)  (17, 2.66667)  (18, -0.333333)  (19, -0.333333)  (20, -0.333333)  (27, -0.333333)  (28, -0.333333)  (29, -0.333333)  (30, 0.)  (31, 0.)  (34, 0.) 
row 18: (15, 0.)  (16, 0.)  (17, 0.)  (18, 1.)  (19, 0.)  (20, 0.)  (30, 0.)  (31, 0.) 
row 19: (17, 0.)  (18, 0.)  (19, 1.)  (20, 0.)  (30, 0.)  (31, 0.) 
row 20: (17, -0.333333)  (18, -0.333333)  (19, -0.166667)  (20, 1.33333)  (27, -0.166667)  (28, -0.333333)  (30, 0.)  (31, 0.25)  (34, 1.38778e-17) 
row 21: (21, 1.)  (22, 0.)  (23, 0.)  (24, 0.) 
row 22: (21, 0.)  (22, 1.)  (23, 0.)  (24, 0.)  (25, 0.)  (26, 0.)  (32, 0.)  (33, 0.) 
row 23: (21, -0.333333)  (22, -0.333333)  (23, 2.66667)  (24, -0.333333)  (25, -0.333333)  (26, -0.333333)  (27, -0.333333)  (28, -0.333333)  (29, -0.333333)  (32, 0.)  (33, 0.)  (34, 0.) 
row 24: (21, 0.)  (22, 0.)  (23, 0.)  (24, 1.)  (28, 0.)  (29, 0.) 
row 25: (22, 0.)  (23, 0.)  (25, 1.)  (26, 0.)  (32, 0.)  (33, 0.) 
row 26: (22, -0.333333)  (23, -0.333333)  (25, -0.166667)  (26, 1.33333)  (27, -0.166667)  (28, -0.333333)  (32, 0.25)  (33, 1.38778e-17)  (34, 0.) 
row 27: (17, -0.333333)  (20, -0.166667)  (23, -0.333333)  (26, -0.166667)  (27, 1.33333)  (28, -0.333333)  (31, 0.)  (32, 1.38778e-17)  (34, 0.25) 
row 28: (16, -0.333333)  (17, -0.333333)  (20, -0.333333)  (23, -0.333333)  (24, -0.333333)  (26, -0.333333)  (27, -0.333333)  (28, 2.66667)  (29, -0.333333)  (31, 0.)  (32, 0.)  (34, 0.) 
row 29: (16, 0.)  (17, 0.)  (23, 0.)  (24, 0.)  (28, 0.)  (29, 1.) 
row 30: (0, 0.)  (3, 0.)  (4, -0.125)  (5, 0.)  (17, 0.)  (18, 0.)  (19, 0.125)  (20, 0.)  (30, 0.)  (31, 0.) 
row 31: (0, 0.)  (3, 0.)  (4, -1.38778e-17)  (5, -0.25)  (12, 0.)  (13, 0.)  (17, 0.)  (18, 0.)  (19, 1.38778e-17)  (20, 0.25)  (27, 0.)  (28, 0.)  (30, 0.)  (31, 0.)  (34, 0.) 
row 32: (6, 0.)  (7, 0.)  (10, -0.25)  (11, 0.)  (12, 0.)  (13, -1.38778e-17)  (22, 0.)  (23, 0.)  (25, 0.)  (26, 0.25)  (27, 1.38778e-17)  (28, 0.)  (32, 0.)  (33, 0.)  (34, 0.) 
row 33: (6, 0.)  (7, 0.)  (10, -1.38778e-17)  (11, -0.125)  (22, 0.)  (23, 0.)  (25, 0.125)  (26, 1.38778e-17)  (32, 0.)  (33, 0.) 
row 34: (0, 0.)  (5, -1.38778e-17)  (6, 0.)  (10, 0.)  (12, 0.)  (13, -0.25)  (17, 0.)  (20, 1.38778e-17)  (23, 0.)  (26, 0.)  (27, 0.25)  (28, 0.)  (31, 0.)  (32, 0.)  (34, 0.)

lindsayad · 2020-05-19T06:14:18Z

If I make the primal variable u second order, and I make the lm non-dual (I got a segmentation fault when running with Mesh/second_order=true and Variables/lm/use_dual=true), then the problem is no longer singular. What's interesting is that an equal-order discretization on the non-conforming mesh is nonsingular (both for dual and non-dual LM)

jwpeterson · 2020-05-19T13:38:55Z

Hmm... I don't think I ever tested a fully "conforming mesh" case, but I should have. I did have a few meshes where some of the nodes just happened to align, in order to make sure that was working OK. Recall that for the non-dual, equal-order Lagrange multiplier formulation, a sufficient condition for LBB-stability is:

h_lambda >= C * h_u

with C > 1, where h_lambda and h_u are the Lagrange multiplier and primal variable mesh sizes, respectively. In the conforming mesh case you are right at the C=1 limit so there's no guarantee it would produce stable solutions, but I'm also a bit surprised that it would produce a singular matrix. In the dual Lagrange multiplier formulation, are you "condensing" out the DOFs for the LM variable and still getting a singular matrix? Or does the current solver retain all those DOFs like in the non-dual LM formulation?

dewenyushu · 2020-05-19T14:09:50Z

Lol I can't get this to fail locally. Are you using conda libmesh? Oh I was in the wrong directory lol. non-conforming one looks great!

I am using my own libmesh. Yeah, it is also confusing the non-conforming case works fine!

jwpeterson · 2020-05-19T14:18:19Z

@lindsayad I didn't look at the mesh used for this specific problem but if it's still singular for say a 2x2 conforming mesh of Quad4's with the interface through the center, it might be easier to see the linear combinations and/or visualize the null space.

I am not sure why the solver fails using 2cores

@dewenyushu are you saying that the test works OK in serial but not parallel? @lindsayad mentioned using -pc_type svd -pc_svd_monitor so IIRC that's a serial solver, so I'm guessing the test doesn't work in serial or parallel...

dewenyushu · 2020-05-19T14:20:15Z

Recall that the non-dual, equal-order Lagrange multiplier formulation is not LBB-stable unless

h_lambda >= C * h_u

with C > 1, where h_lambda and h_u are the Lagrange multiplier and primal variable mesh sizes, respectively. In the conforming mesh case you are right at the C=1 limit so I would not expect it to produce stable solutions

Thanks for the explanation! I have been having the question of the stability of the schemes while not been able to find a good answer.

In the dual Lagrange multiplier formulation, are you "condensing" out the DOFs for the LM variable and still getting a singular matrix?

No I am yet to condense out the LM DOFs. I am starting to work on the "condensing" in another PR. I doubt if the condensing process would change the singularity of the matrix. I will have a try!

dewenyushu · 2020-05-19T14:21:43Z

does the current solver retain all those DOFs like in the non-dual LM formulation

Yes, the solver retains all DOFs

jwpeterson · 2020-05-19T14:27:05Z

I have been having the question of the stability of the schemes while not been able to find a good answer.

It's a surprisingly complicated issue. There is a Tech. Report which has some basic information in Sections 3 and 4 on stability and error estimates for the non-dual method that you may find interesting.

dewenyushu · 2020-05-19T14:37:23Z

@lindsayad I didn't look at the mesh used for this specific problem but if it's still singular for say a 2x2 conforming mesh of Quad4's with the interface through the center, it might be easier to see the linear combinations and/or visualize the null space.

I am not sure why the solver fails using 2cores

@dewenyushu are you saying that the test works OK in serial but not parallel? @lindsayad mentioned using -pc_type svd -pc_svd_monitor so IIRC that's a serial solver, so I'm guessing the test doesn't work in serial or parallel...

Yes, the test works fine in serial on my local machine. But I am using -pc_type lu. emm I just tried -pc_type svd -pc_svd_monitor it works fine both in serial and using 2 cores for me now...

dewenyushu · 2020-05-19T14:41:52Z

@lindsayad I didn't look at the mesh used for this specific problem but if it's still singular for say a 2x2 conforming mesh of Quad4's with the interface through the center, it might be easier to see the linear combinations and/or visualize the null space.

I am not sure why the solver fails using 2cores

@dewenyushu are you saying that the test works OK in serial but not parallel? @lindsayad mentioned using -pc_type svd -pc_svd_monitor so IIRC that's a serial solver, so I'm guessing the test doesn't work in serial or parallel...

Yes, the test works fine in serial on my local machine. But I am using -pc_type lu. emm I just tried -pc_type svd -pc_svd_monitor it works fine both in serial and using 2 cores for me now...

But I do see 2 singular values for both serial and parallel cases. Same as Alex reported

jwpeterson · 2020-05-19T14:44:54Z

I just tried -pc_type svd -pc_svd_monitor it works fine both in serial and using 2 cores for me now...

Hmm... according to this it looks like -pc_type svd works in parallel by making a separate copy of the matrix on each process.

But I think what we are really interested in is the output of -pc_svd_monitor, it should print some information about singular values, etc.? Also, before debugging too much I guess it would be good if you and @lindsayad can agree on a test case that fails in the same way for both of you.

lindsayad · 2020-05-19T15:03:28Z

But I think what we are really interested in is the output of -pc_svd_monitor, it should print some information about singular values, etc.? Also, before debugging too much I guess it would be good if you and @lindsayad can agree on a test case that fails in the same way for both of you.

Correct. The conforming case runs "fine" in serial but it's still singular. Just getting lucky.

here's no guarantee it would produce stable solutions, but I'm also a bit surprised that it would produce a singular matrix.

This is why I tried second order for the primal variable. All the tests I added originally for mortar used a mixed order discretization for saddle point problems. But yea I also didn't expect the theorized instability to necessarily translate to a singular matrix. And the non-conforming case with equal-order is not singular (in the "non-conforming" case the interface end-points are perfectly aligned, but those are the only points of alignment; it's your old 2 element on one side, 3 elements on another mesh).

The 2x2 mesh suggestion is a good one. I'll give that shot and report back.

dewenyushu · 2020-05-19T15:15:03Z

But I think what we are really interested in is the output of -pc_svd_monitor, it should print some information about singular values, etc.? Also, before debugging too much I guess it would be good if you and @lindsayad can agree on a test case that fails in the same way for both of you.

Yes, I agree. I think @lindsayad and I both observe the same performance now for this test moose/test/tests/mortar/continuity-2d-conforming/dual-conforming.i. The output of -pc_svd_monitor is:

SVD: condition number 1.535552857557e+18, 2 of 35 singular values are (nearly) zero
      SVD: smallest singular values: 2.075340939175e-18 7.243352805199e-17 1.573140358615e-02 4.374676748655e-02 1.096517322962e-01
      SVD: largest singular values : 2.755787629229e+00 2.819886851105e+00 2.820940606716e+00 3.186721454198e+00 3.186795709555e+00

lindsayad · 2020-05-19T15:54:33Z

Ok I still can't see the linear combinations in this guy, but I know what's happening

ow 0: (0, 1.)  (1, 0.)  (2, 0.)  (3, 0.)  (12, 0.)  (13, 0.) 
row 1: (0, 0.)  (1, 1.)  (2, 0.)  (3, 0.)  (12, 0.)  (13, 0.) 
row 2: (0, -0.416667)  (1, -0.583333)  (2, 1.66667)  (3, 0.333333)  (4, -0.583333)  (5, -0.416667)  (12, -0.0833333)  (13, -0.333333)  (14, -0.0833333) 
row 3: (0, 0.)  (1, 0.)  (2, 0.)  (3, 1.)  (4, 0.)  (5, 0.)  (12, 0.)  (13, 0.)  (14, 0.) 
row 4: (2, 0.)  (3, 0.)  (4, 1.)  (5, 0.)  (13, 0.)  (14, 0.) 
row 5: (2, 0.)  (3, 0.)  (4, 0.)  (5, 1.)  (13, 0.)  (14, 0.) 
row 6: (6, 1.)  (7, 0.)  (8, 0.)  (9, 0.)  (12, 0.)  (13, 0.) 
row 7: (6, 0.)  (7, 1.)  (8, 0.)  (9, 0.)  (12, 0.)  (13, 0.) 
row 8: (6, 0.)  (7, 0.)  (8, 1.)  (9, 0.)  (10, 0.)  (11, 0.)  (12, 0.)  (13, 0.)  (14, 0.) 
row 9: (6, -0.583333)  (7, -0.416667)  (8, 0.333333)  (9, 1.66667)  (10, -0.416667)  (11, -0.583333)  (12, 0.0833333)  (13, 0.333333)  (14, 0.0833333) 
row 10: (8, 0.)  (9, 0.)  (10, 1.)  (11, 0.)  (13, 0.)  (14, 0.) 
row 11: (8, 0.)  (9, 0.)  (10, 0.)  (11, 1.)  (13, 0.)  (14, 0.) 
row 12: (0, 0.)  (1, -0.166667)  (2, -0.0833333)  (3, 0.)  (6, 0.166667)  (7, 0.)  (8, 0.)  (9, 0.0833333)  (12, 0.)  (13, 0.) 
row 13: (0, 0.)  (1, -0.0833333)  (2, -0.333333)  (3, 0.)  (4, -0.0833333)  (5, 0.)  (6, 0.0833333)  (7, 0.)  (8, 0.)  (9, 0.333333)  (10, 0.)  (11, 0.0833333)  (12, 0.)  (13, 0.)  (14, 0.) 
row 14: (2, -0.0833333)  (3, 0.)  (4, -0.166667)  (5, 0.)  (8, 0.)  (9, 0.0833333)  (10, 0.)  (11, 0.166667)  (13, 0.)  (14, 0.)

This is two elements on the left, two on the right. Two singular values. The singularity doesn't have anything to do with mortar per se. The problem is that at interface end points we are simultaneously specifying Dirichlet boundary conditions and solution continuity through a Lagrange Multiplier. So we have something like these equations:

u0_slave = f0_slave (Dirichlet)
u0_master = f0_master (Dirichlet)
u0_slave = u0_master (LM)

u1_slave = f1_slave (Dirichlet)
u1_master = f1_master (Dirichlet)
u1_slave = u1_master (LM)

There ya go, two singularities. The moment I remove the competing Dirichlet conditions (and implicitly replace them with natural boundary conditions), the singularities go away. Similarly, if the primal variable is made to be second order, then the two end-point LMs that previously corresponded to the two singularities have another primal dof to contribute to, and the singularities also disappear.

For the "non-conforming" input file with aligning end-points, we don't see this same issue because @dewenyushu is applying Neumann conditions instead of competing Dirichlet ones.

jwpeterson · 2020-05-19T16:03:58Z

OK, that's an interesting result. I guess you'll need to be careful to avoid this issue in general problems as well since it seems like it could be a pretty common thing to have Dirichlet BCs on one or both "ends" of a mortar region? Possibly it could be detected when the lower-dimensional elements are initially added, or something along those lines...

dewenyushu · 2020-05-19T16:17:31Z

Ok I still can't see the linear combinations in this guy, but I know what's happening
ow 0: (0, 1.)  (1, 0.)  (2, 0.)  (3, 0.)  (12, 0.)  (13, 0.) 
row 1: (0, 0.)  (1, 1.)  (2, 0.)  (3, 0.)  (12, 0.)  (13, 0.) 
row 2: (0, -0.416667)  (1, -0.583333)  (2, 1.66667)  (3, 0.333333)  (4, -0.583333)  (5, -0.416667)  (12, -0.0833333)  (13, -0.333333)  (14, -0.0833333) 
row 3: (0, 0.)  (1, 0.)  (2, 0.)  (3, 1.)  (4, 0.)  (5, 0.)  (12, 0.)  (13, 0.)  (14, 0.) 
row 4: (2, 0.)  (3, 0.)  (4, 1.)  (5, 0.)  (13, 0.)  (14, 0.) 
row 5: (2, 0.)  (3, 0.)  (4, 0.)  (5, 1.)  (13, 0.)  (14, 0.) 
row 6: (6, 1.)  (7, 0.)  (8, 0.)  (9, 0.)  (12, 0.)  (13, 0.) 
row 7: (6, 0.)  (7, 1.)  (8, 0.)  (9, 0.)  (12, 0.)  (13, 0.) 
row 8: (6, 0.)  (7, 0.)  (8, 1.)  (9, 0.)  (10, 0.)  (11, 0.)  (12, 0.)  (13, 0.)  (14, 0.) 
row 9: (6, -0.583333)  (7, -0.416667)  (8, 0.333333)  (9, 1.66667)  (10, -0.416667)  (11, -0.583333)  (12, 0.0833333)  (13, 0.333333)  (14, 0.0833333) 
row 10: (8, 0.)  (9, 0.)  (10, 1.)  (11, 0.)  (13, 0.)  (14, 0.) 
row 11: (8, 0.)  (9, 0.)  (10, 0.)  (11, 1.)  (13, 0.)  (14, 0.) 
row 12: (0, 0.)  (1, -0.166667)  (2, -0.0833333)  (3, 0.)  (6, 0.166667)  (7, 0.)  (8, 0.)  (9, 0.0833333)  (12, 0.)  (13, 0.) 
row 13: (0, 0.)  (1, -0.0833333)  (2, -0.333333)  (3, 0.)  (4, -0.0833333)  (5, 0.)  (6, 0.0833333)  (7, 0.)  (8, 0.)  (9, 0.333333)  (10, 0.)  (11, 0.0833333)  (12, 0.)  (13, 0.)  (14, 0.) 
row 14: (2, -0.0833333)  (3, 0.)  (4, -0.166667)  (5, 0.)  (8, 0.)  (9, 0.0833333)  (10, 0.)  (11, 0.166667)  (13, 0.)  (14, 0.)
This is two elements on the left, two on the right. Two singular values. The singularity doesn't have anything to do with mortar per se. The problem is that at interface end points we are simultaneously specifying Dirichlet boundary conditions and solution continuity through a Lagrange Multiplier. So we have something like these equations:

u0_slave = f0_slave (Dirichlet)
u0_master = f0_master (Dirichlet)
u0_slave = u0_master (LM)

u1_slave = f1_slave (Dirichlet)
u1_master = f1_master (Dirichlet)
u1_slave = u1_master (LM)

There ya go, two singularities. The moment I remove the competing Dirichlet conditions (and implicitly replace them with natural boundary conditions), the singularities go away. Similarly, if the primal variable is made to be second order, then the two end-point LMs that previously corresponded to the two singularities have another primal dof to contribute to, and the singularities also disappear.

For the "non-conforming" input file with aligning end-points, we don't see this same issue because @dewenyushu is applying Neumann conditions instead of competing Dirichlet ones.

Nice findings! I got it that there are competing Dirichlet BC on the slave at the two ends. However, it is not clear to me how that causes a singularity, and how 'replacing it with a natural BC' will fix the issue.. @lindsayad would you mind explaining more in this regard? Or maybe we can talk offline if that is easier

Get conforming and nonconforming cases to converge, checked Jacobian pattern Fix broken tests, code restructure, and add tests Address issue idaholab#15215

Remove test for the conforming mesh Rebase against devel Address issue idaholab#15215

dewenyushu · 2020-05-20T01:04:26Z

@lindsayad this PR is finally ready for review!

lindsayad

I think this is a really nice start, but I think we can refine the design in Assembly

lindsayad · 2020-05-20T16:40:49Z

framework/include/constraints/MortarConstraintBase.h

@@ -151,6 +159,7 @@ class MortarConstraintBase : public Constraint,
  using ADMortarConstraint<compute_stage>::_phys_points_slave;                                     \
  using ADMortarConstraint<compute_stage>::_phys_points_master;                                    \
  using ADMortarConstraint<compute_stage>::_has_master;                                            \
+  using ADMortarConstraint<compute_stage>::_use_dual;                                              \


It looks like I missed removing this macro. compute_stage doesn't exist anymore. Can you remove this?

lindsayad · 2020-05-20T16:50:33Z

framework/src/base/Assembly.C

+    FEShapeData * fesd = _fe_shape_data_dual_lower[fe_type];
+
+    fe_lower->set_calculate_dual(true);


I ran into a segmentation fault the other day when trying to use a mixed-order discretization with u/order=SECOND and lm/order=FIRST and lm/use_dual=true. ~~I suspect this is what caused it.~~ (I actually don't think that's what caused it because I don't think you should have been creating any non-dual _fe_lowers... 🤔) You don't necessarily have a non-NULL _fe_shape_data_dual_lower for every _fe_lower do you? An _fe_lower is created anytime you call either buildLowerDDualFE or buildLowerDFE, but an _fe_shape_data_dual_lower is only created when you call buildLowerDDualFE. I can imagine an (admittedly very unlikely) scenario in which someone requests a standard LM and a dual LM of different orders and then you would be attempting to dereference a nullptr. I think that you should do something like this to dodge the issue:

if (!fesd) continue;

You don't necessarily have a non-NULL _fe_shape_data_dual_lower for every _fe_lower do you? An _fe_lower is created anytime you call either buildLowerDDualFE or buildLowerDFE, but an _fe_shape_data_dual_lower is only created when you call buildLowerDDualFE.

That is correct. The current implementation only gives non-NULL _fe_shape_data_dual_lower only when the user choose to use the dual method.

I ran into a segmentation fault the other day when trying to use a mixed-order discretization with u/order=SECOND and lm/order=FIRST and lm/use_dual=true.

That is right. The dual method requires lm/order = u/order use_dual=true, otherwise there would be an issue. Where do you suggest I may add a warning about this, just in case someone uses different order for LM and u?

Do we even need to use this API at all? If I had noticed this in the libMesh PR, I probably would've asked for it to be removed. There's no precedent for it in the FE hierarchy. We should know whether to calculate_dual based on the prerequests in buildLowerDDualFE. This is how I think a single reinitLowerDElemRef method should look:

void Assembly::reinitLowerDElemRef(const Elem * elem, const std::vector<Point> * const pts, const std::vector<Real> * const weights) { mooseAssert(pts->size(), "Currently reinitialization of lower d elements is only supported with custom " "quadrature points; there is no fall-back quadrature rule. Consequently make sure " "you never try to use JxW coming from a fe_lower object unless you are also passing " "a weights argument"); _current_lower_d_elem = elem; unsigned int elem_dim = elem->dim(); for (const auto & it : _fe_lower[elem_dim]) { FEBase * fe_lower = it.second; FEType fe_type = it.first; fe_lower->reinit(elem, pts, weights); if (FEShapeData * fesd = _fe_shape_data_lower[fe_type]) { fesd->_phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe_lower->get_phi())); fesd->_grad_phi.shallowCopy( const_cast<std::vector<std::vector<RealGradient>> &>(fe_lower->get_dphi())); if (_need_second_derivative_neighbor.find(fe_type) != _need_second_derivative_neighbor.end()) fesd->_second_phi.shallowCopy( const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_lower->get_d2phi())); } if (FEShapeData * fesd = _fe_shape_data_dual_lower[fe_type]) { fesd->_phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe_lower->get_dual_phi())); fesd->_grad_phi.shallowCopy( const_cast<std::vector<std::vector<RealGradient>> &>(fe_lower->get_dual_dphi())); if (_need_second_derivative_neighbor.find(fe_type) != _need_second_derivative_neighbor.end()) fesd->_second_phi.shallowCopy( const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_lower->get_dual_d2phi())); } } }

My concern was whether to reinitialize both dual and non-dual shape functions, as I imagine only one of them will be used in one simulation. If yes, then yeah I agree this would be a better way to combine both functions and remove this API

You need to compute the non-dual shape functions anyways in order to compute the dual shape functions. If you have a completely non-dual simulation then get_dual_phi will never get called, and fe_lower->reinit won't do any dual related computation. I don't see any scenario in the code I posted where you will be wasting work, but maybe I'm missing something...

I just checked after making the changes. You are right: for a completely non-dual simulation, get_dual_phi is never called. So no additional computation will be caused due to the change. Very helpful comments! Appreciate it!

lindsayad · 2020-05-20T17:04:24Z

framework/src/loops/ComputeMortarFunctor.C

+    // Whether to use dual mortar. Same for all motrat constraints
+    _use_dual = mc->use_dual();


Hmm why do we have to do this? Now that I'm thinking about it, I can actually fairly easily imagine non-dual and dual LM calculations within the same simulation and with different orders for the LMs. Any simulation with thermal and mechanical contact. Thermal contact does not have any saddle point issues and hence does not have a clear reason for using dual; moreover, you get optimal convergence in the primal variable if you use an order for the thermal LM of (primal_order - 1). I think that we need to figure out a way not to do what you're doing here. (Also _use_dual in ComputeMortarFunctor would just end up corresponding to whatever the last mortar constraint's _use_dual value was based on what you're doing here)

Hmm why do we have to do this? Now that I'm thinking about it, I can actually fairly easily imagine non-dual and dual LM calculations within the same simulation and with different orders for the LMs. Any simulation with thermal and mechanical contact. Thermal contact does not have any saddle point issues and hence does not have a clear reason for using dual; moreover, you get optimal convergence in the primal variable if you use an order for the thermal LM of (primal_order - 1). I think that we need to figure out a way not to do what you're doing here. (Also _use_dual in ComputeMortarFunctor would just end up corresponding to whatever the last mortar constraint's _use_dual value was based on what you're doing here)

The purpose was to pass in the use_dual input parameter in order to initialize the dual/no-dual lowerD elements accordingly. I have not met cases where different orders for the LMs exist for one simulation (and not sure if that would necessarily be a common practice). So I just left a comment there and assumed sameuse_dual value for all mortar constraints.

Not sure how to change it at this moment, but I probably can address issue while dealing with the ones you pointed out later in this PR.

I have not met cases where different orders for the LMs exist for one simulation (and not sure if that would necessarily be a common practice).

See #13080 (comment). If using first order displacements, you will want to use a non-dual CONSTANT MONOMIAL for the thermal contact Lagrange multiplier. If using second order displacements, you would want to use a LAGRANGE FIRST for the thermal contact LM. Based on what you told me, you will be wanting to use an order equal to the displacement order for the mechanical contact LM. In any realistic BISON simulation you will have both thermal and mechanical contact and so you will have mixing of dual and non-dual and multiple LM orders.

lindsayad · 2020-05-20T17:06:10Z

framework/src/loops/ComputeMortarFunctor.C

+    if (_use_dual)
+      _subproblem.reinitLowerDElemDualRef(slave_face_elem, &custom_xi1_pts);
+    else
+      _subproblem.reinitLowerDElemRef(slave_face_elem, &custom_xi1_pts);


It would be nice to just have a single reinitLowerDElemRef call that would initialize both "standard" and "dual" data structures. I feel like we could do that...

It would be nice to just have a single reinitLowerDElemRef call that would initialize both "standard" and "dual" data structures. I feel like we could do that...

(Uh sorry I did not entirely get your point earlier...) I am not sure if we will need to initialize both of them as we probably will (mostly) only be using one of them..

But I do not mind initializing both as that would be a easier implementation :)

lindsayad · 2020-05-20T17:07:07Z

framework/src/problems/SubProblem.C

+void
+SubProblem::reinitLowerDElemDualRef(const Elem * elem,
+                                    const std::vector<Point> * const pts,
+                                    const std::vector<Real> * const weights,
+                                    THREAD_ID tid)
+{
+  // - Set our _current_lower_d_elem for proper dof index getting in the moose variables
+  // - Reinitialize all of our lower-d FE objects so we have current phi, dphi, etc. data
+  assembly(tid).reinitLowerDElemDualRef(elem, pts, weights);
+
+  // Actually get the dof indices in the moose variables
+  systemBaseNonlinear().prepareLowerD(tid);
+  systemBaseAuxiliary().prepareLowerD(tid);
+
+  // With the dof indices set in the moose variables, now let's properly size
+  // our local residuals/Jacobians
+  assembly(tid).prepareLowerD();
+
+  // Let's finally compute our variable values!
+  systemBaseNonlinear().reinitLowerD(tid);
+  systemBaseAuxiliary().reinitLowerD(tid);
+}
+


Yea there's a lot of code duplication here. Again, I think we can design this such that a call to Assembly::reinitLowerDElemRef initializes both standard and dual data structures

moosebuild · 2020-05-20T19:36:00Z

Job Precheck on 0650de2 wanted to post the following:

Your code requires style changes.

A patch was auto generated and copied here
You can directly apply the patch by running, in the top level of your repository:

curl -s https://mooseframework.inl.gov/docs/PRs/15216/style.patch | git apply -v

Alternatively, with your repository up to date and in the top level of your repository:

git clang-format 6fd9c7f846a43385430795ac9e971716d2218d1c

…mputeMortarFunctor classes, issue idaholab#15215

lindsayad

Ah it looks so clean and beautiful! I love it. Great work!

andrsd · 2020-05-21T14:36:42Z

@dewenyushu Looks like your PR caused some valgrind problems. Would you mind to go see if you can fix them?

dewenyushu · 2020-05-21T14:43:14Z

@dewenyushu Looks like your PR caused some valgrind problems. Would you mind to go see if you can fix them?

Will do. Looking into the problem now.

dewenyushu · 2020-05-21T15:10:55Z

@dewenyushu Looks like your PR caused some valgrind problems. Would you mind to go see if you can fix them?

@andrsd Emm...I thought it would be an easy fix but it turns out that I am not familiar with the valgrind testing system and have no clue where this memory issue comes from... would you mind giving me a hint on this? Thanks!

dschwen · 2020-05-21T15:50:03Z

~~I think the issue was introduced in your libMesh PR!~~

andrsd · 2020-05-21T15:52:05Z

framework/src/base/Assembly.C

+Assembly::buildLowerDDualFE(FEType type) const
+{
+  if (!_fe_shape_data_dual_lower[type])
+    _fe_shape_data_dual_lower[type] = new FEShapeData;


This needs to be freed in Assembly::~Assembly().

Can't we rather use smart pointers?

nice catch!

andrsd · 2020-05-21T15:55:56Z

framework/src/base/Assembly.C

+Assembly::buildVectorDualLowerDFE(FEType type) const
+{
+  if (!_vector_fe_shape_data_dual_lower[type])
+    _vector_fe_shape_data_dual_lower[type] = new VectorFEShapeData;


And also this one.

dewenyushu · 2020-05-21T15:56:22Z

I think the issue was introduced in your libMesh PR!

Thanks for trying to help!

dewenyushu force-pushed the dual_mortar_15215 branch from 742caa2 to c6b76ff Compare May 4, 2020 21:35

dewenyushu changed the title ~~Enable usage of dual shape functions for the Lagrange multipliers~~ WIP: Enable usage of dual shape functions for the Lagrange multipliers May 4, 2020

dewenyushu force-pushed the dual_mortar_15215 branch from c6b76ff to 3eca7e2 Compare May 4, 2020 21:44

dewenyushu changed the title ~~WIP: Enable usage of dual shape functions for the Lagrange multipliers~~ Enable usage of dual shape functions for the Lagrange multipliers May 5, 2020

dewenyushu commented May 11, 2020

View reviewed changes

dewenyushu mentioned this pull request May 11, 2020

Add dual shape functions libMesh/libmesh#2552

Merged

dewenyushu force-pushed the dual_mortar_15215 branch from 3eca7e2 to b3a9dbe Compare May 13, 2020 18:21

dewenyushu force-pushed the dual_mortar_15215 branch from b3a9dbe to d10f8fc Compare May 18, 2020 18:28

dewenyushu closed this May 19, 2020

dewenyushu reopened this May 19, 2020

Enable 'use_dual' for the Lagrange multipliers

0785347

Get conforming and nonconforming cases to converge, checked Jacobian pattern Fix broken tests, code restructure, and add tests Address issue idaholab#15215

dewenyushu mentioned this pull request May 19, 2020

Dual mortar #15215

Open

4 tasks

Add test for non-conforming mesh using 2nd order elements

96d5ef8

Remove test for the conforming mesh Rebase against devel Address issue idaholab#15215

dewenyushu force-pushed the dual_mortar_15215 branch from d10f8fc to 96d5ef8 Compare May 19, 2020 22:33

moosebuild added the PR: Failed but allowed label May 20, 2020

lindsayad requested changes May 20, 2020

View reviewed changes

moosebuild removed the PR: Failed but allowed label May 20, 2020

Code simplification and restructuring in Assembly, SubProblem, and Co…

be92a2c

…mputeMortarFunctor classes, issue idaholab#15215

dewenyushu force-pushed the dual_mortar_15215 branch from 0650de2 to be92a2c Compare May 20, 2020 19:38

dewenyushu requested a review from lindsayad May 20, 2020 22:49

moosebuild added the PR: Failed but allowed label May 20, 2020

lindsayad approved these changes May 20, 2020

View reviewed changes

lindsayad merged commit 8e2c8ac into idaholab:next May 20, 2020

andrsd reviewed May 21, 2020

View reviewed changes

This was referenced May 21, 2020

Fix PR 15216 #15332

Closed

Fix PR 15216 #15336

Merged

		FEShapeData * fesd = _fe_shape_data_dual_lower[fe_type];

		fe_lower->set_calculate_dual(true);

		// Whether to use dual mortar. Same for all motrat constraints
		_use_dual = mc->use_dual();

Enable usage of dual shape functions for the Lagrange multipliers #15216

Enable usage of dual shape functions for the Lagrange multipliers #15216

Conversation

dewenyushu commented May 4, 2020

dewenyushu commented May 4, 2020

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

lindsayad May 11, 2020 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dewenyushu May 11, 2020 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dewenyushu May 11, 2020 • edited Loading

Choose a reason for hiding this comment

dewenyushu commented May 13, 2020

lindsayad commented May 13, 2020

dewenyushu commented May 13, 2020

lindsayad commented May 13, 2020 via email

moosebuild commented May 18, 2020 • edited Loading

dewenyushu commented May 18, 2020

lindsayad commented May 19, 2020 • edited Loading

lindsayad commented May 19, 2020

lindsayad commented May 19, 2020

jwpeterson commented May 19, 2020 • edited Loading

dewenyushu commented May 19, 2020

jwpeterson commented May 19, 2020

dewenyushu commented May 19, 2020

dewenyushu commented May 19, 2020

jwpeterson commented May 19, 2020

dewenyushu commented May 19, 2020 • edited Loading

dewenyushu commented May 19, 2020 • edited Loading

jwpeterson commented May 19, 2020

lindsayad commented May 19, 2020

dewenyushu commented May 19, 2020

lindsayad commented May 19, 2020

jwpeterson commented May 19, 2020

dewenyushu commented May 19, 2020

dewenyushu commented May 20, 2020

lindsayad left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dewenyushu May 20, 2020 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dewenyushu May 20, 2020 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dewenyushu May 20, 2020 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

moosebuild commented May 20, 2020

lindsayad left a comment

Choose a reason for hiding this comment

andrsd commented May 21, 2020

dewenyushu commented May 21, 2020

dewenyushu commented May 21, 2020 • edited Loading

dschwen commented May 21, 2020 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dewenyushu commented May 21, 2020

lindsayad May 11, 2020 •

edited

Loading

dewenyushu May 11, 2020 •

edited

Loading

dewenyushu May 11, 2020 •

edited

Loading

moosebuild commented May 18, 2020 •

edited

Loading

lindsayad commented May 19, 2020 •

edited

Loading

jwpeterson commented May 19, 2020 •

edited

Loading

dewenyushu commented May 19, 2020 •

edited

Loading

dewenyushu commented May 19, 2020 •

edited

Loading

dewenyushu May 20, 2020 •

edited

Loading

dewenyushu May 20, 2020 •

edited

Loading

dewenyushu May 20, 2020 •

edited

Loading

dewenyushu commented May 21, 2020 •

edited

Loading

dschwen commented May 21, 2020 •

edited

Loading