erdc / proteus

A computational methods and simulation toolkit
http://proteustoolkit.org
MIT License
88 stars 56 forks source link

Heavy duplication of code between Newton and NewtonNS classes in NonlinearSolvers.py #119

Open ahmadia opened 10 years ago

ahmadia commented 10 years ago

It looks to me like the only differences between the two are that NewtonNS contains a specialized norm. Each of these classes is about 300 lines of code, and I see at least one fix that's been committed to NewtonNS that has not made it back in to Newton. It looks like we could generalize the Newton class slightly and sharply reduce the number of lines in NewtonNS (and probably Picard/SSPRK as well, but I haven't looked at them carefully yet).

Any objections to a PR refactoring some of this functionality?

cekees commented 10 years ago

I think that would be good. The original Newton class was written to allow passing in other norms and just defaults to the l_2 norm. Unfortunately the default convergence test it provides isn't general enough to support splitting tolerances for different components of the PDE unless you were to go to the trouble of writing a special weighted norm.

I would really like to do a consistent refactoring of the way we handle tolerances and feed them to temporal/spatial/nonlinear/linear error norms. There's a pretty decent approach using weighted norms that is used in many adaptive time integration codes. The basic idea is to allow vector absolute and relative tolerances of the same length as the measured quantity. For example you can write a weighted l_2 norm as

| x |_{w2} = \sqrt(\sum (w_j x_j)^2)

where the weight is

w_j = 1/(atol_j + y_j*rtol_j)

and y is a vector of reference values (e.g. the solution at the last time step). In that norm | x |_{w2} < 1 implies that in the mean

x_j < atol_j + s_j rtol_j

You could basically replace NewtonNS with Newton using a weighted norm like that and an array of tolerances that are constant in two blocks (one block for the continuity equation and one for the momentum equations). You could also add scaling with respect to the element diameter so the norm would behave nicely with respect to mesh refinement (and more closely approximate the L_2 norm).

ahmadia commented 10 years ago

As noted in https://github.com/erdc-cm/proteus/commit/62521f83526ff91e26e08062f4a8fe2bf854f378#commitcomment-7909712, there appears to be a bug in the current implementation of NewtonNS's solve method, so this is another reason for refactoring the method out.

adimako commented 5 years ago

@cekees Not sure if still is the case. I have seen that Aron has done some refactoring after this issue was opened. If we do not plan to update these classess (apart from maintenance e.g. moving to python 3) and we are happy that these work because they pass tests, maybe not an issue any more