Incorrect boundary enforcement of DirectDirichletBC

[TheGreatCid]

Bug Description

In the process of merging the DirectDirichletBC type, the calculation of the residual in the base class was changed from

  Real resid = (computeQpValue() - _sys.solutionOld()(dofnum)) / (avg_dt * _dt) -
               (*_sys.solutionUDotOld())(dofnum) / avg_dt;

to

  Real resid = (computeQpValue() - _u_old[_qp]) / (avg_dt * _dt) - _u_dot_old[_qp] / avg_dt;

where _u_old= _var.slnOld() and _u_dot_old = _var.uDotOld()

These should be equivalent, and indeed the tests did not catch any changes in these cases. However I found that there a sort of compounding error in the solution at the boundaries as number of time steps increases. This error eventually leads to the incorrect enforcement of the boundary condition. This error is not present in the original method of indexing of the solution vector.

I found that in the DirectDirichletBCBase class, _u_old[_qp] != _sys.solutionOld()(dofnum). This behavior is puzzling and fixing this should fix the class.

Steps to Reproduce

Running this file will show the issue. The BC at the bottom of the plate will be incorrectly enforced.

[Problem]
    extra_tag_matrices = 'mass'
[]
[Mesh]
    [plate]
        type = GeneratedMeshGenerator
        dim = 2
        nx = 5
        ny = 5
        xmax = 10
        xmin = -10
        ymin = -10
        ymax = 10
        boundary_id_offset = 10
        boundary_name_prefix = plate
        subdomain_ids = 2
    []
    [plate_transform]
        type = TransformGenerator
        input = plate
        transform = TRANSLATE
        vector_value = '100 100 0'
    []
[]

[GlobalParams]
    displacements = 'disp_x disp_y'
[]

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

[AuxVariables]

    [vel_x]
    []
    [vel_y]
    []
    [vel_z]
    []
    [accel_x]
    []
    [accel_y]
    []
    [accel_x]
    []
[]

[Kernels]
    [Mass_x]
        type = MassMatrix
        variable = disp_x
        matrix_tags = 'mass'
        density = density
    []
    [Mass_y]
        type = MassMatrix
        variable = disp_y
        matrix_tags = 'mass'
        density = density
    []
    [DynamicTensorMechanics]
        displacements = 'disp_x disp_y'
        stiffness_damping_coefficient = 1e-8
        generate_output = 'stress_zz strain_zz'
    []
[]

[AuxKernels]
    [accel_x]
        type = TestNewmarkTI
        variable = accel_x
        displacement = disp_x
        first = false
        execute_on = 'LINEAR TIMESTEP_BEGIN TIMESTEP_END'
    []
    [vel_x]
        type = TestNewmarkTI
        variable = vel_x
        displacement = disp_x
        execute_on = 'LINEAR TIMESTEP_BEGIN TIMESTEP_END'
    []
    [accel_y]
        type = TestNewmarkTI
        variable = accel_y
        displacement = disp_y
        first = false
        execute_on = 'LINEAR TIMESTEP_BEGIN TIMESTEP_END'
    []
    [vel_y]
        type = TestNewmarkTI
        variable = vel_y
        displacement = disp_y
        execute_on = 'LINEAR TIMESTEP_BEGIN TIMESTEP_END'
    []
[]
[Materials]
    [density_two]
        type = GenericConstantMaterial
        prop_names = density
        prop_values = 2710
        outputs = 'exodus'
        output_properties = 'density'
    []
    [wave_speed]
        type = WaveSpeed
        outputs = 'exodus'
        output_properties = 'wave_speed'
    []

    [elasticity_rubber]
        type = ComputeIsotropicElasticityTensor
        youngs_modulus = 1e6
        poissons_ratio = 0.33
    []
    [strain_block]
        type = ComputeSmallStrain # ComputeIncrementalSmallStrain
        displacements = 'disp_x disp_y'
        implicit = false
    []
    [stress_block]
        type = ComputeLinearElasticStress
    []
[]

[BCs]
    [fix_x]
        type = DirectDirichletBC
        boundary = 'plate_bottom' # 'plate_back plate_left plate_right'
        value = 0
        variable = 'disp_x'
    []
    [fix_y]
        type = DirectDirichletBC
        boundary = 'plate_bottom' # 'plate_back plate_left plate_right'
        value = 0
        variable = 'disp_y'
    []
    [Pull]
        type = DirectFunctionDirichletBC
        boundary = plate_top
        function = 't'
        variable = disp_y
    []
[]

[Postprocessors]
    [critical_time_step]
        type = CriticalTimeStep
    []
[]

[Executioner]
    type = Transient
    dt = 2e-3

    end_time = 0.2
    [TimeIntegrator]
        type = DirectCentralDifference
        mass_matrix_tag = 'mass'
    []
[]

[Outputs]
    [exodus]
        type = Exodus
        file_base = '2D'
        additional_execute_on = FAILED
    []
    print_linear_residuals = false
    interval = 1
[]

Impact

This gives incorrect values and should be fixed ASAP.

[Optional] Diagnostics

No response

[TheGreatCid]

@lindsayad @hugary1995 @dschwen

[hugary1995]

One of the reinit methods must have not been called (at the correct time).

[TheGreatCid]

@dschwen and I found that using _u_old(_var.nodalValueOld()) instead of _u_old(_var.slnOld()) returns the correct value. So there is certainly something odd going on here.

[TheGreatCid]

@lindsayad Should I open a PR using _u_old(_var.nodalValueOld())? and just ignore the issue with _u_old(_var.slnOld()) for now?

[lindsayad]

You should probably do that since nodalValueOld() is probably the best named API to use for a nodal BC. In which case you can either ref this issue and leave it open or mark this as closed and we can open a new one for the framework issue