Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

LCM: Convergence of Schwarz coupling does not match independent implementation. #53

Closed
lxmota opened this issue Feb 2, 2017 · 11 comments
Assignees
Labels
Projects

Comments

@lxmota
Copy link
Contributor

lxmota commented Feb 2, 2017

The Modified Schwarz implementation in Albany displays oscillatory and slow convergence compared to a Matlab implementation that shows monotonic and fast convergence. The root cause is not known.

@lxmota lxmota self-assigned this Feb 2, 2017
@lxmota
Copy link
Contributor Author

lxmota commented Feb 8, 2017

@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
Copy link
Contributor

calleman21 commented Feb 9, 2017

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
Copy link
Contributor Author

lxmota commented Feb 9, 2017

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
Copy link
Contributor Author

lxmota commented Feb 9, 2017

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 added a commit that referenced this issue Feb 10, 2017
and its time derivaties in the application for convenient
retrieval. (#53)
lxmota added a commit that referenced this issue Feb 21, 2017
@lxmota
Copy link
Contributor Author

lxmota commented Feb 21, 2017

@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
Copy link
Contributor

jwfoulk commented Feb 21, 2017

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
Copy link
Contributor Author

lxmota commented Feb 21, 2017

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

@jwfoulk
Copy link
Contributor

jwfoulk commented Feb 21, 2017

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
Copy link
Contributor Author

lxmota commented Feb 22, 2017

@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
Copy link
Collaborator

ikalash commented Feb 22, 2017

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
Copy link
Contributor Author

lxmota commented Feb 22, 2017

@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.

lxmota added a commit that referenced this issue Mar 8, 2017
lxmota added a commit that referenced this issue Mar 8, 2017
lxmota added a commit that referenced this issue Mar 8, 2017
lxmota added a commit that referenced this issue Mar 8, 2017
Remove T suffix to denote Tpetra from variable and functions names
as it clutters notation and is unnecessary here. (#22, #75)
lxmota added a commit that referenced this issue Mar 8, 2017
alternating version of Schwarz. It seems that continuation
will be done using NOX's Anderson acceleration instead of LOCA
for this case, so this feature will need to be reworked once
continuation is done. (#53)
lxmota added a commit that referenced this issue Mar 9, 2017
lxmota added a commit that referenced this issue Mar 12, 2017
to be an easier path to support continuation. Remove left-over
code from Schwarz_Coupled. Enable acceleration support by default
for later implementation of dynamics. (#53)
lxmota added a commit that referenced this issue Mar 14, 2017
how to implement the convergence criterion of the Schwarz loop
that matches what is in the paper and the Matlab code. (#53)
lxmota added a commit that referenced this issue Mar 14, 2017
lxmota added a commit that referenced this issue Mar 15, 2017
for validating parameter list. Remove unneeded code. (#53)
lxmota added a commit that referenced this issue Mar 16, 2017
lxmota added a commit that referenced this issue Mar 16, 2017
Remove unneeded member functions. Use clang-format as we touch other
source files.

[#22, #38, #53]
lxmota added a commit that referenced this issue Mar 16, 2017
apparently not useful because Albany::ModelEvaluator complains
regarldess that solution is null. [#53]

Run through clang-format and add \dotdot{x} support now that
it exists in Thyra. [#22, #38]
lxmota added a commit that referenced this issue Mar 17, 2017
in_args instead of application by unpacking the ugly blob
that comes in in_args. Still get null pointer. [#53]
lxmota added a commit that referenced this issue Mar 20, 2017
lxmota added a commit that referenced this issue Mar 21, 2017
is due to evalModel for each subdomain called from within evalModel
of the controlling Schwarz loop. In a normal solve, NOX sets up
the Tpetra vectors and stuffs them into inArgs before calling evalModel.
This is missing here because the controlling Schwarz loop is not called
from within NOX. Trying to fix this. Slight change to input syntax to
reflect the above. [#53]
lxmota added a commit that referenced this issue Mar 22, 2017
just models in SchwarzAlternating class. Rework the class and Schwarz
loop accordingly. This version appears to work in performing full Schwarz
on the cuboids problem. Needs much more testing. Issues:
- Work on convergence criterion. Right now is hard-coded on maximum
number of iterations.
- Read convergence criterion from YAML/XML input.
- The algorithm writes the results of each Schwarz iteration to the
exodus output, which is useful for debugging but not for production.
- Main_SolveT expects responses for the main solver, but Schwarz
alternating has none, so work around that.
- Runs to completion, but Albany complains in the end of RCP errors.
May be related to responses. Need to debug.

[#53]
lxmota added a commit that referenced this issue Mar 23, 2017
number of Schwarz iterations, absolute and relative tolerances
for the Schwarz loop. [#53]
lxmota added a commit that referenced this issue Mar 23, 2017
number of Schwarz iterations, absolute and relative tolerances
for the Schwarz loop. [#53]
@ibaned ibaned added this to Backlog in Albany Apr 13, 2017
lxmota added a commit that referenced this issue Apr 19, 2017
PrePostOperator. Fix handling of parameter lists. [#53]
lxmota added a commit that referenced this issue Apr 27, 2017
to handle also the Schwarz termination criterion.

It appears that there can be only one NOX pre-post operator, thus
there was a conflict between the model evaluator status flag operator
and the Schwarz termination criterion operator.

So, for now, add the functionality needed for Schwarz to the status
flag. This requires performing some operations to compute norms of
the solution at various times and adding some storage to keep these
values for later use in the Schwarz loop. Thus, this does not affect
a non-Schwarz analysis.

Still, it is probably a good idea to come up with a plan to handle
more than one NOX pre-post operator, so that this situation does
not repeat itself. And also to avoid out-of-control growth of
this class if more pre-post functionality is needed in the future.
lxmota added a commit that referenced this issue Apr 27, 2017
lxmota added a commit that referenced this issue Apr 27, 2017
functionality has been folded into the NOXSolverPrePostOperator class.
lxmota added a commit that referenced this issue Apr 28, 2017
Schwarz paper. Issues fixed:

- Conflict between Schwarz convergence criterion and status test
NOX pre/post operators resolved with folding both into the same
class. In addition, when trying to create one or the other, check
for the existence of a pre/post operator. If one exists already,
skip the creation and use the existing one as it supports both the
status test and Schwarz convergence criterion. Still need to discuss
at large what to do in the long term to avoid out-of-control
bloat in case more functionality is needed.

- The NOX solver resets the solution to zeros before each non-linear
solve, but keeps a copy of the old solution in a previous solution
group. Use this group and the current group to compute the norms that
define the Schwarz convergence criterion. Store this info in
the pre/post operator class.

- Keep an array of the pre/post operators for each subdomain when
the solvers are created. These operators are the only way to access
the solutions and their corresponding norms in the Schwarz loop
that comes much later.
lxmota added a commit that referenced this issue Apr 28, 2017
alternating method:

- Signal that this is a Schwarz alternating method run.
- Track each Schwarz iteration.
- Track solution method for each subdomain.
- Print norms for each subdomain and composite norms that
determine the Schwarz convergence criterion after each Schwarz iteration.
- Signal when the Schwarz iteration has converged.
lxmota added a commit that referenced this issue May 1, 2017
lxmota added a commit that referenced this issue May 2, 2017
Things to do:

- Find a way to avoid exodus output at every Schwarz iteration.
This requires taking control of the output away from the nonlinear
solver and give it to the Schwarz method.

- Continuation. It appears that doing this through LOCA would be
too difficult. A suggested way to do this is through Anderson
acceleration, a variant of fixed point methods. But now that Tempus
is enabled, perhaps it would be easier to do continuation through
dynamics. Then quasi-static problems would be run as truly
quasi-static: very slow dynamic problems.

Sice we now the root cause of the unmatched convergence for the
old Schwarz implementation, I'm closing this issue and opening
new issues for the above.
@lxmota lxmota closed this as completed May 2, 2017
@lxmota lxmota moved this from Backlog to Done in Albany May 2, 2017
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
Albany
  
Done
Development

No branches or pull requests

4 participants