A modified vector_fe_ex4 failed to converge

[qianhanhua]

Hi, dear maintainers,

I had a discussion [https://github.com//discussions/3243] about using libmesh to reproduce results of Mie scattering (Plane electromagnetic wave scattered by dielectric sphere). I was modifying vector_fe_ex1 for such model. With your suggestions I was able to get some meaningful results, but there were always big differences around the surface of the sphere, as compared with the analytical solution.

So I decided to try out with the curl curl example since the governing equation is very similar. I didn't change the main function much but mainly the curl_curl_system.c file, as attached here.
scatter_system.c.txt
However, after compilation and executing, the solver failed to converge:

Assembling the System
Nonlinear Residual: 2.79877e+08
Linear solve starting, tolerance 1e-12
Linear solve finished, step 5000, residual 0.0542179
Trying full Newton step
  Current Residual: 2.98125e+08
  Shrinking Newton step to 0.5
  Current Residual: 2.84622e+08
  Shrinking Newton step to 0.25
  Current Residual: 2.81108e+08
  Shrinking Newton step to 0.125
  Current Residual: 2.80204e+08
  Shrinking Newton step to 0.0625
  Current Residual: 2.79968e+08
  Shrinking Newton step to 0.03125
  Current Residual: 2.79904e+08
  Shrinking Newton step to 0.015625
  Current Residual: 2.79886e+08
  Shrinking Newton step to 0.0078125
  Current Residual: 2.7988e+08
  Shrinking Newton step to 0.00390625
  Current Residual: 2.79878e+08
  Shrinking Newton step to 0.00195312
  Current Residual: 2.79877e+08
  Shrinking Newton step to 0.000976562
  Current Residual: 2.79877e+08
  Shrinking Newton step to 0.000488281
  Current Residual: 2.79877e+08
  Shrinking Newton step to 0.000244141
  Current Residual: 2.79877e+08
  Shrinking Newton step to 0.00012207
  Current Residual: 2.79877e+08
  Shrinking Newton step to 6.10352e-05
  Current Residual: 2.79877e+08
  Shrinking Newton step to 3.05176e-05
  Current Residual: 2.79877e+08
  Shrinking Newton step to 1.52588e-05
  Current Residual: 2.79877e+08
Inexact Newton step FAILED at step 0

May I ask if there is any way to improve so the solver could produce a converged result? Thank you!

[roystgnr]

No time to check out your specific modifications right now, but in general the second step after "I tried a new formulation but my iterative solver failed to converge" (looks like 5000 steps didn't even give you a descent direction for the first quasi-Newton step here!) is "try a small enough problem with it that a direct solver is usable". If -pc_type lu (or more sophisticated, depending on your PETSc config) doesn't work either then the third step is generally to find the problem with the formulation, if it does work then the third step is generally to look for a better preconditioner.

[jwpeterson]

So should I look for a good preconditioner in Petsc's pc_type database after this?

No, I think you should first try a direct solver. If that has trouble solving, then it points to a bug of some kind in your implementation, so I would then focus on debugging with print statements, small test problems, etc. I would also confirm whether the problem is properly constrained with boundary conditions... in the previous discussion you mentioned using penalty BCs, presumably you are doing that here as well? When the residual stays extremely high like it is in your case, it can sometimes indicate an issue with constraints.

Finally, just FYI if I understand correctly, you are solving a linear problem, so while it's fine that you are using a Newton solver, you should not need to take more than one "nonlinear" iteration to solve the problem.

[qianhanhua]

Hi, @jwpeterson , I am not exactly sure how to switch to direct solver since -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type superlu doesn't work here. I think that's because solver was specified in the main function:

  // This example only implements the steady-state problem
  system.time_solver = libmesh_make_unique<SteadySolver>(system);

  // Initialize the system
  equation_systems.init();

  // And the nonlinear solver options
  DiffSolver & solver = *(system.time_solver->diff_solver().get());
  solver.quiet = infile("solver_quiet", true);
  solver.verbose = !solver.quiet;
  solver.max_nonlinear_iterations = infile("max_nonlinear_iterations", 50);
  solver.relative_step_tolerance = infile("relative_step_tolerance", 1.e-3);
  solver.relative_residual_tolerance = infile("relative_residual_tolerance", 1.0e-13);
  solver.absolute_residual_tolerance = infile("absolute_residual_tolerance", 1e-8);
  // And the linear solver options
  solver.max_linear_iterations = infile("max_linear_iterations", 50000);
  solver.initial_linear_tolerance = infile("initial_linear_tolerance", 1.e-10);

For BCs, I am using penalty BCs as in the original curl_curl_system.c, except that I provided the analytical solution of Mie theory.
p.s., I also tried with BCs that forces the boundary unknowns to be zero using penalty method. Same convergence issue appeared.

[jwpeterson]

Just use -ksp_type preonly -pc_type lu and drop the part about superlu. That only works if you enabled superlu in your PETSc build. The command line solver options should then override any solver options that are set in the code.

[qianhanhua]

Tried that but it still went to the Newton steps (and failed) though the residual is 0 after the Linear solve step.

[qianhanhua]

Just use -ksp_type preonly -pc_type lu and drop the part about superlu. That only works if you enabled superlu in your PETSc build. The command line solver options should then override any solver options that are set in the code.

My "doesn't work" meant that the command line options was not overriding code specified solver. With the superlu option the program starts fine but failed for same reason and output.

[vikramvgarg]

Hi qianhanhua, What happens if you change the example in a more incremental fashion ? For example, change only the boundary conditions to the one that you want and see if the solver converges. That way we can better understand if the constraints are an issue like @jwpeterson mentioned. Vikram Garg

…

On Thu, May 12, 2022 at 9:24 AM qianhanhua ***@***.***> wrote: Just use -ksp_type preonly -pc_type lu and drop the part about superlu. That only works if you enabled superlu in your PETSc build. The command line solver options should then override any solver options that are set in the code. My "doesn't work" meant that the command line options was not overriding code specified solver. with the superlu the program starts fine but failed for same reason and output. — Reply to this email directly, view it on GitHub <#3286 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABGNEEKWDVQVZI2MO2LYWCTVJUWFPANCNFSM5VYOGD6A> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[qianhanhua]

Hi, all, after changing BCs to simply zero at the boundary (model itself also added a layer of PML), the solver was able to give a result using direct solver. @jwpeterson Please ignore my previous comments on -ksp_type preonly -pc_type lu not overriding code specified solver. I was wrong due to lack of knowledge. The direct solver didn't work because of previous BCs.

@vikramvgarg I think we can confirm that it was due to the BCs. I was trying to truncate the problem with the analytical solutions at some virtual surface that encloses the problem domain. Seems more work needs to be done to make this idea working.

Now the modified curl curl example give the same results as the modified vector Poisson example (vector_fe_ex1). They both have large error near the sphere's outer surface. So I am facing the same problem to improve the model's accuracy then.

Thanks to all for the help on fixing the solver problem!

[vikramvgarg]

You are welcome.