chop_restart.i

# ------------------------ 

# slope of the bottom boundary (in degrees)
bed_slope = 15 

# change coordinate system to add a slope
gravity_x = ${fparse
  	      cos((90 - bed_slope) / 180 * pi) * 9.81
              } 
gravity_z = ${fparse
	      cos(bed_slope / 180 * pi) * 9.81
              } 

# account for slope in BC on displacement  
# bed_factor_x = ${fparse
#              	 cos(bed_slope / 180 * pi)
#                 }
# bed_factor_z = ${fparse
# 	         sin(bed_slope / 180 * pi)
#                 }

# ------------------------

[Mesh]
  [channel]      
  type = FileMeshGenerator
    file = chop_init_out.e
  []

  [chop_area]
    type = SubdomainBoundingBoxGenerator
    input = channel
    block_id = 2
    bottom_left = '29500  -100 0'
    top_right = '30000 2100 1000'
  []

  [chop_event]
    type = BlockDeletionGenerator
    input = chop_area
    block = 2
    new_boundary = downstream
  []

[]

[Variables]
  [disp_x]
    order = FIRST
    family = LAGRANGE
  []
  [disp_y]
    order = FIRST
    family = LAGRANGE
  []
  [disp_z]
    order = FIRST
    family = LAGRANGE
  []
[]

[AuxVariables]
  [vel_x]
  []
  [accel_x]
  []
  [vel_y]
  []
  [accel_y]
  []
  [vel_z]
  []
  [accel_z]
  []
  [stress_xx]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_xy]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_xz]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_yx]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_yz]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_zx]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_zy]
    order = CONSTANT
    family = MONOMIAL
  []
  [stress_zz]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Functions]
  [weight]
    type = ParsedFunction
    value = '-8829*(1000-z)'    # initial stress that should result from the weight force
  []
  [upstream_dirichlet]
    type = ParsedFunction
    value = '0'
  []
  [frontal_pressure]
    type = ParsedFunction
    value = '8829*(1000-z)'   
  []
  [solution_disp_x]
    type = SolutionFunction
    from_variable = disp_x
    solution = init_solution
  []
  [solution_disp_y]
    type = SolutionFunction
    from_variable = disp_y
    solution = init_solution
  []
  [solution_disp_z]
    type = SolutionFunction
    from_variable = disp_z
    solution = init_solution
  []
[]

[Kernels]
  [gravity_x]
    type = Gravity
    variable = disp_x
    value= ${gravity_x}
  []
  [gravity_y]
    type = Gravity
    variable = disp_y
    value = 0.
  []
  [gravity_z]
    type = Gravity
    variable = disp_z
    value = -${gravity_z}
  []
  [DynamicTensorMechanics]
    stiffness_damping_coefficient = 0.02
    mass_damping_coefficient = 0.02
    displacements = 'disp_x disp_y disp_z'
  []
  [inertia_x]
    type = InertialForce
    variable = disp_x
    velocity = vel_x
    acceleration = accel_x
    beta = 0.25
    gamma = 0.5
  []
  [inertia_y]
    type = InertialForce
    variable = disp_y
    velocity = vel_y
    acceleration = accel_y
    beta = 0.25
    gamma = 0.5
  []
  [inertia_z]
    type = InertialForce
    variable = disp_z
    velocity = vel_z
    acceleration = accel_z
    beta = 0.25
    gamma = 0.5
  []
[]

[AuxKernels]
  [accel_x]
    type = SolutionAux
    variable = accel_x
    from_variable = accel_x
    solution = init_solution
  []
  [vel_x]
    type = SolutionAux
    variable = vel_x
    from_variable = vel_x
    solution = init_solution
  []
  [accel_y]
    type = SolutionAux
    variable = accel_y
    from_variable = accel_y
    solution = init_solution
  []
  [vel_y]
    type = SolutionAux
    variable = vel_y
    from_variable = vel_y
    solution = init_solution
  []
  [accel_z]
    type = SolutionAux
    variable = accel_z
    from_variable = accel_z
    solution = init_solution
  []
  [vel_z]
    type = SolutionAux
    variable = vel_z
    from_variable = vel_z
    solution = init_solution
  []
  [stress_xx]
    type = SolutionAux
    variable = stress_xx
    from_variable = stress_xx
    solution = init_solution
  []
  [stress_xy]
    type = SolutionAux
    variable = stress_xy
    from_variable = stress_xy
    solution = init_solution
  []
  [stress_xz]
    type = SolutionAux
    variable = stress_xz
    from_variable = stress_xz
    solution = init_solution
  []
  [stress_yx]
    type = SolutionAux
    variable = stress_yx
    from_variable = stress_yx
    solution = init_solution
  []
  [stress_yy]
    type = SolutionAux
    variable = stress_yy
    from_variable = stress_yy
    solution = init_solution
  []
  [stress_yz]
    type = SolutionAux
    variable = stress_yz
    from_variable = stress_yz
    solution = init_solution
  []
  [stress_zx]
    type = SolutionAux
    variable = stress_zx
    from_variable = stress_zx
    solution = init_solution
  []
  [stress_zy]
    type = SolutionAux
    variable = stress_zy
    from_variable = stress_zy
    solution = init_solution
  []
  [stress_zz]
    type = SolutionAux
    variable = stress_zz
    from_variable = stress_zz
    solution = init_solution
  []
[]

[Materials]
  [ice_elasticity]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 8.7e9 #Pa
    poissons_ratio = 0.31
  []
  [strain]
    type = ComputeIncrementalSmallStrain
    displacements = 'disp_x disp_y disp_z'
  []
  [density]
    type = GenericConstantMaterial
    prop_names = density
    prop_values = 900 #kg/m3
  []
  [stress]
    type = ComputeFiniteStrainElasticStress
  []
  [von_mises]
    type = RankTwoInvariant
    invariant = 'VonMisesStress'
    property_name = von_mises
    rank_two_tensor = stress
    outputs = exodus
  []
[]

[UserObjects]
  [init_solution]
    type = SolutionUserObject
    mesh = chop_init_out.e
    timestep = LATEST
  []
[]

[ICs]
  [init_disp_x]
    type = FunctionIC
    variable = disp_x
    function = solution_disp_x
  []
  [init_disp_y]
    type = FunctionIC
    variable = disp_y
    function = solution_disp_y
  []
  [init_disp_z]
    type = FunctionIC
    variable = disp_z
    function = solution_disp_z
  []
[]

[BCs]
  [upstream_dirichlet]
    type = DirichletBC                                               
    boundary = 'upstream'
    variable = disp_x
    value    = 0.0
  []  
  [Pressure]
    [downstream_pressure]  
    boundary = downstream
    function = frontal_pressure
    displacements = 'disp_x disp_y disp_z'
    []
  []
  [anchor_bottom_x]
    type = DirichletBC
    variable = disp_x
    boundary = 'bottom'
    value = 0.0
  []  
  [anchor_botom_y]
    type = DirichletBC
    variable = disp_y
    boundary = 'bottom'
    value = 0.0
  []
  [anchor_bottom_z]
    type = DirichletBC
    variable = disp_z
    boundary = 'bottom'
    value = 0.0
  []

  [anchor_sides_x]
    type = DirichletBC
    variable = disp_x
    boundary = 'left right'
    value = 0.0
  []
  [anchor_sides_y]
    type = DirichletBC
    variable = disp_y
    boundary = 'left right'
    value = 0.
  []
  [anchor_sides_z]
    type = DirichletBC
    variable = disp_z
    boundary = 'left right'
    value = 0.0
  []

[]

[Preconditioning]
  [andy]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
  petsc_options = '-ksp_snes_ew'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu       superlu_dist'
  solve_type = 'NEWTON'
  nl_rel_tol = 1e-7
  nl_abs_tol = 1e-12
  dt = 0.02
  end_time = 50.
  timestep_tolerance = 1e-6
  automatic_scaling = true
  [TimeIntegrator]
    type = NewmarkBeta
    beta = 0.25
    gamma = 0.5
  []
[]

[Outputs]
  exodus = true
  perf_graph = true
[]