jacobian block for implicit solver with non explicit flux function

naffrancois · February 15, 2018, 15:45

Hello,

I am coding implicit finite volume. I have a doubt on how to construct a jacobian block associated to a face boundary which does not have an explicit form. The flux is of the form F(Q*) where Q* is the solution to a non linear equation solved numerically with Newton and which depends on Qint. I need dF/dQint. I'd like to avoid local finite differencing and stick to formal differentiation. Anyone has an idea ? I guess this situation appears as well for the differentiation of the Godunov flux with exact Riemann solver, if someone ever tried to do this.

Thanks !

sbaffini · February 16, 2018, 10:49

If you think about wall functions, I guess, they are not much different from what you want to do.

Only, note that, typically, one does not differentiate with respect to Qint (assuming with Qint you refer to the variable on the face), but with respect to the variable in the adjacent fluid cells (that is, with respect to your actual independent variables).

In general, with WF, you would tipically express them as a classical, no WF, linear contribution + a correction with respect to it, and differentiate with respect to the linear contribution only.

This general approach is know as deferred correction and is also used for convection schemes, anything of order higher than 1 (including diffusion), whatever you would be otherwise forced to discretize explicitly.

Obviously, this only makes sense if the part you account for implicitly is somehow still related to the whole term.

Note that, when implemented within a straight newton linearization, what you actually do is:

- just use the full intended term in the explicit residual on the RHS
- use your best approximation in the jacobian on the LHS

and at convergence, whatever the approximation in the jacobian was, it has no more effect.

naffrancois · February 16, 2018, 16:04

Many thanks Sbaffini for your explanations. I am going to reformulate, to check if I understood.

Let's assume a cell centered 1st order FV to make things simpler. At the face boundary I have the constant average state Qint coming from the adjacent cell in the computational domain and known independant boundary data Qbc. The boundary condition consists of non linear relations, which give me a face vector Qf through the solution of a problem G(Qf,Qint,Qbc)=0. Then the final flux is computed F(Qf) and added to the RHS.

Then you suggest linearizing Qf about Qint such that Qf=Qint+other terms. Then F(Qf) is approximately F(Qint) and take dF(Qint)/dQint as the Jacobian block contribution. This is a great complexity reduction but I am wondering how much information I loose in the process.

Thank you for your time

praveen · February 16, 2018, 23:54

If

and

is obtained from $Q_{int}$ by solving $N(Q^*,Q_{int})=0$ then

$\frac{\partial F}{\partial Q_{int}} = \frac{\partial F}{\partial Q^*} \frac{\partial Q^*}{\partial Q_{int}}$

and

$\frac{\partial N}{\partial Q_{int}} + \frac{\partial N}{\partial Q^*} \frac{\partial Q^*}{\partial Q_{int}} = 0$

Hence

$\frac{\partial F}{\partial Q_{int}} = -\frac{\partial F}{\partial Q^*} \left[ \frac{\partial N}{\partial Q^*} \right]^{-1} \frac{\partial N}{\partial Q_{int}}$

You need a formula for

to compute its partial derivatives. If its too complicated, you can use automatic differentiation.

sbaffini · February 17, 2018, 06:05

As suggested by praveen, you actually need dF/dQint itself as Jacobian block contribution, NOT dF/dQf (Q* in praveen answer), and this block contribution goes to the cell relative to Qint.

Still, while the form used by praveen is the mathematically correct one, I guess you should take control of dN/dQ*, in case it might become 0.

Remember that the whole point in making things implicit is to then have a nice matrix, even if not accurate.

naffrancois · February 18, 2018, 17:38

Many thanks for your answers !

As a matter of fact, I use automatic differentiation to propagate derivatives with respect to the conservative variables. As I implemented this procedure with the goal of not worrying about complex flux functions or eos, I was having a doubt about its usability in this precise scenario where the function is not explicit. In this case, the derivatives are propagated through the whole Newton iterative scheme and I am wondering if I get at convergence (based on the norm of N) something even close to the analytic expression given by praveen. I guess I would have to derive by hand these expressions, or more safely use formal calculus software. I don't feel very confident in using finite differencing for checking purpose would you ?

Sbaffini: if dN/dQ*=0 then the Newton scheme may stall to a local extremum, in that case I switch to a bracketed dichotomy algorithm

sbaffini · February 19, 2018, 05:52

What I meant is that, in the solution provided by praveen for dF/dQint, it might happen that your solution for N=0 ends up having also dN/dQ* = 0. You can use a bisection method to find the solution for N=0, but what about (dN/dQ*)^-1 in the formula for dF/dQint if that is the actual correct solution?

