Re: [firedrake] Understanding solver parameters in NonlinearVariationalSolver
I think I mistakenly clicked reply instead of reply all. Attached is my response: -- Thanks Matt, Thanks for the swift response. These are the logs with basic, backtracking, cp and nleqerr. If I don't specify 'snes_linesearch_type' it defaults to 'basic'. At least that's what I could gather from the log. - `cp` seemingly diverges the fastest. Of course I haven't touched the other options. - `nleqerr` goes the furthest I believe, but not that far from `basic` or `bt`. My concern was that since the forms are essentially the same, and I am passing "seemingly" the same parameters in FEniCS and Firedrake (both of which are using PETSc's SNES), what is it that is different which causes the solver to fail in one but not the other. Edit: Just saw Ed's message. Indeed it defaults to `basic` For definiteness this was the dictionary of solver parameters I am passing to NLVS in Firedrake solver_params={ "snes_type": "newtonls", "ksp_type": "preonly", "pc_type": "lu", "snes_max_it": 100, # "snes_stol": 1.e-6, # "snes_rtol": 1.e-6, # "snes_atol": 1.e-6, "snes_view": None, # "ksp_atol": 1.e-10, # "ksp_rtol":1.e-8, "ksp_view": None, "snes_linesearch_type": "bt", "snes_converged_reason": None, "snes_monitor": None, "ksp_converged_reason": None, "ksp_monitor_true_residual": None } On Sun, 27 Jun 2021 at 11:22, Bhavesh Shrimali <bhavesh.shrimali@gmail.com> wrote:
Thanks Matt,
Thanks for the swift response. These are the logs with basic, backtracking, cp and nleqerr. If I don't specify 'snes_linesearch_type' it defaults to 'basic'. At least that's what I could gather from the log.
- `cp` seemingly diverges the fastest. Of course I haven't touched the other options. - `nleqerr` goes the furthest I believe, but not that far from `basic` or `bt`.
My concern was that since the forms are essentially the same, and I am passing "seemingly" the same parameters in FEniCS and Firedrake (both of which are using PETSc's SNES), what is it that is different which causes the solver to fail in one but not the other.
Edit: Just saw Ed's message. Indeed it defaults to `basic`
For definiteness this was the dictionary of solver parameters I am passing to NLVS in Firedrake
solver_params={ "snes_type": "newtonls", "ksp_type": "preonly", "pc_type": "lu", "snes_max_it": 100, # "snes_stol": 1.e-6, # "snes_rtol": 1.e-6, # "snes_atol": 1.e-6, "snes_view": None, # "ksp_atol": 1.e-10, # "ksp_rtol":1.e-8, "ksp_view": None, "snes_linesearch_type": "bt", "snes_converged_reason": None, "snes_monitor": None, "ksp_converged_reason": None, "ksp_monitor_true_residual": None }
On Sun, 27 Jun 2021 at 06:00, Matthew Knepley <knepley@gmail.com> wrote:
Give
-snes_view -snes_converged_reason -snes_monitor -ksp_converged_reason -ksp_monitor_true_residual
and send all the output. That way I might be able to see what is happening when you do not converge. Also, line search 'basic' is not a good choice usually. I would either use the default, or 'cp', or maybe 'nleqerr' that Patrick added.
Thanks,
Matt
On Sat, Jun 26, 2021 at 10:14 PM Bhavesh Shrimali < bhavesh.shrimali@gmail.com> wrote:
Hi everyone,
I am trying to understand the usage of solver parameters that are passed to the NonlinearVariationalSolver in firedrake. In order to do some tests, I decided to reproduce some of my previous calculations done in FEniCS for nonlinear viscoelasticity. In old dolfin, I had the following for solver parameters to be passed to PETSc's SNES:
snes_solver_parameters = {"nonlinear_solver": "snes", "snes_solver": {"linear_solver": "lu", "maximum_iterations": 50, "report": True, "line_search": 'basic', "error_on_nonconvergence": False, "relative_tolerance": 1e-7, "absolute_tolerance": 1e-8}} ... ... solver = NonlinearVariationalSolver(problem) solver.parameters.update(snes_solver_parameters)
Now as I understand, in Firedrake the options passed to the solver are consistent with the corresponding names in PETSc (correct me if I am wrong), and therefore I did:
solver_params={ "snes_type": "newtonls", "ksp_type": "preonly", "pc_type": "lu", "snes_max_it": 100, "snes_rtol": 1.e-7, "snes_atol": 1.e-8, "snes_view": None, "ksp_atol": 1.e-10, # not needed? "ksp_rtol":1.e-8, "ksp_view": None, "snes_linesearch_type": "basic", "snes_no_convergence_test": 1 # this has no effect } solver = NonlinearVariationalSolver(problem, solver_parameters=solver_params)
However, for some reason, I am unable to get the solver to converge when trying to solve the exact same problem as FEniCS. I experimented with different quadrature degrees as well as tweaked the solver parameters for SNES to best match what I am doing in dolfin. I understand that the form compiler is different hence the precise code that is generated may be different (I am not an expert on code generation) even though the python scripts may look similar due to UFL syntax.
In fact the weird part is that the solver seems to converge fine till a certain point in my input (loading). You can see this in the plot below where I plot the final stress-stretch relation obtained at a point in my domain from dolfin and firedrake..
[image: firedrake_vs_fenics.png]
I have marked the abscissa beyond which the solver fails to converge in Firedrake. The complete codes are here: https://github.com/bhaveshshrimali/comparingFiredrakeFEniCS
I tried setting the option `snes_no_convergence_test` to 1/True but that had no effect as the SNES solver still diverges with `DIVERGED_DTOL`. In the end I put the call to solve inside a try/catch block to at least be able to see the solution till that point. And it surprisingly is accurate till that point.
Would anyone have any clue of what could be going wrong? I wouldn't expect that you go through the entire code(s) as that could be fairly long, but any pointers would be helpful :)
I am using firedrake via Docker on windows.
Thanks
--
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener
https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
--
*Bhavesh Shrimali* *Ph.D. Student, **Civil and Environmental Engineering,*
*University of Illinois Urbana-Champaign*
*Champaign, IL, United States*
*Email : bshrima2@illinois.edu <bshrima2@illinois.edu>* *Web : **https://bhaveshshrimali.github.io/ <https://bhaveshshrimali.github.io/>*
-- *Bhavesh Shrimali* *Ph.D. Student, **Civil and Environmental Engineering,* *University of Illinois Urbana-Champaign* *Champaign, IL, United States* *Email : bshrima2@illinois.edu <bshrima2@illinois.edu>* *Web : **https://bhaveshshrimali.github.io/ <https://bhaveshshrimali.github.io/>*
Something is wrong with the Newton system in the last step. You can see that the LU residual is steadily increasing. At the last step it is 10^-6. The iteration appears to be approaching a singular Jacobian. This seems like a formulation problem. Thanks, Matt On Sun, Jun 27, 2021 at 11:24 AM Bhavesh Shrimali < bhavesh.shrimali@gmail.com> wrote:
I think I mistakenly clicked reply instead of reply all. Attached is my response:
--
Thanks Matt,
Thanks for the swift response. These are the logs with basic, backtracking, cp and nleqerr. If I don't specify 'snes_linesearch_type' it defaults to 'basic'. At least that's what I could gather from the log.
- `cp` seemingly diverges the fastest. Of course I haven't touched the other options. - `nleqerr` goes the furthest I believe, but not that far from `basic` or `bt`.
My concern was that since the forms are essentially the same, and I am passing "seemingly" the same parameters in FEniCS and Firedrake (both of which are using PETSc's SNES), what is it that is different which causes the solver to fail in one but not the other.
Edit: Just saw Ed's message. Indeed it defaults to `basic`
For definiteness this was the dictionary of solver parameters I am passing to NLVS in Firedrake
solver_params={ "snes_type": "newtonls", "ksp_type": "preonly", "pc_type": "lu", "snes_max_it": 100, # "snes_stol": 1.e-6, # "snes_rtol": 1.e-6, # "snes_atol": 1.e-6, "snes_view": None, # "ksp_atol": 1.e-10, # "ksp_rtol":1.e-8, "ksp_view": None, "snes_linesearch_type": "bt", "snes_converged_reason": None, "snes_monitor": None, "ksp_converged_reason": None, "ksp_monitor_true_residual": None }
On Sun, 27 Jun 2021 at 11:22, Bhavesh Shrimali <bhavesh.shrimali@gmail.com> wrote:
Thanks Matt,
Thanks for the swift response. These are the logs with basic, backtracking, cp and nleqerr. If I don't specify 'snes_linesearch_type' it defaults to 'basic'. At least that's what I could gather from the log.
- `cp` seemingly diverges the fastest. Of course I haven't touched the other options. - `nleqerr` goes the furthest I believe, but not that far from `basic` or `bt`.
My concern was that since the forms are essentially the same, and I am passing "seemingly" the same parameters in FEniCS and Firedrake (both of which are using PETSc's SNES), what is it that is different which causes the solver to fail in one but not the other.
Edit: Just saw Ed's message. Indeed it defaults to `basic`
For definiteness this was the dictionary of solver parameters I am passing to NLVS in Firedrake
solver_params={ "snes_type": "newtonls", "ksp_type": "preonly", "pc_type": "lu", "snes_max_it": 100, # "snes_stol": 1.e-6, # "snes_rtol": 1.e-6, # "snes_atol": 1.e-6, "snes_view": None, # "ksp_atol": 1.e-10, # "ksp_rtol":1.e-8, "ksp_view": None, "snes_linesearch_type": "bt", "snes_converged_reason": None, "snes_monitor": None, "ksp_converged_reason": None, "ksp_monitor_true_residual": None }
On Sun, 27 Jun 2021 at 06:00, Matthew Knepley <knepley@gmail.com> wrote:
Give
-snes_view -snes_converged_reason -snes_monitor -ksp_converged_reason -ksp_monitor_true_residual
and send all the output. That way I might be able to see what is happening when you do not converge. Also, line search 'basic' is not a good choice usually. I would either use the default, or 'cp', or maybe 'nleqerr' that Patrick added.
Thanks,
Matt
On Sat, Jun 26, 2021 at 10:14 PM Bhavesh Shrimali < bhavesh.shrimali@gmail.com> wrote:
Hi everyone,
I am trying to understand the usage of solver parameters that are passed to the NonlinearVariationalSolver in firedrake. In order to do some tests, I decided to reproduce some of my previous calculations done in FEniCS for nonlinear viscoelasticity. In old dolfin, I had the following for solver parameters to be passed to PETSc's SNES:
snes_solver_parameters = {"nonlinear_solver": "snes", "snes_solver": {"linear_solver": "lu", "maximum_iterations": 50, "report": True, "line_search": 'basic', "error_on_nonconvergence": False, "relative_tolerance": 1e-7, "absolute_tolerance": 1e-8}} ... ... solver = NonlinearVariationalSolver(problem) solver.parameters.update(snes_solver_parameters)
Now as I understand, in Firedrake the options passed to the solver are consistent with the corresponding names in PETSc (correct me if I am wrong), and therefore I did:
solver_params={ "snes_type": "newtonls", "ksp_type": "preonly", "pc_type": "lu", "snes_max_it": 100, "snes_rtol": 1.e-7, "snes_atol": 1.e-8, "snes_view": None, "ksp_atol": 1.e-10, # not needed? "ksp_rtol":1.e-8, "ksp_view": None, "snes_linesearch_type": "basic", "snes_no_convergence_test": 1 # this has no effect } solver = NonlinearVariationalSolver(problem, solver_parameters=solver_params)
However, for some reason, I am unable to get the solver to converge when trying to solve the exact same problem as FEniCS. I experimented with different quadrature degrees as well as tweaked the solver parameters for SNES to best match what I am doing in dolfin. I understand that the form compiler is different hence the precise code that is generated may be different (I am not an expert on code generation) even though the python scripts may look similar due to UFL syntax.
In fact the weird part is that the solver seems to converge fine till a certain point in my input (loading). You can see this in the plot below where I plot the final stress-stretch relation obtained at a point in my domain from dolfin and firedrake..
[image: firedrake_vs_fenics.png]
I have marked the abscissa beyond which the solver fails to converge in Firedrake. The complete codes are here: https://github.com/bhaveshshrimali/comparingFiredrakeFEniCS
I tried setting the option `snes_no_convergence_test` to 1/True but that had no effect as the SNES solver still diverges with ` DIVERGED_DTOL`. In the end I put the call to solve inside a try/catch block to at least be able to see the solution till that point. And it surprisingly is accurate till that point.
Would anyone have any clue of what could be going wrong? I wouldn't expect that you go through the entire code(s) as that could be fairly long, but any pointers would be helpful :)
I am using firedrake via Docker on windows.
Thanks
--
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener
https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
--
*Bhavesh Shrimali* *Ph.D. Student, **Civil and Environmental Engineering,*
*University of Illinois Urbana-Champaign*
*Champaign, IL, United States*
*Email : bshrima2@illinois.edu <bshrima2@illinois.edu>* *Web : **https://bhaveshshrimali.github.io/ <https://bhaveshshrimali.github.io/>*
--
*Bhavesh Shrimali* *Ph.D. Student, **Civil and Environmental Engineering,*
*University of Illinois Urbana-Champaign*
*Champaign, IL, United States*
*Email : bshrima2@illinois.edu <bshrima2@illinois.edu>* *Web : **https://bhaveshshrimali.github.io/ <https://bhaveshshrimali.github.io/>* _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
Thanks Matt, The only difference I see in the forms is that I am using `project` in Firedrake which I believe does a cell-wise projection for a DG-space whereas I use the `LocalSolver` in dolfin to do a cell-wise projection myself. Even the SNES function norms that I am seeing aren't the same between FEniCS and Firedrake. Possible that I made some mistake in the formulation. Thanks for pointing that out. I'll investigate this further today Thanks, Bhavesh On Sun, 27 Jun 2021 at 13:10, Matthew Knepley <knepley@gmail.com> wrote:
Something is wrong with the Newton system in the last step. You can see that the LU residual is steadily increasing. At the last step it is 10^-6. The iteration appears to be approaching a singular Jacobian. This seems like a formulation problem.
Thanks,
Matt
On Sun, Jun 27, 2021 at 11:24 AM Bhavesh Shrimali < bhavesh.shrimali@gmail.com> wrote:
I think I mistakenly clicked reply instead of reply all. Attached is my response:
--
Thanks Matt,
Thanks for the swift response. These are the logs with basic, backtracking, cp and nleqerr. If I don't specify 'snes_linesearch_type' it defaults to 'basic'. At least that's what I could gather from the log.
- `cp` seemingly diverges the fastest. Of course I haven't touched the other options. - `nleqerr` goes the furthest I believe, but not that far from `basic` or `bt`.
My concern was that since the forms are essentially the same, and I am passing "seemingly" the same parameters in FEniCS and Firedrake (both of which are using PETSc's SNES), what is it that is different which causes the solver to fail in one but not the other.
Edit: Just saw Ed's message. Indeed it defaults to `basic`
For definiteness this was the dictionary of solver parameters I am passing to NLVS in Firedrake
solver_params={ "snes_type": "newtonls", "ksp_type": "preonly", "pc_type": "lu", "snes_max_it": 100, # "snes_stol": 1.e-6, # "snes_rtol": 1.e-6, # "snes_atol": 1.e-6, "snes_view": None, # "ksp_atol": 1.e-10, # "ksp_rtol":1.e-8, "ksp_view": None, "snes_linesearch_type": "bt", "snes_converged_reason": None, "snes_monitor": None, "ksp_converged_reason": None, "ksp_monitor_true_residual": None }
On Sun, 27 Jun 2021 at 11:22, Bhavesh Shrimali < bhavesh.shrimali@gmail.com> wrote:
Thanks Matt,
Thanks for the swift response. These are the logs with basic, backtracking, cp and nleqerr. If I don't specify 'snes_linesearch_type' it defaults to 'basic'. At least that's what I could gather from the log.
- `cp` seemingly diverges the fastest. Of course I haven't touched the other options. - `nleqerr` goes the furthest I believe, but not that far from `basic` or `bt`.
My concern was that since the forms are essentially the same, and I am passing "seemingly" the same parameters in FEniCS and Firedrake (both of which are using PETSc's SNES), what is it that is different which causes the solver to fail in one but not the other.
Edit: Just saw Ed's message. Indeed it defaults to `basic`
For definiteness this was the dictionary of solver parameters I am passing to NLVS in Firedrake
solver_params={ "snes_type": "newtonls", "ksp_type": "preonly", "pc_type": "lu", "snes_max_it": 100, # "snes_stol": 1.e-6, # "snes_rtol": 1.e-6, # "snes_atol": 1.e-6, "snes_view": None, # "ksp_atol": 1.e-10, # "ksp_rtol":1.e-8, "ksp_view": None, "snes_linesearch_type": "bt", "snes_converged_reason": None, "snes_monitor": None, "ksp_converged_reason": None, "ksp_monitor_true_residual": None }
On Sun, 27 Jun 2021 at 06:00, Matthew Knepley <knepley@gmail.com> wrote:
Give
-snes_view -snes_converged_reason -snes_monitor -ksp_converged_reason -ksp_monitor_true_residual
and send all the output. That way I might be able to see what is happening when you do not converge. Also, line search 'basic' is not a good choice usually. I would either use the default, or 'cp', or maybe 'nleqerr' that Patrick added.
Thanks,
Matt
On Sat, Jun 26, 2021 at 10:14 PM Bhavesh Shrimali < bhavesh.shrimali@gmail.com> wrote:
Hi everyone,
I am trying to understand the usage of solver parameters that are passed to the NonlinearVariationalSolver in firedrake. In order to do some tests, I decided to reproduce some of my previous calculations done in FEniCS for nonlinear viscoelasticity. In old dolfin, I had the following for solver parameters to be passed to PETSc's SNES:
snes_solver_parameters = {"nonlinear_solver": "snes", "snes_solver": {"linear_solver": "lu", "maximum_iterations": 50, "report": True, "line_search": 'basic', "error_on_nonconvergence": False, "relative_tolerance": 1e-7, "absolute_tolerance": 1e-8}} ... ... solver = NonlinearVariationalSolver(problem) solver.parameters.update(snes_solver_parameters)
Now as I understand, in Firedrake the options passed to the solver are consistent with the corresponding names in PETSc (correct me if I am wrong), and therefore I did:
solver_params={ "snes_type": "newtonls", "ksp_type": "preonly", "pc_type": "lu", "snes_max_it": 100, "snes_rtol": 1.e-7, "snes_atol": 1.e-8, "snes_view": None, "ksp_atol": 1.e-10, # not needed? "ksp_rtol":1.e-8, "ksp_view": None, "snes_linesearch_type": "basic", "snes_no_convergence_test": 1 # this has no effect } solver = NonlinearVariationalSolver(problem, solver_parameters=solver_params)
However, for some reason, I am unable to get the solver to converge when trying to solve the exact same problem as FEniCS. I experimented with different quadrature degrees as well as tweaked the solver parameters for SNES to best match what I am doing in dolfin. I understand that the form compiler is different hence the precise code that is generated may be different (I am not an expert on code generation) even though the python scripts may look similar due to UFL syntax.
In fact the weird part is that the solver seems to converge fine till a certain point in my input (loading). You can see this in the plot below where I plot the final stress-stretch relation obtained at a point in my domain from dolfin and firedrake..
[image: firedrake_vs_fenics.png]
I have marked the abscissa beyond which the solver fails to converge in Firedrake. The complete codes are here: https://github.com/bhaveshshrimali/comparingFiredrakeFEniCS
I tried setting the option `snes_no_convergence_test` to 1/True but that had no effect as the SNES solver still diverges with ` DIVERGED_DTOL`. In the end I put the call to solve inside a try/catch block to at least be able to see the solution till that point. And it surprisingly is accurate till that point.
Would anyone have any clue of what could be going wrong? I wouldn't expect that you go through the entire code(s) as that could be fairly long, but any pointers would be helpful :)
I am using firedrake via Docker on windows.
Thanks
--
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener
https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
--
*Bhavesh Shrimali* *Ph.D. Student, **Civil and Environmental Engineering,*
*University of Illinois Urbana-Champaign*
*Champaign, IL, United States*
*Email : bshrima2@illinois.edu <bshrima2@illinois.edu>* *Web : **https://bhaveshshrimali.github.io/ <https://bhaveshshrimali.github.io/>*
--
*Bhavesh Shrimali* *Ph.D. Student, **Civil and Environmental Engineering,*
*University of Illinois Urbana-Champaign*
*Champaign, IL, United States*
*Email : bshrima2@illinois.edu <bshrima2@illinois.edu>* *Web : **https://bhaveshshrimali.github.io/ <https://bhaveshshrimali.github.io/>* _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener
https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
-- *Bhavesh Shrimali* *Ph.D. Student, **Civil and Environmental Engineering,* *University of Illinois Urbana-Champaign* *Champaign, IL, United States* *Email : bshrima2@illinois.edu <bshrima2@illinois.edu>* *Web : **https://bhaveshshrimali.github.io/ <https://bhaveshshrimali.github.io/>*
After some more digging I am still unable to get the same solver setup using the same set of parameters. This time trying with tsfc for the form compiler (an old fork though: git+ https://github.com/blechta/tsfc.git@2018.1.0). Even with the same quadrature degree and the same form I get different SNES function norms in the very first solving step. Note that there are no convergence issues in the first step but eventually there are later. Though with different PETSc versions should I expect similar convergence behavior ? There is very little difference in my formulation since most of it is in UFL syntax. I will keep looking ... In the meantime, here are the outputs -- Firedrake: 0 SNES Function norm 4.630756653490e+03 1 SNES Function norm 6.610067614593e+02 2 SNES Function norm 8.595689678029e+00 3 SNES Function norm 7.177260814570e-03 4 SNES Function norm 7.839549826428e-09 Nonlinear firedrake_0_ solve converged due to CONVERGED_FNORM_RELATIVE iterations 4 Dolfin: Solving nonlinear variational problem. 0 SNES Function norm 4.048156097083e+03 1 SNES Function norm 7.701194025455e+01 2 SNES Function norm 9.032060723182e-02 3 SNES Function norm 2.188411552504e-08 PETSc SNES solver converged in 3 iterations with convergence reason CONVERGED_FNORM_RELATIVE. and the codes (I will try to reduce these further as they still are far from minimal (~ 30 lines)) Dolfin (petsc4py version: 3.12.0): import numpy as np from dolfin import * parameters["form_compiler"]["cpp_optimize"] = True parameters["form_compiler"]["quadrature_degree"] = 5 parameters["form_compiler"]["representation"] = "tsfc" parameters["linear_algebra_backend"] = "PETSc" snes_solver_parameters = {"nonlinear_solver": "snes", "snes_solver": {"linear_solver": "lu", "maximum_iterations": 50, "report": True, "line_search": 'basic', "error_on_nonconvergence": True} } def Psi(u, p, Cv): I = Identity(3) F = I + grad(u) C = F.T * F Cv_inv = inv(Cv) I1 = tr(C) JF = det(F) Ie1 = inner(C, Cv_inv) Je = JF / sqrt(det(Cv)) Jstar = (lam1 + p + sqrt((p + lam1)**2 + 4 * lam1 * (mu1 + mu2))) / (2 * lam1) psi_Eq = ((3**(1 - alph1)) / (2 * alph1)) * mu1 * (I1**alph1 - 3**alph1) +\ ((3**(1 - alph2)) / (2 * alph2)) * mu2 * (I1**alph2 - 3**alph2) -\ (mu1 + mu2) * ln(Jstar) + 0.5 * lam1 * (Jstar - 1)**2 psi_Neq = ((3**(1 - a1)) / (2 * a1)) * m1 * (Ie1**a1 - 3**a1) +\ ((3**(1 - a2)) / (2 * a2)) * m2 * (Ie1**a2 - 3**a2) -\ (m1 + m2) * ln(Je) W_hat = (psi_Eq + p * (JF - Jstar) + psi_Neq) * dx return W_hat mesh = UnitCubeMesh(3, 3, 3) W1 = VectorElement('CG', mesh.ufl_cell(), 2) W2 = FiniteElement('CG', mesh.ufl_cell(), 1) W3 = TensorElement("DG", mesh.ufl_cell(), 0) V = FunctionSpace(mesh, W1 * W2) V_Cv = FunctionSpace(mesh, W3) du, dp = TrialFunctions(V) dw = TrialFunction(V) v, ph = TestFunctions(V) wh = TestFunction(V) w = Function(V) u, p = split(w) Cv = Function(V_Cv) Cv.interpolate(Constant(np.eye(3))) left_fc = CompiledSubDomain( "near(x[0], 0.) && on_boundary" ) right_fc = CompiledSubDomain( "near(x[0], 1.) && on_boundary" ) bottom_fc = CompiledSubDomain( "near(x[1], 0.) && on_boundary" ) back_fc = CompiledSubDomain( "near(x[2], 0.) && on_boundary" ) strtch = Constant(0.0) bc_right = DirichletBC(V.sub(0).sub(0), strtch, right_fc) bc_left = DirichletBC(V.sub(0).sub(0), Constant(0.0), left_fc) bc_bottom = DirichletBC(V.sub(0).sub(1), Constant(0.0), bottom_fc) bc_back = DirichletBC(V.sub(0).sub(2), Constant(0.0), back_fc) bcs = [bc_right, bc_left, bc_bottom, bc_back] # Material Properties VHB4910 mu1 = Constant(13.54 * 1e3) mu2 = Constant(1.08 * 1e3) lam1 = Constant((mu1 + mu2) * 10**3) alph1 = Constant(1.) alph2 = Constant(-2.474) m1 = Constant(5.42 * 1e3) m2 = Constant(20.78 * 1e3) a1 = Constant(-10.) a2 = Constant(1.948) K1 = Constant(3507 * 1e3) K2 = Constant(1.0 * 1.e-12) bta1 = Constant(1.852) bta2 = Constant(0.26) eta0 = Constant(7014 * 1e3) etaInf = Constant(0.1 * 1e3) S = derivative(Psi(u, p, Cv), w, wh) Jac = derivative(S, w, dw) problem = NonlinearVariationalProblem(S, w, bcs, J=Jac) solver = NonlinearVariationalSolver(problem) solver.parameters.update(snes_solver_parameters) lamdot = 0.01 # Loading Rate Tfinal = 2 * (2 / lamdot) n_count = 101 dt = Tfinal / (n_count - 1) timeVals = np.linspace(0, Tfinal, n_count) stretchVals = np.hstack((lamdot * timeVals[:len(timeVals) // 2], lamdot * (-timeVals[len(timeVals) // 2:] + 2 * timeVals[len(timeVals) // 2]))) strtch.assign(stretchVals[1]) solver.solve() u, p = w.split() Firedrake (petsc4py version: 3.15.1) import numpy as np from firedrake import * parameters["form_compiler"]["quadrature_degree"] = 5 snes_solver_parameters ={"snes_type": "newtonls", "ksp_type": "preonly", "pc_type": "lu", "snes_max_it": 50, "snes_linesearch_type": "basic", "snes_monitor": None, "snes_converged_reason": None} def Psi(u, p, Cv): I = Identity(3) F = I + grad(u) C = F.T * F Cv_inv = inv(Cv) I1 = tr(C) JF = det(F) Ie1 = inner(C, Cv_inv) Je = JF / sqrt(det(Cv)) Jstar = (lam1 + p + sqrt((p + lam1)**2 + 4 * lam1 * (mu1 + mu2))) / (2 * lam1) psi_Eq = ((3**(1 - alph1)) / (2 * alph1)) * mu1 * (I1**alph1 - 3**alph1) +\ ((3**(1 - alph2)) / (2 * alph2)) * mu2 * (I1**alph2 - 3**alph2) -\ (mu1 + mu2) * ln(Jstar) + 0.5 * lam1 * (Jstar - 1)**2 psi_Neq = ((3**(1 - a1)) / (2 * a1)) * m1 * (Ie1**a1 - 3**a1) +\ ((3**(1 - a2)) / (2 * a2)) * m2 * (Ie1**a2 - 3**a2) -\ (m1 + m2) * ln(Je) W_hat = (psi_Eq + p * (JF - Jstar) + psi_Neq) * dx return W_hat mesh = UnitCubeMesh(3, 3, 3) W1 = VectorElement('CG', mesh.ufl_cell(), 2) W2 = FiniteElement('CG', mesh.ufl_cell(), 1) W3 = TensorElement("DG", mesh.ufl_cell(), 0) V = FunctionSpace(mesh, W1 * W2) V_Cv = FunctionSpace(mesh, W3) du, dp = TrialFunctions(V) dw = TrialFunction(V) v, ph = TestFunctions(V) wh = TestFunction(V) w = Function(V) u, p = split(w) Cv = Function(V_Cv) Cv.interpolate(Constant(np.eye(3))) left_fc = 1 right_fc = 2 bottom_fc = 3 back_fc = 5 strtch = Constant(0.0) bc_right = DirichletBC(V.sub(0).sub(0), strtch, right_fc) bc_left = DirichletBC(V.sub(0).sub(0), Constant(0.0), left_fc) bc_bottom = DirichletBC(V.sub(0).sub(1), Constant(0.0), bottom_fc) bc_back = DirichletBC(V.sub(0).sub(2), Constant(0.0), back_fc) bcs = [bc_right, bc_left, bc_bottom, bc_back] # Material Properties VHB4910 mu1 = Constant(13.54 * 1e3) mu2 = Constant(1.08 * 1e3) lam1 = Constant((mu1 + mu2) * 10**3) alph1 = Constant(1.) alph2 = Constant(-2.474) m1 = Constant(5.42 * 1e3) m2 = Constant(20.78 * 1e3) a1 = Constant(-10.) a2 = Constant(1.948) K1 = Constant(3507 * 1e3) K2 = Constant(1.0 * 1.e-12) bta1 = Constant(1.852) bta2 = Constant(0.26) eta0 = Constant(7014 * 1e3) etaInf = Constant(0.1 * 1e3) S = derivative(Psi(u, p, Cv), w, wh) Jac = derivative(S, w, dw) problem = NonlinearVariationalProblem(S, w, bcs, J=Jac) solver = NonlinearVariationalSolver(problem, solver_parameters=snes_solver_parameters) lamdot = 0.01 # Loading Rate Tfinal = 2 * (2 / lamdot) n_count = 101 dt = Tfinal / (n_count - 1) timeVals = np.linspace(0, Tfinal, n_count) stretchVals = np.hstack((lamdot * timeVals[:len(timeVals) // 2], lamdot * (-timeVals[len(timeVals) // 2:] + 2 * timeVals[len(timeVals) // 2]))) strtch.assign(stretchVals[1]) solver.solve() u, p = w.split() Thanks for your help! On Sun, 27 Jun 2021 at 13:21, Bhavesh Shrimali <bhavesh.shrimali@gmail.com> wrote:
Thanks Matt,
The only difference I see in the forms is that I am using `project` in Firedrake which I believe does a cell-wise projection for a DG-space whereas I use the `LocalSolver` in dolfin to do a cell-wise projection myself.
Even the SNES function norms that I am seeing aren't the same between FEniCS and Firedrake. Possible that I made some mistake in the formulation.
Thanks for pointing that out. I'll investigate this further today
Thanks, Bhavesh
On Sun, 27 Jun 2021 at 13:10, Matthew Knepley <knepley@gmail.com> wrote:
Something is wrong with the Newton system in the last step. You can see that the LU residual is steadily increasing. At the last step it is 10^-6. The iteration appears to be approaching a singular Jacobian. This seems like a formulation problem.
Thanks,
Matt
On Sun, Jun 27, 2021 at 11:24 AM Bhavesh Shrimali < bhavesh.shrimali@gmail.com> wrote:
I think I mistakenly clicked reply instead of reply all. Attached is my response:
--
Thanks Matt,
Thanks for the swift response. These are the logs with basic, backtracking, cp and nleqerr. If I don't specify 'snes_linesearch_type' it defaults to 'basic'. At least that's what I could gather from the log.
- `cp` seemingly diverges the fastest. Of course I haven't touched the other options. - `nleqerr` goes the furthest I believe, but not that far from `basic` or `bt`.
My concern was that since the forms are essentially the same, and I am passing "seemingly" the same parameters in FEniCS and Firedrake (both of which are using PETSc's SNES), what is it that is different which causes the solver to fail in one but not the other.
Edit: Just saw Ed's message. Indeed it defaults to `basic`
For definiteness this was the dictionary of solver parameters I am passing to NLVS in Firedrake
solver_params={ "snes_type": "newtonls", "ksp_type": "preonly", "pc_type": "lu", "snes_max_it": 100, # "snes_stol": 1.e-6, # "snes_rtol": 1.e-6, # "snes_atol": 1.e-6, "snes_view": None, # "ksp_atol": 1.e-10, # "ksp_rtol":1.e-8, "ksp_view": None, "snes_linesearch_type": "bt", "snes_converged_reason": None, "snes_monitor": None, "ksp_converged_reason": None, "ksp_monitor_true_residual": None }
On Sun, 27 Jun 2021 at 11:22, Bhavesh Shrimali < bhavesh.shrimali@gmail.com> wrote:
Thanks Matt,
Thanks for the swift response. These are the logs with basic, backtracking, cp and nleqerr. If I don't specify 'snes_linesearch_type' it defaults to 'basic'. At least that's what I could gather from the log.
- `cp` seemingly diverges the fastest. Of course I haven't touched the other options. - `nleqerr` goes the furthest I believe, but not that far from `basic` or `bt`.
My concern was that since the forms are essentially the same, and I am passing "seemingly" the same parameters in FEniCS and Firedrake (both of which are using PETSc's SNES), what is it that is different which causes the solver to fail in one but not the other.
Edit: Just saw Ed's message. Indeed it defaults to `basic`
For definiteness this was the dictionary of solver parameters I am passing to NLVS in Firedrake
solver_params={ "snes_type": "newtonls", "ksp_type": "preonly", "pc_type": "lu", "snes_max_it": 100, # "snes_stol": 1.e-6, # "snes_rtol": 1.e-6, # "snes_atol": 1.e-6, "snes_view": None, # "ksp_atol": 1.e-10, # "ksp_rtol":1.e-8, "ksp_view": None, "snes_linesearch_type": "bt", "snes_converged_reason": None, "snes_monitor": None, "ksp_converged_reason": None, "ksp_monitor_true_residual": None }
On Sun, 27 Jun 2021 at 06:00, Matthew Knepley <knepley@gmail.com> wrote:
Give
-snes_view -snes_converged_reason -snes_monitor -ksp_converged_reason -ksp_monitor_true_residual
and send all the output. That way I might be able to see what is happening when you do not converge. Also, line search 'basic' is not a good choice usually. I would either use the default, or 'cp', or maybe 'nleqerr' that Patrick added.
Thanks,
Matt
On Sat, Jun 26, 2021 at 10:14 PM Bhavesh Shrimali < bhavesh.shrimali@gmail.com> wrote:
Hi everyone,
I am trying to understand the usage of solver parameters that are passed to the NonlinearVariationalSolver in firedrake. In order to do some tests, I decided to reproduce some of my previous calculations done in FEniCS for nonlinear viscoelasticity. In old dolfin, I had the following for solver parameters to be passed to PETSc's SNES:
snes_solver_parameters = {"nonlinear_solver": "snes", "snes_solver": {"linear_solver": "lu", "maximum_iterations": 50, "report": True, "line_search": 'basic', "error_on_nonconvergence": False, "relative_tolerance": 1e-7, "absolute_tolerance": 1e-8}} ... ... solver = NonlinearVariationalSolver(problem) solver.parameters.update(snes_solver_parameters)
Now as I understand, in Firedrake the options passed to the solver are consistent with the corresponding names in PETSc (correct me if I am wrong), and therefore I did:
solver_params={ "snes_type": "newtonls", "ksp_type": "preonly", "pc_type": "lu", "snes_max_it": 100, "snes_rtol": 1.e-7, "snes_atol": 1.e-8, "snes_view": None, "ksp_atol": 1.e-10, # not needed? "ksp_rtol":1.e-8, "ksp_view": None, "snes_linesearch_type": "basic", "snes_no_convergence_test": 1 # this has no effect } solver = NonlinearVariationalSolver(problem, solver_parameters=solver_params)
However, for some reason, I am unable to get the solver to converge when trying to solve the exact same problem as FEniCS. I experimented with different quadrature degrees as well as tweaked the solver parameters for SNES to best match what I am doing in dolfin. I understand that the form compiler is different hence the precise code that is generated may be different (I am not an expert on code generation) even though the python scripts may look similar due to UFL syntax.
In fact the weird part is that the solver seems to converge fine till a certain point in my input (loading). You can see this in the plot below where I plot the final stress-stretch relation obtained at a point in my domain from dolfin and firedrake..
[image: firedrake_vs_fenics.png]
I have marked the abscissa beyond which the solver fails to converge in Firedrake. The complete codes are here: https://github.com/bhaveshshrimali/comparingFiredrakeFEniCS
I tried setting the option `snes_no_convergence_test` to 1/True but that had no effect as the SNES solver still diverges with ` DIVERGED_DTOL`. In the end I put the call to solve inside a try/catch block to at least be able to see the solution till that point. And it surprisingly is accurate till that point.
Would anyone have any clue of what could be going wrong? I wouldn't expect that you go through the entire code(s) as that could be fairly long, but any pointers would be helpful :)
I am using firedrake via Docker on windows.
Thanks
--
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener
https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
--
*Bhavesh Shrimali* *Ph.D. Student, **Civil and Environmental Engineering,*
*University of Illinois Urbana-Champaign*
*Champaign, IL, United States*
*Email : bshrima2@illinois.edu <bshrima2@illinois.edu>* *Web : **https://bhaveshshrimali.github.io/ <https://bhaveshshrimali.github.io/>*
--
*Bhavesh Shrimali* *Ph.D. Student, **Civil and Environmental Engineering,*
*University of Illinois Urbana-Champaign*
*Champaign, IL, United States*
*Email : bshrima2@illinois.edu <bshrima2@illinois.edu>* *Web : **https://bhaveshshrimali.github.io/ <https://bhaveshshrimali.github.io/>* _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener
https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
--
*Bhavesh Shrimali* *Ph.D. Student, **Civil and Environmental Engineering,*
*University of Illinois Urbana-Champaign*
*Champaign, IL, United States*
*Email : bshrima2@illinois.edu <bshrima2@illinois.edu>* *Web : **https://bhaveshshrimali.github.io/ <https://bhaveshshrimali.github.io/>*
-- *Bhavesh Shrimali* *Ph.D. Student, **Civil and Environmental Engineering,* *University of Illinois Urbana-Champaign* *Champaign, IL, United States* *Email : bshrima2@illinois.edu <bshrima2@illinois.edu>* *Web : **https://bhaveshshrimali.github.io/ <https://bhaveshshrimali.github.io/>*
On Sun, Jun 27, 2021 at 8:14 PM Bhavesh Shrimali <bhavesh.shrimali@gmail.com> wrote:
After some more digging I am still unable to get the same solver setup using the same set of parameters. This time trying with tsfc for the form compiler (an old fork though: git+ https://github.com/blechta/tsfc.git@2018.1.0). Even with the same quadrature degree and the same form I get different SNES function norms in the very first solving step. Note that there are no convergence issues in the first step but eventually there are later.
Though with different PETSc versions should I expect similar convergence behavior ?
The first SNES residual is just the norm of F(x) where x is the input vector. So if the input vectors are the same, then the equations are not. This is not yet connected with PETSc. Thanks, Matt
There is very little difference in my formulation since most of it is in UFL syntax. I will keep looking ... In the meantime, here are the outputs --
Firedrake:
0 SNES Function norm 4.630756653490e+03 1 SNES Function norm 6.610067614593e+02 2 SNES Function norm 8.595689678029e+00 3 SNES Function norm 7.177260814570e-03 4 SNES Function norm 7.839549826428e-09 Nonlinear firedrake_0_ solve converged due to CONVERGED_FNORM_RELATIVE iterations 4
Dolfin:
Solving nonlinear variational problem. 0 SNES Function norm 4.048156097083e+03 1 SNES Function norm 7.701194025455e+01 2 SNES Function norm 9.032060723182e-02 3 SNES Function norm 2.188411552504e-08 PETSc SNES solver converged in 3 iterations with convergence reason CONVERGED_FNORM_RELATIVE.
and the codes (I will try to reduce these further as they still are far from minimal (~ 30 lines)) Dolfin (petsc4py version: 3.12.0):
import numpy as np from dolfin import *
parameters["form_compiler"]["cpp_optimize"] = True parameters["form_compiler"]["quadrature_degree"] = 5 parameters["form_compiler"]["representation"] = "tsfc" parameters["linear_algebra_backend"] = "PETSc"
snes_solver_parameters = {"nonlinear_solver": "snes", "snes_solver": {"linear_solver": "lu", "maximum_iterations": 50, "report": True, "line_search": 'basic', "error_on_nonconvergence": True} }
def Psi(u, p, Cv): I = Identity(3) F = I + grad(u) C = F.T * F Cv_inv = inv(Cv) I1 = tr(C) JF = det(F)
Ie1 = inner(C, Cv_inv) Je = JF / sqrt(det(Cv))
Jstar = (lam1 + p + sqrt((p + lam1)**2 + 4 * lam1 * (mu1 + mu2))) / (2 * lam1)
psi_Eq = ((3**(1 - alph1)) / (2 * alph1)) * mu1 * (I1**alph1 - 3**alph1) +\ ((3**(1 - alph2)) / (2 * alph2)) * mu2 * (I1**alph2 - 3**alph2) -\ (mu1 + mu2) * ln(Jstar) + 0.5 * lam1 * (Jstar - 1)**2
psi_Neq = ((3**(1 - a1)) / (2 * a1)) * m1 * (Ie1**a1 - 3**a1) +\ ((3**(1 - a2)) / (2 * a2)) * m2 * (Ie1**a2 - 3**a2) -\ (m1 + m2) * ln(Je)
W_hat = (psi_Eq + p * (JF - Jstar) + psi_Neq) * dx
return W_hat
mesh = UnitCubeMesh(3, 3, 3) W1 = VectorElement('CG', mesh.ufl_cell(), 2) W2 = FiniteElement('CG', mesh.ufl_cell(), 1) W3 = TensorElement("DG", mesh.ufl_cell(), 0)
V = FunctionSpace(mesh, W1 * W2) V_Cv = FunctionSpace(mesh, W3) du, dp = TrialFunctions(V) dw = TrialFunction(V) v, ph = TestFunctions(V) wh = TestFunction(V) w = Function(V) u, p = split(w) Cv = Function(V_Cv) Cv.interpolate(Constant(np.eye(3)))
left_fc = CompiledSubDomain( "near(x[0], 0.) && on_boundary" ) right_fc = CompiledSubDomain( "near(x[0], 1.) && on_boundary" ) bottom_fc = CompiledSubDomain( "near(x[1], 0.) && on_boundary" ) back_fc = CompiledSubDomain( "near(x[2], 0.) && on_boundary" )
strtch = Constant(0.0)
bc_right = DirichletBC(V.sub(0).sub(0), strtch, right_fc) bc_left = DirichletBC(V.sub(0).sub(0), Constant(0.0), left_fc) bc_bottom = DirichletBC(V.sub(0).sub(1), Constant(0.0), bottom_fc) bc_back = DirichletBC(V.sub(0).sub(2), Constant(0.0), back_fc)
bcs = [bc_right, bc_left, bc_bottom, bc_back] # Material Properties VHB4910 mu1 = Constant(13.54 * 1e3) mu2 = Constant(1.08 * 1e3) lam1 = Constant((mu1 + mu2) * 10**3) alph1 = Constant(1.) alph2 = Constant(-2.474) m1 = Constant(5.42 * 1e3) m2 = Constant(20.78 * 1e3) a1 = Constant(-10.) a2 = Constant(1.948) K1 = Constant(3507 * 1e3) K2 = Constant(1.0 * 1.e-12) bta1 = Constant(1.852) bta2 = Constant(0.26) eta0 = Constant(7014 * 1e3) etaInf = Constant(0.1 * 1e3)
S = derivative(Psi(u, p, Cv), w, wh) Jac = derivative(S, w, dw) problem = NonlinearVariationalProblem(S, w, bcs, J=Jac) solver = NonlinearVariationalSolver(problem) solver.parameters.update(snes_solver_parameters)
lamdot = 0.01 # Loading Rate Tfinal = 2 * (2 / lamdot) n_count = 101 dt = Tfinal / (n_count - 1)
timeVals = np.linspace(0, Tfinal, n_count)
stretchVals = np.hstack((lamdot * timeVals[:len(timeVals) // 2], lamdot * (-timeVals[len(timeVals) // 2:] + 2 * timeVals[len(timeVals) // 2]))) strtch.assign(stretchVals[1]) solver.solve() u, p = w.split()
Firedrake (petsc4py version: 3.15.1)
import numpy as np from firedrake import *
parameters["form_compiler"]["quadrature_degree"] = 5
snes_solver_parameters ={"snes_type": "newtonls", "ksp_type": "preonly", "pc_type": "lu", "snes_max_it": 50, "snes_linesearch_type": "basic", "snes_monitor": None, "snes_converged_reason": None}
def Psi(u, p, Cv): I = Identity(3) F = I + grad(u) C = F.T * F Cv_inv = inv(Cv) I1 = tr(C) JF = det(F)
Ie1 = inner(C, Cv_inv) Je = JF / sqrt(det(Cv))
Jstar = (lam1 + p + sqrt((p + lam1)**2 + 4 * lam1 * (mu1 + mu2))) / (2 * lam1)
psi_Eq = ((3**(1 - alph1)) / (2 * alph1)) * mu1 * (I1**alph1 - 3**alph1) +\ ((3**(1 - alph2)) / (2 * alph2)) * mu2 * (I1**alph2 - 3**alph2) -\ (mu1 + mu2) * ln(Jstar) + 0.5 * lam1 * (Jstar - 1)**2
psi_Neq = ((3**(1 - a1)) / (2 * a1)) * m1 * (Ie1**a1 - 3**a1) +\ ((3**(1 - a2)) / (2 * a2)) * m2 * (Ie1**a2 - 3**a2) -\ (m1 + m2) * ln(Je)
W_hat = (psi_Eq + p * (JF - Jstar) + psi_Neq) * dx
return W_hat
mesh = UnitCubeMesh(3, 3, 3) W1 = VectorElement('CG', mesh.ufl_cell(), 2) W2 = FiniteElement('CG', mesh.ufl_cell(), 1) W3 = TensorElement("DG", mesh.ufl_cell(), 0)
V = FunctionSpace(mesh, W1 * W2) V_Cv = FunctionSpace(mesh, W3) du, dp = TrialFunctions(V) dw = TrialFunction(V) v, ph = TestFunctions(V) wh = TestFunction(V) w = Function(V) u, p = split(w) Cv = Function(V_Cv) Cv.interpolate(Constant(np.eye(3)))
left_fc = 1 right_fc = 2 bottom_fc = 3 back_fc = 5
strtch = Constant(0.0)
bc_right = DirichletBC(V.sub(0).sub(0), strtch, right_fc) bc_left = DirichletBC(V.sub(0).sub(0), Constant(0.0), left_fc) bc_bottom = DirichletBC(V.sub(0).sub(1), Constant(0.0), bottom_fc) bc_back = DirichletBC(V.sub(0).sub(2), Constant(0.0), back_fc)
bcs = [bc_right, bc_left, bc_bottom, bc_back] # Material Properties VHB4910 mu1 = Constant(13.54 * 1e3) mu2 = Constant(1.08 * 1e3) lam1 = Constant((mu1 + mu2) * 10**3) alph1 = Constant(1.) alph2 = Constant(-2.474) m1 = Constant(5.42 * 1e3) m2 = Constant(20.78 * 1e3) a1 = Constant(-10.) a2 = Constant(1.948) K1 = Constant(3507 * 1e3) K2 = Constant(1.0 * 1.e-12) bta1 = Constant(1.852) bta2 = Constant(0.26) eta0 = Constant(7014 * 1e3) etaInf = Constant(0.1 * 1e3)
S = derivative(Psi(u, p, Cv), w, wh) Jac = derivative(S, w, dw) problem = NonlinearVariationalProblem(S, w, bcs, J=Jac) solver = NonlinearVariationalSolver(problem, solver_parameters=snes_solver_parameters) lamdot = 0.01 # Loading Rate Tfinal = 2 * (2 / lamdot) n_count = 101 dt = Tfinal / (n_count - 1)
timeVals = np.linspace(0, Tfinal, n_count)
stretchVals = np.hstack((lamdot * timeVals[:len(timeVals) // 2], lamdot * (-timeVals[len(timeVals) // 2:] + 2 * timeVals[len(timeVals) // 2]))) strtch.assign(stretchVals[1]) solver.solve() u, p = w.split()
Thanks for your help!
On Sun, 27 Jun 2021 at 13:21, Bhavesh Shrimali <bhavesh.shrimali@gmail.com> wrote:
Thanks Matt,
The only difference I see in the forms is that I am using `project` in Firedrake which I believe does a cell-wise projection for a DG-space whereas I use the `LocalSolver` in dolfin to do a cell-wise projection myself.
Even the SNES function norms that I am seeing aren't the same between FEniCS and Firedrake. Possible that I made some mistake in the formulation.
Thanks for pointing that out. I'll investigate this further today
Thanks, Bhavesh
On Sun, 27 Jun 2021 at 13:10, Matthew Knepley <knepley@gmail.com> wrote:
Something is wrong with the Newton system in the last step. You can see that the LU residual is steadily increasing. At the last step it is 10^-6. The iteration appears to be approaching a singular Jacobian. This seems like a formulation problem.
Thanks,
Matt
On Sun, Jun 27, 2021 at 11:24 AM Bhavesh Shrimali < bhavesh.shrimali@gmail.com> wrote:
I think I mistakenly clicked reply instead of reply all. Attached is my response:
--
Thanks Matt,
Thanks for the swift response. These are the logs with basic, backtracking, cp and nleqerr. If I don't specify 'snes_linesearch_type' it defaults to 'basic'. At least that's what I could gather from the log.
- `cp` seemingly diverges the fastest. Of course I haven't touched the other options. - `nleqerr` goes the furthest I believe, but not that far from `basic` or `bt`.
My concern was that since the forms are essentially the same, and I am passing "seemingly" the same parameters in FEniCS and Firedrake (both of which are using PETSc's SNES), what is it that is different which causes the solver to fail in one but not the other.
Edit: Just saw Ed's message. Indeed it defaults to `basic`
For definiteness this was the dictionary of solver parameters I am passing to NLVS in Firedrake
solver_params={ "snes_type": "newtonls", "ksp_type": "preonly", "pc_type": "lu", "snes_max_it": 100, # "snes_stol": 1.e-6, # "snes_rtol": 1.e-6, # "snes_atol": 1.e-6, "snes_view": None, # "ksp_atol": 1.e-10, # "ksp_rtol":1.e-8, "ksp_view": None, "snes_linesearch_type": "bt", "snes_converged_reason": None, "snes_monitor": None, "ksp_converged_reason": None, "ksp_monitor_true_residual": None }
On Sun, 27 Jun 2021 at 11:22, Bhavesh Shrimali < bhavesh.shrimali@gmail.com> wrote:
Thanks Matt,
Thanks for the swift response. These are the logs with basic, backtracking, cp and nleqerr. If I don't specify 'snes_linesearch_type' it defaults to 'basic'. At least that's what I could gather from the log.
- `cp` seemingly diverges the fastest. Of course I haven't touched the other options. - `nleqerr` goes the furthest I believe, but not that far from `basic` or `bt`.
My concern was that since the forms are essentially the same, and I am passing "seemingly" the same parameters in FEniCS and Firedrake (both of which are using PETSc's SNES), what is it that is different which causes the solver to fail in one but not the other.
Edit: Just saw Ed's message. Indeed it defaults to `basic`
For definiteness this was the dictionary of solver parameters I am passing to NLVS in Firedrake
solver_params={ "snes_type": "newtonls", "ksp_type": "preonly", "pc_type": "lu", "snes_max_it": 100, # "snes_stol": 1.e-6, # "snes_rtol": 1.e-6, # "snes_atol": 1.e-6, "snes_view": None, # "ksp_atol": 1.e-10, # "ksp_rtol":1.e-8, "ksp_view": None, "snes_linesearch_type": "bt", "snes_converged_reason": None, "snes_monitor": None, "ksp_converged_reason": None, "ksp_monitor_true_residual": None }
On Sun, 27 Jun 2021 at 06:00, Matthew Knepley <knepley@gmail.com> wrote:
Give
-snes_view -snes_converged_reason -snes_monitor -ksp_converged_reason -ksp_monitor_true_residual
and send all the output. That way I might be able to see what is happening when you do not converge. Also, line search 'basic' is not a good choice usually. I would either use the default, or 'cp', or maybe 'nleqerr' that Patrick added.
Thanks,
Matt
On Sat, Jun 26, 2021 at 10:14 PM Bhavesh Shrimali < bhavesh.shrimali@gmail.com> wrote:
> Hi everyone, > > I am trying to understand the usage of solver parameters that are > passed to the NonlinearVariationalSolver in firedrake. In order to do some > tests, I decided to reproduce some of my previous calculations done in > FEniCS for nonlinear viscoelasticity. In old dolfin, I had the following > for solver parameters to be passed to PETSc's SNES: > > snes_solver_parameters = {"nonlinear_solver": "snes", > "snes_solver": {"linear_solver": "lu", > "maximum_iterations": 50, > "report": True, > "line_search": 'basic', > "error_on_nonconvergence": False, > "relative_tolerance": 1e-7, > "absolute_tolerance": 1e-8}} > ... > ... > solver = NonlinearVariationalSolver(problem) > solver.parameters.update(snes_solver_parameters) > > Now as I understand, in Firedrake the options passed to the solver > are consistent with the corresponding names in PETSc (correct me if I am > wrong), and therefore I did: > > solver_params={ > "snes_type": "newtonls", > "ksp_type": "preonly", > "pc_type": "lu", > "snes_max_it": 100, > "snes_rtol": 1.e-7, > "snes_atol": 1.e-8, > "snes_view": None, > "ksp_atol": 1.e-10, # not needed? > "ksp_rtol":1.e-8, > "ksp_view": None, > "snes_linesearch_type": "basic", > "snes_no_convergence_test": 1 # this has no effect > } > solver = NonlinearVariationalSolver(problem, > solver_parameters=solver_params) > > However, for some reason, I am unable to get the solver to converge > when trying to solve the exact same problem as FEniCS. I experimented with > different quadrature degrees as well as tweaked the solver parameters for > SNES to best match what I am doing in dolfin. I understand that the form > compiler is different hence the precise code that is generated may be > different (I am not an expert on code generation) even though the python > scripts may look similar due to UFL syntax. > > In fact the weird part is that the solver seems to converge fine > till a certain point in my input (loading). You can see this in the plot > below where I plot the final stress-stretch relation obtained at a point in > my domain from dolfin and firedrake.. > > [image: firedrake_vs_fenics.png] > > > I have marked the abscissa beyond which the solver fails to converge > in Firedrake. The complete codes are here: > https://github.com/bhaveshshrimali/comparingFiredrakeFEniCS > > I tried setting the option `snes_no_convergence_test` to 1/True but > that had no effect as the SNES solver still diverges with ` > DIVERGED_DTOL`. In the end I put the call to solve inside a > try/catch block to at least be able to see the solution till that point. > And it surprisingly is accurate till that point. > > Would anyone have any clue of what could be going wrong? I wouldn't > expect that you go through the entire code(s) as that could be fairly long, > but any pointers would be helpful :) > > I am using firedrake via Docker on windows. > > Thanks > > -- > > _______________________________________________ > firedrake mailing list > firedrake@imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake >
-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener
https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
--
*Bhavesh Shrimali* *Ph.D. Student, **Civil and Environmental Engineering,*
*University of Illinois Urbana-Champaign*
*Champaign, IL, United States*
*Email : bshrima2@illinois.edu <bshrima2@illinois.edu>* *Web : **https://bhaveshshrimali.github.io/ <https://bhaveshshrimali.github.io/>*
--
*Bhavesh Shrimali* *Ph.D. Student, **Civil and Environmental Engineering,*
*University of Illinois Urbana-Champaign*
*Champaign, IL, United States*
*Email : bshrima2@illinois.edu <bshrima2@illinois.edu>* *Web : **https://bhaveshshrimali.github.io/ <https://bhaveshshrimali.github.io/>* _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener
https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
--
*Bhavesh Shrimali* *Ph.D. Student, **Civil and Environmental Engineering,*
*University of Illinois Urbana-Champaign*
*Champaign, IL, United States*
*Email : bshrima2@illinois.edu <bshrima2@illinois.edu>* *Web : **https://bhaveshshrimali.github.io/ <https://bhaveshshrimali.github.io/>*
--
*Bhavesh Shrimali* *Ph.D. Student, **Civil and Environmental Engineering,*
*University of Illinois Urbana-Champaign*
*Champaign, IL, United States*
*Email : bshrima2@illinois.edu <bshrima2@illinois.edu>* *Web : **https://bhaveshshrimali.github.io/ <https://bhaveshshrimali.github.io/>*
-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
Thanks Matt, That makes sense. The first SNES residual is close but not the same for dolfin/firedrake but the meshes (and function spaces) are the same, so I would guess that the input vectors are the same. Note that if I don't prescribe anything, and the system is in the initial state I get the same norms. But then there is nothing to solve for -- hence I should get the exact same norms. 0 SNES Function norm 3.751493749134e+03 1 SNES Function norm 1.865848147286e-06 Nonlinear firedrake_0_ solve converged due to CONVERGED_FNORM_RELATIVE iterations 1 Solving nonlinear variational problem. 0 SNES Function norm 3.751493749134e+03 1 SNES Function norm 1.865848023631e-06 PETSc SNES solver converged in 1 iterations with convergence reason CONVERGED_FNORM_RELATIVE. All this may not even be relevant for the ultimate goal but I just wanted to make sure if I have implemented everything as it should be -- more of as a sanity check. And it looks like that might not be true. On Sun, 27 Jun 2021 at 21:07, Matthew Knepley <knepley@gmail.com> wrote:
On Sun, Jun 27, 2021 at 8:14 PM Bhavesh Shrimali < bhavesh.shrimali@gmail.com> wrote:
After some more digging I am still unable to get the same solver setup using the same set of parameters. This time trying with tsfc for the form compiler (an old fork though: git+ https://github.com/blechta/tsfc.git@2018.1.0). Even with the same quadrature degree and the same form I get different SNES function norms in the very first solving step. Note that there are no convergence issues in the first step but eventually there are later.
Though with different PETSc versions should I expect similar convergence behavior ?
The first SNES residual is just the norm of F(x) where x is the input vector. So if the input vectors are the same, then the equations are not. This is not yet connected with PETSc.
Thanks,
Matt
There is very little difference in my formulation since most of it is in UFL syntax. I will keep looking ... In the meantime, here are the outputs --
Firedrake:
0 SNES Function norm 4.630756653490e+03 1 SNES Function norm 6.610067614593e+02 2 SNES Function norm 8.595689678029e+00 3 SNES Function norm 7.177260814570e-03 4 SNES Function norm 7.839549826428e-09 Nonlinear firedrake_0_ solve converged due to CONVERGED_FNORM_RELATIVE iterations 4
Dolfin:
Solving nonlinear variational problem. 0 SNES Function norm 4.048156097083e+03 1 SNES Function norm 7.701194025455e+01 2 SNES Function norm 9.032060723182e-02 3 SNES Function norm 2.188411552504e-08 PETSc SNES solver converged in 3 iterations with convergence reason CONVERGED_FNORM_RELATIVE.
and the codes (I will try to reduce these further as they still are far from minimal (~ 30 lines)) Dolfin (petsc4py version: 3.12.0):
import numpy as np from dolfin import *
parameters["form_compiler"]["cpp_optimize"] = True parameters["form_compiler"]["quadrature_degree"] = 5 parameters["form_compiler"]["representation"] = "tsfc" parameters["linear_algebra_backend"] = "PETSc"
snes_solver_parameters = {"nonlinear_solver": "snes", "snes_solver": {"linear_solver": "lu", "maximum_iterations": 50, "report": True, "line_search": 'basic', "error_on_nonconvergence": True} }
def Psi(u, p, Cv): I = Identity(3) F = I + grad(u) C = F.T * F Cv_inv = inv(Cv) I1 = tr(C) JF = det(F)
Ie1 = inner(C, Cv_inv) Je = JF / sqrt(det(Cv))
Jstar = (lam1 + p + sqrt((p + lam1)**2 + 4 * lam1 * (mu1 + mu2))) / (2 * lam1)
psi_Eq = ((3**(1 - alph1)) / (2 * alph1)) * mu1 * (I1**alph1 - 3**alph1) +\ ((3**(1 - alph2)) / (2 * alph2)) * mu2 * (I1**alph2 - 3**alph2) -\ (mu1 + mu2) * ln(Jstar) + 0.5 * lam1 * (Jstar - 1)**2
psi_Neq = ((3**(1 - a1)) / (2 * a1)) * m1 * (Ie1**a1 - 3**a1) +\ ((3**(1 - a2)) / (2 * a2)) * m2 * (Ie1**a2 - 3**a2) -\ (m1 + m2) * ln(Je)
W_hat = (psi_Eq + p * (JF - Jstar) + psi_Neq) * dx
return W_hat
mesh = UnitCubeMesh(3, 3, 3) W1 = VectorElement('CG', mesh.ufl_cell(), 2) W2 = FiniteElement('CG', mesh.ufl_cell(), 1) W3 = TensorElement("DG", mesh.ufl_cell(), 0)
V = FunctionSpace(mesh, W1 * W2) V_Cv = FunctionSpace(mesh, W3) du, dp = TrialFunctions(V) dw = TrialFunction(V) v, ph = TestFunctions(V) wh = TestFunction(V) w = Function(V) u, p = split(w) Cv = Function(V_Cv) Cv.interpolate(Constant(np.eye(3)))
left_fc = CompiledSubDomain( "near(x[0], 0.) && on_boundary" ) right_fc = CompiledSubDomain( "near(x[0], 1.) && on_boundary" ) bottom_fc = CompiledSubDomain( "near(x[1], 0.) && on_boundary" ) back_fc = CompiledSubDomain( "near(x[2], 0.) && on_boundary" )
strtch = Constant(0.0)
bc_right = DirichletBC(V.sub(0).sub(0), strtch, right_fc) bc_left = DirichletBC(V.sub(0).sub(0), Constant(0.0), left_fc) bc_bottom = DirichletBC(V.sub(0).sub(1), Constant(0.0), bottom_fc) bc_back = DirichletBC(V.sub(0).sub(2), Constant(0.0), back_fc)
bcs = [bc_right, bc_left, bc_bottom, bc_back] # Material Properties VHB4910 mu1 = Constant(13.54 * 1e3) mu2 = Constant(1.08 * 1e3) lam1 = Constant((mu1 + mu2) * 10**3) alph1 = Constant(1.) alph2 = Constant(-2.474) m1 = Constant(5.42 * 1e3) m2 = Constant(20.78 * 1e3) a1 = Constant(-10.) a2 = Constant(1.948) K1 = Constant(3507 * 1e3) K2 = Constant(1.0 * 1.e-12) bta1 = Constant(1.852) bta2 = Constant(0.26) eta0 = Constant(7014 * 1e3) etaInf = Constant(0.1 * 1e3)
S = derivative(Psi(u, p, Cv), w, wh) Jac = derivative(S, w, dw) problem = NonlinearVariationalProblem(S, w, bcs, J=Jac) solver = NonlinearVariationalSolver(problem) solver.parameters.update(snes_solver_parameters)
lamdot = 0.01 # Loading Rate Tfinal = 2 * (2 / lamdot) n_count = 101 dt = Tfinal / (n_count - 1)
timeVals = np.linspace(0, Tfinal, n_count)
stretchVals = np.hstack((lamdot * timeVals[:len(timeVals) // 2], lamdot * (-timeVals[len(timeVals) // 2:] + 2 * timeVals[len(timeVals) // 2]))) strtch.assign(stretchVals[1]) solver.solve() u, p = w.split()
Firedrake (petsc4py version: 3.15.1)
import numpy as np from firedrake import *
parameters["form_compiler"]["quadrature_degree"] = 5
snes_solver_parameters ={"snes_type": "newtonls", "ksp_type": "preonly", "pc_type": "lu", "snes_max_it": 50, "snes_linesearch_type": "basic", "snes_monitor": None, "snes_converged_reason": None}
def Psi(u, p, Cv): I = Identity(3) F = I + grad(u) C = F.T * F Cv_inv = inv(Cv) I1 = tr(C) JF = det(F)
Ie1 = inner(C, Cv_inv) Je = JF / sqrt(det(Cv))
Jstar = (lam1 + p + sqrt((p + lam1)**2 + 4 * lam1 * (mu1 + mu2))) / (2 * lam1)
psi_Eq = ((3**(1 - alph1)) / (2 * alph1)) * mu1 * (I1**alph1 - 3**alph1) +\ ((3**(1 - alph2)) / (2 * alph2)) * mu2 * (I1**alph2 - 3**alph2) -\ (mu1 + mu2) * ln(Jstar) + 0.5 * lam1 * (Jstar - 1)**2
psi_Neq = ((3**(1 - a1)) / (2 * a1)) * m1 * (Ie1**a1 - 3**a1) +\ ((3**(1 - a2)) / (2 * a2)) * m2 * (Ie1**a2 - 3**a2) -\ (m1 + m2) * ln(Je)
W_hat = (psi_Eq + p * (JF - Jstar) + psi_Neq) * dx
return W_hat
mesh = UnitCubeMesh(3, 3, 3) W1 = VectorElement('CG', mesh.ufl_cell(), 2) W2 = FiniteElement('CG', mesh.ufl_cell(), 1) W3 = TensorElement("DG", mesh.ufl_cell(), 0)
V = FunctionSpace(mesh, W1 * W2) V_Cv = FunctionSpace(mesh, W3) du, dp = TrialFunctions(V) dw = TrialFunction(V) v, ph = TestFunctions(V) wh = TestFunction(V) w = Function(V) u, p = split(w) Cv = Function(V_Cv) Cv.interpolate(Constant(np.eye(3)))
left_fc = 1 right_fc = 2 bottom_fc = 3 back_fc = 5
strtch = Constant(0.0)
bc_right = DirichletBC(V.sub(0).sub(0), strtch, right_fc) bc_left = DirichletBC(V.sub(0).sub(0), Constant(0.0), left_fc) bc_bottom = DirichletBC(V.sub(0).sub(1), Constant(0.0), bottom_fc) bc_back = DirichletBC(V.sub(0).sub(2), Constant(0.0), back_fc)
bcs = [bc_right, bc_left, bc_bottom, bc_back] # Material Properties VHB4910 mu1 = Constant(13.54 * 1e3) mu2 = Constant(1.08 * 1e3) lam1 = Constant((mu1 + mu2) * 10**3) alph1 = Constant(1.) alph2 = Constant(-2.474) m1 = Constant(5.42 * 1e3) m2 = Constant(20.78 * 1e3) a1 = Constant(-10.) a2 = Constant(1.948) K1 = Constant(3507 * 1e3) K2 = Constant(1.0 * 1.e-12) bta1 = Constant(1.852) bta2 = Constant(0.26) eta0 = Constant(7014 * 1e3) etaInf = Constant(0.1 * 1e3)
S = derivative(Psi(u, p, Cv), w, wh) Jac = derivative(S, w, dw) problem = NonlinearVariationalProblem(S, w, bcs, J=Jac) solver = NonlinearVariationalSolver(problem, solver_parameters=snes_solver_parameters) lamdot = 0.01 # Loading Rate Tfinal = 2 * (2 / lamdot) n_count = 101 dt = Tfinal / (n_count - 1)
timeVals = np.linspace(0, Tfinal, n_count)
stretchVals = np.hstack((lamdot * timeVals[:len(timeVals) // 2], lamdot * (-timeVals[len(timeVals) // 2:] + 2 * timeVals[len(timeVals) // 2]))) strtch.assign(stretchVals[1]) solver.solve() u, p = w.split()
Thanks for your help!
On Sun, 27 Jun 2021 at 13:21, Bhavesh Shrimali < bhavesh.shrimali@gmail.com> wrote:
Thanks Matt,
The only difference I see in the forms is that I am using `project` in Firedrake which I believe does a cell-wise projection for a DG-space whereas I use the `LocalSolver` in dolfin to do a cell-wise projection myself.
Even the SNES function norms that I am seeing aren't the same between FEniCS and Firedrake. Possible that I made some mistake in the formulation.
Thanks for pointing that out. I'll investigate this further today
Thanks, Bhavesh
On Sun, 27 Jun 2021 at 13:10, Matthew Knepley <knepley@gmail.com> wrote:
Something is wrong with the Newton system in the last step. You can see that the LU residual is steadily increasing. At the last step it is 10^-6. The iteration appears to be approaching a singular Jacobian. This seems like a formulation problem.
Thanks,
Matt
On Sun, Jun 27, 2021 at 11:24 AM Bhavesh Shrimali < bhavesh.shrimali@gmail.com> wrote:
I think I mistakenly clicked reply instead of reply all. Attached is my response:
--
Thanks Matt,
Thanks for the swift response. These are the logs with basic, backtracking, cp and nleqerr. If I don't specify 'snes_linesearch_type' it defaults to 'basic'. At least that's what I could gather from the log.
- `cp` seemingly diverges the fastest. Of course I haven't touched the other options. - `nleqerr` goes the furthest I believe, but not that far from `basic` or `bt`.
My concern was that since the forms are essentially the same, and I am passing "seemingly" the same parameters in FEniCS and Firedrake (both of which are using PETSc's SNES), what is it that is different which causes the solver to fail in one but not the other.
Edit: Just saw Ed's message. Indeed it defaults to `basic`
For definiteness this was the dictionary of solver parameters I am passing to NLVS in Firedrake
solver_params={ "snes_type": "newtonls", "ksp_type": "preonly", "pc_type": "lu", "snes_max_it": 100, # "snes_stol": 1.e-6, # "snes_rtol": 1.e-6, # "snes_atol": 1.e-6, "snes_view": None, # "ksp_atol": 1.e-10, # "ksp_rtol":1.e-8, "ksp_view": None, "snes_linesearch_type": "bt", "snes_converged_reason": None, "snes_monitor": None, "ksp_converged_reason": None, "ksp_monitor_true_residual": None }
On Sun, 27 Jun 2021 at 11:22, Bhavesh Shrimali < bhavesh.shrimali@gmail.com> wrote:
Thanks Matt,
Thanks for the swift response. These are the logs with basic, backtracking, cp and nleqerr. If I don't specify 'snes_linesearch_type' it defaults to 'basic'. At least that's what I could gather from the log.
- `cp` seemingly diverges the fastest. Of course I haven't touched the other options. - `nleqerr` goes the furthest I believe, but not that far from `basic` or `bt`.
My concern was that since the forms are essentially the same, and I am passing "seemingly" the same parameters in FEniCS and Firedrake (both of which are using PETSc's SNES), what is it that is different which causes the solver to fail in one but not the other.
Edit: Just saw Ed's message. Indeed it defaults to `basic`
For definiteness this was the dictionary of solver parameters I am passing to NLVS in Firedrake
solver_params={ "snes_type": "newtonls", "ksp_type": "preonly", "pc_type": "lu", "snes_max_it": 100, # "snes_stol": 1.e-6, # "snes_rtol": 1.e-6, # "snes_atol": 1.e-6, "snes_view": None, # "ksp_atol": 1.e-10, # "ksp_rtol":1.e-8, "ksp_view": None, "snes_linesearch_type": "bt", "snes_converged_reason": None, "snes_monitor": None, "ksp_converged_reason": None, "ksp_monitor_true_residual": None }
On Sun, 27 Jun 2021 at 06:00, Matthew Knepley <knepley@gmail.com> wrote:
> Give > > -snes_view -snes_converged_reason -snes_monitor > -ksp_converged_reason -ksp_monitor_true_residual > > and send all the output. That way I might be able to see what is > happening when you do not converge. Also, > line search 'basic' is not a good choice usually. I would either use > the default, or 'cp', or maybe 'nleqerr' that > Patrick added. > > Thanks, > > Matt > > On Sat, Jun 26, 2021 at 10:14 PM Bhavesh Shrimali < > bhavesh.shrimali@gmail.com> wrote: > >> Hi everyone, >> >> I am trying to understand the usage of solver parameters that are >> passed to the NonlinearVariationalSolver in firedrake. In order to do some >> tests, I decided to reproduce some of my previous calculations done in >> FEniCS for nonlinear viscoelasticity. In old dolfin, I had the following >> for solver parameters to be passed to PETSc's SNES: >> >> snes_solver_parameters = {"nonlinear_solver": "snes", >> "snes_solver": {"linear_solver": "lu", >> "maximum_iterations": 50, >> "report": True, >> "line_search": 'basic', >> "error_on_nonconvergence": False, >> "relative_tolerance": 1e-7, >> "absolute_tolerance": 1e-8}} >> ... >> ... >> solver = NonlinearVariationalSolver(problem) >> solver.parameters.update(snes_solver_parameters) >> >> Now as I understand, in Firedrake the options passed to the solver >> are consistent with the corresponding names in PETSc (correct me if I am >> wrong), and therefore I did: >> >> solver_params={ >> "snes_type": "newtonls", >> "ksp_type": "preonly", >> "pc_type": "lu", >> "snes_max_it": 100, >> "snes_rtol": 1.e-7, >> "snes_atol": 1.e-8, >> "snes_view": None, >> "ksp_atol": 1.e-10, # not needed? >> "ksp_rtol":1.e-8, >> "ksp_view": None, >> "snes_linesearch_type": "basic", >> "snes_no_convergence_test": 1 # this has no effect >> } >> solver = NonlinearVariationalSolver(problem, >> solver_parameters=solver_params) >> >> However, for some reason, I am unable to get the solver to converge >> when trying to solve the exact same problem as FEniCS. I experimented with >> different quadrature degrees as well as tweaked the solver parameters for >> SNES to best match what I am doing in dolfin. I understand that the form >> compiler is different hence the precise code that is generated may be >> different (I am not an expert on code generation) even though the python >> scripts may look similar due to UFL syntax. >> >> In fact the weird part is that the solver seems to converge fine >> till a certain point in my input (loading). You can see this in the plot >> below where I plot the final stress-stretch relation obtained at a point in >> my domain from dolfin and firedrake.. >> >> [image: firedrake_vs_fenics.png] >> >> >> I have marked the abscissa beyond which the solver fails to >> converge in Firedrake. The complete codes are here: >> https://github.com/bhaveshshrimali/comparingFiredrakeFEniCS >> >> I tried setting the option `snes_no_convergence_test` to 1/True >> but that had no effect as the SNES solver still diverges with ` >> DIVERGED_DTOL`. In the end I put the call to solve inside a >> try/catch block to at least be able to see the solution till that point. >> And it surprisingly is accurate till that point. >> >> Would anyone have any clue of what could be going wrong? I wouldn't >> expect that you go through the entire code(s) as that could be fairly long, >> but any pointers would be helpful :) >> >> I am using firedrake via Docker on windows. >> >> Thanks >> >> -- >> >> _______________________________________________ >> firedrake mailing list >> firedrake@imperial.ac.uk >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >> > > > -- > What most experimenters take for granted before they begin their > experiments is infinitely more interesting than any results to which their > experiments lead. > -- Norbert Wiener > > https://www.cse.buffalo.edu/~knepley/ > <http://www.cse.buffalo.edu/~knepley/> >
--
*Bhavesh Shrimali* *Ph.D. Student, **Civil and Environmental Engineering,*
*University of Illinois Urbana-Champaign*
*Champaign, IL, United States*
*Email : bshrima2@illinois.edu <bshrima2@illinois.edu>* *Web : **https://bhaveshshrimali.github.io/ <https://bhaveshshrimali.github.io/>*
--
*Bhavesh Shrimali* *Ph.D. Student, **Civil and Environmental Engineering,*
*University of Illinois Urbana-Champaign*
*Champaign, IL, United States*
*Email : bshrima2@illinois.edu <bshrima2@illinois.edu>* *Web : **https://bhaveshshrimali.github.io/ <https://bhaveshshrimali.github.io/>* _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener
https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
--
*Bhavesh Shrimali* *Ph.D. Student, **Civil and Environmental Engineering,*
*University of Illinois Urbana-Champaign*
*Champaign, IL, United States*
*Email : bshrima2@illinois.edu <bshrima2@illinois.edu>* *Web : **https://bhaveshshrimali.github.io/ <https://bhaveshshrimali.github.io/>*
--
*Bhavesh Shrimali* *Ph.D. Student, **Civil and Environmental Engineering,*
*University of Illinois Urbana-Champaign*
*Champaign, IL, United States*
*Email : bshrima2@illinois.edu <bshrima2@illinois.edu>* *Web : **https://bhaveshshrimali.github.io/ <https://bhaveshshrimali.github.io/>*
-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener
https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
-- *Bhavesh Shrimali* *Ph.D. Student, **Civil and Environmental Engineering,* *University of Illinois Urbana-Champaign* *Champaign, IL, United States* *Email : bshrima2@illinois.edu <bshrima2@illinois.edu>* *Web : **https://bhaveshshrimali.github.io/ <https://bhaveshshrimali.github.io/>*
So after exploring a bit more, I am inclined to think that it is indeed the final equations (the nonlinear residual) that are different between dolfin and Firedrake and nothing wrong with the solver per se. Moreover I am also noticing weird errors with respect to function evaluations. I am getting NaN's in places where I shouldn't. Maybe there are floating point errors but I am not sure where. Consider, for instance, the function: def Eveq_Cv(u, Cv): I = Identity(3) # Identity Tensor F = I + grad(u) C = F.T * F Cv_inv = inv(Cv) Iv1 = tr(Cv) Ie1 = inner(C, Cv_inv) Ie2 = 0.5 * (Ie1**2. - inner(Cv_inv * C, C * Cv_inv)) c_neq = m1 * (Ie1 / 3.)**(a1 - 1.) + m2 * (Ie1 / 3.)**(a2 - 1) J2Neq = (Ie1**2 / 3. - Ie2) * c_neq**2. etaK = etaInf + (eta0 - etaInf + K1 * (Iv1**bta1 - 3.**bta1)) / (1 + (K2 * J2Neq)**bta2) G = (c_neq / etaK) * (C - Ie1 * Cv / 3) G = project(G, V_Cv) print(f"Int Ie1: {assemble(Ie1*dx)}") print(f"Int Ie2: {assemble(Ie2*dx)}") print(f"Int Iv1: {assemble(Iv1*dx)}") print(f"Int J2: {assemble(J2Neq*dx)}") print(f"Int etaK: {assemble(etaK*dx)}") print(f"Norm G: {norm(G)}") return G in my code. For initial conditions (u=(0,0,0), Cv=Identity matrix) I should get I1 = Ie1 = Iv1 = 3, J2 = 0.0, etaK = eta0 and G = 0 (matrix). However I get the following outputs (etaK and G are NaN) 0 SNES Function norm 3.751493749134e+03 1 SNES Function norm 1.865847686176e-06 ||u||: 4.060936555993801e-10, ||p||: 14619.999995628163, ||E||: -7.361361589262569e-12 Int Ie1: 2.999999998010551 # 3.0 Int Ie2: 2.9999999960211134 # 3.0 Int Iv1: 2.999999999999996 # 3.0 Int J2: 1.8368572690218592e-07 # 0.0 Int etaK: nan # eta0 = 7014 * 1e3 Norm G: nan # 0.0 whereas the same operations when I carry out in dolfin I get the following output which seems more reasonable.. ||u||: 4.0609362990701997e-10, ||p||: 14619.99999562816, ||E||: 6.166784849648951e-13 Int Ie1: 2.9999999980105647 Int Ie2: 2.999999996021099 Int Iv1: 3.0000000000000075 Int J2: 0.0 Int etaK: 7014000.000000014 Norm G: 3.0338331632858357e-12 The complete code for reproducing this: import numpy as np from firedrake import * parameters["form_compiler"]["quadrature_degree"]=7 def Eveq_Cv(u, Cv): I = Identity(3) # Identity Tensor F = I + grad(u) C = F.T * F Cv_inv = inv(Cv) Iv1 = tr(Cv) Ie1 = inner(C, Cv_inv) Ie2 = 0.5 * (Ie1**2. - inner(Cv_inv * C, C * Cv_inv)) c_neq = m1 * (Ie1 / 3.)**(a1 - 1.) + m2 * (Ie1 / 3.)**(a2 - 1) J2Neq = (Ie1**2 / 3. - Ie2) * c_neq**2. etaK = etaInf + (eta0 - etaInf + K1 * (Iv1**bta1 - 3.**bta1)) / (1 + (K2 * J2Neq)**bta2) G = (c_neq / etaK) * (C - Ie1 * Cv / 3) G = project(G, V_Cv) print(f"Int Ie1: {assemble(Ie1*dx)}") print(f"Int Ie2: {assemble(Ie2*dx)}") print(f"Int Iv1: {assemble(Iv1*dx)}") print(f"Int J2: {assemble(J2Neq*dx)}") print(f"Int etaK: {assemble(etaK*dx)}") print(f"Norm G: {norm(G)}") return G def Psi(u, p, Cv): # Kinematics I = Identity(3) # Identity Tensor F = I + grad(u) C = F.T * F Cv_inv = inv(Cv) I1 = tr(C) JF = det(F) Ie1 = inner(C, Cv_inv) Je = JF / sqrt(det(Cv)) Jstar = (lam1 + p + sqrt((p + lam1)**2 + 4 *\ lam1 * (mu1 + mu2))) / (2 * lam1) psi_Eq = ((3**(1 - alph1)) / (2 * alph1)) * mu1 * (I1**alph1 - 3**alph1) +\ ((3**(1 - alph2)) / (2 * alph2)) * mu2 * (I1**alph2 - 3**alph2) -\ (mu1 + mu2) * ln(Jstar) + 0.5 * lam1 * (Jstar - 1)**2 psi_Neq = ((3**(1 - a1)) / (2 * a1)) * m1 * (Ie1**a1 - 3**a1) +\ ((3**(1 - a2)) / (2 * a2)) * m2 * (Ie1**a2 - 3**a2) -\ (m1 + m2) * ln(Je) W_hat = (psi_Eq + p * (JF - Jstar) + psi_Neq) * dx return W_hat mesh = UnitCubeMesh(3, 3, 3) W1 = VectorElement('CG', mesh.ufl_cell(), 2) W2 = FiniteElement('CG', mesh.ufl_cell(), 1) W3 = TensorElement("DG", mesh.ufl_cell(), 0) # for Cv V = FunctionSpace(mesh, W1 * W2) V_Cv = FunctionSpace(mesh, W3) w = Function(V) u, p = split(w) wh = TestFunction(V) dw = TrialFunction(V) Cv = Function(V_Cv) # u.interpolate(Constant([0., 0., 0.])) Cv.interpolate(Constant(np.eye(3))) mu1 = Constant(13.54 * 1e3) mu2 = Constant(1.08 * 1e3) lam1 = Constant((mu1 + mu2) * 10**3) # make this value very high alph1 = Constant(1.) alph2 = Constant(-2.474) m1 = Constant(5.42 * 1e3) m2 = Constant(20.78 * 1e3) a1 = Constant(-10.) a2 = Constant(1.948) K1 = Constant(3507 * 1e3) K2 = Constant(1.0 * 1.e-12) bta1 = Constant(1.852) bta2 = Constant(0.26) eta0 = Constant(7014 * 1e3) etaInf = Constant(0.1 * 1e3) lamdot = 0.01 Tfinal = 2 * (2. / lamdot) n_count = 101 dt = Tfinal / (n_count - 1) timeVals = np.linspace(0, Tfinal, n_count) stretchVals = np.hstack((lamdot * timeVals[:len(timeVals) // 2], lamdot * (-timeVals[len(timeVals) // 2:] + 2 * timeVals[len(timeVals) // 2]))) left_fc = 1 right_fc = 2 bottom_fc = 3 back_fc = 5 strtch = Constant(0.0) bc_right = DirichletBC(V.sub(0).sub(0), strtch, right_fc) bc_left = DirichletBC(V.sub(0).sub(0), Constant(0.0), left_fc) bc_bottom = DirichletBC(V.sub(0).sub(1), Constant(0.0), bottom_fc) bc_back = DirichletBC(V.sub(0).sub(2), Constant(0.0), back_fc) bcs = [bc_right, bc_left, bc_bottom, bc_back] S = derivative(Psi(u, p, Cv), w, wh) Jac = derivative(S, w, dw) problem = NonlinearVariationalProblem(S, w, bcs, J=Jac) solver_params={ "snes_type": "newtonls", "ksp_type": "preonly", "pc_type": "lu", "snes_max_it": 100, "snes_rtol": 1.e-7, "snes_atol": 1.e-8, "snes_monitor": None, "snes_linesearch_type": "basic" } solver = NonlinearVariationalSolver(problem, solver_parameters=solver_params) solver.solve() u, p = w.split() print( f"||u||: {norm(u)},\ ||p||: {norm(p)},\ ||E||: {assemble(Psi(u, p, Cv))}" ) Eveq_Cv(u, Cv) On Sun, 27 Jun 2021 at 21:28, Bhavesh Shrimali <bhavesh.shrimali@gmail.com> wrote:
Thanks Matt,
That makes sense. The first SNES residual is close but not the same for dolfin/firedrake but the meshes (and function spaces) are the same, so I would guess that the input vectors are the same.
Note that if I don't prescribe anything, and the system is in the initial state I get the same norms. But then there is nothing to solve for -- hence I should get the exact same norms.
0 SNES Function norm 3.751493749134e+03 1 SNES Function norm 1.865848147286e-06 Nonlinear firedrake_0_ solve converged due to CONVERGED_FNORM_RELATIVE iterations 1
Solving nonlinear variational problem. 0 SNES Function norm 3.751493749134e+03 1 SNES Function norm 1.865848023631e-06 PETSc SNES solver converged in 1 iterations with convergence reason CONVERGED_FNORM_RELATIVE.
All this may not even be relevant for the ultimate goal but I just wanted to make sure if I have implemented everything as it should be -- more of as a sanity check. And it looks like that might not be true.
On Sun, 27 Jun 2021 at 21:07, Matthew Knepley <knepley@gmail.com> wrote:
On Sun, Jun 27, 2021 at 8:14 PM Bhavesh Shrimali < bhavesh.shrimali@gmail.com> wrote:
After some more digging I am still unable to get the same solver setup using the same set of parameters. This time trying with tsfc for the form compiler (an old fork though: git+ https://github.com/blechta/tsfc.git@2018.1.0). Even with the same quadrature degree and the same form I get different SNES function norms in the very first solving step. Note that there are no convergence issues in the first step but eventually there are later.
Though with different PETSc versions should I expect similar convergence behavior ?
The first SNES residual is just the norm of F(x) where x is the input vector. So if the input vectors are the same, then the equations are not. This is not yet connected with PETSc.
Thanks,
Matt
There is very little difference in my formulation since most of it is in UFL syntax. I will keep looking ... In the meantime, here are the outputs --
Firedrake:
0 SNES Function norm 4.630756653490e+03 1 SNES Function norm 6.610067614593e+02 2 SNES Function norm 8.595689678029e+00 3 SNES Function norm 7.177260814570e-03 4 SNES Function norm 7.839549826428e-09 Nonlinear firedrake_0_ solve converged due to CONVERGED_FNORM_RELATIVE iterations 4
Dolfin:
Solving nonlinear variational problem. 0 SNES Function norm 4.048156097083e+03 1 SNES Function norm 7.701194025455e+01 2 SNES Function norm 9.032060723182e-02 3 SNES Function norm 2.188411552504e-08 PETSc SNES solver converged in 3 iterations with convergence reason CONVERGED_FNORM_RELATIVE.
and the codes (I will try to reduce these further as they still are far from minimal (~ 30 lines)) Dolfin (petsc4py version: 3.12.0):
import numpy as np from dolfin import *
parameters["form_compiler"]["cpp_optimize"] = True parameters["form_compiler"]["quadrature_degree"] = 5 parameters["form_compiler"]["representation"] = "tsfc" parameters["linear_algebra_backend"] = "PETSc"
snes_solver_parameters = {"nonlinear_solver": "snes", "snes_solver": {"linear_solver": "lu", "maximum_iterations": 50, "report": True, "line_search": 'basic', "error_on_nonconvergence": True} }
def Psi(u, p, Cv): I = Identity(3) F = I + grad(u) C = F.T * F Cv_inv = inv(Cv) I1 = tr(C) JF = det(F)
Ie1 = inner(C, Cv_inv) Je = JF / sqrt(det(Cv))
Jstar = (lam1 + p + sqrt((p + lam1)**2 + 4 * lam1 * (mu1 + mu2))) / (2 * lam1)
psi_Eq = ((3**(1 - alph1)) / (2 * alph1)) * mu1 * (I1**alph1 - 3**alph1) +\ ((3**(1 - alph2)) / (2 * alph2)) * mu2 * (I1**alph2 - 3**alph2) -\ (mu1 + mu2) * ln(Jstar) + 0.5 * lam1 * (Jstar - 1)**2
psi_Neq = ((3**(1 - a1)) / (2 * a1)) * m1 * (Ie1**a1 - 3**a1) +\ ((3**(1 - a2)) / (2 * a2)) * m2 * (Ie1**a2 - 3**a2) -\ (m1 + m2) * ln(Je)
W_hat = (psi_Eq + p * (JF - Jstar) + psi_Neq) * dx
return W_hat
mesh = UnitCubeMesh(3, 3, 3) W1 = VectorElement('CG', mesh.ufl_cell(), 2) W2 = FiniteElement('CG', mesh.ufl_cell(), 1) W3 = TensorElement("DG", mesh.ufl_cell(), 0)
V = FunctionSpace(mesh, W1 * W2) V_Cv = FunctionSpace(mesh, W3) du, dp = TrialFunctions(V) dw = TrialFunction(V) v, ph = TestFunctions(V) wh = TestFunction(V) w = Function(V) u, p = split(w) Cv = Function(V_Cv) Cv.interpolate(Constant(np.eye(3)))
left_fc = CompiledSubDomain( "near(x[0], 0.) && on_boundary" ) right_fc = CompiledSubDomain( "near(x[0], 1.) && on_boundary" ) bottom_fc = CompiledSubDomain( "near(x[1], 0.) && on_boundary" ) back_fc = CompiledSubDomain( "near(x[2], 0.) && on_boundary" )
strtch = Constant(0.0)
bc_right = DirichletBC(V.sub(0).sub(0), strtch, right_fc) bc_left = DirichletBC(V.sub(0).sub(0), Constant(0.0), left_fc) bc_bottom = DirichletBC(V.sub(0).sub(1), Constant(0.0), bottom_fc) bc_back = DirichletBC(V.sub(0).sub(2), Constant(0.0), back_fc)
bcs = [bc_right, bc_left, bc_bottom, bc_back] # Material Properties VHB4910 mu1 = Constant(13.54 * 1e3) mu2 = Constant(1.08 * 1e3) lam1 = Constant((mu1 + mu2) * 10**3) alph1 = Constant(1.) alph2 = Constant(-2.474) m1 = Constant(5.42 * 1e3) m2 = Constant(20.78 * 1e3) a1 = Constant(-10.) a2 = Constant(1.948) K1 = Constant(3507 * 1e3) K2 = Constant(1.0 * 1.e-12) bta1 = Constant(1.852) bta2 = Constant(0.26) eta0 = Constant(7014 * 1e3) etaInf = Constant(0.1 * 1e3)
S = derivative(Psi(u, p, Cv), w, wh) Jac = derivative(S, w, dw) problem = NonlinearVariationalProblem(S, w, bcs, J=Jac) solver = NonlinearVariationalSolver(problem) solver.parameters.update(snes_solver_parameters)
lamdot = 0.01 # Loading Rate Tfinal = 2 * (2 / lamdot) n_count = 101 dt = Tfinal / (n_count - 1)
timeVals = np.linspace(0, Tfinal, n_count)
stretchVals = np.hstack((lamdot * timeVals[:len(timeVals) // 2], lamdot * (-timeVals[len(timeVals) // 2:] + 2 * timeVals[len(timeVals) // 2]))) strtch.assign(stretchVals[1]) solver.solve() u, p = w.split()
Firedrake (petsc4py version: 3.15.1)
import numpy as np from firedrake import *
parameters["form_compiler"]["quadrature_degree"] = 5
snes_solver_parameters ={"snes_type": "newtonls", "ksp_type": "preonly", "pc_type": "lu", "snes_max_it": 50, "snes_linesearch_type": "basic", "snes_monitor": None, "snes_converged_reason": None}
def Psi(u, p, Cv): I = Identity(3) F = I + grad(u) C = F.T * F Cv_inv = inv(Cv) I1 = tr(C) JF = det(F)
Ie1 = inner(C, Cv_inv) Je = JF / sqrt(det(Cv))
Jstar = (lam1 + p + sqrt((p + lam1)**2 + 4 * lam1 * (mu1 + mu2))) / (2 * lam1)
psi_Eq = ((3**(1 - alph1)) / (2 * alph1)) * mu1 * (I1**alph1 - 3**alph1) +\ ((3**(1 - alph2)) / (2 * alph2)) * mu2 * (I1**alph2 - 3**alph2) -\ (mu1 + mu2) * ln(Jstar) + 0.5 * lam1 * (Jstar - 1)**2
psi_Neq = ((3**(1 - a1)) / (2 * a1)) * m1 * (Ie1**a1 - 3**a1) +\ ((3**(1 - a2)) / (2 * a2)) * m2 * (Ie1**a2 - 3**a2) -\ (m1 + m2) * ln(Je)
W_hat = (psi_Eq + p * (JF - Jstar) + psi_Neq) * dx
return W_hat
mesh = UnitCubeMesh(3, 3, 3) W1 = VectorElement('CG', mesh.ufl_cell(), 2) W2 = FiniteElement('CG', mesh.ufl_cell(), 1) W3 = TensorElement("DG", mesh.ufl_cell(), 0)
V = FunctionSpace(mesh, W1 * W2) V_Cv = FunctionSpace(mesh, W3) du, dp = TrialFunctions(V) dw = TrialFunction(V) v, ph = TestFunctions(V) wh = TestFunction(V) w = Function(V) u, p = split(w) Cv = Function(V_Cv) Cv.interpolate(Constant(np.eye(3)))
left_fc = 1 right_fc = 2 bottom_fc = 3 back_fc = 5
strtch = Constant(0.0)
bc_right = DirichletBC(V.sub(0).sub(0), strtch, right_fc) bc_left = DirichletBC(V.sub(0).sub(0), Constant(0.0), left_fc) bc_bottom = DirichletBC(V.sub(0).sub(1), Constant(0.0), bottom_fc) bc_back = DirichletBC(V.sub(0).sub(2), Constant(0.0), back_fc)
bcs = [bc_right, bc_left, bc_bottom, bc_back] # Material Properties VHB4910 mu1 = Constant(13.54 * 1e3) mu2 = Constant(1.08 * 1e3) lam1 = Constant((mu1 + mu2) * 10**3) alph1 = Constant(1.) alph2 = Constant(-2.474) m1 = Constant(5.42 * 1e3) m2 = Constant(20.78 * 1e3) a1 = Constant(-10.) a2 = Constant(1.948) K1 = Constant(3507 * 1e3) K2 = Constant(1.0 * 1.e-12) bta1 = Constant(1.852) bta2 = Constant(0.26) eta0 = Constant(7014 * 1e3) etaInf = Constant(0.1 * 1e3)
S = derivative(Psi(u, p, Cv), w, wh) Jac = derivative(S, w, dw) problem = NonlinearVariationalProblem(S, w, bcs, J=Jac) solver = NonlinearVariationalSolver(problem, solver_parameters=snes_solver_parameters) lamdot = 0.01 # Loading Rate Tfinal = 2 * (2 / lamdot) n_count = 101 dt = Tfinal / (n_count - 1)
timeVals = np.linspace(0, Tfinal, n_count)
stretchVals = np.hstack((lamdot * timeVals[:len(timeVals) // 2], lamdot * (-timeVals[len(timeVals) // 2:] + 2 * timeVals[len(timeVals) // 2]))) strtch.assign(stretchVals[1]) solver.solve() u, p = w.split()
Thanks for your help!
On Sun, 27 Jun 2021 at 13:21, Bhavesh Shrimali < bhavesh.shrimali@gmail.com> wrote:
Thanks Matt,
The only difference I see in the forms is that I am using `project` in Firedrake which I believe does a cell-wise projection for a DG-space whereas I use the `LocalSolver` in dolfin to do a cell-wise projection myself.
Even the SNES function norms that I am seeing aren't the same between FEniCS and Firedrake. Possible that I made some mistake in the formulation.
Thanks for pointing that out. I'll investigate this further today
Thanks, Bhavesh
On Sun, 27 Jun 2021 at 13:10, Matthew Knepley <knepley@gmail.com> wrote:
Something is wrong with the Newton system in the last step. You can see that the LU residual is steadily increasing. At the last step it is 10^-6. The iteration appears to be approaching a singular Jacobian. This seems like a formulation problem.
Thanks,
Matt
On Sun, Jun 27, 2021 at 11:24 AM Bhavesh Shrimali < bhavesh.shrimali@gmail.com> wrote:
I think I mistakenly clicked reply instead of reply all. Attached is my response:
--
Thanks Matt,
Thanks for the swift response. These are the logs with basic, backtracking, cp and nleqerr. If I don't specify 'snes_linesearch_type' it defaults to 'basic'. At least that's what I could gather from the log.
- `cp` seemingly diverges the fastest. Of course I haven't touched the other options. - `nleqerr` goes the furthest I believe, but not that far from `basic` or `bt`.
My concern was that since the forms are essentially the same, and I am passing "seemingly" the same parameters in FEniCS and Firedrake (both of which are using PETSc's SNES), what is it that is different which causes the solver to fail in one but not the other.
Edit: Just saw Ed's message. Indeed it defaults to `basic`
For definiteness this was the dictionary of solver parameters I am passing to NLVS in Firedrake
solver_params={ "snes_type": "newtonls", "ksp_type": "preonly", "pc_type": "lu", "snes_max_it": 100, # "snes_stol": 1.e-6, # "snes_rtol": 1.e-6, # "snes_atol": 1.e-6, "snes_view": None, # "ksp_atol": 1.e-10, # "ksp_rtol":1.e-8, "ksp_view": None, "snes_linesearch_type": "bt", "snes_converged_reason": None, "snes_monitor": None, "ksp_converged_reason": None, "ksp_monitor_true_residual": None }
On Sun, 27 Jun 2021 at 11:22, Bhavesh Shrimali < bhavesh.shrimali@gmail.com> wrote:
> Thanks Matt, > > > Thanks for the swift response. > These are the logs with basic, backtracking, cp and nleqerr. If I > don't specify 'snes_linesearch_type' it defaults to 'basic'. At least > that's what I could gather from the log. > > > - `cp` seemingly diverges the fastest. Of course I haven't > touched the other options. > - `nleqerr` goes the furthest I believe, but not that far from > `basic` or `bt`. > > My concern was that since the forms are essentially the same, and I > am passing "seemingly" the same parameters in FEniCS and Firedrake (both of > which are using PETSc's SNES), what is it that is different which causes > the solver to fail in one but not the other. > > Edit: Just saw Ed's message. Indeed it defaults to `basic` > > > For definiteness this was the dictionary of solver parameters I am > passing to NLVS in Firedrake > > solver_params={ > "snes_type": "newtonls", > "ksp_type": "preonly", > "pc_type": "lu", > "snes_max_it": 100, > # "snes_stol": 1.e-6, > # "snes_rtol": 1.e-6, > # "snes_atol": 1.e-6, > "snes_view": None, > # "ksp_atol": 1.e-10, > # "ksp_rtol":1.e-8, > "ksp_view": None, > "snes_linesearch_type": "bt", > "snes_converged_reason": None, > "snes_monitor": None, > "ksp_converged_reason": None, > "ksp_monitor_true_residual": None > } > > > On Sun, 27 Jun 2021 at 06:00, Matthew Knepley <knepley@gmail.com> > wrote: > >> Give >> >> -snes_view -snes_converged_reason -snes_monitor >> -ksp_converged_reason -ksp_monitor_true_residual >> >> and send all the output. That way I might be able to see what is >> happening when you do not converge. Also, >> line search 'basic' is not a good choice usually. I would either >> use the default, or 'cp', or maybe 'nleqerr' that >> Patrick added. >> >> Thanks, >> >> Matt >> >> On Sat, Jun 26, 2021 at 10:14 PM Bhavesh Shrimali < >> bhavesh.shrimali@gmail.com> wrote: >> >>> Hi everyone, >>> >>> I am trying to understand the usage of solver parameters that are >>> passed to the NonlinearVariationalSolver in firedrake. In order to do some >>> tests, I decided to reproduce some of my previous calculations done in >>> FEniCS for nonlinear viscoelasticity. In old dolfin, I had the following >>> for solver parameters to be passed to PETSc's SNES: >>> >>> snes_solver_parameters = {"nonlinear_solver": "snes", >>> "snes_solver": {"linear_solver": "lu", >>> "maximum_iterations": 50, >>> "report": True, >>> "line_search": 'basic', >>> "error_on_nonconvergence": False, >>> "relative_tolerance": 1e-7, >>> "absolute_tolerance": 1e-8}} >>> ... >>> ... >>> solver = NonlinearVariationalSolver(problem) >>> solver.parameters.update(snes_solver_parameters) >>> >>> Now as I understand, in Firedrake the options passed to the solver >>> are consistent with the corresponding names in PETSc (correct me if I am >>> wrong), and therefore I did: >>> >>> solver_params={ >>> "snes_type": "newtonls", >>> "ksp_type": "preonly", >>> "pc_type": "lu", >>> "snes_max_it": 100, >>> "snes_rtol": 1.e-7, >>> "snes_atol": 1.e-8, >>> "snes_view": None, >>> "ksp_atol": 1.e-10, # not needed? >>> "ksp_rtol":1.e-8, >>> "ksp_view": None, >>> "snes_linesearch_type": "basic", >>> "snes_no_convergence_test": 1 # this has no effect >>> } >>> solver = NonlinearVariationalSolver(problem, >>> solver_parameters=solver_params) >>> >>> However, for some reason, I am unable to get the solver to >>> converge when trying to solve the exact same problem as FEniCS. I >>> experimented with different quadrature degrees as well as tweaked >>> the solver parameters for SNES to best match what I am doing in dolfin. I >>> understand that the form compiler is different hence the precise code that >>> is generated may be different (I am not an expert on code generation) even >>> though the python scripts may look similar due to UFL syntax. >>> >>> In fact the weird part is that the solver seems to converge fine >>> till a certain point in my input (loading). You can see this in the plot >>> below where I plot the final stress-stretch relation obtained at a point in >>> my domain from dolfin and firedrake.. >>> >>> [image: firedrake_vs_fenics.png] >>> >>> >>> I have marked the abscissa beyond which the solver fails to >>> converge in Firedrake. The complete codes are here: >>> https://github.com/bhaveshshrimali/comparingFiredrakeFEniCS >>> >>> I tried setting the option `snes_no_convergence_test` to 1/True >>> but that had no effect as the SNES solver still diverges with ` >>> DIVERGED_DTOL`. In the end I put the call to solve inside a >>> try/catch block to at least be able to see the solution till that point. >>> And it surprisingly is accurate till that point. >>> >>> Would anyone have any clue of what could be going wrong? I >>> wouldn't expect that you go through the entire code(s) as that could be >>> fairly long, but any pointers would be helpful :) >>> >>> I am using firedrake via Docker on windows. >>> >>> Thanks >>> >>> -- >>> >>> _______________________________________________ >>> firedrake mailing list >>> firedrake@imperial.ac.uk >>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >>> >> >> >> -- >> What most experimenters take for granted before they begin their >> experiments is infinitely more interesting than any results to which their >> experiments lead. >> -- Norbert Wiener >> >> https://www.cse.buffalo.edu/~knepley/ >> <http://www.cse.buffalo.edu/~knepley/> >> > > > -- > > *Bhavesh Shrimali* > *Ph.D. Student, **Civil and Environmental Engineering,* > > *University of Illinois Urbana-Champaign* > > > *Champaign, IL, United States* > > *Email : bshrima2@illinois.edu <bshrima2@illinois.edu>* > *Web : **https://bhaveshshrimali.github.io/ > <https://bhaveshshrimali.github.io/>* >
--
*Bhavesh Shrimali* *Ph.D. Student, **Civil and Environmental Engineering,*
*University of Illinois Urbana-Champaign*
*Champaign, IL, United States*
*Email : bshrima2@illinois.edu <bshrima2@illinois.edu>* *Web : **https://bhaveshshrimali.github.io/ <https://bhaveshshrimali.github.io/>* _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener
https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
--
*Bhavesh Shrimali* *Ph.D. Student, **Civil and Environmental Engineering,*
*University of Illinois Urbana-Champaign*
*Champaign, IL, United States*
*Email : bshrima2@illinois.edu <bshrima2@illinois.edu>* *Web : **https://bhaveshshrimali.github.io/ <https://bhaveshshrimali.github.io/>*
--
*Bhavesh Shrimali* *Ph.D. Student, **Civil and Environmental Engineering,*
*University of Illinois Urbana-Champaign*
*Champaign, IL, United States*
*Email : bshrima2@illinois.edu <bshrima2@illinois.edu>* *Web : **https://bhaveshshrimali.github.io/ <https://bhaveshshrimali.github.io/>*
-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener
https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
--
*Bhavesh Shrimali* *Ph.D. Student, **Civil and Environmental Engineering,*
*University of Illinois Urbana-Champaign*
*Champaign, IL, United States*
*Email : bshrima2@illinois.edu <bshrima2@illinois.edu>* *Web : **https://bhaveshshrimali.github.io/ <https://bhaveshshrimali.github.io/>*
-- *Bhavesh Shrimali* *Ph.D. Student, **Civil and Environmental Engineering,* *University of Illinois Urbana-Champaign* *Champaign, IL, United States* *Email : bshrima2@illinois.edu <bshrima2@illinois.edu>* *Web : **https://bhaveshshrimali.github.io/ <https://bhaveshshrimali.github.io/>*
participants (2)
-
Bhavesh Shrimali
-
Matthew Knepley