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
Compute a better alpha factor in the implementation of the SPD-ness of the Newton method. #5555
Comments
Thank @bangerth for raising this issue. It may not be a good time to discuss about it at this busy moment, but I would be grateful if someone can take a look and check my derivations in his spare time. My point is that the top-left block of the Jacobian matrix for most of the homogeneous rheology models we use are non-negative definite. The detailed derivations are presented in the following 4-page pdf: If my conclusions are correct, then we can simply scale the Newton step by a user-specified number I have implemented a test code in https://github.com/YiminJin/aspect/tree/newton-solver. In this branch, I set The reason that the Newton solver is unstable without stabilization in the original code is complex. To make a proper modification, one may have to change the structure of MaterialModel. I think it would be more efficient to have an online meeting about this issue when most of us are free. |
I found some mistakes in the pdf. Here is the revised version: |
For reference, the analysis in @YiminJin's paper only applies if you have that |
Also for reference, if we wanted to address this, the function that does the work in question is in
|
I agree with @bangerth that my conclusion only applies to homogeneous models with moderate strain weakening. If the strain weakening is extremely strong (for example, the stress exponent in the power-law model is The biggest trouble in modifying the code is in the elastic rheology. For viscoelastic models, we have |
Do I understand correctly that this was fixed by the linked PR? If so, feel free to close this issue. |
@YiminJin pointed out in an email to @MFraters that we have a mistake in the paper that describes the Newton method. Specifically, in the derivation of eq. (18), we introduce a vector
that for different values must include the two eigenvectors of the operator
E^sym = a \otimes b + b \otimes a
. But this is not correct, as he points out: The formula forv
is only true ifa
andb
are parallel.For vectors (actually, symmetric tensors)
a
,b
that are not perpendicular, one should replaceb
byc = b - b:a/(||a|| ||b||)
, i.e., the component ofb
perpendicular toa
. Assumingc
is not (numerically) zero, the formulas below remain true. Forc=0
, i.e., ifb
is parallel toa
,E^sym
has one zero eigenvalue and one that is easily computable as2 a:b
(which can be both positive and negative).All of this could probably be used to compute an
alpha
factor in the Newton implementation that is closer to one and that perhaps leads to faster convergence.The text was updated successfully, but these errors were encountered: