SegFault using LinearNodalConstraint with Distributed Mesh

[TheBEllis]

Hi MOOSE team,
I've been trying to use LinearNodalConstraint to replicate an RBE3-esque constraint on a problem that will need a distributed mesh, and I am running into a SegFault. I have put together a small reproducible example. I have uploaded the mesh here https://drive.google.com/file/d/1ZlKYpy96iZ3vGojwCULTtBbGSVbFBUxG/view?usp=sharing.

The independent (primary) nodes in this situation are all the nodes on the cylindrical extrude face, and the dependent (secondary) node is located at 0,0,0 at the centre of the mesh. All the weightings are equal. I have pre-split the mesh into 2 parts. Using a debug executable I immediately get this error, which points to libmesh_assert_valid_neighbors.

Using an opt executable and running a debugger for each MPI process, the backtrace looks like this,

Also I am running b9eefd8, so aaedfb9 is included.

Any help is much appreciated:) The problem converges fine without --use-split, but as soon as the distributed mesh is used the segfault appears.

Edit: After a bit of digging it seems the problem is the elemPtr() call at line 459 of SystemBase.C, in augmentSendList(). Should this be a queryElemPtr()?

penalty = 1e13

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

[Mesh]
  type = FileMesh
  file = '../mesh/rbetest_no_elems.e'
[] 
  
[Variables]
  [temperature]
    initial_condition = 0 # Start at room temperature
  []
[]
  
[Kernels]
  # Heat conduction kernels
  [heat_conduction]
    type = ADHeatConduction
    variable = temperature
  []
[]

[Modules/TensorMechanics/Master]
  [all]
    use_automatic_differentiation = true
    volumetric_locking_correction = true
    strain = SMALL
    add_variables = true
    incremental = false
    generate_output = 'strain_xx strain_yy strain_zz vonmises_stress'
    eigenstrain_names = eigenstrain
  []
[]
  
[Materials]
  [CCZ_elasticity]
    type = ADComputeIsotropicElasticityTensor
    youngs_modulus = 100000
    poissons_ratio = 0.33
  []
  [CCZ_heat]
    type = ADHeatConductionMaterial
    temp = temperature
    thermal_conductivity = 400.0
    specific_heat = 330.0
  []

  [CCZ_thermal_expansion]
    type = ADComputeThermalExpansionEigenstrain
    stress_free_temperature = 10
    thermal_expansion_coeff = 5e-04
    temperature = temperature
    eigenstrain_name = eigenstrain
    block = 1
  []

  [rbe_thermal_expansion]
    type = ADComputeThermalExpansionEigenstrain
    stress_free_temperature = 0
    thermal_expansion_coeff = 0
    temperature = temperature
    eigenstrain_name = eigenstrain
    block = 2
  []

  [stress]
    type = ADComputeLinearElasticStress
  []
[]

[Constraints]
    [rbe3_x]
        type = LinearNodalConstraint
        primary = '7 1 149 8 190 191 192 193 13 194 195 196 197 198 199 200 14 15 16 17 18 19 20 21 22 23 2 64 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922'
        secondary_node_set = 'secondary'
        variable = disp_x
        penalty = ${penalty}
        weights = '0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005'
    []

    [rbe3_y]
        type = LinearNodalConstraint
        primary = '7 1 149 8 190 191 192 193 13 194 195 196 197 198 199 200 14 15 16 17 18 19 20 21 22 23 2 64 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922'
        secondary_node_set = 'secondary'
        variable = disp_y
        penalty = ${penalty}
        weights = '0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005'
    []

    [rbe3_z]
        type = LinearNodalConstraint
        primary = '7 1 149 8 190 191 192 193 13 194 195 196 197 198 199 200 14 15 16 17 18 19 20 21 22 23 2 64 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922'
        secondary_node_set = 'secondary'
        variable = disp_z
        penalty = ${penalty}
        weights = '0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005 0.005'
    []
[]

[BCs]
  [temp_high]
    type = ADDirichletBC
    variable = temperature
    boundary = 'hot'
    value = 100
  []
  [temp_low]
    type = ADDirichletBC
    variable = temperature
    boundary = 'cold'
    value = 0
  []
[]

[Postprocessors]
  [average_x]
    type = AverageNodalVariableValue
    variable = disp_x
    boundary = 'primary'
  []
  [average_y]
    type = AverageNodalVariableValue
    variable = disp_y
    boundary = 'primary'
  []
  [average_z]
    type = AverageNodalVariableValue
    variable = disp_z
    boundary = 'primary'
  []
[]
  
#[Preconditioning]
#  [./SMP]
#    #Creates the entire Jacobian, for the Newton solve
#    type = SMP
#    full = true
#  [../]
#[]
  
[Executioner]
  type = Steady
  automatic_scaling = true
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre    boomeramg'
  nl_rel_tol = 1e-20
  nl_abs_tol = 1e-5
  l_tol = 1e-8
  l_abs_tol = 1e-10
  line_search = none
[]

[Outputs]
  exodus = true
  perf_graph = true
[]

[GiudGiud]

hello

Can you try this instead:

• generate the split with --split-mesh as usual
• load the split mesh using a FileMeshGenerator (FileMesh -- like you have now -- can work too probably)
• do not use --use-split

use-split follows a different code path that does not do a good job at keeping track of neighbors iirc

Guillaume

[TheBEllis]

Hi Guillaume,
This worked, thank you lots! I should have paid more attention to this.

Note for anyone else who is trying to do this, the only other thing I needed to get it working was to set skip_partitioning to true in FileMeshGenerator, as my mesh was pre-split and it is set to False by default (and use the --distributed-mesh command line argument when running).
Again, many thanks,
Bill