用于解决 2D 热棒问题的 Fortran 代码

问：

在我使用之前，我只是使用转置来转置T_l。代码可以成功 mpirun。但是，在我使用您提供的代码后，代码无法成功 mpirun。以下是我使用后的代码以及终端给我的错误消息。你能通读它们并给我建议来修复错误吗？您可以丢弃子例程，这绝对是正确的，因为唯一更改的部分是 ghost swap 部分，如果我简单地使用转置，不会发生错误。MPI_Type_vectorMPI_Type_vectorMPI_Type_vector

我的代码：

program asn
    use mpi

    ! Parameter declaration
    implicit none
    integer :: ierr, processID, totalProcesses, nx, ny, i, j, r_start, r_end, numRows, iterationCount, &
    itemax, m_start, m_end, comm_count, lowProcessID, upProcessID, upSendTag, lowSendTag, dims(1), &
    request(4),rowDimensions(2),assembleDimensions(2),globalDimensions(2),startPoint(2), rowtype
    real(8) :: a, b, dx, dy, maxResidual, localMaxResidual, globalMaxResidual, pi, L
    real(8), allocatable :: x(:), y(:), U(:,:), V(:,:), U_l(:,:), V_l(:,:), T(:,:), T_l(:,:),&
    T_assemble(:,:), T_l_transpose(:,:), residual(:,:)
    logical :: periodicBoundary(1), reorderProcesses
    integer :: resized_send_type,resized_receive_type,cart_comm, L_row, typeSize, send_type, recv_type
    integer, allocatable :: rowsPerProcess(:), displacement(:)
    INTEGER(KIND=MPI_ADDRESS_KIND) :: startMemory, strideMemory

    ! MPI Initialization:
    ! Start the MPI environment, get the rank of the current process and the total number of processes.
    call MPI_INIT(ierr)
    call MPI_COMM_RANK(MPI_COMM_WORLD, processID, ierr)
    call MPI_COMM_SIZE(MPI_COMM_WORLD, totalProcesses, ierr)

    ! Solution Array Initialization:
    ! Define the dimensions of the solution array and allocate memory for it. Initialize all elements to zero.
    nx = 100
    ny = 100
    allocate(T(ny, nx))
    T = 0.0d0

    ! Physical and Numerical Parameters Initialization:
    ! Set the domain length, coefficients 'a' and 'b', and calculate grid spacings 'dx' and 'dy'.
    L = 1d0
    a = 80
    b = 0.5d0
    pi = 3.14159260
    dx = L/(nx-1)
    dy = L/(ny-1)

    ! Position Vector Initialization:
    ! Allocate memory for position vectors 'x' and 'y' and compute their values based on grid spacings.
    allocate(x(nx))
    allocate(y(ny))
    do i =1,nx
        x(i) = (i-1)*dx
        y(i) = (i-1)*dy
    end do

    ! Boundary Conditions Setup:
    ! Set the boundary values for the solution array 'T' based on a given function.
    do i = 1, nx
        T(1,i) = -cos(2 * pi * x(i) / L)**2
        T(ny,i) = -T(1,i)
    end do

    ! Velocity Field Initialization:
    ! Allocate memory for velocity fields 'U' and 'V' and compute their values based on given functions.
    allocate(U(ny,nx))
    allocate(V(ny,nx))
    do i = 1, ny
        do j = 1, nx
            U(i,j) = sin(pi * x(i) / L) * cos(pi * y(j) / L) &
                     +sin(2 * pi * x(i) / L) * cos(2 * pi * y(j) / L)
            V(i,j) = -cos(pi * x(i) / L) * sin(pi * y(j) / L) &    
                     -cos(2 * pi * x(i) / L) * sin(2 * pi * y(j) / L)
        end do
    end do 

    ! MPI Topology Creation:
    ! Create a one-dimensional Cartesian topology for processes and determine neighboring processes.
    dims(1) = totalProcesses
    periodicBoundary(1) = .false.
    reorderProcesses = .true.
    call MPI_CART_CREATE(MPI_COMM_WORLD, 1, dims, periodicBoundary, reorderProcesses, cart_comm, ierr)
    call MPI_CART_SHIFT(cart_comm, 0, 1, lowProcessID, upProcessID, ierr)

    ! Domain Decomposition:
    ! Decompose the computational domain among the available processes based on their rank.
    call DomainDecomposition(ny, totalProcesses, processID, r_start, r_end, m_start, m_end)
    
    ! Local Solution Array Initialization:
    ! Allocate memory for local solution arrays and copy the relevant sections from the global solution array.
    numRows = r_end-r_start+1
    allocate(T_l(numRows,nx))
    allocate(T_l_transpose(nx,numRows))
    T_l(:,:) = T(r_start:r_end,:)

    ! Local Velocity Field Initialization:
    ! Allocate memory for local velocity fields and copy the relevant sections from the global velocity fields.
    allocate(U_l(numRows,nx))
    allocate(V_l(numRows,nx))
    U_l(:,:) = U(r_start:r_end,:)
    V_l(:,:) = V(r_start:r_end,:)

    ! Ghost Cell Communication Initialization:
    ! Define tags for MPI communication to identify the data being sent/received.
    upSendTag = 1
    lowSendTag = 2

    ! Convergence Criteria Initialization:
    ! Allocate memory for the residual array and set the maximum allowed residual and maximum iteration count.
    allocate(residual(numRows,nx))
    iterationCount = 0
    maxResidual = 1d-5
    itemax = 20000

!                           2. calculate and update the local domain
!============================================================================================================
    ! Loop until not converged after large iteration
    do while (iterationCount < itemax)
!------------------------------------------------------------------------------------------------------------
    ! For red nodes
    ! The red nodes are updated based on the parity (odd or even) of 'r_start + 1'.
    ! Depending on the parity, different starting points and steps are used for the update.
    
    ! Check the parity of 'r_start + 1'
    if (mod(r_start + 1, 2) == 1) then
        ! If 'r_start + 1' is odd:
        
        ! Update the matrix for odd rows starting from the second row and third column.
        ! The updates are done in steps of 2 for both rows and columns.
        call updateMatrix(2, 3, 2, 2, T_l, U_l, V_l, a, b, dx, dy, numRows, nx)
        ! Update the matrix for even rows starting from the third row and second column.
        call updateMatrix(3, 2, 2, 2, T_l, U_l, V_l, a, b, dx, dy, numRows, nx)
    else
        ! If 'r_start + 1' is even:
        
        ! Update the matrix for even rows starting from the second row and second column.
        call updateMatrix(2, 2, 2, 2, T_l, U_l, V_l, a, b, dx, dy, numRows, nx)
        ! Update the matrix for odd rows starting from the third row and third column.
        call updateMatrix(3, 3, 2, 2, T_l, U_l, V_l, a, b, dx, dy, numRows, nx)
    end if

!------------------------------------------------------------------------------------------------------------
! Ghost row swap:
! This section of the code deals with the communication of ghost rows between adjacent processes.
! Ghost rows are boundary rows of a local domain that are used to exchange information with neighboring domains.
! Using MPI_Type_vector, we can directly send and receive rows without needing to transpose the matrix.

! Define a new MPI data type for a row of the matrix:
call MPI_Type_vector(1, nx, nx, MPI_DOUBLE_PRECISION, rowtype, ierr)
call MPI_Type_commit(rowtype, ierr)

! Initialize communication count:
comm_count = 0

! Receive ghost row from the lower local domain (if not the first process):
if (processID > 0) then
    call MPI_IRECV(T_l(1,:), 1, rowtype, lowProcessID, upSendTag, &
                   cart_comm, request(comm_count+1), ierr)
    comm_count = comm_count + 1
end if

! Send ghost row to the upper local domain (if not the last process):
if (processID < totalProcesses-1) then
    call MPI_ISEND(T_l(numRows-1,:), 1, rowtype, upProcessID, upSendTag, &
                   cart_comm, request(comm_count+1), ierr)
    comm_count = comm_count + 1
end if

! Receive ghost row from the upper local domain (if not the last process):
if (processID < totalProcesses-1) then
    call MPI_IRECV(T_l(numRows,:), 1, rowtype, upProcessID, lowSendTag, &
                   cart_comm, request(comm_count+1), ierr)
    comm_count = comm_count + 1
end if

! Send ghost row to the lower local domain (if not the first process):
if (processID > 0) then
    call MPI_ISEND(T_l(2,:), 1, rowtype, lowProcessID, lowSendTag, &
                   cart_comm, request(comm_count+1), ierr)
    comm_count = comm_count + 1
end if

! Wait for all non-blocking communication operations to complete:
if (comm_count > 0) then
    call MPI_WAITALL(comm_count, request, MPI_STATUSES_IGNORE, ierr)
end if

! Free the custom MPI data type:
call MPI_Type_free(rowtype, ierr)


!-------------------------------------------------------------------------------------------------------------
    ! For black nodes
    ! The black nodes are updated similarly to the red nodes but with different starting points.
    ! The parity of 'r_start + 1' determines the starting points for the update.

    ! Check the parity of 'r_start + 1'
    if (mod(r_start + 1, 2) == 1) then
        ! If 'r_start + 1' is odd:
        
        ! Update the matrix for even rows starting from the second row and second column.
        call updateMatrix(2, 2, 2, 2, T_l, U_l, V_l, a, b, dx, dy, numRows, nx)
        ! Update the matrix for odd rows starting from the third row and third column.
        call updateMatrix(3, 3, 2, 2, T_l, U_l, V_l, a, b, dx, dy, numRows, nx)
    else
        ! If 'r_start + 1' is even:
        
        ! Update the matrix for odd rows starting from the second row and third column.
        call updateMatrix(2, 3, 2, 2, T_l, U_l, V_l, a, b, dx, dy, numRows, nx)
        ! Update the matrix for even rows starting from the third row and second column.
        call updateMatrix(3, 2, 2, 2, T_l, U_l, V_l, a, b, dx, dy, numRows, nx)
    end if

!-------------------------------------------------------------------------------------------------------------
! Ghost row swap:
! This section of the code deals with the communication of ghost rows between adjacent processes.
! Ghost rows are boundary rows of a local domain that are used to exchange information with neighboring domains.
! Using MPI_Type_vector, we can directly send and receive rows without needing to transpose the matrix.

! Define a new MPI data type for a row of the matrix:
call MPI_Type_vector(1, nx, nx, MPI_DOUBLE_PRECISION, rowtype, ierr)
call MPI_Type_commit(rowtype, ierr)

! Initialize communication count:
comm_count = 0

! Receive ghost row from the lower local domain (if not the first process):
if (processID > 0) then
    call MPI_IRECV(T_l(1,:), 1, rowtype, lowProcessID, upSendTag, &
                   cart_comm, request(comm_count+1), ierr)
    comm_count = comm_count + 1
end if

! Send ghost row to the upper local domain (if not the last process):
if (processID < totalProcesses-1) then
    call MPI_ISEND(T_l(numRows-1,:), 1, rowtype, upProcessID, upSendTag, &
                   cart_comm, request(comm_count+1), ierr)
    comm_count = comm_count + 1
end if

! Receive ghost row from the upper local domain (if not the last process):
if (processID < totalProcesses-1) then
    call MPI_IRECV(T_l(numRows,:), 1, rowtype, upProcessID, lowSendTag, &
                   cart_comm, request(comm_count+1), ierr)
    comm_count = comm_count + 1
end if

! Send ghost row to the lower local domain (if not the first process):
if (processID > 0) then
    call MPI_ISEND(T_l(2,:), 1, rowtype, lowProcessID, lowSendTag, &
                   cart_comm, request(comm_count+1), ierr)
    comm_count = comm_count + 1
end if

! Wait for all non-blocking communication operations to complete:
if (comm_count > 0) then
    call MPI_WAITALL(comm_count, request, MPI_STATUSES_IGNORE, ierr)
end if

! Free the custom MPI data type:
call MPI_Type_free(rowtype, ierr)


!------------------------------------------------------------------------------------------------------------
! update and check stopping criteria

        ! update iteration notation
        iterationCount = iterationCount+1
        ! calculate the maximum residual on the global domain 
        ! this calculate residual on each local node
        do i = 2, numRows-1
            do j = 2, ny-1
            ! residual = B-Ax, where B is RHS of equation, Ax is LHS of equation
                residual(i,j) = 0 - a * U_l(i,j) * ((T_l(i+1,j) - T_l(i-1,j)) / (2 * dx)) &
                              - a * V_l(i,j) * ((T_l(i,j+1) - T_l(i,j-1)) / (2* dy)) &
                              + ((T_l(i+1,j) - 2 * T_l(i,j) + T_l(i-1,j)) / dx**2) &
                              + ((T_l(i,j+1) - 2 * T_l(i,j) + T_l(i,j-1)) / dy**2) &
                              + b * T_l(i,j)
            end do
        end do
        
        localMaxResidual = maxval(abs(residual)) ! calculate the maximum residual over the local domain
        call MPI_ALLREDUCE(localMaxResidual, globalMaxResidual, 1, MPI_DOUBLE_PRECISION,&
         MPI_MAX, cart_comm, ierr) !calculate the maximum residual over all domain

        ! check if the maximum residual is smaller than setted
        if (globalMaxResidual < maxResidual) then
            exit
        end if
!------------------------------------------------------------------------------------------------------------
    end do

我的错误：

 Program received signal SIGABRT: Process abort signal.

Backtrace for this error:

Could not print backtrace: executable file is not an executable

Primary job  terminated normally, but 1 process returned
a non-zero exit code. Per user-direction, the job has been aborted.
--------------------------------------------------------------------------
--------------------------------------------------------------------------
mpirun noticed that process rank 1 with PID 0 on node Daniels-Air exited on signal 6 (Abort trap: 6).

在我使用之前，我只是使用转置来转置T_l。代码可以成功 mpirun。但是，在我使用您提供的代码后，代码无法成功 mpirun。以下是我使用后的代码以及终端给我的错误消息。你能通读它们并给我建议来修复错误吗？您可以丢弃子例程，这绝对是正确的，因为唯一更改的部分是 ghost swap 部分，如果我简单地使用转置，不会发生错误。MPI_Type_vectorMPI_Type_vectorMPI_Type_vector