postprocessed variable values are zero at all timesteps, while solving a dynamic phase field fracture problem.

[zuckarm]

hi all,
i am trying to solve a dynamic phase field fracture problem.
geometry is a unit square plate, with right and left edges having roller supports. bottom edge is fixed, and a load is applied on the top edge. the load is varying with time.
i am running the files, and whatever i do the solution (postprocessed) values are zero at all timesteps, right from step one.
following are my input files:

main file:

E = 1000
nu = 0.3
K = '${fparse E/3/(1-2*nu)}'
G = '${fparse E/2/(1+nu)}'
Gc = 0.025
l = 0.1

[MultiApps]
    [fracture]
        type = TransientMultiApp
        input_files = prob3_fracture.i
        cli_args = 'Gc=${Gc};l=${l}'
        execute_on = 'TIMESTEP_END'
    []
[]

[Transfers]
    [from_d]
        type = MultiAppCopyTransfer
        from_multi_app = 'fracture'
        variable = d
        source_variable = d
    []
    [to_psie_active]
        type = MultiAppCopyTransfer
        to_multi_app = 'fracture'
        variable = psie_active
        source_variable = psie_active
    []
[]

[Mesh]
	type = FileMesh
	file = prob3.msh
[]

[GlobalParams]
  volumetric_locking_correction = false
  use_displaced_mesh = true
  displacements = 'disp_x disp_y'  
  order = FIRST
  family = LAGRANGE
[]

[Variables]
	[./disp_x]
        order = FIRST
        family = LAGRANGE
  	[../]
  	[./disp_y]
        order = FIRST
        family = LAGRANGE
  	[../]
[]

[Kernels]
  [dynamic_x]
	type = ADDynamicStressDivergenceTensors
    variable = disp_x
    component = 0
	zeta = 0
  []
  [dynamic_y]
	type = ADDynamicStressDivergenceTensors
    variable = disp_y
    component = 1
	zeta = 0
  []
  [inertia_x]
    type = ADInertialForce
    variable = disp_x
	density = density
	eta = 0
  []
  [inertia_y]
    type = ADInertialForce
    variable = disp_y
	density = density
	eta = 0
  []
[]

[AuxVariables] 
    [fx]
    []
    [fy]
    []
    [d]
    []
  	[./accel_x]
  	[../]
  	[./vel_x]
  	[../]
  	[./accel_y]
  	[../]
  	[./vel_y]
  	[../]
[]

[AuxKernels]
	[./acceleration_x]
    type = TestNewmarkTI
    variable = accel_x
    displacement = disp_x
    first = false
  [../]
  [./velocity_x]
    type = TestNewmarkTI
    variable = vel_x
    displacement = disp_x
  [../]
  [./acceleration_y]
    type = TestNewmarkTI
    variable = accel_y
    displacement = disp_y
    first = false
  [../]
  [./velocity_y]
    type = TestNewmarkTI
    variable = vel_y
    displacement = disp_y
  [../]
[]		

[Functions]
	[./load]
		type = PiecewiseLinear
		x = '0 	0.05 0.1'
		y = '0 	1    1'
	[../]
[]
  
[BCs]
	[./roller_left]
		type = ADDirichletBC
		variable = disp_x
		boundary = left
		value = 0
	[../]
    [./roller_right]
		type = ADDirichletBC
		variable = disp_x
		boundary = right
		value = 0
	[../]			
	[.fix_y]
		type = ADDirichletBC
		variable = disp_y
		boundary = bottom
		value = 0
	[../]
	[.fix_x]
		type = ADDirichletBC
		variable = disp_x
		boundary = bottom
		value = 0
	[../]
	[./traction]
		type = ADPressure
		variable = 'disp_y'
		boundary = 'top'
		function = load
		factor = -1
	[../]
[]

[Materials]
    [bulk]
        type = ADGenericConstantMaterial
        prop_names = 'K G'
        prop_values = '${K} ${G}'
    []
    [./density]
		type = ADGenericConstantMaterial
		prop_names = density
		prop_values = 1
	[../]
    [degradation]
        type = PowerDegradationFunction
        f_name = g
        function = (1-d)^p*(1-eta)+eta
        phase_field = d
        parameter_names = 'p eta'
        parameter_values = '2 0'
    []
    [strain]
        type = ADComputeSmallStrain
    []
    [elasticity]
        type = SmallDeformationIsotropicElasticity
        bulk_modulus = K
        shear_modulus = G
        phase_field = d
        degradation_function = g
        decomposition = NONE
        output_properties = 'elastic_strain psie_active'
        outputs = exodus
    []
    [stress]
        type = ComputeSmallDeformationStress
        elasticity_model = elasticity
        output_properties = 'stress'
        outputs = exodus
    []
[]

