Question about ConstantBounds and VariableOldValueBounds not working #27645

Yasunaoka · 2024-05-16T13:46:14Z

Yasunaoka
May 16, 2024

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

Hello MOOSE community!
I've recently been working on the phase field fracture problem, where the phase field parameters need to be constrained to boundaries of 0 and 1 and monotonically increasing. I found that ConstantBounds and VariableOldValueBounds met my requirements well, and I made my input file settings as follows:

[Bounds]
  [irreversibility]
    type = VariableOldValueBounds
    variable = bounds_dummy
    bounded_variable = c
    bound_type = lower
  []
  [upper]
    type = ConstantBounds
    variable = bounds_dummy
    bounded_variable = c
    bound_type = upper
    bound_value = 1.0
  []
[]

[Executioner]
   type = Transient
   petsc_options_iname = '-pc_type -pc_factor_shift  -snes_type'
   petsc_options_value = 'lu NONZERO vinewtonrsls'

  [TimeIntegrator]
    type = ActuallyExplicitEuler
    solve_type = lumped
  []
[]

However, from what I've observed, they don't seem to be working correctly:

Time Step 21987, time = 11.8948, dt = 0.0005
 Solve Converged!
fracture0: 
fracture0: Time Step 21987, time = 11.8948, dt = 0.0005
fracture0:  Solve Converged!

Postprocessor Values:
+----------------+----------------+----------------+----------------+----------------+----------------+
| time           | max_c          | max_von        | min_degration  | numelemqps     | time_step      |
+----------------+----------------+----------------+----------------+----------------+----------------+
:                :                :                :                :                :                :
|   1.188778e+01 |   7.880643e-01 |   9.987684e+02 |   1.086355e-01 |   1.451520e+06 |   5.000000e-04 |
|   1.188828e+01 |   7.910732e-01 |   9.985850e+02 |   1.060929e-01 |   1.451520e+06 |   5.000000e-04 |
|   1.188878e+01 |   7.944909e-01 |   1.000007e+03 |   1.032356e-01 |   1.451520e+06 |   5.000000e-04 |
|   1.188928e+01 |   7.983761e-01 |   1.001331e+03 |   1.000241e-01 |   1.451520e+06 |   5.000000e-04 |
|   1.188978e+01 |   8.027761e-01 |   1.000273e+03 |   9.643148e-02 |   1.451520e+06 |   5.000000e-04 |
|   1.189028e+01 |   8.077593e-01 |   9.986019e+02 |   9.241737e-02 |   1.451520e+06 |   5.000000e-04 |
|   1.189078e+01 |   8.134319e-01 |   9.969938e+02 |   8.791847e-02 |   1.451520e+06 |   5.000000e-04 |
|   1.189128e+01 |   8.198848e-01 |   9.980568e+02 |   8.287774e-02 |   1.451520e+06 |   5.000000e-04 |
|   1.189178e+01 |   8.268397e-01 |   1.000820e+03 |   7.748472e-02 |   1.451520e+06 |   5.000000e-04 |
|   1.189228e+01 |   8.329732e-01 |   1.001274e+03 |   7.268414e-02 |   1.451520e+06 |   5.000000e-04 |
|   1.189278e+01 |   8.390025e-01 |   1.002774e+03 |   6.800265e-02 |   1.451520e+06 |   5.000000e-04 |
|   1.189328e+01 |   8.471114e-01 |   1.002910e+03 |   6.219615e-02 |   1.451520e+06 |   5.060536e-04 |
|   1.189378e+01 |   1.495091e+00 |   2.950757e+03 |   8.821367e-05 |   1.451520e+06 |   5.186652e-04 |
|   1.189430e+01 |   8.508068e-01 |   2.133388e+42 |   1.090975e-01 |   1.451520e+06 |   5.000000e-04 |
|   1.189480e+01 |   1.003630e+00 |            inf |   1.348090e-02 |   1.451520e+06 |   5.180453e-04 |
+----------------+----------------+----------------+----------------+----------------+----------------+


Time Step 21988, time = 11.8953, dt = 0.000518045

Here I output the global maximum value of the phase fieldmax_c, and I can see that it not only exceeds 1 (violating ConstantBounds), but also decreases (violating VariableOldValueBounds)

