Unexpected behavior when using GridGenerator::subdivided_hyper_rectangle in parallel

[Jimmy Ho]

Hi All,

I am using the Step 40 tutorial to build a parallel program using MPI. The code runs but generates different results when using one and multiple processors. After stripping it down to the bare minimum, it appears that when the mesh is built using GridGenerator::subdivided_hyper_rectangle without any subsequent refinement, the compress() function does not work properly, leading to a different global stiffness when using multiple processors.

A minimum example to reproduce this is attached. When the mesh is built using GridGenerator::hyper_cube or GridGenerator::subdivided_hyper_rectangle with subsequent refinement, the program works as expected. When the same mesh is generated using GridGenerator::subdivided_hyper_rectangle without any subsequent refinement, some entries in the global stiffness matrix (the nodes on the right hand side of the mesh, in this case, using 2 processors) do not get updated. The example output of stiffness matrices when using one and two processors are also attached for reference.

So my question is, is this the expected behavior of the function? If so, why is that the case?

Thanks a lot for any answers! Your inputs are highly appreciated!

Best,

Jimmy

[Wolfgang Bangerth]

Jimmy,

> A minimum example to reproduce this is attached. When the mesh is built using
> GridGenerator::hyper_cube or GridGenerator::subdivided_hyper_rectangle with
> subsequent refinement, the program works as expected. When the same mesh is
> generated using GridGenerator::subdivided_hyper_rectangle without any
> subsequent refinement, some entries in the global stiffness matrix (the nodes
> on the right hand side of the mesh, in this case, using 2 processors) do not
> get updated. The example output of stiffness matrices when using one and two
> processors are also attached for reference.
>
> So my question is, is this the expected behavior of the function? If so, why
> is that the case?

There is still a lot of stuff in the program that could be remove, including
all of the comments, to make it substantially smaller and easier to understand.

I don't know whether there is a bug, but here is a suggestion: The finite
element solution u_h(x) that results from the linear system should be the same
independent of the partitioning. But the order of degrees of freedom may be
different, and consequently the matrix may not be the same -- it should only
be the same up to some column and row permutation. Have you verified that the
*solution function* (not the solution vector) that results is the same
independent of the number of processors?

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Jimmy Ho]

Hi Dr. Bangerth,

Thanks a lot for your guidance! I compared the solution in the vtu files using the minimal example above, they are nearly identical. Looking back into the code, I am outputting the system matrix from processor 0, which probably only printed the part that it locally owns, hence there's a difference in the matrices when running in serial and parallel. I guess I rushed to conclusion too quickly.

Best,

Jimmy

[Jimmy Ho]

Hi Dr. Bangerth,

As a follow-up question, upon calling compress(), will the local copy of the system matrix on a specific processor get updated to contain information from all other processors? In other words, if I print out the system matrix from a particular processor after calling compress(), is that the same global system matrix that the linear solver is solving?

Thanks a lot for clarifying!

Best,

Jimmy

On Wednesday, July 29, 2020 at 9:55:38 AM UTC-5, Wolfgang Bangerth wrote:

[Wolfgang Bangerth]

On 7/30/20 10:11 AM, Jimmy Ho wrote:
>
> As a follow-up question, upon calling compress(), will the local copy of the
> system matrix on a specific processor get updated to contain information from
> all other processors? In other words, if I print out the system matrix from a
> particular processor after calling compress(), is that the same global system
> matrix that the linear solver is solving?

Each processor only stores a subset of the rows of the matrix. During
assembly, each processor computes a number of matrix entries, but it can not
compute all entries for the rows of the matrix it owns -- some need
contributions from other processes. The call to compress() makes sure the
contributions from these other processes are sent to the one that owns these
rows of the matrix.

In any case, after compress(), the rows each processor owns are correct -- but
each processor doesn't know anything about the rows of the matrix it doesn't own.

[Jimmy Ho]

Hi Dr. Bangerth,

Thanks a lot for the clarification! They are really helpful!

Best,

Jimmy