Skip to content

Waq/bugfix/delwaq 622 stable numbering boundaries #374

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
arjenmarkusdeltaresnl wants to merge 4 commits into
base: main
Choose a base branch
from
Open

Waq/bugfix/delwaq 622 stable numbering boundaries #374

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
1fc19c8
First step to implement partition-independent numbering
arjenmarkusdeltaresnl Nov 7, 2025
bfe0cc3
Correct use of the various indices
arjenmarkusdeltaresnl Nov 10, 2025
4b18252
DELWAQ-622: Remove debug output
arjenmarkusdeltaresnl Nov 10, 2025
e68930f
Count the number of boundary cells
arjenmarkusdeltaresnl Nov 11, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions src/tools_gpl/waq_tools/waqmerge/waqmerge.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@ program waqmerge
use waqmerge_version_module, only: getfullversionstring_waqmerge
use m_date_time_utils_external, only: fill_in_date_time
use m_file_path_utils, only: extract_file_extension
use m_merge_domains

implicit none

Expand Down
212 changes: 195 additions & 17 deletions src/utils_lgpl/io_hyd/packages/io_hyd/src/merge_domains.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,11 @@
!-------------------------------------------------------------------------------
!
!
module m_merge_domains
implicit none

contains

subroutine merge_domains(hyd, domain_hyd_coll)

! function : merge the domains, make pointers to final domain
Expand All @@ -40,15 +45,15 @@ subroutine merge_domains(hyd, domain_hyd_coll)
! declaration of the arguments

type(t_hydrodynamics) :: hyd ! description of the hydrodynamics
type(t_hydrodynamics_collection) :: domain_hyd_coll ! description of all domain hydrodynamics
type(t_hydrodynamics_collection) :: domain_hyd_coll ! description of all domain hydrodynamics

! local declarations

integer :: n_domain ! number of domains
integer :: i_domain ! index in collection
integer :: idmn ! flow like domain index (0:n_domain-1)
type(t_hydrodynamics), pointer :: d_hyd ! description of one domain hydrodynamics
type(t_hydrodynamics), pointer :: l_hyd ! description of a linked domain hydrodynamics
type(t_hydrodynamics), pointer :: d_hyd ! description of one domain hydrodynamics
type(t_hydrodynamics), pointer :: l_hyd ! description of a linked domain hydrodynamics
integer :: nosegl ! total number of segments per layer
integer :: num_boundary_conditions ! total number of boundaries
integer :: nobndl ! total number of boundaries per layer
Expand All @@ -68,9 +73,9 @@ subroutine merge_domains(hyd, domain_hyd_coll)
integer :: isect ! index of section
integer :: no_bnd ! number of boundaries in section
integer :: i_bnd ! index of boundary
type(t_openbnd_section), pointer :: openbndsect ! single section
type(t_open_boundary_line),pointer :: openbndlin ! single open boundary lin
type(t_openbnd_section) :: new_sect ! single section new
type(t_openbnd_section), pointer :: openbndsect ! single section
type(t_open_boundary_line),pointer :: openbndlin ! single open boundary lin
type(t_openbnd_section) :: new_sect ! single section new
logical :: bnd_active ! if a boundary is active
integer :: iret ! return value
integer :: ik ! node counter
Expand Down Expand Up @@ -113,7 +118,7 @@ subroutine merge_domains(hyd, domain_hyd_coll)
hyd%crs = d_hyd%crs
hyd%conv_type = d_hyd%conv_type
hyd%conv_version = d_hyd%conv_version

! copy waqgeom layer information from first domain
hyd%waqgeom%num_layers = d_hyd%waqgeom%num_layers
hyd%waqgeom%numtopsig = d_hyd%waqgeom%numtopsig
Expand Down Expand Up @@ -453,6 +458,9 @@ subroutine merge_domains(hyd, domain_hyd_coll)
end do
end do

! determine the stable numbering for the boundary cells
call sort_and_number_boundaries( domain_hyd_coll%hyd_pnts(1:domain_hyd_coll%current_size), num_boundary_conditions )

! exchanges
! gather the exchanges per layer, and in such a way that the internal exchanges come first, and overlap with the edge face table
num_exchanges_u_dir = 0
Expand Down Expand Up @@ -489,6 +497,7 @@ subroutine merge_domains(hyd, domain_hyd_coll)
do iq = 1, d_hyd%num_exchanges_u_dir
ip1 = d_hyd%ipoint(1,iq)
ip2 = d_hyd%ipoint(2,iq)

