Problem with serial decomposition of extremely big meshes, use PT-Scotch (the parallel version of SCOTCH) in SPECFEM3D instead of (or as an option in addition to) the serial version of SCOTCH

[pbasini]

There is a sort of bug that appears when dealing with big meshes, in my case something like 65 millions of elements: when running xdecompose_mesh I get:

total number of nodes:
   nnodes =     66076729
 total number of spectral elements:
   nspec =     65270709
 materials:
   num_mat =            1
   defined =            1 undefined =            0
 absorbing boundaries:
   nspec2D_xmin =       170235
   nspec2D_xmax =        60379
   nspec2D_ymin =        27992
   nspec2D_ymax =        19563
   nspec2D_bottom =       481552
   nspec2D_top =       829214
 Par_file_faults not found: assume no faults
 node valence:
   min =            1  max =           12
   nsize =           12 sup_neighbour =           70
forrtl: severe (408): fort: (3): Subscript #1 of the array ADJNCY has value -401485046 which is less than the lower bound of 0

Image              PC                Routine            Line        Source            
xdecompose_mesh    00000000004F3B5A  Unknown               Unknown  Unknown
xdecompose_mesh    00000000004F2656  Unknown               Unknown  Unknown
xdecompose_mesh    00000000004AD2A0  Unknown               Unknown  Unknown
xdecompose_mesh    0000000000475FCE  Unknown               Unknown  Unknown
xdecompose_mesh    00000000004764F1  Unknown               Unknown  Unknown
xdecompose_mesh    000000000042403F  part_decompose_me         130  part_decompose_mesh.f90
xdecompose_mesh    00000000004371B5  decompose_mesh_mp         808  decompose_mesh.F90
xdecompose_mesh    000000000045A302  MAIN__                     93  program_decompose_mesh.f90
xdecompose_mesh    0000000000403A3C  Unknown               Unknown  Unknown
libc.so.6          00002B499DBE076D  Unknown               Unknown  Unknown
xdecompose_mesh    0000000000403939  Unknown               Unknown  Unknown

the problem is in the term nodes_elmnts(l+j_nsize)_sup_neighbour, in line 130 of part_decompose_mesh.f90, because it exceeds the limit of 2147483646 for a 4byte integer.

With Elliot we found a temporary solution to this problem by modifying decompose_mesh.F90 and part_decompose_mesh.f90 as follows:

decompose_mesh.F90

integer(kind=8), dimension(:), allocatable  :: xadj
integer(kind=8) :: nb_edges
integer(kind=8) :: max_neighbour   ! real maximum number of neighbours per element

part_decompose_mesh.f90

integer(kind=8), dimension(0:nspec)  :: xadj
integer(kind=8), intent(out) :: max_neighbour
integer(kind=8)  :: i, j, k, l, m, nb_edges
integer(kind=8) :: tmp

As you noticed we simply redefined the kind for those variables we need, then the variable tmp will be defined three times as follows

tmp=int(nodes_elmnts(k+j*nsize),kind=8)*sup_neighbour+m
tmp=int(nodes_elmnts(k+j*nsize),kind=8)*sup_neighbour &
                          + xadj(nodes_elmnts(k+j*nsize))
tmp = int(nodes_elmnts(l+j*nsize), kind=8) * sup_neighbour &
                       + xadj(nodes_elmnts(l+j*nsize))

and will be used as the index for adjncy.

Everything seems to be fine, the code runs through this part without a problem but it's not yet time for a happy ending. In fact I get the following error:

total number of nodes:
   nnodes =     66076729
 total number of spectral elements:
   nspec =     65270709
 materials:
   num_mat =            1
   defined =            1 undefined =            0
 absorbing boundaries:
   nspec2D_xmin =       170235
   nspec2D_xmax =        60379
   nspec2D_ymin =        27992
   nspec2D_ymax =        19563
   nspec2D_bottom =       481552
   nspec2D_top =       829214
 Par_file_faults not found: assume no faults
 node valence:
   min =            1  max =           12
   nsize =           12 sup_neighbour =           70
 mesh2dual:
   max_neighbour =                     38
 1
ERROR: graphCheck: arc data do not match
ERROR : MAIN : Invalid check

This happens right after call scotchfgraphcheck (scotchgraph (1), ier) in decompose_mesh.F90. I think that this can be caused by the fact that in the previous call:

call scotchfgraphbuild (scotchgraph (1), 0, nspec, &
                          xadj (1), xadj (1), &
                          elmnts_load (1), xadj (1), &
                          nb_edges, adjncy (1), &
                          adjncy (1), ier)

