Theta method for DAEs

[oriolcg]

Hi @santiagobadia, there is a problem when working with solutions of a DAE system solved using the theta-method. The values of the constraining field at n+1 are wrong when 0 < \theta < 1.

For instance for a Navier-Stokes problem, currently we have the following approach:

1. Solve for (u_\theta,p_\theta)
2. Update the solution (u_np1,p_np1) = 1/\theta * (u_\theta,\p_theta) + (1-\theta)/\theta * (u_n,p_n)

Since p is an instantaneous field that balances the equations at n+\theta, this approach leads to values of p_{n+1} that don't satisfy the momentum equation at n+1. This is not a big deal if you only want to work with the velocities, but if you want to output forces, e.g. drag and lift, the stress at n+1 computed from (u_{n+1},p_{n+1}) is wrong.

One solution would be to update p at n+1 solving the pressure Poison equation at the driver level (expensive). Another way would be to output the solution where it satisfies the equations, i.e. at n+\theta. What do you think? Would it be possible to return the values at n+\theta and {n+1}? In that way you could compute the quantity of interest at n+\theta.

We could also consider sending a functional to be evaluated inside GridapODEs, but I don't know if this is something contemplated in the GridapODEs philosophy...

[santiagobadia]

Hi @oriolcg

I don't think this is an error. I would say this is pretty common, the theta method return the solution at the end of the step, not intermediate stages.

I don't think that the drag for u^n+1,p^n+1 is wrong, I think that you can only expect first order accuracy for the pressure, and thus, it is as good approximation p^n+1 and p^n+\theta.

I would not solve for a pressure Poisson equation, it introduces artificial boundary conditions, etc. It is a mess, I would skip this. It is inherently wrong at the Dirichlet boundaries, and 100% sure it is much worse that p^n+1 when computing drags.

We probably agree that using the theta method and returning the value of the solution at t^n+\theta is not standard. People will expect the value at n+1, which you need anyway at the next step.

My philosophy is to keep Gridap ODEs as simple and clean as possible. If people want to do non-standard things, I would keep it at the driver level.

So, I would consider two possible options:

1. The easiest and probably best option is to do this at the driver level, you just need to store pn, get pn+1 from the solver and compute p^{n+\theta} = \theta*p^{n+1} + (1-\theta)*p^{n}. I would even consider some syntactic sugar, and define +,- binary operations over FEFunctions and the unary operation * a scalar. These are three line functions. (To be done in Gridap, if not already there, I should check.)

2. Add a kwargs-based constructor for the ThetaMethod solver that provides information like whether you want to return the solution at a step different from n+1 (using n+1 as default). But I would probably explore 1. and if not convinced go for 2. Again, it is 3 lines of code in GridapODEs.

What do you think?

[oriolcg]

Hi @santiagobadia, thanks for the comments. I didn't mean to say that the theta method is wrong, just that for this particular case you can have artificial pressure oscillations in time due to this linear extrapolation. I agree with you that returning the results at n+1 is, and should be, the standard of practice.

Actually, I've been using the option 1 in all my drivers (NSI and FSI) and it works well. As you say, it's only a couple of lines of code since all the operators +,-,* can operate with FEFunctions. I just wanted to see if there is a more efficient way, but given that this issue does not appear for other methods (BDF, RK, ...), it's cleaner to keep everything at the driver level.

Closing issue.