Closed gdmcbain closed 3 years ago
kinnala
commented
3 years ago 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
commented
3 years ago 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
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
commented
3 years ago 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-systemis 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, uAs 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.:
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
commented
3 years ago 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 − uold}/h) + (1 − θ) f (t − h, uold, {u − uold}/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 fu = K and fv = 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
commented
3 years ago 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
commented
3 years ago 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:
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
commented
3 years ago 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.
gdmcbain
commented
3 years ago 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
commented
3 years ago 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
commented
3 years ago 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
commented
3 years ago In Python this is naturally a generator, like
ex19.evolvewithout thewhile-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 criteriaThen various parameters to evolve could pick the time discretization method and its parameters?
nikosmatsa
commented
3 years ago gdmcbain
commented
3 years ago 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 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
commented
3 years ago 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.:
breakEdit: I suppose it doesn't matter because
for (t, x) in DynamicalSystem(...).evolve(...):
if t > 1.:
breakis just as simple and maybe even more descriptive.
gdmcbain
commented
3 years ago Sure, why not.
Later the class might have other methods alongside evolve like:
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, or unsubscribe.
gdmcbain
commented
3 years ago 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
commented
3 years ago A default steady state method can be obtained from a backward Euler step of infinite length.
gdmcbain
commented
3 years ago 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
commented
3 years ago 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
commented
3 years ago Besides inputs, another important feature is outputs. Common uses:
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
commented
3 years ago 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 abovesksparse.cholmod is long deprecated and getting difficult to installI'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
nikosmatsa
commented
3 years ago gdmcbain
commented
3 years ago 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
commented
3 years ago Based on ex19 what the "u_init" will be? 0?
gdmcbain
commented
3 years ago 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 https://github.com/kinnala/scikit-fem/blob/6e31b46e1d11d12d93798405d38c3764a18a8f72/docs/examples/ex13.py#L42.
nikosmatsa
commented
3 years ago 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
commented
3 years ago Here's a minimal modification of ex19 to include steady uniform volumetric (areal) generation:
https://gist.github.com/gdmcbain/253b75c25af1daf286b2e66be21c9d8a
nikosmatsa
commented
3 years ago i can change in s = asm(unit_load, basis)
the unit_load to loading
def loading(v, w): return(10 w.x[0] v)
nikosmatsa
commented
3 years ago s represents f(x,t)?
nikosmatsa
commented
3 years ago gdmcbain
commented
3 years ago i can change in
s = asm(unit_load, basis)the unit_load to loading
def loading(v, w): return(10 * w.x[0] * v)
Yes, if you like, but don't forget to decorate it as a @LinearForm, as in the original:
gdmcbain
commented
3 years ago s represents f(x,t)?
Yes, s for source.
gdmcbain
commented
3 years ago 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())` 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()`and making a 3d plot?
nikosmatsa
commented
3 years ago 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 "
gdmcbain
commented
3 years ago gdmcbain
commented
3 years ago NameError: name 'backsolve' is not defined
Yeah, backsolve is pretty important, you need that.
nikosmatsa
commented
3 years ago `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() `plt3(what, what)
gdmcbain
commented
3 years ago 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
commented
3 years ago `NameError Traceback (most recent call last)
nikosmatsa
commented
3 years ago 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
commented
3 years ago 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
commented
3 years ago gdmcbain
commented
3 years ago 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.
nikosmatsa
commented
3 years ago I didn't put it right.I don't want theta rule
gdmcbain
commented
3 years ago 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 / 2or
A = M + L * dt
B = Mrespectively.
From #529 at 2020-12-24: