3.3 Control Section¶

The control section is defined by the Control element. This section defines all parameters that are used to control the evolution of the solution as well as parameters for the nonlinear solution procedure. These parameters are defined as child elements of the Control element and depend somewhat on the analysis as defined by the Module section. The control section defines the control parameters, the solver parameters, and optionally, the time_stepper parameters.

Note that for multistep analyses (see chapter 4) the Control section is specified for each step separately and defined as a child element of each step section.

Example:

<Control>
    <analysis>STATIC</analysis>
    <time_steps>50</time_steps>
    <step_size>0.02</step_size>
    <solver>
        <max_refs>25</max_refs>
        <diverge_reform>1</diverge_reform>
        <reform_each_time_step>1</reform_each_time_step>
        <dtol>0.001</dtol>
        <etol>0.01</etol>
        <rtol>0</rtol>
        <lstol>0.9</lstol>
        <min_residual>1e-20</min_residual>
        <rhoi>-2</rhoi>
        <qn_method type="BFGS">
            <max_ups>10</max_ups>
        </qn_method>
    </solver>
    <time_stepper>
        <dtmin>0.0002</dtmin>
        <dtmax>0.02</dtmax>
        <max_retries>5</max_retries>
        <opt_iter>6</opt_iter>
    </time_stepper>
</Control>

Control Parameters¶

The following parameters are common for all analysis types. If not specified they are assigned default values, which are found in the last column. An asterisk (*) after the name indicates a required parameter. The numbers behind the description refer to the comments following the table.

Parameter	Description	Default
analysis	Sets the analysis type (1)	static
time_steps*	Total number of time steps. (=_ ntime_)(2)	(none)
step_size*	The initial time step size. (= dt) (2)	(none)
plot_level	Sets the level of state dumps to the plot file (3)	PLOT_MAJOR_ITRS
plot_range	Set the range of the states that will be stored to the plot file (4)	0,-1
plot_stride	Set the stride of the states that will be stored to the plot file (4)	1
plot_zero_state	Flag that controls whether the “zero” state will be written to the plot file, even if it is not defined in the range (4)	0 (false)
output_level	Controls when to output data to file (5)	OUTPUT_MAJOR_ITRS
output_stride	Set the output stride for writing out	1
adaptor_resolve	re-solve time step after mesh adaptations	1 (yes)

Comments:

The analysis element sets the analysis type. The text value of this element is determined by the module, but the numeric value should be 0 for (quasi-)static or steady-state analysis and 1 for dynamic or transient analysis. In a quasi-static analysis, inertial effects are ignored and an equilibrium solution is sought. Note that in this analysis mode it is still possible to simulate time-dependent effects such as viscoelasticity. In a dynamic analysis the inertial effects are included.

Value	Description
static	(quasi-) static analysis
steady-state	steady-state response of a transient (quasi-static) biphasic, biphasic-solute, multiphasic, fluid or fluid-FSI analysis
dynamic	dynamic analysis for solid, fluid and fluid-FSI analyses

The total simulation time of the analysis is determined by ntime * dt. Note that when the auto-time stepper is enabled (see below), the actual number of time steps and time step size may be different than specified in the input file. However, the total simulation time will always be determined by ntime * dt.
The _plot_level _allows the user to control exactly when the solution is to be saved to the plot file. The following values are allowed:

Value	Description
PLOT_NEVER	Don't save the solution to the plot file
PLOT_MAJOR_ITRS	Save the solution after each converged timestep
PLOT_MINOR_ITRS	Save the solution for every quasi-Newton iteration
PLOT_MUST_POINTS	Only save the solution at the must points
PLOT_FINAL	Only store the final converged state.
PLOT_STEP_FINAL	Store only the final converged state of each analysis step.
PLOT_AUGMENTATIONS	Store states at the start of each augmentation.

The PLOT_MUST_POINTS option must be used in conjunction of a must-point curve. See the comments on the dtmax parameter for more information on must-point curves. When the _plot_level _option is set to PLOT_MUST_POINTS, only the time-points defined in the must-point curve are stored to the plotfile.

When using the fixed time stepper, several parameters control which time steps are stored to the plot file.

The plot_range parameter sets the range of the states that will be stored to the plot file. The range is defined by two values that specify the first and last time step that will be stored to the plot file. The value “0” refers to the initial time step, usually time zero, and negative values count backwards from the final time step (as defined by the time_steps parameter). For instance, the default values,

