skfem.ivp

[gdmcbain]

From #529 at 2020-12-24:

Note that we don't have any helper functions for time integration yet. That would be a welcome improvement though, e.g., skfem.ivp or something.

[kinnala]

It's fine if most of the functionality are wrappers to stuff like https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.RK45.html but I think having a consistent interface (easy to experiment with the methods and parameters) which is specifically designed for finite element calculations (e.g. use of mass matrices in Mu'(t) = Au(t) and especially sparse matrices) would be beneficial.

[gdmcbain]

This is quite a big issue. I think getting the right abstraction could be very helpful but I fear that it's hard; not many packages get this right. Some of the big questions:

• The entire scipy.integrate.solve_ivp framework is useless here since it assumes the form y ' = f (y, t) which doesn't cover finite element systems with unlumped mass matrices.

• Even then, I find it restrictive that it requires t_span, the ‘interval of integration’; very often one would like to integrate until some condition is met, e.g. a steady state is obtained.

• An initial value problem framework should work together with other frameworks, for related problems like

  • stationary nonlinear problems (cf. scipy.optimize.root as a nonlinear solver #119)
    • continuation
      • natural parameter
      • pseudotransient
  • eigenproblems
  • ...

All of the above should really be outside the scope of scikit-fem but although I've been looking for a couple of decades, I haven't really found anything suitable. There's PETSc but it's a bit heavy. There was GarlicSim but it had already been abandoned before I found it.

[gdmcbain]

On the t_span question, I'm very much informed by the example at the end of the Revised Report on the Algorithmic Language Scheme, in particular that:

The value returned by integrate-system is an infinite stream of system states.

In Python this is naturally a generator, like ex19.evolve without the while-termination criterion; i.e. something like:

def evolve(t: float,
           u: np.ndarray) -> Iterator[Tuple[float, np.ndarray]]:
        t, u = step(t, u)
        yield t, u

As mentioned in #529, this is a lot like the iteration of the (problem-specific) function step which turns a continuous-time dynamical system into a discrete dynamical system. This kind of iteration is implemented in, e.g.:

• more_itertools.iterate
• toolz.itertoolz.iterate

I'm not sure why this isn't in itertools in the Python Standard Library; maybe because the implementation is literally just

    while True:
        yield x
        x = func(x)

(Here using this would involve packing the time t and the state u into x as a tuple or similar compound.)

[gdmcbain]

On the question of implicit dynamical systems, if one begins with the general implicit system defined by the function f(t, u, v) = 0 with v being the rate of change of u. Then the stationary problem is obtained by setting v=0 and a θ-approximation, e.g., by putting v = (u − u _old) / h and g (u) = θf (t, u, {u − u_old}/h) + (1 − θ) f (t − h, u_old, {u − u_old}/h).

This can be solved howsoever (#119). This often involves the jacobian j(u) = g'(u) which is

+
+

For linear systems, f(t, u, v) = M @ v + K @ u + s (t) this is just f_u = K and f_v = M.

In general the idea would be to define a dynamical system by the function f (the residual) and its two partial derivatives with respect to its second and third arguments; i.e. the stiffness and mass matrices, K and M.

[gdmcbain]

Solid mechanical problems often involve the second derivative with respect to time. These can always be reduced to larger first-order systems by appending the first derivative to the state. In practice, it is common though to treat the second-order system directly; e.g. as in the Newmark-beta method. I'm not sure whether the latter really offer much benefit. If so, the previous could be augmented by defining f(t, u, v, w), &c. In my own work I did implement things like Newmark but then switched to internally converting second-order systems to first since this simplified the code-base.

[gdmcbain]

In pacopy 0.1.2, as used in exx 23 and 27, a steady nonlinear system is defined by an object with .F and .jacobian_solver methods. The specification of the latter rather than just the jacobian is very useful because it enables:

• different problem-specific techniques to be specified for the solution (direct, iterative, preconditioned)
• cheating, as in not actually solving the exact Jacobian, as in the damping of the Newton iteration in ex10.

The closeness of this .F and .jacobian_solver methods to what I'd previously coded for nonlinear dynamical systems makes me wonder whether the interfaces can be combined.

Generally though I prefer functions to methods and I've had to replace pacopy with something taking residual and jacobian (or jacobian-solver) functions as arguments rather than an object with those methods; otherwise I find myself always having to construct temporary objects on the fly.

[gdmcbain]

wrappers to stuff like https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.RK45.html

I don't think this'll work will it? Does it assume dividing through by the mass matrix? As in converting Mu'(t) = Au(t) to u'(t) = M⁻¹ Au(t)? This is a possibility, but not if the mass matrix is singular, as in incompressible Stokes flow.

[gdmcbain]

Another consideration: coupling. Can one drive one system with the output of another? Can two or more systems be strongly coupled? The latter might lead to the concept of hierarchical systems.

[kinnala]

wrappers to stuff like https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.RK45.html

I don't think this'll work will it? Does it assume dividing through by the mass matrix? As in converting Mu'(t) = Au(t) to u'(t) = M−1 Au(t)? This is a possibility, but not if the mass matrix is singular, as in incompressible Stokes flow.

Alright, that may not make sense. Would it work by caching the LU-decomposition of M and write the RHS function handle using that or is it too slow? How about lumping the mass matrix, i.e. summing all off-diagonal entries to the diagonal after which the inverse is trivial? I'm not very familiar with time discretization as you may notice. :-)

[kinnala]

All of the above should really be outside the scope of scikit-fem but although I've been looking for a couple of decades, I haven't really found anything suitable. There's PETSc but it's a bit heavy. There was GarlicSim but it had already been abandoned before I found it.

If we are able to make interesting progress then we could eventually spin it off into a separate package.

[kinnala]

In Python this is naturally a generator, like ex19.evolve without the while-termination criterion; i.e. something like:

def evolve(t: float,
           u: np.ndarray) -> Iterator[Tuple[float, np.ndarray]]:
        t, u = step(t, u)
        yield t, u

So based on this, would the simple use of skfem.ivp be something like:

from skfem.ivp import evolve

for (k, u) in evolve(M, A, I=I, x=x, dt=1e-2):
    if k > 10:
        break
    # possibly plot u or check convergence criteria

Then various parameters to evolve could pick the time discretization method and its parameters?

[nikosmatsa]

https://github.com/nikosmatsa/Finite-Differences-Methods-for-PDE-/blob/main/parabolic%20fem%20.ipynb

[gdmcbain]

A brief comment before getting back to this properly in the new year.

Then various parameters to evolve could pick the time discretization method and its parameters?

Actually I generally leave the evolve function very simple. It's really just the same as more-itertools.iterate. All the machinery goes in the step function and the step function is passed to evolve. step can be a closure or a method. One implementation looks roughly like:

class DynamicalSystem:
  def step(h: float, t: float, x: ndarray) -> Tuple[float, ndarray]:
    raise NotImplementedError
  def evolve(h: float, t: float, x: ndarray):
     while True:
        yield t, x
        t, x = self.step(h, t, x)

class LinearDynamicalSystem(DynamicalSystem):
   M, L: spmatrix
   theta: float = 0.5

   def step(h, t, x):
       x = ... # take a generalized trapezoidal step as in ex19
       return t + h, x       

[gdmcbain]

One case in which the evolve function gets slightly more complicated is hybrid simulation, involving both continuous and discrete evolution. In this case, the tuple (t, x) is augmented to (t, x, d) where d is a dict or other object containing discrete dynamical variables which only change in instantaneous jumps at certain events, which might be scheduled in time or occur when the evolving continuous state x causes some bool predicate to flip. I have a lot of that sort of thing but we probably don't need to worry about it here.

[kinnala]

class DynamicalSystem:
  def step(h: float, t: float, x: ndarray) -> Tuple[float, ndarray]:
    raise NotImplementedError
  def evolve(h: float, t: float, x: ndarray):
     while True:
        yield t, x
        t, x = self.step(h, t, x)

class LinearDynamicalSystem(DynamicalSystem):
   M, L: spmatrix
   theta: float = 0.5

   def step(h, t, x):
       x = ... # take a generalized trapezoidal step as in ex19
       return t + h, x       

Should we have DynamicalSystem.__call__ = DynamicalSystem.evolve? So that one could use it like this:

for (t, x) in DynamicalSystem(...):
    if t > 1.:
        break

Edit: I suppose it doesn't matter because

for (t, x) in DynamicalSystem(...).evolve(...):
    if t > 1.:
        break

is just as simple and maybe even more descriptive.

[gdmcbain]

NameError: name 'backsolve' is not defined

Yeah, backsolve is pretty important, you need that.

scikit-fem/docs/examples/ex19.py

Line 73 in 6e31b46

backsolve = splu(condense(A, D=boundary, expand=False).T).solve

[nikosmatsa]

`if name == 'main':

from argparse import ArgumentParser
from pathlib import Path

from matplotlib.animation import FuncAnimation
import matplotlib.pyplot as plt

from skfem.visuals.matplotlib import plot


ax = plot(m, u_init, shading='gouraud')
title = ax.set_title('t = 0.00')
field = ax.get_children()[0]  # vertex-based temperature-colour
fig = ax.get_figure()
fig.colorbar(field)

def update(event):
    t, u = event

    u0 = {'skfem': basis.interpolator(u)(np.zeros((2, 1)))[0]}
   

    title.set_text(f'$t$ = {t:.2f}')
    field.set_array(u)

animation = FuncAnimation(fig, update, evolve(0., u_init), repeat=False)
plt.show()


from skfem.visuals.matplotlib import plot3, show 
plot3(?,?) 
show() `

[nikosmatsa]

plt3(what, what)

[gdmcbain]

As in ex10 cited above, the two arguments are the mesh and the solution. (Or in this unsteady case a snapshot of the solution.)

[nikosmatsa]

`NameError Traceback (most recent call last)
in
30
31 from skfem.visuals.matplotlib import plot3, show
---> 32 plot3(m,u)
33 show()

NameError: name 'u' is not defined
`

[nikosmatsa]

How can I change the theta rule method in order to implement to the CN or backward
Euler method for the time derivative in the code :

Ba^{n}+\frac{1}{2}t Aa^{n}=Ba^{n-1}-\frac{1}{2}t Aa^{n-1}+tb^{n-1/2}

or

Ba^{n} =& Ba^{n-1}-t A^{n} + t b^{n}

[gdmcbain]

NameError: name 'u' is not defined

In ex10 the mesh is called m and the unknown deflection x; in ex19 the mesh is called mesh and the unknown temperature at any given time step is called u. The names are more or less arbitrary, as usual in mathematics and computation.

[nikosmatsa]

https://github.com/nikosmatsa/Finite-Differences-Methods-for-PDE-/blob/main/parabolic%20fem%20.ipynb

[gdmcbain]

How can I change the theta rule method in order to implement to the CN or backward
Euler method for the time derivative

The theta method is already implemented in ex19; theta = 1/2 gives Crank–Nicholson, theta = 1 gives backward Euler.

scikit-fem/docs/examples/ex19.py

Lines 65 to 67 in 6e31b46

	theta = 0.5 # Crank–Nicolson
	A = M + theta * L * dt
	B = M - (1 - theta) * L * dt

[nikosmatsa]

I didn't put it right.I don't want theta rule

[gdmcbain]

Crank–Nicolson and backward Euler are just special cases of the θ-rule. If you want to hard-code them, you could; e.g.

A = M + L * dt / 2
B = M - L * dt / 2

or

A = M + L * dt
B = M

respectively.

[nikosmatsa]

`if name == 'main':

import matplotlib.pyplot as plt
from skfem.visuals.matplotlib import plot3, show 


def update(event):
    t, u = event

for t, u in evolve(0.0, u_init):
    plot3(m,u)
    plt.show() `

plots me all the graphs but I want only the last one

[gdmcbain]

Try moving the plotting outside the for-loop:

for t, u in evolve(0.0, u_init):
    pass

plot3(m,u)
plt.show()

[gdmcbain]

Oh, I don't think you want plt.show since you've already imported show from skfem.visuals.matplotlib, so just

show()

[nikosmatsa]

close the issue.thnka for everything.you are awesome