Bayesian Inversion with Open-Source Codes for Various One-Dimensional Model Problems in Computational Mechanics

The complexity of many problems in computational mechanics calls for reliable programming codes and accurate simulation systems. Typically, simulation responses strongly depend on material and model parameters, where one distinguishes between backward and forward models. Providing reliable information for the material/model parameters, enables us to calibrate the forward model (e.g., a system of PDEs). Markov chain Monte Carlo methods are efficient computational techniques to estimate the posterior density of the parameters. In the present study, we employ Bayesian inversion for several mechanical problems and study its applicability to enhance the model accuracy. Seven different boundary value problems in coupled multi-field (and multi-physics) systems are presented. To provide a comprehensive study, both rate-dependent and rate-independent equations are considered. Moreover, open source codes (https://doi.org/10.5281/zenodo.6451942) are provided, constituting a convenient platform for future developments for, e.g., multi-field coupled problems. The developed package is written in MATLAB and provides useful information about mechanical model problems and the backward Bayesian inversion setting.


Introduction
Bayesian inversion has been used in different applied and engineering problems to identify specific material parameters that either cannot be measured with usual techniques, or for which significant experimental efforts are needed to provide a reasonable estimation. In mechanical systems, the uncertainty arises from the spatial fluctuation (varying properties within space), uncertainty in the dimensional measurements, model structural errors, or misalignment of the specimen or measurement device [1]. Therefore, considering the effect of uncertainties is crucial for having a reliable model.
Markov chain Monte Carlo (MCMC) methods are among the most effective sampling techniques for obtaining knowledge and estimating posterior densities in Bayesian inversion. These methods are effective since they can be implemented easily for different scientific and engineering problems [2]. However, traditional methods are computationally inefficient in high-dimensional parameter spaces. There are effective improvements to enhance MCMC convergence to the target density. An efficient remedy is the adaptation of the proposal, which can be done locally and globally [3,4]. In this work, we introduce some common MCMC techniques and discuss their computational features. Then, they will be used to identify the unknown material parameters in different mechanical problems.
In computational mechanics, a Bayesian framework, employing the Metropolis-Hastings algorithm was introduced in [5] to accurately capture the mechanical behavior of solids. In the context of continuum thermodynamics, a fundamental study for model calibration in solid mechanics was given in [6]. The MCMC based on the Metropolis-Hastings algorithm has been employed to facilitate the Bayesian system identification of a nonlinear dynamical system from experimentally observed acceleration time histories [7]. The delayed acceptance Metropolis-Hastings has been used in geomechanics [8,9] to estimate the posterior density, where the forward model describes the flow in porous media with or without fractures as well as coupled hydro-mechanical processes. The Bayesian parameter identification of viscoplastic-damaging models was designed in [10] by developing a polynomial chaos expansion version of Kalman filter. This has been further evaluated with a CT-test and by considering the model error [11]. Subsequently, different Bayesian inversion models including Transitional Markov Chain Monte Carlo method and Gauss-Markov Kalman filter approach were used to compare and identify the material parameters of viscoplastic damaging materials [12]. An overview on Bayesian inversion in solid mechanics to incorporate several uncertainty sources in the model problem is provided in [1].
Due to the complexity of these mechanical problems, an efficient parameter estimation setting provides means to enhance the accuracy and reliability of the computations. Employing Bayesian inversion enables us to calibrate the forward model and provides an accurate estimation of the model parameters. In [50], the authors use MCMC methods in ductile fracture, where a Bayesian setting is developed to infer the mechanical characteristics. Here, we use the local and global proposal adaptation to improve the MCMC efficiency and employ it in different mechanical and thermomechanical systems. The approach provides informative knowledge about the influence of the mechanical characteristics, i.e., which parameter has more effect on the system. Furthermore, the convergence study (hereR-statistics) monitors the fast convergence of the MCMC chains to the target density and guarantees that all parameters converge accurately at the same time.
In mechanical problems, uncertainties in parameters arise from several sources, e.g., spatial variability and thermal diffusivity. As an example, the material stiffness usually presents strong spatial fluctuation. Direct estimation of these parameters is often difficult and gives rise to significant experimental efforts. Therefore, providing a robust and efficient computational framework to estimate reliable information of the mechanical parameter is of utmost importance. In [51], we presented a probabilistic method to identify the parameters in brittle fracture using the Metropolis-Hastings algorithm. Later, we extended the parameter identification to a more complicated PDE-based model for hydraulic fracture in the isotropic and anisotropic setting [9]. The efficient MCMC approach enabled us to identify several parameters at the same time. Recently, the authors extended the idea of parameter estimation to ductile materials in the context of a unified phase-field framework [50], where experimental data was considered for parameter identification.
An important part of this paper is the accompanying open-source programming codes published on zenodo 1 . In recent years, we have published several studies in computational mechanics and phase-field fracture without explicitly providing codes 2 . However, we strongly believe and agree with the reasons for open-source research software initiatives outlined in [53][54][55], that is, that codes must be available at some point. In order to demonstrate the problem statements, numerical experiments, implementation and documentation, we concentrate in this work on one-dimensional model problems. The purpose is two-fold: firstly, to verify the results of the current work, and secondly, to provide prototype codes of our previous published studies [49][50][51]. Together with the algorithms published in these papers, the current codes allow for extensions towards two-and three-dimensional problems.
In summary, the objectives of this work are: 1.
To provide compact open-source codes 2 for the Bayesian inversion method for one-dimensional analysis in the context of: i Elasticity, ii Heat equation with convection, iii Elastoplasticity, iv Phase-field fracture for brittle materials, v Phase-field fracture for ductile materials (gradient-extended plasticity), vi Phase-field fracture for thermoelasticity, vii Phase-field fatigue fracture.

To investigate convergence performance of different
MCMC techniques in the aforementioned boundary value problems.
The paper is structured as follows. In Sect. 2, we introduce different MCMC techniques and compare their computational features. In Sect. 3, we introduce seven mechanical problems and employ the presented probabilistic inverse technique to precisely estimate the effective material parameters for the associated boundary value problems. Specifically, in this section, we highlight the efficiency of the MCMC techniques. An overview of the codes including convergence analysis is given in Sect. 4. Finally, some concluding remarks are elaborated in Sect. 5.

Parameter Estimation Based on Bayesian Inference
In this section, we review the how efficient MCMC techniques can be used to provide accurate information of material parameters. As a probabilistic model, Bayesian inversion is used to update the available information (prior knowledge) and provide more accurate data (i.e., the posterior density). As a result, a very good agreement between the measurement and the model response (as an aim of the inference) can be achieved. First of all, we define the following statistical model where S denotes an n-dimensional measurement/reference vector, L refers to the model response affected by the k-dimensional model parameters h which belong to the random field H and the position x, and e is the observation error, e.g., a Gaussian independent and identically distributed error e $ N ð0; r 2 IÞ, including the parameter r 2 .
In Bayesian estimation, a parametric forward model (e.g., a PDE-based model) is used to update the available data (considered as random variables) based on the available information (denoting the prior knowledge). Having a measurement S ¼ obs, the conditional density has the following form: Considering the parameter estimation setting, by a specific observation m, we strive to estimate the posterior distribution pðhjmÞ, where prior information is available. Employing the statistical model (1), the likelihood function can be estimated as pðSjhÞ ¼ Lðh; r 2 jSÞ ¼ 1 where is the sum of square errors. The random walk Metropolis (RWM) algorithm was introduced by Metropolis et al. [56]. The algorithm starts with a given initial guess h 0 according to the prior density. Then, based on the proposal distribution (the current state of the chain is perturbed randomly) a new candidate that depends on the chain (the previous candidates) h H is proposed. The acceptance rate is computed by where the likelihood function is given in (4). Either the new candidate is accepted, or the MCMC chain follows in the current way. Here, a symmetric proposal density is used. A generalized form of the algorithm using a nonsymmetric was introduced by Hastings [57], where the probability of the forward jump is not equal to the backward one. This algorithm has been widely used to estimate the desired distribution of the parameters, specifically when the problem dimensionality is low. In contrast, in the highdimensional case, the rejection rate increases dramatically. Furthermore, since the proposal density has no tunning, the convergence is significantly slow. It will be more pronounced when an unsuitable initial guess is used. Although the algorithm is easy to implement and versatile, its computational drawbacks (mostly in multi-dimensional domains) motivate to develop improvements to the Metropolis-Hastings (MH) algorithm.
In the algorithmic treatment, we first determine a range for each parameter (based on the prior information). Considering N iterations, we assign an initial guess (prior knowledge) h 0 for the proposals. During the iterations (j ¼ 2; . . .; N) the next proposal h H is estimated as where COV denotes the Cholesky decomposition of the given proposal function and randn(k,1) is a k-dimensional normally distributed random variable. The function controls whether the new candidate is in the determined range. If each element of the proposal exceeds the maximum value, this element is replaced by the upper bound. If the value is less than the minimum of the range, it is replaced by the minimum. In the MH algorithm, the fixed (not adopted) covariance matrix is used during the MCMC chain. If the acceptance rate k is greater than a produced random variable between 0 and 1 (rand), the proposal is accepted, and the chain is followed with the new proposal. Otherwise, we continue with the previous candidate h jÀ1 . A summary of the algorithm is given in Algorithm 1.

Algorithm 1 The Metropolis-Hastings algorithm.
Initialization: set prior data θ 0 and number of samples N . while j < N

Adaptive Metropolis (AM)
The adaptive Metropolis algorithm [3] is based on tuning the proposal density. This is done based on the obtained information from the previous samples (path of the MCMC chain) by computation of the covariance function. According to the available information, we can adapt the proposal distribution as Here, the parameter is chosen close to zero to guarantee that the covariance is positive definite, and K p ¼ 2:38 2 j is chosen according to [2]. The covariance function is given by [2]. To enhance the efficiency of the method, the proposal can be adapted after a specific number of iterations, e.g., 100. We introduce a new proposal employing the Cholesky decomposition of the covariance function (9): where I j is the identity matrix. In contrast to (6), we adopt the proposal iteratively (e.g., each 100 iterations). By using the updated proposal, the chain will have sufficient information from the accepted candidates which enhances the method efficiency. In this step we employ the code introduced in [4], named covupd. Finally, the acceptance/rejection probability is estimated by

Delayed Rejection (DR)
As previously mentioned, the MH algorithm suffers from computational disadvantages, mostly in a high-dimensional parameter space. The delayed rejection [58] is based on giving another chance to the rejected candidate by using a second stage move. For this, we change the proposal density to increase the acceptance probability of the rejected candidate. The essential parameter is c 2 , which will be selected in the range c 2 \1 such that the next stage has a narrower proposal function (normally, c 2 ¼ 1=5 is chosen). The alternative proposal h ÃÃ is chosen using the proposal function where V j is the covariance matrix estimated by the adaptive algorithm [58]. We use the following acceptance ratio The DR process can be carried out once, or at random/fixed number of iterations.

Delayed Rejection Adaptive Metropolis (DRAM)
In the DRAM algorithm, we follow two steps in order to reduce the rejected candidates and enhance the convergence. Globally, we adapt the proposal density based on the already obtained information from the chain (AM). Locally, we pursue the DR strategy that in each iteration, the rejected candidate will have at least one additional chance. In the case that the initial values are too far from the desired density (target distribution), the adaptive strategy may not be started appropriately. This is due to the large variance of a near singular covariance matrix. As a solution, the DR strategy reduces the variance in different delayed stages, so that some accepted proposals are available, enabling the AM to start appropriately. Therefore, the DRAM technique can effectively overcome the drawbacks of the MH algorithm [4]. A summary of the process is given in Algorithm 2.
The algorithm includes two important aspects • Increasing the acceptance rate by a local adaptation (DR). • Increasing MCMC chain accuracy by a global adaptation (AM).

MCMC with Ensemble-Kalman Filter (EnKF-MCMC)
As previously mentioned, a good detection of the proposal density will enhance the Markov chain movement to the target density. EnKF represents the error covariance matrix by a large stochastic ensemble of model realizations [59]. In this method, to have an accurate estimation of posterior density, a Kalman gain (similar to a weighting matrix in particle filter) is calculated employing the mean and the covariance of the prior distribution and the cross-covariance between parameters and observations. Here, we introduce another proposal distribution detection technique using an ensemble-Kalman filter [60], where one has where Dh indicates the jump of Kalman-inspired proposal. For updating the proposed candidates, we can rewrite (13) as Here K points out the so-called Kalman gain, where C hM denotes the covariance matrix between the identified unknowns and the PDE-based model, C MM is the covariance matrix of the model solution (i.e., the PDE response), and R indicates the measurement noise covariance matrix [61]. Furthermore, y jÀ1 indicates the residual of candidates with respect to the model and s jÀ1 $ N ð0; RÞ refers to the density of measurement. Denoting obs as an observation, y jÀ1 ¼ obs À f ðh jÀ1 Þ. In order to enhance the method efficiency, i.e., reduce the effect of the initial guess, different MCMC chains instead of a single chain can be employed.

Algorithm 2
The DRAM algorithm.
while j < N 1. Propose a new candidate θ * = θ j−1 + R j Z j where R j is the Cholesky decomposition of V j and Z j ∼ Uniform (0, I j ) where I j denotes the identity matrix.
After the estimation of h H , we check its bounds, and calculate its acceptance rate k that determines if the new proposal is rejected or if the chain should be continued with h jÀ1 (the previous candidate). We summarized the procedure in Algorithm 3. The main advantage of EnKF-MCMC lies in using the covariance between the model response and between parameters and the PDE-response. This approach increases the possibility of acceptance of the candidate noticeably. Table 1 provides an informative summary of the computational characteristics of the discussed MCMC techniques. The MH algorithm yields a more efficient performance when compared to RWM since an asymmetric proposal density is used. Also, the implementation of the method is very simple. AM overcomes the drawback of MH by adapting the proposal based on the obtained information. The DR method improves AM, specifically when the initial guess is not suitable and a local adaptation is used. DRAM takes advantage of the benefits of both AM and DR (global and local chain adjustment). EnKF-MCMC has an effective proposal improvement and leads to significant improvements for large number of unknowns. The theoretical results show that DRAM and EnKF perform effectively compared to other methods. In the sequel, we will validate this discussion with practical examples.

One-Dimensional Examples
In this section, we aim at describing how Bayesian inference can be employed to identify material parameters in computational mechanics. Seven different boundary value problems including complex coupled multi-field (and multi-physics) systems are considered. Moreover, compact open-source codes are provided.
In order to introduce Bayesian inference as simply as possible, one-dimensional analyses of uniaxial tests are performed. This can be extended to a multidimensional setting in a straightforward manner. Additionally, Bayesian inference codes to reproduce all the examples are given in detail.

Example 1: Linear Elasticity
This example provides a brief illustrative parameter estimation based on Bayesian inference for one-dimensional linear elasticity. Specifically, we estimate Young's modulus as the single uncertain parameter in the given problem. The observation, i.e., the reference response, is taken as the analytical solution of the governing second-order ordinary differential equation.

State Variables and Kinematics
Let us consider the one-dimensional boundary value problem (BVP) shown in Fig. 1a. An isothermal elastic rod is completely fixed on the left side and subjected to an imposed traction per unit area t on the right side. The external body force qA bðxÞ has the sinusoidal form qA bðxÞ ¼ qðxÞ ¼ q 0 sinð2px=lÞ, as shown in Fig. 1b, where b is a body force per unit mass. Moreover, the deformation process is assumed to be quasi-static, such that inertial effects can be neglected.
The material is assumed to obey linear elasticity in the small strain setting. Consequently, letting B :¼ ½0; L, the displacement field u : B ! R constitutes the only primary variable. To elaborate the corresponding variational formulation, we employ the function space for the displacement field: In particular, we denote W u 0 when we have zero boundary condition (i.e., u ¼ 0). Additionally, the strain field is obtained as The stress field then follows from Hooke's law as where E is Young's modulus, assumed to be constant throughout the bar for the sake of simplicity. With these along with the Dirichlet boundary condition u ¼ 0 at x ¼ 0 and the Neumann boundary condition du=dx ¼ t=E at x ¼ L. Here, A is the cross-sectional area of the bar, assumed to be constant for simplicity. The second-order ordinary differential Eq. (19) has the exact solution for the displacement field: The analytical solution is taken here as the reference response. For the numerical solution, the finite element method is employed, for which the variational formulation of the BVP is presented below.

Energy Quantities and Variational Principles
The strain energy density for one-dimensional linear elasticity reads The total energy functional EðuÞ is then given by

Update the Kalman gain
4. Accepted/rejected the forecast.

Set
where the external forces evaluated at the Neumann boundary o N B fLg have been used (Fig. 1a). Minimization of the energy functional (22) with respect to the displacement field leads to the following weak formulation.
Here, the aim of the Bayesian inversion is to determine Young's modulus E.

Bayesian Inference
In the present example, we consider a bar of unitary length L ¼ 1 with x 2 B :¼ ½0; 1 that is initially unstretched, q 0 ¼ 14 N/m, and a unit cross-sectional area throughout the bar The statistical problem consists of identifying Young's modulus E. Bayesian inference is then applied to determine the target value, while no specific information (other than its positivity) is available a priori.
From now on, in order to estimate the model response, the function is employed. Then, response is used to calculate the likelihood of the given proposal. The main script for the present example is named elasticity 1D. We assume Young's modulus is between E ¼ 0:1 GPa and E ¼ 20 GPa, the initial guess is h 0 ¼ 3, and we use r 2 ¼ 0:1. As a synthetic measurement, we employ E ¼ 4:259 to estimate the displacement. As a one-dimensional Bayesian inversion (k ¼ 1), we plan to compare the performance of all methods, including the convergence of the MCMC chain, elapsed CUP time, and the convergence rate. Figure 2 shows the evolution of the MCMC chain using different Bayesian techniques. All chains converged effectively and EnKF shows a significantly lower fluctuation. As a result, as Fig. 3 illustrates that a narrow posterior distribution has been obtained using the MCMC techniques. Again, EnKF-MCMC shows a better performance (a narrower density) compared to others. The posterior density provides informative knowledge about the parameter. We extract the median of the accepted candidates as the inferred solution. Table 2 shows the computational features of our Bayesian inversion. Although all methods estimate the Young's modulus with a reasonable accuracy, using DRAM and EnKF results in a more precise inference. As expected, the acceptance rate of MH is minimum and it increases by a global adaptation (adaptive Metropolis). The delayed rejection increases significantly the acceptance ratio; however, it also raises the computational time. Combining AM and DR again increases the acceptance ratio by 10%. Since we adopt the proposal every 100 iterations, AM does not noticeably increase the computational time. Regarding EnKF-MCMC, the CPU time is much higher; however, using Kalman again leads to an excellent ratio and a very low uncertainty in MCMC; see

Example 2: Heat Equation with Convection Contribution
The heat equation is a convection-diffusion equation used to model heat transfer from regions of higher temperature to regions of lower temperature. Convective transport is a physical process that takes place along the streamlines in the flow field. In the second example, we apply Bayesian inversion to a one-dimensional convection-diffusion equation. Herein, the goal is to estimate thermal diffusivity and thermal velocity.

State Variables and Kinematics
In the context of the one-dimensional bar treated in the previous example, we now consider a thermal process occurring in a time interval T :¼ ½0; T. In particular, at material points x 2 B and time t 2 T, the BVP solution provides the temperature field T : B Â T ! R. To elaborate the corresponding variational formulation, we employ the function space for the temperature field In particular, we use W T 0 when we have zero boundary condition ( T ¼ 0). In the present one-dimensional problem, the Dirichlet boundary is set at the left boundary and the ending points (right boundary), and the problem includes the following boundary conditions: along with Neumann boundary conditions where Q is the thermal volume flux and q is the thermal transport on the Neumann boundary, where n ¼ 1 for x ¼ L and n ¼ À1 for x ¼ 0, and T is the imposed temperature on the Dirichlet boundary. Moreover, the heat flux in (26) can be described by the negative direction of the thermal heat gradient T 0 (denoting o ox T ) through Fourier's law of heat conduction as where K is the thermal conductivity.

Time Discretization of the Temporal Variables
In the following, we shortly comment on the time-discretization of the temporal variables. To do so, let the interval T :¼ ðt 0 ; TÞ be discretized using the discrete time (loading) points with the end time value T [ 0. The parameter t 2 T denotes for rate-dependent problems the time, and for rateindependent problems an incremental loading parameter. In the following, to formulate the incremental problems, we denote the current time t :¼ t nþ1 [ 0 with known fields at t n . In order to advance the solution within a specific time step, we focus on the finite time increment ½t n ; t nþ1 , where denotes the time step length.

Strong Formulation and Variational Principles
The one-dimensional diffusion equation with convection can be written as where c is the heat capacity, q is the heat density, and j indicates the thermal velocity. Typically, to identify the relative importance between convection and diffusion in a heat equation, the dimensionless quantity called Péclet number is estimated by Pe :¼ ðj hÞ=ð2 KÞ, where h is the element size. This value expresses the ratio of convective to diffusive transport in a heat flow problem.
For test functions dT 2 W T 0 (see (24) with zero boundary condition), the heat problem given in (30) has the following variational form: The standard variational equation (initiated from a Galerkin formulation) for the heat equation including convection and diffusion effects where convection is dominant compared with diffusion typically leads to unstable results, and thus numerical stabilization has to be used [62,63]. More specifically, without stabilization techniques, some spurious oscillations (or wiggles) appear, which prevent the numerical solution from converging to a physical response. It can be observed that stronger convective effects lead to larger oscillations in time. It is shown for the value Pe [ 1, where convective response is dominating, that the solution obtained through standard FEM presents spurious oscillations [64,65]. To remedy these shortcomings, we consider here the modified test function where b opt is given by and coth indicates the hyperbolic cotangent. Observe that for the linear interpolation basis functions used for the finite element discretization, we have such that the modified weak formulation reads: This variational equation typically refers to the Streamline-Upwind (SU) formulation. An algorithm for the update of the temperature field T in the increment ½t n ; t nþ1 results in time-discrete field variables in an incremental setting of the convection-diffusion equation. Equation (36) then becomes in terms of the x-algorithmic expressions given by such that GðtÞ is a temporal weighting function. The scalar value x for different variations of GðtÞ is given in Table 3.
Regarding the analytical solution, we will consider (25). Thus, the solution has the following form: where The reader is referred to [67] for further details.
The unknown material parameters for the convection-diffusion problem denoted as F2 then follow from the probabilistic Bayesian inversion P BI À F2ðuÞ Á with u :¼ T and P :¼ ðj; KÞ:

Bayesian Inference
For the simulation, we consider a bar of unitary length h is imposed on the left end, i.e., x ¼ 0; see Fig. 4b. Temperature boundary conditions are set according to (25) such that T 1 ¼ 1 C and T 2 ¼ 0 C. Additionally, we set initial temperature as The Bayesian inversion aims to determine two thermal parameters, namely (j, K). We set c ¼ 1:2 J/K and q ¼ 1:6 W=m 2 and use the online MATLAB codes presented in [68]. We assume j is between 5 and 20 m=s with K varying from 0.5 to 2 W=mK uniformly. The script heat 1D has a comparison between all methods considering a two-dimensional Bayesian inversion (simultaneous interference of j and K). Finally, we set h 0 ¼ ð1:2; 15Þ and r 2 ¼ 0:05. In order to provide a reference observation, we use K ¼ 0:9 and j ¼ 15, which can be found in reference:mat.
We plot the evolution of the MCMC chain of K and j as shown respectively in Figs. 5 and 6. MH shows a significant uncertainty in both parameters. However, EnKF and DRAM estimated K and j with a low fluctuation. Figure 7 depicts a comparison between the posterior densities using all MCMC techniques. As shown, the MH method failed to estimate the parameters effectively since a wide distribution is obtained, and the estimated values are relatively far from the true values. The proposal adaptation enhanced the performance of MH; however, the results are similar to a bimodal distribution. The local adaptation (DR) gives rise to a better inference and the peak point of both distributions was located very close to the true values. In general, EnKF and DRAM showed an excellent performance.
The computational features of the methods are listed in Table 4. The median of the posterior densities of all methods are close to the true values. Here, MH and AM have a low (less than 10%) acceptance rate and EnKF has an excellent ratio. Finally, the evolution of the temperature using the estimated parameters from DRAM and EnKF is shown in Fig. 8. The results are compared with the analytical solution, where a very good agreement is observed.

Example 3: Elastoplasticity
As a third example, we apply Bayesian inversion to a onedimensional elastoplasticity problem. Herein, the goal is to estimate Young's modulus E, in addition to the material parameters characterizing the elastoplastic response, namely, the plastic yield strength r Y and the hardening modulus H. Moreover, due to the dissipative nature of the problem, time evolution is explicitly considered [69].
The plastic strain e p : B Â T ! R is considered as a local internal variable. Moreover, to account for phenomenological hardening effects, a hardening variable a : B Â T ! R þ is introduced as a second local internal variable. For simplicity, a is assumed to coincide with the equivalent plastic strain, such that The hardening variable starts to evolve from the initial condition aðx; 0Þ ¼ 0. The elastoplastic response of the bar is thus characterized by the displacement field u, the plastic strain e p , and the hardening variable a.
With the purpose of stating variational principles, we further introduce the function spaces for the plasticity variables W p :¼ L 2 ðBÞ; W a a n ; q :¼ fa 2 H 1 ðBÞ : a ¼ a n þ jqj; The hardening law (42) is then enforced in incremental form as a 2 W a a n ; e p Àe p n and a n 2 W a a nÀ1 ; e p Àe p nÀ1 . In particular, we use W u 0 , and W a 0;q when we have zero boundary condition (i.e., u ¼ 0, and a n ¼ 0 respectively).

Energy Quantities and Variational Principles
Let us define the following pseudo-energy density per unit volume: Wðe; e p ; aÞ :¼ W elas ðe; e p Þ þ W plas ðaÞ: The elastic contribution is assumed to have the simple quadratic form W elas ðe; e p Þ : Accordingly, the plastic contribution in (44) is assumed to have the form with the initial yield stress r Y ! and the isotropic hardening modulus H ! 0. Note that W elas can be viewed as a stored energy density, while W plas can be viewed as a dissipated energy density, defined here as a state function. From the strain energy density, the constitutive stressstrain relation and the generalized stress follow as r ¼ oW elas oe ¼ Eðe À e p Þ and s p ¼ À oW elas oe p r: ð47Þ The elastoplasticity model additionally requires evolution equations for the plastic strain. To this end, the yield function is first defined in terms of the generalized stress s p and the resisting force R p as With the yield function at hand, an associative flow rule is adopted for the plastic strain, such that where k p ! 0 is the plastic multiplier. The flow rule is complemented with the set of Kuhn-Tucker conditions b 0; k p ! 0; and b k p ¼ 0: To derive the weak form of the evolution problem, the total energy functional for the elastoplasticity model is defined as Minimization of the energy functional (22) with respect to the displacement field leads to the following weak formulation.
where e p ðx; tÞ evolves according to the evolution Eqs. (49) and (50), while aðx; tÞ evolves according to the hardening law (42). The unknown material parameters for elastoplasticity denoted as F3 then follow from the probabilistic Bayesian inversion P BI À F3ðuÞ Á with u :¼ ðu; e p ; aÞ and P :¼ ðE; H; r Y Þ: One can show that the evolution Eqs. (42), (49) and (50), which completely determine the evolution of the plasticity variables ðe p ; aÞ, can also be obtained in weak from through minimization of the energy functional (51) with respect to ðe p ; aÞ. Due to the local nature of the plasticity variables in the present example, only the weak form of the mechanical balance equation is included in Formulation 3.3.

Bayesian Inference
In this numerical example, we consider a bar of unitary length L ¼ 1 with x 2 B :¼ ½0; 1 that is initially unstretched and unplasticized. Its left end is fixed, i.e., x ¼ 0, while on its right end, i.e., x ¼ 1, a monotonic displacement increment D u ¼ 5 Â 10 À4 mm is applied for 151 time steps; see Fig. 4a. We set ða; e p Þ ¼ ð0; 0Þ at t ¼ 0. Regarding the finite element mesh size, 100 elements are used. The main script is EP_1D.m and the code forwardmodel.m computes the solution with respect to the given candidates.
This problem is a good example to compare EnKF-MCMC and DRAM. The unknown ranges are mentioned in Table 5. The true values are used to produce the reference observation (reference.mat). We use h 0 ¼ ðE 0 ; H 0 ; r 0 Y Þ ¼ ð60 000; 200; 300Þ. In order to implement the Bayesian inversion, we require the available range for the parameters (the prior density), for which we adopt the uniform distribution summarized in Table 5.
The inferred values were used to estimate the load-displacement curve shown in Fig. 10. Here, we use 8 000 samples to estimate the posterior density. As Fig. 9 shows, the chains estimated by the DRAM algorithm have a low uncertainty for E and r Y ; however, the variation for H is notable. On the other hand, the fluctuation of all unknowns in EnKF is negligible. Of course, this is an important advantage of using the Kalman filter compared to DRAM. The posterior densities verify this statement, where a wide distribution is achieved for H. For the yield stress, both methods performed similarly and have a relatively similar distribution. Also, the acceptance rate of EnKF is 98.4% and DRAM is 6.9%. The medians are summarized in Table 5, indicating the accuracy of both techniques (there is an excellent agreement between the simulations and the true values). In total, due to the higher acceptance ratio and significantly lower uncertainty (specifically of hardening), we conclude that EnKF is more suitable for the elastoplasticity example. Finally, we use the estimated values of EnKF and DRAM to plot the load-displacement curve during different time-steps; see Fig. 10.

Example 4: Phase-Field Fracture for Brittle Materials
The fourth example addresses the variational phase-field modeling of brittle fracture. Starting with the work of [13,15] and the subsequent developments [16][17][18][19][21][22][23], this framework has gained significant popularity in the past decade due to its ability to naturally handle the complex behavior of fractured solids. A Bayesian inversion framework was applied to phasefield modeling of brittle fracture in [51]. In the present illustrative example, we consider a one-dimensional setting, where the objective is to estimate Young's modulus E and the fracture toughness G c for a fixed length scale parameter (although the length scale can also be estimated).

State Variables and Kinematics
The response of a fracturing solid obeying a brittle phasefield model is governed by a two-field system of PDEs. In particular, at material points x 2 B and time t 2 T, the solution of the coupled system provides the displacement field u : B Â T ! R and the phase-field fracture variable d : B Â T ! ½0; 1. Here, dðx; tÞ ¼ 0 and dðx; tÞ ¼ 1 characterize an undamaged and a completely fractured material point, respectively. With the purpose of stating variational principles, we introduce the function spaces In particular, we use W u 0 and W d 0 when we have zero boundary condition (i.e., u ¼ 0 and d ¼ 0, respectively). Here, o D B f0; Lg denotes the Dirichlet boundary, with uð0; tÞ ¼ 0. On the other hand, W d d n is a non-empty, closed, and convex subset of H 1 ðBÞ. This function space introduces the evolutionary character of the phase-field, incorporating an irreversibility condition in incremental form, where d n is the damage value at a previous time instant.

Energy Quantities and Variational Principles
The energy density per unit volume in the present model is additively decomposed into stored elastic strain energy and crack-surface dissipated energy. In particular, one has Wðe; dÞ :¼ W elas ðe; dÞ þ W frac ðd; d 0 Þ: The elastic contribution is assumed to have the simple quadratic form W elas ðe; dÞ : where g(d) is a degradation function for which the usual quadratic form is adopted. On the other hand, the fracture contribution represents the regularized crack-surface energy density, expressed in terms of the current damage state as well as its spatial gradient, and modulated by a length scale parameter. In particular, we adopt the so-called AT-2 version of [28], where In (56), G c [ 0 is the fracture toughness (more specifically, Griffith's critical energy release rate) and l f [ 0 is the fracture length scale parameter that governs the regularization.
With the above definitions, a global energy functional is defined as Minimization of the energy functional (57) with respect to the displacement field and the crack phase-field leads to the following weak formulation.
The unknown material parameters for phase-field fracture denoted as F4 then follow from the probabilistic Bayesian inversion P BI À F4ðuÞ Á with u :¼ ðu; dÞ and P :¼ ðE; G c Þ: The inequality in () 2 is a consequence of the irreversibility condition represented by the constraint d n 2 W d d nÀ1 . To resolve this issue and solve () 2 as an equality, the history field approach, proposed in [20] and widely employed in the phase-field literature, is adopted in this work.

Bayesian Inference
In this numerical example, we consider a bar of unitary length L ¼ 1 with x 2 B :¼ ½0; 1 that is initially unstretched and undamaged. Its left end is fixed, i.e., x ¼ 0, while on its right end, i.e., x ¼ 1, a monotonic displacement increment D u ¼ 1 Â 10 À5 mm is applied for 151 time steps. The example setup is shown in Fig. 4a. Regarding the finite element mesh size, 300 elements are used.
We use a two-dimensional Bayesian inversion to identify (E, G c ). The reference solution is a 1 Â 151 vector computed by E ¼ 70 500 MPa and G c ¼ 0:027 kN/mm 2 (the true values). For the Bayesian inference, we assume E is between 50 000 and 100 000 and G c is between 0.02 and 0.03. As the initial guess, we set h 0 ¼ ð60 000; 0:025Þ. The script name is brittle 1D:m.
Using EnKF-MCMC with 1 000 iterations leads to the estimation of MCMC and also the posterior density. Figure 11 depicts the MCMC evolution for both parameters. A narrow posterior density for both unknowns points out the low uncertainty and fast convergence obtained employing EnKF-MCMC. The acceptance ratio of the chain is 12.2%. The estimated medians are E ¼ 70493 and G c ¼ 0:027, which highlights the method accuracy. Finally, we use these values to estimate the load-displacement curve, and the results are shown in Fig. 12.

Example 5: Phase-Field Fracture for Gradient-Extended Ductile Materials
In this example, the phase-field model for brittle fracture explored above is extended to the ductile case. To this end, the evolution of the phase-field variable is coupled to elastoplasticity, allowing for the initiation and propagation of ductile cracks. Initial works on this topic include [32][33][34][35][36]38] (see [29,50] for an overview). Moreover, to achieve a regularized plastic response, the extension to gradient-extended plasticity [27,39] is considered, which has shown to ensure mesh-objective and physically meaningful responses in the post-critical stage. The present Bayesian inversion framework was applied to the phase-field modeling of ductile fracture in our recent work [50], where three dimensional simulations were performed and compared to experimental observations. Herein, we present an illustrative one-dimensional example, where the objective is to estimate the parameters involved in both the elatoplasticity example (Sect. 3.2) and the brittle fracture example (Sect. 3.4).

State Variables and Kinematics
The response of a fracturing solid obeying the present ductile phase-field model with gradient plasticity is governed by a three-field system of PDEs and a local plastic evolution equation. In particular, at material points x 2 B and time t 2 T, the solution of the coupled system provides the displacement field u : B Â T ! R, the phase-field fracture variable d : B Â T ! ½0; 1 introduced in Sect. 3.4, and the plasticity variables e p : B Â T ! R and a : B Â T ! R þ introduced in Sect. 3.2. Being a phasefield model, as in Sect. 3.4, the constitutive state includes the gradient of Notethe crack phase-field d 0 ðx; tÞ. Moreover, in contrast with Sect. 3.2, we now render the equivalent plastic strain a a non-local internal variable to regularize the plastic response in post-critical stages. As such, the constitutive state further includes the spatial gradient a 0 ðx; tÞ.
Note, with the purpose of stating variational principles, we consider the function spaces (52) for the displacement field, the crack phase-field, and further for plasticity variables by (43).

Energy Quantities and Variational Principles
The energy density per unit volume in the present model is additively decomposed into stored elastic strain energy and the sum of plastic and crack-surface dissipation. In particular, one has Wðe; e p ; a; d; a 0 ; d 0 Þ :¼ W elas ðe; e p ; dÞ þ W plas ða; d; a 0 Þ þ W frac ðd; d 0 Þ: In agreement with (45) and (54), the elastic contribution reads W elas ðu; e p ; dÞ : where g(d) is the quadratic degradation function defined in (55). The stress-strain relation then reads The plasticity contribution is defined as a gradient-extended damaged version of (46): where, for simplicity, and as commonly done in the literature, the same quadratic degradation function (55) is employed to degrade the plastic energy density. Finally, the fracture contribution reads Note that, as opposed to (56), the fracture energy density is here defined as a linear function of d, representing a socalled AT-1 type formulation, where an elastic stage is included in the deformation process [29]. In (56), w 0 [ 0 is a critical fracture energy density, which can be directly related to the fracture toughness in the case of brittle fracture, while l d [ 0 is the fracture length scale parameter. It is worth mentioning that (61) and (62) imply that the dissipated energy in the present model is a state function. With the above definitions, a global energy functional is defined as Eðu; e p ; a; dÞ ¼ Minimization of the energy functional (63) with respect to the displacement field, the crack phase-field, and the plasticity variables leads to the following weak formulation.
The unknown material parameters for ductile fracture denoted as F5 then follow from the probabilistic Bayesian inversion P BI À F5ðuÞ Á with u :¼ ðu; e p ; a; dÞ and P :¼ ðE; H; r Y ; w 0 ; l p Þ: As in Sect. 3.4, the history field approach [20] is employed to solve the variational inequality () 3 . On the other hand, in () 2 , one makes use of the relation da ¼ jde p j and the incremental hardening law a n ¼ a n þ je p n À e p nÀ1 j. This problem can be solved in a convenient way by eliminating e p through the incremental flow rule e p n ¼ e p nÀ1 þ ða n À a nÀ1 Þ signðr trial n Þ, where r trial n :¼ rðe n ; e p nÀ1 ; d n Þ, with r given in (60); see, e.g., [50]. For loading paths without plastic flow reversals, e.g., monotonic loading with possible elastic unloading, the 1D problem is further simplified by setting e p a.

Bayesian Inference
In this numerical example, we consider a bar of unitary length L ¼ 1 with x 2 B :¼ ½0; 1 that is initially unstretched, unplasticized and undamaged. Its left end is fixed, i.e., x ¼ 0, while on its right end, i.e., x ¼ 1, a monotonic displacement increment D u ¼ 5 Â 10 À4 mm is applied for 151 time steps; see Fig. 4a. For the spatial discretization, we use 100 elements. The main script is named ductile_1D.m.
In this example we use a four-dimensional Bayesian inversion to identify ðE; r Y ; w 0 ; HÞ simultaneously. To this end, EnKF is chosen and we set an initial guess of E 0 ¼ 60 000, r 0 Y ¼ 300, H 0 ¼ 300, w 0 0 ¼ 20. The parameter ranges are listed in Table 6 including the true values. Figure 13 shows the MCMC chains where they converge very fast. The posterior densities are depicted in Fig. 14, where for all unknowns, a narrow distribution is obtained. The medians of the MCMC chains are mentioned in Table 6, showing a very good accuracy. These values are used to plot the load-displacement curve shown in Fig. 15. In this example, following the new proposed the stepwise Bayesian inversion for ductile phase-field fracture [50], we consider four stages of Bayesian inversion to identify the parameters E (elastic stage), r Y (perfectly plastic stage), H (hardening stage), and w 0 (damaging stage). The uniform prior densities including the true values are listed in Table 6. Such a step-wise Bayesian framework is computationally effective since a single parameter is estimated at each stage (and used as input in subsequent stages). The posterior densities (using the last 4 000 iterations) are illustrated in Fig. 13. TheR-convergence of the parameters is depicted in Fig. 14 indicating the rapid convergence of all unknowns. Finally, the loaddisplacement curve obtained using the estimated parameters is shown in Fig. 15, indicating a typical ductile fracture process.

Example 6: Phase-Field Fracture for Thermoelastic Materials
Crack propagation due to the thermal effect occurs in different engineering problems such laser fracture, additive manufacturing, drying shrinking cracking of concrete.
Here, the main challenge is the estimation of the displacements, the damage field, and temperature. In [42,70], a fully coupled phase-field approach for solving thermomechanical systems was presented. This technique overcame the computational problems related to realization of sharp crack discontinuities, specifically in complex crack topologies. The authors in [71], introduced a parallel algorithm to simulate the dynamic and quasi-static brittle fracture of thermoelastic materials including source codes. Moreover, in order to enhance the efficiency of the monolithic solver, in [70], a monolithic BFGS algorithm was implemented to solve the multi-field discretized equations.

State Variables and Kinematics
In this example, the given BVP is a coupled multi-field system for fracturing solids with thermal effects. Such a model can be formulated based on a coupled three-field system. In particular, at material points x 2 B and time t 2 T, the BVP solution provides the displacement field u : B Â T ! R, the temperature field T : B Â T ! R, and the crack phase-field d : B Â T ! ½0; 1. Adopting the small strain hypothesis, the strain tensor is additively decomposed into an elastic part e e and a thermal part e T , that is, For materials that expand isotropically when the temperature is increased, the thermal part can be computed as where T 0 is the initial temperature and c is the thermal expansion coefficient.
In the present one-dimensional problem, the thermal Dirichlet and Neumann boundary conditions are described by where Q is the thermal volume flux and q is the thermal transport on the Neumann boundary, where n ¼ 1 for x ¼ L and n ¼ À1 for x ¼ 0, and T is the imposed temperature on the Dirichlet boundary. Moreover, the heat flux in (66) can be described by the negative direction of the thermal heat gradient T 0 through Fourier's law of heat conduction as Here, K d ¼ gðdÞK is the degraded thermal conductivity with the same quadratic degradation function given in (55).
To elaborate the corresponding variational formulation, we employ the function spaces (52) for the displacement and the crack phase-field, and further consider the function space (24) for the temperature field.
Moreover, the energy balance statement based on the evolution of temperature yields the following strong form initialized with Additionally, the global evolution equation for the crack phase-field in the rate-dependent form yields the following local statement: together with a homogeneous Neumann condition for the crack surface: Here, we introduce a positive crack driving force D in t 2 ½0; t n . To define D, in agreement with (45) and (54), the elastic strain energy affected by temperature reads W elas ðu; e T ; dÞ : such that a positive crack driving force D can be defined as where g(d) is a degradation function. The weak forms of the local statements (68), (69), and (71) induce the following formulation: q; E; G c ; l f ; K; q; CÞ be given with the initial conditions u 0 ¼ uðx; 0Þ, T 0 ¼ T ðx; 0Þ, and d 0 ¼ dðx; 0Þ. For the loading increments n ¼ 1; 2; . . .; N, find u :¼ u n 2 W u u , T :¼ T n 2 W T T , and d :¼ d n 2 W d d n satisfying the following Euler-Lagrange equations in weak form: The unknown material parameters for phase-field fracture in thermoelastic materials denoted as F6 then follow from the probabilistic Bayesian inversion P BI À F6ðuÞ Á with u :¼ ðu; T ; dÞ and P :¼ ðE; G c ; K; cÞ:

Bayesian Inference
The range of the parameters and the true values are listed in Table 7. The script named 1D thermoelastic:m with EnKF is used to identify the unknown parameters. We use the load-displacement curve (estimated by the true values) as the reference observation for a time interval T ¼ 1 Â 10 5 seconds.
Consider a bar of unitary length L ¼ 1 with x 2 B :¼ ½0; 1 that is initially unstretched and undamaged. Its left end is fixed, i.e., x ¼ 0, while on its right end, i.e., x ¼ 1, a monotonic displacement increment D u ¼ 1 Â 10 À3 mm is applied for 155 time steps; see Fig. 4c. We use 99 elements for the spatial discretization. The initial values are G 0 c ¼ 0:05, E 0 ¼ 300, K 0 ¼ 5 Â 10 À2 , and c 0 ¼ 1 Â 10 À3 . For temperature boundary condition, the one dimensional bar is restricted according to (25) such that T 1 ¼ 1 C and T 2 ¼ 0 C. The reference simulation is performed for t ¼ 1 Â 10 5 s.
The evolutions of the MCMC chains for all parameters are shown in Fig. 16. The unknowns show a negligible uncertainty and converge rapidly to the true values with an acceptance rate of 71.6%. The posterior densities are illustrated in Fig. 17. The distribution of G c and E is symmetric, and the true values are in the center (the median). The medians are shown in Table 7. We compare the load-displacement diagram using inferred and true values, as shown in Fig. 18. An excellent agreement between the curves points out the accuracy of the method.
As previously mentioned, the reference observation is chosen at the last time step. We show the displacement u, the phase-field d, and the temperature T for 100 time steps in Fig. 19.

Example 7: Phase-Field Approach to Predict Fatigue in Brittle Materials
The last example considers the phase-field approach to predict fatigue, a topic that has gained attention in the recent literature. Most works on this subject focus on brittle materials [43,44,[72][73][74], although models for rubber [75], hydrogen-assisted fatigue [76], ductile materials [45], and ductile materials with plastic coupling [46,49] have also been proposed. Here, we focus on the 1D formulation for the brittle case, initially addressed in [43].
The brittle fracture problem presented in Sect. 3.4 is taken as a point of departure. However, the displacements are now imposed cyclically. To account for the effect of cyclic loading, a fatigue degradation function is introduced, characterized in the simplest case by two material parameters which are hereby estimated along with Young's modulus and the fracture toughness.

State Variables and Kinematics
As in Sect. 3.4, the primary fields in the present model consist of the displacement field u : B Â T ! R and the non-local phase-field fracture variable d : B Â T ! ½0; 1, with the corresponding function spaces given in (52). Consequently, the fatigue mechanism does not entail the introduction of a new internal variable, but rather, of a history-dependent parameter, as will become clear below.
subsubsectionEnergy Quantities and Variational Principles The standard phase-field model for brittle fracture, i.e., based on the regularized crack surface energy density (56), Fig. 16 Example 6. The evolution of the MCMC chains using EnKF for four different material parameters in thermoelastic example predicts crack propagation at critical stress amplitudes. As such, the model is not suitable to describe fatigue fracture under cyclic loading, which may occur at repetitive load amplitudes well below the ultimate strength of the material. To account for fracture under such conditions, Alessi et al. [43] proposed to progressively degrade the crack resisting force by means of a suitable accumulation variable, for which a measure of accumulated energy was later established [44] and adopted in subsequent works.
In the context of the above-mentioned framework, the definition of energy quantities that comply with the second law of thermodynamics, i.e., non-decreasing energy dissipation, requires a delicate treatment. Specifically, the definition of a state function representing the regularized dissipated energy due to fracture, as done in (56) and (62), is no longer possible. Instead, we must rely on a pathdependent quantity, with a fatigue degradation mechanism acting on the rate of crack-surface energy dissipation. In particular, we may define a total energy density of the form Wðe; d; cÞ :¼ W elas ðe; dÞ þ Because the crack driving force is not modified, the stored elastic strain energy (54) is again employed in the present model. On the other hand, the time derivative of the regularized fracture energy density (56) is scaled by a nonnegative degradation function f ðcÞ, to be defined below, which renders the dissipated energy a path-dependent quantity. Note that in (75), the second law of thermodynamics is satisfied by the non-negativity of the dissipation rate, i.e., f ðcÞ d dt W frac ðd; d 0 Þ ! 0. In (75), we have introduced a fatigue accumulation variable c : B Â T ! R þ and a fatigue degradation function c7 !f ðcÞ 2 ½0; 1, endowed with the properties f ðc c c Þ ¼ 1, f ðc [ c c Þ 2 ½0; 1 and f 0 ðdÞ ! 0. In particular, we adopt the following form [44]: where k [ 0 is a material parameter that controls the rate of (logarithmic) decay of the fatigue degradation function, and c c ! 0 is a material threshold parameter, under which no fatigue effects are triggered. The fatigue variable c is defined as an accumulation of strain energy density during loading stages.
where H is the Heaviside function. In view of the path-dependent energy density (75), a total energy functional can only be approximated in incremental form. In agreement with previous incremental treatments for state-dependent dissipation potentials [77,78], we consider the form where the current time step is denoted with a subscript n þ 1. Note that the pseudo-time integral in (75) has been evaluated in incremental form by fixing the fatigue variable c at the previous time step. Minimization of the energy functional (78) with respect to the displacement field and the crack phase-field leads to the following weak formulation.
The unknown material parameters for phase-field approach to fatigue denoted as F7 then follow from the probabilistic Bayesian inversion P BI À F7ðuÞ Á with u :¼ ðu; dÞ and P :¼ ðE; G c ; c c ; kÞ: Once more, the history field approach [20] is employed here to handle the inequality constraint in (79) 2 .

Bayesian Inference
In this numerical example, we consider a bar of unitary length L ¼ 1 with x 2 B :¼ ½0; 1 that is initially unstretched and undamaged. We set c ¼ 0 at t ¼ 0. The left end of the bar is fixed, i.e., x ¼ 0, while cyclic displacements are imposed on the right end , i.e., x ¼ 1, in increments D u ¼ 2 Â 10 À4 mm. We consider 30 loading cycles with maximum and minimum values of 5 mm and -5 mm, respectively; see Fig. 4a. Regarding the spatial discretization, we use 300 elements. The cyclic displacement loading applied to the right-edge of the bar is shown in Fig. 20.
In order to perform the Bayesian inversion, we estimate simultaneously ðE; k; G c ; c c Þ (four-dimensional Bayesian inversion). The assumed ranges and true values are given in Table 8. The main script is 1D fatigue:m and we employ the DRAM algorithm to infer the materials unknowns. Figure 21 shows the effect of different parameters on the load-displacement diagram. For materials with low Young's modulus (e.g., ductile materials), several load cycles are needed for fracture. However, for too stiff materials (e.g., brittle materials), fewer load cycles lead to a full fracture. Regarding higher G c , as more energy is needed for fracture, more load cycles are necessary. Hence, e.g., G c ¼ 0:04 MPa requires more time (i.e., higher cycle number), and the difference between the peak points is more pronounced. Inversely, less G c facilitates full fracture. As shown, the same pattern is observed for c c . Finally, k is a phenomenological parameter used to calibrate the fatigue degradation process once the threshold value c c is overcome. Figure 22 shows the homogeneous response for a uniaxial test under the displacement loading displayed in Fig.4c. In order to enhance the Bayesian inference efficiency, we only select some critical points. In fact, as Fig. 21 depicts, the effect of different material unknowns on the load-displacement is more tangible on the peak points and after the fracture. Therefore, instead of all points, only these points are used for the Bayesian framework. These points are illustrated in Fig. 22. The evolutions of the MCMC chains for all parameters are shown in Fig. 23.
We use the DRAM algorithm to identify the parameters. Fig. 24 shows the MCMC chains of the parameters where the acceptance rate is 3.3% and all parameters converged. The initial values are E 0 ¼ 60 000, k 0 ¼ 0:4, G 0 c ¼ 0:3, and c 0 c ¼ 0:0075. The respective posterior densities are shown in Fig. 24, and the medians are summarized in Table 8. It seems that E and k are inferred correctly. However, the estimation of G c and c c is not sufficiently accurate. This probably due to the similar effect of the energy parameters on the load-displacement curve, e.g., causing an overestimation of G c and an underestimation of c c that may balanced each other. We use the inferred values (Table 8)    with initial values (h 0 ) and the reference values. As Fig. 25 shows, excluding the last peak point, all points are estimated acceptably, showing the relative accuracy of the model.

Overview Software Package
Bayesian inversion is a robust, simple, and efficient computational technique to provide reliable information about the model parameters and uncertainties. The mechanical models given in this work constitute relevant model problems that can be further considered in further works. The provided codes can be introduced/discussed in academic classes and can be extended to two-or three-dimensional settings. Moreover, the given examples have been considered in recent fundamental and applied studies. Of course, there are still several open questions in these fields, and we hope that the codes provided herein can be a starting point for further researches In this section, we present an overview of the codes and discuss their performance as the forward and backward models. Here, we strive to provide a comprehensive package to review the common MCMC techniques.   Table 9 shows the MCMC codes including their description. We should note that in all examples the function forwardmodel.m is the forward model and considers the response of the model problem to the given candidates. An overview of the mechanical codes is given in Table 10.
We observed that for the problems with one unknown, all MCMC techniques worked efficiently. Similar results have been achieved for a two-dimensional probabilistic problem, in case that the parameters do not correlate (heat problem). For more complicated examples, MH, AMH, and DR were not effective; therefore, DRAM and EnKF have been employed. Due to the more efficient proposal adaptation when using the Kalman filter, the acceptance rate is noticeably higher compared to DRAM. EnKF works well for phase-field problems, including brittle and ductile fracture. The obtained results for fatigue fracture were not as good as for other examples, likely due to the high correlation of the critical fracture and fatigue parameters.
Additionally, to evaluate the stability of our implementation, let us consider Example 3.5, i.e., phase-field fracture for gradient-extended ductile materials, and investigate the convergence behavior of the model problem. More specifically, we consider the residual of the staggered solution procedure. Therefore, we assume the solutions at the converged state as U k :¼ ðu k ; e p;k ; a k ; d k Þ at iteration k, such that two options for checking the staggered residual can be defined: 1. Balance of equilibrium equations (weak formulation) at the converged state: Res k Stag :¼ jE u ðU k ; duÞj þ jE a ðU k ; de p ; daÞj þ jE d ðU k ; ddÞj TOL Stag ; 2. Subsequent solutions at the converged state: Res k Stag :¼ ju k À u kÀ1 j þ ja k À a kÀ1 j þ jd k À d kÀ1 j TOL Stag : We then define the staggered residual in logarithmic space as residual :¼ log 10 ðRes Stag Þ: ð82Þ We examine the ductile phase-field fracture model through the posterior density of the unknowns obtained with EnKF in Table 6. We set TOL Stag ¼ 1 Â 10 À6 . The convergence performance based on (80) is shown in Fig. 26 (left), while the alternative approximation based on (81) is shown in Fig. 26 (right). The results show that based on the user criteria, the code reaches the desired convergence state, thus suggesting a stable implementation.

Conclusion
The implementation of Bayesian inversion in different benchmark mechanical problems has been studied in this work. We compared the performance of different common MCMC techniques. Random-walk MCMC and MH can be easily used; however, they could not infer the unknowns in   (9), where the code is taken from [4] prior.m The generation of the uniform prior distribution

EES.m
The computation of square error (4) multi-dimensional cases. Using proposal adaptation and delayed rejection (global and local adaptation) improved the Bayesian inversion. Combing both methods gives rise to an effective MCMC technique. The Kalman filter in Bayesian inference leads to an accurate estimation of the desired parameters. As observed, different material parameters can be estimated simultaneously using the probabilistic MCMC techniques. We presented a MATLAB package including all codes alongside the manuscript.
In mechanical problems, the simultaneous and accurate inference of the parameters is crucial; otherwise, the estimation is not reliable. We conclude that, despite their higher complexity, DRAM and EnKF are more suitable MCMC techniques for parameter estimation in the present examples when compared to the other presented methods. Moreover, EnKF leads to less fluctuation of MCMC in comparison with DRAM and a higher acceptance rate. Therefore, the posterior densities are relatively narrower, indicating more reliable results. Funding Open Access funding enabled and organized by Projekt DEAL.

Declarations
Conflict of interest The authors declare no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.