Of course, your problem can be such that a correct solution to N=0 never happen to also have dN/dQ*=0, but noting that you have a backup algorithm, I suspect that might not be the case.

Also, not an expert here, but I guess that automatic differentiation algorithm would exactly replicate any unboundedness in the original jacobian; do you have a backup also for that?

praveen · February 19, 2018, 06:01

You can solve for

from $N(Q^*, Q_{int})=0$ only if $\frac{\partial N}{\partial Q^*}$ is invertible. This is impicit function theorem

https://en.wikipedia.org/wiki/Implicit_function_theorem

If that jacobian is singular, that means your problem

is ill-posed and you have to fix it at a mathematical/modeling level. It is not a numerical issue.

sbaffini · February 19, 2018, 07:23

Well, ill-posedness happens nonetheless. Also, quoting from the link, the implicit function theorem gives a sufficient condition to ensure that there is such a function.

It is plenty of nonlinear problems whose solutions are in points of null derivatives. I don't think you can arbitrarily ban them from the physics just because of that. I think it is our math that should put them in a well posed manner if possible at all, or try to exclude those ill-posed solutions from the set of reachable ones.

February 15, 2018, 15:45	jacobian block for implicit solver with non explicit flux function	#1
naffrancois Senior Member Join Date: Oct 2011 Posts: 242 Rep Power: 17	Hello, I am coding implicit finite volume. I have a doubt on how to construct a jacobian block associated to a face boundary which does not have an explicit form. The flux is of the form F(Q) where Q is the solution to a non linear equation solved numerically with Newton and which depends on Qint. I need dF/dQint. I'd like to avoid local finite differencing and stick to formal differentiation. Anyone has an idea ? I guess this situation appears as well for the differentiation of the Godunov flux with exact Riemann solver, if someone ever tried to do this. Thanks !

February 16, 2018, 23:54		#4
praveen Super Moderator Praveen. C Join Date: Mar 2009 Location: Bangalore Posts: 343 Blog Entries: 6 Rep Power: 18	If and is obtained from $Q_{int}$ by solving $N(Q^,Q_{int})=0$ then $\frac{\partial F}{\partial Q_{int}} = \frac{\partial F}{\partial Q^} \frac{\partial Q^}{\partial Q_{int}}$ and $\frac{\partial N}{\partial Q_{int}} + \frac{\partial N}{\partial Q^} \frac{\partial Q^}{\partial Q_{int}} = 0$ Hence $\frac{\partial F}{\partial Q_{int}} = -\frac{\partial F}{\partial Q^} \left[ \frac{\partial N}{\partial Q^*} \right]^{-1} \frac{\partial N}{\partial Q_{int}}$ You need a formula for to compute its partial derivatives. If its too complicated, you can use automatic differentiation. sbaffini likes this. __________________ http://cpraveen.github.io http://twitter.com/cfdlab http://github.com/cpraveen

February 18, 2018, 17:38		#6
naffrancois Senior Member Join Date: Oct 2011 Posts: 242 Rep Power: 17	Many thanks for your answers ! As a matter of fact, I use automatic differentiation to propagate derivatives with respect to the conservative variables. As I implemented this procedure with the goal of not worrying about complex flux functions or eos, I was having a doubt about its usability in this precise scenario where the function is not explicit. In this case, the derivatives are propagated through the whole Newton iterative scheme and I am wondering if I get at convergence (based on the norm of N) something even close to the analytic expression given by praveen. I guess I would have to derive by hand these expressions, or more safely use formal calculus software. I don't feel very confident in using finite differencing for checking purpose would you ? Sbaffini: if dN/dQ*=0 then the Newton scheme may stall to a local extremum, in that case I switch to a bracketed dichotomy algorithm sbaffini likes this.

February 19, 2018, 06:01		#8
praveen Super Moderator Praveen. C Join Date: Mar 2009 Location: Bangalore Posts: 343 Blog Entries: 6 Rep Power: 18	You can solve for from $N(Q^, Q_{int})=0$ only if $\frac{\partial N}{\partial Q^}$ is invertible. This is impicit function theorem https://en.wikipedia.org/wiki/Implicit_function_theorem If that jacobian is singular, that means your problem is ill-posed and you have to fix it at a mathematical/modeling level. It is not a numerical issue. __________________ http://cpraveen.github.io http://twitter.com/cfdlab http://github.com/cpraveen

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
[Other] mesh airfoil NACA0012	anand_30	OpenFOAM Meshing & Mesh Conversion	13	March 7, 2022 18:22
[ANSYS Meshing] Help with element size	sandri_92	ANSYS Meshing & Geometry	14	November 14, 2018 08:54
[blockMesh] non-orthogonal faces and incorrect orientation?	nennbs	OpenFOAM Meshing & Mesh Conversion	7	April 17, 2013 06:42
[blockMesh] error message with modeling a cube with a hold at the center	hsingtzu	OpenFOAM Meshing & Mesh Conversion	2	March 14, 2012 10:56
Droplet Evaporation	Christian	Main CFD Forum	2	February 27, 2007 07:27

