# linear_elasticity/elastodynamic_identification.py¶

Description

The linear elastodynamics solution of an iron plate impact problem with identification of material parameters from simulated measurement data.

Find such that:

where

## Usage Examples¶

• Run without the identification:

sfepy-run sfepy/examples/linear_elasticity/elastodynamic_identification.py
sfepy-view output/edi/user_block.h5 -f u:wu:f1e3:p0 1:vw:p0

• Get help:

python3 sfepy/examples/linear_elasticity/elastodynamic_identification.py -h

• Run the identification with default parameters, show live plot of convergence and launch ipython shell after the computation:

python3 sfepy/examples/linear_elasticity/elastodynamic_identification.py --plot-log --shell


Result figures are in output/edi, if not changed using –output-dir option.

• Check the Jacobian matrix by finite differences:

python3 sfepy/examples/linear_elasticity/elastodynamic_identification.py --opt-conf=max_nfev=1 --check-jac --shell

• Identify also the damping parameters (zero by default):

python3 sfepy/examples/linear_elasticity/elastodynamic_identification.py --par-names=young,poisson,density,alpha,beta --plot-log --shell


source code

r"""
The linear elastodynamics solution of an iron plate impact problem with
identification of material parameters from simulated measurement data.

Find :math:\ul{u} such that:

.. math::
\int_{\Omega} \rho \ul{v} \pddiff{\ul{u}}{t}
+ \int_{\Omega} D_{ijkl}\ e_{ij}(\ul{v}) e_{kl}(\ul{u})
= 0

where

.. math::
D_{ijkl} = \mu (\delta_{ik} \delta_{jl}+\delta_{il} \delta_{jk}) +
\lambda \ \delta_{ij} \delta_{kl} \;,

\lambda = E \nu / ((1 + \nu)(1 - 2\nu)), \\ \mu = E / 2(1 + \nu)
\;.

Usage Examples
--------------

- Run without the identification::

sfepy-run sfepy/examples/linear_elasticity/elastodynamic_identification.py
sfepy-view output/edi/user_block.h5 -f u:wu:f1e3:p0 1:vw:p0

- Get help::

python3 sfepy/examples/linear_elasticity/elastodynamic_identification.py -h

- Run the identification with default parameters, show live plot of
convergence and launch ipython shell after the computation::

python3 sfepy/examples/linear_elasticity/elastodynamic_identification.py --plot-log --shell

Result figures are in output/edi, if not changed using --output-dir option.

- Check the Jacobian matrix by finite differences::

python3 sfepy/examples/linear_elasticity/elastodynamic_identification.py --opt-conf=max_nfev=1 --check-jac --shell

- Identify also the damping parameters (zero by default)::

python3 sfepy/examples/linear_elasticity/elastodynamic_identification.py --par-names=young,poisson,density,alpha,beta --plot-log --shell

See also :ref:linear_elasticity-elastodynamic.
"""
from argparse import ArgumentParser, RawDescriptionHelpFormatter
import sys
import os.path as op
from functools import partial

import numpy as nm
import scipy.sparse as sps
import matplotlib.pyplot as plt

from sfepy.base.base import output, Struct, IndexedStruct
from sfepy.base.conf import dict_from_string as parse_as_dict
from sfepy.discrete.fem.meshio import UserMeshIO
from sfepy.mesh.mesh_generators import gen_block_mesh
from sfepy.homogenization.utils import define_box_regions
from sfepy.solvers import register_solver, Solver
from sfepy.base.timing import Timer
import sfepy.mechanics.matcoefs as mc

def define(
young=200e9, poisson=0.3, density=7800,
alpha=0.0, beta=0.0,
plane='strain',
dims=(1e-2, 2.5e-3, 2.5e-3),
shape=(21, 6, 6),
v0=1.0,
ct1=1.5,
dt=None,
active_only=False,
save_times=20,
output_dir='output/edi',
**kwargs,
):
"""
Parameters
----------
young, poisson, density: material parameters
plane: plane strain or stress hypothesis
dims: physical dimensions of the block (L, d, x)
shape: numbers of mesh vertices along each axis
v0: initial impact velocity
ct1: final time in L / "longitudinal wave speed" units
dt: time step (None means automatic)
save_times: number of time steps to save
output_dir: output directory
"""
dim = len(dims)

