Open GiudGiud opened 3 years ago
At a high level, I thought PicardSolve just provided a mechanism for performing an arbitrary number of "inner" iterations that terminate based on some custom criterion. It sounds like what you want here is maybe a way to customize the contents of PicardSolve::solveStep. Really this is getting back to design discussion we've have a number of times over the years about a more flexible, extendable executioner system - which you've referenced a bit in the tagged issues here. Regardless - I'd prefer any additions/mods related to this either be:
Rather than a new layer - I'd prefer to see a way to swap the picard layer in/out for these alternate inner iteration loops.
I want to modify solve as well, as the relaxation of the subapps (for picard, whichever other method for the others) is done there.
Making the PicardSolve general enough might make the code a bit hard to read imo. Picard is pretty simple, just one step and iterate. Some methods use half steps, others use a perturbation step, maybe in the future we'll have methods that work on the derivatives. Might be a lot of code for one file
What would you be thinking though if not a single class below SolveObject that each solve method is inheriting from, have the executioner hold a solveObject and cast it to the right type with a switch statement? I'm not sure how to do this without using a base class holding the shared routines.
The current executioner system has been flexible and extendable. Creating separate SolveObject
s will be the solution. Of course we can still use c++ inheritance to reduce potential code duplication. Note that there are other solve objects derived from SolveObject
. For example I added FixedPoint
in test. Griffin has another one. TimeStepper should be made as one of the solve objects. When I created SolveObject
, my main goal was to clean up Picard iteration so I did not finish it. Another thing I did not finish is to clean up TransientMultiApp
.
Having other options for solving a fully-coupled multiphysics system is great.
Yeah - I think I'm with Yaqi on this - those other solve types should each be their own "SolveObject" subclasses like PicardSolve (i.e. they should be sibling classes of PicardSolve).
I agree they should be sibling classes with PicardSolve. My question is now mostly, is there enough shared code that adding a layer, IterativeSolve, makes sense. (or if there is enough shared code, can I just add a layer?) Semantically it certainly does. A PicardSolve and a SecantSolve are much closer on paper than FEProblemSolve
I will start coding and believe that we will soon realize if creating another layer like SolveObject
=>MultiAppIteration
=>PicardSolve
or something else is needed.
@YaqiWang I was thinking of doing that
Yeah - I'd say start without a new shared layer. If as you go along, you see lots of duplicate code, we can revisit. If it's only a bit of duplicate code, I'd say its not worth adding the extra layer.
Reason
Currently, multiphysics couplings can be handled in a single application (full coupling) or using multiapps (loose and tight coupling). Tight coupling is achieved by iterating each physics and exchanging information until convergence. This is done using a Picard iteration process.
There are other options, namely Newton (or Householder)-like methods, bisection (or false position)-like methods. The former have a better asymptotic convergence rate (super-linear, quadratic for most) and the latter can be guaranteed to converge. I'm sure there are other more modern root finding methods that can apply as well, using stochastic processes or machine learning for example.
We should add these methods as an option to user, so they can converge the coupling faster. This requires some infrastructure change, as the Picard method is hardwired in TimeStepper and Executioner, Steady, Transient and Eigenvalue.
Design
I'm thinking a class between SolveObject and all the iteration algorithms (PicardSolve, SecantSolve etc) that would host as much of the shared parameters and solve routines as possible. It would also be called from the timestepper and initialized by the executioner.
Impact
Access to quadratic convergence for multiapp tight coupling Between 4 and 8x faster convergence for the 1D-2D primary loop TH coupling