#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------
#
# _____ _ _ _ _____ _____
# | ____|__| | ___| |_ _____(_)___ ___| ___| ____|
# | _| / _` |/ _ \ \ \ /\ / / _ \ / __/ __| |_ | _|
# | |__| (_| | __/ |\ V V / __/ \__ \__ \ _| | |___
# |_____\__,_|\___|_| \_/\_/ \___|_|___/___/_| |_____|
#
#
# Unit of Strength of Materials and Structural Analysis
# University of Innsbruck,
# 2017 - today
#
# Matthias Neuner matthias.neuner@uibk.ac.at
#
# This file is part of EdelweissFE.
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
# License as published by the Free Software Foundation; either
# version 2.1 of the License, or (at your option) any later version.
#
# The full text of the license can be found in the file LICENSE.md at
# the top level directory of EdelweissFE.
# ---------------------------------------------------------------------
"""The nonlinear implicit dynamic solver: Newmark-beta time integration on top of the Newton loop
of :class:`~edelweissfe.solvers.nonlinearimplicitstatic.NIST`.
An increment solves the semi-discrete equation of motion
.. math::
\\boldsymbol{M} \\ddot{\\boldsymbol{u}}_{n+1} + \\boldsymbol{C} \\dot{\\boldsymbol{u}}_{n+1}
+ \\boldsymbol{P}_\\mathrm{int}(\\boldsymbol{u}_{n+1}) = \\boldsymbol{P}_\\mathrm{ext}(t_{n+1})
for the displacement at the end of the increment, with the velocity and the acceleration expressed
through the displacement increment :math:`\\Delta \\boldsymbol{u} = \\boldsymbol{u}_{n+1} - \\boldsymbol{u}_n`
by the two Newmark relations (Newmark 1959; Hughes, *The Finite Element Method*, ch. 9)
.. math::
\\boldsymbol{u}_{n+1} &= \\boldsymbol{u}_n + \\Delta t \\, \\dot{\\boldsymbol{u}}_n
+ \\tfrac{\\Delta t^2}{2} \\left[ (1 - 2\\beta) \\, \\ddot{\\boldsymbol{u}}_n
+ 2\\beta \\, \\ddot{\\boldsymbol{u}}_{n+1} \\right], \\\\
\\dot{\\boldsymbol{u}}_{n+1} &= \\dot{\\boldsymbol{u}}_n
+ \\Delta t \\left[ (1 - \\gamma) \\, \\ddot{\\boldsymbol{u}}_n + \\gamma \\, \\ddot{\\boldsymbol{u}}_{n+1} \\right].
Solved for the end-of-increment kinematics these read
.. math::
\\ddot{\\boldsymbol{u}}_{n+1} &= \\frac{\\Delta \\boldsymbol{u} - \\Delta t \\, \\dot{\\boldsymbol{u}}_n
- \\Delta t^2 \\left( \\tfrac{1}{2} - \\beta \\right) \\ddot{\\boldsymbol{u}}_n}{\\beta \\, \\Delta t^2}, \\\\
\\dot{\\boldsymbol{u}}_{n+1} &= \\dot{\\boldsymbol{u}}_n
+ \\Delta t \\left[ (1 - \\gamma) \\, \\ddot{\\boldsymbol{u}}_n + \\gamma \\, \\ddot{\\boldsymbol{u}}_{n+1} \\right],
so that the effective residual and tangent handed to the Newton loop are
.. math::
\\boldsymbol{R} &= \\boldsymbol{f}_\\mathrm{ext} - \\boldsymbol{f}_\\mathrm{int}
- \\boldsymbol{M} \\ddot{\\boldsymbol{u}}_{n+1} - \\boldsymbol{C} \\dot{\\boldsymbol{u}}_{n+1}, \\\\
\\boldsymbol{K}_\\mathrm{eff} &= \\boldsymbol{K} + \\frac{1}{\\beta \\, \\Delta t^2} \\boldsymbol{M}
+ \\frac{\\gamma}{\\beta \\, \\Delta t} \\boldsymbol{C}.
Everything else -- element evaluation, loads, constraints, multi-point-constraint condensation,
Dirichlet handling, the convergence test, cutbacks and the linear solver -- is the parent's: this
solver runs :meth:`~edelweissfe.solvers.nonlinearimplicitstatic.NIST.solveIncrement` itself and
enters it only through its three hooks -- :meth:`NonlinearImplicitDynamic.initializeIncrement`
(mass, damping, dynamic stiffness, initial acceleration), :meth:`NonlinearImplicitDynamic.assembleAdditionalTerms`
(two terms in the residual and two in the tangent, where the parent forms ``R = f_ext - f_int``
and before the tangent is converted to CSR) and :meth:`NonlinearImplicitDynamic.finalizeIncrement`
(the converged velocity and acceleration become the state).
**Parameters.** The default :math:`\\beta = 1/4`, :math:`\\gamma = 1/2` is the average-acceleration
(trapezoidal) rule: second-order accurate, unconditionally stable for linear problems and free of
numerical damping, so the discrete energy of an undamped linear system is conserved exactly and
any drift in it points at a bug rather than at the scheme. :math:`\\gamma > 1/2` adds algorithmic
damping at the price of first-order accuracy; :math:`\\gamma \\geq 1/2` together with
:math:`\\beta \\geq (\\gamma + 1/2)^2 / 4` keeps the scheme unconditionally stable. A choice outside
that region is accepted with a warning, :math:`\\beta = 0` (the explicit central-difference member
of the family) is refused, because the formulation above divides by it.
**Mass.** The mass is the CONSISTENT one, :math:`\\int_\\Omega \\rho N^T N \\, dV`, assembled from
:meth:`~edelweissfe.elements.base.baseelement.BaseElement.computeConsistentInertia` into the very
sparsity pattern the stiffness uses -- both are scattered through the same VIJ layout of the
:class:`~edelweissfe.numerics.dofmanager.DofManager`, so adding :math:`\\boldsymbol{M}/(\\beta \\Delta t^2)`
to the tangent is an entry-wise scaled addition of two value vectors, before the parent's in-place CSR update, its
multi-point-constraint condensation and its Dirichlet row replacement, all of which therefore act
on the effective matrix without knowing it is one. For the same reason the mass is summed into
CSR form by the stiffness' own CSR generator. The damping :math:`\\boldsymbol{C}` is the diagonal
each element reports through ``computeLumpedDamping``; being diagonal, it is kept only as its
nonzero entries and their positions in the same layout, not as a full value vector. Both are
assembled in one loop over the elements (in chunks of elements on the thread pool, see
:doc:`/documentation/parallelization`) when the equation system is (re)built -- at a step's start
and after a topology change -- not per Newton iteration. A rebuild caused only by a constraint changing
its connectivity (a contact candidate list, which can change on every increment) leaves every
element, and with it every element's slot in the layout, where it was: the operators are then
reused in the new layout rather than reassembled -- see
:meth:`NonlinearImplicitDynamic._reuseMassAndDamping`. The one exception is a step that changes a
material property mid-step (``>>changematerialproperty``), which invalidates the density and makes
them reassemble every increment.
**Checks on the mass.** The assembled mass and damping are checked at runtime for non-finite
entries, the damping for negative entries, and every time-integrated degree of freedom for a
nonzero mass; the total mass of each field is reported. The *symmetry* of the mass is not checked
at runtime -- on the assembled matrix that costs a transposed copy of it, gigabytes on a large
model -- but pinned once per element type in ``tests/test_nid_element_mass.py``, together with the
total mass :math:`\\rho V`. It does matter: a non-symmetric mass still converges, to a wrong response.
That test covers only the element types it lists -- the Marmot hexahedra ``C3D8``, ``C3D20``,
``C3D20R`` and the gradient-enhanced ``GC3D8``, ``GC3D20R``; any other element type (the plane
elements, for one) is *not* guarded anywhere, so add it to that list before using it with this
solver.
**Which fields.** Only the fields whose inertia is a mass
(:func:`~edelweissfe.config.phenomena.carriesLinearMomentum`, i.e. the displacement) are integrated
in time. Rows and columns of the mass and the damping belonging to any other field are zeroed, and
those fields keep the parent's quasi-static treatment. A nonzero inertia an element reports there --
a rotational inertia, or the micro-inertia of a gradient-enhanced element -- is discarded with a
warning, once per field. A model with no such field is refused.
**Zero-length and negligible increments.** The time stepper proposes a zero-length increment before
the first real one of every step. A quasi-static solver equilibrates it; this solver skips it with
the state kept, because in zero time nothing can move: displacement and velocity are continuous,
and a load appearing at that instant is answered by the acceleration below, not by a displacement.
Solving it statically would put a suddenly loaded model at its static deflection before the first
real increment -- the whole dynamic response skipped. An increment shorter than :math:`10^{-10}` of
the time elapsed in the step is skipped for the same reason: it is the round-off remainder the
stepper's progress accumulation leaves at the end of a step (400 increments of 0.0025 leave 2e-14),
and the Newmark update would divide by its square. A skipped increment evaluates no force, so it
leaves the node field result ``P`` as the last accepted increment wrote it.
**Initial acceleration.** At the start of a step the acceleration is computed from equilibrium,
:math:`\\boldsymbol{M} \\ddot{\\boldsymbol{u}}_0 = \\boldsymbol{P}_\\mathrm{ext}(t_0) - \\boldsymbol{P}_\\mathrm{int}(\\boldsymbol{u}_0)
- \\boldsymbol{C} \\dot{\\boldsymbol{u}}_0`, with the loads evaluated at the step's own start (an
``f(t)`` amplitude at :math:`t = 0`). Without it, a load applied suddenly at the start of a step
would enter the trapezoidal rule as though it were zero at :math:`t_0`, which costs one first-order
error in the momentum, and the scheme's second order is lost. It is computed in the first real
increment of the step. A step that continues from a restart checkpoint -- a resumed run, or a step
that falls back to one -- therefore computes it only if the checkpoint was written before that,
after the zero increment: every later checkpoint carries the acceleration, and recomputing it would
replace a consistent state with one that agrees only to solver tolerance. Switch it off with
``computeInitialAcceleration=False`` to continue a step from whatever acceleration was carried over.
**State and restart.** Velocity and acceleration are ordinary node-field entries, ``V`` and ``A``,
read from the model whenever the equation system is (re)built and written at the acceptance of every
increment alongside the parent's ``U``/``P``/``dU``, both through the model's layout of the equation
system -- the same mechanism the explicit solver uses for its velocity. That makes them part of every checkpoint the
``*output, type=restart`` manager writes and of every ``*fieldOutput`` (``result=V``,
``result=A``), with no solver-specific restart code beyond the marker of whether the initial
acceleration is still to be computed. The force written as ``P`` is the internal force with the
inertial and damping forces of the converged state, :math:`\\boldsymbol{f}_\\mathrm{int} +
\\boldsymbol{M} \\ddot{\\boldsymbol{u}}_{n+1} + \\boldsymbol{C} \\dot{\\boldsymbol{u}}_{n+1}`, as the
documentation page :doc:`/documentation/nodefieldresults` defines it: at a support without load, the
reaction.
**Live h-adaptivity.** A ``*modelModifier, type=hAdaptivity`` refining the mesh mid-step is
supported. Three things make that work, and only the third is this solver's own:
* Being node-field entries, ``V`` and ``A`` are carried onto the nodes a refinement creates by the
same isoparametric interpolation from the parent element that already carries ``U``
(:data:`~edelweissfe.modelmodifiers.adaptivity.hadaptivity.WARM_STARTED_NODE_FIELD_ENTRIES`).
Nothing here interpolates anything itself.
* The mass and the damping are not transferred at all: they are reassembled from the elements that
now exist, because the parent rebuilds its
:class:`~edelweissfe.numerics.dofmanager.DofManager` on a topology change and everything bound to
it is rebuilt with it. There is therefore no re-lumping error of the kind an explicit solver has
to account for -- the total mass is an assembly of the new mesh, not a redistribution of the old.
* The interpolated acceleration is only as good as an interpolation of the *previous* mesh's
acceleration, and it is not in equilibrium with the operators that were just reassembled. So the
increment after a topology change **recomputes the acceleration from equilibrium**, with the
cold-start machinery a step start already uses (``computeInitialAcceleration``, which also
switches this off). The displacement and the velocity are kept: they are the kinematic state, and
interpolation is the best statement of them the new mesh admits. Only a change that the topology
pipeline actually recorded re-arms it -- a mere change of a constraint's connectivity (a contact
candidate list) rebuilds the equation system without moving a node, and must not restart the
acceleration of every degree of freedom in the model.
What is exactly conserved across a refinement, and what is not, is reported at every event and
follows from the same partition-of-unity argument the explicit solver's own live-adaptivity notes
make: the assembled total mass is a geometric identity (children tile their parent at the same
density) and is *asserted*; the linear momentum and the kinetic energy are exact for a spatially
uniform velocity field and differ at second order in the velocity gradient across a refined parent
otherwise, which is discretisation error rather than a defect, and are reported only. Unlike an
explicit run, the increment that follows re-establishes equilibrium by Newton iteration, so an
imbalance the transfer leaves behind is absorbed rather than radiated.
**Prescribed motion.** On a degree of freedom with a Dirichlet condition the acceleration is the
one the Newmark relations imply for the prescribed displacement history, starting from zero at the
step's start. For a prescribed displacement that is not smooth in time -- a velocity ramp switched
on suddenly -- the average-acceleration rule answers with an acceleration that oscillates from
increment to increment on those degrees of freedom, and the consistent mass, having off-diagonal
entries, passes the corresponding inertia forces on to the free neighbours. A smooth amplitude, or
some algorithmic damping (``newmarkGamma > 1/2``), keeps that out of the response.
**Not supported in this version.** Rayleigh or any other damping model beyond what the elements
report as their lumped damping; elements without a consistent mass (the pure-Python element
library), which are refused at the first increment.
"""
from dataclasses import dataclass
import numpy as np
from scipy.sparse import coo_matrix, csr_matrix, diags
import edelweissfe.utils.performancetiming as performancetiming
from edelweissfe.config.phenomena import carriesLinearMomentum
from edelweissfe.models.femodel import FEModel
from edelweissfe.numerics.dofmanager import DofManager, DofVector, VIJSystemMatrix
from edelweissfe.numerics.parallelizationutilities import (
forEachSliceInParallel,
mapChunksInParallel,
)
from edelweissfe.numerics.systemmatrixassembly import (
SystemMatrixAssembly,
VIJSystemMatrixAssembly,
)
from edelweissfe.outputmanagers.base.outputmanagerbase import OutputManagerBase
from edelweissfe.solvers.base.conservationchecks import (
CONSERVATION_TOLERANCE,
ConservationCheck,
formatMomentumAndKineticEnergy,
linearMomentum,
)
from edelweissfe.solvers.nonlinearimplicitstatic import NIST, NISTSchema
from edelweissfe.timesteppers.timestep import TimeStep
from edelweissfe.utils.fieldoutput import FieldOutputController
from edelweissfe.utils.schema import schemaField
#: An increment shorter than this fraction of the elapsed step time is the round-off remainder of the
#: time stepper's progress accumulation and is skipped with the state kept; see
#: :meth:`NonlinearImplicitDynamic.solveIncrement`.
_NEGLIGIBLE_INCREMENT_FRACTION = 1e-10
#: Rows per block of the mass matrix, in which the inertia forces are computed on the thread pool (see
#: :meth:`NonlinearImplicitDynamic._rowBlocks`). A smaller model is a single block, computed serially.
_MASS_ROWS_PER_BLOCK = 1 << 14
[docs]
@dataclass(frozen=True)
class NIDSchema(NISTSchema):
""":class:`~edelweissfe.solvers.nonlinearimplicitstatic.NISTSchema`'s options, plus the
Newmark parameters and the initial-acceleration switch. Mirrors
:attr:`NonlinearImplicitDynamic.SolverSpecificOptions` one-for-one; that dict remains the
runtime source of truth, for the reasons given on the parent schema.
"""
newmarkBeta: float | None = schemaField(
description=(
"Newmark parameter beta, weighting the end-of-increment acceleration in the displacement "
"update. Must be > 0. 1/4 with gamma = 1/2 is the average-acceleration (trapezoidal) rule."
),
dtype=float,
default=0.25,
)
newmarkGamma: float | None = schemaField(
description=(
"Newmark parameter gamma, weighting the end-of-increment acceleration in the velocity "
"update. 1/2 gives second-order accuracy and no algorithmic damping; larger values damp."
),
dtype=float,
default=0.5,
)
computeInitialAcceleration: bool | None = schemaField(
description=(
"Compute the acceleration at the start of every step from equilibrium with the loads at "
"the step's start, instead of continuing from the carried-over acceleration. Skipped "
"automatically when the step resumes from a restart checkpoint. Also re-armed for the "
"increment following a topology change (h-adaptivity), whose interpolated acceleration "
"is not in equilibrium with the operators reassembled on the refined mesh."
),
dtype=bool,
default=True,
)
massConservationTolerance: float | None = schemaField(
description=(
"Relative-change tolerance for the per-topology-change check on a dynamic field's total "
"assembled mass. Children of a refined element tile it and carry the same density, so "
"this is conserved geometrically; the default is already generous relative to floating-"
"point precision and exists to catch a genuinely wrong refinement, not to absorb "
"ordinary quadrature noise. The counterpart of the explicit solver's "
"lumped-quantity-conservation-tolerance."
),
dtype=float,
default=CONSERVATION_TOLERANCE,
)
@dataclass
class _NewmarkSystem:
"""Everything the Newmark increment reads that is sized by the current equation system.
Bound to the :class:`~edelweissfe.numerics.dofmanager.DofManager` it was assembled for: the
moment the parent rebuilds that manager, every entry here is stale together, so the whole
bundle is replaced rather than any part of it patched.
Parameters
----------
dofManager
The manager these operators and vectors are indexed by.
dynamicDofs
The degrees of freedom integrated in time -- those of the fields carrying a mass.
dynamicFields
The names of those fields, in the model's field order.
Mvij
The consistent mass as a VIJ value vector, in the stiffness' layout, ready to be scaled
and added onto the tangent.
dampingVIJIndices, dampingVIJValues
The diagonal damping, only where it is nonzero: each element's lumped damping at its own
diagonal slots of the VIJ layout. Diagonal by construction, so a full VIJ value vector for
it would be all zeros but for a few entries per element.
M, C
The mass and the damping as CSR matrices, for the residual's matrix-vector products.
MRowBlocks, MRowsPerBlock
``M`` in blocks of ``MRowsPerBlock`` consecutive rows (:data:`_MASS_ROWS_PER_BLOCK`), sharing
its arrays, for the products with the mass on the thread pool (:meth:`massTimes`).
massByField
Each dynamic field's total assembled mass (:meth:`NonlinearImplicitDynamic._totalMassByField`),
a property of ``M`` alone, so computed once here rather than per accepted increment.
V, A
The velocity and acceleration at the last converged increment: the solver's own vectors of the
current equation system, read from the model when it was built.
"""
dofManager: DofManager
dynamicDofs: np.ndarray
dynamicFields: list
Mvij: VIJSystemMatrix
dampingVIJIndices: np.ndarray
dampingVIJValues: np.ndarray
M: csr_matrix
C: csr_matrix
MRowBlocks: list[csr_matrix]
MRowsPerBlock: int
massByField: dict
V: DofVector
A: DofVector
def massTimes(self, vector) -> np.ndarray:
"""``M @ vector``, one row block of :attr:`MRowBlocks` per call on the thread pool -- bitwise
the product with the whole of ``M`` (see :meth:`NonlinearImplicitDynamic._rowBlocks`)."""
vector = np.asarray(vector)
product = np.empty_like(vector)
def multiplyRows(rows: slice):
product[rows] = self.MRowBlocks[rows.start // self.MRowsPerBlock] @ vector
forEachSliceInParallel(multiplyRows, len(product), self.MRowsPerBlock)
return product
@dataclass
class _NewmarkIncrement:
"""What the Newton iterations of one increment share: the system, the increment's Newmark
coefficients and tangent terms, and the trial kinematics of the current iteration.
Parameters
----------
system
The Newmark system of the increment.
dT
The time increment.
beta, gamma
The Newmark parameters.
massFactor, dampingFactor
:math:`1 / (\\beta \\Delta t^2)` and :math:`\\gamma / (\\beta \\Delta t)`: the derivatives of the
trial acceleration and velocity with respect to the displacement increment, i.e. the factors
the mass matrix and the damping enter the effective tangent with.
V_np, A_np
The trial velocity and acceleration of the current iteration.
"""
system: _NewmarkSystem
dT: float
beta: float
gamma: float
massFactor: float
dampingFactor: float
V_np: DofVector
A_np: DofVector
@dataclass(frozen=True)
class _ConservedQuantities:
"""What a topology change ought to leave alone, measured on one :class:`_NewmarkSystem`.
Parameters
----------
massByField
Each dynamic field's own total assembled mass. Per field, never summed across them: a
violation in a numerically small field would otherwise hide inside a large one.
momentum
The linear momentum :math:`M v`, per spatial component. Per component, not summed over
them: adding a momentum's x, y and z contributions produces a number with no physical
meaning and would hide a component-wise error behind a cancellation.
kineticEnergy
:math:`\\tfrac{1}{2} v^T M v`.
"""
massByField: dict
momentum: np.ndarray
kineticEnergy: float
#: Restored with the solver that measured it; see :attr:`NonlinearImplicitDynamic.checkpointedState`.
checkpointedState = {"massByField": dict, "momentum": np.ndarray, "kineticEnergy": float}
[docs]
class NonlinearImplicitDynamic(NIST):
"""This is the Nonlinear Implicit Dynamic -- solver (``NID``), Newmark-beta time integration on
top of the Newton loop of :class:`~edelweissfe.solvers.nonlinearimplicitstatic.NIST`.
Parameters
----------
jobInfo
A dictionary containing the job information.
journal
The journal instance for logging.
"""
identification = "NID"
#: Live h-adaptivity is supported; see the module docstring.
supportsModelModifiers = True
#: Option schema for this solver, per OptionSchemaProvider.
schema = NIDSchema
#: Only the VIJ assembly: the effective tangent adds the mass and the damping entry by entry to the VIJ value array
#: (see :meth:`assembleAdditionalTerms`), and the mass is gathered into CSR by the generator's gather map.
systemMatrixAssemblies = {"vij": VIJSystemMatrixAssembly}
#: The predictor's state, as in the parent, the drift the conservation checks accumulate over the
#: step, the conserved totals after the last accepted increment, which a topology change is
#: compared against, and whether the next increment computes the acceleration from equilibrium
#: first, which a checkpoint written after the zero increment of a step still has to do. The
#: velocity and acceleration are node-field entries, restored with the model.
checkpointedState = NIST.checkpointedState | {
"_conservationCheck": ConservationCheck,
"_conservedAtLastAccept": _ConservedQuantities,
"_initialAccelerationPending": bool,
}
SolverSpecificOptions = NIST.SolverSpecificOptions | {
"newmarkBeta": 0.25,
"newmarkGamma": 0.5,
"computeInitialAcceleration": True,
"massConservationTolerance": CONSERVATION_TOLERANCE,
}
def __init__(self, jobInfo, journal, **kwargs):
super().__init__(jobInfo, journal, **kwargs)
#: The operators and kinematic state of the current equation system; None until the
#: first increment of a step has built it. See :class:`_NewmarkSystem`.
self._newmarkSystem = None
#: Whether the next increment has to compute the acceleration from equilibrium first.
self._initialAccelerationPending = False
#: Length of ``model.topology.history`` when the current equation system was assembled; tells
#: a mesh change from a constraint re-reporting its connectivity (see :meth:`_updateNewmarkSystem`).
self._topologyRecordsAtLastBuild = 0
#: The per-topology-change mass check and the drift it accumulates over a step; see
#: :class:`~edelweissfe.solvers.base.conservationchecks.ConservationCheck`.
self._conservationCheck = ConservationCheck(journal, self.identification)
#: The state shared by the Newton iterations of the current increment; see
#: :class:`_NewmarkIncrement`. Set by :meth:`initializeIncrement`.
self._currentIncrement = None
#: The non-time-integrated fields a discarded inertia has already been warned about.
self._fieldsWarnedAboutDiscardedInertia = set()
#: The velocity and the acceleration of the current equation system, read from the model when it is
#: (re)built (see :meth:`readStateFromModel`), and the state of the last converged increment from then on.
self._V = None
self._A = None
#: The names of the fields of the current equation system that carry a mass, whose velocity and acceleration
#: this solver reads from and writes to the model.
self._dynamicFields = None
self._validateNewmarkParameters()
def _validateNewmarkParameters(self):
"""Refuse a beta the implicit formulation divides by zero with, and warn about a pair
outside the unconditionally stable region."""
beta = self.options["newmarkBeta"]
gamma = self.options["newmarkGamma"]
if beta <= 0.0:
raise ValueError(
"newmarkBeta must be > 0 (got {:}): the implicit Newmark update divides by beta*dt^2; "
"beta = 0 is the explicit central-difference member of the family, which this solver "
"does not integrate -- use NED for that.".format(beta)
)
if gamma < 0.5 or beta < 0.25 * (gamma + 0.5) ** 2:
self.journal.message(
"Newmark parameters beta={:}, gamma={:} lie outside the unconditionally stable region "
"(gamma >= 1/2 and beta >= (gamma + 1/2)^2 / 4); the time increment is then bounded by "
"the highest frequency of the mesh, which nothing here checks.".format(beta, gamma),
self.identification,
0,
)
elif gamma != 0.5:
self.journal.message(
"Newmark gamma={:} differs from 1/2: the scheme is numerically damped and first-order "
"accurate in time.".format(gamma),
self.identification,
1,
)
[docs]
def beginStep(
self,
step,
model: FEModel,
fieldOutputController: FieldOutputController,
outputmanagers: dict[str, OutputManagerBase],
):
"""Start a step; see the parent.
Arms the initial-acceleration computation and drops the operators of the previous step -- the
parent rebuilds its equation system at the start of every step, and the operators follow that
system. A step that continues from a restart checkpoint takes over, with the checkpointed state
restored after this, whether the initial acceleration is still to be computed: the first real
increment of the step computes it, and a checkpoint written after the zero increment comes
before that. The velocity and acceleration are not dropped: they
are re-read from the ``V``/``A`` node-field entries the previous step (or the checkpoint)
left, which is what carries the kinematic state across steps.
Parameters
----------
step
The step to be solved.
model
The model tree.
fieldOutputController
The field output controller.
outputmanagers
The output managers.
"""
# >>options blocks may have changed them since construction.
self._validateNewmarkParameters()
self._newmarkSystem = None
self._conservationCheck.reset()
#: Mass, momentum and kinetic energy after the last accepted increment of this step, or None.
self._conservedAtLastAccept = None
# The mesh the step starts on, restored or not: a topology change is a change from this.
self._topologyRecordsAtLastBuild = len(model.topology.history)
resumed = not step.timeStepper.isAtStepStart()
# A resumed step overwrites it with the checkpointed one.
self._initialAccelerationPending = bool(self.options["computeInitialAcceleration"])
self.journal.message(
"Newmark-beta time integration, beta={:}, gamma={:}; consistent mass{:}".format(
self.options["newmarkBeta"],
self.options["newmarkGamma"],
(
"; continuing from the restart checkpoint"
if resumed
else (
"; initial acceleration from equilibrium"
if self._initialAccelerationPending
else "; carrying the acceleration over"
)
),
),
self.identification,
1,
)
return super().beginStep(step, model, fieldOutputController, outputmanagers)
[docs]
def solveIncrement(
self,
U_n: DofVector,
dU: DofVector,
P: DofVector,
K: SystemMatrixAssembly,
stepActions: list,
model: FEModel,
timeStep: TimeStep,
prevTimeStep: TimeStep,
extrapolation: str,
maxIter: int,
maxGrowingIter: int,
) -> tuple[DofVector, DofVector, DofVector, int, dict]:
"""Newton-Raphson scheme on the Newmark-effective equation of motion of an increment.
The parent's loop, into which the dynamics enter through its hooks
(:meth:`initializeIncrement`, :meth:`assembleAdditionalTerms`, :meth:`finalizeIncrement`). A
zero-length increment -- the one the time stepper
yields before the first real increment of a step -- has no equation of motion and is
skipped with the state kept: displacement and velocity are continuous in time, and a load
appearing at that instant is answered by the initial acceleration, not by a displacement.
Parameters
----------
U_n
The old solution vector.
dU
The old solution increment.
P
The old reaction vector.
K
The system matrix to be used.
stepActions
The list of active step actions.
model
The model tree.
timeStep
The time step.
prevTimeStep
The previous time step.
extrapolation
The type of extrapolation to be used.
maxIter
The maximum number of iterations to be used.
maxGrowingIter
The maximum number of growing residuals until the Newton-Raphson is terminated.
Returns
-------
tuple[DofVector,DofVector,DofVector,int,dict]
A tuple containing
- the new solution vector
- the solution increment
- the new reaction vector
- the number of required iterations
- the history of residuals per field
"""
dT = timeStep.timeIncrement
if dT <= _NEGLIGIBLE_INCREMENT_FRACTION * abs(timeStep.stepTime):
# A zero-length increment has no equation of motion: displacement and velocity cannot change
# in zero time, and a load appearing now appears as an acceleration (the initial acceleration
# solve), not as a static displacement. The time stepper proposes one before every step.
# A negligible round-off increment is skipped too: the Newmark update divides by beta*dt^2
# and would amplify the round-off of the linear solve into the acceleration and velocity.
self.applyStepActionsAtIncrementStart(model, timeStep, stepActions)
self.journal.message(
"{:} increment ({:e}): no equation of motion to integrate, state kept".format(
"zero-length" if dT <= 0.0 else "negligible", dT
),
self.identification,
2,
)
dU[:] = 0.0
# No Newmark increment, so no force of this state is evaluated: the model keeps the force it holds (see
# writeResultsToModel).
self._currentIncrement = None
return U_n, dU, P, 0, {field: [] for field in self.theDofManager.idcsOfFieldsInDofVector}
return super().solveIncrement(
U_n, dU, P, K, stepActions, model, timeStep, prevTimeStep, extrapolation, maxIter, maxGrowingIter
)
[docs]
def initializeIncrement(self, U_n: DofVector, stepActions: dict, model: FEModel, timeStep: TimeStep):
"""Bring the Newmark operators up to date with the equation system, form the two terms the
dynamics add to the tangent, and -- at a cold step start or after a topology change --
compute the acceleration from equilibrium; see the parent's hook.
After the step actions, which the parent has applied just before: a
``>>changematerialproperty`` acting on this increment changes the density the mass is
assembled from.
Parameters
----------
U_n
The solution at the start of the increment.
stepActions
The active step actions.
model
The model tree.
timeStep
The time step.
"""
system = self._updateNewmarkSystem(model, stepActions)
beta = self.options["newmarkBeta"]
gamma = self.options["newmarkGamma"]
dT = timeStep.timeIncrement
self._currentIncrement = _NewmarkIncrement(
system=system,
dT=dT,
beta=beta,
gamma=gamma,
massFactor=1.0 / (beta * dT * dT),
dampingFactor=gamma / (beta * dT),
V_np=self.theDofManager.constructDofVector(),
A_np=self.theDofManager.constructDofVector(),
)
if self._initialAccelerationPending:
self._initialAccelerationPending = False
self._computeInitialAcceleration(system, U_n, stepActions, model, timeStep)
[docs]
@performancetiming.timeit("inertia and damping terms")
def assembleAdditionalTerms(
self, dU: DofVector, R: DofVector, F: DofVector, K: SystemMatrixAssembly, timeStep: TimeStep
):
"""Add the inertia and damping forces to the residual and their derivatives to the
tangent; see the parent's hook. The trial displacement increment fixes the trial kinematics.
Parameters
----------
dU
The current trial solution increment.
R
The residual, augmented in place.
F
The reference flux scale of the convergence test, augmented in place.
K
The tangent being assembled; its VIJ value array is augmented in place.
timeStep
The time step.
"""
increment = self._currentIncrement
system = increment.system
with performancetiming.timeit("newmark kinematics"):
self._newmarkVelocityAndAcceleration(
dU,
system.V,
system.A,
increment.dT,
increment.beta,
increment.gamma,
system.dynamicDofs,
increment.V_np,
increment.A_np,
)
with performancetiming.timeit("inertia and damping forces"):
# The damping is diagonal, so C @ V_np costs next to nothing.
PInertia = system.massTimes(increment.A_np)
PDamping = system.C @ np.asarray(increment.V_np)
R -= PInertia
R -= PDamping
# Inertia and damping forces enter the reference scale of the relative flux tolerance too
# (a body in free flight has no internal force at all).
F += np.abs(PInertia)
F += np.abs(PDamping)
# The effective tangent K + M / (beta dT^2) + gamma C / (beta dT). Same VIJ layout as K, so it
# is an entry-wise sum, before the parent's CSR conversion, MPC condensation and Dirichlet row
# replacement. The mass is added slice by slice on the thread pool; the damping is diagonal.
with performancetiming.timeit("effective tangent"):
Kvij = K.vijArray
# A slice-wise sum would silently add only a prefix of a longer mass vector.
if len(system.Mvij) != len(Kvij):
raise ValueError(
"The mass ({:} entries) and the tangent ({:} entries) differ in their VIJ layout.".format(
len(system.Mvij), len(Kvij)
)
)
def addScaledMass(s: slice):
Kvij[s] += increment.massFactor * system.Mvij[s]
forEachSliceInParallel(addScaledMass, len(Kvij))
Kvij[system.dampingVIJIndices] += increment.dampingFactor * system.dampingVIJValues
[docs]
@performancetiming.timeit("finalize increment")
def finalizeIncrement(self, model: FEModel):
"""The converged trial kinematics become the state of the solver, which the model receives at
the acceptance of the increment (see :meth:`writeResultsToModel`); see the parent's hook. Only a
converged increment gets here -- a failing one leaves by exception and keeps the previous
velocity and acceleration for the cutback.
Parameters
----------
model
The model tree.
"""
increment = self._currentIncrement
system = increment.system
system.V[:] = increment.V_np
system.A[:] = increment.A_np
# What a topology change before the next increment is compared against.
self._conservedAtLastAccept = self._conservedQuantities(system, model)
self.journal.message(
"kinetic energy {:e}".format(self._conservedAtLastAccept.kineticEnergy), self.identification, 2
)
@staticmethod
def _newmarkVelocityAndAcceleration(
dU: DofVector,
V_n: DofVector,
A_n: DofVector,
dT: float,
beta: float,
gamma: float,
dynamicDofs: np.ndarray,
V_np: DofVector,
A_np: DofVector,
):
"""The end-of-increment velocity and acceleration implied by a displacement increment.
The two Newmark relations solved for the new acceleration and velocity (see the module
docstring), on the dynamic degrees of freedom; every other entry of the outputs is zero,
matching the zero rows of the mass and the damping there.
Parameters
----------
dU
The trial displacement increment.
V_n, A_n
The velocity and acceleration at the start of the increment.
dT
The time increment.
beta, gamma
The Newmark parameters.
dynamicDofs
The degrees of freedom integrated in time.
V_np, A_np
The vectors to write the trial velocity and acceleration into.
"""
dUDynamic = np.asarray(dU)[dynamicDofs]
VDynamic = np.asarray(V_n)[dynamicDofs]
ADynamic = np.asarray(A_n)[dynamicDofs]
ANew = (dUDynamic - dT * VDynamic - dT * dT * (0.5 - beta) * ADynamic) / (beta * dT * dT)
VNew = VDynamic + dT * ((1.0 - gamma) * ADynamic + gamma * ANew)
A_np[:] = 0.0
V_np[:] = 0.0
A_np[dynamicDofs] = ANew
V_np[dynamicDofs] = VNew
[docs]
def readStateFromModel(self, U: DofVector):
"""Read the state the next increment starts from out of the model, into the vectors of the equation system just
(re)built: the solution, as in the parent, and the velocity and the acceleration of the fields carrying a mass,
where the last accepted increment (or the checkpoint, or nothing -- a start from rest) left them, and where a
refinement has meanwhile interpolated them onto the nodes it created.
Parameters
----------
U
The solution vector of the equation system, zero, overwritten.
"""
_, self._dynamicFields = self._locateDynamicDofs()
self._V = self.theDofManager.constructDofVector()
self._A = self.theDofManager.constructDofVector()
self.equationSystemLayout.readStateFromModel(U, self._V, self._A, self._dynamicFields)
[docs]
def writeResultsToModel(self, U: DofVector, dU: DofVector, P: DofVector):
"""Write the results of a converged increment to the model, as the parent does, and the velocity and the
acceleration of the fields carrying a mass with them: what field outputs read, what a checkpoint stores, and
what the next equation system reads back.
The force written is the internal force with the inertial and damping forces of the converged state,
:math:`\\boldsymbol{f}_\\mathrm{int} + \\boldsymbol{M} \\ddot{\\boldsymbol{u}}_{n+1} + \\boldsymbol{C}
\\dot{\\boldsymbol{u}}_{n+1}` (see the documentation page *Node field results*): at a support that no load
acts on, the reaction, the inertia of the support's own mass included. A skipped (zero-length or negligible)
increment evaluates no force; the model keeps the force the last accepted increment wrote, of this step or of
the step before.
Parameters
----------
U
The solution.
dU
The solution increment.
P
The internal force, the elements' residual contributions.
"""
increment = self._currentIncrement
if increment is None:
internalInertialAndDampingForce = None
else:
system = increment.system
internalInertialAndDampingForce = self.theDofManager.constructDofVector()
internalInertialAndDampingForce[:] = P + system.massTimes(self._A) + system.C @ np.asarray(self._V)
super().writeResultsToModel(U, dU, internalInertialAndDampingForce, self._V, self._A, self._dynamicFields)
def _updateNewmarkSystem(self, model: FEModel, stepActions: dict) -> _NewmarkSystem:
"""The Newmark operators and state for the parent's current equation system, assembled if
that system is new and reused otherwise.
The parent rebuilds its :class:`~edelweissfe.numerics.dofmanager.DofManager` at the start
of a step and whenever the topology or a constraint's connectivity changes; the operators
are bound to that manager by identity and follow it. The velocity and acceleration are the
ones the rebuild has read from the model (see :meth:`readStateFromModel`). A step changing a
material property mid-step invalidates the density, so the operators are reassembled every
increment of such a step.
A rebuild that follows a recorded **topology change** -- a ``model.topology.history`` grown
since the last build, or since the step started -- additionally re-arms the
initial-acceleration solve and reports what the change did to the conserved quantities (if
there is a previous system to compare with) -- see the module docstring. A step boundary
rebuilds the manager too, but changes no topology, so it is not such a rebuild; and the
first increment of a resumed step that refines is, exactly as it was in the uninterrupted
run.
It takes a grown history, not merely a rebuild: a rebuild triggered by a constraint
re-reporting its connectivity (a contact candidate list, which can tick on any increment)
has moved no node and interpolated nothing, so there is no stale acceleration to replace and
no conservation statement to make. Restarting the acceleration of the whole model on it
would be both wrong and, repeated per increment, expensive.
Parameters
----------
model
The model tree.
stepActions
The step's actions.
Returns
-------
_NewmarkSystem
The system for the current DofManager.
"""
system = self._newmarkSystem
dofManagerChanged = system is None or system.dofManager is not self.theDofManager
if not dofManagerChanged and not stepActions["changematerialproperty"]:
return system
# Safe to advance here: a topology change always rebuilds the manager, so it is never skipped above.
topologyRecords = len(model.topology.history)
topologyChanged = topologyRecords != self._topologyRecordsAtLastBuild
self._topologyRecordsAtLastBuild = topologyRecords
with performancetiming.timeit("newmark system setup"):
dynamicDofs, dynamicFields = self._locateDynamicDofs()
reusedMassAndDamping = None
if (
dofManagerChanged
and system is not None
and not topologyChanged
and not stepActions["changematerialproperty"]
):
reusedMassAndDamping = self._reuseMassAndDamping(system, model)
if reusedMassAndDamping is not None:
Mvij, damping, M, C, massByField = reusedMassAndDamping
else:
Mvij, damping, M, C, massByField = self._assembleMassAndDamping(model, dynamicDofs, dynamicFields)
self._newmarkSystem = _NewmarkSystem(
dofManager=self.theDofManager,
dynamicDofs=dynamicDofs,
dynamicFields=dynamicFields,
Mvij=Mvij,
dampingVIJIndices=damping[0],
dampingVIJValues=damping[1],
M=M,
C=C,
MRowBlocks=self._rowBlocks(M, _MASS_ROWS_PER_BLOCK),
MRowsPerBlock=_MASS_ROWS_PER_BLOCK,
massByField=massByField,
V=self._V,
A=self._A,
)
if dofManagerChanged and topologyChanged:
if self._conservedAtLastAccept is not None:
with performancetiming.timeit("conservation report"):
self.reportTopologyChangeConservation(
self._conservedAtLastAccept,
self._conservedQuantities(self._newmarkSystem, model),
)
self._initialAccelerationPending = bool(self.options["computeInitialAcceleration"])
self.journal.message(
(
"acceleration will be recomputed from equilibrium on the refined mesh"
if self._initialAccelerationPending
else "keeping the interpolated acceleration (computeInitialAcceleration is off)"
),
self.identification,
2,
)
return self._newmarkSystem
@performancetiming.timeit("reuse mass and damping")
def _reuseMassAndDamping(
self, system: _NewmarkSystem, model: FEModel
) -> tuple[VIJSystemMatrix, tuple[np.ndarray, np.ndarray], csr_matrix, csr_matrix, dict] | None:
"""The mass and the damping of ``system``, moved into the layout of the current
:class:`~edelweissfe.numerics.dofmanager.DofManager` without reassembling them -- or None if
that manager differs from the old one in anything the operators depend on.
The parent builds a new manager whenever a constraint reports a change of its connectivity,
which a contact candidate list does on almost every increment. Such a change moves no node
and touches no element: the degree-of-freedom numbering, the elements, their order and so
their slots in the VIJ layout -- which come first, ahead of every constraint's -- are all
where they were, and only the constraints' slots behind them change. Mass and damping live
entirely in the elements' slots (a constraint carries none), so the CSR operators are
unchanged and the VIJ value vectors only need their element part copied into a vector of the
new length. That replaces a loop over every element and the checks on the assembled mass.
Nothing about that is assumed: the numbering, the element set, every element's slot and the
row and column indices of all element slots are compared, and anything that differs falls
back to a full assembly.
Parameters
----------
system
The Newmark system of the previous manager.
model
The model tree.
Returns
-------
tuple[VIJSystemMatrix, tuple[np.ndarray, np.ndarray], csr_matrix, csr_matrix, dict] | None
The mass as a VIJ value vector of the current manager, the damping's VIJ indices and
values (all in the unchanged element slots), both as CSR matrices, and the mass by field
(unchanged, as the fields' slots are) -- or None if they have to be reassembled.
"""
old = system.dofManager
new = self.theDofManager
if new.nDof != old.nDof or new.idcsOfFieldsInDofVector != old.idcsOfFieldsInDofVector:
return None
if new.idcsOfElementsInDofVector.keys() != old.idcsOfElementsInDofVector.keys():
return None
nElementVIJ = 0
for el in model.elements.values():
start = new.idcsOfHigherOrderEntitiesInVIJ[el]
if start != old.idcsOfHigherOrderEntitiesInVIJ[el]:
return None
nElementVIJ = max(nElementVIJ, start + el.getVIJContributionSize())
if not (
np.array_equal(new.I[:nElementVIJ], old.I[:nElementVIJ])
and np.array_equal(new.J[:nElementVIJ], old.J[:nElementVIJ])
):
return None
# Behind the elements' slots lie only constraints', which carry no mass and no damping.
if np.any(np.asarray(system.Mvij)[nElementVIJ:]) or np.any(system.dampingVIJIndices >= nElementVIJ):
return None
Mvij = new.constructVIJSystemMatrix()
Mvij[:] = 0.0
Mvij[:nElementVIJ] = system.Mvij[:nElementVIJ]
self.journal.message(
"constraint connectivity changed only: mass and damping reused, not reassembled",
self.identification,
2,
)
return Mvij, (system.dampingVIJIndices, system.dampingVIJValues), system.M, system.C, system.massByField
def _conservedQuantities(self, system: _NewmarkSystem, model: FEModel) -> _ConservedQuantities:
"""Total mass, linear momentum and kinetic energy of one Newmark system.
Read off that system's own :class:`~edelweissfe.numerics.dofmanager.DofManager`, so the
same method measures the state before and after a rebuild without either one having to
know about the other.
A field occupies a contiguous slice of the dof vector, node-major with the component
innermost -- what ``writeNodeFieldToDofVector``'s flatten establishes -- so reshaping the
slice recovers the per-node vectors, and the row sums of the mass over the slice count
every field component once, hence the division by the field's dimension.
Parameters
----------
system
The system to measure.
model
The model tree, for the fields' spatial dimension.
Returns
-------
_ConservedQuantities
The three quantities.
Raises
------
ValueError
If two dynamic fields differ in spatial dimension, so that their momenta have no
common components to add.
"""
V = np.asarray(system.V)
MV = system.massTimes(V)
return _ConservedQuantities(
massByField=system.massByField,
momentum=linearMomentum(MV, system.dofManager, system.dynamicFields, model),
kineticEnergy=0.5 * float(V @ MV),
)
@staticmethod
def _totalMassByField(M: csr_matrix, dofManager: DofManager, fieldNames: list, model: FEModel) -> dict:
"""Each field's own total assembled mass, never summed across fields.
The row sums of the consistent mass over a field's slice count every spatial component of
the field once, hence the division by the field's dimension.
Parameters
----------
M
The consistent mass.
dofManager
The manager it is indexed by.
fieldNames
The fields whose inertia is a mass.
model
The model tree, for the fields' spatial dimension.
Returns
-------
dict
Field name to its total mass.
"""
rowSums = np.asarray(M.sum(axis=1)).ravel()
return {
fieldName: float(np.sum(rowSums[dofManager.idcsOfFieldsInDofVector[fieldName]]))
/ model.nodeFields[fieldName].dimension
for fieldName in fieldNames
}
@staticmethod
def _rowBlocks(M: csr_matrix, rowsPerBlock: int) -> list[csr_matrix]:
"""``M`` cut into CSR matrices of ``rowsPerBlock`` consecutive rows (the last one may be
shorter), whose products with a vector can run on separate threads.
A block's column indices and values are views of ``M``'s, so the blocks cost only their row
pointers. SciPy computes every row of ``M @ x`` on its own, from zero, adding the row's
products in storage order; a block holds the same entries in the same order, with indices of
the same type, so its rows of ``block @ x`` are bitwise those of ``M @ x``.
Parameters
----------
M
The matrix.
rowsPerBlock
The number of rows of a block.
Returns
-------
list[csr_matrix]
The blocks, in row order.
"""
nRows, nColumns = M.shape
blocks = []
for start in range(0, nRows, rowsPerBlock):
stop = min(start + rowsPerBlock, nRows)
entries = slice(M.indptr[start], M.indptr[stop])
# Built empty and given the views afterwards: the constructor would copy them (SciPy's prune).
block = csr_matrix((stop - start, nColumns), dtype=M.dtype)
block.indptr = M.indptr[start : stop + 1] - M.indptr[start]
block.indices = M.indices[entries]
block.data = M.data[entries]
blocks.append(block)
return blocks
[docs]
def reportTopologyChangeConservation(self, before: _ConservedQuantities, after: _ConservedQuantities):
"""Report what a topology change did to the quantities that ought to survive it: every
dynamic field's total mass is checked (and raises when violated), momentum and kinetic
energy are reported -- see :mod:`~edelweissfe.solvers.base.conservationchecks`
for which of them is exact when. The mass is reassembled here rather than transferred, so
its conservation is an identity of the assembly and the cheapest correctness check on the
whole refinement.
Parameters
----------
before, after
The quantities measured on the system before and on the one after the change.
Raises
------
RuntimeError
If any dynamic field's total mass changed by more than ``massConservationTolerance``.
"""
tolerance = self.options["massConservationTolerance"]
worstField, worstRelativeChange = "", 0.0
for fieldName, massBefore in before.massByField.items():
relativeChange = self._conservationCheck.check(
"mass of field '{:}'".format(fieldName), massBefore, after.massByField.get(fieldName, 0.0), tolerance
)
if relativeChange >= worstRelativeChange:
worstField, worstRelativeChange = fieldName, relativeChange
self.journal.message(
"Topology change: worst-conserved mass, field '{:}', to {:.1e} relative; {:}".format(
worstField,
worstRelativeChange,
formatMomentumAndKineticEnergy(
before.momentum, after.momentum, before.kineticEnergy, after.kineticEnergy
),
),
self.identification,
2,
)
def _locateDynamicDofs(self) -> tuple[np.ndarray, list]:
"""The degrees of freedom integrated in time: those of every field whose inertia is a mass.
Returns
-------
tuple[np.ndarray, list]
Their indices in the dof vector, and the names of their fields in the model's order.
Raises
------
ValueError
If the model carries no such field.
"""
isDynamic = np.zeros(self.theDofManager.nDof, dtype=bool)
dynamicFields = []
for fieldName, indices in self.theDofManager.idcsOfFieldsInDofVector.items():
if carriesLinearMomentum(fieldName):
isDynamic[indices] = True
dynamicFields.append(fieldName)
if not dynamicFields:
raise ValueError(
"No field of this model carries a mass (see phenomena.inertiaKind), so there is nothing "
"for {:} to integrate in time; it carries: {:}.".format(
self.identification, ", ".join(self.theDofManager.idcsOfFieldsInDofVector) or "(no field)"
)
)
return np.flatnonzero(isDynamic), dynamicFields
@performancetiming.timeit("assemble mass and damping")
def _assembleMassAndDamping(
self, model: FEModel, dynamicDofs: np.ndarray, dynamicFields: list
) -> tuple[VIJSystemMatrix, tuple[np.ndarray, np.ndarray], csr_matrix, csr_matrix]:
"""Assemble the consistent mass and the diagonal damping, in the stiffness' VIJ layout and
as CSR matrices.
Entries coupling to a degree of freedom that is not integrated in time are zeroed in both:
the mass an element reports on a non-mechanical block (a micro-inertia) belongs to a field
this solver keeps quasi-static. Elements without kernels -- contact facets -- carry geometry,
not a material, and are skipped.
Parameters
----------
model
The model tree.
dynamicDofs
The degrees of freedom integrated in time.
dynamicFields
Their fields, for the per-field mass report.
Returns
-------
tuple[VIJSystemMatrix, tuple[np.ndarray, np.ndarray], csr_matrix, csr_matrix, dict]
The mass as a VIJ value vector, the damping as the VIJ indices and values of its
nonzero diagonal entries (see :meth:`_lumpedDampingOfElement`), the same two as CSR
matrices, and each dynamic field's total mass.
Raises
------
ValueError
If the assembled mass or damping is not finite, the damping is negative, or the mass
leaves a dynamic degree of freedom without any mass.
"""
Mvij = self.theDofManager.constructVIJSystemMatrix()
# The elements are independent: each writes its mass into its own slice of Mvij and reports its damping, so
# they are done in chunks on the thread pool. The chunks come back in the order of the elements, which makes
# the concatenated damping the one of a serial loop.
def assembleChunk(elements):
dampingIndices, dampingValues = [], []
for el in elements:
try:
el.computeConsistentInertia(Mvij[el])
except NotImplementedError as error:
raise NotImplementedError(
"{:} needs a consistent mass matrix, which element {:} ({:}) does not provide; "
"use elements that implement computeConsistentInertia (the Marmot elements), or "
"the explicit solver NED, which needs only a lumped mass.".format(
self.identification, el.elNumber, type(el).__name__
)
) from error
indices, values = self._lumpedDampingOfElement(el)
dampingIndices.append(indices)
dampingValues.append(values)
return dampingIndices, dampingValues
chunks = mapChunksInParallel([el for el in model.elements.values() if el.hasKernels], assembleChunk)
# Seeded with empty arrays, so that a model without any damping still concatenates to
# empty arrays of the right dtype.
dampingVIJIndices = np.concatenate([np.zeros(0, dtype=np.int64)] + [i for chunk in chunks for i in chunk[0]])
dampingVIJValues = np.concatenate([np.zeros(0)] + [v for chunk in chunks for v in chunk[1]])
nDof = self.theDofManager.nDof
I = self.theDofManager.I # noqa: E741
J = self.theDofManager.J
isDynamic = np.zeros(nDof, dtype=bool)
isDynamic[dynamicDofs] = True
# The three passes over the VIJ entries below are as long as the system matrix has entries (hundreds of
# millions on a large model): in slices on the thread pool, each slice as the pass over the whole would do it.
couplesDynamicOnly = np.empty(len(I), dtype=bool)
def markCouplingOfDynamicDofs(s: slice):
couplesDynamicOnly[s] = isDynamic[I[s]] & isDynamic[J[s]]
forEachSliceInParallel(markCouplingOfDynamicDofs, len(I))
self._warnAboutDiscardedInertia(Mvij, dampingVIJIndices, couplesDynamicOnly)
# Mass and damping on the quasi-static fields are discarded, see the warning above.
massValues = np.asarray(Mvij)
def zeroMassOfQuasiStaticDofs(s: slice):
massValues[s][~couplesDynamicOnly[s]] = 0.0
forEachSliceInParallel(zeroMassOfQuasiStaticDofs, len(I))
keptDamping = couplesDynamicOnly[dampingVIJIndices]
dampingVIJIndices = dampingVIJIndices[keptDamping]
dampingVIJValues = dampingVIJValues[keptDamping]
# The mass matrix has the stiffness' VIJ layout, so K's CSR generator sums it into K's pattern.
M = self.csrGenerator.updateCSR(np.asarray(Mvij))
C = coo_matrix((dampingVIJValues, (I[dampingVIJIndices], J[dampingVIJIndices])), shape=(nDof, nDof)).tocsr()
# Symmetry is checked per element type in the tests; see the module documentation.
if not np.all(np.isfinite(M.data)) or not np.all(np.isfinite(C.data)):
raise ValueError("The assembled consistent mass or damping contains non-finite entries.")
if np.any(C.data < 0.0):
raise ValueError("The assembled damping has negative entries; a damping must dissipate.")
# A dynamic degree of freedom without mass (zero density, or only mass-less entities) would be
# integrated as though it had inertia.
diagonal = M.diagonal()
massless = dynamicDofs[diagonal[dynamicDofs] <= 0.0]
if massless.size:
raise ValueError(
"{:} of {:} time-integrated degrees of freedom received no mass (first: dof {:}). Every "
"element carrying a mass field must report a positive density.".format(
massless.size, dynamicDofs.size, int(massless[0])
)
)
massByField = self._totalMassByField(M, self.theDofManager, dynamicFields, model)
for fieldName, totalMass in massByField.items():
self.journal.message(
"consistent mass assembled: field '{:}' carries a total mass of {:e}".format(fieldName, totalMass),
self.identification,
2,
)
return Mvij, (dampingVIJIndices, dampingVIJValues), M, C, massByField
def _lumpedDampingOfElement(self, el) -> tuple[np.ndarray, np.ndarray]:
"""The lumped (diagonal) damping one element reports, placed on the diagonal of its own
block in the VIJ layout -- only the nonzero entries, as their VIJ indices and values.
An element's block occupies ``nDof * nDof`` consecutive VIJ entries in row-major order, so
its ``k``-th diagonal entry is at ``start + k * (nDof + 1)``. Each slot belongs to one
element, so indices collected over all elements are unique.
Parameters
----------
el
The element.
Returns
-------
tuple[np.ndarray, np.ndarray]
The VIJ indices of the element's nonzero damping entries, and the damping values there.
"""
Ce = np.zeros(el.nDof)
el.computeLumpedDamping(Ce)
nonzero = np.flatnonzero(Ce)
start = self.theDofManager.idcsOfHigherOrderEntitiesInVIJ[el]
return (start + nonzero * (el.nDof + 1)).astype(np.int64), Ce[nonzero]
def _warnAboutDiscardedInertia(
self, Mvij: VIJSystemMatrix, dampingVIJIndices: np.ndarray, couplesDynamicOnly: np.ndarray
):
"""Warn, once per field, when an element reported a nonzero inertia or damping on a field
this solver keeps quasi-static -- a rotational inertia, or a gradient-enhanced element's
micro-inertia -- which is about to be discarded.
Parameters
----------
Mvij
The assembled mass, before the entries are zeroed.
dampingVIJIndices
The VIJ indices of the nonzero damping, before those are dropped.
couplesDynamicOnly
Per VIJ entry, whether it couples two time-integrated degrees of freedom.
"""
massValues = np.asarray(Mvij)
def discardedMassInSlice(s: slice) -> np.ndarray:
return np.flatnonzero(~couplesDynamicOnly[s] & (massValues[s] != 0.0)) + s.start
discarded = np.concatenate(
forEachSliceInParallel(discardedMassInSlice, len(massValues))
+ [dampingVIJIndices[~couplesDynamicOnly[dampingVIJIndices]]]
)
if not discarded.size:
return
I = self.theDofManager.I # noqa: E741
J = self.theDofManager.J
isAffected = np.zeros(self.theDofManager.nDof, dtype=bool)
isAffected[I[discarded]] = True
isAffected[J[discarded]] = True
for fieldName, indices in self.theDofManager.idcsOfFieldsInDofVector.items():
if carriesLinearMomentum(fieldName) or fieldName in self._fieldsWarnedAboutDiscardedInertia:
continue
if np.any(isAffected[indices]):
self._fieldsWarnedAboutDiscardedInertia.add(fieldName)
self.journal.message(
"WARNING: the elements report a nonzero inertia or damping on field '{:}', which {:} "
"does not integrate in time; it is discarded and the field is treated "
"quasi-statically.".format(fieldName, self.identification),
self.identification,
2,
)
@performancetiming.timeit("initial acceleration")
def _computeInitialAcceleration(
self, system: _NewmarkSystem, U_n: DofVector, stepActions: dict, model: FEModel, timeStep: TimeStep
):
"""Compute the acceleration from equilibrium at the start of an increment, and commit it.
Solves :math:`M a_0 = f_\\mathrm{ext}(t_0) - f_\\mathrm{int}(u_0) - C v_0` on the dynamic,
free degrees of freedom. The loads are evaluated with a synthetic time step at the start of
the increment passed in (zero step-progress increment and zero time increment -- an ``f(t)``
amplitude at that instant), the elements with a zero displacement increment, which leaves
their state untouched. Prescribed degrees of freedom get a zero acceleration; degrees of
freedom this solver does not integrate get a unit diagonal and a zero right-hand side so the
system stays regular there.
Called on the first increment of a cold-started step, where :math:`t_0` is the step's own
start, and on the increment after a topology change, where it is that increment's start and
the acceleration it replaces is the one a refinement interpolated -- see
:meth:`_updateNewmarkSystem`. The two are the same computation; only the instant differs.
Parameters
----------
system
The current Newmark system; its acceleration is overwritten.
U_n
The displacement at the start of the increment.
stepActions
The step's actions.
model
The model tree.
timeStep
The increment about to be solved; the start is derived from it.
"""
startTimeStep = TimeStep(
timeStep.number,
0.0,
timeStep.stepProgress - timeStep.stepProgressIncrement,
0.0,
timeStep.stepTime - timeStep.timeIncrement,
timeStep.totalTime - timeStep.timeIncrement,
)
dirichlets = stepActions["dirichlet"].values()
loads = self.collectLoads(stepActions)
P = self.theDofManager.constructDofVector()
F = self.theDofManager.constructDofVector()
PLoadsAndConstraints = self.theDofManager.constructDofVector()
dU0 = self.theDofManager.constructDofVector()
# the tangent is not needed: assembled into the increment's own system matrix, which the Newton iterations
# assemble anew
K = self._K
K.beginAssembly()
P, K, F = self.assembleInternalForces(U_n, dU0, P, K, F, startTimeStep)
PLoadsAndConstraints, K = self.assembleLoads(loads, U_n, dU0, PLoadsAndConstraints, K, startTimeStep)
PLoadsAndConstraints, K = self.assembleConstraints(
model.constraints, U_n, dU0, PLoadsAndConstraints, K, startTimeStep
)
# f_ext - f_int is the negative of the summed residual contributions, as in solveIncrement
R = self.theDofManager.constructDofVector()
R[:] = -P
R -= PLoadsAndConstraints
R -= system.C @ np.asarray(system.V)
nDof = self.theDofManager.nDof
isStatic = np.ones(nDof, dtype=bool)
isStatic[system.dynamicDofs] = False
R[isStatic] = 0.0
MEff = (system.M + diags(isStatic.astype(float), format="csr")).tocsr()
if self.mpcTransformation is not None:
R[:] = self.mpcTransformation.transformResidual(R, dU0)
# MEff is a matrix of its own, not the CSR generator's, whose pattern holds the condensed one
MEff = self.mpcTransformation.transformSystemMatrix(self.mpcTransformation.withCondensedPattern(MEff))
for dirichlet in dirichlets:
R[dirichlet.constrainedDofIndices] = 0.0
MEff = self.applyDirichletToStiffness(MEff, dirichlets)
A0 = self.linearSolve(MEff, R)
if self.mpcTransformation is not None:
A0 = self.mpcTransformation.recoverSlaveSolution(A0)
system.A[:] = 0.0
system.A[system.dynamicDofs] = A0[system.dynamicDofs]
# It completes the state the increment starts from: a rebuild before the increment is accepted -- a retry
# after a cutback whose contact search changes the connectivity -- reads it back from the model.
self.equationSystemLayout.writeVelocityAndAccelerationToModel(system.V, system.A, system.dynamicFields)
self.journal.message(
"acceleration from equilibrium at t = {:g}: ||a0||inf = {:e}".format(
startTimeStep.stepTime,
float(np.max(np.abs(A0[system.dynamicDofs]))) if system.dynamicDofs.size else 0.0,
),
self.identification,
2,
)