if ( ip1 .lt. 0 .and. ip2 .ge. min_seg .and. ip2 .le. max_seg) then
! ip1 is a boundary and ip2 is in this layer
if(d_hyd%idomain(ip2) .eq. idmn .and. (.not.d_hyd%ispoint_bnd(abs(ip1)))) then
Expand All @@ -497,8 +506,9 @@ subroutine merge_domains(hyd, domain_hyd_coll)
d_hyd%iglobal_link(iq) = num_exchanges_u_dir
if (abs(ip1) .le. d_hyd%nobndl) then
num_boundary_conditions = num_boundary_conditions + 1
d_hyd%iglobal_bnd(-ip1) = -num_boundary_conditions
call renum_bnd(d_hyd%openbndsect_coll,ip1,-num_boundary_conditions)
!d_hyd%iglobal_bnd(-ip1) = -num_boundary_conditions
!call renum_bnd(d_hyd%openbndsect_coll,ip1,-num_boundary_conditions)
call renum_bnd(d_hyd%openbndsect_coll,ip1,d_hyd%iglobal_bnd(-ip1))
end if
end if
else if (ip1 .ge. min_seg .and. ip1 .le. max_seg .and. ip2 .lt. 0) then !ip1 between min +1 and max
Expand All @@ -509,11 +519,13 @@ subroutine merge_domains(hyd, domain_hyd_coll)
d_hyd%iglobal_link(iq) = num_exchanges_u_dir
if (abs(ip2) .le. d_hyd%nobndl) then
num_boundary_conditions = num_boundary_conditions + 1
d_hyd%iglobal_bnd(-ip2) = -num_boundary_conditions
call renum_bnd(d_hyd%openbndsect_coll,ip2,-num_boundary_conditions)
!d_hyd%iglobal_bnd(-ip2) = -num_boundary_conditions
!call renum_bnd(d_hyd%openbndsect_coll,ip2,-num_boundary_conditions)
call renum_bnd(d_hyd%openbndsect_coll,ip2,d_hyd%iglobal_bnd(-ip2))
end if
end if
end if

end do
end do
end do
Expand All @@ -533,8 +545,9 @@ subroutine merge_domains(hyd, domain_hyd_coll)
d_hyd%iglobal_link(iq) = num_exchanges_u_dir
if (abs(ip1) .le. d_hyd%nobndl) then
num_boundary_conditions = num_boundary_conditions + 1
d_hyd%iglobal_bnd(-ip1) = -num_boundary_conditions
call renum_bnd(d_hyd%openbndsect_coll,ip1,-num_boundary_conditions)
!d_hyd%iglobal_bnd(-ip1) = -num_boundary_conditions
!call renum_bnd(d_hyd%openbndsect_coll,ip1,-num_boundary_conditions)
call renum_bnd(d_hyd%openbndsect_coll,ip1,d_hyd%iglobal_bnd(-ip1))
end if
end if
else if (ip1 .gt. 0 .and. ip2 .lt. 0) then
Expand All @@ -545,8 +558,9 @@ subroutine merge_domains(hyd, domain_hyd_coll)
d_hyd%iglobal_link(iq) = num_exchanges_u_dir
if (abs(ip2) .le. d_hyd%nobndl) then
num_boundary_conditions = num_boundary_conditions + 1
d_hyd%iglobal_bnd(-ip2) = -num_boundary_conditions
call renum_bnd(d_hyd%openbndsect_coll,ip2,-num_boundary_conditions)
!d_hyd%iglobal_bnd(-ip2) = -num_boundary_conditions
!call renum_bnd(d_hyd%openbndsect_coll,ip2,-num_boundary_conditions)
call renum_bnd(d_hyd%openbndsect_coll,ip1,d_hyd%iglobal_bnd(-ip2))
end if
end if
end if
Expand Down Expand Up @@ -603,7 +617,7 @@ subroutine merge_domains(hyd, domain_hyd_coll)
endif
enddo
enddo

! For z-layer models, make sure the lower most boundary exist, so that the minimum value of the pointer
! equals the number of boundary conditions per layer times the number of layers
if (minval(hyd%ipoint) /= -hyd%num_boundary_conditions) then
Expand Down Expand Up @@ -1004,7 +1018,7 @@ subroutine merge_domains_old(hyd, domain_hyd_coll)
end if
end do
end do

! For z-layer models, make sure the lower most boundary exist, so that the minimum value of the pointer
! equals the number of boundary conditions per layer times the number of layers
if (minval(hyd%ipoint) /= -hyd%num_boundary_conditions) then
Expand Down Expand Up @@ -1198,3 +1212,167 @@ subroutine merge_domains_old(hyd, domain_hyd_coll)
end do
return
end

subroutine sort_and_number_boundaries( hyd_pnts, num_boundaries )

