How to get only the converged solution for a material property?

[jaYINGLING]

I am writing a material property that updates the elasticity as a function of the highest vonmises_stress in each element.

In the computeQpElasticityTensor() method I have the following logic:

if (_trial_stress[_qp] > _max_trial_stress[_qp])
     _max_trial_stress[_qp] = _trial_stress[_qp];

After which a damage coefficient is calculated to determine the young's modulus and poisson's ratio. _trial_stress is a coupledVar to vonmises_stress and _max_trial_stress is a declared material property.

For some reason using the logic above to specify _max_trial_stress results in some elements getting assigned values of vonmises_stress that were never seen in any previous timestep for a particular element. More specifically, if I use paraview to track the values of vonmises_stress in a particular element, it may never exceed 20MPa. However, the max_trial_stress for that element might be as high as 60MPa.

I have verified that I can correctly assign _max_trial_stress[_qp] = _trial_stress[_qp]; to each element by removing the if statement. However, this logic is essential to describing the elasticity of the material upon relaxation.

It seems like maybe very large values of vonmises_stress are being stored in _max_trial_stress during the iteration process before a solution for the timestep is found. I'm not sure how to get around this problem. I tried using coupledVarOld to use the previously converged solution, but the problem still occurs.

At what point should I assign _max_trial_stress[_qp]=_trial_stress[_qp] so that it only takes the value of the converged solution?

Perhaps a related question is: How is _qp used to write to material properties. Does each thread iterate through a list of _qps or is the material property accessed all at once as a vector?

[aeslaughter]

Material properties are computed on the fly and don't maintain their state like you are expecting. The the _max_trial_stress references memory in the internals of MOOSE that get swapped out at each element before calling the compute methods.

This problem might require the use of old material properties (getMaterialPropertyOld) to use for comparison. Hopefully, someone with more experience with this type of calculation will jump in.

[jaYINGLING]

OK. A solution might be to get the _trial_stress and assign it to a material property _trial_stress_prop. Then use:

if (_trial_stress_prop_old[_qp] > _max_trial_stress[_qp])
     _max_trial_stress[_qp] = _trial_stress_prop_old[_qp];

I'll give this a try. Thank you.

[jaYINGLING]

Unfortunately this does not solve the problem. Here you can see my material that gets strained from left to right. The _max_trial_stress is the same for large groups of elements with sudden transitions that I think are related to the number of threads used to do the solve. In effect, using just one thread would result in all mesh elements taking the same value of _max_trial_stress. The correct behavior would be to see a gradient of _max_trial_stress, largest to smallest, from left to right.

[jaYINGLING]

Perhaps I would be better off writing an AuxKernel that computes _max_trial_stress and then use that to calculate the elasticity?

[aeslaughter]

I wouldn't recommend this there is potential to get interpolation errors because aux variables are computed at differing locations than the material properties. Additionally, automatic differentiation is not supported.

[aeslaughter]

@sapitts @bwspenc Do you have any ideas for solving this problem?

[sapitts]

It might be worthwhile to put the _max_trial_stress[_qp] setting/ calculation in a computeStressFinalize method; however, we only have that type of method within the radial return type methods instead of the ComputeElasticityTensor classes.

You might instead consider storing the maximum von mises stress value (_vonmises_stress_old) from the previous timestep in the first iteration of the current timestep. With this approach your max_trial_stress value would lag by a timestep, but you would know you were working with the converged von mises stress value from that previous timestep. @dschwen helped me put in logic for a different material calculation, which essentially had the form

int number_nonlinear_it =
      _fe_problem.getNonlinearSystemBase().getCurrentNonlinearIterationNumber();
  if (number_nonlinear_it == 0)

in the jacobianSetup() method of my material class. Perhaps this approach would be helpful for your problem too?

[bwspenc]

Sorry, I'm super late to this conversation. It sounds like what you want is a damage model. Rather than updating the elasticity tensor, you would want to update a damage index that is computed based on the strain. We have a fully developed system for damage models in TensorMechanics. Look at ComputeDamageStress. That's used together with the actual damage model, which computes the scalar damage material property. I don't think we have a "real" damage model in MOOSE, but there is one in BlackBear that derives from the ScalarDamageBase class in MOOSE: https://github.com/idaholab/blackbear/blob/devel/src/materials/MazarsDamage.C

[jaYINGLING]

These are good leads. Thank you for your help. This is a side project for me so I apologize for the delayed responses.