Skip to content

Commit

Permalink
[WIP] Create quarter administration (#360)
Browse files Browse the repository at this point in the history
  • Loading branch information
@martijn-siemerink
martijn-siemerink authored Jun 19, 2024
1 parent 8cfa925 commit 9fcf3a9
Show file tree
Hide file tree
Showing 10 changed files with 519 additions and 4 deletions.
2 changes: 1 addition & 1 deletion CHANGES.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ Changelog of threedigrid-builder
1.15.1 (unreleased)
-------------------

- Nothing changed yet.
- Create quarter administration.


1.15.0 (2024-05-22)
Expand Down
200 changes: 198 additions & 2 deletions libthreedigrid/cells.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,7 +78,7 @@ subroutine set_2d_computational_nodes_lines(origin, lgrmin, kmax, mmax, nmax, dx
deallocate(area_mask_padded)
endif
write(*,*) '** INFO: Number of 2D nodes is: ', nod - 1
write(*,*) '** INFO: Number of 2D lines is: ', l_u + l_v
write(*,*) '** INFO: Number of 2D lines is: ', l_u + (l_v - l_u)
write(*,*) '** INFO: Done setting 2D calculation cells.'

end subroutine set_2d_computational_nodes_lines
Expand Down Expand Up @@ -227,4 +227,200 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
end subroutine set_2d_computational_lines


end module m_cells
subroutine set_quarter_admin(nodk, nodm, nodn, line, kcu, quarter_line, quarter_neighbour, liutot, livtot, n2dobc)

use parameters, only : LINE_2D_BOUNDARY_EAST, LINE_2D_BOUNDARY_WEST, LINE_2D_BOUNDARY_SOUTH, LINE_2D_BOUNDARY_NORTH

integer, intent(in) :: nodk(:)
integer, intent(in) :: nodm(:)
integer, intent(in) :: nodn(:)
integer, intent(in) :: line(:, :)
integer, intent(in) :: kcu(:)
integer, intent(inout) :: quarter_line(:, :)
integer, intent(inout) :: quarter_neighbour(:, :)
integer, intent(in) :: liutot
integer, intent(in) :: livtot
integer, intent(in) :: n2dobc
integer :: l
integer :: quarter
integer :: nodd
integer :: nodu
integer :: kd
integer :: ku
integer :: md
integer :: mu
integer :: nd
integer :: nu
integer :: nb


!!!!!!!!!!!!!!!!
!! 3 !! 4 !!
!!!!!!!!!!!!!!!!
!! 1 !! 2 !!
!!!!!!!!!!!!!!!!
! Here we try to find all horizontal flow lines associated with a cell quarter.
! For two adjacent cells of samen size both neighbouring quadrants have same flow line.
! When adjacent cells differ in size due to refinement, each quadrant is associated with its own line for the larger cell. The smaller cell will have the same flow line for both quadrants. Hence the "if" structure.
! Same logic applies to finding neighbouring nodes for quarters.
! The cell quadrant index is illustrated above.
do l=1, liutot
nodd = line(l, 1) + 1 ! This is a 0-based python index, therefore +1
nodu = line(l, 2) + 1
kd = nodk(nodd)
ku = nodk(nodu)
md = nodm(nodd)
mu = nodm(nodu)
nd = nodn(nodd)
nu = nodn(nodu)
if (ku==kd) then
quarter_line(get_quarter_idx(nodd, 2), 1) = l - 1 ! This needs to be an index in Python again.
quarter_line(get_quarter_idx(nodd, 4), 1) = l - 1
quarter_line(get_quarter_idx(nodu, 1), 1) = l - 1
quarter_line(get_quarter_idx(nodu, 3), 1) = l - 1
quarter_neighbour(get_quarter_idx(nodd, 2), 1) = nodu - 1
quarter_neighbour(get_quarter_idx(nodd, 4), 1) = nodu - 1
quarter_neighbour(get_quarter_idx(nodu, 1), 1) = nodd - 1
quarter_neighbour(get_quarter_idx(nodu, 3), 1) = nodd - 1
elseif (ku == kd + 1) then
if (nd == 2 * nu - 1) then
quarter_line(get_quarter_idx(nodd, 2), 1) = l - 1
quarter_line(get_quarter_idx(nodd, 4), 1) = l - 1
quarter_line(get_quarter_idx(nodu, 1), 1) = l - 1
quarter_neighbour(get_quarter_idx(nodd, 2), 1) = nodu - 1
quarter_neighbour(get_quarter_idx(nodd, 4), 1) = nodu - 1
quarter_neighbour(get_quarter_idx(nodu, 1), 1) = nodd - 1
endif
if (nd == 2 * nu) then
quarter_line(get_quarter_idx(nodd, 2), 1) = l - 1
quarter_line(get_quarter_idx(nodd, 4), 1) = l - 1
quarter_line(get_quarter_idx(nodu, 3), 1) = l - 1
quarter_neighbour(get_quarter_idx(nodd, 2), 1) = nodu - 1
quarter_neighbour(get_quarter_idx(nodd, 4), 1) = nodu - 1
quarter_neighbour(get_quarter_idx(nodu, 3), 1) = nodd - 1
endif
else if (kd == ku + 1) then
if (nu == 2 * nd - 1) then
quarter_line(get_quarter_idx(nodd, 2), 1) = l - 1
quarter_line(get_quarter_idx(nodu, 1), 1) = l - 1
quarter_line(get_quarter_idx(nodu, 3), 1) = l - 1
quarter_neighbour(get_quarter_idx(nodd, 2), 1) = nodu - 1
quarter_neighbour(get_quarter_idx(nodu, 1), 1) = nodd - 1
quarter_neighbour(get_quarter_idx(nodu, 3), 1) = nodd - 1
elseif (nu == 2 * nd) then
quarter_line(get_quarter_idx(nodd, 4), 1) = l - 1
quarter_line(get_quarter_idx(nodu, 1), 1) = l - 1
quarter_line(get_quarter_idx(nodu, 3), 1) = l - 1
quarter_neighbour(get_quarter_idx(nodd, 4), 1) = nodu - 1
quarter_neighbour(get_quarter_idx(nodu, 1), 1) = nodd - 1
quarter_neighbour(get_quarter_idx(nodu, 3), 1) = nodd - 1
endif
endif
enddo

! Here we try to find all vertical flow lines associated with a cell quarter.
! For two adjacent cells of samen size both neighbouring quadrants have same flow line.
! When adjacent cells differ in size due to refinement, each quadrant is associated with its own line for the larger cell. The smaller cell will have the same flow line for both quadrants. Hence the "if" structure.
! Same logic applies to finding neighbouring nodes for quarters.
! The cell quadrant index is illustrated above.
do l=liutot + 1, liutot+livtot
nodd = line(l, 1) + 1
nodu = line(l, 2) + 1
kd = nodk(nodd)
ku = nodk(nodu)
md = nodm(nodd)
mu = nodm(nodu)
nd = nodn(nodd)
nu = nodn(nodu)
if (ku==kd) then
quarter_line(get_quarter_idx(nodd, 3), 2) = l - 1
quarter_line(get_quarter_idx(nodd, 4), 2) = l - 1
quarter_line(get_quarter_idx(nodu, 1), 2) = l - 1
quarter_line(get_quarter_idx(nodu, 2), 2) = l - 1
quarter_neighbour(get_quarter_idx(nodd, 3), 2) = nodu - 1
quarter_neighbour(get_quarter_idx(nodd, 4), 2) = nodu - 1
quarter_neighbour(get_quarter_idx(nodu, 1), 2) = nodd - 1
quarter_neighbour(get_quarter_idx(nodu, 2), 2) = nodd - 1
elseif (ku == kd + 1) then
if (md == 2 * mu - 1) then
quarter_line(get_quarter_idx(nodd, 3), 2) = l - 1
quarter_line(get_quarter_idx(nodd, 4), 2) = l - 1
quarter_line(get_quarter_idx(nodu, 1), 2) = l - 1
quarter_neighbour(get_quarter_idx(nodd, 3), 2) = nodu - 1
quarter_neighbour(get_quarter_idx(nodd, 4), 2) = nodu - 1
quarter_neighbour(get_quarter_idx(nodu, 1), 2) = nodd - 1
endif
if (md == 2 * mu) then
quarter_line(get_quarter_idx(nodd, 3), 2) = l - 1
quarter_line(get_quarter_idx(nodd, 4), 2) = l - 1
quarter_line(get_quarter_idx(nodu, 2), 2) = l - 1
quarter_neighbour(get_quarter_idx(nodd, 3), 2) = nodu - 1
quarter_neighbour(get_quarter_idx(nodd, 4), 2) = nodu - 1
quarter_neighbour(get_quarter_idx(nodu, 2), 2) = nodd - 1
endif
else if (kd == ku + 1) then
if (mu == 2 * md - 1) then
quarter_line(get_quarter_idx(nodd, 3), 2) = l - 1
quarter_line(get_quarter_idx(nodu, 1), 2) = l - 1
quarter_line(get_quarter_idx(nodu, 2), 2) = l - 1
quarter_neighbour(get_quarter_idx(nodd, 3), 2) = nodu - 1
quarter_neighbour(get_quarter_idx(nodu, 1), 2) = nodd - 1
quarter_neighbour(get_quarter_idx(nodu, 2), 2) = nodd - 1
elseif (mu == 2 * md) then
quarter_line(get_quarter_idx(nodd, 4), 2) = l - 1
quarter_line(get_quarter_idx(nodu, 1), 2) = l - 1
quarter_line(get_quarter_idx(nodu, 2), 2) = l - 1
quarter_neighbour(get_quarter_idx(nodd, 4), 2) = nodu - 1
quarter_neighbour(get_quarter_idx(nodu, 1), 2) = nodd - 1
quarter_neighbour(get_quarter_idx(nodu, 2), 2) = nodd - 1
endif
endif
enddo

do nb=1,n2dobc !nodobc
l = liutot + livtot + nb
nodd = line(l, 1) + 1
nodu = line(l, 2) + 1
select case(kcu(l))
case(LINE_2D_BOUNDARY_WEST)
quarter_line(get_quarter_idx(nodu, 1), 1) = l - 1
quarter_line(get_quarter_idx(nodu, 3), 1) = l
quarter_neighbour(get_quarter_idx(nodu, 1), 1) = nodd - 1
quarter_neighbour(get_quarter_idx(nodu, 3), 1) = nodd - 1
case(LINE_2D_BOUNDARY_EAST)
quarter_line(get_quarter_idx(nodd, 2), 1) = l - 1
quarter_line(get_quarter_idx(nodd, 4), 1) = l - 1
quarter_neighbour(get_quarter_idx(nodd, 2), 1) = nodu - 1
quarter_neighbour(get_quarter_idx(nodd, 4), 1) = nodu - 1
case(LINE_2D_BOUNDARY_SOUTH)
quarter_line(get_quarter_idx(nodu, 1), 2) = l - 1
quarter_line(get_quarter_idx(nodu, 2), 2) = l - 1
quarter_neighbour(get_quarter_idx(nodu, 1), 2) = nodd - 1
quarter_neighbour(get_quarter_idx(nodu, 2), 2) = nodd - 1
case(LINE_2D_BOUNDARY_NORTH)
quarter_line(get_quarter_idx(nodd, 3), 2) = l - 1
quarter_line(get_quarter_idx(nodd, 4), 2) = l - 1
quarter_neighbour(get_quarter_idx(nodd, 3), 2) = nodu - 1
quarter_neighbour(get_quarter_idx(nodd, 4), 2) = nodu - 1
end select
enddo

end subroutine set_quarter_admin


function get_quarter_idx(nod, quadrant) result(idx)
! This function returns the quarter index based on the node and its quadrant index.
!!!!!!!!!!!!!!!!
!! 3 !! 4 !!
!!!!!!!!!!!!!!!!
!! 1 !! 2 !!
!!!!!!!!!!!!!!!!
integer, intent(in) :: nod
integer, intent(in) :: quadrant
integer :: idx

idx = 4 * (nod - 1) + quadrant

end function get_quarter_idx

end module m_cells
26 changes: 26 additions & 0 deletions libthreedigrid/parameters.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,5 +13,31 @@ module parameters
NODE_2D_BOUND = 5,&
NODE_2D_GW_BOUND = 6,&
NODE_1D_BOUND = 7

integer, parameter :: LINE_1D_EMBEDDED = 0,&
LINE_1D_ISOLATED = 1,&
LINE_1D_CONNECTED = 2,&
LINE_1D_LONG_CRESTED = 3,&
LINE_1D_SHORT_CRESTED = 4,&
LINE_1D_DOUBLE_CONNECTED = 5,&
LINE_1D2D_SINGLE_CONNECTED_CLOSED = 51,&
LINE_1D2D_SINGLE_CONNECTED_OPEN_WATER = 52,&
LINE_1D2D_DOUBLE_CONNECTED_CLOSED = 53,&
LINE_1D2D_DOUBLE_CONNECTED_OPEN_WATER = 54,&
LINE_1D2D_POSSIBLE_BREACH = 55,&
LINE_1D2D_ACTIVE_BREACH = 56,&
LINE_1D2D_GROUNDWATER = 57,&
LINE_2D = 100,&
LINE_2D_OBSTACLE = 101,&
LINE_2D_VERTICAL = 150,&
LINE_2D_GROUNDWATER = -150,&
LINE_2D_BOUNDARY_WEST = 200,&
LINE_2D_BOUNDARY_EAST = 300,&
LINE_2D_BOUNDARY_SOUTH = 400,&
LINE_2D_BOUNDARY_NORTH = 500,&
LINE_2D_GROUNDWATER_BOUNDARY_WEST = 600,&
LINE_2D_GROUNDWATER_BOUNDARY_EAST = 700,&
LINE_2D_GROUNDWATER_BOUNDARY_SOUTH = 800,&
LINE_2D_GROUNDWATER_BOUNDARY_NORTH = 900

end module parameters
3 changes: 3 additions & 0 deletions threedigrid_builder/application.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -81,6 +81,8 @@ def _make_gridadmin(
node_id_counter=node_id_counter,
line_id_counter=line_id_counter,
)
grid.set_quarter_administration(quadtree)

if grid.meta.has_groundwater:
grid.add_groundwater(
grid.meta.has_groundwater_flow, node_id_counter, line_id_counter
Expand All @@ -93,6 +95,7 @@ def _make_gridadmin(
node_id_counter,
line_id_counter,
)

dem_average_areas = db.get_dem_average_areas()
grid.set_dem_averaged_cells(dem_average_areas)

Expand Down
1 change: 1 addition & 0 deletions threedigrid_builder/base/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,5 +5,6 @@
from .linestrings import * # NOQA
from .nodes import * # NOQA
from .pumps import * # NOQA
from .quarters import * # NOQA
from .settings import * # NOQA
from .surfaces import * # NOQA
17 changes: 17 additions & 0 deletions threedigrid_builder/base/quarters.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,17 @@
from typing import Tuple

from .array import Array

__all__ = ["Quarters"]


class Quarter:
id: int
line: Tuple[int, int]
neighbour_node: Tuple[int, int]


class Quarters(Array[Quarter]):
"""Calculation quarters."""

pass
13 changes: 13 additions & 0 deletions threedigrid_builder/grid/grid.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@
Nodes,
PointsOnLine,
Pumps,
Quarters,
SurfaceMaps,
Surfaces,
)
Expand Down Expand Up @@ -172,6 +173,7 @@ def __init__(
breaches=None,
meta=None,
quadtree_stats=None,
quarters: Optional[Quarters] = None,
):
if not isinstance(nodes, Nodes):
raise TypeError(f"Expected Nodes instance, got {type(nodes)}")
Expand Down Expand Up @@ -206,8 +208,15 @@ def __init__(
raise TypeError(f"Expected Obstacles instance, got {type(obstacles)}")
if breaches is not None and not isinstance(breaches, PotentialBreaches):
raise TypeError(f"Expected Breaches instance, got {type(breaches)}")

if quarters is None:
quarters = Quarters(id=[])
elif not isinstance(quarters, Quarters):
raise TypeError(f"Expected Quarters instance, got {type(quarters)}")

self.nodes = nodes
self.lines = lines
self.quarters = quarters
self.meta = meta
self.surfaces = surfaces
self.surface_maps = surface_maps
Expand Down Expand Up @@ -605,6 +614,10 @@ def set_boundary_conditions_2d(
self.nodes += nodes
self.lines += lines

def set_quarter_administration(self, quadtree):
"""Sets the quarter cell administration for lines and neighbouring calculation cells."""
self.quarters = quadtree.get_quarters_admin(self.nodes, self.lines)

def set_pumps(self, pumps):
"""Set the pumps on this grid object
Expand Down
40 changes: 39 additions & 1 deletion threedigrid_builder/grid/quadtree.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@

import numpy as np

from threedigrid_builder.base import Lines, Nodes
from threedigrid_builder.base import Lines, Nodes, Quarters
from threedigrid_builder.constants import LineType, NodeType
from threedigrid_builder.exceptions import SchematisationError

Expand Down Expand Up @@ -247,3 +247,41 @@ def get_nodes_lines(self, area_mask, node_id_counter, line_id_counter):
)

return nodes, lines

def get_quarters_admin(self, nodes, lines):
"""Compute 2D openwater quarter administration based on computed Quadtree.
Determine the connecting flowlines for a quarter in u and v direction.
Determine the adjacent calculation nodes for a quarter in u and v direction.
Return:
Quarters object
"""

n_2d_nodes = np.count_nonzero(
np.isin(nodes.node_type, NodeType.NODE_2D_OPEN_WATER)
)
n_2d_bnd_nodes = np.count_nonzero(
np.isin(nodes.node_type, NodeType.NODE_2D_BOUNDARIES)
)

quarter_line = np.full((4 * n_2d_nodes, 2), -9999, dtype=np.int32, order="F")
neighbour_node = np.full((4 * n_2d_nodes, 2), -9999, dtype=np.int32, order="F")

m_cells.set_quarter_admin(
nodes.nodk,
nodes.nodm,
nodes.nodn,
lines.line,
lines.kcu,
quarter_line,
neighbour_node,
self.n_lines_u,
self.n_lines_v,
n_2d_bnd_nodes,
)

return Quarters(
id=np.arange(0, 4 * n_2d_nodes),
line=quarter_line,
neighbour_node=neighbour_node,
)
Loading

0 comments on commit 9fcf3a9

Please sign in to comment.