<div dir="ltr"><div>Hi Barry,</div><div><br></div><div>I believe that what you are referring to is what Jed is referring to in this thread, am I right?</div><div></div><div><a href="https://scicomp.stackexchange.com/questions/3298/appropriate-space-for-weak-solutions-to-an-elliptical-pde-with-mixed-inhomogeneo/3300#3300">https://scicomp.stackexchange.com/questions/3298/appropriate-space-for-weak-solutions-to-an-elliptical-pde-with-mixed-inhomogeneo/3300#3300</a></div><div><br></div><div>We do set the rows/cols of the Jacobian to zero except the diagonal component which is set to one, as you mention. I understand that in the case of only homogeneous Dirichlet BCs it is generally a good idea to scale that diagonal component so that the condition number of the Jacobian improves. I assume that what Jed mentions is the inhomogeneous Dirichlet BC version of this scaling, which also acts on the corresponding indices of the residual, not just the Jacobian. My question is the following, since the case we are encountering problems with is a system with only homogeneous Dirichlet BCs, how does it apply? Also, would this scaling affect the convergence of the NEWTONLS with "bt" line-search? Without any scaling we can solve this example with "basic" (with a very reasonable convergence rate), but not with "bt" line-search.</div><div><br></div><div>Regards,</div><div>Francesc.<br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Tue, Nov 16, 2021 at 6:55 PM Barry Smith <<a href="mailto:bsmith@petsc.dev">bsmith@petsc.dev</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div style="overflow-wrap: break-word;"><div><br></div>   By "scaling" I mean you may need to rescale the residual terms from the Dirichlet boundary conditions to be the same  order  (in terms of the characteristic size  of the elements) as they are for the other residuals.<div><br></div><div>  It seems in your formulation the Jacobian entries of the residual for Dirichlet boundary conditions are 1 on the diagonal while the Jacobian entries of the residual for the PDE entries will have a scaling related to the finite element method you are using, for example in 3d with an operator grad dot grad the Jacobian entry will be order dx = (dx)^3 * (1/dx)^2. Thus you should scale the Dirichlet boundary condition residuals by dx to get the same scaling.</div><div><br></div><div>  Barry</div><div><br><div><br><blockquote type="cite"><div>On Nov 16, 2021, at 1:37 PM, Francesc Levrero-Florencio <<a href="mailto:f.levrero-florencio@onscale.com" target="_blank">f.levrero-florencio@onscale.com</a>> wrote:</div><br><div><div dir="ltr"><div>Hi Barry,</div><div><br></div><div>Thanks for the quick reply. I am not sure what you meant by "scaling". In our code, every time our residual is called we do the following in this order:</div><div>- Apply updated Dirichlet BCs in order to assemble the internal forces, and add them to the residual.<br></div><div>- Assemble the external forces, and subtract them from the residual.</div><div>- Update the residual to take into account Dirichet BCs, so we set the corresponding indices of the residual to - (dirichlet_bc_value - solution_value), the latter obtained directly from the solution vector of PETSc. However, note that in this case all Dirichlet BCs are homogeneous, so this value would be just "solution_value", or 0 throughout the simulation. This ensures that the solution vector coming from PETSc has the Dirichlet BCs applied.<br></div><div><br></div><div>This is basically what we do in the residual. If I am not mistaken, the lambda scaling factor from the line-search would scale the factor "jacobian^{-1} * residual", so it would scale equally both internal and external forces.</div><div><br></div><div>Regards,</div><div>Francesc.<br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Tue, Nov 16, 2021 at 5:38 PM Barry Smith <<a href="mailto:bsmith@petsc.dev" target="_blank">bsmith@petsc.dev</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><br>
  Perhaps this behavior is the result of a "scaling" issue in how various terms affect the residual? In particular perhaps the terms for enforcing boundary conditions are scaled differently than terms for the PDE enforcement?<br>