among the arguments passed to scotchfgraphbuild there are also xadj and nb_edges, that are now defined with kind=8. Temporary solution, but not sure of the results because I still have to generate the database:

decompose_mesh.F90
allocate(tmp_xadj(1:nspec+1),stat=ier)
if( ier /= 0 ) stop 'error allocating array xadj'

then, instead of calling scotchfgraphbuild, I first make sure that nb_edges is smaller than the 4bit limit and then I translate the variables defined as kind=8 into kind=4, and I pass those new variables to scotchfgraphbuild:

if( nb_edges < 2147483646 )then
      tmp_edges=int(nb_edges,kind=4)
      tmp_xadj=int(xadj,kind=4)

      call scotchfgraphbuild (scotchgraph (1), 0, nspec, &
                          tmp_xadj (1), tmp_xadj (1), &
                          elmnts_load (1), tmp_xadj (1), &
                          tmp_edges, adjncy (1), &
                          adjncy (1), ier)

As I said, I am not sure this could fix the whole thing: on the cluster (scinet consortium) even if I ask a node with 64 Gb of RAM I cannot pass through the scotchfgraphpart (scotchgraph (1), nparts, scotchstrat(1),part(1),ier) call because of lack in memory. On my workstation xdecompose_mesh runs and finishes without errors but it creates the proc000000_Database of 5.5 Gb and all the others of 300 bytes.

[komatits]

To summarize, the question is: for very (very!) large meshes, is there a way of calling the SCOTCH domain decomposer with Fortran integer(kind=8) input arrays rather than kind=4?

(in C that would correspond to "long long int" input arrays instead of "long int")

On 25/04/2014 16:24, Piero Basini wrote:

There is a sort of bug that appears when dealing with big meshes, in my case something like 65 millions of elements: when running xdecompose_mesh I get:

|total number of nodes: nnodes = 66076729 total number of spectral elements: nspec = 65270709 materials: num_mat = 1 defined = 1 undefined = 0 absorbing boundaries: nspec2D_xmin = 170235 nspec2D_xmax = 60379 nspec2D_ymin = 27992 nspec2D_ymax = 19563 nspec2D_bottom = 481552 nspec2D_top = 829214 Par_file_faults not found: assume no faults node valence: min = 1 max = 12 nsize = 12 sup_neighbour = 70 forrtl: severe (408): fort: (3): Subscript #1 of the array ADJNCY has value -401485046 which is less than the lower bound of 0

Image PC Routine Line Source xdecompose_mesh 00000000004F3B5A Unknown Unknown Unknown xdecompose_mesh 00000000004F2656 Unknown Unknown Unknown xdecompose_mesh 00000000004AD2A0 Unknown Unknown Unknown xdecompose_mesh 0000000000475FCE Unknown Unknown Unknown xdecompose_mesh 00000000004764F1 Unknown Unknown Unknown xdecompose_mesh 000000000042403F part_decompose_me 130 part_decompose_mesh.f90 xdecompose_mesh 00000000004371B5 decompose_mesh_mp 808 decompose_mesh.F90 xdecompose_mesh 000000000045A302 MAIN__ 93 program_decompose_mesh.f90 xdecompose_mesh 0000000000403A3C Unknown Unknown Unknown libc.so.6 00002B499DBE076D Unknown Unknown Unknown xdecompose_mesh 0000000000403939 Unknown Unknown Unknown |

the problem is in the term nodes_elmnts(l+j_nsize)_sup_neighbour, in line 130 of part_decompose_mesh.f90, because it exceeds the limit of 2147483646 for a 4byte integer.

With Elliot we found a temporary solution to this problem by modifying decompose_mesh.F90 and part_decompose_mesh.f90 as follows:

|decompose_mesh.F90

integer(kind=8), dimension(:), allocatable :: xadj integer(kind=8) :: nb_edges integer(kind=8) :: max_neighbour ! real maximum number of neighbours per element |

|part_decompose_mesh.f90

integer(kind=8), dimension(0:nspec) :: xadj integer(kind=8), intent(out) :: max_neighbour integer(kind=8) :: i, j, k, l, m, nb_edges integer(kind=8) :: tmp |

As you noticed we simply redefined the kind for those variables we need, then the variable tmp will be defined three times as follows

|tmp=int(nodes_elmnts(k+j_nsize),kind=8)_sup_neighbour+m tmp=int(nodes_elmnts(k+j_nsize),kind=8)_sup_neighbour &

• xadj(nodes_elmnts(k+j_nsize)) tmp = int(nodes_elmnts(l+j_nsize), kind=8) * sup_neighbour &
  • xadj(nodes_elmnts(l+j*nsize)) |

and will be used as the index for adjncy.

Everything seems to be fine, the code runs through this part without a problem but it's not yet time for a happy ending. In fact I get the following error:

total number of nodes: nnodes = 66076729 total number of spectral elements: nspec = 65270709 materials: num_mat = 1 defined = 1 undefined = 0 absorbing boundaries: nspec2D_xmin = 170235 nspec2D_xmax = 60379 nspec2D_ymin = 27992 nspec2D_ymax = 19563 nspec2D_bottom = 481552 nspec2D_top = 829214 Par_file_faults not found: assume no faults node valence: min = 1 max = 12 nsize = 12 sup_neighbour = 70 mesh2dual: max_neighbour = 38 1 ERROR: graphCheck: arc data do not match ERROR : MAIN : Invalid check

This happens right after |call scotchfgraphcheck (scotchgraph (1), ier)| in decompose_mesh.F90. I think that this can be caused by the fact that in the previous call:

call scotchfgraphbuild (scotchgraph (1), 0, nspec, & xadj (1), xadj (1), & elmnts_load (1), xadj (1), & nb_edges, adjncy (1), & adjncy (1), ier)

among the arguments passed to scotchfgraphbuild there are also xadj and nb_edges, that are now defined with kind=8. Temporary solution, but not sure of the results because I still have to generate the database:

decompose_mesh.F90 allocate(tmp_xadj(1:nspec+1),stat=ier) if( ier /= 0 ) stop 'error allocating array xadj'

then, instead of calling scotchfgraphbuild, I first make sure that nb_edges is smaller than the 4bit limit and then I translate the variables defined as kind=8 into kind=4, and I pass those new variables to scotchfgraphbuild:

|if( nb_edges < 2147483646 )then tmp_edges=int(nb_edges,kind=4) tmp_xadj=int(xadj,kind=4)

   call scotchfgraphbuild (scotchgraph (1), 0, nspec, &
                       tmp_xadj (1), tmp_xadj (1), &
                       elmnts_load (1), tmp_xadj (1), &
                       tmp_edges, adjncy (1), &
                       adjncy (1), ier)

As I said, I am not sure this could fix the whole thing: on the cluster (scinet consortium) even if I ask a node with 64 Gb of RAM I cannot pass through the |scotchfgraphpart (scotchgraph (1), nparts, scotchstrat(1),part(1),ier)| call because of lack in memory. On my workstation xdecompose_mesh runs and finishes without errors but it creates the proc000000_Database of 5.5 Gb and all the others of 300 bytes.

— Reply to this email directly or view it on GitHub https://github.com/geodynamics/specfem3d/issues/129.

Dimitri Komatitsch CNRS Research Director (DR CNRS), Laboratory of Mechanics and Acoustics, UPR 7051, Marseille, France http://komatitsch.free.fr

[komatits]

Answer from François Pellegrini, the developer of Scotch in Bordeaux (France):

Yes, since about 2010. Just read the INSTALL.txt file at the root of the Scotch tree and compile with -DINTSIZE64 and -DIDXSIZE64.

[komatits]

Piero @pbasini and Elliott @QuLogic, I guess the best way to proceed is then probably to make your SPECFEM patch to change the arrays from integer(4) to integer(8) permanent and then also permanently add -DINTSIZE64 and -DIDXSIZE64 to the Scotch compiler options that we use in our "configure" script, right?

[pbasini]

Dimitri @komatits : I've tried the solution proposed by Pellegrini, and I have also changed the definition of all the integers entering the scotch calls from integer(4) to integer(8), this assuming that by compiling scotch with -DINTSIZE64 all the integers are integer(8). Nevertheless I am still getting the error ERROR: SCOTCH_graphBuild: invalid base parameter. I will further investigate this, to make sure that I did not miss any important variable.

[komatits]

Hi Piero @pbasini,

Thanks! Did you try with both -DINTSIZE64 and -DIDXSIZE64?

On 28/04/2014 18:27, Piero Basini wrote:

Dimitri @komatits https://github.com/komatits : I've tried the solution proposed by Pellegrini, and I have also changed the definition of all the integers entering the scotch calls from integer(4) to integer(8), this assuming that by compiling scotch with -DINTSIZE64 all the integers are integer(8). Nevertheless I am still getting the error ERROR: SCOTCH_graphBuild: invalid base parameter. I will further investigate this, to make sure that I did not miss any important variable.

