SNES solver: Added a PID controller to update the timestep

[malamast]

The user sets target_its a desired number of nonlinear iterations and the timestep is updated based on the ratio of that number and the actual number on nonlinear iteration of each Newton step. The exponent kI can be tuned to the problem needs. A smaller exponent results in smaller changes of dt.

It has been tested with the 1D-threshold example and it gave a small speed up over the old method. There is good agreement between the different solvers for most of the variables. It seems that velocities are more sensitive to the accuracy of the solver. See also the attached figures where I compare the solution at the last timestep. Some figures show the time evolution of the variables at the middle of the 1D domain.

I have also added an option to let the user set the maximum number of function evaluations maxf . One might need to increase maxf, if large values of output_step are set.

How to use:
nout = 50
timestep = 95788 * 0.05 # 95788 = 1 milisecond

[solver]
type = beuler # Backward Euler steady-state solver
snes_type = newtonls # Newton iterations
ksp_type = gmres # GMRES method
pc_type = hypre # Preconditioner type
pc_hypre_type = euclid # Hypre preconditioner type
max_nonlinear_iterations = 16 # Max Newton iterations per timestep
lag_jacobian = 6 # How long to wait before Jacobian recalculation per Newton iteration

atol = 1e-12 # Absolute tolerance
rtol = 1e-6 # Relative tolerance
stol = 1e-10 # Convergence tolerance
maxf = 20000 # Maximum number of function evaluations
maxl = 120 # default: 20
#diagnose = true
#diagnose_failures = true

pidController = true
target_its = 5
kP = 0.7
kI = 0.3
kD = 0.2

[github-actions]

clang-tidy made some suggestions

[ZedThree]

I've deleted the clang-tidy comments as the clang-format commit rearranged things and confused it, and put them back on the right places.

It would be nice to get a test on this. We have an integrated test for all the solvers, but it's not currently set up to handle different options on the same solver. I think I have an idea on how to handle that -- but it might be best to do that in a separate PR?

[ZedThree]

We can drop the static_cast here

Suggested change

	.withDefault(static_cast<int>(7))),
	.withDefault(7)),

[ZedThree]

These can all be const:

Suggested change

	BoutReal facP = std::pow(double(target_its) / double(nl_its), kP);
	BoutReal facI = std::pow(double(nl_its_prev) / double(nl_its), kI);
	BoutReal facD = std::pow(double(nl_its_prev) * double(nl_its_prev) / double(nl_its)
	/ double(nl_its_prev2),
	kD);
	const BoutReal facP = std::pow(double(target_its) / double(nl_its), kP);
	const BoutReal facI = std::pow(double(nl_its_prev) / double(nl_its), kI);
	const BoutReal facD = std::pow(double(nl_its_prev) * double(nl_its_prev) / double(nl_its)
	/ double(nl_its_prev2),
	kD);

[ZedThree]

We can use std::clamp from <algorithm> to make this const:

Suggested change

	// clamp groth factor to avoid huge changes
	BoutReal fac = facP * facI * facD;
	if (fac < 0.2)
	fac = 0.2;
	else if (fac > 5.0)
	fac = 5.0;
	// clamp growth factor to avoid huge changes
	const BoutReal fac = std::clamp(facP * facI * facD, 0.2, 5.0);

[ZedThree]

We can use std::min to clamp to a max value and make this const:

Suggested change

	BoutReal dt_new = timestep * fac;
	if (dt_new > max_timestep) {
	dt_new = max_timestep;
	}
	const BoutReal dt_new = std::min(timestep * fac, max_timestep);

[ZedThree]

Suggested change

	nl_its_prev = static_cast<int>(nl_its);
	nl_its_prev = nl_its;

[ZedThree]

This can be:

Suggested change

	if (timestep > max_timestep) {
	timestep = max_timestep;
	}
	timestep = std::min(timestep, max_timestep);

[mikekryjak]

@malamast you need to run the simulations for longer - 50 timesteps is not enough to reach full steady state. It was chosen in the example to keep the solution time reasonable for both CVODE and SNES simultaneously. Maybe 200 timesteps would be enough to get all of the oscillations out for sure, and should hopefully result in all three converging to the same profiles.

The main question is how is the time evolution affected. It would be good to have plots of upstream density and target temperature as a function of time if you can. Your results so far definitely show that there are some differences - this is fine as we want a fast steady state solution, but it would be good to quantify it a bit.

Also, what was your motivation for implementing the maxf setting? Is it to prevent the solver from getting stuck in some new way?

P.S. it makes sense that velocity is the most different. The plasma sloshes around left and right for a while because this example doesn't have viscosity enabled.