Hello,
Thanks for the detailed answer!
I've just re-run the executable taking off the -snes_mf_operator and LU
options, but it doesn't look any better... :
Running ./example-opt -snes_type ls -snes_linesearch_type basic -ksp_rtol 1e-4
NL step 0, |residual|_2 = 1.936492e-05
NL step 1, |residual|_2 = 1.938856e-05
NL step 2, |residual|_2 = 1.941222e-05
NL step 3, |residual|_2 = 1.943592e-05
NL step 4, |residual|_2 = 1.945964e-05
Then I stopped it because it didn't look like converging, after one hour or so.
Monitoring the linear iterations showed:
Running ./example-opt -ksp_rtol 1e-3 -ksp_monitor
NL step 0, |residual|_2 = 1.936492e-05
0 KSP Residual norm 7.859921803229e-07
1 KSP Residual norm 3.815524538428e-07
2 KSP Residual norm 3.364144277753e-07
3 KSP Residual norm 2.909576661282e-07
4 KSP Residual norm 2.772044053816e-07
...
1654 KSP Residual norm 8.626050416153e-10
1655 KSP Residual norm 8.145093822144e-10
1656 KSP Residual norm 6.346659238618e-10
The same happens with the GAMG options (either -pc_gamg_type agg or
-pc_gamg_agg_nsmooths 1, from PETSc manual). But I could still reconfigure
PETSc and libMesh with the hypre option.
I still haven't asked my FEAP colleague about the specifications she's using,
but she told me it takes 4 or 5 nonlinear iterations to converge, depending on
the loading.
Thanks!
Lorenzo
On Nov 22, 2013, at 7:10 PM, Derek Gaston wrote:
> -snes_mf_operator is a PETSc option specifying that you want to do
> preconditioned Jacobian Free Newton Krylov (JFNK). Is the guy using FEAP
> doing JFNK?
>
> JFNK means that every linear iteration inside a Newton step you must
> recompute your residual. That means a full sweep over the mesh,
> re-evaluating your material properties and residual statements at every
> quadrature point to assemble a full residual vector. This is expensive.
>
> If you are just doing a single physics problem (and it sounds like you are)
> and you have the ability to compute the exact Jacobian (which it sounds like
> you do) then using -snes_mf_operator will more than likely be MUCH slower
> than just solving using Newton's method like normal (where you assemble a
> matrix and a RHS just _once_ per Newton step - then throw a linear solver at
> it).
>
> To use "regular" Newton just leave that option off.
>
> Further, you have specified "-pc_type lu"... which is generally a bad idea
> for performance. With LU you are doing a direct inversion of your Jacobian
> matrix (old-school style!) and using it as a preconditioner. It is generally
> much better to use an "inexact" Newton formulation where you don't perfectly
> invert your Jacobian matrix and instead let your Krylov solver solve to some
> (fairly loose) tolerance (ie use -ksp_rtol 1e-4 or larger) so that you are
> not "oversolving" your linear problem inside each Newton step.
>
> Instead of using LU I highly recommend using an algebraic multigrid
> preconditioner for solid mechanics. Look into using GAMG in PETSc or use
> --download-hypre when configuring PETSc and use: -pc_type hypre
> -pc_hypre_type boomeramg.
>
>
> Basically: Your solver options are non-optimal. Make sure you are solving
> using _exactly_ the same solver options between libMesh and FEAP before doing
> any comparisons. You should probably run with -ksp_monitor to show you how
> many linear iterations you're taking and then make sure that both your
> libMesh and FEAP implementations are taking the same (or VERY similar) number
> of both nonlinear _and_ linear iterations
>
>
> Derek
>
>
>
> On Fri, Nov 22, 2013 at 10:08 AM, Lorenzo Zanon <[email protected]>
> wrote:
> Hello,
>
> I have an example based on miscellaneous_ex3.C, where I implemented a
> nonlinear elastic problem with StVenant stress-strain law on a 3D domain
> 5x1x1, blocked at x=0 and with an applied load at x=5. The problem is, the
> computation of the displacement (along x y and z) is really slow (hours), on
> our cluster it goes out of CPU time already with a mesh of 64x8x1 (the
> loading acts only on the y-direction, so no more elements are needed along
> z). 4 or 5 Newton steps should be enough for solving the problem...
>
> A colleague of mine implemented the same problem on the software called FEAP,
> and it takes only a few minutes there. I think my implementation is correct
> because the results on a very coarse mesh (10x2x1) are roughly similar to the
> FEAP ones on a finer mesh. I don't have any problems for the 2D case also.
>
> Is there anything I can do? I'm running in opt mode with the following
> options concerning the Newton solver:
>
> -snes_type ls -snes_linesearch_type basic -snes_mf_operator -pc_type lu
> -pc_factor_nonzeros_along_diagonal
>
> Thanks!
> Lorenzo
> ------------------------------------------------------------------------------
> Shape the Mobile Experience: Free Subscription
> Software experts and developers: Be at the forefront of tech innovation.
> Intel(R) Software Adrenaline delivers strategic insight and game-changing
> conversations that shape the rapidly evolving mobile landscape. Sign up now.
> http://pubads.g.doubleclick.net/gampad/clk?id=63431311&iu=/4140/ostg.clktrk
> _______________________________________________
> Libmesh-users mailing list
> [email protected]
> https://lists.sourceforge.net/lists/listinfo/libmesh-users
>
------------------------------------------------------------------------------
Shape the Mobile Experience: Free Subscription
Software experts and developers: Be at the forefront of tech innovation.
Intel(R) Software Adrenaline delivers strategic insight and game-changing
conversations that shape the rapidly evolving mobile landscape. Sign up now.
http://pubads.g.doubleclick.net/gampad/clk?id=63431311&iu=/4140/ostg.clktrk
_______________________________________________
Libmesh-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users
