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]

Sure, why not. Later the class might have other methods alongside `evolve` like: * steady state * harmonic response * eigenanalysis but I think the transient response is the fundamental one so it's reasonable to make it the default. Envoyé depuis ProtonMail mobile

…

-------- Message d'origine -------- Le 27 déc. 2020 à 22:24, T. Gustafsson a écrit :

> 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 — You are receiving this because you authored the thread. Reply to this email directly, [view it on GitHub](#531 (comment)), or [unsubscribe](https://github.com/notifications/unsubscribe-auth/AAMD5U6PTXH4MBRMPBHW5GLSW4KNNANCNFSM4VHSAUKQ).

[gdmcbain]

I suppose if a time limit is the most common way to terminate a simulation, it might be reasonable to give evolve an optional argument tmax: float: np.inf. Charging while True to while t<tmax doesn't greatly increase the cyclomatic complexity.

[gdmcbain]

A default steady state method can be obtained from a backward Euler step of infinite length.

[gdmcbain]

One design decision is how to handle forcing. Is it part of the system or an input? It was implicitly left in the system above.but having it separate is often more physical and is more flexible.

[gdmcbain]

The whole nonlinear framework above does stumble at the moment because it relies on #119 having being solved which it hasn't despite considerable effort.

[gdmcbain]

Besides inputs, another important feature is outputs. Common uses:

• uncondensing: putting the condensed Dirichlet degrees of freedom back in
• untransforming: if the system has been transformed into the canonical form, it's nice to transform it back
  • higher, like second, order systems to first

These outputs and transformations should be composable. I have used toolz.functoolz.compose a lot for this, but like iterate, it's trivial to copy.

(I wonder whether one day Python's @ will grow up to be like Haskell's . operator. It looks the part. If matrices are thought of, as they should be, as linear operators, A @ B is already the composition.)

[gdmcbain]

I remembered over the break that back on 2019-05-07 I'd solved a textbook example exactly like the one asked about by @nikosmatsa . Here's the old source code, verbatim. N. B.:—This won't work as is!

"""

Reproducing Carslaw & Jaeger's 1959 figure 20: 'Temperatures in the slab -l < x < l with constant heat production at the rate A0 per unit volume and zero surface temperatures.  The numbers on the curves are values of kappa t / l**2.'

"""

from matplotlib.pyplot import subplots, pause
import numpy as np
from scipy.sparse import csr_matrix

from sksparse.cholmod import cholesky

import skfem
from skfem.models.poisson import laplace, mass, unit_load

from dysys import SparseDySys

Q_in = 1.0                      # mW
length = 911.0                  # mm
area = 17.0

n_lumps = 2**6 // 2

thermal_conductivity = 2.0      # mW / K mm
density = 3.0                   # mg / uL
specific_heat = 5.0             # uJ / mg K

mesh = skfem.MeshLine(np.linspace(0., length / 2, n_lumps + 1))
basis = skfem.InteriorBasis(mesh, skfem.ElementLineP1())

A0 = Q_in / area / length
f1 = skfem.asm(unit_load, basis)

def heating(*_) -> np.ndarray:
    return A0 * f1

sys = SparseDySys(
    density * specific_heat * skfem.asm(mass, basis),
    thermal_conductivity * skfem.asm(laplace, basis),
    heating).constrain([-1])

RC = (length / thermal_conductivity) * length * density * specific_heat
T_norm = A0 * length**2 / 8 / thermal_conductivity

t_final = 1.0 * RC / 4.0
dt = t_final / 1e2

def record(_, t, u, __):
    """put the nodal values in the dict and tag with time"""
    return u, {'time': t, 'temperature': sys.reconstitute(u)}

samples = [0.02, 0.06, 0.1, 0.2, 0.3, 0.4, 0.6, 0.8, 1.0]

temperature = {}
for t, _, d in sys.march_till(t_final + dt / 2, dt, d={},
                              events=[(tstar * RC / 4, record)
                                      for tstar in samples]):
    if d.get('time') == t:
        temperature[4 * t / RC] = d['temperature'] / T_norm

def exact(t, x, l, A0, k, alpha):
    return (1 - (x / l)**2 - 32 / np.pi**3 *
            sum((-1)**i / (2 * i + 1)**3 *
                np.cos((2 * i + 1) * np.pi * x / 2 / l) *
                np.exp(-alpha * (2 * i + 1)**2 * np.pi**2 * t / 4 / l**2)
                for i in range(100)))

fig, ax = subplots()
fig.suptitle('constant heat load to slab, by scikit-fem')
x_cj = np.linspace(0, 1) * length / 2
for t in temperature:
    line = ax.plot(2 * mesh.p[0] / length,
                   temperature[t],
                   label=f'{t:.2g}',
                   marker='o', linestyle='None')
    ax.plot(x_cj / length * 2,
            exact(t * RC / 4, x_cj,
                  0.5 * length, A0, thermal_conductivity,
                  thermal_conductivity / density / specific_heat),
            color=line[0].get_color())

ax.legend(loc=1)
ax.set_xlim((0, 1))
ax.set_xlabel(r'$x/\ell$')
ax.set_ylim((0, 1))
ax.set_ylabel(r'$2 k A T / \ell Q$')
fig.savefig('slab.png')

It won't work because:

• dysys is an unpublished in-house module, like what I've been talking about above
• sksparse.cholmod is long deprecated and getting difficult to install

I'll revise this to work with current master after I finish reviewing #530, but it might give some ideas. Here's the plot, like the one in

• Carslaw, H. S., Jaeger, J. C. (1959). Conduction of Heat in Solids. Oxford University Press

[nikosmatsa]

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

[gdmcbain]

This line

x = solve(*condense(L+M, b, I=m.interior_nodes()))

needs to be changed; skfem.utils.solve only solves a single linear algebraic system, it doesn't generate a trajectory for an initial value problem.

What's to be solved to take a step Δ t then is La = M a _old + b Δ t, so the first argument of condense should be L and the second the right-hand side, which includes both the actual generation term b and the thermal inertia terms from the old value of a. Since in your example (like ex19) the Dirichlet conditions and source term (absent in ex19) are independent of time, the condensation doesn't need to be done at every step. Anyway the solve needs to be called repeatedly with updated right-hand side; i.e. it needs to be iterated.

[nikosmatsa]

Based on ex19 what the "u_init" will be? 0?

[gdmcbain]

Is your initial condition zero? If so, you need an array of zeros, one for each degree of freedom. The number is given by the .N attribute of the InteriorBasis, as in

scikit-fem/docs/examples/ex13.py

Line 42 in 6e31b46

u = np.zeros(basis.N)

.

[nikosmatsa]

u = np.zeros(basis.N)
and my f(x,t)=10 will remain unchanged as i have defined it previously.
I will check it tomorrow then.

[gdmcbain]

Here's a minimal modification of ex19 to include steady uniform volumetric (areal) generation:

https://gist.github.com/gdmcbain/253b75c25af1daf286b2e66be21c9d8a

[nikosmatsa]

i can change in
s = asm(unit_load, basis)

the unit_load to loading

def loading(v, w):
return(10 * w.x[0] * v)

[gdmcbain]

s represents f(x,t)?

Yes, s for source.

[gdmcbain]

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

You do have to call evolve or step at some point. In ex19 this is done inside FuncAnimation. If you don't want to use that, use a for-loop, for example:

    for t, u in evolve(0.0, u_init):
        print(t, u.max())

[nikosmatsa]

`
def update(event):
t, u = event

    u0 = {'skfem': basis.interpolator(u)(np.zeros((2, 1)))[0]}
 
for t, u in evolve(0.0, u_init):
    plt.plot(t, u.max())
    plt.show()`

[nikosmatsa]

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

[nikosmatsa]

and making a 3d plot?

[nikosmatsa]

Traceback (most recent call last):
File "/opt/anaconda3/lib/python3.8/site-packages/matplotlib/cbook/init.py", line 196, in process
func(*args, **kwargs)
File "/opt/anaconda3/lib/python3.8/site-packages/matplotlib/animation.py", line 951, in _start
self._init_draw()
File "/opt/anaconda3/lib/python3.8/site-packages/matplotlib/animation.py", line 1743, in _init_draw
self._draw_frame(next(self.new_frame_seq()))
File "", line 66, in evolve
t, u = step(t, u)
File "", line 58, in step
u_new[interior] = backsolve(b1)
NameError: name 'backsolve' is not defined

[gdmcbain]

and making a 3d plot?

scikit-fem/docs/examples/ex10.py

Lines 47 to 49 in 6e31b46

	from skfem.visuals.matplotlib import plot3, show
	plot3(m, x)
	show()

[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