lam, mu = mc.lame_from_youngpoisson(young, poisson, plane=plane)
# Longitudinal and shear wave propagation speeds.
cl = nm.sqrt((lam + 2.0 * mu) / density)
cs = nm.sqrt(mu / density)

# Element size.
L, d = dims[:2]
H = L / (nm.max(shape) - 1)

# Time-stepping parameters.
if dt is None:
# For implicit schemes, dt based on the Courant number C0 = dt * cl / H
# equal to 1.
dt = H / cl # C0 = 1

t1 = ct1 * L / cl

def mesh_hook(mesh, mode):
"""
Generate the block mesh.
"""
mesh = gen_block_mesh(dims, shape, 0.5 * nm.array(dims),
name='user_block', verbose=False)
return mesh

elif mode == 'write':
pass

filename_mesh = UserMeshIO(mesh_hook)

bbox = [[0] * dim, dims]
regions = define_box_regions(dim, bbox[0], bbox[1], 1e-5)

dx = dims[0] / (shape[0] - 1)
sx = int(0.2 * (shape[0] - 1)) * dx
regions.update({
'Omega' : 'all',
'Plane' : ('vertices in (x > %.12f) & (x < %.12f)'
% (0.99999 * sx, 1.00001 * sx),
'facet'),
})
if dim == 3:
regions.update({
'Sensor' : ('r.Plane *v r.Top *v r.Far', 'vertex'),
})

else:
regions.update({
'Sensor' : ('r.Plane *v r.Top', 'vertex'),
})

# Iron.
materials = {
'm' : ({
'young' : young,
'poisson' : poisson,
'density': density,
},),
}

fields = {
'displacement': ('real', 'vector', 'Omega', 1),
}

integrals = {
'i' : 2,
}

# Notes:
# 1. The order of the variables in the solution vector is specified here
#    (3rd tuple member), since that specific order is expected by the
#    elastodynamic time-stepping solvers.
# 2. For the same reason, we won't explicitly define below the equations
#    du = du/dt and ddu = ddu/dt - these are implicitly defined by
#    the time-stepping solver. see the step() method of the solvers.
variables = {
'u' : ('unknown field', 'displacement', 0),
'du' : ('unknown field', 'displacement', 1),
'ddu' : ('unknown field', 'displacement', 2),
'v' : ('test field', 'displacement', 'u'),
'dv' : ('test field', 'displacement', 'du'),
'ddv' : ('test field', 'displacement', 'ddu'),
}
# The mapping of variables for the elastodynamics solvers - keys are given,
# values correspond to the names of the actual variables.
var_names = {'u' : 'u', 'du' : 'du', 'ddu' : 'ddu'}

ebcs = {
'Impact' : ('Left', {'u.0' : 0.0, 'du.0' : 0.0, 'ddu.0' : 0.0}),
}
if dim == 3:
ebcs.update({
'Symmtery-y' : ('Near',
{'u.1' : 0.0, 'du.1' : 0.0, 'ddu.1' : 0.0}),
'Symmetry-z' : ('Bottom',
{'u.2' : 0.0, 'du.2' : 0.0, 'ddu.2' : 0.0}),
})

def get_ic(coor, ic, mode='u'):
val = nm.zeros_like(coor)
if mode == 'u':
val[:, 0] = 0.0

elif mode == 'du':
val[:, 0] = -1.0

return val

functions = {
'get_ic_u' : (get_ic,),
'get_ic_du' : (lambda coor, ic: get_ic(coor, None, mode='du'),),
}

ics = {
'ic' : ('Omega', {'u.all' : 'get_ic_u', 'du.all' : 'get_ic_du'}),
}

equations = {
'eq' :
f"""
+ {alpha} * dw_mass_ad.i.Omega(m.density, dv, du)
+ {beta} * dw_lin_elastic_yp_ad.i.Omega(m.young, m.poisson, dv, du)
+ dw_lin_elastic_yp_ad.i.Omega(m.young, m.poisson, v, u) = 0""",
}

solvers = {
'lsd' : ('ls.auto_direct', {
# 'lsd' : ('ls.scipy_direct', {
# Reuse the factorized linear system from the first time step.
'use_presolve' : True,
# Speed up the above by omitting the matrix digest check used
# normally for verification that the current matrix corresponds to
# the factorized matrix stored in the solver instance. Use with
# care!
'use_mtx_digest' : False,
# Increase when getting MUMPS error -9.
'memory_relaxation' : 50,
}),
'newton' : ('nls.newton', {
'i_max'      : 1,
'eps_a'      : 1e-6,
'eps_r'      : 1e-6,
}),
'tsn' : ('ts.newmark', {
't0' : 0.0,
't1' : t1,
'dt' : dt,
'n_step' : None,

'is_linear'  : True,

'beta' : 0.25,
'gamma' : 0.5,

'var_names' : var_names,
'verbose' : 1,
}),
}