use m_hydmod

implicit none

! Use the properties of the boundary sections to sort
! the individual boundary cells in such a way that a
! numbering results that is independent of the partitioning.
! Note: D-Flow FM does not provide such a numbering itself,
! unlike with the internal cells.

type(t_hydrodynamics), dimension(:), intent(inout) :: hyd_pnts
integer, intent(out) :: num_boundaries

type :: t_cell_properties
character(len=40) :: name
real(kind=kind(1.0d0)), dimension(4) :: coordinates
integer :: domain_index
integer :: section_index
integer :: bound_index
integer :: original_index
integer :: sorted_index = -1
end type t_cell_properties

type(t_cell_properties), allocatable, dimension(:) :: cell_list
type(t_cell_properties), allocatable, dimension(:) :: new_cell_list
type(t_cell_properties) :: new_cell
type(t_openbnd_section) :: section
type(t_open_boundary_line) :: bound
integer :: b, c, d, e, s

!
! Gather the information per domain into a single list
!
allocate( cell_list(0) )

e = 0
do d = 1,size(hyd_pnts)
do s = 1,hyd_pnts(d)%openbndsect_coll%current_size
section = hyd_pnts(d)%openbndsect_coll%openbndsect_pnts(s)

do b = 1,section%openbndlin_coll%current_size
bound = section%openbndlin_coll%openbndlin_pnts(b)
e = e + 1

new_cell%name = section%name
new_cell%coordinates = [bound%x1, bound%y1, bound%x2, bound%y2]
new_cell%domain_index = d
new_cell%section_index = s
new_cell%bound_index = e
new_cell%original_index = bound%ibnd

cell_list = [cell_list, new_cell]
enddo
enddo
enddo

!
! Sort the list:
! - First by section name
! - Then by coordinates (one by one)
!
do d = 1,size(hyd_pnts)
do s = 1,hyd_pnts(d)%openbndsect_coll%current_size
section = hyd_pnts(d)%openbndsect_coll%openbndsect_pnts(s)

do b = 1,section%openbndlin_coll%current_size
bound = section%openbndlin_coll%openbndlin_pnts(b)
enddo
enddo
enddo

call sort_lexicographically( cell_list )

call identify_duplicates( cell_list, new_cell_list )

!
! Fill the new numbering in in the domain list
!
num_boundaries = size(new_cell_list)

do c = 1,size(new_cell_list)
d = new_cell_list(c)%domain_index
s = new_cell_list(c)%section_index
b = abs( new_cell_list(c)%original_index )

!hyd_pnts(d)%openbndsect_coll%openbndsect_pnts(s)%openbndlin_coll%openbndlin_pnts(b)%ibnd_new = &
hyd_pnts(d)%iglobal_bnd(b) = -new_cell_list(c)%sorted_index
enddo

contains


subroutine sort_lexicographically( cell_list )

! Use a simple bubblesort to sort the boundary segments in a persistent way

type(t_cell_properties), dimension(:), intent(inout) :: cell_list

type(t_cell_properties) :: bstore
integer :: i, c
logical :: changed, swap

changed = .true.

do while (changed)
changed = .false.

do i = 2,size(cell_list)
swap = .false.
if ( cell_list(i)%name < cell_list(i-1)%name ) then
swap = .true.
elseif ( cell_list(i)%name == cell_list(i-1)%name ) then
do c = 1,4
if ( cell_list(i)%coordinates(c) < cell_list(i-1)%coordinates(c) ) then
swap = .true.
elseif ( cell_list(i)%coordinates(c) == cell_list(i-1)%coordinates(c) ) then
cycle
else
exit
endif
enddo
endif
if ( swap ) then
changed = .true.
bstore = cell_list(i)
cell_list(i) = cell_list(i-1)
cell_list(i-1) = bstore
endif
enddo
enddo
end subroutine sort_lexicographically


subroutine identify_duplicates( cell_list, new_cell_list )

! Go through the sorted list and remove the duplicates

type(t_cell_properties), dimension(:), intent(in) :: cell_list
type(t_cell_properties), dimension(:), intent(out), allocatable :: new_cell_list

integer :: i, j, k

new_cell_list = cell_list
j = 1
k = 1
new_cell_list(1)%sorted_index = k

do i = 2,size(cell_list)
if ( new_cell_list(i)%name /= new_cell_list(j)%name .or. &
any( new_cell_list(i)%coordinates /= new_cell_list(j)%coordinates ) ) then
j = i ! Register this entry
k = k + 1
endif
new_cell_list(i)%sorted_index = k
enddo

end subroutine identify_duplicates

end subroutine sort_and_number_boundaries

end module m_merge_domains