From 4941bdde72bfa61d3d3ad6884ad65f14f226fd07 Mon Sep 17 00:00:00 2001
From: Chloe Mimeau <Chloe.Mimeau@imag.fr>
Date: Thu, 25 Jul 2013 10:35:19 +0000
Subject: [PATCH] update advec_line
---
.../{advecX_line.F90 => advecX_line.f90} | 27 +-
.../particles/advec_line/advecY_line.f90 | 27 +-
.../particles/advec_line/advecZ_line.f90 | 27 +-
.../advec_line/advec_common_line.f90 | 115 +-
.../advec_line/advec_remesh_formula_line.f90 | 9 +-
.../advec_line/advec_remesh_line.F90 | 1312 +++++++++++++++++
6 files changed, 1417 insertions(+), 100 deletions(-)
rename HySoP/src/scalesInterface/particles/advec_line/{advecX_line.F90 => advecX_line.f90} (90%)
create mode 100644 HySoP/src/scalesInterface/particles/advec_line/advec_remesh_line.F90
diff --git a/HySoP/src/scalesInterface/particles/advec_line/advecX_line.F90 b/HySoP/src/scalesInterface/particles/advec_line/advecX_line.f90
similarity index 90%
rename from HySoP/src/scalesInterface/particles/advec_line/advecX_line.F90
rename to HySoP/src/scalesInterface/particles/advec_line/advecX_line.f90
index 80aabe3a2..e311ef87d 100644
--- a/HySoP/src/scalesInterface/particles/advec_line/advecX_line.F90
+++ b/HySoP/src/scalesInterface/particles/advec_line/advecX_line.f90
@@ -1,3 +1,4 @@
+!USEFORTEST advec
!> @addtogroup part
!------------------------------------------------------------------------------
@@ -77,22 +78,22 @@ subroutine advecX_calc_line(dt,Vx,scal3D)
! Input/Output
real(WP), intent(in) :: dt
- real(WP), dimension(N_proc(1), N_proc(2), N_proc(3)), intent(in) :: Vx
- real(WP), dimension(N_proc(1), N_proc(2), N_proc(3)), intent(inout) :: scal3D
+ real(WP), dimension(mesh_sc%N_proc(1), mesh_sc%N_proc(2), mesh_sc%N_proc(3)), intent(in) :: Vx
+ real(WP), dimension(mesh_sc%N_proc(1), mesh_sc%N_proc(2), mesh_sc%N_proc(3)), intent(inout) :: scal3D
! Other local variables
integer :: j,k ! indice of the currend mesh point
integer, dimension(2) :: ind_group ! indice of the currend group of line (=(i,k) by default)
integer :: direction=1 ! current direction = along Y
- real(WP), dimension(N_proc(1)) :: p_pos_adim ! adimensionned particles position
- real(WP), dimension(N_proc(1)) :: p_V ! particles velocity
+ real(WP), dimension(mesh_sc%N_proc(1)) :: p_pos_adim ! adimensionned particles position
+ real(WP), dimension(mesh_sc%N_proc(1)) :: p_V ! particles velocity
logical, dimension(bl_nb(1)+1) :: bl_type ! is the particle block a center block or a left one ?
logical, dimension(bl_nb(1)) :: bl_tag ! indice of tagged particles
ind_group = 0
- do k = 1, N_proc(3)
+ do k = 1, mesh_sc%N_proc(3)
ind_group(2) = ind_group(2) + 1
ind_group(1) = 0
- do j = 1, N_proc(2)
+ do j = 1, mesh_sc%N_proc(2)
ind_group(1) = ind_group(1) + 1
! ===== Init particles =====
@@ -102,7 +103,7 @@ subroutine advecX_calc_line(dt,Vx,scal3D)
! -- Compute velocity (with a RK2 scheme) --
call AC_particle_velocity_line(dt, direction, ind_group, p_pos_adim, p_V)
! -- Advec particles --
- p_pos_adim = p_pos_adim + dt*p_V/d_sc(direction)
+ p_pos_adim = p_pos_adim + dt*p_V/mesh_sc%dx(direction)
! ===== Remeshing =====
! -- Pre-Remeshing: Determine blocks type and tag particles --
@@ -147,7 +148,7 @@ subroutine Xremesh_O2_line(ind_group, p_pos_adim, bl_type, bl_tag,j,k,scal)
logical, dimension(:), intent(in) :: bl_type
logical, dimension(:), intent(in) :: bl_tag
real(WP), dimension(:), intent(in) :: p_pos_adim
- real(WP), dimension(N_proc(1), N_proc(2), N_proc(3)), intent(inout) :: scal
+ real(WP), dimension(mesh_sc%N_proc(1), mesh_sc%N_proc(2), mesh_sc%N_proc(3)), intent(inout) :: scal
! Other local variables
! Variable used to remesh particles in a buffer
real(WP),dimension(:),allocatable :: send_buffer ! buffer use to remesh the scalar before to send it to the right subdomain
@@ -165,12 +166,12 @@ subroutine Xremesh_O2_line(ind_group, p_pos_adim, bl_type, bl_tag,j,k,scal)
! First particle is a left one
send_j_min = floor(p_pos_adim(1))-1
end if
- if (bl_type(N_proc(direction)/bl_size +1)) then
+ if (bl_type(mesh_sc%N_proc(direction)/bl_size +1)) then
! Last particle is a centered one
- send_j_max = nint(p_pos_adim(N_proc(direction)))+1
+ send_j_max = nint(p_pos_adim(mesh_sc%N_proc(direction)))+1
else
! Last particle is a left one
- send_j_max = floor(p_pos_adim(N_proc(direction)))+1
+ send_j_max = floor(p_pos_adim(mesh_sc%N_proc(direction)))+1
end if
! -- Determine the communication needed : who will communicate whit who ? (ie compute sender and recer) --
@@ -210,12 +211,12 @@ subroutine advecX_init_line(Vx, j, k, p_pos_adim, p_V)
! Input/Output
integer, intent(in) :: j,k
- real(WP), dimension(N_proc(direction)), intent(out) :: p_pos_adim, p_V
+ real(WP), dimension(mesh_sc%N_proc(direction)), intent(out) :: p_pos_adim, p_V
real(WP), dimension(:,:,:), intent(in) :: Vx
! Other local variables
integer :: ind ! indice
- do ind = 1, N_proc(direction)
+ do ind = 1, mesh_sc%N_proc(direction)
p_pos_adim(ind) = ind
p_V(ind) = Vx(ind,j,k)
end do
diff --git a/HySoP/src/scalesInterface/particles/advec_line/advecY_line.f90 b/HySoP/src/scalesInterface/particles/advec_line/advecY_line.f90
index f050f2474..2e7d1c4f0 100644
--- a/HySoP/src/scalesInterface/particles/advec_line/advecY_line.f90
+++ b/HySoP/src/scalesInterface/particles/advec_line/advecY_line.f90
@@ -1,3 +1,4 @@
+!USEFORTEST advec
!> @addtogroup part
!! @{
@@ -74,21 +75,21 @@ subroutine advecY_calc_line(dt,Vy,scal3D)
! input/output
real(WP), intent(in) :: dt
- real(WP), dimension(N_proc(1), N_proc(2), N_proc(3)), intent(in) :: Vy
- real(WP), dimension(N_proc(1), N_proc(2), N_proc(3)), intent(inout) :: scal3D
+ real(WP), dimension(mesh_sc%N_proc(1), mesh_sc%N_proc(2), mesh_sc%N_proc(3)), intent(in) :: Vy
+ real(WP), dimension(mesh_sc%N_proc(1), mesh_sc%N_proc(2), mesh_sc%N_proc(3)), intent(inout) :: scal3D
! other local variables
integer :: i,k ! indice of the currend mesh point
integer, dimension(direction) :: ind_group ! indice of the currend group of line ((i,k) by default)
- real(WP), dimension(N_proc(direction)) :: p_pos_adim ! adimensionned particles position
- real(WP), dimension(N_proc(direction)) :: p_V ! particles velocity
+ real(WP), dimension(mesh_sc%N_proc(direction)) :: p_pos_adim ! adimensionned particles position
+ real(WP), dimension(mesh_sc%N_proc(direction)) :: p_V ! particles velocity
logical, dimension(bl_nb(direction)+1) :: bl_type ! is the particle block a center block or a left one ?
logical, dimension(bl_nb(direction)) :: bl_tag ! indice of tagged particles
ind_group = 0
- do k = 1, N_proc(3)
+ do k = 1, mesh_sc%N_proc(3)
ind_group(2) = ind_group(2) + 1
ind_group(1) = 0
- do i = 1, N_proc(1)
+ do i = 1, mesh_sc%N_proc(1)
ind_group(1) = ind_group(1) + 1
! ===== Init particles =====
@@ -98,7 +99,7 @@ subroutine advecY_calc_line(dt,Vy,scal3D)
! -- Compute velocity (with a RK2 scheme) --
call AC_particle_velocity_line(dt, direction, ind_group, p_pos_adim, p_V)
! -- Advec particles --
- p_pos_adim = p_pos_adim + dt*p_V/d_sc(direction)
+ p_pos_adim = p_pos_adim + dt*p_V/mesh_sc%dx(direction)
! ===== Remeshing =====
! -- Pre-Remeshing: Determine blocks type and tag particles --
@@ -142,7 +143,7 @@ subroutine Yremesh_O2_line(ind_group, p_pos_adim, bl_type, bl_tag,i,k,scal)
logical, dimension(:), intent(in) :: bl_type
logical, dimension(:), intent(in) :: bl_tag
real(WP), dimension(:), intent(in) :: p_pos_adim
- real(WP), dimension(N_proc(1), N_proc(2), N_proc(3)), intent(inout) :: scal
+ real(WP), dimension(mesh_sc%N_proc(1), mesh_sc%N_proc(2), mesh_sc%N_proc(3)), intent(inout) :: scal
! Other local variables
! Variable used to remesh particles in a buffer
real(WP),dimension(:),allocatable :: send_buffer ! buffer use to remesh the scalar before to send it
@@ -160,12 +161,12 @@ subroutine Yremesh_O2_line(ind_group, p_pos_adim, bl_type, bl_tag,i,k,scal)
! First particle is a left one
send_j_min = floor(p_pos_adim(1))-1
end if
- if (bl_type(N_proc(direction)/bl_size +1)) then
+ if (bl_type(mesh_sc%N_proc(direction)/bl_size +1)) then
! Last particle is a centered one
- send_j_max = nint(p_pos_adim(N_proc(direction)))+1
+ send_j_max = nint(p_pos_adim(mesh_sc%N_proc(direction)))+1
else
! Last particle is a left one
- send_j_max = floor(p_pos_adim(N_proc(direction)))+1
+ send_j_max = floor(p_pos_adim(mesh_sc%N_proc(direction)))+1
end if
! -- Determine the communication needed : who will communicate whit who ? (ie compute sender and recer) --
@@ -205,12 +206,12 @@ subroutine advecY_init_line(Vy, i, k, p_pos_adim, p_V)
! input/output
integer, intent(in) :: i,k
- real(WP), dimension(N_proc(direction)), intent(out) :: p_pos_adim, p_V
+ real(WP), dimension(mesh_sc%N_proc(direction)), intent(out) :: p_pos_adim, p_V
real(WP), dimension(:,:,:), intent(in) :: Vy
! Other local variables
integer :: ind ! indice
- do ind = 1, N_proc(direction)
+ do ind = 1, mesh_sc%N_proc(direction)
p_pos_adim(ind) = ind
p_V(ind) = Vy(i,ind,k)
end do
diff --git a/HySoP/src/scalesInterface/particles/advec_line/advecZ_line.f90 b/HySoP/src/scalesInterface/particles/advec_line/advecZ_line.f90
index 26d2e867a..3a5e0a014 100644
--- a/HySoP/src/scalesInterface/particles/advec_line/advecZ_line.f90
+++ b/HySoP/src/scalesInterface/particles/advec_line/advecZ_line.f90
@@ -1,3 +1,4 @@
+!USEFORTEST advec
!> @addtogroup part
!! @{
@@ -72,21 +73,21 @@ subroutine advecZ_calc_line(dt,Vz,scal3D)
! Input/Output
real(WP), intent(in) :: dt
- real(WP), dimension(N_proc(1), N_proc(2), N_proc(3)), intent(in) :: Vz
- real(WP), dimension(N_proc(1), N_proc(2), N_proc(3)), intent(inout) :: scal3D
+ real(WP), dimension(mesh_sc%N_proc(1), mesh_sc%N_proc(2), mesh_sc%N_proc(3)), intent(in) :: Vz
+ real(WP), dimension(mesh_sc%N_proc(1), mesh_sc%N_proc(2), mesh_sc%N_proc(3)), intent(inout) :: scal3D
! Other local variables
integer :: i,j ! indice of the currend mesh point
integer, dimension(2) :: ind_group ! indice of the currend group of line (=(i,k) by default)
- real(WP), dimension(N_proc(direction)) :: p_pos_adim ! adimensionned particles position
- real(WP), dimension(N_proc(direction)) :: p_V ! particles velocity
+ real(WP), dimension(mesh_sc%N_proc(direction)) :: p_pos_adim ! adimensionned particles position
+ real(WP), dimension(mesh_sc%N_proc(direction)) :: p_V ! particles velocity
logical, dimension(bl_nb(direction)+1) :: bl_type ! is the particle block a center block or a left one ?
logical, dimension(bl_nb(direction)) :: bl_tag ! indice of tagged particles
ind_group = 0
- do j = 1, N_proc(2)
+ do j = 1, mesh_sc%N_proc(2)
ind_group(2) = ind_group(2) + 1
ind_group(1) = 0
- do i = 1, N_proc(1)
+ do i = 1, mesh_sc%N_proc(1)
ind_group(1) = ind_group(1) + 1
! ===== Init particles =====
@@ -96,7 +97,7 @@ subroutine advecZ_calc_line(dt,Vz,scal3D)
! -- Compute velocity (with a RK2 scheme) --
call AC_particle_velocity_line(dt, direction, ind_group, p_pos_adim, p_V)
! -- Advec particles --
- p_pos_adim = p_pos_adim + dt*p_V/d_sc(direction)
+ p_pos_adim = p_pos_adim + dt*p_V/mesh_sc%dx(direction)
! ===== Remeshing =====
! -- Pre-Remeshing: Determine blocks type and tag particles --
@@ -140,7 +141,7 @@ subroutine Zremesh_O2_line(ind_group, p_pos_adim, bl_type, bl_tag,i,j,scal)
logical, dimension(:), intent(in) :: bl_type
logical, dimension(:), intent(in) :: bl_tag
real(WP), dimension(:), intent(in) :: p_pos_adim
- real(WP), dimension(N_proc(1), N_proc(2), N_proc(3)), intent(inout) :: scal
+ real(WP), dimension(mesh_sc%N_proc(1), mesh_sc%N_proc(2), mesh_sc%N_proc(3)), intent(inout) :: scal
! Other local variables
! Variable used to remesh particles in a buffer
real(WP),dimension(:),allocatable :: send_buffer ! buffer use to remesh the scalar before to send it
@@ -159,12 +160,12 @@ subroutine Zremesh_O2_line(ind_group, p_pos_adim, bl_type, bl_tag,i,j,scal)
! First particle is a left one
send_j_min = floor(p_pos_adim(1))-1
end if
- if (bl_type(N_proc(direction)/bl_size +1)) then
+ if (bl_type(mesh_sc%N_proc(direction)/bl_size +1)) then
! Last particle is a centered one
- send_j_max = nint(p_pos_adim(N_proc(direction)))+1
+ send_j_max = nint(p_pos_adim(mesh_sc%N_proc(direction)))+1
else
! Last particle is a left one
- send_j_max = floor(p_pos_adim(N_proc(direction)))+1
+ send_j_max = floor(p_pos_adim(mesh_sc%N_proc(direction)))+1
end if
! -- Determine the communication needed : who will communicate whit who ? (ie compute sender and recer) --
@@ -204,12 +205,12 @@ subroutine advecZ_init_line(Vz, i, j, p_pos_adim, p_V)
! Input/Output
integer, intent(in) :: i,j
- real(WP), dimension(N_proc(direction)), intent(out) :: p_pos_adim, p_V
+ real(WP), dimension(mesh_sc%N_proc(direction)), intent(out) :: p_pos_adim, p_V
real(WP), dimension(:,:,:), intent(in) :: Vz
! Other local variables
integer :: ind ! indice
- do ind = 1, N_proc(direction)
+ do ind = 1, mesh_sc%N_proc(direction)
p_pos_adim(ind) = ind
p_V(ind) = Vz(i,j,ind)
end do
diff --git a/HySoP/src/scalesInterface/particles/advec_line/advec_common_line.f90 b/HySoP/src/scalesInterface/particles/advec_line/advec_common_line.f90
index d33264d6a..7b4325354 100644
--- a/HySoP/src/scalesInterface/particles/advec_line/advec_common_line.f90
+++ b/HySoP/src/scalesInterface/particles/advec_line/advec_common_line.f90
@@ -1,3 +1,4 @@
+!USEFORTEST advec
!> @addtogroup part
!! @{
!------------------------------------------------------------------------------
@@ -39,7 +40,6 @@
module advec_common_line
use precision_tools
- use structure_tools
implicit none
@@ -105,11 +105,11 @@ subroutine AC_obtain_receivers_line(direction, ind_group, rece_ind_min, rece_ind
tag_min = 5
tag_max = 6
- send_gap = 3*N(direction)
+ send_gap = 3*mesh_sc%N(direction)
+
+ rece_gap(1) = floor(real(rece_ind_min-1)/mesh_sc%N_proc(direction))
+ rece_gap(2) = floor(real(rece_ind_max-1)/mesh_sc%N_proc(direction))
- rece_gap(1) = floor(real(rece_ind_min-1)/N_proc(direction))
- rece_gap(2) = floor(real(rece_ind_max-1)/N_proc(direction))
-
! ===== Communicate with my neigbors -> obtain ghost ! ====
! Compute their rank
call mpi_cart_shift(D_comm(direction), 0, 1, rankP, rankN, ierr)
@@ -152,10 +152,10 @@ subroutine AC_obtain_receivers_line(direction, ind_group, rece_ind_min, rece_ind
! ===== Receive information form the first and the last processus which need a part of my local velocity field =====
- if (send_gap(1) == 3*N(direction)) then
+ if (send_gap(1) == 3*mesh_sc%N(direction)) then
call mpi_recv(send_gap(1), 1, MPI_INTEGER, MPI_ANY_SOURCE, tag_table(1), D_comm(direction), statut, ierr)
end if
- if (send_gap(2) == 3*N(direction)) then
+ if (send_gap(2) == 3*mesh_sc%N(direction)) then
call mpi_recv(send_gap(2), 1, MPI_INTEGER, MPI_ANY_SOURCE, tag_table(2), D_comm(direction), statut, ierr)
end if
@@ -191,6 +191,7 @@ subroutine AC_particle_velocity_line(dt, direction, ind_group, p_pos_adim, p_V)
! a memory copy is still needed to send velocity field to other processus : mpi send contiguous memory values
use mpi
+ use structure_tools
use cart_topology ! info about mesh and mpi topology
use advec_variables ! contains info about solver parameters and others.
@@ -201,10 +202,10 @@ subroutine AC_particle_velocity_line(dt, direction, ind_group, p_pos_adim, p_V)
real(WP), dimension(:), intent(in) :: p_pos_adim
real(WP), dimension(:), intent(inout) :: p_V
! Others, local
- real(WP), dimension(N_proc(direction)) :: p_pos_bis ! adimensionned position of the middle point
- real(WP), dimension(N_proc(direction)), target :: p_V_bis ! velocity of the middle point
- real(WP), dimension(N_proc(direction)) :: weight ! interpolation weight
- type(real_pter), dimension(N_proc(direction)) :: Vp, Vm ! Velocity on previous and next mesh point
+ real(WP), dimension(mesh_sc%N_proc(direction)) :: p_pos_bis ! adimensionned position of the middle point
+ real(WP), dimension(mesh_sc%N_proc(direction)), target :: p_V_bis ! velocity of the middle point
+ real(WP), dimension(mesh_sc%N_proc(direction)) :: weight ! interpolation weight
+ type(real_pter), dimension(mesh_sc%N_proc(direction)) :: Vp, Vm ! Velocity on previous and next mesh point
real(WP), dimension(:), allocatable, target :: V_buffer ! Velocity buffer for postion outside of the local subdomain
integer :: size_buffer ! buffer size
integer :: rece_ind_min ! the minimal indice used in velocity interpolation
@@ -230,20 +231,20 @@ subroutine AC_particle_velocity_line(dt, direction, ind_group, p_pos_adim, p_V)
! -- Initialisation --
ind_com = 0
- do ind = 1, N_proc(direction)
+ do ind = 1, mesh_sc%N_proc(direction)
nullify(Vp(ind)%pter)
nullify(Vm(ind)%pter)
end do
! Compute the midlle point
- p_pos_bis = p_pos_adim + (dt/2.0)*p_V/d_sc(direction)
+ p_pos_bis = p_pos_adim + (dt/2.0)*p_V/mesh_sc%dx(direction)
p_V_bis = p_V
! Compute range of the set of point where I need the velocity value
rece_ind_min = floor(p_pos_bis(1))
- rece_ind_max = floor(p_pos_bis(N_proc(direction))) + 1
+ rece_ind_max = floor(p_pos_bis(mesh_sc%N_proc(direction))) + 1
! Allocate the buffer
- ! If rece_ind_min and rece_ind_max are not in [N_proc(direction);1] then it will change the number of communication
- ! size_buffer = max(temp - N_proc(direction), 0) - min(0, temp)
- !size_buffer = - max(temp - N_proc(direction), 0) - min(0, temp)
+ ! If rece_ind_min and rece_ind_max are not in [mesh_sc%N_proc(direction);1] then it will change the number of communication
+ ! size_buffer = max(temp - mesh_sc%N_proc(direction), 0) - min(0, temp)
+ !size_buffer = - max(temp - mesh_sc%N_proc(direction), 0) - min(0, temp)
! It must work, but for first test we prefer compute size_buffer more simply
size_buffer = 0
@@ -257,10 +258,10 @@ subroutine AC_particle_velocity_line(dt, direction, ind_group, p_pos_adim, p_V)
call mpi_cart_shift(D_comm(direction), 0, proc_gap, rankP, rece_rank(proc_gap), ierr)
if (rece_rank(proc_gap) /= D_rank(direction)) then
! Range I want
- gap = proc_gap*N_proc(direction)
+ gap = proc_gap*mesh_sc%N_proc(direction)
rece_range(1) = max(rece_ind_min, gap+1) ! fortran => indice start from 0
- rece_range(2) = min(rece_ind_max, gap+N_proc(direction))
- ! Tag = concatenation of (rank+1), ind_group(1), ind_group(2), direction et unique Id.
+ rece_range(2) = min(rece_ind_max, gap+mesh_sc%N_proc(direction))
+ ! Tag = concatenation of (rank+1), ind_group(1), ind_group(2), direction et unique Id.
tag = compute_tag(ind_group, tag_velo_range, direction, proc_gap)
! Send message
size_buffer = size_buffer + (rece_range(2)-rece_range(1)) + 1
@@ -280,7 +281,7 @@ subroutine AC_particle_velocity_line(dt, direction, ind_group, p_pos_adim, p_V)
tag = compute_tag(ind_group, tag_velo_range, direction, -proc_gap)
! Ib - Receive the message
call mpi_recv(send_range(1), 2, MPI_INTEGER, send_rank, tag, D_comm(direction), rece_status, ierr)
- send_range = send_range + proc_gap*N_proc(direction)
+ send_range = send_range + proc_gap*mesh_sc%N_proc(direction)
! II - Send it
! IIa - Compute send tag
tag = compute_tag(ind_group, tag_velo_V, direction, proc_gap)
@@ -298,9 +299,9 @@ subroutine AC_particle_velocity_line(dt, direction, ind_group, p_pos_adim, p_V)
! IIa - Compute reception tag
tag = compute_tag(ind_group, tag_velo_V, direction, -proc_gap)
! IIb - Receive message
- gap = proc_gap*N_proc(direction)
+ gap = proc_gap*mesh_sc%N_proc(direction)
rece_range(1) = max(rece_ind_min, gap+1) ! fortran => indice start from 0
- rece_range(2) = min(rece_ind_max, gap+N_proc(direction))
+ rece_range(2) = min(rece_ind_max, gap+mesh_sc%N_proc(direction))
msg_size = rece_range(2)-rece_range(1)+1
call mpi_Irecv(V_buffer(ind), msg_size, MPI_DOUBLE_PRECISION, rece_rank(proc_gap), tag, D_comm(direction), &
& rece_request(proc_gap), ierr)
@@ -315,19 +316,19 @@ subroutine AC_particle_velocity_line(dt, direction, ind_group, p_pos_adim, p_V)
pos = floor(p_pos_bis(ind))
weight(ind) = p_pos_bis(ind)-pos
! Vm = V(pos)
- proc_gap = floor(real(pos-1)/N_proc(direction))
+ proc_gap = floor(real(pos-1)/mesh_sc%N_proc(direction))
call mpi_cart_shift(D_comm(direction), 0, proc_gap, rankP, send_rank, ierr)
if (send_rank == D_rank(direction)) then
- Vm(ind)%pter => p_V_bis(pos-proc_gap*N_proc(direction))
+ Vm(ind)%pter => p_V_bis(pos-proc_gap*mesh_sc%N_proc(direction))
else
ind_com = ind_com + 1
Vm(ind)%pter => V_buffer(ind_com)
end if
! Vp = V(pos+1)
- proc_gap = floor(real(pos+1-1)/N_proc(direction))
+ proc_gap = floor(real(pos+1-1)/mesh_sc%N_proc(direction))
call mpi_cart_shift(D_comm(direction), 0, proc_gap, rankP, send_rank, ierr)
if (send_rank == D_rank(direction)) then
- Vp(ind)%pter => p_V_bis(pos+1-proc_gap*N_proc(direction))
+ Vp(ind)%pter => p_V_bis(pos+1-proc_gap*mesh_sc%N_proc(direction))
else
ind_com = ind_com + 1
Vp(ind)%pter => V_buffer(ind_com)
@@ -335,7 +336,7 @@ subroutine AC_particle_velocity_line(dt, direction, ind_group, p_pos_adim, p_V)
pos_old = pos
! Following indice : we use previous work (already done)
- do ind = 2, N_proc(direction)
+ do ind = 2, mesh_sc%N_proc(direction)
pos = floor(p_pos_bis(ind))
weight(ind) = p_pos_bis(ind)-pos
select case(pos-pos_old)
@@ -347,10 +348,10 @@ subroutine AC_particle_velocity_line(dt, direction, ind_group, p_pos_adim, p_V)
! The particle follows the previous one
Vm(ind)%pter => Vp(ind-1)%pter
! Vp = V(pos+1)
- proc_gap = floor(real(pos+1-1)/N_proc(direction)) ! fortran -> indice starts from 1
+ proc_gap = floor(real(pos+1-1)/mesh_sc%N_proc(direction)) ! fortran -> indice starts from 1
call mpi_cart_shift(D_comm(direction), 0, proc_gap, rankP, send_rank, ierr)
if (send_rank == D_rank(direction)) then
- Vp(ind)%pter => p_V_bis(pos+1-proc_gap*N_proc(direction))
+ Vp(ind)%pter => p_V_bis(pos+1-proc_gap*mesh_sc%N_proc(direction))
else
ind_com = ind_com + 1
Vp(ind)%pter => V_buffer(ind_com)
@@ -358,19 +359,19 @@ subroutine AC_particle_velocity_line(dt, direction, ind_group, p_pos_adim, p_V)
case(2)
! pos = pos_old +2, wich correspond to "extention"
! Vm = V(pos)
- proc_gap = floor(real(pos-1)/N_proc(direction))
+ proc_gap = floor(real(pos-1)/mesh_sc%N_proc(direction))
call mpi_cart_shift(D_comm(direction), 0, proc_gap, rankP, send_rank, ierr)
if (send_rank == D_rank(direction)) then
- Vm(ind)%pter => p_V_bis(pos-proc_gap*N_proc(direction))
+ Vm(ind)%pter => p_V_bis(pos-proc_gap*mesh_sc%N_proc(direction))
else
ind_com = ind_com + 1
Vm(ind)%pter => V_buffer(ind_com)
end if
! Vp = V(pos+1)
- proc_gap = floor(real(pos+1-1)/N_proc(direction))
+ proc_gap = floor(real(pos+1-1)/mesh_sc%N_proc(direction))
call mpi_cart_shift(D_comm(direction), 0, proc_gap, rankP, send_rank, ierr)
if (send_rank == D_rank(direction)) then
- Vp(ind)%pter => p_V_bis(pos+1-proc_gap*N_proc(direction))
+ Vp(ind)%pter => p_V_bis(pos+1-proc_gap*mesh_sc%N_proc(direction))
else
ind_com = ind_com + 1
Vp(ind)%pter => V_buffer(ind_com)
@@ -391,11 +392,11 @@ subroutine AC_particle_velocity_line(dt, direction, ind_group, p_pos_adim, p_V)
! Then compute the field
- do ind = 1, N_proc(direction)
+ do ind = 1, mesh_sc%N_proc(direction)
p_V(ind) = weight(ind)*Vp(ind)%pter + (1-weight(ind))*Vm(ind)%pter
end do
- do ind = 1, N_proc(direction)
+ do ind = 1, mesh_sc%N_proc(direction)
nullify(Vp(ind)%pter)
nullify(Vm(ind)%pter)
end do
@@ -457,11 +458,11 @@ subroutine AC_type_and_block_line(dt, direction, ind_group, p_V, &
logical, dimension(bl_nb(direction)+1), intent(out) :: bl_type ! is the particle block a center block or a left one ?
logical, dimension(bl_nb(direction)), intent(out) :: bl_tag ! indice of tagged particles
! Local variables
- real(WP), dimension(bl_nb(direction)+1) :: bl_lambdaMin ! for a particle, lamda = V*dt/dx ; bl_lambdaMin = min of
+ real(WP), dimension(bl_nb(direction)+1) :: bl_lambdaMin ! for a particle, lamda = V*dt/dx ; bl_lambdaMin = min of
! lambda on a block (take also into account first following particle)
real(WP) :: lambP, lambN ! buffer to exchange some lambda min with other processus
integer, dimension(bl_nb(direction)+1) :: bl_ind ! block index : integer as lambda in (bl_ind,bl_ind+1) for a left block
- ! and lambda in (bl_ind-1/2, bl_ind+1/2) for a right block
+ ! and lambda in (bl_ind-1/2, bl_ind+1/2) for a right block
integer :: ind, i_p ! some indices
real(WP) :: cfl ! = d_sc
integer :: rankP, rankN ! processus rank for shift (P= previous, N = next)
@@ -473,7 +474,7 @@ subroutine AC_type_and_block_line(dt, direction, ind_group, p_V, &
integer :: ierr ! mpi error code
! ===== Initialisation =====
- cfl = dt/d_sc(direction)
+ cfl = dt/mesh_sc%dx(direction)
! ===== Compute bl_lambdaMin =====
! -- Compute rank of my neighbor --
@@ -491,7 +492,7 @@ subroutine AC_type_and_block_line(dt, direction, ind_group, p_V, &
! -- For the last block (1/2) --
! The processus contains only its first half => exchange ghost with the next processus
ind = bl_nb(direction) + 1
- bl_lambdaMin(ind) = minval(p_V(N_proc(direction)-(bl_size/2)+1:N_proc(direction)))*cfl
+ bl_lambdaMin(ind) = minval(p_V(mesh_sc%N_proc(direction)-(bl_size/2)+1:mesh_sc%N_proc(direction)))*cfl
! Send message
call mpi_Isend(bl_lambdaMin(ind), 1, MPI_DOUBLE_PRECISION, rankN, tag_table(2), D_comm(direction), send_request(2), ierr)
! Receive it
@@ -587,14 +588,14 @@ subroutine AC_obtain_senders_line(send_i_min, send_i_max, direction, ind_group,
tag_table = compute_tag(ind_group, tag_obtsend_NP, direction)
- rece_proc = 3*N(direction)
+ rece_proc = 3*mesh_sc%N(direction)
- proc_min = floor(real(send_i_min-1)/N_proc(direction))
- proc_max = floor(real(send_i_max-1)/N_proc(direction))
+ proc_min = floor(real(send_i_min-1)/mesh_sc%N_proc(direction))
+ proc_max = floor(real(send_i_max-1)/mesh_sc%N_proc(direction))
allocate(send_request(proc_min:proc_max,3))
send_request(:,3) = 0
-
+
! Send
do proc_gap = proc_min, proc_max
! Compute the rank of the target processus
@@ -602,10 +603,10 @@ subroutine AC_obtain_senders_line(send_i_min, send_i_max, direction, ind_group,
! Determine if I am the the first or the last processes (considering my
! coordinate along the current directory) to send information to
! one of these processes.
- ! Note that local indice go from 1 to N_proc (fortran).
+ ! Note that local indice go from 1 to mesh_sc%N_proc (fortran).
! I am the first ?
- if ((send_i_min< +1-2*bl_bound_size + proc_gap*N_proc(direction)+1).AND. &
- & (send_i_max>= proc_gap*N_proc(direction))) then
+ if ((send_i_min< +1-2*bl_bound_size + proc_gap*mesh_sc%N_proc(direction)+1).AND. &
+ & (send_i_max>= proc_gap*mesh_sc%N_proc(direction))) then
if(rankN /= D_rank(direction)) then
call mpi_Isend(-proc_gap, 1, MPI_INTEGER, rankN, tag_table(1), D_comm(direction), &
& send_request(proc_gap,1), ierr)
@@ -615,8 +616,8 @@ subroutine AC_obtain_senders_line(send_i_min, send_i_max, direction, ind_group,
end if
end if
! I am the last ?
- if ((send_i_max > -1+2*bl_bound_size + (proc_gap+1)*N_proc(direction)) &
- & .AND.(send_i_min<= (proc_gap+1)*N_proc(direction))) then
+ if ((send_i_max > -1+2*bl_bound_size + (proc_gap+1)*mesh_sc%N_proc(direction)) &
+ & .AND.(send_i_min<= (proc_gap+1)*mesh_sc%N_proc(direction))) then
if(rankN /= D_rank(direction)) then
call mpi_Isend(-proc_gap, 1, MPI_INTEGER, rankN, tag_table(2), D_comm(direction), &
& send_request(proc_gap,2), ierr)
@@ -629,10 +630,10 @@ subroutine AC_obtain_senders_line(send_i_min, send_i_max, direction, ind_group,
! Receive
- if (rece_proc(1) == 3*N(direction)) then
+ if (rece_proc(1) == 3*mesh_sc%N(direction)) then
call mpi_recv(rece_proc(1), 1, MPI_INTEGER, MPI_ANY_SOURCE, tag_table(1), D_comm(direction), statut, ierr)
end if
- if (rece_proc(2) == 3*N(direction)) then
+ if (rece_proc(2) == 3**mesh_sc%N(direction)) then
call mpi_recv(rece_proc(2), 1, MPI_INTEGER, MPI_ANY_SOURCE, tag_table(2), D_comm(direction), statut, ierr)
end if
@@ -684,7 +685,7 @@ subroutine AC_bufferToScalar_line(direction, ind_group, send_i_min, send_i_max,
integer, intent(in) :: proc_min ! smaller gap between me and the processes to where I send data
integer, intent(in) :: proc_max ! smaller gap between me and the processes to where I send data
real(WP), dimension(send_i_min:send_i_max), intent(in) :: send_buffer
- real(WP), dimension(N_proc(direction)), intent(inout) :: scal1D
+ real(WP), dimension(mesh_sc%N_proc(direction)), intent(inout) :: scal1D
! Variables used to communicate between subdomains. A variable prefixed by "send_"(resp "rece")
! design something I send (resp. I receive).
@@ -730,9 +731,9 @@ subroutine AC_bufferToScalar_line(direction, ind_group, send_i_min, send_i_max,
do proc_gap = proc_min, proc_max
! Compute the rank of the target processus
call mpi_cart_shift(D_comm(direction), 0, proc_gap, rankP, send_rank, ierr)
- send_gap = proc_gap*N_proc(direction)
+ send_gap = proc_gap*mesh_sc%N_proc(direction)
send_range(1, proc_gap) = max(send_i_min, send_gap+1) ! fortran => indice start from 0
- send_range(2, proc_gap) = min(send_i_max, send_gap+N_proc(direction))
+ send_range(2, proc_gap) = min(send_i_max, send_gap+mesh_sc%N_proc(direction))
if (send_rank/=D_rank(direction)) then
! Determine quantity of information to send
comm_size = send_range(2, proc_gap)-send_range(1, proc_gap)+1
@@ -774,18 +775,18 @@ subroutine AC_bufferToScalar_line(direction, ind_group, send_i_min, send_i_max,
comm_size=(rece_i_max-rece_i_min+1)
allocate(rece_buffer(rece_i_min:rece_i_max)) ! XXX possible optimisation
! by allocating one time to the max size, note that the range use in
- ! this allocation instruction is include in (1, N_proc(2))
+ ! this allocation instruction is include in (1, mesh_sc%N_proc(2))
tag = compute_tag(ind_group, tag_bufToScal_buffer, direction, -proc_gap)
call mpi_recv(rece_buffer(rece_i_min), comm_size, MPI_DOUBLE_PRECISION, &
& rece_rank(proc_gap), tag, D_comm(direction), rece_status, ierr)
! Update the scalar field
- send_gap = proc_gap*N_proc(direction)
+ send_gap = proc_gap*mesh_sc%N_proc(direction)
scal1D(rece_i_min+send_gap:rece_i_max+send_gap) = scal1D(rece_i_min+send_gap:rece_i_max+send_gap) &
& + rece_buffer(rece_i_min:rece_i_max)
deallocate(rece_buffer)
end if
end do
-
+
! Free Isend buffer
do proc_gap = proc_min, proc_max
diff --git a/HySoP/src/scalesInterface/particles/advec_line/advec_remesh_formula_line.f90 b/HySoP/src/scalesInterface/particles/advec_line/advec_remesh_formula_line.f90
index 6afd42c7c..2cca2e257 100644
--- a/HySoP/src/scalesInterface/particles/advec_line/advec_remesh_formula_line.f90
+++ b/HySoP/src/scalesInterface/particles/advec_line/advec_remesh_formula_line.f90
@@ -1,3 +1,4 @@
+!USEFORTEST advec
!> @addtogroup part
!! @{
!------------------------------------------------------------------------------
@@ -63,7 +64,7 @@ subroutine AC_remesh_lambda2corrected_basic(direction, p_pos_adim, scal1D, bl_ty
! Input/Output
integer, intent(in) :: direction
real(WP), dimension(:), intent(in) :: p_pos_adim
- real(WP), dimension(N_proc(direction)), intent(in) :: scal1D
+ real(WP), dimension(mesh_sc%N_proc(direction)), intent(in) :: scal1D
logical, dimension(:), intent(in) :: bl_type
logical, dimension(:), intent(in) :: bl_tag
integer, intent(in) :: ind_min, ind_max
@@ -75,7 +76,7 @@ subroutine AC_remesh_lambda2corrected_basic(direction, p_pos_adim, scal1D, bl_ty
send_j_min = ind_min
send_j_max = ind_max
- do p_ind = 1, N_proc(direction), bl_size
+ do p_ind = 1, mesh_sc%N_proc(direction), bl_size
bl_ind = p_ind/bl_size + 1
if (bl_tag(bl_ind)) then
! Tag case
@@ -142,7 +143,7 @@ subroutine AC_remesh_lambda4corrected_basic(direction, p_pos_adim, scal1D, bl_ty
! Input/Output
integer, intent(in) :: direction
real(WP), dimension(:), intent(in) :: p_pos_adim
- real(WP), dimension(N_proc(direction)), intent(in) :: scal1D
+ real(WP), dimension(mesh_sc%N_proc(direction)), intent(in) :: scal1D
logical, dimension(:), intent(in) :: bl_type
logical, dimension(:), intent(in) :: bl_tag
integer, intent(in) :: ind_min, ind_max
@@ -154,7 +155,7 @@ subroutine AC_remesh_lambda4corrected_basic(direction, p_pos_adim, scal1D, bl_ty
send_j_min = ind_min
send_j_max = ind_max
- do p_ind = 1, N_proc(direction), bl_size
+ do p_ind = 1, mesh_sc%N_proc(direction), bl_size
bl_ind = p_ind/bl_size + 1
if (bl_tag(bl_ind)) then
! Tagged case
diff --git a/HySoP/src/scalesInterface/particles/advec_line/advec_remesh_line.F90 b/HySoP/src/scalesInterface/particles/advec_line/advec_remesh_line.F90
new file mode 100644
index 000000000..6331ae32d
--- /dev/null
+++ b/HySoP/src/scalesInterface/particles/advec_line/advec_remesh_line.F90
@@ -0,0 +1,1312 @@
+!USEFORTEST advec
+!> @addtogroup part
+
+!------------------------------------------------------------------------------
+!
+! MODULE: advec_remesh_line
+!
+!
+! DESCRIPTION:
+!> The module advec_remesh_line contains different semi-optimized remeshing
+!! procedure. They are here for debugging/test/comparaison purpose and will
+!! be deleted in "not to far" future (after adding optimized M'6, having a lot
+!! of validation and having performed benchmarks).
+!
+!! The module "test_advec" can be used in order to validate the procedures
+!! embedded in this module.
+!
+!> @author
+!! Jean-Baptiste Lagaert, LEGI
+!
+!------------------------------------------------------------------------------
+
+module advec_remesh_line
+
+ use precision_tools
+ use advec_abstract_proc
+ use advec_correction
+
+ implicit none
+
+ ! ===== Public procedures =====
+ !----- (corrected) lambda 2 Remeshing method -----
+ public :: Xremesh_O2 ! order 2
+ public :: Yremesh_O2 ! order 2
+ public :: Zremesh_O2 ! order 2
+ !----- (corrected) lambda 4 Remeshing method -----
+ public :: Xremesh_O4 ! order 4
+ public :: Yremesh_O4 ! order 4
+ public :: Zremesh_O4 ! order 4
+ !----- M'6 remeshing method -----
+ public :: Xremesh_Mprime6
+ public :: Yremesh_Mprime6
+ public :: Zremesh_Mprime6
+
+
+ ! ===== Private variable ====
+
+contains
+
+! #####################################################################################
+! ##### #####
+! ##### Public procedure #####
+! ##### #####
+! #####################################################################################
+
+! ============================================================================
+! ==================== Remeshing along X subroutines ====================
+! ============================================================================
+
+!> remeshing along Xwith an order 2 method, corrected to allow large CFL number - group version
+!! @param[in] direction = current direction
+!! @param[in] ind_group = coordinate of the current group of lines
+!! @param[in] gs = size of groups (along X direction)
+!! @param[in] p_pos_adim = adimensionned particles position
+!! @param[in] p_V = particles velocity (needed for tag and type)
+!! @param[in] j,k = indice of of the current line (x-coordinate and z-coordinate)
+!! @param[in] dt = time step (needed for tag and type)
+!! @param[in,out] scal = scalar field to advect
+subroutine Xremesh_O2(direction, ind_group, gs, p_pos_adim, p_V, j, k, scal, dt)
+
+ use advec_common ! Some procedures common to advection along all directions
+ use advec_common_line ! Some procedures common to advection along all directions
+ use advec_remeshing_line ! Remeshing formula
+ use advec_variables ! contains info about solver parameters and others.
+ use cart_topology ! Description of mesh and of mpi topology
+
+ ! Input/Output
+ integer, intent(in) :: direction
+ integer, dimension(2), intent(in) :: ind_group
+ integer, dimension(2), intent(in) :: gs
+ integer, intent(in) :: j, k
+ real(WP), dimension(:,:,:), intent(in) :: p_pos_adim ! adimensionned particles position
+ real(WP), dimension(:,:,:), intent(in) :: p_V ! particles velocity
+ real(WP), dimension(:,:,:), intent(inout) :: scal
+ real(WP), intent(in) :: dt
+ ! Other local variables
+ ! To type and tag particles
+ logical, dimension(bl_nb(direction)+1,gs(1),gs(2)) :: bl_type ! is the particle block a center block or a left one ?
+ logical, dimension(bl_nb(direction),gs(1),gs(2)) :: bl_tag ! indice of tagged particles
+ ! To compute recquired communications
+ integer, dimension(gs(1), gs(2)) :: send_group_min ! distance between me and processus wich send me information
+ integer, dimension(gs(1), gs(2)) :: send_group_max ! distance between me and processus wich send me information
+ ! Variable used to remesh particles in a buffer
+ real(WP),dimension(:),allocatable, target:: send_buffer ! buffer use to remesh the scalar before to send it to the right subdomain
+ ! sorted by receivers and not by coordinate.
+ integer, dimension(2,gs(1),gs(2)) :: rece_proc ! minimal and maximal gap between my Y-coordinate and the one from which
+ ! I will receive data
+ integer, dimension(gs(1),gs(2)) :: proc_min ! smaller gap between me and the processes to where I send data
+ integer, dimension(gs(1),gs(2)) :: proc_max ! smaller gap between me and the processes to where I send data
+
+ integer :: i1, i2 ! indice of a line into the group
+
+ ! -- Pre-Remeshing: Determine blocks type and tag particles --
+ call AC_type_and_block_group(dt, direction, gs, ind_group, p_V, bl_type, bl_tag)
+
+ ! -- Compute ranges --
+ where (bl_type(1,:,:))
+ ! First particle is a centered one
+ send_group_min = nint(p_pos_adim(1,:,:))-1
+ elsewhere
+ ! First particle is a left one
+ send_group_min = floor(p_pos_adim(1,:,:))-1
+ end where
+ where (bl_type(mesh_sc%N_proc(direction)/bl_size +1,:,:))
+ ! Last particle is a centered one
+ send_group_max = nint(p_pos_adim(mesh_sc%N_proc(direction),:,:))+1
+ elsewhere
+ ! Last particle is a left one
+ send_group_max = floor(p_pos_adim(mesh_sc%N_proc(direction),:,:))+1
+ end where
+
+ ! -- Determine the communication needed : who will communicate whit who ? (ie compute sender and recer) --
+ call AC_obtain_senders_group(direction, gs, ind_group, &
+ & send_group_min, send_group_max, proc_min, proc_max, rece_proc)
+
+ do i2 = 1, gs(2)
+ do i1 = 1, gs(1)
+ send_j_min = send_group_min(i1,i2)
+ send_j_max = send_group_max(i1,i2)
+
+ ! -- Allocate buffer for remeshing of local particles --
+ allocate(send_buffer(send_j_min:send_j_max))
+ send_buffer = 0.0;
+
+ ! -- Remesh the particles in the buffer --
+ call AC_remesh_lambda2corrected_basic(direction, p_pos_adim(:,i1,i2), scal(:,j+i1-1,k+i2-1), &
+ & bl_type(:,i1,i2), bl_tag(:,i1,i2), send_j_min, send_j_max, send_buffer)
+
+ ! -- Send the buffer to the matching processus and update the scalar field --
+ scal(:,j+i1-1,k+i2-1) = 0
+ call AC_bufferToScalar_line(direction, ind_group , send_j_min, send_j_max, proc_min(i1,i2), proc_max(i1,i2), &
+ & rece_proc(:,i1,i2), send_buffer, scal(:,j+i1-1,k+i2-1))
+
+ ! Deallocate all field
+ deallocate(send_buffer)
+
+ end do
+ end do
+
+end subroutine Xremesh_O2
+
+
+!> remeshing along X with an order 4 method, corrected to allow large CFL number - untagged particles
+!! @param[in] direction = current direction
+!! @param[in] ind_group = coordinate of the current group of lines
+!! @param[in] gs = size of groups (along X direction)
+!! @param[in] p_pos_adim = adimensionned particles position
+!! @param[in] p_V = particles velocity (needed for tag and type)
+!! @param[in] j,k = indice of of the current line (x-coordinate and z-coordinate)
+!! @param[in,out] scal = scalar field to advect
+!! @param[in] dt = time step (needed for tag and type)
+subroutine Xremesh_O4(direction, ind_group, gs, p_pos_adim, p_V, j,k, scal, dt)
+
+ use advec_common ! Some procedures common to advection along all directions
+ use advec_common_line ! Some procedures common to advection along all directions
+ use advec_remeshing_line ! Remeshing formula
+ use advec_variables ! contains info about solver parameters and others.
+ use cart_topology ! Description of mesh and of mpi topology
+
+ ! Input/Output
+ integer, intent(in) :: direction
+ integer, dimension(2), intent(in) :: ind_group
+ integer, dimension(2), intent(in) :: gs
+ integer, intent(in) :: j, k
+ real(WP), dimension(:,:,:), intent(in) :: p_pos_adim ! adimensionned particles position
+ real(WP), dimension(:,:,:), intent(in) :: p_V ! particles velocity
+ real(WP), dimension(:,:,:), intent(inout) :: scal
+ real(WP), intent(in) :: dt
+ ! Other local variables
+ ! To compute recquired communications
+ integer, dimension(gs(1), gs(2)) :: send_group_min ! distance between me and processus wich send me information
+ integer, dimension(gs(1), gs(2)) :: send_group_max ! distance between me and processus wich send me information
+ ! To type and tag particles
+ logical, dimension(bl_nb(direction)+1,gs(1),gs(2)) :: bl_type ! is the particle block a center block or a left one ?
+ logical, dimension(bl_nb(direction),gs(1),gs(2)) :: bl_tag ! indice of tagged particles
+ ! Variables used to remesh particles ...
+ ! ... and to communicate between subdomains. A variable prefixed by "send_"(resp "rece")
+ ! designes something I send (resp. I receive).
+ real(WP),dimension(:),allocatable :: send_buffer ! buffer use to remesh the scalar before to send it to the right subdomain
+ integer, dimension(2,gs(1),gs(2)) :: rece_proc ! minimal and maximal gap between my Y-coordinate and the one from which
+ ! I will receive data
+ integer, dimension(gs(1),gs(2)) :: proc_min ! smaller gap between me and the processes to where I send data
+ integer, dimension(gs(1),gs(2)) :: proc_max ! smaller gap between me and the processes to where I send data
+ integer :: i1, i2 ! indice of a line into the group
+
+ ! -- Pre-Remeshing: Determine blocks type and tag particles --
+ call AC_type_and_block_group(dt, direction, gs, ind_group, p_V, bl_type, bl_tag)
+
+ ! -- Compute ranges --
+ where (bl_type(1,:,:))
+ ! First particle is a centered one
+ send_group_min = nint(p_pos_adim(1,:,:))-2
+ elsewhere
+ ! First particle is a left one
+ send_group_min = floor(p_pos_adim(1,:,:))-2
+ end where
+ where (bl_type(mesh_sc%N_proc(direction)/bl_size +1,:,:))
+ ! Last particle is a centered one
+ send_group_max = nint(p_pos_adim(mesh_sc%N_proc(direction),:,:))+2
+ elsewhere
+ ! Last particle is a left one
+ send_group_max = floor(p_pos_adim(mesh_sc%N_proc(direction),:,:))+2
+ end where
+
+ ! -- Determine the communication needed : who will communicate whit who ? (ie compute sender and recer) --
+ call AC_obtain_senders_com(direction, gs, ind_group, send_group_min, &
+ & send_group_max, proc_min, proc_max, rece_proc)
+
+ do i2 = 1, gs(2)
+ do i1 = 1, gs(1)
+ send_j_min = send_group_min(i1,i2)
+ send_j_max = send_group_max(i1,i2)
+
+ ! -- Allocate buffer for remeshing of local particles --
+ allocate(send_buffer(send_j_min:send_j_max))
+ send_buffer = 0.0;
+
+ ! -- Remesh the particles in the buffer --
+ call AC_remesh_lambda4corrected_basic(direction, p_pos_adim(:,i1,i2), scal(:,j+i1-1,k+i2-1), &
+ & bl_type(:,i1,i2), bl_tag(:,i1,i2), send_j_min, send_j_max, send_buffer)
+
+ ! -- Send the buffer to the matching processus and update the scalar field --
+ scal(:,j+i1-1,k+i2-1) = 0
+ call AC_bufferToScalar_line(direction, ind_group, send_j_min, send_j_max, proc_min(i1,i2), proc_max(i1,i2), &
+ & rece_proc(:,i1,i2), send_buffer, scal(:,j+i1-1,k+i2-1))
+
+ ! Deallocate all field
+ deallocate(send_buffer)
+
+ end do
+ end do
+
+end subroutine Xremesh_O4
+
+
+!> remeshing along X with M'6 formula - No tag neither correction for large time steps.
+!! @param[in] direction = current direction
+!! @param[in] ind_group = coordinate of the current group of lines
+!! @param[in] gs = size of groups (along X direction)
+!! @param[in] p_pos_adim = adimensionned particles position
+!! @param[in] p_V = particles velocity (only to have the same profile
+!! then other remeshing procedures)
+!! @param[in] j,k = indice of of the current line (y-coordinate and z-coordinate)
+!! @param[in] dt = time step (only to have the same profile
+!! then other remeshing procedures)
+!! @param[in,out] scal = scalar field to advect
+subroutine Xremesh_Mprime6(direction, ind_group, gs, p_pos_adim, p_V, j,k,scal, dt)
+
+ use advec_common_line ! Some procedures common to advection along all directions
+ use advec_remeshing_line ! Remeshing formula
+ use advec_variables ! contains info about solver parameters and others.
+ use cart_topology ! Description of mesh and of mpi topology
+
+ ! Input/outpu
+ integer, intent(in) :: direction
+ integer, dimension(2), intent(in) :: ind_group
+ integer, dimension(2), intent(in) :: gs
+ integer, intent(in) :: j, k
+ real(WP), dimension(:,:,:), intent(in) :: p_pos_adim ! adimensionned particles position
+ real(WP), dimension(:,:,:), intent(in) :: p_V ! particles velocity
+ real(WP), dimension(:,:,:), intent(inout) :: scal
+ real(WP), intent(in) :: dt
+ ! Other local variables
+ ! To compute recquired communications
+ integer, dimension(gs(1), gs(2)) :: send_group_min ! distance between me and processus wich send me information
+ integer, dimension(gs(1), gs(2)) :: send_group_max ! distance between me and processus wich send me information
+ ! Variables used to remesh particles ...
+ ! ... and to communicate between subdomains. A variable prefixed by "send_"(resp "rece")
+ ! designes something I send (resp. I receive).
+ real(WP),dimension(:),allocatable :: send_buffer ! buffer use to remesh the scalar before to send it to the right subdomain
+ integer, dimension(2,gs(1),gs(2)) :: rece_proc ! minimal and maximal gap between my Y-coordinate and the one from which
+ ! I will receive data
+ integer, dimension(gs(1),gs(2)) :: proc_min ! smaller gap between me and the processes to where I send data
+ integer, dimension(gs(1),gs(2)) :: proc_max ! smaller gap between me and the processes to where I send data
+ integer :: i1, i2 ! indice of a line into the group
+ integer :: i ! indice of the current particle
+
+ ! -- Compute the remeshing domain --
+ send_group_min = floor(p_pos_adim(1,:,:)-2)
+ send_group_max = floor(p_pos_adim(mesh_sc%N_proc(direction),:,:)+3)
+
+ ! -- Determine the communication needed : who will communicate whit who ? (ie compute sender and recer) --
+ call AC_obtain_senders_com(direction, gs, ind_group, send_group_min, &
+ & send_group_max, proc_min, proc_max, rece_proc)
+
+ do i2 = 1, gs(2)
+ do i1 = 1, gs(1)
+ send_j_min = send_group_min(i1,i2)
+ send_j_max = send_group_max(i1,i2)
+
+ ! -- Allocate buffer for remeshing of local particles --
+ allocate(send_buffer(send_j_min:send_j_max))
+ send_buffer = 0.0;
+
+ ! -- Remesh the particles in the buffer --
+ do i = 1, mesh_sc%N_proc(direction), 1
+ call AC_remesh_Mprime6(p_pos_adim(i,i1,i2),scal(i,j+i1-1,k+i2-1), send_buffer)
+ end do
+
+ ! -- Send the buffer to the matching processus and update the scalar field --
+ scal(:,j+i1-1,k+i2-1) = 0
+ call AC_bufferToScalar_line(direction, ind_group, send_j_min, send_j_max, proc_min(i1,i2), proc_max(i1,i2), &
+ & rece_proc(:,i1,i2), send_buffer, scal(:,j+i1-1,k+i2-1))
+
+ ! Deallocate all field
+ deallocate(send_buffer)
+
+ end do
+ end do
+
+end subroutine Xremesh_Mprime6
+
+
+! ============================================================================
+! ==================== Remeshing along Y subroutines ====================
+! ============================================================================
+
+!> remeshing along Y with an order 2 method, corrected to allow large CFL number - group version
+!! @param[in] direction = current direction
+!! @param[in] ind_group = coordinate of the current group of lines
+!! @param[in] gs = size of groups (along Y direction)
+!! @param[in] p_pos_adim = adimensionned particles position
+!! @param[in] p_V = particles velocity (needed for tag and type)
+!! @param[in] i,k = indice of of the current line (y-coordinate and z-coordinate)
+!! @param[in] dt = time step (needed for tag and type)
+!! @param[in,out] scal = scalar field to advect
+subroutine Yremesh_O2(direction, ind_group, gs, p_pos_adim, P_V,i,k,scal, dt)
+
+ use advec_common ! Some procedures common to advection along all directions
+ use advec_common_line ! Some procedures common to advection along all directions
+ use advec_remeshing_line ! Remeshing formula
+ use advec_variables ! contains info about solver parameters and others.
+ use cart_topology ! Description of mesh and of mpi topology
+
+ ! Input/Output
+ integer, intent(in) :: direction
+ integer, dimension(2), intent(in) :: ind_group
+ integer, dimension(2), intent(in) :: gs
+ integer, intent(in) :: i, k
+ real(WP), dimension(:,:,:), intent(in) :: p_pos_adim ! adimensionned particles position
+ real(WP), dimension(:,:,:), intent(in) :: p_V ! particles velocity
+ real(WP), dimension(:,:,:), intent(inout) :: scal
+ real(WP), intent(in) :: dt
+ ! Other local variables
+ ! To type and tag particles
+ logical, dimension(bl_nb(direction)+1,gs(1),gs(2)) :: bl_type ! is the particle block a center block or a left one ?
+ logical, dimension(bl_nb(direction),gs(1),gs(2)) :: bl_tag ! indice of tagged particles
+ ! To compute recquired communications
+ integer, dimension(gs(1), gs(2)) :: send_group_min ! distance between me and processus wich send me information
+ integer, dimension(gs(1), gs(2)) :: send_group_max ! distance between me and processus wich send me information
+ ! Variable used to remesh particles in a buffer
+ real(WP),dimension(:),allocatable :: send_buffer ! buffer use to remesh the scalar before to send it to the right subdomain
+ integer, dimension(2,gs(1),gs(2)) :: rece_proc ! minimal and maximal gap between my Y-coordinate and the one from which
+ ! I will receive data
+ integer, dimension(gs(1),gs(2)) :: proc_min ! smaller gap between me and the processes to where I send data
+ integer, dimension(gs(1),gs(2)) :: proc_max ! smaller gap between me and the processes to where I send data
+
+ integer :: i1, i2 ! indice of a line into the group
+
+ ! -- Pre-Remeshing: Determine blocks type and tag particles --
+ call AC_type_and_block_group(dt, direction, gs, ind_group, p_V, bl_type, bl_tag)
+
+ ! -- Compute ranges --
+ where (bl_type(1,:,:))
+ ! First particle is a centered one
+ send_group_min = nint(p_pos_adim(1,:,:))-1
+ elsewhere
+ ! First particle is a left one
+ send_group_min = floor(p_pos_adim(1,:,:))-1
+ end where
+ where (bl_type(mesh_sc%N_proc(direction)/bl_size +1,:,:))
+ ! Last particle is a centered one
+ send_group_max = nint(p_pos_adim(mesh_sc%N_proc(direction),:,:))+1
+ elsewhere
+ ! Last particle is a left one
+ send_group_max = floor(p_pos_adim(mesh_sc%N_proc(direction),:,:))+1
+ end where
+
+ ! -- Determine the communication needed : who will communicate whit who ? (ie compute sender and recer) --
+ call AC_obtain_senders_group(direction, gs, ind_group, send_group_min, &
+ & send_group_max, proc_min, proc_max, rece_proc)
+
+ do i2 = 1, gs(2)
+ do i1 = 1, gs(1)
+ send_j_min = send_group_min(i1,i2)
+ send_j_max = send_group_max(i1,i2)
+
+ ! -- Allocate buffer for remeshing of local particles --
+ allocate(send_buffer(send_j_min:send_j_max))
+ send_buffer = 0.0;
+
+ ! -- Remesh the particles in the buffer --
+ call AC_remesh_lambda2corrected_basic(direction, p_pos_adim(:,i1,i2), scal(i+i1-1,:,k+i2-1), &
+ & bl_type(:,i1,i2), bl_tag(:,i1,i2), send_j_min, send_j_max, send_buffer)
+
+ ! -- Send the buffer to the matching processus and update the scalar field --
+ scal(i+i1-1,:,k+i2-1) = 0
+ call AC_bufferToScalar_line(direction, ind_group , send_j_min, send_j_max, proc_min(i1,i2), proc_max(i1,i2), &
+ & rece_proc(:,i1,i2), send_buffer, scal(i+i1-1,:,k+i2-1))
+
+ ! Deallocate all field
+ deallocate(send_buffer)
+
+ end do
+ end do
+
+end subroutine Yremesh_O2
+
+
+!> remeshing along Y with an order 4 method, corrected to allow large CFL number - untagged particles
+!! @param[in] direction = current direction
+!! @param[in] ind_group = coordinate of the current group of lines
+!! @param[in] gs = size of groups (along Y direction)
+!! @param[in] p_pos_adim = adimensionned particles position
+!! @param[in] p_V = particles velocity (needed for tag and type)
+!! @param[in] bl_tag = contains information about block (is it tagged ?)
+!! @param[in] i,k = indice of of the current line (x-coordinate and z-coordinate)
+!! @param[in] bl_type = equal 0 (resp 1) if the block is left (resp centered)
+!! @param[in,out] scal = scalar field to advect
+!! @param[in] dt = time step (needed for tag and type)
+subroutine Yremesh_O4(direction, ind_group, gs, p_pos_adim, p_V, i,k,scal, dt)
+
+ use advec_common ! Some procedures common to advection along all directions
+ use advec_common_line ! Some procedures common to advection along all directions
+ use advec_remeshing_line ! Remeshing formula
+ use advec_variables ! contains info about solver parameters and others.
+ use cart_topology ! Description of mesh and of mpi topology
+
+ ! input/output
+ integer, intent(in) :: direction
+ integer, dimension(2), intent(in) :: ind_group
+ integer, dimension(2), intent(in) :: gs
+ integer, intent(in) :: i, k
+ real(WP), dimension(:,:,:), intent(in) :: p_pos_adim ! adimensionned particles position
+ real(WP), dimension(:,:,:), intent(in) :: p_V ! particles velocity
+ real(WP), dimension(:,:,:), intent(inout) :: scal
+ real(WP), intent(in) :: dt
+ ! Other local variables
+ ! To compute recquired communications
+ integer, dimension(gs(1), gs(2)) :: send_group_min ! distance between me and processus wich send me information
+ integer, dimension(gs(1), gs(2)) :: send_group_max ! distance between me and processus wich send me information
+ ! To type and tag particles
+ logical, dimension(bl_nb(direction)+1,gs(1),gs(2)) :: bl_type ! is the particle block a center block or a left one ?
+ logical, dimension(bl_nb(direction),gs(1),gs(2)) :: bl_tag ! indice of tagged particles
+ ! Variables used to remesh particles ...
+ ! ... and to communicate between subdomains. A variable prefixed by "send_"(resp "rece")
+ ! designes something I send (resp. I receive).
+ real(WP),dimension(:),allocatable :: send_buffer ! buffer use to remesh the scalar before to send it to the right subdomain
+ integer, dimension(2,gs(1),gs(2)) :: rece_proc ! minimal and maximal gap between my Y-coordinate and the one from which
+ ! I will receive data
+ integer, dimension(gs(1),gs(2)) :: proc_min ! smaller gap between me and the processes to where I send data
+ integer, dimension(gs(1),gs(2)) :: proc_max ! smaller gap between me and the processes to where I send data
+ integer :: i1, i2 ! indice of a line into the group
+
+ ! -- Pre-Remeshing: Determine blocks type and tag particles --
+ call AC_type_and_block_group(dt, direction, gs, ind_group, p_V, bl_type, bl_tag)
+
+ ! -- Compute ranges --
+ where (bl_type(1,:,:))
+ ! First particle is a centered one
+ send_group_min = nint(p_pos_adim(1,:,:))-2
+ elsewhere
+ ! First particle is a left one
+ send_group_min = floor(p_pos_adim(1,:,:))-2
+ end where
+ where (bl_type(mesh_sc%N_proc(direction)/bl_size +1,:,:))
+ ! Last particle is a centered one
+ send_group_max = nint(p_pos_adim(mesh_sc%N_proc(direction),:,:))+2
+ elsewhere
+ ! Last particle is a left one
+ send_group_max = floor(p_pos_adim(mesh_sc%N_proc(direction),:,:))+2
+ end where
+
+ ! -- Determine the communication needed : who will communicate whit who ? (ie compute sender and recer) --
+ call AC_obtain_senders_group(direction, gs, ind_group, send_group_min, &
+ & send_group_max, proc_min, proc_max, rece_proc)
+
+ do i2 = 1, gs(2)
+ do i1 = 1, gs(1)
+ send_j_min = send_group_min(i1,i2)
+ send_j_max = send_group_max(i1,i2)
+
+ ! -- Allocate buffer for remeshing of local particles --
+ allocate(send_buffer(send_j_min:send_j_max))
+ send_buffer = 0.0;
+
+ ! -- Remesh the particles in the buffer --
+ call AC_remesh_lambda4corrected_basic(direction, p_pos_adim(:,i1,i2), scal(i+i1-1,:,k+i2-1), &
+ & bl_type(:,i1,i2), bl_tag(:,i1,i2), send_j_min, send_j_max, send_buffer)
+
+ ! -- Send the buffer to the matching processus and update the scalar field --
+ scal(i+i1-1,:,k+i2-1) = 0
+ call AC_bufferToScalar_line(direction, ind_group, send_j_min, send_j_max, proc_min(i1,i2), proc_max(i1,i2), &
+ & rece_proc(:,i1,i2), send_buffer, scal(i+i1-1,:,k+i2-1))
+
+ ! Deallocate all field
+ deallocate(send_buffer)
+
+ end do
+ end do
+
+end subroutine Yremesh_O4
+
+
+!> remeshing along Y with M'6 formula - No tag neither correction for large time steps.
+!! @param[in] direction = current direction
+!! @param[in] ind_group = coordinate of the current group of lines
+!! @param[in] gs = size of groups (along Y direction)
+!! @param[in] p_pos_adim = adimensionned particles position
+!! @param[in] p_V = particles velocity (only to have the same profile
+!! then other remeshing procedures)
+!! @param[in] i,k = indice of of the current line (x-coordinate and z-coordinate)
+!! @param[in,out] scal = scalar field to advect
+!! @param[in] dt = time step (only to have the same profile
+!! then other remeshing procedures)
+subroutine Yremesh_Mprime6(direction, ind_group, gs, p_pos_adim, p_V, i,k,scal, dt)
+
+ use advec_common_line ! Some procedures common to advection along all directions
+ use advec_remeshing_line ! Remeshing formula
+ use advec_variables ! contains info about solver parameters and others.
+ use cart_topology ! Description of mesh and of mpi topology
+
+ ! input/output
+ integer, intent(in) :: direction
+ integer, dimension(2), intent(in) :: ind_group
+ integer, dimension(2), intent(in) :: gs
+ integer, intent(in) :: i, k
+ real(WP), dimension(:,:,:), intent(in) :: p_pos_adim ! adimensionned particles position
+ real(WP), dimension(:,:,:), intent(in) :: p_V ! particles velocity
+ real(WP), dimension(:,:,:), intent(inout) :: scal
+ real(WP), intent(in) :: dt
+ ! Other local variables
+ ! To compute recquired communications
+ integer, dimension(gs(1), gs(2)) :: send_group_min ! distance between me and processus wich send me information
+ integer, dimension(gs(1), gs(2)) :: send_group_max ! distance between me and processus wich send me information
+ ! Variables used to remesh particles ...
+ ! ... and to communicate between subdomains. A variable prefixed by "send_"(resp "rece")
+ ! designes something I send (resp. I receive).
+ real(WP),dimension(:),allocatable :: send_buffer ! buffer use to remesh the scalar before to send it to the right subdomain
+ integer, dimension(2,gs(1),gs(2)) :: rece_proc ! minimal and maximal gap between my Y-coordinate and the one from which
+ ! I will receive data
+ integer, dimension(gs(1),gs(2)) :: proc_min ! smaller gap between me and the processes to where I send data
+ integer, dimension(gs(1),gs(2)) :: proc_max ! smaller gap between me and the processes to where I send data
+ integer :: i1, i2 ! indice of a line into the group
+ integer :: ind_p ! indice of the current particle
+
+ ! -- Compute the remeshing domain --
+ send_group_min = floor(p_pos_adim(1,:,:)-2)
+ send_group_max = floor(p_pos_adim(mesh_sc%N_proc(direction),:,:)+3)
+
+ ! -- Determine the communication needed : who will communicate whit who ? (ie compute sender and recer) --
+ call AC_obtain_senders_com(direction, gs, ind_group, send_group_min, &
+ & send_group_max, proc_min, proc_max, rece_proc)
+
+ do i2 = 1, gs(2)
+ do i1 = 1, gs(1)
+ send_j_min = send_group_min(i1,i2)
+ send_j_max = send_group_max(i1,i2)
+
+ ! -- Allocate and initialize the buffer --
+ allocate(send_buffer(send_j_min:send_j_max))
+ send_buffer = 0.0;
+
+ ! -- Remesh the particles in the buffer --
+ do ind_p = 1, mesh_sc%N_proc(direction), 1
+ call AC_remesh_Mprime6(p_pos_adim(ind_p,i1,i2),scal(i+i1-1,ind_p,k+i2-1), send_buffer)
+ end do
+
+ ! -- Send the buffer to the matching processus and update the scalar field --
+ scal(i+i1-1,:,k+i2-1) = 0
+ call AC_bufferToScalar_line(direction, ind_group, send_j_min, send_j_max, proc_min(i1,i2), proc_max(i1,i2), &
+ & rece_proc(:,i1,i2), send_buffer, scal(i+i1-1,:,k+i2-1))
+
+ ! Deallocate all field
+ deallocate(send_buffer)
+
+ end do
+ end do
+
+end subroutine Yremesh_Mprime6
+
+
+! ============================================================================
+! ==================== Remeshing along Z subroutines ====================
+! ============================================================================
+
+!> remeshing along Z with an order 2 method, corrected to allow large CFL number - group version
+!! @param[in] direction = current direction
+!! @param[in] ind_group = coordinate of the current group of lines
+!! @param[in] gs = size of groups (along Z direction)
+!! @param[in] p_pos_adim = adimensionned particles position
+!! @param[in] p_V = particles velocity (needed for tag and type)
+!! @param[in] i,j = indice of of the current line (x-coordinate and z-coordinate)
+!! @param[in,out] scal = scalar field to advect
+!! @param[in] dt = time step (needed for tag and type)
+subroutine Zremesh_O2(direction, ind_group, gs, p_pos_adim, p_V,i,j,scal, dt)
+
+ use advec_common ! Some procedures common to advection along all directions
+ use advec_common_line ! Some procedures common to advection along all directions
+ use advec_remeshing_line ! Remeshing formula
+ use advec_variables ! contains info about solver parameters and others.
+ use cart_topology ! Description of mesh and of mpi topology
+
+ ! Input/Output
+ integer, intent(in) :: direction
+ integer, dimension(2), intent(in) :: ind_group
+ integer, dimension(2), intent(in) :: gs
+ integer, intent(in) :: i, j
+ real(WP), dimension(:,:,:), intent(in) :: p_pos_adim ! adimensionned particles position
+ real(WP), dimension(:,:,:), intent(in) :: p_V ! particles velocity
+ real(WP), dimension(:,:,:), intent(inout) :: scal
+ real(WP), intent(in) :: dt
+ ! Other local variables
+ ! To compute recquired communications
+ integer, dimension(gs(1), gs(2)) :: send_group_min ! distance between me and processus wich send me information
+ integer, dimension(gs(1), gs(2)) :: send_group_max ! distance between me and processus wich send me information
+ ! To type and tag particles
+ logical, dimension(bl_nb(direction)+1,gs(1),gs(2)) :: bl_type ! is the particle block a center block or a left one ?
+ logical, dimension(bl_nb(direction),gs(1),gs(2)) :: bl_tag ! indice of tagged particles
+ ! Variable used to remesh particles in a buffer
+ real(WP),dimension(:),allocatable :: send_buffer ! buffer use to remesh the scalar before to send it to the right subdomain
+ integer, dimension(2,gs(1),gs(2)) :: rece_proc ! minimal and maximal gap between my Y-coordinate and the one from which
+ ! I will receive data
+ integer, dimension(gs(1),gs(2)) :: proc_min ! smaller gap between me and the processes to where I send data
+ integer, dimension(gs(1),gs(2)) :: proc_max ! smaller gap between me and the processes to where I send data
+
+ integer :: i1, i2 ! indice of a line into the group
+
+ ! -- Pre-Remeshing: Determine blocks type and tag particles --
+ call AC_type_and_block_group(dt, direction, gs, ind_group, p_V, bl_type, bl_tag)
+
+ ! -- Compute ranges --
+ where (bl_type(1,:,:))
+ ! First particle is a centered one
+ send_group_min = nint(p_pos_adim(1,:,:))-1
+ elsewhere
+ ! First particle is a left one
+ send_group_min = floor(p_pos_adim(1,:,:))-1
+ end where
+ where (bl_type(mesh_sc%N_proc(direction)/bl_size +1,:,:))
+ ! Last particle is a centered one
+ send_group_max = nint(p_pos_adim(mesh_sc%N_proc(direction),:,:))+1
+ elsewhere
+ ! Last particle is a left one
+ send_group_max = floor(p_pos_adim(mesh_sc%N_proc(direction),:,:))+1
+ end where
+
+ ! -- Determine the communication needed : who will communicate whit who ? (ie compute sender and recer) --
+ call AC_obtain_senders_group(direction, gs, ind_group, send_group_min, &
+ & send_group_max, proc_min, proc_max, rece_proc)
+
+ do i2 = 1, gs(2)
+ do i1 = 1, gs(1)
+ send_j_min = send_group_min(i1,i2)
+ send_j_max = send_group_max(i1,i2)
+
+ ! -- Allocate buffer for remeshing of local particles --
+ allocate(send_buffer(send_j_min:send_j_max))
+ send_buffer = 0.0;
+
+ ! -- Remesh the particles in the buffer --
+ call AC_remesh_lambda2corrected_basic(direction, p_pos_adim(:,i1,i2), scal(i+i1-1,j+i2-1,:), &
+ & bl_type(:,i1,i2), bl_tag(:,i1,i2), send_j_min, send_j_max, send_buffer)
+
+ ! -- Send the buffer to the matching processus and update the scalar field --
+ scal(i+i1-1,j+i2-1,:) = 0
+ call AC_bufferToScalar_line(direction, ind_group , send_j_min, send_j_max, proc_min(i1,i2), proc_max(i1,i2), &
+ & rece_proc(:,i1,i2), send_buffer, scal(i+i1-1,j+i2-1,:))
+
+ ! Deallocate all field
+ deallocate(send_buffer)
+
+ end do
+ end do
+
+end subroutine Zremesh_O2
+
+
+!> remeshing along Z with an order 4 method, corrected to allow large CFL number - untagged particles
+!! @param[in] direction = current direction
+!! @param[in] ind_group = coordinate of the current group of lines
+!! @param[in] gs = size of groups (along Z direction)
+!! @param[in] p_pos_adim = adimensionned particles position
+!! @param[in] p_V = particles velocity (needed for tag and type)
+!! @param[in] i,j = indice of of the current line (x-coordinate and z-coordinate)
+!! @param[in,out] scal = scalar field to advect
+!! @param[in] dt = time step (needed for tag and type)
+subroutine Zremesh_O4(direction, ind_group, gs, p_pos_adim, p_V,i,j,scal, dt)
+
+ use advec_common ! Some procedures common to advection along all directions
+ use advec_common_line ! Some procedures common to advection along all directions
+ use advec_remeshing_line ! Remeshing formula
+ use advec_variables ! contains info about solver parameters and others.
+ use cart_topology ! Description of mesh and of mpi topology
+
+ ! Input/Output
+ integer, intent(in) :: direction
+ integer, dimension(2), intent(in) :: ind_group
+ integer, dimension(2), intent(in) :: gs
+ integer, intent(in) :: i, j
+ real(WP), dimension(:,:,:), intent(in) :: p_pos_adim ! adimensionned particles position
+ real(WP), dimension(:,:,:), intent(in) :: p_V ! particles velocity
+ real(WP), dimension(:,:,:), intent(inout) :: scal
+ real(WP), intent(in) :: dt
+ ! Other local variables
+ ! To type and tag particles
+ logical, dimension(bl_nb(direction)+1,gs(1),gs(2)) :: bl_type ! is the particle block a center block or a left one ?
+ logical, dimension(bl_nb(direction),gs(1),gs(2)) :: bl_tag ! indice of tagged particles
+ ! To compute recquired communications
+ integer, dimension(gs(1), gs(2)) :: send_group_min ! distance between me and processus wich send me information
+ integer, dimension(gs(1), gs(2)) :: send_group_max ! distance between me and processus wich send me information
+ ! Variables used to remesh particles ...
+ ! ... and to communicate between subdomains. A variable prefixed by "send_"(resp "rece")
+ ! designes something I send (resp. I receive).
+ real(WP),dimension(:),allocatable :: send_buffer ! buffer use to remesh the scalar before to send it to the right subdomain
+ integer, dimension(2,gs(1),gs(2)) :: rece_proc ! minimal and maximal gap between my Y-coordinate and the one from which
+ ! I will receive data
+ integer, dimension(gs(1),gs(2)) :: proc_min ! smaller gap between me and the processes to where I send data
+ integer, dimension(gs(1),gs(2)) :: proc_max ! smaller gap between me and the processes to where I send data
+ integer :: i1, i2 ! indice of a line into the group
+
+ ! -- Pre-Remeshing: Determine blocks type and tag particles --
+ call AC_type_and_block_group(dt, direction, gs, ind_group, p_V, bl_type, bl_tag)
+
+ ! -- Compute ranges --
+ where (bl_type(1,:,:))
+ ! First particle is a centered one
+ send_group_min = nint(p_pos_adim(1,:,:))-2
+ elsewhere
+ ! First particle is a left one
+ send_group_min = floor(p_pos_adim(1,:,:))-2
+ end where
+ where (bl_type(mesh_sc%N_proc(direction)/bl_size +1,:,:))
+ ! Last particle is a centered one
+ send_group_max = nint(p_pos_adim(mesh_sc%N_proc(direction),:,:))+2
+ elsewhere
+ ! Last particle is a left one
+ send_group_max = floor(p_pos_adim(mesh_sc%N_proc(direction),:,:))+2
+ end where
+
+ ! -- Determine the communication needed : who will communicate whit who ? (ie compute sender and recer) --
+ call AC_obtain_senders_group(direction, gs, ind_group, send_group_min, &
+ & send_group_max, proc_min, proc_max, rece_proc)
+
+ do i2 = 1, gs(2)
+ do i1 = 1, gs(1)
+ send_j_min = send_group_min(i1,i2)
+ send_j_max = send_group_max(i1,i2)
+
+ ! -- Allocate buffer for remeshing of local particles --
+ allocate(send_buffer(send_j_min:send_j_max))
+ send_buffer = 0.0;
+
+ ! -- Remesh the particles in the buffer --
+ call AC_remesh_lambda4corrected_basic(direction, p_pos_adim(:,i1,i2), scal(i+i1-1,j+i2-1,:), &
+ & bl_type(:,i1,i2), bl_tag(:,i1,i2), send_j_min, send_j_max, send_buffer)
+
+ ! -- Send the buffer to the matching processus and update the scalar field --
+ scal(i+i1-1,j+i2-1,:) = 0
+ call AC_bufferToScalar_line(direction, ind_group , send_j_min, send_j_max, proc_min(i1,i2), proc_max(i1,i2), &
+ & rece_proc(:,i1,i2), send_buffer, scal(i+i1-1,j+i2-1,:))
+
+ ! Deallocate all field
+ deallocate(send_buffer)
+
+ end do
+ end do
+
+end subroutine Zremesh_O4
+
+
+!> remeshing along Z with M'6 formula - No tag neither correction for large time steps.
+!! @param[in] direction = current direction
+!! @param[in] ind_group = coordinate of the current group of lines
+!! @param[in] gs = size of groups (along Z direction)
+!! @param[in] p_pos_adim = adimensionned particles position
+!! @param[in] p_V = particles velocity (only to have the same profile
+!! then other remeshing procedures)
+!! @param[in] i,j = indice of of the current line (x-coordinate and y-coordinate)
+!! @param[in,out] scal = scalar field to advect
+!! @param[in] dt = time step (only to have the same profile
+!! then other remeshing procedures)
+subroutine Zremesh_Mprime6(direction, ind_group, gs, p_pos_adim, p_V, i,j,scal, dt)
+
+ use advec_common_line ! Some procedures common to advection along all directions
+ use advec_remeshing_line ! Remeshing formula
+ use advec_variables ! contains info about solver parameters and others.
+ use cart_topology ! Description of mesh and of mpi topology
+
+ ! input/output
+ integer, intent(in) :: direction
+ integer, dimension(2), intent(in) :: ind_group
+ integer, dimension(2), intent(in) :: gs
+ integer, intent(in) :: i, j
+ real(WP), dimension(:,:,:), intent(in) :: p_pos_adim ! adimensionned particles position
+ real(WP), dimension(:,:,:), intent(in) :: p_V ! particles velocity
+ real(WP), dimension(:,:,:), intent(inout) :: scal
+ real(WP), intent(in) :: dt
+ ! Other local variables
+ ! To compute recquired communications
+ integer, dimension(gs(1), gs(2)) :: send_group_min ! distance between me and processus wich send me information
+ integer, dimension(gs(1), gs(2)) :: send_group_max ! distance between me and processus wich send me information
+ ! Variables used to remesh particles ...
+ ! ... and to communicate between subdomains. A variable prefixed by "send_"(resp "rece")
+ ! designes something I send (resp. I receive).
+ real(WP),dimension(:),allocatable :: send_buffer ! buffer use to remesh the scalar before to send it to the right subdomain
+ integer, dimension(2,gs(1),gs(2)) :: rece_proc ! minimal and maximal gap between my Y-coordinate and the one from which
+ ! I will receive data
+ integer, dimension(gs(1),gs(2)) :: proc_min ! smaller gap between me and the processes to where I send data
+ integer, dimension(gs(1),gs(2)) :: proc_max ! smaller gap between me and the processes to where I send data
+ integer :: i1, i2 ! indice of a line into the group
+ integer :: ind_p ! indice of the current particle
+
+ ! -- Compute the remeshing domain --
+ send_group_min = floor(p_pos_adim(1,:,:)-2)
+ send_group_max = floor(p_pos_adim(mesh_sc%N_proc(direction),:,:)+3)
+
+ ! -- Determine the communication needed : who will communicate whit who ? (ie compute sender and recer) --
+ call AC_obtain_senders_com(direction, gs, ind_group, send_group_min, &
+ & send_group_max, proc_min, proc_max, rece_proc)
+
+ do i2 = 1, gs(2)
+ do i1 = 1, gs(1)
+ send_j_min = send_group_min(i1,i2)
+ send_j_max = send_group_max(i1,i2)
+
+ ! -- Allocate and initialize the buffer --
+ allocate(send_buffer(send_j_min:send_j_max))
+ send_buffer = 0.0;
+
+ ! -- Remesh the particles in the buffer --
+ do ind_p = 1, mesh_sc%N_proc(direction), 1
+ call AC_remesh_Mprime6(p_pos_adim(ind_p,i1,i2),scal(i+i1-1,j+i2-1, ind_p), send_buffer)
+ end do
+
+ ! -- Send the buffer to the matching processus and update the scalar field --
+ scal(i+i1-1,j+i2-1,:) = 0
+ call AC_bufferToScalar_line(direction, ind_group, send_j_min, send_j_max, proc_min(i1,i2), proc_max(i1,i2), &
+ & rece_proc(:,i1,i2), send_buffer, scal(i+i1-1,j+i2-1,:))
+
+ ! Deallocate all field
+ deallocate(send_buffer)
+
+ end do
+ end do
+
+end subroutine Zremesh_Mprime6
+
+
+! #####################################################################################
+! ##### #####
+! ##### Private procedure #####
+! ##### #####
+! #####################################################################################
+
+! =====================================================================================
+! ==================== Remeshing tool to determine comunications ====================
+! =====================================================================================
+
+!> Determine the set of processes wich will send me information during the remeshing
+!! and compute for each of these processes the range of wanted data. Use implicit
+!! computation rather than communication (only possible if particle are gather by
+!! block whith contrainst on velocity variation - as corrected lambda formula.) -
+!! work directly on a group of particles lines.
+! @param[in] send_group_min = minimal indice of mesh involved in remeshing particles (of the particles in my local subdomains)
+! @param[in] send_group_max = maximal indice of mesh involved in remeshing particles (of the particles in my local subdomains)
+!! @param[in] direction = current direction (1 = along X, 2 = along Y, 3 = along Z)
+!! @param[in] ind_group = coordinate of the current group of lines
+!! @param[out] send_min = minimal indice of mesh involved in remeshing particles (of the particles in my local subdomains)
+!! @param[out] send_max = maximal indice of mesh involved in remeshing particles (of the particles in my local subdomains)
+!! @param[out] proc_min = gap between my coordinate and the processes of minimal coordinate which will receive information from me
+!! @param[out] proc_max = gap between my coordinate and the processes of maximal coordinate which will receive information from me
+!! @param[out] rece_proc = coordinate range of processes which will send me information during the remeshing.
+!! @param[in] gp_s = size of group of line along the current direction
+!! @details
+!! Work on a group of line of size gs(1) x gs(2))
+!! Obtain the list of processts which are associated to sub-domain where my partticles
+!! will be remeshed and the list of processes wich contains particles which
+!! have to be remeshed in my sub-domain. This way, this procedure determine
+!! which processus need to communicate together in order to proceed to the
+!! remeshing (as in a parrallel context the real space is subdivised and each
+!! processus contains a part of it)
+!! In the same time, it computes for each processus with which I will
+!! communicate, the range of mesh point involved for each line of particles
+!! inside the group and it stores it by using some sparse matrix technics
+!! (see cartography defined in the algorithm documentation)
+!! This routine does not involve any computation to determine if
+!! a processus is the first or the last processes (considering its coordinate along
+!! the current directory) to send remeshing information to a given processes.
+!! It directly compute it using contraints on velocity (as in corrected lambda
+!! scheme) When possible use it rather than AC_obtain_senders_com
+subroutine AC_obtain_senders_group(direction, gp_s, ind_group, send_min, send_max, proc_min, proc_max, rece_proc)
+! XXX Work only for periodic condition. For dirichlet conditions : it is
+! possible to not receive either rece_proc(1), either rece_proc(2) or none of
+! these two => detect it (track the first and the last particles) and deal with it.
+
+ use cart_topology ! info about mesh and mpi topology
+ use advec_variables ! contains info about solver parameters and others.
+ use mpi
+
+ ! Input/output
+ integer, intent(in) :: direction
+ integer, dimension(2), intent(in) :: ind_group
+ integer(kind=4), dimension(:,:), intent(out) :: proc_min, proc_max
+ integer, dimension(:,:,:), intent(out) :: rece_proc
+ integer, dimension(2), intent(in) :: gp_s
+ integer, dimension(:,:), intent(in) :: send_min ! distance between me and processus wich send me information
+ integer, dimension(:,:), intent(in) :: send_max ! distance between me and processus wich send me information
+ ! Other local variable
+ integer(kind=4) :: proc_gap ! gap between a processus coordinate (along the current
+ ! direction) into the mpi-topology and my coordinate
+ integer :: rankP, rankN ! processus rank for shift (P= previous, N = next)
+ integer, dimension(2) :: tag_table ! mpi message tag (for communicate rece_proc(1) and rece_proc(2))
+ integer :: proc_max_abs ! maximum of proc_max array
+ integer :: proc_min_abs ! minimum of proc_min array
+ integer, dimension(:,:), allocatable :: first, last ! Storage processus to which I will be the first (or the last) to send
+ ! remeshed particles
+ integer, dimension(2) :: first_condition ! allowed range of value of proc_min and proc_max for being the first
+ integer, dimension(2) :: last_condition ! allowed range of value of proc_min and proc_max for being the last
+ integer, dimension(:,:),allocatable :: send_request ! mpi status of nonblocking send
+ integer :: ierr ! mpi error code
+ integer, dimension(MPI_STATUS_SIZE) :: statut ! mpi status
+ integer :: ind1, ind2 ! indice of the current line inside the group
+ integer :: min_size ! begin indice in first and last to stock indice along first dimension of the group line
+ integer :: max_size ! maximum size of first/last along the first direction
+ integer :: indice ! internal indice
+ integer, dimension(1 + gp_s(2)*(2+gp_s(1))) :: rece_buffer ! buffer for reception of rece_max
+
+ rece_proc = 3*mesh_sc%N(direction)
+ max_size = size(rece_buffer) + 1
+
+ proc_min = floor(real(send_min-1)/mesh_sc%N_proc(direction))
+ proc_max = floor(real(send_max-1)/mesh_sc%N_proc(direction))
+ proc_min_abs = minval(proc_min)
+ proc_max_abs = maxval(proc_max)
+
+ allocate(send_request(proc_min_abs:proc_max_abs,3))
+ send_request(:,3) = 0
+
+ ! -- Determine if I am the first or the last to send information to a given
+ ! processus and sort line by target processes for which I am the first and
+ ! for which I am the last. --
+ tag_table = compute_tag(ind_group, tag_obtsend_NP, direction)
+ min_size = 2 + gp_s(2)
+ allocate(first(max_size,proc_min_abs:proc_max_abs))
+ first = 0
+ first(1,:) = min_size
+ allocate(last(max_size,proc_min_abs:proc_max_abs))
+ last = 0
+ last(1,:) = min_size
+ do proc_gap = proc_min_abs, proc_max_abs
+ first(2,proc_gap) = -proc_gap
+ last(2,proc_gap) = -proc_gap
+ first_condition(2) = proc_gap*mesh_sc%N_proc(direction)+1
+ first_condition(1) = 1-2*bl_bound_size + first_condition(2)
+ last_condition(2) = (proc_gap+1)*mesh_sc%N_proc(direction)
+ last_condition(1) = -1+2*bl_bound_size + last_condition(2)
+ do ind2 = 1, gp_s(2)
+ first(2+ind2,proc_gap) = 0
+ last(2+ind2,proc_gap) = 0
+ do ind1 = 1, gp_s(1)
+ ! Compute if I am the first.
+ if ((send_min(ind1,ind2)< first_condition(1)).AND. &
+ & (send_max(ind1,ind2)>= first_condition(2))) then
+ first(2+ind2,proc_gap) = first(2+ind2,proc_gap)+1
+ first(1,proc_gap) = first(1,proc_gap) + 1
+ first(first(1,proc_gap),proc_gap) = ind1
+ end if
+ ! Compute if I am the last.
+ if ((send_max(ind1,ind2) > last_condition(1)) &
+ & .AND.(send_min(ind1,ind2)<= last_condition(2))) then
+ last(2+ind2,proc_gap) = last(2+ind2,proc_gap)+1
+ last(1,proc_gap) = last(1,proc_gap) + 1
+ last(last(1,proc_gap),proc_gap) = ind1
+ end if
+ end do
+ end do
+ end do
+
+#ifdef PART_DEBUG
+ do proc_gap = proc_min_abs, proc_max_abs
+ if (first(1,proc_gap)>max_size) then
+ print*, 'too big array on proc = ', cart_rank, ' - proc_gap = ', proc_gap
+ print*, 'it occurs on AC_obtain_senders_group - array concerned : "first"'
+ print*, 'first = ', first(1,proc_gap)
+ end if
+ if (last(1,proc_gap)>max_size) then
+ print*, 'too big array on proc = ', cart_rank, ' - proc_gap = ', proc_gap
+ print*, 'it occurs on AC_obtain_senders_group - array concerned : "last"'
+ print*, 'last = ', last(1,proc_gap)
+ end if
+ end do
+#endif
+
+ ! -- Send information if I am the first or the last --
+ do proc_gap = proc_min_abs, proc_max_abs
+ ! I am the first ?
+ if (first(1,proc_gap)>min_size) then
+ ! Compute the rank of the target processus
+ call mpi_cart_shift(D_comm(direction), 0, proc_gap, rankP, rankN, ierr)
+ if(rankN /= D_rank(direction)) then
+ call mpi_Isend(first(2,proc_gap), first(1,proc_gap)-1, MPI_INTEGER, rankN, tag_table(1), D_comm(direction), &
+ & send_request(proc_gap,1), ierr)
+ send_request(proc_gap,3) = 1
+ else
+ indice = min_size
+ do ind2 = 1, gp_s(2)
+ do ind1 = 1, first(2+ind2,proc_gap)
+ indice = indice+1
+ rece_proc(1,first(indice,proc_gap),ind2) = -proc_gap
+ end do
+ end do
+ end if
+ end if
+ ! I am the last ?
+ if (last(1,proc_gap)>min_size) then
+ ! Compute the rank of the target processus
+ call mpi_cart_shift(D_comm(direction), 0, proc_gap, rankP, rankN, ierr)
+ if(rankN /= D_rank(direction)) then
+ call mpi_Isend(last(2,proc_gap), last(1,proc_gap)-1, MPI_INTEGER, rankN, tag_table(2), D_comm(direction), &
+ & send_request(proc_gap,2), ierr)
+ send_request(proc_gap,3) = send_request(proc_gap, 3) + 2
+ else
+ indice = min_size
+ do ind2 = 1, gp_s(2)
+ do ind1 = 1, last(2+ind2,proc_gap)
+ indice = indice+1
+ rece_proc(2,last(indice,proc_gap),ind2) = -proc_gap
+ end do
+ end do
+ end if
+ end if
+ end do
+
+ ! -- Receive it --
+ ! size_max = size(rece_buffer) ! 2 + 2*gp_s(1)*gp_s(2)
+ max_size = max_size-1
+ do while(any(rece_proc(1,:,:) == 3*mesh_sc%N(direction)))
+ call mpi_recv(rece_buffer(1), max_size, MPI_INTEGER, MPI_ANY_SOURCE, tag_table(1), D_comm(direction), statut, ierr)
+ indice = min_size-1
+ do ind2 = 1, gp_s(2)
+ do ind1 = 1, rece_buffer(1+ind2)
+ indice = indice+1
+ rece_proc(1,rece_buffer(indice),ind2) = rece_buffer(1)
+ end do
+ end do
+ end do
+ do while(any(rece_proc(2,:,:) == 3*mesh_sc%N(direction)))
+ call mpi_recv(rece_buffer(1), max_size, MPI_INTEGER, MPI_ANY_SOURCE, tag_table(2), D_comm(direction), statut, ierr)
+ indice = min_size-1
+ do ind2 = 1, gp_s(2)
+ do ind1 = 1, rece_buffer(1+ind2)
+ indice = indice+1
+ rece_proc(2,rece_buffer(indice),ind2) = rece_buffer(1)
+ end do
+ end do
+ end do
+
+ ! -- Free Isend buffer --
+ do proc_gap = proc_min_abs, proc_max_abs
+ select case (send_request(proc_gap,3))
+ case (3)
+ call mpi_wait(send_request(proc_gap,1), statut, ierr)
+ call mpi_wait(send_request(proc_gap,2), statut, ierr)
+ case (2)
+ call mpi_wait(send_request(proc_gap,2), statut, ierr)
+ case (1)
+ call mpi_wait(send_request(proc_gap,1), statut, ierr)
+ end select
+ end do
+
+ deallocate(first)
+ deallocate(last)
+ deallocate(send_request)
+
+end subroutine AC_obtain_senders_group
+
+
+!> Determine the set of processes wich will send me information during the
+!! scalar remeshing by explicit (and exensive) way : communications !
+!! @param[in] direction = current direction (1 = along X, 2 = along Y, 3 = along Z)
+!! @param[in] ind_group = coordinate of the current group of lines
+!! @param[out] send_min = minimal indice of mesh involved in remeshing particles (of the particles in my local subdomains)
+!! @param[out] send_max = maximal indice of mesh involved in remeshing particles (of the particles in my local subdomains)
+!! @param[out] proc_min = gap between my coordinate and the processes of minimal coordinate which will receive information from me
+!! @param[out] proc_max = gap between my coordinate and the processes of maximal coordinate which will receive information from me
+!! @param[out] rece_proc = coordinate range of processes which will send me information during the remeshing.
+!! @param[in] gp_s = size of group of line along the current direction
+!! @param[in] com = integer used to distinguish this function from AC_obtain_senders_group.
+!! @details
+!! Obtain the list of processus which contains some particles which belong to
+!! my subdomains after their advection (and thus which will be remeshing into
+!! my subdomain). This result is return as an interval [send_min; send_max].
+!! All the processus whose coordinate (into the current direction) belong to
+!! this segment are involved into scalar remeshing into the current
+!! subdomains. Use this method when the sender are not predictable without
+!! communication, as in M'6 schemes for instance. More precisly, it
+!! correspond do scheme without bloc of particle involving velocity variation
+!! contrainsts to avoid that the distance between to particle grows (or dimishes)
+!! too much.
+subroutine AC_obtain_senders_com(direction, gp_s, ind_group, send_min, send_max, proc_min, proc_max, rece_proc)
+! XXX Work only for periodic condition. See AC_obtain_senders. Adapt it for
+! other condition must be more easy.
+
+ use cart_topology ! info about mesh and mpi topology
+ use advec_variables ! contains info about solver parameters and others.
+ use mpi
+
+ ! Input/output
+ integer, intent(in) :: direction
+ integer, dimension(2), intent(in) :: ind_group
+ integer(kind=4), dimension(:,:), intent(out) :: proc_min, proc_max
+ integer, dimension(:,:,:), intent(out) :: rece_proc
+ integer, dimension(2), intent(in) :: gp_s
+ integer, dimension(:,:), intent(in) :: send_min ! distance between me and processus wich send me information
+ integer, dimension(:,:), intent(in) :: send_max ! distance between me and processus wich send me information
+ ! Other local variable
+ integer(kind=4) :: proc_gap ! gap between a processus coordinate (along the current
+ ! direction) into the mpi-topology and my coordinate
+ integer :: rankP, rankN ! processus rank for shift (P= previous, N = next)
+ integer, dimension(2) :: tag_table ! mpi message tag (for communicate rece_proc(1) and rece_proc(2))
+ integer, dimension(gp_s(1), gp_s(2)) :: proc_max_prev ! maximum gap between previous processus and the receivers of its remeshing buffer
+ integer, dimension(gp_s(1), gp_s(2)) :: proc_min_next ! minimum gap between next processus and the receivers of its remeshing buffer
+ integer :: proc_max_abs ! maximum of proc_max array
+ integer :: proc_min_abs ! minimum of proc_min array
+ integer, dimension(:,:), allocatable :: first, last ! Storage processus to which I will be the first (or the last) to send
+ ! remeshed particles
+ integer, dimension(:,:),allocatable :: send_request ! mpi status of nonblocking send
+ integer :: ierr ! mpi error code
+ integer, dimension(MPI_STATUS_SIZE) :: statut ! mpi status
+ integer :: ind1, ind2 ! indice of the current line inside the group
+ integer :: min_size ! begin indice in first and last to stock indice along first dimension of the group line
+ integer :: max_size ! maximum size of first/last along the first direction
+ integer :: indice ! internal indice
+ integer, dimension(1 + gp_s(2)*(2+gp_s(1))) :: rece_buffer ! buffer for reception of rece_max
+
+
+ rece_proc = 3*mesh_sc%N(direction)
+ max_size = size(rece_buffer) + 1
+
+ proc_min = floor(real(send_min-1)/mesh_sc%N_proc(direction))
+ proc_max = floor(real(send_max-1)/mesh_sc%N_proc(direction))
+ proc_min_abs = minval(proc_min)
+ proc_max_abs = maxval(proc_max)
+
+ allocate(send_request(proc_min_abs:proc_max_abs,3))
+ send_request(:,3) = 0
+
+ ! -- Exchange send_block_min and send_block_max to determine if I am the first
+ ! or the last to send information to a given target processus. --
+ min_size = gp_s(1)*gp_s(2)
+ ! Compute message tag - we re-use tag_part_tag_NP id as using this procedure
+ ! suppose not using "AC_type_and_block"
+ tag_table = compute_tag(ind_group, tag_part_tag_NP, direction)
+ ! Exchange "ghost"
+ call mpi_Sendrecv(proc_min(1,1), min_size, MPI_INTEGER, neighbors(direction,1), tag_table(1), &
+ & proc_min_next(1,1), min_size, MPI_INTEGER, neighbors(direction,2), tag_table(1), &
+ & D_comm(direction), statut, ierr)
+ call mpi_Sendrecv(proc_max(1,1), min_size, MPI_INTEGER, neighbors(direction,2), tag_table(2), &
+ & proc_max_prev(1,1), min_size, MPI_INTEGER, neighbors(direction,1), tag_table(2), &
+ & D_comm(direction), statut, ierr)
+
+ ! -- Determine if I am the first or the last to send information to a given
+ ! processus and sort line by target processes for which I am the first and
+ ! for which I am the last. --
+ tag_table = compute_tag(ind_group, tag_obtsend_NP, direction)
+ min_size = 2 + gp_s(2)
+ allocate(first(max_size,proc_min_abs:proc_max_abs))
+ first = 0
+ first(1,:) = min_size
+ allocate(last(max_size,proc_min_abs:proc_max_abs))
+ last = 0
+ last(1,:) = min_size
+ do proc_gap = proc_min_abs, proc_max_abs
+ first(2,proc_gap) = -proc_gap
+ last(2,proc_gap) = -proc_gap
+ end do
+ do ind2 = 1, gp_s(2)
+ first(2+ind2,:) = 0
+ last(2+ind2,:) = 0
+ do ind1 = 1, gp_s(1)
+ ! Compute if I am the first, ie if:
+ ! a - proc_min <= proc_gap <= proc_max,
+ ! b - proc_gap > proc_max_prev -1.
+ do proc_gap = max(proc_min(ind1,ind2), proc_max_prev(ind1,ind2)), proc_max(ind1,ind2)
+ first(2+ind2,proc_gap) = first(2+ind2,proc_gap)+1
+ first(1,proc_gap) = first(1,proc_gap) + 1
+ first(first(1,proc_gap),proc_gap) = ind1
+ end do
+ ! Compute if I am the last, ie if:
+ ! a - proc_min <= proc_gap <= proc_max,
+ ! b - proc_gap < proc_min_next+1.
+ do proc_gap = proc_min(ind1,ind2), min(proc_min_next(ind1,ind2), proc_max(ind1,ind2))
+ last(2+ind2,proc_gap) = last(2+ind2,proc_gap)+1
+ last(1,proc_gap) = last(1,proc_gap) + 1
+ last(last(1,proc_gap),proc_gap) = ind1
+ end do
+ end do
+ end do
+
+#ifdef PART_DEBUG
+ do proc_gap = proc_min_abs, proc_max_abs
+ if (first(1,proc_gap)>max_size) then
+ print*, 'too big array on proc = ', cart_rank, ' - proc_gap = ', proc_gap
+ print*, 'it occurs on AC_obtain_senders_group - array concerned : "first"'
+ print*, 'first = ', first(1,proc_gap)
+ end if
+ if (last(1,proc_gap)>max_size) then
+ print*, 'too big array on proc = ', cart_rank, ' - proc_gap = ', proc_gap
+ print*, 'it occurs on AC_obtain_senders_group - array concerned : "last"'
+ print*, 'last = ', last(1,proc_gap)
+ end if
+ end do
+#endif
+
+ ! -- Send information if I am the first or the last --
+ do proc_gap = proc_min_abs, proc_max_abs
+ ! I am the first ?
+ if (first(1,proc_gap)>min_size) then
+ ! Compute the rank of the target processus
+ call mpi_cart_shift(D_comm(direction), 0, proc_gap, rankP, rankN, ierr)
+ if(rankN /= D_rank(direction)) then
+ call mpi_Isend(first(2,proc_gap), first(1,proc_gap)-1, MPI_INTEGER, rankN, tag_table(1), D_comm(direction), &
+ & send_request(proc_gap,1), ierr)
+ send_request(proc_gap,3) = 1
+ else
+ indice = min_size
+ do ind2 = 1, gp_s(2)
+ do ind1 = 1, first(2+ind2,proc_gap)
+ indice = indice+1
+ rece_proc(1,first(indice,proc_gap),ind2) = -proc_gap
+ end do
+ end do
+ end if
+ end if
+ ! I am the last ?
+ if (last(1,proc_gap)>min_size) then
+ ! Compute the rank of the target processus
+ call mpi_cart_shift(D_comm(direction), 0, proc_gap, rankP, rankN, ierr)
+ if(rankN /= D_rank(direction)) then
+ call mpi_Isend(last(2,proc_gap), last(1,proc_gap)-1, MPI_INTEGER, rankN, tag_table(2), D_comm(direction), &
+ & send_request(proc_gap,2), ierr)
+ send_request(proc_gap,3) = send_request(proc_gap, 3) + 2
+ else
+ indice = min_size
+ do ind2 = 1, gp_s(2)
+ do ind1 = 1, last(2+ind2,proc_gap)
+ indice = indice+1
+ rece_proc(2,last(indice,proc_gap),ind2) = -proc_gap
+ end do
+ end do
+ end if
+ end if
+ end do
+
+
+ ! -- Receive it --
+ ! size_max = size(rece_buffer) ! 2 + 2*gp_s(1)*gp_s(2)
+ max_size = max_size-1
+ do while(any(rece_proc(1,:,:) == 3*mesh_sc%N(direction)))
+ call mpi_recv(rece_buffer(1), max_size, MPI_INTEGER, MPI_ANY_SOURCE, tag_table(1), D_comm(direction), statut, ierr)
+ indice = min_size-1
+ do ind2 = 1, gp_s(2)
+ do ind1 = 1, rece_buffer(1+ind2)
+ indice = indice+1
+ rece_proc(1,rece_buffer(indice),ind2) = rece_buffer(1)
+ end do
+ end do
+ end do
+ do while(any(rece_proc(2,:,:) == 3*mesh_sc%N(direction)))
+ call mpi_recv(rece_buffer(1), max_size, MPI_INTEGER, MPI_ANY_SOURCE, tag_table(2), D_comm(direction), statut, ierr)
+ indice = min_size-1
+ do ind2 = 1, gp_s(2)
+ do ind1 = 1, rece_buffer(1+ind2)
+ indice = indice+1
+ rece_proc(2,rece_buffer(indice),ind2) = rece_buffer(1)
+ end do
+ end do
+ end do
+
+ ! -- Free Isend buffer --
+ do proc_gap = proc_min_abs, proc_max_abs
+ select case (send_request(proc_gap,3))
+ case (3)
+ call mpi_wait(send_request(proc_gap,1), statut, ierr)
+ call mpi_wait(send_request(proc_gap,2), statut, ierr)
+ case (2)
+ call mpi_wait(send_request(proc_gap,2), statut, ierr)
+ case (1)
+ call mpi_wait(send_request(proc_gap,1), statut, ierr)
+ end select
+ end do
+
+ deallocate(first)
+ deallocate(last)
+ deallocate(send_request)
+
+end subroutine AC_obtain_senders_com
+
+
+end module advec_remesh_line
+!> @}
--
GitLab