options = {
'ts' : 'tsn',
'tsc' : None,
'nls' : 'newton',
'ls' : 'lsd',

'save_times' : save_times,

'active_only' : active_only,
'auto_transform_equations' : False,

'output_format' : 'h5',
'output_dir' : output_dir,
}

return locals()

from sfepy.solvers.ts_solvers import NewmarkTS, _cache

class NewmarkSATS(NewmarkTS):
r"""
Solve elastodynamics problems by the Newmark method.

The method was introduced in [1]. Common settings [2]:

==================== ======== ==== ===== ==========
name                 kind     beta gamma Omega_crit
==================== ======== ==== ===== ==========
trapezoidal rule:    implicit 1/4  1/2   unconditional
linear acceleration: implicit 1/6  1/2   :math:2\sqrt{3}
Fox-Goodwin:         implicit 1/12 1/2   :math:\sqrt{6}
central difference:  explicit 0    1/2   2
==================== ======== ==== ===== ==========

All of these methods are 2-order of accuracy.

[1] Newmark, N. M. (1959) A method of computation for structural dynamics.
Journal of Engineering Mechanics, ASCE, 85 (EM3) 67-94.

[2] Arnaud Delaplace, David Ryckelynck: Solvers for Computational Mechanics
"""
name = 'ts.newmark_sa'

def create_nlst(self, nls, dt, gamma, beta, u0, e0, v0, a0, pack, unpack):
dt2 = dt**2
iue, iu, ie, iv, ia = pack.indices

if iue != iu:
raise ValueError('extra variables not supported!')

cc0 = (1.0 - gamma) * dt
cc = gamma * dt
ck0 = (0.5 - beta) * dt2
ck = beta * dt2

def v(a):
return v0 + cc0 * a0 + cc * a

def u(a):
return u0 + dt * v0 + ck0 * a0 + ck * a

def fun(at):
vec = nm.r_[u(at), v(at), at]

aux = nls.fun(vec)

rt = aux[iu] + aux[iv] + aux[ia]
return rt

@_cache(self, 'matrix', self.conf.is_linear)
vec = None if self.conf.is_linear else nm.r_[u(at), v(at), at]
M, C, K = self.get_matrices(nls, vec, unpack)

Kt = M + cc * C + ck * K
return Kt

@_cache(self, 'matrix_sa', self.conf.is_linear and (self.ts.step > 0))
vec = None if self.conf.is_linear else nm.r_[u(at), v(at), at]
M, C, K = self.get_matrices(nls, vec, unpack)

size = at.shape[0]
I = sps.eye(size, dtype=nm.float64)
zz = sps.dia_matrix((size, size), dtype=nm.float64)
if self.ts.step > 0:
duda = -ck * I
dvda = -cc * I

else:
duda = dvda = zz
# The initial state might depend on parameters - it does not in
# this example as u0 = 0, v0 = const and a0 = 0, but the
# following code should work in the general case if the next
# line is commented out.
K = C = zz

return sps.bmat([[K, C, M],
[I, zz, duda],
[zz, I, dvda]])

@_cache(self, 'matrix_sa0', self.conf.is_linear and (self.ts.step > 0))
if self.ts.step > 0:
size = at.shape[0]
I = sps.eye(size, dtype=nm.float64)
zz = sps.dia_matrix((size, size), dtype=nm.float64)
duda0 = -ck0 * I
dudv0 = -dt * I
dudu0 = -I
dvda0 = -cc0 * I
dvdv0 = -I
return sps.bmat([[zz, zz, zz],
[dudu0, dudv0, duda0],
[zz, dvdv0, dvda0]])

else:
return sps.csr_matrix((self.sa_info.n_dof, self.sa_info.n_dof),
dtype=nm.float64)

mtx = nm.zeros((self.sa_info.n_dof, self.sa_info.n_par),
dtype=nm.float64, order='F')
mtx = self.context.equations.evaluate(
mode='weak', dw_mode='sensitivity',
diff_vars=self.sa_info.par_names, asm_obj=mtx,
)
for ic, name in enumerate(self.sa_info.par_names):
if name == 'alpha': # = M v
info = self.sa_info.par_info['alpha']

