Pressure projection methods

General Discussion
21

Thank you, professor. but there might be something wrong.

About comment 2. We can replace \bar{p} = 2*p^n - p^{n-1} with p^n, and this format is also formally second-order, but it is unstable. There is currently no proof of stability. There is also a rotational pressure correction method to avoid artificial Neumann BC maybe we can use it, but I am trying to update this code. Page 5 An overview of projection methods for incompressible flows.pdf (943.4 KB)

About comment 3. I understand what you mean. But I’m not sure if I can theoretically use the Dirichelt BC instead of the artificial Neumann BC. For the Neumann BC. I only need to add int1d(Th)(1e-8*pnp1*pp), but the Dirichlet BC requires on(1,2,3,4,pnp1 = pe(t)).

About comment 4. I understand roughly what you mean. I’m still trying to figure out how to fix or avoid this problem.

22

Ok, this is possible that the use of \bar p is unstable in general, even if it performs well on this test.

The Dirichlet BC is not usable in general since we do not know pe(t). Here it was just to check that it gives the same result as Neumann BC, so that we conclude that on this test there is no boundary layer due to the Neumann BC. Therefore this validates the Neumann BC.

About comment 4, there is no problem, the scheme is OK. This is just to explain that one has to be cautious on the procedure of evaluation of accuracy. One should not take into account the results for too small timestep, because when the timestep tends to 0 the error does not tend to 0, it tends to the error due to mesh size.
For a rigorous testing one should diminish timestep and mesh size simultaneously (but this can be very costly in computational time).

23

About comment 4. Does this require limiting the spatial step size? Professor. I always use Th = square(200,200), is this enough?

24

The right way is that each time you divide the timestep by 2, you also multiply the number of space points by 2.
But you can do this only when the H1 error attains its floor value.
Thus an appropriate experiment could be
dt=1/4, 1/8, 1/16 → nx=100
dt=1/32 → nx=200
dt=1/64 -->nx=400
dt=1/128–>nx=800
The last line is probably too expensive, thus you can stop at dt=1/64, nx=400.

25

I can understand that refining the time step and space step simultaneously can produce better results, but how to calculate the convergence order (time accuracy or spatiotemporal accuracy)?

At a fixed spatial step size, the way I calculate the temporal accuracy is as follows

if(i > 0){
        ConL2u[i] = log(errL2u(i-1)/errL2u(i))/log(dt(i-1)/dt(i));
        ConH1u[i] = log(errH1u(i-1)/errH1u(i))/log(dt(i-1)/dt(i));
        ConL2p[i] = log(errL2p(i-1)/errL2p(i))/log(dt(i-1)/dt(i));
    }

dt is the time step size.

26

If we assume that the error is comparable to C_t dt^{\alpha}+C_x h^{\beta}, we see that dividing simultaneously dt and h by 2, the error is divided by 2^{\min(\alpha,\beta)}.
It follows that evaluating the error by dividing the time and space steps, if we obtain a rate of \gamma, it implies that \alpha\geq\gamma and \beta\geq\gamma.
Thus you prove simultaneously the order at least \gamma for time and space.

If you want to prove the same order \gamma for both time and space this is done.
However if for example you want to prove order 1 for time and order 2 for space,
then you have to divide h by 2 and dt by 4, and verify that the error is divided by 4.

If you take h fixed and divide dt by 2, you see that the error C_t dt^{\alpha}+C_x h^{\beta}, decays only if C_t dt^{\alpha}>>C_x h^{\beta}, that is if dt is not too small. Otherwise if C_t dt^{\alpha}<<C_x h^{\beta} it will remain constant equal to C_x h^{\beta}. This is what you observe.