Initialize old Material Property

[delpierrena]

Hello everyone,

I am currently trying to implement the particle finite element method using the Moose FEM solver.
The application could be landslide, earthen dam failures, footing etc.. The basic idea of the method is to remesh when necessary to allow for modelling large deformations problems. The remeshing algorithm is working fine but I am struggling harder with reinitializing the moose simulation from a past solution.
Since our problem involves large deformations, we need to remap a lot of quantities from one mesh to another. So we run a simulation, retrieve the plastic strain, internal parameter, stresses, velocities and accelerations and we want to remap these into the new displaced mesh.
The issue we face is that we are unable to initialize the old material properties before starting the simulation. We modified the initQpStatefulProperties function from the ComputeMultiPlasticityStress class in order to initialize the plastic_strain, the stress and the intnl parameter but this is obviously not the way to go.
Example : We performed a slope failure simulation. We output some quantities at a given point in the slope after 0.1s of simulation. Here the elastic strain xx (es_xx) the plastic strain xx (ps_xx) the strain xx (s_xx) and the stress yy.

Postprocessor Values:
+----------------+----------------+----------------+----------------+----------------+
| time           | es_xx          | ps_xx          | s_xx           | stress_yy      |
+----------------+----------------+----------------+----------------+----------------+
|   0.000000e+00 |   0.000000e+00 |   0.000000e+00 |   0.000000e+00 |   0.000000e+00 |
|   1.000000e-02 |  -1.467515e-05 |   1.252200e-05 |  -2.153149e-06 |  -3.078785e+04 |
|   2.000000e-02 |  -1.678404e-05 |   1.252200e-05 |  -4.262041e-06 |  -3.115457e+04 |
|   3.000000e-02 |  -4.448685e-06 |   1.252200e-05 |   8.073315e-06 |  -3.157687e+04 |
|   4.000000e-02 |   2.029382e-05 |   1.252200e-05 |   3.281582e-05 |  -3.271324e+04 |
|   5.000000e-02 |   4.244592e-05 |   1.252200e-05 |   5.496792e-05 |  -3.546752e+04 |
|   6.000000e-02 |   7.036422e-05 |   1.252200e-05 |   8.288622e-05 |  -3.904428e+04 |
|   7.000000e-02 |   1.045518e-04 |   1.252200e-05 |   1.170738e-04 |  -4.171049e+04 |
|   8.000000e-02 |   1.193600e-04 |   1.252200e-05 |   1.318820e-04 |  -4.339706e+04 |
|   9.000000e-02 |   1.225850e-04 |   2.543216e-05 |   1.480172e-04 |  -4.536374e+04 |
|   1.000000e-01 |   1.144755e-04 |   8.649673e-05 |   2.009722e-04 |  -4.688521e+04 |
+----------------+----------------+----------------+----------------+----------------+

Then, we want to restart the simulation based on our modified version of InitQpStatefulProperties but the results are clearly wrong since the strain is not the sum of the elastic and of the plastic contribution and this, since the first step :
Postprocessor Values: +----------------+----------------+----------------+----------------+----------------+ | time | es_xx | ps_xx | s_xx | stress_yy | +----------------+----------------+----------------+----------------+----------------+ | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | | 1.000000e-02 | 3.105196e-05 | 1.682521e-04 | 2.885585e-05 | -3.673527e+04 | +----------------+----------------+----------------+----------------+----------------+

We believe that the issue is coming from the reinitialization of some fields (most probably the plastic strain ). In particular, when reinitializing the finite element problem, we cannot initialize the old properties needed for the return map algorithm. Do you have any idea of how we could initialize these old material properties ?

[GiudGiud]

Hello

Things have improved recently on this front, with the new stateful material properties objects merged this year

see https://mooseframework.inl.gov/source/actions/ProjectedStatefulMaterialStorageAction.html

Did you try to rely on those for your work?
How are you reinitializing the problem? Is this using restart/recover and checkpoints? Or exodus? Or some special problem?

[maxnezdyur]

I don't know off hand what you should do. One thing to consider is that you need an Auxkernel for each component of the vector but I don't know if the storage action knows how many components there are where the action happens. I feel that's the reason it is not on the supported types.

[GiudGiud]

it should use an array auxkernel most likely

[GiudGiud]

So in addition of putting std::vector in the TypeList

did you manage to do that? Seems to be some pretty difficult template meta programming
I think this will be easier, though it s not finished
https://github.com/GiudGiud/moose/tree/PR_vector_projected

[GiudGiud]

I just need to add a few more objects for interpolated stateful properties. I ll work on this tomorrow

[delpierrena]

I was trying to do it but indeed the templating is pretty difficult to understand from scratch

did you manage to do that? Seems to be some pretty difficult template meta programming

[delpierrena]

Just to clarify the full process, I show you different important blocks of my input file:
I declared all the different projected properties in the following block :