elif name == 'beta': # = K v
info = self.sa_info.par_info['beta']

else:
continue

term = self.context.equations[info[1]].terms[info[2]]
val, iels, status = term.evaluate(mode='weak',
diff_var=None,
standalone=False,
ret_status=True)
val /= term.sign
term.assemble_to(mtx[:, ic], val, iels)

Vt = mtx[iu] + mtx[iv] + mtx[ia]
zz = nm.zeros_like(Vt, dtype=nm.float64)
return nm.block([[Vt],
[zz],
[zz]])

nlst = nls.copy()
nlst.fun = fun
nlst.u = u
nlst.v = v

self.nlst = nlst

return nlst

register_solver(NewmarkSATS)

def apply_sensor(pb, ts, state, out):
us = state['u'].get_state_in_region(pb.domain.regions['Sensor'])
output('sensor', ts.step, us)
out.append(us)

def update_pars(materials, equations, pars, par_names, par_info):
"""
Materials and options are updated in place.
"""
done = set()
for key in par_names:
if key in done: continue

info = par_info[key]
ip = par_names.index(key)
if info[0] == 'term': # Update term coefficient.
equations[info[1]].terms[info[2]].sign = pars[ip]

else: # Update materials.
mat = materials[info[0]]
mkeys = mat.get_keys(info[1])
for mkey in mkeys:
if key == info[2]:
val = pars[ip]

else:
raise ValueError

mat.datas[mkey][info[2]][:, :] = val

output('updated:', par_names)

def eval_fun(pars, data, pb, options, par_names, par_info, opt_data, plog,
inodir):
opt_data.tfun.start()

materials = pb.get_materials()
update_pars(materials, pb.equations, pars, par_names, par_info)

pb.ts.set_step() # Reset ts.
pb.get_solver().clear_lin_solver() # No digest -> clear manually.

out = []
pb.solve(save_results=False, step_hook=partial(apply_sensor, out=out))
usens = nm.concatenate(out)

sim_data = usens[:, -1]
res = sim_data - data

plog(*(tuple(pars) + (nm.linalg.norm(res),)))

times = pb.ts.times
fig, ax = plt.subplots()
ax.plot(times, sim_data, label='sim')
ax.plot(times, data, ls='', marker='o', mew=0, label='exp')
ax.plot(times, opt_data.true_data, ls=':', label='true')
ax.legend()
plt.tight_layout()
fig.savefig(inodir(f'res{opt_data.n_fun:05d}.png'), bbox_inches='tight')

opt_data.n_fun += 1
opt_data.usens = usens

dt = opt_data.tfun.stop()
output(f'function evaluation done in {dt} s.')

return res

def eval_jac_fd(pars, *args, **kwargs):
v0 = eval_fun(pars, *args, **kwargs)

eps = 1e-6
jac = []
for ip, par in enumerate(pars):
pars1 = pars.copy()
dp = eps * par
pars1[ip] += dp

v1 = eval_fun(pars1, *args, **kwargs)
jac.append((v1 - v0) / dp)

jac = nm.array(jac).T

return jac

def eval_jac(pars, data, pb, options, par_names, par_info, opt_data, plog,
inodir):
opt_data.tjac.start()

materials = pb.get_materials()
update_pars(materials, pb.equations, pars, par_names, par_info)

pb.ts.set_step() # Reset ts.
pb.get_solver().clear_lin_solver() # No digest -> clear manually.

variables = pb.get_initial_state()
n_par = len(par_names)

_tss = pb.get_solver()
conf = _tss.conf.copy(name='tss_sa')
conf.kind = 'ts.newmark_sa'

tss = Solver.any_from_conf(conf, nls=_tss.nls, tsc=_tss.tsc, context=pb)
tss.sa_info = Struct(
par_info=par_info, par_names=par_names, n_par=n_par, n_dof=n_dof,
)
tss.matrix_sa = tss.matrix_sa0 = tss.matrix_pars = None

ls = Solver.any_from_conf(tss.nls.lin_solver.conf)

_init_fun, prestep_fun, _poststep_fun = pb.get_tss_functions(
update_bcs=True, update_materials=True,
save_results=False,
step_hook=None, post_process_hook=None,
)

