Multiapp with keep_solution_during_restore=true problem with history dependent subapp materials

[lynnmunday]

Check these boxes if you have followed the posting rules.

• Q&A General is the most appropriate section for my question
• I have consulted the posting Guidelines on the Discussions front page
• I have searched the Discussions forum and my question has not been asked before
• I have searched the MOOSE website and the documentation does not answer my question
• I have formatted my post following the posting guidelines (screenshots as a last resort, triple back quotes around pasted text)

Question

I made the below input files as a minimum working example to show the problem I'm having with a multiapp and keep_solution_during_restore=true with picard iterations. My main app is a transient simulation that calls a full solve multiapp that is doing a thermal mechanical transient solve. The thermomechanical solve has history dependent material properties (incremental deformation). I run a transient for the thermomechanical solve until it reaches equilibrium. During the picard solves, I would like to use keep_solution_during_restore=true so I will be closer to the equilibration of the subapp. However, the displacement keeps growing by the same amount for every picard solve, see the below main app output for the postprocessor.

Main app solve:
 0 Nonlinear |R| = 0.000000e+00
 Solve Converged!
 1 fixed point disp_x_max = 5.27730465e-03
 2 fixed point disp_x_max = 1.05686819e-02
 3 fixed point disp_x_max = 1.58743749e-02
 4 fixed point disp_x_max = 2.11946292e-02

. use keep_solution_during_restore=false gives me the correct behavior. I didn't really look into the issue and using the new mechanics kernel system solves the problem but I wanted to see if anyone knows the problem or if I should create an issue. Maybe it should error out if you have history dependent material properties.

main.i:

[Mesh]
  [generated]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 10
    ny = 10
    xmax = 2
    ymax = 1
  []
  [pin]
    type = ExtraNodesetGenerator
    input = generated
    new_boundary = pin
    coord = '0 0 0'
  []
[]

[Variables]
  [diffused]
    order = FIRST
    family = LAGRANGE
  []
[]
[Kernels]
  [diff]
    type = Diffusion
    variable = diffused
  []
[]
[AuxVariables]
  [disp_x]
    order = FIRST
    family = LAGRANGE
    initial_condition = 0.0
  []
[]
[Executioner]
  type = Steady
  solve_type = 'PJFNK'
  fixed_point_algorithm = picard
  custom_pp = disp_x_max
  custom_abs_tol = 1e-6
  fixed_point_min_its = 4
  fixed_point_max_its = 5
[]

[MultiApps]
  [sub_app]
    type = FullSolveMultiApp
    input_files = 'subapp.i'
    execute_on = 'timestep_begin'
    keep_solution_during_restore = true # uses previous iterate as initial guess for the next iteration
  []
[]

[Transfers]
  [disp_z_avg_clad]
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = sub_app
    source_variable = 'disp_x'
    variable = disp_x
    displaced_source_mesh = false
    displaced_target_mesh = false
    search_value_conflicts = false
    greedy_search = true
  []
[]
[Postprocessors]
  [disp_x_max]
    type = ElementExtremeValue
    variable = 'disp_x'
    value_type = max
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]
[Outputs]
  csv = true
[]

subapp.i:

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  [generated]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 10
    ny = 10
    xmax = 2
    ymax = 1
  []
  [pin]
    type = ExtraNodesetGenerator
    input = generated
    new_boundary = pin
    coord = '0 0 0'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [all]
    add_variables = true
    strain = FINITE
    incremental=true
  []
[]

[Materials]
  [elasticity]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1e9
    poissons_ratio = 0.3
  []
  [stress]
    type = ComputeFiniteStrainElasticStress
  []
[]

[BCs]
  [pull_x]
    type = NeumannBC
    variable = disp_x
    boundary = right
    value = 1e6
  []
  [pin_x]
    type = DirichletBC
    variable = disp_x
    boundary = pin
    value = 0
  []
  [bottom_y]
    type = DirichletBC
    variable = disp_y
    boundary = bottom
    value = 0
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  num_steps =1
  dt = 1
[]

[Outputs]
[]

[lynnmunday]

@bwspenc @dschwen @hugary1995
Someone asked me about this yesterday for a coupled griffin bison model.

[hugary1995]

The issue is here https://github.com/idaholab/moose/blob/next/modules/solid_mechanics/src/materials/ComputeFiniteStrain.C#L57-L59

It calculates old deformation gradient using old displacement gradients. This saves us some memory apparently.

[hugary1995]

Though to be fair this is an issue with keep_solution_during_restore = true. Even if we keep the solution, we should still restore the old solution.

[hugary1995]

I just read the source code for keep_solution_during_restore, and honestly I can't see anything wrong with it... Perhaps somewhere else we are calling advanceState()? You could add two print statements, one inside MultiApp::restore() after line 848, the other one inside copyOldSolutions() to first confirm whether we are somehow copying solutions into old solutions after we manually "unrestored" the current solution that should have been restored.

[bwspenc]

I just looked at this a bit with Lynn. He's going to look at what's going on with that old state. It seems like regardless of how you have keep_solution_during_restore set it shouldn't be changing the old state of the variables. While we could adjust that code in our material model to make it work, it seems like we should fix keep_solution_during_restore.

Now that I think about this, this rings a bell about another issue we were seeing with Picard iterations with xfem. The Picard iterations were being treated as new time steps, so timestepSetup() was getting called on things when we didn't anticipate it. Maybe the state is getting advanced somewhere. Probably the easiest thing to do is print out the new and old displacement gradients and infer from that what is happening.

[hugary1995]

Did you try the two print statements I suggested above? That way we can at least first confirm this is caused by old solution being copied back. That'll also serve as a good break point later when we run it through a debugger.

[hugary1995]

Okay, I just did that and can confirm we are indeed copying old solutions after unrestoring the current solution. So that's the root cause. Let me step through the debugger to figure out how to fix this.

[GiudGiud]

Taking a look.
The problem seems to come from interference between FEProblemBase::advanceState() and MultiApp::restore as Gary said

[GiudGiud]

can you guys try using a TransientMultiApp instead?

[MultiApps]
  [sub_app]
    type = TransientMultiApp
    input_files = 'subapp.i'
    execute_on = 'timestep_begin'
    keep_solution_during_restore = true # uses previous iterate as initial guess for the next iteration
    catch_up = true
  []
[]

there is custom logic in trnasientMultiapp for this case.

Note in the subapp I also did this to get it to converge and accepts steps after the first

[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  num_steps =1
  start_time = -1
  dt = 1
  nl_abs_tol = 1e-6
[]