February 16, 2018, 10:49		#2
sbaffini Senior Member Paolo Lampitella Join Date: Mar 2009 Location: Italy Posts: 2,195 Blog Entries: 29 Rep Power: 39	If you think about wall functions, I guess, they are not much different from what you want to do. Only, note that, typically, one does not differentiate with respect to Qint (assuming with Qint you refer to the variable on the face), but with respect to the variable in the adjacent fluid cells (that is, with respect to your actual independent variables). In general, with WF, you would tipically express them as a classical, no WF, linear contribution + a correction with respect to it, and differentiate with respect to the linear contribution only. This general approach is know as deferred correction and is also used for convection schemes, anything of order higher than 1 (including diffusion), whatever you would be otherwise forced to discretize explicitly. Obviously, this only makes sense if the part you account for implicitly is somehow still related to the whole term. Note that, when implemented within a straight newton linearization, what you actually do is: - just use the full intended term in the explicit residual on the RHS - use your best approximation in the jacobian on the LHS and at convergence, whatever the approximation in the jacobian was, it has no more effect.

February 16, 2018, 16:04		#3
naffrancois Senior Member Join Date: Oct 2011 Posts: 242 Rep Power: 17	Many thanks Sbaffini for your explanations. I am going to reformulate, to check if I understood. Let's assume a cell centered 1st order FV to make things simpler. At the face boundary I have the constant average state Qint coming from the adjacent cell in the computational domain and known independant boundary data Qbc. The boundary condition consists of non linear relations, which give me a face vector Qf through the solution of a problem G(Qf,Qint,Qbc)=0. Then the final flux is computed F(Qf) and added to the RHS. Then you suggest linearizing Qf about Qint such that Qf=Qint+other terms. Then F(Qf) is approximately F(Qint) and take dF(Qint)/dQint as the Jacobian block contribution. This is a great complexity reduction but I am wondering how much information I loose in the process. Thank you for your time

February 17, 2018, 06:05		#5
sbaffini Senior Member Paolo Lampitella Join Date: Mar 2009 Location: Italy Posts: 2,195 Blog Entries: 29 Rep Power: 39	As suggested by praveen, you actually need dF/dQint itself as Jacobian block contribution, NOT dF/dQf (Q* in praveen answer), and this block contribution goes to the cell relative to Qint. Still, while the form used by praveen is the mathematically correct one, I guess you should take control of dN/dQ*, in case it might become 0. Remember that the whole point in making things implicit is to then have a nice matrix, even if not accurate.

February 19, 2018, 05:52		#7
sbaffini Senior Member Paolo Lampitella Join Date: Mar 2009 Location: Italy Posts: 2,195 Blog Entries: 29 Rep Power: 39	What I meant is that, in the solution provided by praveen for dF/dQint, it might happen that your solution for N=0 ends up having also dN/dQ* = 0. You can use a bisection method to find the solution for N=0, but what about (dN/dQ)^-1 in the formula for dF/dQint if that is the actual correct solution? Of course, your problem can be such that a correct solution to N=0 never happen to also have dN/dQ=0, but noting that you have a backup algorithm, I suspect that might not be the case. Also, not an expert here, but I guess that automatic differentiation algorithm would exactly replicate any unboundedness in the original jacobian; do you have a backup also for that?

February 19, 2018, 07:23		#9
sbaffini Senior Member Paolo Lampitella Join Date: Mar 2009 Location: Italy Posts: 2,195 Blog Entries: 29 Rep Power: 39	Well, ill-posedness happens nonetheless. Also, quoting from the link, the implicit function theorem gives a sufficient condition to ensure that there is such a function. It is plenty of nonlinear problems whose solutions are in points of null derivatives. I don't think you can arbitrarily ban them from the physics just because of that. I think it is our math that should put them in a well posed manner if possible at all, or try to exclude those ill-posed solutions from the set of reachable ones.