<plot_range>0,-1</plot_range>

store all time steps to the plot file, including the initial “zero” time step. As another example, to store only the last five time steps, set

<plot_range>-5,-1</plot_range>

By default, all time steps within the range will be stored to the plot file. Time steps can be skipped using the _plot_stride _parameter. For instance, to store only every 10 steps, set

<plot_stride>10</plot_stride>

Note that the first and last time step defined by the range will always be stored, regardless of the plot stride.

The “zero” time step refers to the initial state of the model, before any calculations are done. This state will only be stored to the plot file if the minimal value of the plot range is set to zero. To force storing this state to the plot file, set the plot_zero_state parameter to one.

<plot_zero_state>1</plot_zero_state>

This will store the zero time step to the plot file, even when it is not specified inside the plot range.

Again, the plot_range,_ plot_stride_, and plot_zero_state parameters are only used by the fixed time stepper. Currently, these parameters are not used with the auto time stepper.

The output_level can be used to control when FEBio outputs the data files. The following values are supported.

Value	Description
OUTPUT_NEVER	Don't generate any output
OUTPUT_MUST_POINTS	Only output at must points
OUTPUT_MAJOR_ITRS	Output at end of each time step
OUTPUT_MINOR_ITRS	Output at each iteration
OUTPUT_FINAL	Only output the data at the last converged time step.

Time Stepper parameters¶

The optional _time_stepper _element sets the parameters that control the FEBio auto-time stepper. This auto-time stepper will adjust the time step size depending on the convergence stats of the previous time step. It defines the following parameters.

Value	Description	Default
dtmin	Minimum time step size (1)
dtmax	Maximum time step size (3)
opt_iter	Optimal, or desired, number of iterations per time step (2)	10
max_retries	Maximum number of times a time step is restarted. (2)	5
aggressiveness	Algorithm for cutting the time step size after a failed time step (4)	0
cutback	Time step scale factor, used when aggressiveness = 1	0.5

Comments:

The dtmin and dtmax values are used to constrain the range of possible time step values. The _opt_iter _defines the estimated optimal number of quasi-Newton iterations. If the actual number of iterations is less than or equal to this value the time step size is increased, otherwise it is decreased, as detailed further below.
When a time step fails (e.g. due to a negative Jacobian), FEBio will retry the time step with a smaller time step size. The max_retries parameter determines the maximum number of times a timestep may be retried before FEBio error terminates. The new time step size is determined by the ratio of the time step size before restarts started and max_retries+1. For example, if the time step size is 0.1 and _max_retries is set to 4, then the time step size is adjusted by 0.02: The first retry will have a step size of 0.08; the next will be 0.06, and so on. Specifically, the new time step size is determined as follows. (This is only for the case when the aggressiveness _parameter is set to zero.)

\[ \Delta t^{k+1}=\Delta t^{k}-\Delta t^{0}/\left(m_{\mathrm{max}}+1\right) \]

Here, \(\Delta t^{k}\)is the current time step size, \(\Delta t^{k+1}\) is the new time step size, \(\Delta t^{0}\) is the time step size before retries started, and \(m_{\mathrm{max}}\)is the max number of retries (i.e. max__retries).
When the time step succeeds, the time step size will be adjusted by comparing the opt_iter_ parameter to the actual number of iterations it took to converge. If the actual number of iterations is larger than opt_iter, the time step size will be decreased.

\[ \Delta t^{k+1}=\Delta t^{k}-\left(\Delta t^{k}-\Delta t_{\mathrm{min}}\right)\left(1-r\right) \]

Here, \(\Delta t^{k}\)is the current time step size, \(\Delta t^{k+1}\) is the new time step size, \(\Delta t_{\mathrm{min}}\)is the minimum time step size and \(r=\sqrt{n_{\mathrm{opt}}/n}\), where \(n\) is the number of iterations and \(n_{\mathrm{opt}}\)is the optimal number of iterations as specified by the user (i.e. opt_iter).
On the other hand, if the actual number of iterations is less than_ opt_iter_, then the time step size will be increased. The exact equation follows.

\[ \Delta t^{k+1}=\Delta t^{k}+\left(\Delta t_{\mathrm{max}}-\Delta t^{k}\right)\mathrm{min}\left(0.2,r-1\right) \]