<br>
<br>
<br>
> On Nov 16, 2021, at 11:19 AM, Francesc Levrero-Florencio <<a href="mailto:f.levrero-florencio@onscale.com" target="_blank">f.levrero-florencio@onscale.com</a>> wrote:<br>
> <br>
> Dear PETSc team and users,<br>
> <br>
> We are running a simple cantilever beam bending, where the profile of the beam is I-shaped, where we apply a bending force on one end and fully constrained displacements on the other end. The formulation is a large strain formulation in Total Lagrangian form, where the material of the beam is a Saint Venant-Kirchhoff hyperelastic material that uses the same constants as steel (200E9 GPa Young’s modulus and 0.3 Poisson’s ratio).<br>
> <br>
> The simulation finishes in the requested number of time steps by using the “basic” type of line-search in the SNES (i.e. traditional Newton method without line-search) in a reasonable number of Newton iterations per time step (3 or 4 iterations). However, when using the “bt” (or “l2”, and in fact no type of line-search ends up yielding convergence) line-search type, the convergence never happens within the SNES default maximum number of iterations of 50.<br>
> <br>
> During solving with traditional Newton, the general behaviour of each time step is that the norm of the residual increases on the second call to the residual function, but then hugely decreases in the following one or two, up to the point where convergence is achieved. Using “bt” line-search, the line-search discards the step at lambda=1 immediately because the norm of the residual is higher than that produced in the first call to the residual function, cutting down the value of lambda to a value significantly lower than 1. The simulation then progresses in following Newton iterations in a similar fashion, the line-search step at lambda=1 is always discarded, and then smaller steps are taken but convergence never occurs, for even the first time step.<br>
> <br>
> Here are a few values of the relevant norms (using traditional Newton) in the first time step:<br>
> <br>
> BASIC NEWTON LS<br>
> Norm of the internal forces is 0<br>
> Norm of the external forces is 1374.49<br>
> Norm of the residual is 1374.49<br>
> Norm of the solution with Dirichlet BCs is 0<br>
> Number of SNES iteration is 0<br>
> ---------------------------------------------------------------------<br>
> Norm of the internal forces is 113498<br>
> Norm of the external forces is 1374.49<br>
> Norm of the residual is 105053<br>
> Norm of the solution with Dirichlet BCs is 0.441466<br>
> Number of SNES iteration is 0<br>
> ---------------------------------------------------------------------<br>
> Norm of the internal forces is 42953.5<br>
> Norm of the external forces is 1374.49<br>
> Norm of the residual is 11.3734<br>
> Norm of the solution with Dirichlet BCs is 0.441438<br>
> Number of SNES iteration is 1<br>
> <br>
> Here are a few values of the relevant norms (using “bt” line-search) in the first time step:<br>
> <br>
> BT NEWTONLS<br>
> Norm of the internal forces is 0<br>
> Norm of the external forces is 1374.49<br>
> Norm of the residual is 1374.49<br>
> Norm of the solution with Dirichlet BCs is 0<br>
> Number of SNES iteration is 0<br>
> ---------------------------------------------------------------------<br>
> Norm of the internal forces is 113498<br>
> Norm of the external forces is 1374.49<br>
> Norm of the residual is 105053<br>
> Norm of the solution with Dirichlet BCs is 0.441466<br>
> Number of SNES iteration is 0<br>
> (I assume this is the first try at lambda=1)<br>
> ---------------------------------------------------------------------<br>
> Norm of the internal forces is 4422.12<br>
> Norm of the external forces is 1374.49<br>
> Norm of the residual is 1622.74<br>
> Norm of the solution with Dirichlet BCs is 0.0441466<br>
> Number of SNES iteration is 0<br>
> Line search: gnorm after quadratic fit 1.622742343614e+03<br>
> (I assume that in this line-search iteration 0.05 < lambda < 1, but the corresponding residual is not smaller than the one in the first call, 1622 > 1374)<br>
> ---------------------------------------------------------------------<br>
> Norm of the internal forces is 2163.76<br>
> Norm of the external forces is 1374.49<br>
> Norm of the residual is 1331.88<br>
> Norm of the solution with Dirichlet BCs is 0.0220733<br>
> Number of SNES iteration is 0<br>
> Line search: Cubically determined step, current gnorm 1.331884625811e+03 lambda=5.0000000000000003e-02<br>
> (This is the accepted lambda for the current Newton iteration)<br>
> ---------------------------------------------------------------------<br>
> Norm of the internal forces is 104020<br>
> Norm of the external forces is 1374.49<br>
> Norm of the residual is 94739<br>
> Norm of the solution with Dirichlet BCs is 0.441323<br>
> Number of SNES iteration is 1<br>
> <br>
> My question would be, any idea on how to deal with this situation? I would imagine a “hack” would be to bypass the first residual norm, and have the line-search use the following one as the “base norm” to check its reduction in further iterations, but we are open to any ideas.<br>
> <br>
> Thanks for your help in advance and please keep up the good work!<br>
> <br>
> Regards,<br>
> Francesc.<br>
<br>
</blockquote></div>
</div></blockquote></div><br></div></div></blockquote></div>