It should be noted that I have two input files that are linked by MultiApps, both of which are solved by Explicit TimeIntegrator, and I don't know if either of these features would bypass this limitation.

Is there something wrong with my settings? Any help is greatly appreciated!

Finally here is my input file

parent inputfile:

#use miug mius mm
[MultiApps]
  [fracture]
    type = TransientMultiApp
    input_files = fracture.i
    execute_on = TIMESTEP_END
    clone_parent_mesh = true
  []
[]

[Transfers]
  [from_c]
    type = MultiAppCopyTransfer
    from_multi_app = fracture
    variable = c
    source_variable = c
  []
  [to_shearband_driving_force]
    type = MultiAppCopyTransfer
    to_multi_app = fracture
    variable = shearband_driving_force
    source_variable = shearband_driving_force
  []
[]

[Mesh]
  allow_renumbering = false
  [TWC]
    type = AnnularMeshGenerator
    rmax = 5
    rmin = 3.5
    nt = 960
    nr = 378
  []
[]

[Problem]
  #restart_file_base = displacement_t_out_cp/LATEST  # You may also use a specific number here
[]

[Variables]
  [disp_x]
  []
  [disp_y]
  []
[]

[AuxVariables]
  [c]
  []
  [shearband_driving_force]
    order = CONSTANT
    family = MONOMIAL
    [AuxKernel]
      type = MaterialRealAux
      property = local_energy_dev
    []
  []
  [vonmises]
    order = CONSTANT
    family = MONOMIAL
    [AuxKernel]
      type = RankTwoScalarAux
      rank_two_tensor = stress
      variable = vonmises
      scalar_type = VonMisesStress
    []
  []
  [pressure]
    order = CONSTANT
    family = MONOMIAL
    [AuxKernel]
      type = RankTwoScalarAux
      rank_two_tensor = stress
      variable = pressure
      scalar_type = Hydrostatic
    []
  []
  [temp]
    order = CONSTANT
    family = MONOMIAL
    [AuxKernel]
      type = MaterialRealAux
      property = T
      variable = temp
    []
  []  
  [EPS]
    order = CONSTANT
    family = MONOMIAL
    [AuxKernel]
      type = MaterialRealAux
      property = eqv_plastic_strain
      variable = EPS
    []
  []
  [yield_factor]
    order = CONSTANT
    family = MONOMIAL
    [AuxKernel]
      type = MaterialRealAux
      property = yield_factor
      variable = yield_factor
    []
  []
  [accel_x]
  []
  [vel_x]
  []
  [accel_y]
  []
  [vel_y]
  []
[]


[AuxKernels]
  [accel_x]
      type = TestNewmarkTI
      variable = accel_x
      displacement = disp_x
      first = false
  []
  [vel_x]
      type = TestNewmarkTI
      variable = vel_x
      displacement = disp_x
  []
  [accel_y]
      type = TestNewmarkTI
      variable = accel_y
      displacement = disp_y
      first = false
  []
  [vel_y]
      type = TestNewmarkTI
      variable = vel_y
      displacement = disp_y
  []
[]

[Kernels]
  [DynamicSolidMechanics]
    displacements = 'disp_x disp_y'
    strain = FINITE
    incremental = true
    out_of_plane_direction = z
    planar_formulation = PLANE_STRAIN
    decomposition_method = HughesWinget
  []
  [inertia_x]
    type = InertialForce
    variable = disp_x
    density = density
    density_scaling = density_scaling
  []
  [inertia_y]
    type = InertialForce
    variable = disp_y
    density = density
    density_scaling = density_scaling
  []   
[]

[BCs]
  [Pressure]
      [outter_surface]
          boundary = 'rmax'
          function = current_pulse
          displacements = 'disp_x disp_y'
      []
  []
[]