In addition, the new time step size is enforced to be less than 5 times the old time step size, to prevent the time step size from growing too fast.
After calculating the new time step size, it will be ensured that the new timestep size falls within the range defined by dtmin and dtmax and that must-points are respected.
The user can specify a load curve for the dtmax parameter. This load curve is referred to as the must-point curve and serves two purposes. Firstly, it defines the value of the dtmax parameter as a function of time. This can be useful, for instance, to enforce smaller time steps during rapid loading or larger time steps when approaching steady-state in a transient analysis. Secondly, the time points of the dtmax loadcurve define so-called must-points. A must-point is a time point where FEBio must pass through. This is useful for synchronizing the auto-time stepper with the loading scenario. For instance, when loading starts at time 0.5, adding a must-point at this time will guarantee that the timestepper evaluates the model at that time. In conjunction with the PLOT_MUST_POINT value of the _plot_level _parameter, this option can also be used to only write results to the plotfile at the specified time points. Consider the following example.

<dtmax lc="1">0.1</dtmax>
...
<loadcontroller id="1" type="loadcurve">
  <interpolate>STEP</interpolate>
    <points>
      <point>0,0</point>
      <point>0.5, 0.1</point>
      <point>1.0, 0.2</point>
    </points>
</loadcontroller>

This example defines load curve 1 as the must-point curve. This curve defines three points where FEBio will pass through (namely 0, 0.5 and 1.0). The values of each time point is the value of the maximum time-step size (dtmax). Since the curve is defined as a step-function, each value is valid up to the corresponding time-point. Thus, between time 0 and time 0.5, the maximum time step value is 0.1. Between 0.5 and 1.0 the maximum time step value is 0.2. If the _plot_level _parameter is set to PLOT_MUST_POINTS, then only the three defined time points will be stored to the plotfile.

Note that when specifying a loadcurve for the dtmax parameter, the value of this parameter will be ignored.

When FEBio fails to converge a time step, the time step size is reduced and then FEBio tries to solve the time step again. The algorithm for determining the reduced time step size, depends on the aggressiveness parameter.
aggressiveness = 0 (default): \(\Delta t^{k+1}=\Delta t^{k}-\Delta t^{0}/\left(m_{\mathrm{max}}+1\right)\), (see comment 2 for more details).
aggressiveness = 1: \(\Delta t^{k+1}=\gamma\Delta t^{k}\), where \(\gamma\) is the cutback factor (i.e. the cutback parameter).

Common Solver Parameters¶

Many of the solvers define the same parameters. They are listed in the table below. Additional parameters that are only defined for a specific type of analysis are listed in the following sections.

Parameter	Description	Default
etol	Convergence tolerance on energy (1)	0.01
rtol	Convergence tolerance on residual (1)	0 (disabled)
lstol	Convergence tolerance on line search (2)	0.9
lsmin	minimum line search step	0.01
lsiter	Maximum number of line search iterations	5
ls_check_jacobians	Check for negative Jacobians during linesearch	0 (disabled)
max_refs	Max number of stiffness reformations (3)	15
qn_method	Quasi-Newton update method (3)	BFGS
max_ups	Max number of BFGS/Broyden stiffness updates (3)	10
diverge_reform	Flag for reforming stiffness matrix when the solution diverges (3)	1
min_residual	Sets minimal value for residual norm (5)	1e-20
max_residual	Sets the max value for residual norm	0 (disabled)
symmetric_stiffness	Use symmetric stiffness matrix flag (6)	1
equation_scheme	Set the equation allocation scheme	0
equation_order	Set the equation allocation ordering	0
optimize_bw	Optimize bandwidth of stiffness matrix (4)	0
linear_solver	Set the preferred linear solver to use. (7)	(default)
check_zero_diagonal	Check for zero diagonals in the stiffness matrix	0 (disabled)
zero_diagonal_tol	Tolerance value for checking zero diagonals	0.0
force_partition	For a partition in the stiffness matrix	0 (disabled)
reform_each_time_step	Reform the stiffness matrix at the start of each time step	1 (enabled)
reform_augment	Reform stiffness matrix after each augmentation	0 (disabled)

Comments:

FEBio determines convergence of a time step based on three convergence criteria: displacement, residual and energy (that is, residual multiplied by displacement). Each of these criteria requires a tolerance value that will determine convergence when the relative change will drop below this value. For example, a displacement tolerance of \(\varepsilon\) means that the ratio of the displacement increment (i.e. the solution of the linearized FE equations, \(\Delta\mathbf{U}=-\mathbf{K}^{-1}_{k}\mathbf{R}_{k}\)) norm at the current iteration \(k+\)1 to the norm of the total displacement (\(\mathbf{U}_{k+1}=\mathbf{U}_{k}+\Delta\mathbf{U}\)) must be less than \(\varepsilon\):

\[ \frac{\left\Vert \Delta\mathbf{U}\right\Vert }{\left\Vert \mathbf{U}_{k+1}\right\Vert }<\varepsilon \]

For the residual and energy norms, it is the ratio of the current residual norm (resp. energy norm) to the initial one that needs to be less than the specified convergence tolerance.

To disable a specific convergence criterion, set the corresponding tolerance to zero. For example, by default, the residual tolerance is zero, so that this convergence criterion is not used.

The lstol parameter controls the scaling of the vector direction obtained from the line search. A line search method is used to improve the convergence of the nonlinear (quasi-) Newton solution algorithm. After each quasi-Newton iteration, this algorithm searches in the direction of the displacement increment for a solution that has less energy (that is, residual multiplied with the displacement increment) than the previous iteration. In many problems this will automatically be the case. However, in some problems that are very nonlinear (e.g. contact), the line search can improve convergence significantly. The line search can be disabled by setting the lstol parameter to zero, although this is not recommended.
The _qn_method property determines the quasi-Newton update method. The particular method is set via the type attribute. Additional parameters are then specified as child elements of this tag. _See section Quasi-Newton Solver Parameters for more details.
The optimize_bw parameter enables bandwidth minimization for the global stiffness matrix. This can drastically decrease the memory requirements and running times when using the skyline solver. It is highly recommended when using the skyline solver.

<optimize_bw>1</optimize_bw>

When using a different linear solver (e.g., pardiso or SuperLU), the bandwidth optimization can still be performed if so desired. However, for these solvers there will be little or no effect since these solvers are not as sensitive to the bandwidth as the skyline solver.

If no force is acting on the model, then convergence might be problematic due to numerical noise in the calculations. For example, this can happen in a displacement driven contact problem where one of the contacting bodies is moved before initial contact is made. When this happens, the residual norm will be very small. When it drops below the tolerance set by min_residual, FEBio will assume that there is no force acting on the system and will converge the time step.
The _symmetric_stiffness _flag is used to set the symmetry of the global stiffness matrix. Only specify this flag if you want to override the default symmetry, which depends somewhat on the analysis type. Forcing a symmetric matrix when the problem is non-symmetric may have negative impact on convergence. The values for this parameter are:
- 0 = unsymmetic
- 1 = symmetric
- 2 = structurally symmetric
- 3 = preferred
If the option is set to 3 (preferred), then one of the other formats will be chosen depending on the model components in the model. For instance, if a non-symmetric contact formulation is used, then a non-symmetric stiffness matrix will be used.
The default linear solver used by FEBio is usually configured in the FEBio configuration file. However, it can also be specified directly in the FEBio input file as part of the solver parameters using the linear_solver tag. See Configuring Linear Solvers for more details on configuring linear solvers.

<linear_solver type="pardiso"/>

Quasi-Newton Solver Parameters¶

Most solvers in FEBio use the quasi-Newton method for solving the nonlinear FE equations. Several different QN algorithms are implemented as described below. The specific QN method is selected using the qn_method property of the solver. The following QN methods are supported.

BFGS The Broyden-Fletcher-G-S method. (default if qn_method is not specified). This method assumes the global stiffness matrix is symmetric. It can still be used with non–symmetric matrices, but the convergence might be compromised.

Broyden The Broyden method is similar to BFGS, but works with both symmetric and non-symmetric matrices. Therefore, it is the recommended method for problems that generate a non-symmetric stiffness matrix, such as (some) contact formulations, biphasic, multi-phasic, fluid, fluid FSI.

JFNK This is an implementation of the Jacobian-Free-Newton-Krylov method. When using this option, you must also use an iterative linear solver such as FGMRES.

