method of characteristics

[gdmcbain]

The method of characteristics is tailored to solve unsteady advection equations

(Ern & Guermond 2004, Theory and Practice of Finite Elements, §6.3.3 ‘The method of characteristics’).

(This is a stub which I'll fill out shortly, but for now to refer to from the discussion of the evaluation matrix #593. I expect it to be a heavy load on that or probes #572.)

[kinnala]

Should it be solved like this: 1) find the points that map to the quadrature points through the advection map, 2) find out which of them are outside of the domain, 3) evluate rest of them through probes.

Then we could have a separate function which checks if a point is inside the domain or not.

[gdmcbain]

Following the discussion in #646, in particular your question

Is the time step so large that the point can move over several elements?

which is called in the finite difference context the Courant number u Δ t / Δ x or CFL condition after Courant, Friedrich, & Lewy (1928, 1967), I'm thinking that (3) should use a modified version of probes that is already told which cell each point is in and so doesn't have to call element_finder as in

scikit-fem/skfem/assembly/basis/cell_basis.py

Line 146 in d08af2a

cells = self.mesh.element_finder(mapping=self.mapping)(*x)

To begin with, one might check using the local-to-canonical map whether the point is in the guessed cell. For those that aren't, try the adjacent cell. I haven't looked at the details of what the local-to-canonical map looks like for a point not in the cell to see whether that might guide which facet has been crossed and so which of the adjacent cells is to be taken as the next guess, but if there's something there, it be better to use it than to fall back on an exhaustive search of the neighbours.

The mesh's facet connectivity array already indicates which facets are on the boundary of the domain, via the sentinel value −1.

scikit-fem/skfem/mesh/mesh.py

Line 670 in 0855dae

inverse[1, np.nonzero(inverse[0] == inverse[1])[0]] = -1

—

• Courant, R., Friedrichs, K. & Lewy, H. (1967). On the Partial Difference Equations of Mathematical Physics. IBM Journal, 11, 215--234.

[gdmcbain]

On the other hand, in examples like

• exx 25, 28 (simple channels)
• ex 27 (backwardfacing step)

it's easy to know which feet are upstream of the inlet because the mathematical description of the inlet is so simple (x=const.).

That shortcut might arise often enough to be worth exploiting. That would take care of

2. find out which of them are outside of the domain

after which simple application of (3) to the remainder should be valid. (If not optimal, but of course premature optimization is to be avoided.)

[gdmcbain]

A very simple one-dimensional, uniform advecting velocity, zero diffusion, example is discussed in ‘Is this a good way to implement the method of characteristics for advection equation?’ (FEniCS Q&A, 2017-03-19).

[kinnala]

There is also an efficient "in_polygon" implementation in matplotlib (Path class). There are also algorithms for 3D meshes, e.g. https://github.com/marmakoide/inside-3d-mesh