Create simple MOOSE nonlinear solver

[WilkAndy]

Description of the enhancement or error report

Create a utility that can solve small nonlinear systems.

Rationale for the enhancement or information for reproducing the error

In plasticity in tensor mechanics we have to solve a small nonlinear system at each quadpoint in order to compute the stress. The system has between 1 and 20 dofs, and is typically around 4 dofs. Currently i have my own solver, that employs the standard Newton-Raphson + cubic-line-search approach. Possibly (probably?) others have their own equivalents buried in the code somewhere.

It'd be great if there was a central solver that could do this for us.

When i implemented my own, people advised me from using PETSc, but now we've got @fdkong, perhaps PETSc is the thing to use.

I imagine the interface might be something like this:

void nonlinearSolver(void *residual_calculator, void *jacobian_calculator, const std::vector<Real> & initial_guess, std::vector<Real> & solution, bool use_linesearch)

It should be efficient as this thing will get used (by me) at every quadpoint for every MOOSE Residual calculation.

Identified impact

(i.e. Internal object changes, limited interface changes, public API change, or a list of specific applications impacted)

Uniformize MOOSE, improve maintainability, improve efficiency.

I don't have time myself to do this, but am hoping it'll interest someone!

[fdkong]

@WilkAndy, could you show me an example? It is interesting, but I want to understand how you use this.

[YaqiWang]

We did this. It will be nice to have a uniform utility function for doing this. But PETSc might be too heavy.

[fdkong]

We did this. It will be nice to have a uniform utility function for doing this. But PETSc might be too heavy.

Yes, you are right. But if we could sweep over the quadrature points using the same SNES, it may be not too bad.

[WilkAndy]

@YaqiWang has done this?! Where is it (in MOOSE)? I guess this just reinforces my request - it's dumb if we all have our own newton-raphson-plus-line-search routines inside MOOSE!

@fdkong , see tensor_mechanics/src/materials/PQPlasticModel. The core methods are nrStep and lineSearch, and the iterative loop lies between line 284 and line 306.

a

[YaqiWang]

It is in Rattlesnake, our internal radiation transport application.

[WilkAndy]

@YaqiWang , if you can bring it out of Rattlesnake and into framework/utils, can you have a quick look at the code i mentioned above and see if your Rattlesnake version could accommodate what i want to do? Then perhaps it's an easy job to move it to framework?

[YaqiWang]

Just looked. Mine is much simpler, only one variable, from an adiabatic temperature model, where heat capacity depends on temperature as a Function.:

rho C_p(T) dT/dt = power.

Probably whatever you do will be directly useable by me. So feel free to contribute. I will patch if later I find out it is not enough to me.

[WilkAndy]

OK, good to know, thanks @YaqiWang ,

a

[jessecarterMOOSE]

FYI I'm all for this. We could definitely make use of it.

[friedmud]

As far as I know - there is no third party package that does a dense Newton solve that is optimized for tiny matrices. The one that I can see that's the closest is Ceres: http://ceres-solver.org/index.html

@WilkAndy it looks like you're using LAPACK directly for your linear solver.

My suggestion would be to utilize Eigen matrices and vectors and build up a nice little object-oriented nonlinear solver. This could actually be done as a separate project from MOOSE and then we could import it into either libMesh/contrib or moose/contrib so it's available for all apps.

I think the wider computational science community would appreciate an open-source, hardcore, tiny nonlinear solver library built using Eigen. Could generate quite a bit of interest and contributions.

(If you're not familiar with Eigen: http://eigen.tuxfamily.org . We already have it available in libMesh/contrib and a lot of the libMesh DenseMatrix and DenseVector routines can fall back to using it. I've used it pretty extensively and find that it is incredibly fast.)

[WilkAndy]

Thanks for all those remarks, @friedmud . That Eigen looks nice. It's truly amazing that a C++ NR+LineSearch doesn't yet exist!

@fdkong , are you interested in doing this?

[friedmud]

@WilkAndy I think the reason it doesn't is that most people want to "code their own" because at this size you can really squeeze out large computational gains by tying your solver a little closer to your problem.

For instance: if you know your matrix is always symmetric you can optimize storage and matvecs. If you know that the Jacobian isn't changing (or is varying slowly) you can lag it's re-evaluation. If you know the Jacobian is pieced together by several non-changing matrix plus a changing one... then you can take advantage of that... etc. There are so many tiny little optimizations (many of which you are actually using yourself in the TensorMechanics stuff) at this level that trying to write a fast general solver may not be possible.

That said: the "MOOSE way" is generally to ignore the loss of optimality in order to create a better dev environment... so I'm actually for making an attempt at this!

[WilkAndy]

Yea, if someone wants to optimise the core functionality then they're welcome to! For a start, I'd just like something that works.

[fdkong]

I will try a PETSc snes, and if it is much slower than the one from @WilkAndy, we possibly consider Eigen (the name is misleading). For a small system, we could choose to have an analytical inverse in PETSc, and it costs not much.

[WilkAndy]

Sounds great, thanks @fdkong ,

a