fu = pb.fields['displacement']
ii = fu.get_dofs_in_region(pb.domain.regions['Sensor'])[0]
idof = int(fu.n_components * ii + 2 + variables.adi.indx['u'].start)

def init_fun(tss, ts, vec0):
vec0 = _init_fun(ts, vec0)
# Previous step dy/da.
tss.gy0 = nm.zeros((tss.sa_info.n_dof, tss.sa_info.n_par),
dtype=nm.float64, order='F')
return vec0

jac = []
def poststep_fun(tss, ts, vec):
# EBCs need to be applied here because of algebraically computed
# variables in elastodynamics solvers.
variables.set_state(vec, pb.active_only, apply_ebc=True)

unpack = tss.unpack
utp, etp, vtp, atp = unpack(vec)
if ts.step == 0:
nlst = tss.create_nlst(tss.nls, ts.dt, tss.conf.gamma, tss.conf.beta,
utp, etp, vtp, atp, tss.pack, unpack)

else:
nlst = tss.nlst

ry = -(Vtp + Wtp @ tss.gy0)

if options.multi_rhs:
gy = ls(ry, mtx=Atp)

else:
gy = nm.empty_like(ry)
for ic in range(gy.shape[1]):
gy[:, ic] = ls(ry[:, ic], mtx=Atp)

if ts.step == 0:
# No digest -> clear manually.
ls.clear()

jac.append(gy[idof])

tss.gy0 = gy

return vec

status = IndexedStruct()
tss(variables.get_state(pb.active_only, force=True),
init_fun=init_fun,
prestep_fun=prestep_fun,
poststep_fun=poststep_fun,
status=status)

jac = nm.array(jac)

if options.check_jac:
jac_fd = eval_jac_fd(pars, data, pb, options, par_names, par_info,
opt_data, plog, inodir)
output('Jacobian:')
output(jac)
output('Jacobian using finite differences:')
output(jac_fd)
output('Relative difference:')
output((jac_fd - jac) / nm.where(jac_fd, jac_fd, 1.0))

opt_data.jac = jac.copy()
opt_data.jac_fd = jac_fd

opt_data.n_jac += 1

dt = opt_data.tjac.stop()
output(f'jacobian evaluation done in {dt} s.')

return jac

def parse_args(args=None):
helps = {}

sdefault_opt_conf = """ftol=1e-08, xtol=1e-4, gtol=None, method='trf',
f_scale=3e-8, loss='soft_l1', verbose=2
"""
opts = dict(
par_names = ('density,young,poisson', 'parameters to be identified'),
opt_conf = (sdefault_opt_conf, 'optimization solver options'),
jac = (True, 'do not use the semi-analytical Jacobian calculation'),
check_jac = (False,
'check the Jacobian using finite differences'),
multi_rhs = (False, 'solve all rhs of sensitivity analysis in one call'),
young = (200e9, """Young's modulus"""),
poisson = (0.3, """Poisson's ratio"""),
density = (7800.0, 'density'),
alpha = (0.0, 'proportional damping coefficient (M)'),
beta = (0.0, 'proportional damping coefficient (K)'),
dims = ('1e-2,2.5e-3,2.5e-3',
'physical dimensions of the block (L, d, x)'),
shape = ('21,6,6',
'numbers of mesh vertices along each axis'),
v0 = (1.0, 'initial impact velocity'),
ct1 = (1.5, 'final time in L / "longitudinal wave speed" units'),
dt = (None, 'time step (None means automatic)'),
active_only = (False,
'vectors and matrices contain only active DOFs'),
save_times = (20, 'number of time steps to save'),
output_dir = ('output/edi', 'output directory'),
plot_log = (False, 'show live log figures'),
shell = (False, 'run ipython shell after all computations'),
debug = (False,
'automatically start debugger when an exception is raised'),
)

parser = ArgumentParser(description=__doc__.rstrip(),
formatter_class=RawDescriptionHelpFormatter)

dhelp = ' [default: %(default)s]'
for key, (val, msg) in opts.items():
helps[key] = msg
action = 'store'
vtype = type(val)
choices = None
option = key
if val is True:
action = 'store_false'
option = 'no_' + key

elif val is False:
action = 'store_true'

elif isinstance(val, tuple):
choices = val
vtype = type(val[0])
val = val[0]

elif isinstance(val, list):
vtype = type(val[1])
val = val[0]