[Functions]
  [current_pulse]
      #4Gpa=4e9n/m^2=4e9 kg/(m*s^2)=4e9 kg/(e3e12mm*mius^2)=4e3 miug/(mm*mius^2)
      type = PiecewiseLinear
      x = '0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1	1.1	1.2	1.3	1.4	1.5	1.6	1.7	1.8	1.9	2	2.1	2.2	2.3	2.4	2.5	2.6	2.7	2.8	2.9	3	3.1	3.2	3.3	3.4	3.5	3.6	3.7	3.8	3.9	4	4.1	4.2	4.3	4.4	4.5	4.6	4.7	4.8	4.9	5	5.1	5.2	5.3	5.4	5.5	5.6	5.7	5.8	5.9	6	6.1	6.2	6.3	6.4	6.5	6.6	6.7	6.8	6.9	7	7.1	7.2	7.3	7.4	7.5	7.6	7.7	7.8	7.9	8	8.1	8.2	8.3	8.4	8.5	8.6	8.7	8.8	8.9	9	9.1	9.2	9.3	9.4	9.5	9.6	9.7	9.8	9.9	10	10.1	10.2	10.3	10.4	10.5	10.6	10.7	10.8	10.9	11	11.1	11.2	11.3	11.4	11.5	11.6	11.7	11.8	11.9	12	12.1	12.2	12.3	12.4	12.5	12.6	12.7	12.8	12.9	13	13.1	13.2	13.3	13.4	13.5	13.6	13.7	13.8	13.9	14	14.1	14.2	14.3	14.4	14.5	14.6	14.7	14.8	14.9	15	15.1	15.2	15.3	15.4	15.5	15.6	15.7	15.8	15.9	16	16.1	16.2	16.3	16.4	16.5	16.6	16.7	16.8	16.9	17	17.1	17.2	17.3	17.4	17.5	17.6	17.7	17.8	17.9	18	18.1	18.2	18.3	18.4	18.5	18.6	18.7	18.8	18.9	19	19.1	19.2	19.3	19.4	19.5	19.6	19.7	19.8	19.9	20'
      y = '4.02E-02	4.03E-02	3.96E-02	3.93E-02	2.76E-02	-8.19E-03	-8.96E-02	-2.25E-01	-4.23E-01	-6.42E-01	-8.92E-01	-1.06E+00	-1.00E+00	-5.29E-01	4.22E-01	2.00E+00	4.27E+00	7.34E+00	1.12E+01	1.59E+01	2.15E+01	2.80E+01	3.54E+01	4.36E+01	5.25E+01	6.20E+01	7.22E+01	8.31E+01	9.45E+01	1.07E+02	1.19E+02	1.33E+02	1.47E+02	1.61E+02	1.76E+02	1.91E+02	2.07E+02	2.23E+02	2.40E+02	2.56E+02	2.73E+02	2.90E+02	3.07E+02	3.24E+02	3.41E+02	3.59E+02	3.76E+02	3.94E+02	4.11E+02	4.28E+02	4.46E+02	4.63E+02	4.81E+02	4.99E+02	5.17E+02	5.34E+02	5.52E+02	5.69E+02	5.86E+02	6.02E+02	6.19E+02	6.35E+02	6.51E+02	6.67E+02	6.83E+02	6.98E+02	7.13E+02	7.27E+02	7.41E+02	7.55E+02	7.69E+02	7.82E+02	7.95E+02	8.07E+02	8.18E+02	8.29E+02	8.39E+02	8.48E+02	8.56E+02	8.64E+02	8.70E+02	8.75E+02	8.80E+02	8.84E+02	8.87E+02	8.89E+02	8.90E+02	8.90E+02	8.89E+02	8.87E+02	8.84E+02	8.80E+02	8.75E+02	8.68E+02	8.60E+02	8.52E+02	8.43E+02	8.32E+02	8.21E+02	8.09E+02	7.96E+02	7.83E+02	7.69E+02	7.55E+02	7.40E+02	7.26E+02	7.12E+02	6.97E+02	6.82E+02	6.67E+02	6.50E+02	6.33E+02	6.15E+02	5.97E+02	5.79E+02	5.61E+02	5.42E+02	5.23E+02	5.03E+02	4.83E+02	4.64E+02	4.45E+02	4.26E+02	4.07E+02	3.88E+02	3.70E+02	3.52E+02	3.35E+02	3.17E+02	3.00E+02	2.82E+02	2.65E+02	2.49E+02	2.33E+02	2.17E+02	2.02E+02	1.87E+02	1.72E+02	1.59E+02	1.45E+02	1.32E+02	1.20E+02	1.08E+02	9.71E+01	8.63E+01	7.61E+01	6.63E+01	5.69E+01	4.83E+01	4.03E+01	3.31E+01	2.66E+01	2.09E+01	1.59E+01	1.15E+01	7.91E+00	4.95E+00	2.66E+00	9.90E-01	4.47E-02	-2.36E-01	2.15E-01	1.34E+00	3.22E+00	5.72E+00	8.89E+00	1.27E+01	1.72E+01	2.23E+01	2.81E+01	3.47E+01	4.18E+01	4.97E+01	5.81E+01	6.73E+01	7.71E+01	8.72E+01	9.79E+01	1.09E+02	1.20E+02	1.32E+02	1.44E+02	1.57E+02	1.70E+02	1.83E+02	1.96E+02	2.10E+02	2.24E+02	2.38E+02	2.52E+02	2.67E+02	2.82E+02	2.97E+02	3.12E+02	3.27E+02	3.42E+02	3.56E+02	3.71E+02	3.85E+02	3.99E+02	4.13E+02'
  []
