simulkade / FVTool

Finite volume toolbox for Matlab/Octave
http://fvt.simulkade.com
BSD 2-Clause "Simplified" License
97 stars 56 forks source link

Implicit vs explicit Euler #16

Closed aao24 closed 6 years ago

aao24 commented 6 years ago

Hi Ali,

I was wondering why, while using backward Euler in time, the convection term obtained by TVD usees the pervious iterate (phi_old)... Would it be more natural to use explicit Euler?

[See e.g. ConvectionTVDBurgers and diffconv_variable_diff in Examples]

Thanks, Anna

simulkade commented 6 years ago

Hi Anna, The implicit Euler and TVD scheme make the discretized convection term nonlinear. Therefore, I calculate the nonlinear term at the previous time step (as an initial estimate) and move it to the right hand side, then iterate a few times. You can also use it with the explicit Euler scheme. I'll add an example later this week. Most of the time, I need to have large time steps; therefore I prefer to use the implicit scheme.

aao24 commented 6 years ago

Thanks for you promt reply, Ali!

I understand. I guess this is the common practice for the flux limiting convection schemes to use this semi-implicit version? Is this method has the same stability and accuracy properties as the fully implicit version?

Assuming the time step is small enough, do both methods still have the same accuracy? I attached a test code where I use backward and forward Euler for a small enough time step to ensure stability (dt=1e-03)...takes some time! BE.gif and FE.gif are the realizations, and relative_abs_diff.gif is the abs difference between them (devided by 2-norm of one of the solutions)...

fe be relative_abs_diff

ConvDiff_where_u_depends_on_t.txt

simulkade commented 6 years ago

Hi Anna,

What happens if you iterate in the implicit scheme:

ci_temp = ci_old;
for i = 1:3
  [Mc_i, RHSc_i] = convectionTvdTerm2D(W_i, ci_temp, FL);
  Mi= Mc_i+Mt_i+Mbc-Mdiff;
  RHSi = RHSt_i+RHSbc+RHSc_i;
  %adding the point source
  RHSi(xiX+(xiY-1)*(Nx+2))=RHSi(xiX+(xiY-1)*(Nx+2))+a;
  %solving, implicity
  ci_temp = solvePDE(m, Mi, RHSi);
end
ci = ci_temp
aao24 commented 6 years ago

I get the following relative error with the inner loop relative_abs_diff2

mm... looks too similar... so here is the difference of two "implicit scheme" solutions (scaled by the explicit solution norm as above)

bevsbe2

In the plot below $\phi_1$ is the solution with the implicit Euler without the inner loop, and $\phi_2$ with the inner loop, $\phi_e$ corresponds to the explicit Euler. norms

simulkade commented 6 years ago

It is indeed common practice to use the semi implicit version. When I was writing the code, I looked at the source code of a big PDE solver project that a friend of mine was involved in, and it turned out that they use the semi-implicit method without even iterating. I have tested in on a few problems, and it seems that 3-5 iterations is enough for the inner loop to converge. However, I must admit there has been a few cases (coupled multi component multi phase flow in porous media with phase equilibrium) that this method was not stable enough; the overshoot was very tiny, but enough to blow up the whole system of equations.
I think the best way to compare the accuracy of these method is to solve a simple convection equation (c_t + u c_x = 0) in 1D and compare the numerical and analytical solutions for different time steps. BTW, the animations you've created look very cool :-)

aao24 commented 6 years ago

Thank you! :-)