— Reply to this email directly or view it on GitHub https://github.com/geodynamics/specfem3d/issues/129#issuecomment-41580031.

Dimitri Komatitsch CNRS Research Director (DR CNRS), Laboratory of Mechanics and Acoustics, UPR 7051, Marseille, France http://komatitsch.free.fr

[pbasini]

Hi Dimitri @komatits

Yes I did, and actually I found the error: in the graphbuild call we were using as base value 0, without defining it as a variable, thus scotch was expecting a 64 integer. What I did was to modify decompose_mesh.F90: I defined the variable basevalue as an integer(8) and I assigned to it the value 0. Now the code runs through the scotchfgraphbuild call. Anyway it is taking a while to finish, and thus I cannot tell you yet that everything went fine.

OK @komatits ! It worked: the decomposition ended and now I have 10 files of different size, while before following Pellegrini's comment I had 1 file of 5.5 Gb and 9 file of 300 bytes. Now the only problem left is that on the cluster I am using for the simulations there could be a memory issue: I will try this new version of the decomposition and let you know. If the problem will remain we will probably have to think about moving toward ptscotch.

I have a question related to the "cluster-issue": could it be possible to decompose the mesh on a different machine and then move the proc*Database files on the cluster? Assuming that xdecompose_mesh has been compiled with the same compiler and options used to compile xgenerate_databases and xspecfem3d on the cluster.

[komatits]

Hi Piero @pbasini and Elliott @QuLogic ,

Thanks! I am glad the problem is solved. You can then commit it to "devel" once it is fully tested I guess.

Please check that the "base value = 0" problem that you mention was specific to the 64-bit integer case (one would think that in case of a bug in the base value the 32-bit case would crash as well).

[pbasini]

@komatits unfortunately the problem is not solved or at least not completely: on the cluster (on a node with 64 Gb of RAM) the code now pass through the 5 calls to scotch but when it is time to allocate the tab_interfaces (line 378 of part_decompose_mesh.f90) it crashes because of a lack of memory. I think the solution now is to migrate toward ptscotch or to find a way to free memore before this call is made.

[komatits]

Use PT-Scotch (the parallel version of SCOTCH) in SPECFEM3D instead of (or as an option in addition to) the serial version of SCOTCH; because decomposing very large 3D meshes on a serial machine (on a single processor core) can be a severe limitation and for some meshes that will be a major bottleneck.

Support ParMetis as well because there is a new, improved version that has been released by Karypis et al recently. To do that, Jeroen sent a Fortran bug fix (a patch to use to call METIS from Fortran) to Matthieu by email on Jan 16, 2013.

What would need to be done is switch to PT-Scotch version 6.0 and/or to ParMETIS (supporting both would be easy, since the subroutine calls are very similar; in the current serial code we support both Scotch and METIS).

To do that, one would need to rewrite "decompose_mesh.F90" in directory src/decompose_mesh. That is the only file that would need to change (in addition to part_decompose_mesh.f90, which is called by decompose_mesh.F90). Unfortunately that file is not very well written nor very well commented. Thus unfortunately that set of routines is not so easy to understand and modify. A bunch of things are hardwired.

Thus, this is not so easy to do (Daniel and I checked about a year ago). It would maybe be easier to rewrite 'decompose_mesh_scotch' from scratch as a new, MPI program instead of trying to add MPI support to the current and not so well written serial version. That program is relatively small, and what it does is relatively straightforward.

[komatits]

For very large meshes the serial domain decomposer runs out of memory because of the size of the connectivity table on a single machine, and that is independent of the operating system.

The clean solution is to use a parallel domain decomposer (PT-Scotch or ParMetis) instead of their serial version (Scotch or Metis); that is planned but will take a few months.

In the meantime something worth testing (and easy to test) is to use Metis instead of Scotch, you can easily replace the call to Scotch with a call to Metis. The call to Metis (off by default) is at line 925 of file src/decompose_mesh/decompose_mesh.F90 . That calls Metis through Scotch (Scotch has an internal way of doing nothing and calling Metis instead), if that fails as well then try calling Metis directly (i.e. download Metis and replace that line with a real call to Metis). If that fails then it is the adjacency table (arrays "xadj" and "adjncy") that is too big and then there is nothing we can do with the serial version.

[komatits]

Vadim Monteiller @vmont is almost done implementing the new parallel mesh decomposer, he will commit it in the next few days.

[komatits]

Done! (parallel mesh decomposer created by Vadim Monteiller @vmont )