[Executioner]
	type = Transient
    solve_type = NEWTON
    petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
    petsc_options_value = 'lu       superlu_dist                 '
    automatic_scaling = true
    nl_rel_tol = 1e-8
    nl_abs_tol = 1e-10
    start_time = 0
    dt = 0.005 
    end_time = 0.1
    fixed_point_max_its = 10
    fixed_point_min_its = 2
    accept_on_max_fixed_point_iteration = true
    fixed_point_rel_tol = 1e-8
    fixed_point_abs_tol = 1e-10
         [TimeIntegrator]
            type = NewmarkBeta
            beta = 0.25
            gamma = 0.5
        []
[]

[Postprocessors]
	[./y_disp]
		type = NormalBoundaryDisplacement
    	#displacements = disp_y
    	boundary = top
    	value_type = absolute_max
  	[../]
    [./phi]
        type = AverageNodalVariableValue
        variable = d
    [../]
[]

[VectorPostprocessors]
	[./y_disp_and_d]
		type = LineValueSampler
		start_point = '0.5 0 0'
		end_point = '0.5 1 0'
		num_points = 51
		variable = 'disp_y d'
		sort_by = y
		contains_complete_history = true
	[../]
[]

fracture file:

[Mesh]
  type = FileMesh
  file = prob3.msh
[]

[Variables]
  [d]
  []
[]

[AuxVariables]
  [bounds_dummy]
  []
  [psie_active]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Bounds]
  [irreversibility]
    type = VariableOldValueBoundsAux
    variable = bounds_dummy
    bounded_variable = d
    bound_type = lower
  []
  [upper]
    type = ConstantBoundsAux
    variable = bounds_dummy
    bounded_variable = d
    bound_type = upper
    bound_value = 1
  []
[]

[Kernels]
  [diff]
    type = ADPFFDiffusion
    variable = d
    fracture_toughness = Gc
    regularization_length = l
    normalization_constant = c0
  []
  [source]
    type = ADPFFSource
    variable = d
    free_energy = psi
  []
[]

[Materials]
  [fracture_properties]
    type = ADGenericConstantMaterial
    prop_names = 'Gc l'
    prop_values = '${Gc} ${l}'
  []
  [./density]
		type = ADGenericConstantMaterial
		prop_names = density
		prop_values = 1
	[../]

  [degradation]
    type = PowerDegradationFunction
    f_name = g
    function = (1-d)^p*(1-eta)+eta
    phase_field = d
    parameter_names = 'p eta'
    parameter_values = '2 0'
  []
  [crack_geometric]
    type = CrackGeometricFunction
    f_name = alpha
    function = 'd^2'
    phase_field = d
  []
  [psi]
    type = ADDerivativeParsedMaterial
    f_name = psi
    function = 'alpha*Gc/c0/l+g*psie_active'
    args = 'd psie_active'
    material_property_names = 'alpha(d) g(d) Gc c0 l'
    derivative_order = 1
  []
[]

[Executioner]
  type = Transient

  solve_type = NEWTON
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -snes_type'
  petsc_options_value = 'lu       superlu_dist                  vinewtonrsls'
  automatic_scaling = true

  nl_rel_tol = 1e-8
  nl_abs_tol = 1e-10
[]

[Outputs]
  print_linear_residuals = true
[]

i have ensured that the length scale is 5 times the mesh element size.

what have i missed?

[GiudGiud]

hello

which postprocesing is 0 at all timesteps?
the postprcoessors or the vectorpostprocessor?

Is the sub-app giving reasonable results on its own?

Guillaume

[GiudGiud]

you can try

• running the subapp input only
  or
• adding outputs to the subapp and checking those

[zuckarm]

i added the following block to my sub app input file:

[Postprocessors]
	[./phi_sub_app]
        type = AverageNodalVariableValue
        variable = d
    [../]
[]

still, it's the same issue. postprocessed values are zero at all timesteps.

[GiudGiud]

ok. And you're confident phi_sub_app should not amount to 0? At the end of the time steps you are examining.
Are you looking at the value in the console or in an output csv file?

[zuckarm]

yes, since i am applying a load on top edge, there should be some evolution of phase field parameter, and there should be gradual increase in displacement too as the load ramps with time. i am looking at both console and csv outputs.

i think there is a problem with the time integrator block. for some reason there is no communication between the main and sub apps. i am not sure what is the reason for this.

[GiudGiud]

you can use Problem/verbose_multiapps = true to see the transfers happen.

The only reasons I could see for this are:

• the solve is failing completely and the fields are not updating. you should see this in the log
• you are looking at the wrong CSV. Either an old one or the one output on initial conditions
• something is overwriting d in the main app, I dont see an auxkernel doing that
• a bug in the VPP, likely because of "contains_complete_history" which is not always implemented