Strong coupling of a non-linear multiphysic problem: failure with Newton Raphson method

Question

I am trying to solve a multiphysic problem using finite elements and a Newton Raphson solution scheme. I have two non-linear subsystems that are coupled bi-directionally.

The first subsystem includes an equation that I solve for $p_f$:

$\nabla u - \nabla a^3 \nabla p_f = \nabla a^3 B$

This equation becomes very non-linear, due to the parameter $a$ with an exponent of 3! $B$ are body forces and $u$ is the solution of the second subsystem that I am trying to solve. This second subsystem includes the Stokes equation that depends on the gradient of $p_f$. Therefore, both system are coupled bi-directionally.

I am using a Newton-Raphson solution scheme for each subsystem. First, I solve for $p_f$. After convergence occurred, I pass this solution to the second subsystem, which I solve using a Newton Raphson scheme, too. Then, the solution $u$ is passed to the first subsystem ... until the solutions of both subsystems do not change anymore. See the flow chart below.

This procedure works very well so far. However, if the values of the parameter $a$ become very low (thus the gradient of $a$ becomes very, very low), there is no convergence anymore! The reason seems to be that the coupling between both subsystem becomes very strong. A small change in one subsystem creates a significant changes in the other one (especially where the gradient of $a$ is very low). I tried a variety of good initial conditions, however, since the first subsystem reacts very sensitiv, convergence has been impossible so far.

In the figure below, you can see how the solution of this approach ("Segregated Approach"-blue line) is different from the "true" solution ("Fully Coupled"-red line). The "true solution" is the result of a benchmark, which I calculated by solving a single system of equations that includes the equations of both subsystems.

The greatest differences between both solutions are located where the gradient of $a$ becomes very small. See the figure below (in this figure the values of $p_f$ and $a$ are normalized to range between 0 and 1):

So: Is a Newton-Raphson solution scheme the wrong method for solving a strongly coupled system, if there are low gradients in a subsystem? Could my method be fixed to solve the problem that I described? If yes, how? Are there alternative solution schemes for this problem?

Answer 1

The issues you're running into now are not a failing of Newton-Raphson, but a question of coupling. You're doing iterated sequential coupling -- solving each equation sequentially and then iterating until (hopeful) convergence.

No solver choice in place of NR is going to fix this lack of convergence, as long as you are doing iterated sequential coupling.

Instead, you want to consider different coupling schemes that use more information about the coupling terms. What @paul suggests is effectively equivalent to doing a quasi-Newton method using an approximation of the Jacobian of the coupled system that is block diagonal. Note that this is NOT Newton-Raphson any more, because you do not have the true Jacobian of the coupled system -- instead you form an approximation of the Jacobian using fewer terms than the true Jacobian and use this instead. This is a good starting point because it requires no additional knowledge of the physics -- you can build this entirely from pieces you already have. Norms for the coupled system become a combination of the norms for the two systems you have. Residuals for the coupled system become the vector of both systems' residuals. And the approximate inverse is block diagonal.

It is possible, but unlikely, that this will help much. The issue is that the approximate (block diagonal) inverse will be too far from the true inverse. One potential option is to try Jacobian-Free Newton Krylov to get a better inverse. This is likely to not work all that much better until you start to add off-diagonal blocks to the approximate inverse.

From a practical standpoint -- your best bet is to bite the bullet and find a way to start to include at least some information about how the equations couple in your approximate inverse, i.e. including off-diagonal blocks in your approximate Jacobian for use in a fully implicit coupling scheme.

Answer 2

Instead of running Newton-Raphson to convergence of each subsystem, try 1 iteration on the first subsystem followed by 1 iteration on the second subsystem. This may keep the subsystems more coupled, not going to "distant" unrelated sub solutions. Repeat this two step iteration and see if it converges.

Answer 3

This is a short answer and mostly a workaround. This could hit you back again later on, but can provide you with a temporary solution. You can try to relax the solution to your Newton-Raphson solver. What I mean by this is do not take the full solution of your non-linear system, but just take a factor of the previous solution and of the new solution (say 0.1 new solution, 0.9 old solution), and use this intermediate solution in your other solver. This means that you would be relaxing your global system of equation (and is similar conceptually to the above answer of using a single Newton-Raphson iteration).

This usually leads to more robust convergence, but it can turn out to be extremely slow.