[]

[Materials]
  [density]
    type = Density
    displacements = 'disp_x disp_y'
    density = 7.9e3 #miug/mm^3
  []
  [density_scaling]
    type = DensityScaling
    density = density
    desired_time_step = 0.0005
    output_properties = density_scaling
    outputs = 'exodus'
    factor = 0.8
  []
  [strain]
    type = ComputeFiniteStrain
    displacements = 'disp_x disp_y'
    implicit = false
    decomposition_method = HughesWinget
  []
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    shear_modulus = 77.5e3  #77.5Gpa
    poissons_ratio = 0.3
  []
  [damage_stress]
    type = JCEOSnlnotemp
    A = 110#Mpa
    B = 1500
    n = 0.36
    C = 0.025
    reference_strain_rate = 1
    m = 0.81
    T_m = 1800#kelvin
    T_r = 298
    T_i = 298
    SpecHeat = 500e-6 #1e-6mm^2/(mius^2*kelvin)=1*m*m/(s*s*kelvin)
    D_i = 7.9e3
    ShearModule = 77.5e3
    C0 = 4569e-3 #1e-3 mm/mius = 1m/s
    S1 = 1.49
    gamma0 = 2.17
    a = 0
    read_prop_user_object = prop_read
    c = c
    D_name = 'degradation'
    shearband_driving_force_name = local_energy_dev
  []
  [degradation]#a1
    type = DerivativeParsedMaterial
    property_name = degradation
    coupled_variables = 'c'
    expression = '((1.0-c)^2)/(((1.0-c)^2)+a1*c*(1-0.5*c))'
    constant_names       = 'a1'
    constant_expressions = '0.8'
    derivative_order = 2
  []
[]

[UserObjects]
  [prop_read]
    type = PropertyReadFile
    prop_file_name = 0.02722_0.5606.txt
    nprop = 1
    read_type = element
  []
[]

[Postprocessors]
  [numelemqps]
    type = NumElemQPs
  []
  [max_c]
    type = NodalExtremeValue
    variable = c
  []
  [max_von]
    type = ElementExtremeValue
    variable = vonmises
  []
  [min_degration]
    type = ElementExtremeMaterialProperty
    mat_prop = degradation
    value_type = min
  []
  [time_step]
    type = CriticalTimeStep
    factor = 0.8
    density = density
    density_scaling = density_scaling
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

[Executioner]
  type = Transient
  start_time = 0
  end_time = 20
  [TimeIntegrator]
    type = CentralDifference
    solve_type = lumped
  []
  [TimeStepper]
    type = PostprocessorDT
    postprocessor = time_step
  []
[]

[Outputs]
  [exodus]
    type = Exodus
    start_time = 0
    end_time = 20
    interval = 250
  []
  [exodus_detail]
    type = Exodus
    start_time = 11.7
    end_time = 20
    interval = 2
  []
  [checkpoint]
    type = Checkpoint
    interval = 100
    num_files = 10
  []
[]

sub inputfile

#use miug mius mm
[Mesh]
[]

[AuxVariables]
  [shearband_driving_force]
    order = CONSTANT
    family = MONOMIAL
  []
  [bounds_dummy]
    order = FIRST
    family = LAGRANGE
  []
[]

[Bounds]
  [irreversibility]
    type = VariableOldValueBounds
    variable = bounds_dummy
    bounded_variable = c
    bound_type = lower
  []
  [upper]
    type = ConstantBounds
    variable = bounds_dummy
    bounded_variable = c
    bound_type = upper
    bound_value = 1.0
  []
