#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------
#
# _____ _ _ _ _____ _____
# | ____|__| | ___| |_ _____(_)___ ___| ___| ____|
# | _| / _` |/ _ \ \ \ /\ / / _ \ / __/ __| |_ | _|
# | |__| (_| | __/ |\ V V / __/ \__ \__ \ _| | |___
# |_____\__,_|\___|_| \_/\_/ \___|_|___/___/_| |_____|
#
#
# Unit of Strength of Materials and Structural Analysis
# University of Innsbruck,
# 2017 - today
#
# Alexander Dummer alexander.dummer@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.
# ---------------------------------------------------------------------
from dataclasses import dataclass
import numpy as np
import edelweissfe.utils.performancetiming as performancetiming
from edelweissfe.config.timing import createTimingDict
from edelweissfe.constraints.base.constraintbase import ConstraintBase
from edelweissfe.models.femodel import FEModel
from edelweissfe.numerics.csrgeneratorv2 import CSRGenerator
from edelweissfe.numerics.dofmanager import DofVector
from edelweissfe.numerics.systemmatrixassembly import (
SystemMatrixAssembly,
VIJSystemMatrixAssembly,
)
from edelweissfe.outputmanagers.base.outputmanagerbase import OutputManagerBase
from edelweissfe.solvers.nonlinearimplicitstatic import NIST
from edelweissfe.timesteppers.timestep import TimeStep
from edelweissfe.utils.exceptions import CutbackRequest, IncrementFailed
from edelweissfe.utils.fieldoutput import FieldOutputController
from edelweissfe.utils.schema import schemaField
def getRungeKuttaParameters(rungeKuttaStages: int) -> tuple[dict, dict, dict]:
_alpha = {}
_omega = {}
_lambda = {}
if rungeKuttaStages == 2:
# parameters for stage 2
_alpha[21] = 1
# parameters for increment estimate
_omega[1] = 0.5
_omega[2] = 0.5
# parameters for error control
_lambda[1] = -0.5
_lambda[2] = 0.5
elif rungeKuttaStages == 3:
# parameters for stage 2
_alpha[21] = 2 / 3
# parameters for stage 3
_alpha[31] = 0
_alpha[32] = 2 / 3
# parameters for increment estimate
_omega[1] = 2 / 8
_omega[2] = 3 / 8
_omega[3] = 3 / 8
# parameters for error control
_lambda[1] = 0
_lambda[2] = -3 / 8
_lambda[3] = 3 / 8
elif rungeKuttaStages == 4:
# parameters for stage 2
_alpha[21] = -0.5
# parameters for stage 3
_alpha[31] = 0.6684895833
_alpha[32] = -0.2434895833
# parameters for stage 4
_alpha[41] = -2.323685857
_alpha[42] = 1.125483559
_alpha[43] = 2.198202298
# parameters for increment estimate
_omega[1] = 0.03431372549
_omega[2] = 0.02705627706
_omega[3] = 0.7440130202
_omega[4] = 0.1946169772
# parameters for error control
_lambda[1] = 0.03431372549
_lambda[2] = -0.01262626262
_lambda[3] = -0.0289338397
_lambda[4] = 0.00724637679
elif rungeKuttaStages == 6:
# parameters for stage 2
_alpha[21] = 0.2
# parameters for stage 3
_alpha[31] = 3 / 40
_alpha[32] = 9 / 40
# parameters for stage 4
_alpha[41] = 3 / 10
_alpha[42] = -9 / 10
_alpha[43] = 6 / 5
# parameters for stage 5
_alpha[51] = -11 / 54
_alpha[52] = 5 / 2
_alpha[53] = -70 / 27
_alpha[54] = 35 / 27
# parameters for stage 6
_alpha[61] = 1631 / 55296
_alpha[62] = 175 / 512
_alpha[63] = 575 / 13824
_alpha[64] = 44275 / 110592
_alpha[65] = 253 / 4096
# parameters for increment estimate
_omega[1] = 37 / 378
_omega[2] = 0
_omega[3] = 250 / 621
_omega[4] = 125 / 594
_omega[5] = 0
_omega[6] = 512 / 1771
# parameters for error control
_lambda[1] = -277 / 64512
_lambda[2] = 0
_lambda[3] = 6925 / 370944
_lambda[4] = -6925 / 202752
_lambda[5] = -277 / 14336
_lambda[6] = 277 / 7084
else:
raise NotImplementedError("runge-kutta-stages must be 2, 3, 4 or 6")
return _alpha, _omega, _lambda
[docs]
@dataclass(frozen=True)
class NESTSchema:
"""The options of the ``*solver`` datalines and of an ``>>options`` block routed to this
solver, owned by this module and never mutated from outside it.
Mirrors :attr:`NEST.SolverSpecificOptions` one-for-one; the plain ``self.options`` dict remains
the actual source of truth consulted at runtime (see :class:`NISTSchema` for why). The
``runge-kutta-*`` option names are not valid Python identifiers, hence the ``optionName``
indirection.
"""
rungeKuttaStages: int | None = schemaField(
description="The number of Runge-Kutta stages.", dtype=int, default=2, optionName="runge-kutta-stages"
)
rungeKuttaErrorTolerance: float | None = schemaField(
description="The error tolerance for the Runge-Kutta error control.",
dtype=float,
default=1e-3,
optionName="runge-kutta-error-tolerance",
)
rungeKuttaErrorControl: str | None = schemaField(
description="Activate the Runge-Kutta error control (on|off).",
dtype=str,
default="on",
optionName="runge-kutta-error-control",
)
linsolver: str | None = schemaField(
description=(
"The linear solver: the name of a *linearSolver block, or a type of linear solver with its defaults (e.g. "
"pardiso, superlu, blockamg)."
),
dtype=str,
default="pardiso",
)
# Read by the inherited NIST.applyDirichletToStiffness, so NEST must declare it too.
pruneCondensedMatrixZeros: bool | None = schemaField(
description=(
"Compact explicitly stored zeros out of the multi-point-constraint-condensed system "
"matrix before solving. See NISTSchema for the trade-off; unset (the default) decides "
"automatically from the linear solver in force."
),
dtype=bool,
default=None,
documentedDefault="auto",
)
[docs]
class NEST(NIST):
"""This is the Nonlinear Explicit STatic -- solver.
Parameters
----------
jobInfo
A dictionary containing the job information.
journal
The journal instance for logging.
"""
identification = "NESTSolver"
supportsMPC = False
supportsModelModifiers = False
#: Option schema for this solver, per OptionSchemaProvider.
schema = NESTSchema
SolverSpecificOptions = {
"runge-kutta-stages": 2,
"runge-kutta-error-tolerance": 1e-3,
"runge-kutta-error-control": "on",
"linsolver": "",
"pruneCondensedMatrixZeros": None,
}
#: Not restartable: no restart scenario covers this solver, so it does not claim to be.
checkpointedState = None
def __init__(self, jobInfo, journal, **kwargs):
self.journal = journal
self.options = self.SolverSpecificOptions.copy()
# the *solver datalines belong only to this solver: unknown entries are typos
self._updateOptions(kwargs, journal, strict=True)
#: The ``*linearSolver`` blocks of the input file by name, which linsolver= may name (none without an input file).
self.linearSolverDefinitions = jobInfo.get("linearSolvers", {})
#: The linear solver, and the ``linsolver`` it was made for (see :meth:`selectLinearSolverForStep`).
self.linSolver = None
self.linSolverName = None
[docs]
def validateModelCapabilities(self, model: FEModel):
"""Refuse, on top of the base checks, a constraint that is switched on and needs its connectivity updated in
every increment, i.e. a contact, which searches for what it touches: this solver never performs the update. A
constraint switched off from outside
(:attr:`~edelweissfe.constraints.base.constraintbase.ConstraintBase.switchedOn`) is not refused: it rests,
and takes no part in the step.
Parameters
----------
model
The model tree.
"""
super().validateModelCapabilities(model)
for constraintName, constraint in model.constraints.items():
if constraint.switchedOn and type(constraint).updateConnectivity is not ConstraintBase.updateConnectivity:
raise Exception(
f"Constraint '{constraintName}' requires a dynamic connectivity update "
f"(contact) every increment, which {self.identification} never performs -- "
"contact is not currently supported with this solver."
)
[docs]
def beginStep(
self,
step,
model: FEModel,
fieldOutputController: FieldOutputController,
outputmanagers: dict[str, OutputManagerBase],
):
"""Build the equation system, once per step, and start the step; see
:meth:`~edelweissfe.solvers.base.nonlinearsolverbase.NonlinearSolverBase.beginStep`.
Parameters
----------
step
The step to be solved.
model
The model tree.
fieldOutputController
The field output controller.
outputmanagers
The output managers.
"""
performancetiming.reset()
self.validateModelCapabilities(model)
self.journal.message("Creating monolithic equation system", self.identification, 0)
self.forgetEquationSystem()
self.equationSystemLayout = model.buildEquationSystem()
self.theDofManager = self.equationSystemLayout.dofManager
self.journal.message(
"total size of eq. system: {:}".format(self.theDofManager.nDof),
self.identification,
0,
)
self.journal.printSeperationLine()
self.computationTimes = createTimingDict()
# get parameters for runge kutta scheme
self.rkAlpha, self.rkOmega, self.rkLambda = getRungeKuttaParameters(self.options.get("runge-kutta-stages", 2))
self.rkStages = self.options.get("runge-kutta-stages", 2)
self.tol = self.options.get("runge-kutta-error-tolerance", 1e-3)
self.selectLinearSolverForStep()
# NEST solves a linear system per Runge-Kutta stage (see solveIncrement, via the inherited
# linearSolve()), so its linear solver needs the setModel NIST gives its own: for solvers that
# derive anything beyond the plain (A, b) call -- a field-split solver such as blockamg raises
# without the field structure setModel supplies. NEST's equation system is built once per step here
# rather than being rebuilt on connectivity changes, so one call per step is enough.
self.linSolver.setModel(self.equationSystemLayout.modelForLinearSolver, self.theDofManager)
U = self.theDofManager.constructDofVector()
P = self.theDofManager.constructDofVector()
dU = self.theDofManager.constructDofVector()
self.csrGenerator = CSRGenerator(self.theDofManager)
K = self.createSystemMatrixAssembly(VIJSystemMatrixAssembly)
self.readStateFromModel(U)
self._U, self._dU, self._P, self._K = U, dU, P, K
self.prevTimeStep = None
self.applyStepActionsAtStepStart(model, step)
[docs]
def prepareIncrement(self, step, model: FEModel, isRetry: bool):
"""Nothing to prepare: this solver has neither topology updates nor contact."""
[docs]
def attemptIncrement(self, step, model: FEModel, timeStep: TimeStep):
"""Integrate the increment with the embedded Runge-Kutta scheme; see
:meth:`~edelweissfe.solvers.base.nonlinearsolverbase.NonlinearSolverBase.attemptIncrement`.
Parameters
----------
step
The step being solved.
model
The model tree.
timeStep
The increment.
Raises
------
IncrementFailed
If the scheme requests a smaller increment.
"""
self.incrementStatus = {
"step": step.number,
"inc": timeStep.number,
"iters": None,
"converged": False,
"time inc": timeStep.timeIncrement,
"time end": timeStep.totalTime,
"notes": "",
}
self.journal.printSeperationLine()
self.journal.message(
"increment {:}: {:8f}, {:8f}; time {:10f} to {:10f}".format(
timeStep.number,
timeStep.stepProgressIncrement,
timeStep.stepProgress,
timeStep.totalTime - timeStep.timeIncrement,
timeStep.totalTime,
),
self.identification,
level=1,
)
try:
U, dU, P, incScaleFactor = self.solveIncrement(
self._U,
self._dU,
self._P,
self._K,
step.actions,
model,
timeStep,
self.prevTimeStep,
)
except CutbackRequest as e:
self.prevTimeStep = None
self.incrementStatus["notes"] = str(e)
raise IncrementFailed(str(e), max(e.cutbackSize, 0.25)) from e
self._U, self._dU, self._P = U, dU, P
self._incScaleFactor = incScaleFactor
[docs]
def acceptIncrement(self, step, model: FEModel, timeStep: TimeStep):
"""Commit the increment to the model, and let the error estimate size the next one; see
:meth:`~edelweissfe.solvers.base.nonlinearsolverbase.NonlinearSolverBase.acceptIncrement`.
Parameters
----------
step
The step being solved.
model
The model tree.
timeStep
The solved increment.
"""
U, dU, P = self._U, self._dU, self._P
step.timeStepper.changeIncrementSize(self._incScaleFactor)
self.prevTimeStep = timeStep
self.writeResultsToModel(U, dU, P)
self.updateRigidBodies(model, timeStep)
model.advanceToTime(timeStep.totalTime)
[docs]
def endStep(self, step, model: FEModel):
"""Report the step's performance timing.
Parameters
----------
step
The step that was solved.
model
The model tree.
"""
prettyTable = performancetiming.makePrettyTable()
self.journal.printPrettyTable(prettyTable, self.identification)
performancetiming.reset()
[docs]
def solveIncrement(
self,
U_n: DofVector,
dU: DofVector,
P: DofVector,
K: SystemMatrixAssembly,
stepActions: list,
model: FEModel,
timeStep: TimeStep,
prevTimeStep: TimeStep,
) -> tuple[DofVector, DofVector, DofVector, float]:
"""Explicit Runge-Kutta scheme to solve for an increment.
Parameters
----------
Un
The old solution vector.
dU
The old solution increment.
P
The old reaction vector.
K
The system matrix to be used.
elements
The dictionary containing all elements.
stepActions
The list of active step actions.
model
The model tree.
timeStep
The current time step.
prevTimeStep
The previous time step.
Returns
-------
tuple[DofVector,DofVector,DofVector,int,dict]
A tuple containing
- the new solution vector
- the solution increment
- the new reaction vector
- the scale factor for the next increment size
"""
constraints = model.constraints
R = self.theDofManager.constructDofVector()
F = self.theDofManager.constructDofVector()
PLoadsAndConstraints = self.theDofManager.constructDofVector()
U_np = self.theDofManager.constructDofVector()
dU_ = [self.theDofManager.constructDofVector() for i in range(self.rkStages)]
err = self.theDofManager.constructDofVector()
err[:] = 0.0
dirichlets = stepActions["dirichlet"].values()
loads = self.collectLoads(stepActions)
# Find which global DOFs the Dirichlet BCs constrain, once up front.
self.locateConstrainedDofs(dirichlets)
self.applyStepActionsAtIncrementStart(model, timeStep, stepActions)
for geostatic in stepActions["geostatic"].values():
geostatic.applyAtIterationStart()
dU[:] = 0.0
# iterate over runge kutta stages
for k in range(self.rkStages):
U_np[:] = U_n
dU_[k][:] = 0.0
if k > 0:
for j in range(k + 1 - 1):
a = (k + 1) * 10 + 1 + j
dU_[k] += self.rkAlpha[a] * dU_[j]
U_np += dU_[k]
P[:] = F[:] = PLoadsAndConstraints[:] = 0.0
K.beginAssembly()
P, K, F = self.assembleInternalForces(U_np, dU_[k], P, K, F, timeStep)
PLoadsAndConstraints, K = self.assembleLoads(loads, U_np, dU_[k], PLoadsAndConstraints, K, timeStep)
PLoadsAndConstraints, K = self.assembleConstraints(
constraints, U_np, dU_[k], PLoadsAndConstraints, K, timeStep
)
# the right-hand side is the negative of all residual contributions, R = f_ext - f_int
R[:] = -P
R -= PLoadsAndConstraints
R = self.applyDirichletToResidual(timeStep, R, dirichlets)
K_ = self.assembleStiffnessCSR(K)
# identity rows, and the columns eliminated into R
K_ = self.applyDirichletToStiffness(K_, dirichlets, R)
# solve for increment
dU_[k] = self.linearSolve(K_, R)
# update solution increment
dU += self.rkOmega[k + 1] * dU_[k]
# update error measure
err += self.rkLambda[k + 1] * dU_[k]
# update solution
U_np[:] = U_n
U_np += dU
normErr = np.linalg.norm(err) / (np.linalg.norm(U_np) + 1e-16)
self.journal.message("Estimated error {:5.2e}".format(normErr), self.identification, 2)
if normErr > self.tol and self.options.get("runge-kutta-error-control", "on") != "off":
raise CutbackRequest("Estimated error too high!", np.sqrt(self.tol / normErr) * 0.9)
# compute elements once more to get the final reaction
P[:] = F[:] = 0.0
K.beginAssembly()
P, K, F = self.assembleInternalForces(U_np, dU, P, K, F, timeStep)
# compute new increment size scale factor
if self.options.get("runge-kutta-error-control", "on") == "off":
beta = 1.0
else:
beta = max(0.1, min(np.sqrt(self.tol / (normErr + 1e-15)) * 0.9, 2))
return U_np, dU, P, beta