Both the BFGS and Broyden methods calculate the global stiffness matrix at the beginning of each time step. For each iteration, a matrix update is then done. The maximum number of such updates is set with max_ups. When FEBio reaches this number, or if the solution diverges and_ diverge_reform_ is set to 1, it reforms the global stiffness matrix (that is, it recalculates it) and factorizes it, essentially taking a "full Newton" iteration. Then FEBio continues with BFGS/Broyden iterations. The _max_refs _parameter is used to set the maximum of such reformations FEBio can do, before it fails the timestep. In that case, FEBio will either terminate or, if the auto-time stepper is enabled, retry with a smaller time step size.

The BFGS and Broyden take the following parameters.

Parameter	Description	Default
max_ups	This is the maximum number of updates between stiffness matrix reformations.	10
cmax	Sets the max condition number.	1e5
max_buffer_size	Sets the max buffer size.	0 (ignored)
cycle_buffer	Recycle buffers when max_ups is larger than max_buffer_size	1

Note that when max_ups is set to 0, FEBio effectively will use a Full-Newton method since the stiffness matrix is reformed each iteration. In this case it is recommended to increase the number of_ max_refs_ (to e.g. 50) in the solver settings, since the default value might cause FEBio to terminate prematurely when convergence is slow.

Solver Parameters for a Structural Mechanics Analysis¶

A structural mechanics analysis is defined by using the solid type in Module section. The solver parameters are given in the table below:

Parameter	Description	Default
dtol	convergence tolerance on displacement	0.001
rhoi	Spectral radius parameter \(\rho_{\infty}\) (1)	-2
alpha	generalized alpha method parameter alpha	1
beta	generalized alpha method beta	0.25
gamma	generalized alpha method gamma	0.5
logSolve	Use an acceleration for Newton iterations	0 (off)
arc_length	Specifies arc-length method (2)	0 (off)
arc_length_scale	Set the arc-length scale factor	0

Comments:

There are two time-integrators implemented in FEBio: the classical Newmark integration and the generalized alpha method. Setting rhoi to -2, will use the former, where as setting rhoi to the range [0,1] will use the latter. In the latter case, the values for gamma, beta, and gamma are determined from the spectral radius parameter, rhoi.
When using the arc-length method, there are some important restrictions. Only zero-prescribed boundary conditions are allowed. The loads must be specified via nodal loads and cannot be time-dependent. The time_steps and step_size control parameters are interpreted as arc-length steps and step_size respectively.

When using the “CG-solid” solver, the following parameters are defined.

Parameter	Description	Default
cgmethod	Select solution algorithm (1)	0
preconditioner	Select the preconditioner (2)	0
lsiter	maximum number of line search iterations	10

Comments:

The solution algorithm can be:
0: Hager-Zhang conjugate gradient (This is the default and should be preferred since its faster.)
1: steepest descent
Set 0 for no preconditioner, 1 for diagonal stiffness preconditioner. The diagonal stiffness preconditioner greatly improves convergence when there are different size elements, mixed materials, or quadratic elements with midside nodes.

<solver type="CG-solid">
  <cgmethod>0</cgmethod>
  <preconditioner>1</preconditioner>
  <lsiter>10</lsiter>
</solver>

Solver Parameters for Biphasic Analysis¶

A biphasic analysis is defined by using the biphasic type in Module section. Since a biphasic analysis couples a fluid problem to a solid mechanics problem, all control parameters above can be used in a biphasic analysis. In addition, the following parameters can be defined:

Parameter	Description	Default
ptol	Specify the fluid pressure convergence tolerance	0.01
mixed_formulation	Use a mixed formulation (1)	0

Comment:

The mixed formulation uses discontinuous shape functions for the pressure interpolation. It still uses the full shape functions for displacement.

Solver Parameters for Solute and Multiphasic Analyses¶

When the type attribute of the Module section is set to solute or multiphasic, an analysis is solved that includes solute transport. All parameters for a biphasic analysis can be used (including the ones for a structural mechanics analysis). In addition, the following parameters can be specified:

Parameter	Description	Default
ctol	Specify the concentration convergence tolerance	0.01
force_positive_concentrations	Force the concentration values always to be positive.	1 (enabled)

Solver Parameters for Fluid and Fluid-FSI Analyses¶

A fluid analysis is defined by using the fluid type in Module section (see Section Module Section). In addition to the common parameters, the following parameters can be specified:

Parameter	Description	Default
vtol	convergence tolerance on velocity	0.001
ftol	convergence tolerance on dilatation	0.001
rhoi	Spectral radius parameter \(\rho_{\infty}\)	0
reform_each_time_step	Flag for reforming stiffness matrix at the start of each time step	1
predictor	Sets the predictor method	0
min_volume_ratio	Sets the minimum value for the volume ratio.	0 (disabled)

A fluid-structure interaction analysis is defined by using the fluid-FSI type in Module section. It uses the parameters of solid and fluid analyses.

Transient fluid and fluid-FSI analyses may often run very efficiently using Broyden's method (qnmethod set to 1) with max_ups set to 50; efficiency may be further increased by setting reform_each_time_step to 0, which will postpone reforming the stiffness matrix at subsequent time steps until max_ups updates have been exhausted. When using Broyden's method with fluid and fluid-FSI analyses, set rtol to a non-zero value to ensure an accurate solution, for example rtol=0.001 (see Section Control Parameters).

The spectral radius parameter \(\rho_{\infty}\) determines the time integration scheme: Values in the range \(0\le\rho_{\infty}\le1\) use the generalized \(\alpha-\)method (see FEBio Theory Manual). With \(\rho_{\infty}=0\), damping of higher frequency components (such as those produced by a step increase in velocity) occurs theoretically within a single time step; in contrast, with \(\rho_{\infty}=1\), no damping occurs, potentially leading to instability in the solution process. When solving transient flows that exhibit eddies or shed vortices, a value of \(\rho_{\infty}=0.5\) represents a good balance between too much and too little numerical damping.

Explicit-Solid Solver¶

FEBio supports an explicit-solid solver, which solves dynamic structural mechanics problems using an explicit time-integration scheme. This means that the solution at a given timestep is calculated using only the solution from the previous timestep. The advantage of this approach is that no stiffness matrix needs to be calculated and no linear system of equations needs to be solved (at least, with mass-lumping enabled). As a result, this solver can solve dynamic problems very fast. However, the downside of this approach is that it is only conditionally stable and imposes an upper limit on the time step size. This limit is often much smaller than the time step sizes that can be taken by implicit time-integration schemes (as used by the default Newton-based solvers in FEBio) and thus often will require many more time steps to solve the model. Thus, the explicit solver is most useful in dynamic problems that already require relatively small time step sizes, such as impact-contact models, or models that use non-elastic material behavior that may undergo rapid changes in their constitutive behavior (e.g. damage).

The explicit solver in FEBio uses the midpoint time integration rule and uses mass lumping. Mass lumping is a technique that replaces the consistent mass matrix with a diagonal approximation. In fact, it's this diagonalization of the mass matrix that really prevents the need to solve a linear system of equations. (Since solving a diagonal linear system of equations is trivial.)

The explicit solver can be selected by setting the type attribute of the solver element to explicit-solid.

<solver type="explicit-solid">
...
</solver>

The following parameters can be specified for the explicit-solid solver.

parameter	description	default value
mass_lumping	Chooses the mass lumping strategy. (1)	2
dyn_damping	The damping factor for dynamic damping. (2)	1

Table 1. The parameters of the explicit-solid solver

Comments:

As mentioned, mass lumping is a technique to approximate the mass matrix with a diagonal matrix. Several approaches are supported in FEBio.
mass_lumping = 0: no mass-lumping (This is currently not implemented, but the value is reserved in case support will be added in the future.)
mass_lumping = 1: use row summation approach. Useful for linear elements, but may produce negative nodal masses for quadratic elements and thus should not be used for models containing higher-order elements.
mass_lumping = 2: HRZ lumping technique. Better suited for higher-order elements, but should not be used for shells.
Dynamic damping is a technique to dampen the solution in order to obtain an equilibrium solution (e.g. for false-transient analyses). In FEBio, dynamic damping is applied to the velocity updates and the _dyn_damping _parameter specifies the damping factor. The default value is 1 so no damping is applied. Note that the lower the damping factor, the more damping, but also the longer it may take to find the equilibrium solution.

Important! As a final note, the explicit-solver is still considered an experimental feature in FEBio and may not work correctly with all of the features of FEBio yet. Use at your own risk and please report any problems you may encounter to the FEBio developers.