[]

[Variables]
  [c]
    order = FIRST
    family = LAGRANGE
  []
[]

[Kernels]
  [trans]
    type = TimeDerivative
    variable = c
  []
  [Bulk]
    type = AllenCahn
    variable = c
    f_name = F
    mob_name = L
  []
  [ACInterface]
    type = ACInterface
    variable = c
    kappa_name = kappa_op
    mob_name = L
  []
[]

[Materials]
  [pfbulkmat]
    type = GenericConstantMaterial
    prop_names = 'gc_prop l visco'
    prop_values = '6.666666667	0.003968254	1.26001008' # *2/12 before:/8

    #16.8kj/m^2 = 16.8 1e3 N*m/m^2 = 16.8 1e3 N/m = 16.8 10e3 kg/s^2 = 16.8 miug/mius^2
    #10 mium = 1e-2 mm
    #yita~(zhang) = yita* c_a * l/gc = 0.01 mius, c_a = 2
    #yita(origin) = 0.01*gc/(c_a*l) miug/(mm * mius)
    #visco(moose) = yita(origin)/gc = yita~(zhang)/(c_a*l) = 0.5 mius/mm
  []
  [define_mobility]
    type = ParsedMaterial
    material_property_names = 'gc_prop visco'
    property_name = L
    expression = '1.0/(gc_prop * visco)'
  []
  [define_kappa]
    type = ParsedMaterial
    material_property_names = 'gc_prop l'
    property_name = kappa_op
    expression = 'gc_prop * l'#kappa = 2gc * l/c_a), c_a = 2
  []
  [local_fracture_energy]#c_a = 2, a(d) = d^2, f_frac = (1/c_a) * (gc/l) * a(d)
    type = DerivativeParsedMaterial
    property_name = local_fracture_energy
    coupled_variables = 'c'
    material_property_names = 'gc_prop l'
    expression = '(c^2 * gc_prop)/(2 * l)'
    derivative_order = 2
  []
  [fracture_driving_energy]
    type = DerivativeParsedMaterial
    coupled_variables = 'c shearband_driving_force'
    material_property_names = 'local_fracture_energy(c) degradation(c)'
    expression = 'local_fracture_energy + degradation * shearband_driving_force'
    derivative_order = 2
    property_name = F
  []
  [degradation]#a1
    type = DerivativeParsedMaterial
    property_name = degradation
    coupled_variables = 'c'
    expression = '((1.0-c)^2)/(((1.0-c)^2)+a1*c*(1-0.5*c))'
    constant_names       = 'a1'
    constant_expressions = '0.8'
    derivative_order = 2
  []
[]

[Executioner]
   type = Transient
   petsc_options_iname = '-pc_type -pc_factor_shift  -snes_type'
   petsc_options_value = 'lu NONZERO vinewtonrsls'

  [TimeIntegrator]
    type = ActuallyExplicitEuler
    solve_type = lumped
  []
[]

Answered by GiudGiud

May 16, 2024

Hello

Can you use an implicit time stepper?
The explicit time stepper skips a lot of steps, because of the reduced solve needs. I would not be surprised if bounds are not included in explicit time stepping

Guillaume

View full answer

GiudGiud · 2024-05-16T14:30:40Z

GiudGiud
May 16, 2024
Collaborator

Hello

Can you use an implicit time stepper?
The explicit time stepper skips a lot of steps, because of the reduced solve needs. I would not be surprised if bounds are not included in explicit time stepping

Guillaume

1 reply

Yasunaoka May 18, 2024
Author

Hello
Thank you for your suggestion. After using the implicit time stepper, they worked properly.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Question about ConstantBounds and VariableOldValueBounds not working #27645

{{title}}

{{editor}}'s edit

{{editor}}'s edit

Replies: 1 comment 1 reply

{{title}}

{{title}}

Select a reply

Question about ConstantBounds and VariableOldValueBounds not working #27645

Yasunaoka May 16, 2024

Check these boxes if you have followed the posting rules.

Question

Replies: 1 comment · 1 reply

GiudGiud May 16, 2024 Collaborator

Yasunaoka May 18, 2024 Author

Yasunaoka
May 16, 2024

Replies: 1 comment 1 reply

GiudGiud
May 16, 2024
Collaborator

Yasunaoka May 18, 2024
Author