Explicit Dynamics: Initial Velocity Condition

[ddm42]

Hello! I am trying to run a problem using the "CentralDifference" TimeIntegrator derived from the ActuallyExplicitTimeIntergrator class, and there doesn't seem to be a way to set an initial velocity condition on the body. Is there a way to do this? I tried to implement a quick fix in CentralDifference::computeTimeDerivativeHelper() where I could make temporary u_old and u_older variables on the first time step, but I am unable to add or subtract a number from those classes. Is there a better way to implement this?

[GiudGiud]

Hello

Is velocity a variable?
Or is it computed from the derivative of the displacement?

[ddm42]

It is computed from the derivative of displacement. I used to extract it as an AuxVariable through the implicit dynamics system. And we were also able to set the initial condition through the auxvariable.

[cbolisetti]

@ddm42 can you try using the testnewmarkTI auxkernel as described in this example? Another thing to make sure is that your displacement history (if that's what your input is), is consistent with the velocity and acceleration initial conditions. If no initial conditions are set, moose will assume they are all zero, which is not always the case when you input displacements.
Also, I think @dschwen is working on improving all of this, so he might have more updated information.

[ddm42]

Yes I have also been using the testnewmarkTI auxkernel to extract the velocity and accelertion. This is my input file. The example also did not use ICs. Also, I am in Daniel Schwen's group and he was unable to find a solution.

[Mesh]
  coord_type = RZ

  [rock]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 10
    ny = 80
    # nx = 20
    # ny = 10
    xmin = 0
    xmax = ${units 0.1 in -> m}
    ymin = ${units -10 in -> m}
    ymax = 0
  []
  [rock_block]
    type = SubdomainIDGenerator
    input = rock
    subdomain_id = 1
  []
[]

[GlobalParams]
  displacements = 'disp_x disp_y'
  large_kinematics = true
[]

[Problem]
  kernel_coverage_check = false
[]

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

[AuxVariables]
  [vel_x]
    block = 1
  []
  [vel_y]
    block = 1
  []
  [accel_x]
    block = 1
  []
  [accel_y]
    block = 1
  []
  [T]
    initial_condition = ${units x degC}
  []
  [wall]
  []
[]

[Kernels]
  [sdx]
    type = TotalLagrangianStressDivergenceAxisymmetricCylindrical
    # type = TotalLagrangianStressDivergence
    variable = disp_x
    component = 0
    block = 1
  []
  [sdy]
    type = TotalLagrangianStressDivergenceAxisymmetricCylindrical
    # type = TotalLagrangianStressDivergence
    variable = disp_y
    component = 1
    block = 1
  []

  [ifx]
    type = InertialForce
    variable = disp_x
    # velocity = vel_x
    # acceleration = accel_x
    # alpha = 0 # not implemented in TotalLagrangianStressDivergence
    block = 1
  []
  [ify]
    type = InertialForce
    variable = disp_y
    # velocity = vel_y
    # acceleration = accel_y
    # alpha = 0 # not implemented in TotalLagrangianStressDivergence
    block = 1
  []
[]

[AuxKernels]
  [vel_x]
    type = TestNewmarkTI
    variable = vel_x
    displacement = disp_x
    first = true
    execute_on = 'LINEAR TIMESTEP_BEGIN TIMESTEP_END'
  []
  [vel_y]
    type = TestNewmarkTI
    variable = vel_y
    displacement = disp_y
    first = true
    execute_on = 'LINEAR TIMESTEP_BEGIN TIMESTEP_END'
  []

  [accel_x]
    type = TestNewmarkTI
    variable = accel_x
    displacement = disp_x
    first = false
    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'
  []

  [wall]
    type = ParsedAux
    variable = wall
    coupled_variables = 'vel_y'
    expression = 'if(y>0, -y*1e15-vel_y*1e8, 0)'
    use_xyzt = true
    use_displaced_mesh = true
  []
[]

[Materials]
  [rock_strain]
    type = ComputeLagrangianStrainAxisymmetricCylindrical
    # type = ComputeLagrangianStrain
    block = 1
  []
  [rock_stress]
    type = ComputeLagrangianLinearElasticStress
    block = 1
  []
  [elasticity_rock]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = ${units 20 GPa -> Pa}
    poissons_ratio = 0.3
    block = 1
  []
  [rock_density]
    type = Density
    density = ${units x kg/m^3}
  []
[]

[ICs]
  [vel_y]
    type = ConstantIC
    variable = vel_y
    value = ${units x m/s}
    block = 1
  []
[]

[BCs]
  [wall]
    type = CoupledVarNeumannBC
    variable = disp_y
    boundary = top
    v = wall
    use_displaced_mesh = true
  []
  [symmetry]
    type = DirichletBC
    variable = disp_x
    value = 0
    boundary = left
  []
[]

[Postprocessors]
  [vel_y]
    type = ElementAverageValue
    variable = vel_y
    block = 1
  []
  [gap_y]
    type = SideAverageValue
    variable = disp_y
    boundary = top
  []
[]

[Executioner]
  type = Transient
  [TimeIntegrator]
    type = CentralDifference
    solve_type = lumped
  []
  solve_type = LINEAR
  dt = ${units 1e-8 s}
  num_steps = 2000
[]

[Outputs]
  exodus = true
  csv = true
  print_linear_residuals = false
  perf_graph = true
  time_step_interval = 1
[]

[ddm42]

I am currently trying to make a fix in the initialConditionTempl.C file, where I am trying to set the variable being affected by the IC to an older solution state if specified. Then all the IC objects could work on old solution states, similarly to the current state.

Specifically I have added this line of code in the constructor of "InitialConditionTempl" with the appropriate corresponding changes in the header file and validParams().

_var( getParam("_set_old_soln_state") ? _sys.solutionState(getParam("_old_soln_state_number")) : _sys.getActualFieldVariable(parameters.get<THREAD_ID>("_tid"),parameters.get("variable")) ),

The full constructor is below.

template
InitialConditionTempl::InitialConditionTempl(const InputParameters & parameters)
: InitialConditionBase(parameters),
_fe_problem(*getCheckedPointerParam<FEProblemBase *>("_fe_problem_base")),
_tid(getParam<THREAD_ID>("_tid")),
_t(_fe_problem.time()),
_var((getParam("_set_old_soln_state")) ? _sys.solutionState(getParam("_old_soln_state_number")) : _sys.getActualFieldVariable(parameters.get<THREAD_ID>("_tid"),parameters.get("variable"))),
_fe_var(dynamic_cast<MooseVariableFE *>(&_var)),
_assembly(
_fe_problem.assembly(_tid, _var.kind() == Moose::VAR_SOLVER ? _var.sys().number() : 0)),
_coord_sys(_assembly.coordSystem()),
_current_elem(_var.currentElem()),
_current_elem_volume(_assembly.elemVolume()),
_current_node(nullptr),
_qp(0),
_fe_type(_var.feType())
{
}

The problem is that the old states and current field variables are stored as different data types. How do I access the old state in a way that matches the field variable data type? I need this if it is to be compatible with the current IC system.
@lindsayad

[lindsayad]

You don't need to change the initialization of _var. What you would need to do in order to set "initial conditions" for the old value of _var is to change what's at the end of InitialConditionTempl<T>::compute(). Instead of a call to _var.setDofValue, you would need a call to something like _var.setDofValueOld. That API doesn't currently exist though. If you add it, then I would make it private and make InitialConditionTempl a friend template. We do not want other code to have the ability to modify old values

[ddm42]

Great thank you for the help