Local time stepping for ALECG

[jbakosi]

Note: merge to branch vec. Need to get that one in first. Can you guys take a look at vec? Can we get vec in?

This adds the capability to ALECG to march to steady state using local time stepping.

By default this is off, and can be enabled by, e.g.,

inciter
  ...
  steady_state true
  residual 1.0e-8
  ...
end

Note that if marching to steady state is on, nstep and term are still considered when deciding on the end of time stepping: either of nstep, term, or the L2-norm of the residual, configured by residual above, reaching their respective condition, time stepping will stop.

No regression tests yet, but this is so far tested with vortical flow and the Onera M6 wing, for both this makes a difference in convergence rate when a final steady state can be assumed. Note that we do not deactivate any nodes, but local time stepping allows signals to propagate faster based on local time step sizes.

---

This change is 

[jbakosi]

Here is a convergence plot of the L2-norm of the residual with respect to the number of iterations, comparing global (red) with local (green) time stepping for the vortical flow problem. For this flow, it is not expected to make a large difference, but it still helps, and this also verifies correctness.

[adityakpandare]

Reviewed 30 of 42 files at r1, 3 of 5 files at r2.
Reviewable status: 32 of 42 files reviewed, 2 unresolved discussions (waiting on @adityakpandare, @charest, and @jbakosi)

---

src/Inciter/ALECG.hpp, line 304 at r1 (raw file):

    std::vector< tk::real > m_dtp;
    //! Physical time for each mesh node
    std::vector< tk::real > m_tp;

If it's a steady state calculation, why would we need to keep track of the local physical time?

---

src/Inciter/ALECG.cpp, line 732 at r2 (raw file):

      for (const auto& eq : g_cgpde) {
        auto eqdt = eq.dt( d->Coord(), d->Inpoel(), m_u );
        if (eqdt < mindt) mindt = eqdt;

If my understanding is right, the way to compute dt is the same, regardless if it's local or global. Only difference is, for global dt, we take a minimum. We should think about reusing that common code, if that's the case.

---

src/Inciter/ALECG.cpp, line 971 at r2 (raw file):

        else
          m_tp[i] = term;
      }

I don't think the physical time makes any sense if local time-stepping is used. I vote for removing this element-wise physical time bookkeeping.

[jbakosi]

Thanks for the extra set of eyes.

Reviewable status: 32 of 42 files reviewed, 2 unresolved discussions (waiting on @adityakpandare and @charest)

---

src/Inciter/ALECG.hpp, line 304 at r1 (raw file):

Previously, adityakpandare (Aditya Pandare) wrote…

If it's a steady state calculation, why would we need to keep track of the local physical time?

It's needed to evaluate source terms and Dirichlet BCs. See cg::CompFlow::src() and cg::CompFlow::dirbc().

---

src/Inciter/ALECG.cpp, line 732 at r2 (raw file):

Previously, adityakpandare (Aditya Pandare) wrote…

If my understanding is right, the way to compute dt is the same, regardless if it's local or global. Only difference is, for global dt, we take a minimum. We should think about reusing that common code, if that's the case.

Yeah, in general, you are right. Currently, we have two overloads of cg::CompFlow::dt(): the one that takes a minimum not only takes a minimum but also is written as an element loop, and also defines the local length scale differently (the cubic root of the tet), while the one that computes local time step sizes, is written as a node loop, and defines the local length scales as the volume of the dual mesh (associated to nodes). As a result, the two overloads look so different that I did not try combining them.

---

src/Inciter/ALECG.cpp, line 971 at r2 (raw file):

Previously, adityakpandare (Aditya Pandare) wrote…

I don't think the physical time makes any sense if local time-stepping is used. I vote for removing this element-wise physical time bookkeeping.

See above on why this is needed. However, you are right, that this line could probably be removed without any problems. This is likely still a remnant of me trying to understand what makes this algorithm work.

[adityakpandare]

Reviewable status: 32 of 42 files reviewed, all discussions resolved (waiting on @adityakpandare and @charest)

---

src/Inciter/ALECG.hpp, line 304 at r1 (raw file):

Previously, jbakosi (Jozsef Bakosi) wrote…

It's needed to evaluate source terms and Dirichlet BCs. See cg::CompFlow::src() and cg::CompFlow::dirbc().

Ah, I see. Thanks!

---

src/Inciter/ALECG.cpp, line 732 at r2 (raw file):

Previously, jbakosi (Jozsef Bakosi) wrote…

Yeah, in general, you are right. Currently, we have two overloads of cg::CompFlow::dt(): the one that takes a minimum not only takes a minimum but also is written as an element loop, and also defines the local length scale differently (the cubic root of the tet), while the one that computes local time step sizes, is written as a node loop, and defines the local length scales as the volume of the dual mesh (associated to nodes). As a result, the two overloads look so different that I did not try combining them.

Okay, that makes sense. Thanks for the explanation.

---

src/Inciter/ALECG.cpp, line 971 at r2 (raw file):

Previously, jbakosi (Jozsef Bakosi) wrote…

See above on why this is needed. However, you are right, that this line could probably be removed without any problems. This is likely still a remnant of me trying to understand what makes this algorithm work.

Yes, I understand now. Thanks

[jbakosi]

Thanks for the extra set of eyes!

Reviewable status: 31 of 42 files reviewed, all discussions resolved (waiting on @adityakpandare and @charest)

---

src/Inciter/ALECG.cpp, line 971 at r2 (raw file):

Previously, adityakpandare (Aditya Pandare) wrote…

Yes, I understand now. Thanks

Removed the above.

[adityakpandare]

Not a problem!

Reviewed 1 of 1 files at r3.
Reviewable status: 32 of 42 files reviewed, all discussions resolved (waiting on @adityakpandare and @charest)