[ProjectedStatefulMaterialStorage]
   [all]
     projected_props = 'plastic_strain_interpolated_old stress_interpolated_old total_strain_interpolated_old elastic_strain_interpolated_old plastic_internal_parameter_interpolated_old plastic_transverse_direction_interpolated_old'
     family = MONOMIAL
     order = FIRST
   []
 []

I have a block of function where I retrieve the old state of a previous simulation, or let's say the initial state of the current one:

[Functions]
  [stress_xx_fcn]
    type = PiecewiseConstantFromCSV
    read_type = 'voronoi'
    read_prop_user_object = 'reader_node'
    column_number = '3'
  []
  [stress_xy_fcn]
    type = PiecewiseConstantFromCSV
    read_type = 'voronoi'
    read_prop_user_object = 'reader_node'
    column_number = '4'
  []
......

Then the Materials block :

[Materials]
  [plastic_strain_interpolated_old]
    type = GenericFunctionRankTwoTensor
    tensor_functions = 'plastic_strain_xx_fcn plastic_strain_xy_fcn 0 plastic_strain_xy_fcn plastic_strain_yy_fcn 0  0 0 0'
    tensor_name = plastic_strain_interpolated_old
  []
  [elastic_strain_interpolated_old]
    type = GenericFunctionRankTwoTensor
    tensor_functions = 'elastic_strain_xx_fcn elastic_strain_xy_fcn  0 elastic_strain_xy_fcn elastic_strain_yy_fcn 0  0 0 0'
    tensor_name = elastic_strain_interpolated_old
  []
  .......
   [./mc]
    type = ComputeMultiPlasticityStress
    ep_plastic_tolerance = 1E-9
    plastic_models = mc
    perform_finite_strain_rotations = false
    min_stepsize = 0.01
    ignore_failures = true
    use_interpolated_state = true
  [../]
   ...

[GiudGiud]

in the Materials block, is use_interpolated_state = true in every relevant block?

If you use Debug/show_actions=true,
do you see auxkernels updating every variable that needs to be updated?
the Functions, are they used in AuxKernels? or in FunctionIC (preferrable)

[dschwen]

Guys, if you know the number of elements to store in that vector then just change it into a statically sized std::array<Real, N>. Adding dynamically sized types to the projected material storage will be an ugly mess! There is no enforced guarantee that vector sizes stay the same during a simulation, so all the extra checks and the added parameters for users to specify vector sizes will make the code harder to read and maintain.

[GiudGiud]

I've created #28836 to keep track of this

[GiudGiud]

@delpierrena
is there any reason you use ComputeMultiPlasticityStress ? Its use of std::vector properties is because of the potentially arbitrary number of plasticity models in use.

Could you use another material?
I m going to guess but this one maybe?
https://mooseframework.inl.gov/source/materials/IsotropicPlasticityStressUpdate.html
@dschwen

[GiudGiud]

I think such an auxkernelight be added by the ProjectedStatefulMaterial property action?

[delpierrena]

Because applying or not the initial stress does not change the solution and the output variables at each time step.
I tried to do the same for the elastic_strain_old and initialize them, but it does not apply either:

 [./_var_elastic_strain_4]
    family = MONOMIAL
    order = FIRST
    [./InitialCondition]
      type = FunctionIC
      function = -0.001
    [../]
  [../]

Tomorrow I will add some cout in the code to be sure, but I believe that the solution should at least be different at time 1, influenced by the initial condition.

[maxnezdyur]

When you get a chance try to do both of these for your strain. Without the auxkernel, we don't set the initial condition correctly. The auxkernel is only run on initial. I think this goes back to something I said a while ago, we need the ability to apply ICs on old states of variables. The auxkernel allows us to do that. I tested this with strain_yy and got different results!Let me know if this works for you.

[AuxKernels]
  [harden_initial]
    type = ConstantAux
    variable = _var_hardening_variable_0
    execute_on = 'INITIAL'
    value =  -123456789
  []
  []


  [ICs]
  [harden_ic]
    type = ConstantIC
    variable = _var_hardening_variable_0
    value = -123456789
  []
[]

@GiudGiud
I think we need to add the ability for Material<T>Aux to be able to grab a different material state. It would default to 0 (current) but we could go backwards. This would allow us to check with interpolated state on and off directly.

[GiudGiud]

so this does not do it?
https://mooseframework.inl.gov/source/auxkernels/ProjectedStatefulMaterialAux.html
Reading the documentation it does things on the material property as well. but I would expect the output aux variable to be the old state?

[maxnezdyur]

I don't think it does. I think that projects the current timestep material onto a variable so that for the next timestep, we can use that "old" variable for projecting back.
But even if it did, I think an auxkernel that would get the old material the same way any other material would getMaterialPropertyOld<T>(...), would be useful for debugging in general.