LCM: Convergence of Schwarz coupling does not match independent implementation. about albany

[lxmota]

@ikalash @jwfoulk @calleman21

I see that every time that NOX updates the solution, it calls for a residual evaluation. Then inside SchwarzMultiscale::evalModelImpl we force the writing of the solution to the STK mesh database. The Schwarz BC then has access to the latest solution and it should all be good.

But what happens is that when I compare the residuals produced by YuckySchwarz and Albany, they differ. The residual entries for the Schwarz BCs are the difference between the current solution and the projected solution from the other subdomain. What happens is that in Albany that projection is from the previous iteration, not from the current.

[calleman21]

I have one data point that might be relevant: when I change
<Parameter name="Forcing Term Method" type="string" value="Constant"/>
to
<Parameter name="Forcing Term Method" type="string" value="Type 2"/>
in the NOX->Direction->Newton sublist, the oscillatory behavior is changed, with less oscillations after the first couple of iterations. Also, the residual norm seems to eventually reduce to the limits of numerical precision. @agsalin, do you know what this option does? Does it force re-evaluation of the residual at some additional point in the computation?

[lxmota]

I just tried
<Parameter name="Forcing Term Method" type="string" value="Type 2"/>
on the test problem that I'm using for debugging this.

While previously I had monotonic (although slow) convergence, now the residual oscillates. It also increased the number of iterations significantly.

[lxmota]

From NOX_Description.H:

• "Forcing Term Method" - Method to compute the forcing term, i.e.,
  the tolerance for the linear solver. Defaults to ""
  (nothing). Choices are "Type 1" and "Type 2".

• "Forcing Term Minimum Tolerance" - Minimum acceptable linear
  solver tolerance. Defaults to 1.0e-6.

• "Forcing Term Maximum Tolerance" = Maximum acceptable linear
  solver tolerance. Default to 0.01.

• "Forcing Term Alpha" - Used for the "Type 2" forcing term
  calculation. Defaults to 1.5.

• "Forcing Term Gamma" - Used for the "Type 2" forcing term
  calculation. Defaults to 0.9.

[lxmota]

@ikalash @calleman21 @jtostie @jwfoulk

I have determined the reason why the Albany and Matlab implementations behave differently.

In essence, it is as the difference between the Gauss-Seidel and Jacobi methods for linear systems.

In Albany, the block system is formed and only then an increment for the displacement is calculated. The residuals contain the interpolated displacements between subdomains from the previous iterations. This is similar to the Gauss-Seidel method.

In Matlab, the linear systems are solver sequentially. The first system to be solved contains in its residual interpolated displacements from the other subdomain from the previous iteration. Its solution is then incremented. Then the second system is solved, and its residual contains interpolated displacements from the first subdomain from current iteration. This is similar to the Jacobi method.

This makes the Matlab code converge faster.

No preconditioner helps here, because before any preconditioner is applied, the residuals have already been calculated.

The next step is to figure out how difficult would be to solve the systems sequentially in Albany.

[jwfoulk]

Hi Alejandro - sorry to keep asking this question - please remind me of the matrix-free structure such that it avoids the lag in the interpolated displacements. Thanks again.

[lxmota]

@jwfoulk Sorry Jay, I'm uncertain about your question. Could you elaborate?

[jwfoulk]

Sorry about my brevity. First and foremost - great news in terms of understanding the issue. My questions just stems from the fact that the matrix-free methods do not seem to lag and consequently are probably using the current iteration. I'm just wondering if this indeed true and might this provide some path forward for the "modified" block system given that the matrix-free methods solve a block system.

[lxmota]

@jwfoulk I see. I guess then that the question is whether the matrix-free (or monolithic Schwarz in our parlance) would converge faster sequentially. I would need to implement that in Matlab to answer that question.

[ikalash]

Sorry for the delay in my response. @lxmota : so, am I understanding correctly that for modified Schwarz, we cannot form/solve a block diagonal system as we are doing now; we'd need to solve the system on omega1 using a regular (non-block) linear solver, then solve the system on omega2 (using the solution from omega1) after the system on omega1 is solved using a regular (non-block) linear solver, all within a Newton-Schwarz step?

It would be interesting to see what happens for the matrix-free (monolithic Schwarz) if solved sequentially, but for that method, it is not as critical, since the method already has some of the properties the sequential modified Schwarz has -- e.g., monotonic convergence. I guess things get coupled more tightly than in the modified Schwarz implicitly through the off-diagonal block terms somehow.

[lxmota]

@ikalash I guess you can form the block system for Modified Schwarz, but the price you pay is slower convergence, as we see now in the Albany implementation.

Yes, I have been thinking about whether it is possible to have a sequential version of the monolithic Schwarz. Not sure it is possible. But as you say, at the moment it is not as critical.