Re: [petsc-users] TSBEULER vs TSPSEUDO
Hi Jed, Thanks for the answer. We do have a monolithic arc-length implementation based on the TS/SNES logic, but we are also exploring having a custom SNESSHELL because the arc-length logic is substantially more complex than that of traditional load-controlled continuation methods. It works quite well, the only "issue" is its initiation; we are currently performing load-control (or displacement loading as you mentioned) in the first time increment. Besides load-control and arc-length control, what other continuation methods would you suggest exploring? The test problem we are dealing with assumes plasticity but with small strains so we will not see any snap-throughs, snap-backs or similar. TSBEULER works quite well for this specific case and converges in a few time steps within around 5-10 SNES iterations per time step. What PETSc functions do you suggest exploring for implementing the TS time step extension control you mentioned? Since you mentioned -ts_theta_initial_guess_extrapolate, is it worth using it in highly nonlinear mechanical problems (such as plasticity)? It sounds quite useful if it consistently reduces SNES iterations by one per time step, as each linear solve is quite expensive for large problems. Regards, Francesc. From: Jed Brown Sent: 08 November 2022 17:09 To: Francesc Levrero-Florencio ; petsc-users@mcs.anl.gov Subject: Re: [petsc-users] TSBEULER vs TSPSEUDO [External Sender] First, I believe arc-length continuation is the right approach in this problem domain. I have a branch starting an implementation, but need to revisit it in light of some feedback (and time has been too short lately). My group's nonlinear mechanics solver uses TSBEULER because it's convenient to parametrize loading on T=[0,1]. Unlike arc-length continuation, this can't handle snap-through effects. TSPSEUDO is the usual recommendation if you don't care about time accuracy, though you could register a custom controller for normal TS methods that implements any logic you'd like around automatically extending the time step without using a truncation error estimate. Note that displacement loading (as usually implemented) is really bad (especially for models with plasticity) because increments that are large relative to the mesh size can invert elements or initiate plastic yielding when that would not happen if using smaller increments. Arc-length continuation also helps fix that problem. Note that you can use extrapolation (-ts_theta_initial_guess_extrapolate), though I've found this to be somewhat brittle and only reduce SNES iteration count by about 1 per time step. Francesc Levrero-Florencio writes: > Hi PETSc people, > > We are running highly nonlinear quasi-static (steady-state) mechanical finite > element problems with PETSc, currently using TSBEULER and the basic time > adapt scheme. > > What we do in order to tackle these nonlinear problems is to parametrize the > applied loads with the time in the TS and apply them incrementally. While > this usually works well, we have seen instances in which the adaptor would > reject the time step according to the calculated truncation errors, even if > the SNES converges in a small number of iterations. Another issue that we > have recently observed is that in a sequence of converged time steps the > adaptor decides to start cutting the time step to smaller and smaller values > using the low clip default value of TSAdaptGetClip (again because the > truncation errors are high enough). What can we do in order to avoid these > issues? The first one is avoided by using TSAdaptSetAlwaysAccept, but the > latter remains. We have tried setting the low clip value to its maximum > accepted value of 1, but then the time increment does not increase even if > the SNES always converges in 3 or 4 iterations. Maybe a solution is to > increase the tolerances of the TSAdapt? > > Another potential solution we have recently tried in order to tackle these > issues is using TSPSEUDO (and deparametrizing the applied loads), but > generally find that it takes a much longer time to reach an acceptable > solution compared with TSBEULER. We have mostly used the default KSPONLY > option, but we'd like to explore TSPSEUDO with NEWTONLS. A first question > would be: what happens if the SNES fails to converge, does the solution get > updated somehow in the corresponding time step? We have performed a few tests > with TSPSEUDO and NEWTONLS, setting the maximum number of SNES iterations to > a relatively low number (e.g. 5), and then always setting the SNES as > converged in the poststage function, and found that it performs reasonably > well, at least better than with the default KSPONLY (does this make any > sense?). > > Thanks a lot! > > Regards, > Francesc.
[petsc-users] TSBEULER vs TSPSEUDO
Hi PETSc people, We are running highly nonlinear quasi-static (steady-state) mechanical finite element problems with PETSc, currently using TSBEULER and the basic time adapt scheme. What we do in order to tackle these nonlinear problems is to parametrize the applied loads with the time in the TS and apply them incrementally. While this usually works well, we have seen instances in which the adaptor would reject the time step according to the calculated truncation errors, even if the SNES converges in a small number of iterations. Another issue that we have recently observed is that in a sequence of converged time steps the adaptor decides to start cutting the time step to smaller and smaller values using the low clip default value of TSAdaptGetClip (again because the truncation errors are high enough). What can we do in order to avoid these issues? The first one is avoided by using TSAdaptSetAlwaysAccept, but the latter remains. We have tried setting the low clip value to its maximum accepted value of 1, but then the time increment does not increase even if the SNES always converges in 3 or 4 iterations. Maybe a solution is to increase the tolerances of the TSAdapt? Another potential solution we have recently tried in order to tackle these issues is using TSPSEUDO (and deparametrizing the applied loads), but generally find that it takes a much longer time to reach an acceptable solution compared with TSBEULER. We have mostly used the default KSPONLY option, but we'd like to explore TSPSEUDO with NEWTONLS. A first question would be: what happens if the SNES fails to converge, does the solution get updated somehow in the corresponding time step? We have performed a few tests with TSPSEUDO and NEWTONLS, setting the maximum number of SNES iterations to a relatively low number (e.g. 5), and then always setting the SNES as converged in the poststage function, and found that it performs reasonably well, at least better than with the default KSPONLY (does this make any sense?). Thanks a lot! Regards, Francesc.
Re: [petsc-users] Doubt about BT and BASIC NEWTONLS
Hi Jed,
Thanks for the information.
Yes, we set the rows and cols of the matrix for essential BCs to zero
(except diagonal, to 1). For homogeneous BCs it would be exactly as the
first term you mentioned in the link, and for nonhomogeneous BCs it would
be both terms, but the second one not scale (i.e. alpha set to 1).
I have tried a few of the line searches, "bt" does not work, but "cp" does
progress for a few time steps, with reasonable numbers of Newton iterations
(ends up stalling though). It does look to me that the "bt" line search
seems to "get stuck" on a local minima, while the full Newton, maybe due to
the size of the step, gets to overcome it. I am just a bit "surprised"
because in my experience line searches, particularly of the "bt" kind, tend
to improve, or at least not decrease, the convergence behaviour of the
Newton solver.
Regards,
Francesc.
On Wed, Nov 17, 2021 at 5:49 PM Jed Brown wrote:
> Francesc Levrero-Florencio writes:
>
> > Hi Barry,
> >
> > I believe that what you are referring to is what Jed is referring to in
> > this thread, am I right?
> >
> https://scicomp.stackexchange.com/questions/3298/appropriate-space-for-weak-solutions-to-an-elliptical-pde-with-mixed-inhomogeneo/3300#3300
>
> Yeah, that's the scaling. Are you decoupling the interior in the way I
> described, so the matrix columns for essential boundary conditions are also
> zeroed?
>
> Also note that line searches can prevent a rootfinding method from
> converging, as in this example.
>
> https://scicomp.stackexchange.com/a/2446/119
>
> There is -snes_linesearch_type cp ("critical point"), which has a
> surrogate that looks like aWolfe conditions when your rootfinding problem
> happens to be the first order optimality conditions for a minimization
> problem. There's also SNESSetObjective(), if your problem has an explicit
> objective. In practice, cp usually works well if your problem is "almost"
> coming from a minimization principle, and poorly otherwise.
>
> > 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.
>
>
Re: [petsc-users] Doubt about BT and BASIC NEWTONLS
Hi Barry,
I believe that what you are referring to is what Jed is referring to in
this thread, am I right?
https://scicomp.stackexchange.com/questions/3298/appropriate-space-for-weak-solutions-to-an-elliptical-pde-with-mixed-inhomogeneo/3300#3300
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.
Regards,
Francesc.
On Tue, Nov 16, 2021 at 6:55 PM Barry Smith wrote:
>
>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.
>
> 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.
>
> Barry
>
>
> On Nov 16, 2021, at 1:37 PM, Francesc Levrero-Florencio <
> f.levrero-floren...@onscale.com> wrote:
>
> Hi Barry,
>
> 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:
> - Apply updated Dirichlet BCs in order to assemble the internal forces,
> and add them to the residual.
> - Assemble the external forces, and subtract them from the residual.
> - 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.
>
> 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.
>
> Regards,
> Francesc.
>
> On Tue, Nov 16, 2021 at 5:38 PM Barry Smith wrote:
>
>>
>> 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?
>>
>>
>>
>> > On Nov 16, 2021, at 11:19 AM, Francesc Levrero-Florencio <
>> f.levrero-floren...@onscale.com> wrote:
>> >
>> > Dear PETSc team and users,
>> >
>> > 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).
>> >
>> > 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.
>> >
>> > During solving with traditional Newton, the general behaviour of each
>> time step
Re: [petsc-users] Doubt about BT and BASIC NEWTONLS
Hi Barry,
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:
- Apply updated Dirichlet BCs in order to assemble the internal forces, and
add them to the residual.
- Assemble the external forces, and subtract them from the residual.
- 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.
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.
Regards,
Francesc.
On Tue, Nov 16, 2021 at 5:38 PM Barry Smith wrote:
>
> 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?
>
>
>
> > On Nov 16, 2021, at 11:19 AM, Francesc Levrero-Florencio <
> f.levrero-floren...@onscale.com> wrote:
> >
> > Dear PETSc team and users,
> >
> > 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).
> >
> > 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.
> >
> > 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.
> >
> > Here are a few values of the relevant norms (using traditional Newton)
> in the first time step:
> >
> > BASIC NEWTON LS
> > Norm of the internal forces is 0
> > Norm of the external forces is 1374.49
> > Norm of the residual is 1374.49
> > Norm of the solution with Dirichlet BCs is 0
> > Number of SNES iteration is 0
> > -
> > Norm of the internal forces is 113498
> > Norm of the external forces is 1374.49
> > Norm of the residual is 105053
> > Norm of the solution with Dirichlet BCs is 0.441466
> > Number of SNES iteration is 0
> > -
> > Norm of the internal forces is 42953.5
> > Norm of the external forces is 1374.49
> > Norm of the residual is 11.3734
> > Norm of the solution with Dirichlet BCs is 0.441438
> > Number of SNES iteration is 1
> >
> > Here are a few values of the relevant norms (using “bt” line-search) in
> the first time step:
> >
> > BT NEWTONLS
> > Norm of the internal forces is 0
> > Norm of the external forces is 1374.49
> > Norm of the residual is 1374.49
> > Norm of the solution with Dirichlet BCs is 0
> > Number of SNES iteration is 0
> > -
> > Norm of the internal forces is 113498
> > Norm of the external forces is 1374.49
> > Norm of the residual i
[petsc-users] Doubt about BT and BASIC NEWTONLS
Dear PETSc team and users, 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). 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. 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. Here are a few values of the relevant norms (using traditional Newton) in the first time step: BASIC NEWTON LS Norm of the internal forces is 0 Norm of the external forces is 1374.49 Norm of the residual is 1374.49 Norm of the solution with Dirichlet BCs is 0 Number of SNES iteration is 0 - Norm of the internal forces is 113498 Norm of the external forces is 1374.49 Norm of the residual is 105053 Norm of the solution with Dirichlet BCs is 0.441466 Number of SNES iteration is 0 - Norm of the internal forces is 42953.5 Norm of the external forces is 1374.49 Norm of the residual is 11.3734 Norm of the solution with Dirichlet BCs is 0.441438 Number of SNES iteration is 1 Here are a few values of the relevant norms (using “bt” line-search) in the first time step: BT NEWTONLS Norm of the internal forces is 0 Norm of the external forces is 1374.49 Norm of the residual is 1374.49 Norm of the solution with Dirichlet BCs is 0 Number of SNES iteration is 0 - Norm of the internal forces is 113498 Norm of the external forces is 1374.49 Norm of the residual is 105053 Norm of the solution with Dirichlet BCs is 0.441466 Number of SNES iteration is 0 (I assume this is the first try at lambda=1) - Norm of the internal forces is 4422.12 Norm of the external forces is 1374.49 Norm of the residual is 1622.74 Norm of the solution with Dirichlet BCs is 0.0441466 Number of SNES iteration is 0 Line search: gnorm after quadratic fit 1.622742343614e+03 (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) - Norm of the internal forces is 2163.76 Norm of the external forces is 1374.49 Norm of the residual is 1331.88 Norm of the solution with Dirichlet BCs is 0.0220733 Number of SNES iteration is 0 Line search: Cubically determined step, current gnorm 1.331884625811e+03 lambda=5.0003e-02 (This is the accepted lambda for the current Newton iteration) - Norm of the internal forces is 104020 Norm of the external forces is 1374.49 Norm of the residual is 94739 Norm of the solution with Dirichlet BCs is 0.441323 Number of SNES iteration is 1 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. Thanks for your help in advance and please keep up the good work! Regards, Francesc.
[petsc-users] Fwd: Some guidance on setting an arc-length solver up in PETSc with TS/SNES
Dear PETSc team and users,
I am trying to implement a “non-consistent arc-length method” (i.e.
non-consistent as in the Jacobian from a traditional load-controlled method
is used instead of the “augmented one”, the latter would need an extra/row
column for the constraint terms; the non-consistent method requires two
linear solves, with different RHS, at every nonlinear iteration).
I am trying to do this in the context of TS/SNES. Computation of the
Jacobian remains the same, but I am thinking of doing most of the work
within the function to compute the residual. The rough non-consistent
arc-length algorithm is (delta here refers to the correction per nonlinear
iteration):
- Solve K * delta u^k_1 = - R
- Solve K * delta u^k_2 = - F_ext
- Find the correct root of delta lambda
- lambda^k = lambda^{k-1} + delta lambda^k
- delta u^k = delta u^k_1 + delta lambda^k * delta u^k_2
- u^{k} = u^{k-1} + delta u^k, u here is the solution to the mechanical
problem (displacement).
The rest is the same as a load-controlled scheme (for a more complete
arc-length algorithm, please see “A simple extrapolated predictor for
overcoming the starting and tracking issues in the arc-length method for
nonlinear structural mechanics”, Kadapa, C. (2021), *Eng Struct*). I see
two potential issues with an arc-length implementation by using TS/SNES:
- Needing to modify the solution vector within each nonlinear iteration. I
believe this could be sorted out by allowing TS/SNES to obtain the solution
to the traditional load-controlled problem (u_1 above), and only keep track
of the correction per nonlinear iteration, which would correspond to delta
u^k_1 = u^k_1 (given by the TS at the current iteration) – u^{k-1}_1 (given
by the TS at the previous iteration). I imagine the nonlinear iterative
correction would be correct because the residual (and jacobian) would be
calculated by using an internally stored u (u above, stored outside the TS,
which is not the same as the TS internally stored u_1). I believe this
would be okay in quasistatic simulations since solution derivatives are not
required. How could this be amended for transient simulations? Is there any
other simpler way to achieve this? Can the solution vector be set to a
predictor value at the beginning of a time step (maybe prestep or poststep
related functions)?
- Solving two linear systems per iteration. I imagine the way around this
is to do TSGetKSP(ts, ksp), and then KSPSolve(ksp, - F_ext, delta u^k_2)
within the TS, since the TS would be using the correctly calculated
Jacobian K above. Would this an effective way to achieve this?
Thanks for your help in advance and please keep up the good work!
Regards,
Francesc.
