Minor question about automatic_scaling

[r-williams-2]

Hello,

Just got a minor question about using automatic_scaling=true in the executioner block.

I have been using this quite successfully but I wanted to test individual physics to be sure it is working as intended. As such, I have set up a simple 1D diffusion problem with no-flux boundaries (using the grand potential phase-field method). This has one solution variable and a couple of auxiliary variables.

Without the use of automatic_scaling, the initial residual is massive and reduces easily with Newton iterations. However, when I use automatic_scaling (as used when solving the coupled problems I am working on) the initial nonlinear residual is very small (in the order of e-18) and is assumed to have already converged according to the nl_rel_tol = nl_abs_tol = 1.0e-8.

Is there any reason why automatic_scaling would not work on a problem with a single variable? Is it something to do with the auxiliary variables? I am only concerned that when used in the coupled problem with automatic_scaling, the residual for this variable isn't reducing sufficiently to solve the problem correctly. Is the solution to set a lower absolute tolerance on the nonlinear residuals?

If you require/desire the input file, please let me know and I can provide this.

As always, any help is greatly appreciated!

Kind regards,
Ross.

[aeslaughter]

@lindsayad

[lindsayad]

In general we try to discourage use of nl_abs_tol. I am guessing based on your description that you have Dirichlet boundary conditions? If so are you using preset or not? I would highly recommend using preset Dirichlet conditions if you have Dirichlet conditions.

Anyway one way to ensure that your problem is truly converged is to use ReferenceResidualProblem. Please see this page for more documentation about it.

My first suggestion would be to just try ReferenceResidualProblem and see what happens. If you are unable to converge with nl_rel_tol = 1e-8 (remember to remove nl_abs_tol) perhaps because the initial residual is still too "low", I would try adding resid_vs_jac_scaling_param = 0.5 to your Executioner block. Doing that will ensure that both residual and Jacobian information is considered when determining the scaling factor. See the NonlinearSystemBase for more information. I would also recommend adding verbose = true to your Executioner block just to see what the scaling factor being computed is, both for the default automatic scaling setting and with the added resid_vs_jac_scaling_param = 0.5 parameter.

If you still have issues, please feel free to attach your input file here.

[r-williams-2]

Noted regarding nl_abs_tol. I'd picked it up from an example problem and kept using it. I do not have any Dirichlet BCs. In the simplified problem, I just have grad(mu)=0 on both boundaries and letting a linear chemical potential profile relax. This will be important towards the end of this post as it is similar to the case with no BCs applied in the mechanical case.

ReferenceResidualProblem is also very useful to know about for when I am using other physics in a coupled model. I'm glad I asked about this now. Thanks for that.

So, having tried removing the nl_abs_tol, it seems to work. I think the main issue was the use of nl_abs_tol=1e-8 as the 1D problem only has a single solution variable (so it wasn't other residuals dominating the problem) and there are no fluxes applied on the boundaries (so using the integration of the fluxes instead of the initial residual is probably not useful). When scaled, the initial residual was lower than 1e-8. I think when unscaled, it just wasn't lower than the absolute tolerance and was converging because of the nl_rel_tol. The standard scaling factor was 2.04503e-26 (so you can imagine how even an initial residual of 1e+9 could be small enough to have "converged" when only using automatic_scaling=true and nl_abs_tol=1e-8). When I used resid_vs_jac_scaling_param = 0.5, it's only scaled by 9.99547e-18. While this doesn't solve the problem if nl_abs_tol=1e-8, I think this is a better solution moving forward since, as you say, it takes into account the size of the residual and the Jacobian entries.

Now, coming back to the BCs of grad(mu)=0 and that the problem relaxes to the equilibrium state (so the driving force for evolution vanishes and, consequently, the magnitude of the residual also reduces). I tried initiaising a problem close to equilibrium and did indeed find difficulties in dropping the residual according to nl_rel_tol=1e-8 (managed to get to about 1.1e-8 but increasing the nl_rel_tol only pushes this problem further ahead in simulated time). As per the ReferenceResidualProblem documentation, the only solution I can see for this is to apply a suitably small nl_abs_tol. Is there a good way to identify a physically-reasonable value for this other than separating the physics out like I have done? In this 1D case, I have found that a nl_abs_tol=1e-18 (much smaller than before) works well.

As promised, if you want to play about with it to see what I mean, I have attached the input file. If there is a better way to get this to you, just let me know.

Thanks again!

Diffusion_1D_simplified.zip

[lindsayad]

You have a very interesting input. You have single physics (obviously as you said) and no Dirichlet conditions...so actually your condition number is totally independent of your scaling factor. You can test your condition number by passing on the command line -pc_type svd -pc_svd_monitor. Note that you will not want to do this if your number of dofs is much greater than 1000...it takes too much CPU and memory. So yea regardless of scaling factor your condition number is about 4...which is really good. 1 is perfect. 1e16 is when you start running out of double precision.

And even more unique are the conditions you outlined...you just have two terms, the time derivative and diffusion. As t -> infinity both of those terms will go to zero, so with ReferenceResidualProblem given enough time you will not be able to improve 8 orders of magnitude relative to zero. So you have arrived at one of the cases where nl_abs_tol is truly appropriate.

I have attached an updated input file that uses ReferenceResidualProblem that works well without an absolute tolerance for about 18 time steps, and then it starts running out of precision.
Diffusion_1D_simplified.txt

[r-williams-2]

Fantastic. Thanks for that. I think I understand how everything operates.

The only other questions that remain are:

• Is there a better way of identifying a suitable nl_abs_tol than isolating your physics and using a problem near equilibrium to identify a reasonable value?
• Can I use something similar to ReferenceResidualProblem to set separate nl_abs_tol for the residual for each variable in case they are of different magnitude or should the individual scaling of each residual handle this issue?

[lindsayad]

There is another way to avoid using (and guessing) an nl_abs_tol. We have the Executioner/steady_state_detection boolean parameter. If you set that to true then if the solution changes by less than the tolerance ss_check_tol (default is 1e-8), then your simulation is considered to be at a steady-state and your transient solve stops. This is very handy. However, you may still have to set a nl_abs_tol (or loosen the ss_check_tol) because you may run out of precision for nl_rel_tol as you're running up to that steady state. There is no tried and true way I think for determining nl_abs_tol...you kind of have to run your problem and see.

Unfortunately, it is not currently possible to set an nl_abs_tol for individual variables. I think we may have requests for that feature in the past, but no one has gotten around to implementing it. Feel free to search through our issues, and if we don't have an issue for that yet, please feel free to create one.

[r-williams-2]

Ahh. Yes. I had forgotten about this. I should probably start using that. I will be sure to keep an eye out for this specific scenario.

I will have a look for that now.

Once again, thanks for all your help on this matter.