Skip to content

Albany Issues

Jump to bottom
Brian Granzow edited this page Apr 27, 2015 · 9 revisions

#Individual Issues

Parameter Units

There was a question of whether values were using consistent units. Brian has looked at this, it seems they are consistent.

Convergence for "Null" Creep Model

I will coin the term Null Creep to refer to the current creep model code with material parameters set such that it should be mathematically equivalent to J2. In this case, we would expect convergence rates to be the same. Andrew has confirmed that they are not (Null Creep converges linearly while J2 does not).

Further testing has confirmed that the Null Creep model will converge quadratically when below the yield strength and convergence becomes linear above the yield strength.

In-depth debugging reveals that the output of Null Creep and J2 diverge at the first Jacobian computation where the yield condition is crossed (f > 0).

Convergence depends on Jacobian

We are assuming that the different convergence rate is caused by the Jacobian values being different, likewise we are assuming that Null Creep and J2 do output the exact same residual.

This seems to be confirmed, residuals only diverge and quadratic convergence is only lost after the first Jacobian Creep model outputs diverge.

Jacobian depends on automatic differentiation

Further we assume that the Jacobian is formed entirely using the derivative values obtained by code transformation of the model where every scalar is replaced with an automatic differentiation type.

This has been totally confirmed.

Specifics of automatic differentiation

The output fields "stress", "Fp", and "eqps" are the model's contribution to the residual/Jacobian, and divergence in the derivative values for these fields causes convergence loss. As mentioned above, this occurs when the yield condition is reached, i.e. (f > epsilon).

stress

The diagonal values of the output stress had different derivatives. In both models, stress = p * I + s / J(cell,pt). The values and derivatives for p and J(cell,pt) are the same, I is the identity matrix, so the problem is traced to s. The diagonal values of s have different derivatives, and come from different formulas. In both models, there is an initial value for s and also an N = s / norm(s), which are identical. In Creep, the new value is s = smag_new * N. In J2, the new value is s = s - 2 * mubar * dgam * N. Although both formulas result in the same s values, the derivatives in Creep's smag_new are not aligned with the derivatives in J2's mubar and dgam, such that the resulting s derivatives are different.

Fp

Like stress, Fp had different diagonal derivatives. In both models, Fp = expA * N, and N is the same and expA has different derivatives. Furthermore, expA = exp(A) where A = dgam * N, so the whole issue traces to dgam. Note that dgam is also involved in the stress derivatives for J2. In this case, J2 does have derivatives for dgam while Creep does not. J2 has derivatives because it gets its dgam directly from X[0] after getting X[0] from solver.computeFadInfo (see below). Creep has zero derivatives because it computes dgam = (smag_new - smag) / (-2 * mubar), and in our Null Creep test case smag_new = smag, derivatives and all, so their derivatives cancel.

eqps

In the eqps output, the J2 model showed zero values and zero derivatives, while the Creep model showed non-zero derivatives. This seems to be because the J2 model programmers were careful to separate derivates when using the LocalNonlinearSolver, which has a method called computeFadInfo. Specifically, in J2, eqps = alpha = eqpsold(cell,pt) + sq23 * X[0], and to prevent alpha from having derivatives, X[0] needs to not have them. In the Creep model, alpha = eqpsold(cell,pt) + sq23 * dgam_plastic, and dgam_plastic has derivatives. One less than elegant solution is simply to strip the derivatives from alpha before using it further. This solution did result in identical eqps zero derivatives, but is not sufficient for convergence, stress and Fp must also be resolved.

Resources

The above details were obtained using specially modified versions of both models. Three modifications were made:

  1. a timer was inserted and a conditional check written that trips only during the first divergent computation
  2. overloaded print functions were written that print either normal scalars or value/derivative automatic differentiation types. this function is especially constructed so that a textual diff shows only significant changes in values.
  3. a type-generic function to strip derivatives was written and used on alpha in Creep.

The modified model files can be copied from here:

/lore/dibanez/trilinos/cube_2x2x2/CreepModel_Def.hpp /lore/dibanez/trilinos/cube_2x2x2/J2Model_Def.hpp

Likewise, there is an execution script and log files there:

/lore/dibanez/trilinos/cube_2x2x2/good_run4/Creep_log.txt /lore/dibanez/trilinos/cube_2x2x2/good_run4/J2_log.txt /lore/dibanez/trilinos/cube_2x2x2/run.sh

The recommended workflow is to run the execution script and then use

vimdiff J2_log.txt Creep_log.txt

Performance sensitive to parameters

Currently varying material parameters can greatly vary solution time. this will be revisited once the Null Creep case is resolved.

Condition number sensitive to Essential BCs

Dirichlet boundary conditions are treated as degrees of freedom in every Albany simulation. This means that rows in the Jacobian that are associated with a DBC degree of freedom will be replaced by a row with a one on the diagonal and zeros everywhere else. Every other row in the Jacobian matrix is generally on the order of some material parameter or parameters. For instance, for linear elasticity, the Jacobian gets multiplied through by the elastic modulus E, which is typically on the order of GPa, or 10^9 Pa. The difference between DBC rows ~O(1) and the other rows in the Jacobian ~O(E) will increase the condition number of the matrix by ~O(E).

Left versus Right precoditioning

Consider the Jacobian A formed by Albany in block form.

A = [V B;
     0 D]

where D implements the Dirichlet boundary conditions, V is associated with non-DBC degrees of freedom, and B is the coupling between B and D. We'll assume that the norms of the sub matrices are of the following order

||V|| ~ O(E)
||B|| ~ O(E)
||D|| ~ O(1)

where E >> 1.

Let's assume we are just interested in left vs. right preconditioning to perform a scaling of the matrix A, so that everything is O(1).

Left preconditioning

Consider the block matrix R such that

R = [ R1   0;
      0    R2]

If we are interested in performing scaling by multiplying the Jacobian A by the matrix R on the left, we get

R A = [R1 V   R1 B;
         0    R2 D]

||R1 V|| ~ O(1)
||R1 B|| ~ O(1)
||R2 D|| ~ O(1)

Right preconditioning

Consider the block matrix C such that

C = [ C1   0;
      0    C2]

If we are interested in trying to perform scaling by multiplying the Jacobian A by the matrix C on the right, we get

A C = [V C1    B C2;
         0     D C2]

Clearly, it's not possible for both

||B C2|| ~ O(1)
||D C2|| ~ O(1)

to be satisfied at the same time.

Plan of Action

  1. 04/23/2015 - it is suggested that we build a very small mesh with Null Creep input and print the first available values of the "F" vector in its residual and Jacobian form, comparing these values to the exact same experiment run with J2 instead of Null Creep.
  2. 04/25/2015 - extensive use of this simple 2x2x2 test case is revealing in detail where the problems are. there are a few specific variables that either have derivatives and shouldn't or vice versa. we should ask Andrew Bradley about computeFadInfo and then work out the mathematics to fix the derivatives.
  3. 04/26/2015 - Andrew Bradley has added an option for left preconditioning in Stratimikos. This has been tested using various pre conditioners (ML, MueLu, Ifpack2) and does indeed improve the conditioning of the Jacobian.
Clone this wiki locally