Consistently evaluate Neumann BC in generalized-alpha method

[aabrown100-git]

As I understand the generalized-alpha method (see Chung and Hulbert, 1993) if we have, for example, a time-dependent pressure boundary condition, then this BC contributes to the residual vector using the value of pressure at $t_{n+\alpha_f}$ (or $t_{n+1-\alpha_f}$ depending on the definition of $\alpha_f$). Looking through the code (starting from SETBCNEU()), it appears the value of pressure at $t_{n+1}$ is being used instead. Is this is an error in the code, or else please let me know where I am mistaken.

[vvedula22]

Yes, the paper states that the applied load should be at $n+\alpha_f$ but we compute the load at $n+1$. This variant is shown in Bazilevs et al. (2008) where even the pressure as a state variable is evaluated at $n+1$.

I don't think this makes a noticeable difference in the solution or convergence as the applied load is not a function of the state variable and this term doesn't contribute to LHS. In the case of the follower pressure load, the contribution to LHS is evaluated at $n+\alpha_f$ but the actual load is still at $n+1$. It is more important to have the state variable terms consistently evaluated at $n+\alpha_m$ or $n+\alpha_f$. Nevertheless, you could try and report how it goes.

[aabrown100-git]

Thanks for that paper! I see that pressure as a state variable for the fluid system is evaluated at $n+1$ (eqs. 156-158), but it doesn't seem say when the value of the Neumann BC (like an applied surface pressure) for the solid system should be evaluated. These are two different quantities, right?

I agree that it probably won't make much of a difference. I'm thinking that if this really is an error, then the effect is that we are applying a slightly higher or lower pressure to the solid than we think we are.

[aabrown100-git]

Adding a relevant recent paper A note on the accuracy of the generalized- scheme for the incompressible Navier-Stokes equations. This paper discusses whether to use $p_{n+1}$ or $p_{n+\alpha_f}$ in the NS residual, and also whether to evaluate the Neumann BC at $t_{n+1}$ or $t_{n+\alpha_f}$. See in particular the introduction and Eq. (11) onwards. I believe we are currently performing Scheme 3. The paper recommends Scheme 2, but shows that Scheme 1 (which I believe is what I am suggesting) is better than Scheme 3. Of course, this is for the Navier-Stokes equations, so I'm not sure if these results directly extend to struct and ustruct equations.

[vvedula22]

@aabrown100-git Earlier our pressure integration was based on the Bazilevs 2008 paper which is different from Scheme 1 but I changed it to Scheme 2 later. Reg. the pressure boundary condition, let us change it to evaluate it at $n+\alpha_f$.

[alisonmarsden]

Adding my two cents here, Ingrid and Ju’s paper demonstrates a nice increase in accuracy in their paper (the one that Aaron linked to) and I think it would be great to update svFSI to use this method – this could also be done after the C++ conversion, and could be provided as an option to the user. From: Vijay Vedula ***@***.***> Date: Tuesday, October 4, 2022 at 10:03 AM To: SimVascular/svFSI ***@***.***> Cc: Subscribed ***@***.***> Subject: Re: [SimVascular/svFSI] Error in Neumann BC residual contribution in generalized-alpha method? (Issue #88) @aabrown100-git<https://github.com/aabrown100-git> Earlier our pressure integration was based on the Bazilevs 2008 paper which is different from Scheme 1 but I changed it to Scheme 2 later. Reg. the pressure boundary condition, let us change it to evaluate it at $n+\alpha_f$. — Reply to this email directly, view it on GitHub<#88 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ADJGJZ5N55HDEKIQ7NNXVPDWBPQELANCNFSM6AAAAAAQUHAPEU>. You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[aabrown100-git]

In terms of implementation in svFSI, I see two options. Denoting the Neumann BC value $\mathbf{h}(t)$, we can

1. Evaluate $\mathbf{h}(t_{n+\alpha_f})$ and use that value when constructing the residual
2. Evaluate $\mathbf{h}(t_{n})$ and $\mathbf{h}(t_{n+1})$ and interpolate between these value to get $\mathbf{h}(t_{n+\alpha_f})$, then construct the residual

I'm not sure if there is any meaningful difference between the two, but one may be easier to implement than the other. Option 2 may be more appropriate when the Neumann BC value is obtained from an external LPN code.