if action == 'store':
helps[key] += dhelp
type=vtype,
action=action, dest=key, choices=choices,
default=val, help=helps[key])
else:
action=action, dest=key,
default=val, help=helps[key])
options = parser.parse_args(args=args)

options.dims = [float(ii) for ii in options.dims.split(',')]
options.shape = [int(ii) for ii in options.shape.split(',')]
options.par_names = [ii.strip()  for ii in options.par_names.split(',')]
aux = parse_as_dict(sdefault_opt_conf)
aux.update(parse_as_dict(options.opt_conf))
options.opt_conf = aux

options.plane = 'strain'
options.opt_bounds = {
'density' : [0.5 * options.density, 1.5 * options.density],
'young' : [0.5 * options.young, 1.5 * options.young],
'poisson' : [0.9 * options.poisson, 1.1 * options.poisson],
'alpha' : [0, 1e4],
'beta' : [0, 1e-4],
}

return options, helps

def main():
from scipy.optimize import least_squares

from sfepy.base.base import output, Struct
from sfepy.base.conf import ProblemConf
from sfepy.discrete import Problem
from sfepy.base.log import Log

options, helps = parse_args()
if options.debug:
from sfepy.base.base import debug_on_error; debug_on_error()

conf = ProblemConf.from_dict(define(**vars(options)), sys.modules[__name__])

pb = Problem.from_conf(conf)

out = []
pb.solve(save_results=False, step_hook=partial(apply_sensor, out=out))
usens = nm.concatenate(out)

size = nm.abs(usens[:, -1]).max()
true_data = usens[:, -1]
data = true_data.copy()
i0 = int(0.2 * options.shape[0])
data[i0:] += nm.random.default_rng(seed=42).normal(
0, 0.1 * size, len(data[i0:]),
)
data[max(len(data) // 2, i0)] += 0.5 * size

par_info = {
'density' : ('m', 'Omega', 'density', options.density),
'young' : ('m', 'Omega', 'young', options.young * 0.8),
'poisson' : ('m', 'Omega', 'poisson', options.poisson),
'alpha' : ('term', 'eq', 1, 1000),
'beta' : ('term', 'eq', 2, 1e-5),
}
par_names = options.par_names
all_par_names = set(par_info.keys())
if not set(par_names).issubset(all_par_names):
unknown = set(par_names).difference(all_par_names)
raise ValueError(f'parameters {unknown} not in {all_par_names}!')

labels = [[name] for name in par_names] + [['|res|']]
nlog = len(labels)
formats = [['{:.16e}' for ii in range(len(labs))] for labs in labels]
inodir = partial(op.join, options.output_dir)
plog = Log(labels,
formats=formats,
yscales=['log'] * nlog,
xlabels=['iter'] * nlog,
ylabels=None,
show_legends=True,
is_plot=options.plot_log,
log_filename=inodir('pars.txt'),
aggregate=1000, sleep=0.1)

pars0 = [par_info[name][-1] for name in par_names]
pars0 = nm.array(pars0, dtype=nm.float64)
x_scale = [par if abs(par) > 0 else 1.0 for par in pars0]
bounds = [options.opt_bounds.get(name, (-nm.inf, nm.inf))
for name in par_names]
bounds = tuple(zip(*bounds))
opt_data = Struct(name='opt_data', true_data=true_data, n_fun=0, n_jac=0,
tfun=Timer('fun'), tjac=Timer('jac'))

timer = Timer(start=True)
res = least_squares(
eval_fun, pars0,
jac=eval_jac if options.jac else '2-point',
bounds=bounds,
x_scale=x_scale,
args=(data, pb, options, par_names, par_info, opt_data, plog, inodir),
**options.opt_conf
)
elapsed = timer.stop()
output(f'identification done in {elapsed} s.')
output(res)
od = opt_data
tfi = od.tfun.total / od.n_fun
tji = od.tjac.total / max(od.n_jac, 1)
tfe = tfi * len(pars0)
output(f'evaluations: function: {od.n_fun}, jacobian: {od.n_jac}')
output(f'elapsed: function: {od.tfun.total}, jacobian {od.tjac.total}')
output(f'elapsed/it: function: {tfi}, jacobian {tji}')
output(f'expected finite difference jacobian time: {tfe}')

plog(save_figure=inodir('pars.png'))

if options.shell:
from sfepy.base.base import shell; shell()

plog(finished=True)

if __name__ == '__main__':
main()