MODFLOW 6  version 6.7.0.dev3
USGS Modular Hydrologic Model
gridconnectionmodule Module Reference

Refactoring issues towards parallel: More...

Data Types

type  gridconnectiontype
 This class is used to construct the connections object for the interface model's spatial discretization/grid. More...
 

Functions/Subroutines

subroutine construct (this, model, nrOfPrimaries, connectionName)
 Construct the GridConnection and allocate the data structures for the primary connections. More...
 
subroutine connectprimaryexchange (this, primEx)
 Make connections for the primary exchange. More...
 
subroutine connectcell (this, idx1, v_model1, idx2, v_model2)
 Connect neighboring cells at the interface by storing them in the boundary cell and connected cell arrays. More...
 
subroutine addtoregionalmodels (this, v_model)
 Add a model to a list of all regional models. More...
 
subroutine extendconnection (this)
 Extend the connection topology to deal with higher levels of connectivity (neighbors-of-neighbors, etc.) More...
 
subroutine buildconnections (this)
 Builds a sparse matrix holding all cell connections,. More...
 
recursive subroutine addneighbors (this, cellNbrs, depth, mask, interior)
 
subroutine addremoteneighbors (this, cellNbrs, mask)
 Add cell neighbors across models using the stored exchange data structures. More...
 
subroutine addneighborcell (this, cellNbrs, newNbrIdx, v_nbr_model, mask)
 Add neighboring cell to tree structure. More...
 
recursive subroutine registerinterfacecells (this, cellWithNbrs)
 Recursively set interface cell indexes and. More...
 
subroutine addtoglobalmap (this, ifaceIdx, cell)
 Add entry to lookup table, inflating when necessary. More...
 
subroutine compressglobalmap (this)
 Compress lookup table to get rid of unused entries. More...
 
subroutine sortinterfacegrid (this)
 Soft cell ids in the interface grid such that. More...
 
subroutine makeprimaryconnections (this, sparse)
 Add primary connections to the sparse data structure. More...
 
recursive subroutine connectneighborcells (this, cell, sparse)
 Recursively add higher order connections (from cells neighboring the primarily connected cells) to the. More...
 
subroutine fillconnectiondatainternal (this)
 Fill connection data (ihc, cl1, ...) for. More...
 
subroutine fillconnectiondatafromexchanges (this)
 Fill connection data (ihc, cl1, ...) for. More...
 
subroutine createconnectionmask (this)
 Create the connection masks. More...
 
recursive subroutine maskinternalconnections (this, cell, nbrCell, level)
 Recursively mask connections, increasing the level as we go. More...
 
subroutine setmaskonconnection (this, cell, nbrCell, level)
 Set a mask on the connection from a cell to its neighbor, (and not the transposed!) not overwriting the current level. More...
 
integer(i4b) function getinterfaceindexbycell (this, cell)
 Get interface index from global cell. More...
 
integer(i4b) function getinterfaceindexbyindexmodel (this, index, v_model)
 Get interface index from a model pointer and the local index. More...
 
integer(i4b) function get_regional_offset (this, v_model)
 Get the offset for a regional model. More...
 
subroutine allocatescalars (this)
 Allocate scalar data. More...
 
subroutine getdiscretization (this, disu)
 Sets the discretization (DISU) after all preprocessing by this grid connection has been done,. More...
 
subroutine buildinterfacemap (this)
 Build interface map object for outside use. More...
 
subroutine destroy (this)
 Deallocate grid connection resources. More...
 
logical function isperiodic (this, n, m)
 Test if the connection between nodes within. More...
 

Variables

integer(i4b), parameter initnrneighbors = 7
 

Detailed Description

  • remove camelCase

Function/Subroutine Documentation

◆ addneighborcell()

subroutine gridconnectionmodule::addneighborcell ( class(gridconnectiontype), intent(in)  this,
type(cellwithnbrstype), intent(inout)  cellNbrs,
integer(i4b), intent(in)  newNbrIdx,
class(virtualmodeltype), pointer  v_nbr_model,
type(globalcelltype), optional  mask 
)
private
Parameters
[in]thisthis grid connection instance
[in,out]cellnbrsthe root cell which to add to
[in]newnbridxthe neighboring cell's index
v_nbr_modelthe model where the new neighbor lives
maskdon't add connections to this cell (optional)

Definition at line 483 of file GridConnection.f90.

484  class(GridConnectionType), intent(in) :: this !< this grid connection instance
485  type(CellWithNbrsType), intent(inout) :: cellNbrs !< the root cell which to add to
486  integer(I4B), intent(in) :: newNbrIdx !< the neighboring cell's index
487  class(VirtualModelType), pointer :: v_nbr_model !< the model where the new neighbor lives
488  type(GlobalCellType), optional :: mask !< don't add connections to this cell (optional)
489 
490  if (present(mask)) then
491  if (newnbridx == mask%index .and. mask%v_model == v_nbr_model) then
492  return
493  end if
494  end if
495 
496  call cellnbrs%addNbrCell(newnbridx, v_nbr_model)
497 

◆ addneighbors()

recursive subroutine gridconnectionmodule::addneighbors ( class(gridconnectiontype), intent(inout)  this,
type(cellwithnbrstype), intent(inout), target  cellNbrs,
integer(i4b), intent(inout)  depth,
type(globalcelltype), optional  mask,
logical(lgp)  interior 
)
private
Parameters
[in,out]thisthis grid connection
[in,out]cellnbrscell to add to
[in,out]depthcurrent depth (typically decreases in recursion)
maskmask to excluded back-and-forth connection between cells
interiorwhen true, we are adding from the exchange back into the model

Definition at line 386 of file GridConnection.f90.

387  use simmodule, only: ustop
388  class(GridConnectionType), intent(inout) :: this !< this grid connection
389  type(CellWithNbrsType), intent(inout), target :: cellNbrs !< cell to add to
390  integer(I4B), intent(inout) :: depth !< current depth (typically decreases in recursion)
391  type(GlobalCellType), optional :: mask !< mask to excluded back-and-forth connection between cells
392  logical(LGP) :: interior !< when true, we are adding from the exchange back into the model
393  ! local
394  type(GlobalCellType), pointer :: cell
395  integer(I4B) :: ipos, ipos_start, ipos_end
396  integer(I4B) :: nbrIdx, inbr
397  integer(I4B) :: newDepth
398 
399  ! readability
400  cell => cellnbrs%cell
401 
402  ! if depth == 1, then we are not adding neighbors but use
403  ! the boundary and connected cell only
404  if (depth < 2) then
405  return
406  end if
407  newdepth = depth - 1
408 
409  ! find neighbors local to this cell by looping through grid connections
410  ipos_start = cell%v_model%con_ia%get(cell%index) + 1
411  ipos_end = cell%v_model%con_ia%get(cell%index + 1) - 1
412  do ipos = ipos_start, ipos_end
413  nbridx = cell%v_model%con_ja%get(ipos)
414  call this%addNeighborCell(cellnbrs, nbridx, cellnbrs%cell%v_model, mask)
415  end do
416 
417  ! add remote nbr using the data from the exchanges
418  call this%addRemoteNeighbors(cellnbrs, mask)
419 
420  ! now find nbr-of-nbr
421  do inbr = 1, cellnbrs%nrOfNbrs
422 
423  ! are we leaving the model through another exchange?
424  if (interior .and. cellnbrs%cell%v_model == this%model) then
425  if (.not. cellnbrs%neighbors(inbr)%cell%v_model == this%model) then
426  ! decrement by 1, because the connection we are crossing is not
427  ! calculated by this interface
428  newdepth = newdepth - 1
429  end if
430  end if
431  ! and add neighbors with the new depth
432  call this%addNeighbors(cellnbrs%neighbors(inbr), newdepth, &
433  cellnbrs%cell, interior)
434  end do
435 
This module contains simulation methods.
Definition: Sim.f90:10
subroutine, public ustop(stopmess, ioutlocal)
Stop the simulation.
Definition: Sim.f90:312
Here is the call graph for this function:

◆ addremoteneighbors()

subroutine gridconnectionmodule::addremoteneighbors ( class(gridconnectiontype), intent(inout)  this,
type(cellwithnbrstype), intent(inout)  cellNbrs,
type(globalcelltype), optional  mask 
)
Parameters
[in,out]thisthis grid connection instance
[in,out]cellnbrscell to add to
maska mask to exclude back-and-forth connections

Definition at line 440 of file GridConnection.f90.

441  class(GridConnectionType), intent(inout) :: this !< this grid connection instance
442  type(CellWithNbrsType), intent(inout) :: cellNbrs !< cell to add to
443  type(GlobalCellType), optional :: mask !< a mask to exclude back-and-forth connections
444  ! local
445  integer(I4B) :: ix, iexg
446  class(VirtualExchangeType), pointer :: v_exchange
447  class(VirtualModelType), pointer :: v_m1, v_m2
448 
449  ! loop over all exchanges
450  do ix = 1, this%haloExchanges%size
451 
452  v_exchange => get_virtual_exchange(this%haloExchanges%at(ix))
453  v_m1 => v_exchange%v_model1
454  v_m2 => v_exchange%v_model2
455 
456  ! loop over n-m links in the exchange
457  if (cellnbrs%cell%v_model == v_m1) then
458  do iexg = 1, v_exchange%nexg%get()
459  if (v_exchange%nodem1%get(iexg) == cellnbrs%cell%index) then
460  ! we have a link, now add foreign neighbor
461  call this%addNeighborCell( &
462  cellnbrs, v_exchange%nodem2%get(iexg), v_m2, mask)
463  end if
464  end do
465  end if
466  ! and the reverse
467  if (cellnbrs%cell%v_model == v_m2) then
468  do iexg = 1, v_exchange%nexg%get()
469  if (v_exchange%nodem2%get(iexg) == cellnbrs%cell%index) then
470  ! we have a link, now add foreign neighbor
471  call this%addNeighborCell( &
472  cellnbrs, v_exchange%nodem1%get(iexg), v_m1, mask)
473  end if
474  end do
475  end if
476 
477  end do
478 
Here is the call graph for this function:

◆ addtoglobalmap()

subroutine gridconnectionmodule::addtoglobalmap ( class(gridconnectiontype), intent(inout)  this,
integer(i4b), intent(in)  ifaceIdx,
type(globalcelltype), intent(in)  cell 
)
private
Parameters
[in,out]thisthis grid connection instance
[in]ifaceidxunique idx in the interface grid
[in]cellthe global cell

Definition at line 529 of file GridConnection.f90.

530  class(GridConnectionType), intent(inout) :: this !< this grid connection instance
531  integer(I4B), intent(in) :: ifaceIdx !< unique idx in the interface grid
532  type(GlobalCellType), intent(in) :: cell !< the global cell
533  ! local
534  integer(I4B) :: i, currentSize, newSize
535  type(GlobalCellType), dimension(:), pointer :: tempMap
536 
537  ! inflate?
538  currentsize = size(this%idxToGlobal)
539  if (ifaceidx > currentsize) then
540  newsize = nint(1.5 * currentsize)
541  allocate (tempmap(newsize))
542  do i = 1, currentsize
543  tempmap(i) = this%idxToGlobal(i)
544  end do
545 
546  deallocate (this%idxToGlobal)
547  this%idxToGlobal => tempmap
548  end if
549 
550  this%idxToGlobal(ifaceidx) = cell
551 

◆ addtoregionalmodels()

subroutine gridconnectionmodule::addtoregionalmodels ( class(gridconnectiontype), intent(inout)  this,
class(virtualmodeltype), pointer  v_model 
)
private
Parameters
[in,out]thisthis grid connection
v_modelthe model to add to the region

Definition at line 219 of file GridConnection.f90.

220  class(GridConnectionType), intent(inout) :: this !< this grid connection
221  class(VirtualModelType), pointer :: v_model !< the model to add to the region
222  ! local
223  class(*), pointer :: vm_obj
224 
225  vm_obj => v_model
226  if (.not. this%regionalModels%ContainsObject(vm_obj)) then
227  call this%regionalModels%Add(vm_obj)
228  end if
229 

◆ allocatescalars()

subroutine gridconnectionmodule::allocatescalars ( class(gridconnectiontype this)
private
Parameters
thisthis grid connection instance

Definition at line 981 of file GridConnection.f90.

983  class(GridConnectionType) :: this !< this grid connection instance
984 
985  call mem_allocate(this%nrOfBoundaryCells, 'NRBNDCELLS', this%memoryPath)
986  call mem_allocate(this%indexCount, 'IDXCOUNT', this%memoryPath)
987  call mem_allocate(this%nrOfCells, 'NRCELLS', this%memoryPath)
988 

◆ buildconnections()

subroutine gridconnectionmodule::buildconnections ( class(gridconnectiontype), intent(inout)  this)
private
Parameters
[in,out]thisthis grid connection instance

Definition at line 289 of file GridConnection.f90.

290  class(GridConnectionType), intent(inout) :: this !< this grid connection instance
291  ! local
292  integer(I4B) :: i, icell, iconn
293  integer(I4B), dimension(:), allocatable :: nnz
294  type(SparseMatrix), pointer :: sparse
295  integer(I4B) :: ierror
296  type(ConnectionsType), pointer :: conn
297  class(VirtualModelType), pointer :: vm
298  character(len=LINELENGTH) :: warnmsg
299 
300  ! Recursively generate interface cell indices, fill map to global cells,
301  ! and add to region lookup table
302  this%indexCount = 0
303  do icell = 1, this%nrOfBoundaryCells
304  call this%registerInterfaceCells(this%boundaryCells(icell))
305  end do
306  do icell = 1, this%nrOfBoundaryCells
307  call this%registerInterfaceCells(this%connectedCells(icell))
308  end do
309  this%nrOfCells = this%indexCount
310 
311  ! compress lookup table
312  call this%compressGlobalMap()
313 
314  ! prepare to build vertex and cell2d info, or not
315  do i = 1, this%regionalModels%Count()
316  vm => get_virtual_model_from_list(this%regionalModels, i)
317  if (vm%dis_icondir%get() == 0) this%icondir = 0
318  end do
319 
320  if (this%icondir == 0 .and. this%model%dis%icondir > 0) then
321  write (warnmsg, '(3a)') 'Exchange ', trim(this%primaryExchange%name), &
322  " will not use vertex and/or cell coordinate data because not all &
323  &connected models have those data specified"
324  call store_warning(warnmsg)
325  end if
326 
327  ! sort interface indexes such that 'n > m' means 'n below m'
328  call this%sortInterfaceGrid()
329 
330  ! allocate a map from interface index to global coordinates
331  call mem_allocate(this%idxToGlobalIdx, this%nrOfCells, &
332  'IDXTOGLOBALIDX', this%memoryPath)
333 
334  ! create sparse data structure, to temporarily hold connections
335  allocate (sparse)
336  allocate (nnz(this%nrOfCells))
337  nnz = initnrneighbors + 1
338  call sparse%init(this%nrOfCells, this%nrOfCells, nnz)
339 
340  ! now (recursively) add connections to sparse, start with
341  ! the primary connections (n-m from the exchange files)
342  call this%makePrimaryConnections(sparse)
343  ! then into own domain
344  do icell = 1, this%nrOfBoundaryCells
345  call this%connectNeighborCells(this%boundaryCells(icell), sparse)
346  end do
347  ! and same for the neighbors of connected cells
348  do icell = 1, this%nrOfBoundaryCells
349  call this%connectNeighborCells(this%connectedCells(icell), sparse)
350  end do
351 
352  ! create connections object
353  allocate (this%connections)
354  conn => this%connections
355  call conn%allocate_scalars(this%memoryPath)
356  conn%nodes = this%nrOfCells
357  conn%nja = sparse%nnz
358  conn%njas = (conn%nja - conn%nodes) / 2
359  call conn%allocate_arrays()
360  do iconn = 1, conn%njas
361  conn%anglex(iconn) = -999.
362  end do
363 
364  ! fill connection from sparse
365  call sparse%filliaja(conn%ia, conn%ja, ierror)
366  if (ierror /= 0) then
367  write (*, *) 'Error filling ia/ja in GridConnection: terminating...'
368  call ustop()
369  end if
370  call fillisym(conn%nodes, conn%nja, conn%ia, conn%ja, conn%isym)
371  call filljas(conn%nodes, conn%nja, conn%ia, conn%ja, conn%isym, conn%jas)
372  call sparse%destroy()
373 
374  ! fill connection data (ihc, cl1, cl2, etc.) using data
375  ! from models and exchanges
376  call this%fillConnectionDataInternal()
377  call this%fillConnectionDataFromExchanges()
378 
379  ! set the masks on connections
380  call this%createConnectionMask()
381 
Here is the call graph for this function:

◆ buildinterfacemap()

subroutine gridconnectionmodule::buildinterfacemap ( class(gridconnectiontype this)
Parameters
thisthis grid connection

Definition at line 1068 of file GridConnection.f90.

1070  use stlvecintmodule
1071  class(GridConnectionType) :: this !< this grid connection
1072  ! local
1073  integer(I4B) :: i, j, iloc, jloc
1074  integer(I4B) :: im, ix, mid, n
1075  integer(I4B) :: ipos, ipos_model
1076  type(STLVecInt) :: model_ids
1077  type(STLVecInt) :: src_idx_tmp, tgt_idx_tmp, sign_tmp
1078  class(VirtualExchangeType), pointer :: v_exg
1079  class(VirtualModelType), pointer :: vm
1080  class(VirtualModelType), pointer :: v_model1, v_model2
1081  type(InterfaceMapType), pointer :: imap
1082  integer(I4B), dimension(:), pointer, contiguous :: ia_ptr, ja_ptr
1083 
1084  allocate (this%interfaceMap)
1085  imap => this%interfaceMap
1086 
1087  ! first get the participating models
1088  call model_ids%init()
1089  do i = 1, this%nrOfCells
1090  call model_ids%push_back_unique(this%idxToGlobal(i)%v_model%id)
1091  end do
1092 
1093  ! initialize the map
1094  call imap%init(model_ids%size, this%haloExchanges%size)
1095 
1096  ! for each model part of this interface, ...
1097  do im = 1, model_ids%size
1098  mid = model_ids%at(im)
1099  imap%model_ids(im) = mid
1100  vm => get_virtual_model(mid)
1101  imap%model_names(im) = vm%name
1102  call src_idx_tmp%init()
1103  call tgt_idx_tmp%init()
1104 
1105  ! store the node map for this model
1106  do i = 1, this%nrOfCells
1107  if (mid == this%idxToGlobal(i)%v_model%id) then
1108  call src_idx_tmp%push_back(this%idxToGlobal(i)%index)
1109  call tgt_idx_tmp%push_back(i)
1110  end if
1111  end do
1112 
1113  ! and copy into interface map
1114  allocate (imap%node_maps(im)%src_idx(src_idx_tmp%size))
1115  allocate (imap%node_maps(im)%tgt_idx(tgt_idx_tmp%size))
1116  do i = 1, src_idx_tmp%size
1117  imap%node_maps(im)%src_idx(i) = src_idx_tmp%at(i)
1118  imap%node_maps(im)%tgt_idx(i) = tgt_idx_tmp%at(i)
1119  end do
1120 
1121  call src_idx_tmp%destroy()
1122  call tgt_idx_tmp%destroy()
1123 
1124  ! and for connections
1125  call src_idx_tmp%init()
1126  call tgt_idx_tmp%init()
1127 
1128  ! store the connection map for this model
1129  do i = 1, this%nrOfCells
1130  if (mid /= this%idxToGlobal(i)%v_model%id) cycle
1131  do ipos = this%connections%ia(i), this%connections%ia(i + 1) - 1
1132  j = this%connections%ja(ipos)
1133  if (mid /= this%idxToGlobal(j)%v_model%id) cycle
1134 
1135  ! i and j are now in same model (mid)
1136  iloc = this%idxToGlobal(i)%index
1137  jloc = this%idxToGlobal(j)%index
1138  ia_ptr => this%idxToGlobal(i)%v_model%con_ia%get_array()
1139  ja_ptr => this%idxToGlobal(i)%v_model%con_ja%get_array()
1140  ipos_model = getcsrindex(iloc, jloc, ia_ptr, ja_ptr)
1141  call src_idx_tmp%push_back(ipos_model)
1142  call tgt_idx_tmp%push_back(ipos)
1143  end do
1144  end do
1145 
1146  ! copy into interface map
1147  allocate (imap%conn_maps(im)%src_idx(src_idx_tmp%size))
1148  allocate (imap%conn_maps(im)%tgt_idx(tgt_idx_tmp%size))
1149  do i = 1, src_idx_tmp%size
1150  imap%conn_maps(im)%src_idx(i) = src_idx_tmp%at(i)
1151  imap%conn_maps(im)%tgt_idx(i) = tgt_idx_tmp%at(i)
1152  end do
1153 
1154  call src_idx_tmp%destroy()
1155  call tgt_idx_tmp%destroy()
1156 
1157  end do
1158 
1159  call model_ids%destroy()
1160 
1161  ! for each exchange that is part of this interface
1162  do ix = 1, this%haloExchanges%size
1163 
1164  ! all exchanges in this list should have at
1165  ! least one relevant connection for this map
1166  v_exg => get_virtual_exchange(this%haloExchanges%at(ix))
1167  v_model1 => v_exg%v_model1
1168  v_model2 => v_exg%v_model2
1169 
1170  imap%exchange_ids(ix) = v_exg%id
1171  imap%exchange_names(ix) = v_exg%name
1172 
1173  call src_idx_tmp%init()
1174  call tgt_idx_tmp%init()
1175  call sign_tmp%init()
1176 
1177  do n = 1, v_exg%nexg%get()
1178  i = this%getInterfaceIndex(v_exg%nodem1%get(n), v_model1)
1179  j = this%getInterfaceIndex(v_exg%nodem2%get(n), v_model2)
1180  if (i == -1 .or. j == -1) cycle ! not all exchange nodes are part of the interface
1181  ipos = this%connections%getjaindex(i, j)
1182  if (ipos == 0) then
1183  ! this can typically happen at corner points for a larger
1184  ! stencil (XT3D), when both i and j are included in the
1185  ! interface grid as leaf nodes, but their connection is not.
1186  ! (c.f. 'test_gwf_ifmod_mult_exg.py')
1187  cycle
1188  end if
1189  call src_idx_tmp%push_back(n)
1190  call tgt_idx_tmp%push_back(ipos)
1191  call sign_tmp%push_back(1)
1192 
1193  ! and the reverse connection:
1194  call src_idx_tmp%push_back(n)
1195  call tgt_idx_tmp%push_back(this%connections%isym(ipos))
1196  call sign_tmp%push_back(-1)
1197  end do
1198 
1199  allocate (imap%exchange_maps(ix)%src_idx(src_idx_tmp%size))
1200  allocate (imap%exchange_maps(ix)%tgt_idx(tgt_idx_tmp%size))
1201  allocate (imap%exchange_maps(ix)%sign(sign_tmp%size))
1202  do i = 1, src_idx_tmp%size
1203  imap%exchange_maps(ix)%src_idx(i) = src_idx_tmp%at(i)
1204  imap%exchange_maps(ix)%tgt_idx(i) = tgt_idx_tmp%at(i)
1205  imap%exchange_maps(ix)%sign(i) = sign_tmp%at(i)
1206  end do
1207 
1208  call src_idx_tmp%destroy()
1209  call tgt_idx_tmp%destroy()
1210  call sign_tmp%destroy()
1211 
1212  end do
1213 
1214  ! set the primary exchange idx
1215  ! findloc cannot be used until gfortran 9...
1216  imap%prim_exg_idx = -1
1217  do i = 1, imap%nr_exchanges
1218  if (imap%exchange_names(i) == this%primaryExchange%name) then
1219  imap%prim_exg_idx = i
1220  exit
1221  end if
1222  end do
1223 
1224  ! sanity check
1225  ! do i = 1, interfaceMap%nr_models
1226  ! if (size(interfaceMap%node_map(i)%src_idx) == 0) then
1227  ! write(*,*) 'Error: empty node map in interface for ', &
1228  ! this%primaryExchange%name
1229  ! call ustop()
1230  ! end if
1231  ! if (size(interfaceMap%connection_map(i)%src_idx) == 0) then
1232  ! write(*,*) 'Error: empty connection map in interface for ', &
1233  ! this%primaryExchange%name
1234  ! call ustop()
1235  ! end if
1236  ! end do
1237  ! do i = 1, interfaceMap%nr_exchanges
1238  ! if (size(interfaceMap%exchange_map(i)%src_idx) == 0) then
1239  ! write(*,*) 'Error: empty exchange map in interface for ', &
1240  ! this%primaryExchange%name
1241  ! call ustop()
1242  ! end if
1243  ! end do
1244 
class(basemodeltype) function, pointer, public getbasemodelfromlist(list, idx)
Definition: BaseModel.f90:172
Highest level model type. All models extend this parent type.
Definition: BaseModel.f90:13
Here is the call graph for this function:

◆ compressglobalmap()

subroutine gridconnectionmodule::compressglobalmap ( class(gridconnectiontype), intent(inout)  this)
private
Parameters
[in,out]thisthis grid connection instance

Definition at line 556 of file GridConnection.f90.

557  class(GridConnectionType), intent(inout) :: this !< this grid connection instance
558  ! local
559  type(GlobalCellType), dimension(:), pointer :: tempMap
560 
561  if (size(this%idxToGlobal) > this%nrOfCells) then
562  allocate (tempmap(this%nrOfCells))
563  tempmap(1:this%nrOfCells) = this%idxToGlobal(1:this%nrOfCells)
564  deallocate (this%idxToGlobal)
565  allocate (this%idxToGlobal(this%nrOfCells))
566  this%idxToGlobal(1:this%nrOfCells) = tempmap(1:this%nrOfCells)
567  deallocate (tempmap)
568  end if
569 

◆ connectcell()

subroutine gridconnectionmodule::connectcell ( class(gridconnectiontype), intent(in)  this,
integer(i4b)  idx1,
class(virtualmodeltype), pointer  v_model1,
integer(i4b)  idx2,
class(virtualmodeltype), pointer  v_model2 
)
private
Parameters
[in]thisthis grid connection
idx1local index cell 1
v_model1model of cell 1
idx2local index cell 2
v_model2model of cell 2

Definition at line 182 of file GridConnection.f90.

183  class(GridConnectionType), intent(in) :: this !< this grid connection
184  integer(I4B) :: idx1 !< local index cell 1
185  class(VirtualModelType), pointer :: v_model1 !< model of cell 1
186  integer(I4B) :: idx2 !< local index cell 2
187  class(VirtualModelType), pointer :: v_model2 !< model of cell 2
188  ! local
189  type(GlobalCellType), pointer :: bnd_cell, conn_cell
190 
191  this%nrOfBoundaryCells = this%nrOfBoundaryCells + 1
192  if (this%nrOfBoundaryCells > size(this%boundaryCells)) then
193  write (*, *) 'Error: nr of cell connections exceeds '// &
194  'capacity in grid connection, terminating...'
195  call ustop()
196  end if
197 
198  bnd_cell => this%boundaryCells(this%nrOfBoundaryCells)%cell
199  conn_cell => this%connectedCells(this%nrOfBoundaryCells)%cell
200  if (v_model1 == this%model) then
201  bnd_cell%index = idx1
202  bnd_cell%v_model => v_model1
203  conn_cell%index = idx2
204  conn_cell%v_model => v_model2
205  else if (v_model2 == this%model) then
206  bnd_cell%index = idx2
207  bnd_cell%v_model => v_model2
208  conn_cell%index = idx1
209  conn_cell%v_model => v_model1
210  else
211  write (*, *) 'Error: unable to connect cells outside the model'
212  call ustop()
213  end if
214 
Here is the call graph for this function:

◆ connectneighborcells()

recursive subroutine gridconnectionmodule::connectneighborcells ( class(gridconnectiontype), intent(inout)  this,
type(cellwithnbrstype cell,
type(sparsematrix), pointer  sparse 
)
private
Parameters
[in,out]thisthis grid connection instance
cellthe cell whose connections is to be added
sparsethe sparse data structure to hold the connections

Definition at line 653 of file GridConnection.f90.

654  class(GridConnectionType), intent(inout) :: this !< this grid connection instance
655  type(CellWithNbrsType) :: cell !< the cell whose connections is to be added
656  type(SparseMatrix), pointer :: sparse !< the sparse data structure to hold the connections
657  ! local
658  integer(I4B) :: ifaceIdx, ifaceIdxNbr ! unique idx in the interface grid
659  integer(I4B) :: inbr
660 
661  ifaceidx = this%getInterfaceIndex(cell%cell)
662  do inbr = 1, cell%nrOfNbrs
663  ifaceidxnbr = this%getInterfaceIndex(cell%neighbors(inbr)%cell)
664 
665  call sparse%addconnection(ifaceidxnbr, ifaceidxnbr, 1)
666  call sparse%addconnection(ifaceidx, ifaceidxnbr, 1)
667  call sparse%addconnection(ifaceidxnbr, ifaceidx, 1)
668 
669  ! recurse
670  call this%connectNeighborCells(cell%neighbors(inbr), sparse)
671  end do
672 

◆ connectprimaryexchange()

subroutine gridconnectionmodule::connectprimaryexchange ( class(gridconnectiontype this,
class(disconnexchangetype), pointer  primEx 
)
private
Parameters
thisthis grid connection
primexthe primary exchange for this connection

Definition at line 161 of file GridConnection.f90.

162  class(GridConnectionType) :: this !< this grid connection
163  class(DisConnExchangeType), pointer :: primEx !< the primary exchange for this connection
164  ! local
165  integer(I4B) :: iconn
166 
167  ! store the primary exchange
168  this%primaryExchange => primex
169 
170  ! connect the cells
171  do iconn = 1, primex%nexg
172  call this%connectCell(primex%nodem1(iconn), primex%v_model1, &
173  primex%nodem2(iconn), primex%v_model2)
174  end do
175 

◆ construct()

subroutine gridconnectionmodule::construct ( class(gridconnectiontype), intent(inout)  this,
class(numericalmodeltype), intent(in), pointer  model,
integer(i4b)  nrOfPrimaries,
character(len=*)  connectionName 
)

Definition at line 130 of file GridConnection.f90.

131  class(GridConnectionType), intent(inout) :: this !> this instance
132  class(NumericalModelType), pointer, intent(in) :: model !> the model for which the interface is constructed
133  integer(I4B) :: nrOfPrimaries !> the number of primary connections between the two models
134  character(len=*) :: connectionName !> the name, for memory management mostly
135  ! local
136  class(VirtualModelType), pointer :: v_model
137 
138  this%model => model
139  this%memoryPath = create_mem_path(connectionname, 'GC')
140 
141  call this%allocateScalars()
142 
143  allocate (this%boundaryCells(nrofprimaries))
144  allocate (this%connectedCells(nrofprimaries))
145  allocate (this%idxToGlobal(2 * nrofprimaries))
146 
147  v_model => get_virtual_model(model%id)
148  call this%addToRegionalModels(v_model)
149 
150  this%nrOfBoundaryCells = 0
151 
152  this%internalStencilDepth = 1
153  this%exchangeStencilDepth = 1
154  this%icondir = 1 ! start with assumption that vertex and cell2d data are available in all models
155  this%haloExchanges => null()
156 
Here is the call graph for this function:

◆ createconnectionmask()

subroutine gridconnectionmodule::createconnectionmask ( class(gridconnectiontype), intent(inout)  this)

The level indicates the nr of connections away from the remote neighbor, the diagonal term holds the negated value of their nearest connection. We end with setting

Parameters
[in,out]thisinstance of this grid connection

Definition at line 809 of file GridConnection.f90.

810  class(GridConnectionType), intent(inout) :: this !< instance of this grid connection
811  ! local
812  integer(I4B) :: icell, inbr, n, ipos
813  integer(I4B) :: level, newMask
814  type(CellWithNbrsType), pointer :: cell, nbrCell
815 
816  ! set all masks to zero to begin with
817  do ipos = 1, this%connections%nja
818  call this%connections%set_mask(ipos, 0)
819  end do
820 
821  ! remote connections remain masked
822  ! now set mask for exchange connections (level == 1)
823  level = 1
824  do icell = 1, this%nrOfBoundaryCells
825  call this%setMaskOnConnection(this%boundaryCells(icell), &
826  this%connectedCells(icell), level)
827  ! for cross-boundary connections, we need to apply the mask to both n-m and m-n,
828  ! because if the upper triangular one is disabled, its transposed (lower triangular)
829  ! counter part is skipped in the NPF calculation as well.
830  call this%setMaskOnConnection(this%connectedCells(icell), &
831  this%boundaryCells(icell), level)
832  end do
833 
834  ! now extend mask recursively into the internal domain (level > 1)
835  do icell = 1, this%nrOfBoundaryCells
836  cell => this%boundaryCells(icell)
837  do inbr = 1, cell%nrOfNbrs
838  nbrcell => this%boundaryCells(icell)%neighbors(inbr)
839  level = 2 ! this is incremented within the recursion
840  call this%maskInternalConnections(this%boundaryCells(icell), &
841  this%boundaryCells(icell)% &
842  neighbors(inbr), level)
843  end do
844  end do
845 
846  ! set normalized mask:
847  ! =1 for links with connectivity <= interior stencil depth
848  ! =0 otherwise
849  do n = 1, this%connections%nodes
850  ! set diagonals to zero
851  call this%connections%set_mask(this%connections%ia(n), 0)
852 
853  do ipos = this%connections%ia(n) + 1, this%connections%ia(n + 1) - 1
854  newmask = 0
855  if (this%connections%mask(ipos) > 0) then
856  if (this%connections%mask(ipos) < this%internalStencilDepth + 1) then
857  newmask = 1
858  end if
859  end if
860  ! set mask on off-diag
861  call this%connections%set_mask(ipos, newmask)
862  end do
863  end do
864 

◆ destroy()

subroutine gridconnectionmodule::destroy ( class(gridconnectiontype this)
Parameters
thisthis grid connection instance

Definition at line 1249 of file GridConnection.f90.

1251  class(GridConnectionType) :: this !< this grid connection instance
1252 
1253  call mem_deallocate(this%nrOfBoundaryCells)
1254  call mem_deallocate(this%indexCount)
1255  call mem_deallocate(this%nrOfCells)
1256 
1257  ! arrays
1258  deallocate (this%idxToGlobal)
1259  deallocate (this%boundaryCells)
1260  deallocate (this%connectedCells)
1261 
1262  call mem_deallocate(this%idxToGlobalIdx)
1263 

◆ extendconnection()

subroutine gridconnectionmodule::extendconnection ( class(gridconnectiontype), intent(inout)  this)
private

The following steps are taken:

  1. Recursively add interior neighbors (own model) up to the specified depth
  2. Recursively add exterior neighbors
  3. Allocate a (sparse) mapping table for the region
  4. Build connection object for the interface grid, and the mask
    Parameters
    [in,out]thisthis grid connection

Definition at line 241 of file GridConnection.f90.

242  class(GridConnectionType), intent(inout) :: this !< this grid connection
243  ! local
244  integer(I4B) :: remoteDepth, localDepth
245  integer(I4B) :: icell
246  integer(I4B) :: imod, regionSize, offset
247  class(VirtualModelType), pointer :: v_model
248  !integer(I4B), pointer :: nr_nodes
249 
250  ! we need (stencildepth-1) extra cells for the interior
251  remotedepth = this%exchangeStencilDepth
252  localdepth = 2 * this%internalStencilDepth - 1
253  if (localdepth < remotedepth) then
254  localdepth = remotedepth
255  end if
256 
257  ! first add the neighbors for the interior
258  ! (possibly extending into other models)
259  do icell = 1, this%nrOfBoundaryCells
260  call this%addNeighbors(this%boundaryCells(icell), localdepth, &
261  this%connectedCells(icell)%cell, .true.)
262  end do
263  ! and for the exterior
264  do icell = 1, this%nrOfBoundaryCells
265  call this%addNeighbors(this%connectedCells(icell), remotedepth, &
266  this%boundaryCells(icell)%cell, .false.)
267  end do
268 
269  ! set up mapping for the region (models participating in interface model grid)
270  allocate (this%regionalModelOffset(this%regionalModels%Count()))
271  regionsize = 0
272  offset = 0
273  do imod = 1, this%regionalModels%Count()
274  v_model => get_virtual_model_from_list(this%regionalModels, imod)
275  regionsize = regionsize + v_model%dis_nodes%get()
276  this%regionalModelOffset(imod) = offset
277  offset = offset + v_model%dis_nodes%get()
278  end do
279  ! init to -1, meaning 'interface index was not assigned yet'
280  allocate (this%region_to_iface_map(regionsize))
281  this%region_to_iface_map = -1
282 
283  call this%buildConnections()
284 
Here is the call graph for this function:

◆ fillconnectiondatafromexchanges()

subroutine gridconnectionmodule::fillconnectiondatafromexchanges ( class(gridconnectiontype), intent(inout)  this)
Parameters
[in,out]thisthis grid connection instance

Definition at line 735 of file GridConnection.f90.

736  use constantsmodule, only: dpio180
737  class(GridConnectionType), intent(inout) :: this !< this grid connection instance
738  ! local
739  integer(I4B) :: inx, iexg
740  integer(I4B) :: ipos, isym
741  integer(I4B) :: nOffset, mOffset, nIfaceIdx, mIfaceIdx
742  class(VirtualExchangeType), pointer :: v_exg
743  class(VirtualModelType), pointer :: v_m1, v_m2
744  type(ConnectionsType), pointer :: conn
745 
746  conn => this%connections
747 
748  do inx = 1, this%haloExchanges%size
749  v_exg => get_virtual_exchange(this%haloExchanges%at(inx))
750 
751  v_m1 => v_exg%v_model1
752  v_m2 => v_exg%v_model2
753 
754  if (v_exg%ianglex%get() > 0) then
755  conn%ianglex = 1
756  end if
757 
758  noffset = this%get_regional_offset(v_m1)
759  moffset = this%get_regional_offset(v_m2)
760  do iexg = 1, v_exg%nexg%get()
761  nifaceidx = this%region_to_iface_map(noffset + v_exg%nodem1%get(iexg))
762  mifaceidx = this%region_to_iface_map(moffset + v_exg%nodem2%get(iexg))
763  ! not all nodes from the exchanges are part of the interface grid
764  ! (think of exchanges between neighboring models, and their neighbors)
765  if (nifaceidx == -1 .or. mifaceidx == -1) then
766  cycle
767  end if
768 
769  ipos = conn%getjaindex(nifaceidx, mifaceidx)
770  ! (see prev. remark) sometimes the cells are in the interface grid,
771  ! but the connection isn't. This can happen for leaf nodes of the grid.
772  if (ipos == 0) then
773  ! no match, safely cycle
774  cycle
775  end if
776  isym = conn%jas(ipos)
777 
778  ! note: cl1 equals L_nm: the length from cell n to the shared
779  ! face with cell m (and cl2 analogously for L_mn)
780  if (nifaceidx < mifaceidx) then
781  conn%cl1(isym) = v_exg%cl1%get(iexg)
782  conn%cl2(isym) = v_exg%cl2%get(iexg)
783  if (v_exg%ianglex%get() > 0) then
784  conn%anglex(isym) = &
785  v_exg%auxvar%get(v_exg%ianglex%get(), iexg) * dpio180
786  end if
787  else
788  conn%cl1(isym) = v_exg%cl2%get(iexg)
789  conn%cl2(isym) = v_exg%cl1%get(iexg)
790  if (v_exg%ianglex%get() > 0) then
791  conn%anglex(isym) = mod(v_exg%auxvar%get(v_exg%ianglex%get(), iexg) &
792  + 180.0_dp, 360.0_dp) * dpio180
793  end if
794  end if
795  conn%hwva(isym) = v_exg%hwva%get(iexg)
796  conn%ihc(isym) = v_exg%ihc%get(iexg)
797 
798  end do
799  end do
800 
This module contains simulation constants.
Definition: Constants.f90:9
real(dp), parameter dpio180
real constant
Definition: Constants.f90:130
Here is the call graph for this function:

◆ fillconnectiondatainternal()

subroutine gridconnectionmodule::fillconnectiondatainternal ( class(gridconnectiontype), intent(inout)  this)
private
Parameters
[in,out]thisthis grid connection instance

Definition at line 677 of file GridConnection.f90.

678  use constantsmodule, only: dpi, dtwopi
679  class(GridConnectionType), intent(inout) :: this !< this grid connection instance
680  ! local
681  type(ConnectionsType), pointer :: conn
682  integer(I4B) :: n, m, ipos, isym, ipos_orig, isym_orig
683  type(GlobalCellType), pointer :: ncell, mcell
684  class(VirtualModelType), pointer :: v_m !< pointer to virtual model (readability)
685 
686  conn => this%connections
687 
688  do n = 1, conn%nodes
689  do ipos = conn%ia(n) + 1, conn%ia(n + 1) - 1
690  m = conn%ja(ipos)
691  if (n > m) cycle
692 
693  isym = conn%jas(ipos)
694  ncell => this%idxToGlobal(n)
695  mcell => this%idxToGlobal(m)
696  if (ncell%v_model == mcell%v_model) then
697 
698  ! for readability
699  v_m => ncell%v_model
700 
701  ! within same model, straight copy
702  ipos_orig = getcsrindex(ncell%index, mcell%index, &
703  v_m%con_ia%get_array(), v_m%con_ja%get_array())
704 
705  if (ipos_orig == 0) then
706  ! periodic boundary conditions can add connections between cells in
707  ! the same model, but they are dealt with through the exchange data
708  if (this%isPeriodic(ncell%index, mcell%index)) cycle
709 
710  ! this should not be possible
711  write (*, *) 'Error: cannot find cell connection in model grid'
712  call ustop()
713  end if
714 
715  isym_orig = v_m%con_jas%get(ipos_orig)
716  conn%hwva(isym) = v_m%con_hwva%get(isym_orig)
717  conn%ihc(isym) = v_m%con_ihc%get(isym_orig)
718  if (ncell%index < mcell%index) then
719  conn%cl1(isym) = v_m%con_cl1%get(isym_orig)
720  conn%cl2(isym) = v_m%con_cl2%get(isym_orig)
721  conn%anglex(isym) = v_m%con_anglex%get(isym_orig)
722  else
723  conn%cl1(isym) = v_m%con_cl2%get(isym_orig)
724  conn%cl2(isym) = v_m%con_cl1%get(isym_orig)
725  conn%anglex(isym) = mod(v_m%con_anglex%get(isym_orig) + dpi, dtwopi)
726  end if
727  end if
728  end do
729  end do
730 
real(dp), parameter dtwopi
real constant
Definition: Constants.f90:129
real(dp), parameter dpi
real constant
Definition: Constants.f90:128
Here is the call graph for this function:

◆ get_regional_offset()

integer(i4b) function gridconnectionmodule::get_regional_offset ( class(gridconnectiontype), intent(inout)  this,
class(virtualmodeltype), pointer  v_model 
)
private
Parameters
[in,out]thisthis grid connection instance
v_modelthe model to get the offset for
Returns
the index offset in the regional domain

Definition at line 960 of file GridConnection.f90.

961  class(GridConnectionType), intent(inout) :: this !< this grid connection instance
962  class(VirtualModelType), pointer :: v_model !< the model to get the offset for
963  integer(I4B) :: offset !< the index offset in the regional domain
964  ! local
965  integer(I4B) :: im
966  class(VirtualModelType), pointer :: vm
967 
968  offset = 0
969  do im = 1, this%regionalModels%Count()
970  vm => get_virtual_model_from_list(this%regionalModels, im)
971  if (vm == v_model) then
972  offset = this%regionalModelOffset(im)
973  return
974  end if
975  end do
976 
Here is the call graph for this function:

◆ getdiscretization()

subroutine gridconnectionmodule::getdiscretization ( class(gridconnectiontype this,
class(disutype), pointer  disu 
)
Parameters
thisthe grid connection
disuthe target disu object

Definition at line 994 of file GridConnection.f90.

996  use sparsemodule, only: sparsematrix
997  class(GridConnectionType) :: this !< the grid connection
998  class(DisuType), pointer :: disu !< the target disu object
999  ! local
1000  integer(I4B) :: icell, nrOfCells, idx
1001  class(VirtualModelType), pointer :: v_model
1002  real(DP) :: xglo, yglo
1003 
1004  ! the following is similar to dis_df
1005  nrofcells = this%nrOfCells
1006  disu%nodes = nrofcells
1007  disu%nvert = 0 ! TODO_MJR: vertex reconstruction is not supported yet
1008  disu%icondir = this%icondir
1009  disu%nodesuser = nrofcells
1010  disu%nja = this%connections%nja
1011 
1012  call disu%allocate_arrays()
1013 
1014  ! these are otherwise allocated when sourced
1015  call disu%allocate_arrays_mem()
1016 
1017  ! fill data
1018  do icell = 1, nrofcells
1019  idx = this%idxToGlobal(icell)%index
1020  disu%top(icell) = -huge(1.0_dp)
1021  disu%bot(icell) = -huge(1.0_dp)
1022  disu%area(icell) = -huge(1.0_dp)
1023  end do
1024 
1025  ! grid connections follow from GridConnection:
1026  disu%con => this%connections
1027  disu%njas = disu%con%njas
1028 
1029  if (this%icondir > 0) then
1030  ! copy cell x,y because available in all models
1031  do icell = 1, nrofcells
1032  idx = this%idxToGlobal(icell)%index
1033  v_model => this%idxToGlobal(icell)%v_model
1034 
1035  ! we are merging grids with possibly (likely) different origins,
1036  ! transform to global coordinates:
1037  call dis_transform_xy(v_model%dis_xc%get(idx), &
1038  v_model%dis_yc%get(idx), &
1039  v_model%dis_xorigin%get(), &
1040  v_model%dis_yorigin%get(), &
1041  v_model%dis_angrot%get(), &
1042  xglo, yglo)
1043 
1044  ! NB: usernodes equals internal nodes for interface
1045  disu%cellxy(1, icell) = xglo
1046  disu%xc(icell) = xglo
1047  disu%cellxy(2, icell) = yglo
1048  disu%yc(icell) = yglo
1049  end do
1050  else
1051  ! otherwise compress
1052  call mem_reallocate(disu%xc, 0, 'XC', disu%memoryPath)
1053  call mem_reallocate(disu%yc, 0, 'YC', disu%memoryPath)
1054  end if
1055 
1056  ! if vertices will be needed, it will look like this:
1057  !
1058  ! 1. determine total nr. of verts
1059  ! 2. allocate vertices list
1060  ! 3. create sparse
1061  ! 4. get vertex data per cell, add functions to base
1062  ! 5. add vertex (x,y) to list and connectivity to sparse
1063  ! 6. generate ia/ja from sparse
1064 
Here is the call graph for this function:

◆ getinterfaceindexbycell()

integer(i4b) function gridconnectionmodule::getinterfaceindexbycell ( class(gridconnectiontype), intent(inout)  this,
type(globalcelltype), intent(in)  cell 
)
private
Parameters
[in,out]thisthis grid connection instance
[in]cellthe global cell to get the interface index for
Returns
the index in the interface model

Definition at line 929 of file GridConnection.f90.

930  class(GridConnectionType), intent(inout) :: this !< this grid connection instance
931  type(GlobalCellType), intent(in) :: cell !< the global cell to get the interface index for
932  integer(I4B) :: iface_idx !< the index in the interface model
933  ! local
934  integer(I4B) :: offset, region_idx
935 
936  offset = this%get_regional_offset(cell%v_model)
937  region_idx = offset + cell%index
938  iface_idx = this%region_to_iface_map(region_idx)
939 
Here is the caller graph for this function:

◆ getinterfaceindexbyindexmodel()

integer(i4b) function gridconnectionmodule::getinterfaceindexbyindexmodel ( class(gridconnectiontype), intent(inout)  this,
integer(i4b)  index,
class(virtualmodeltype), pointer  v_model 
)
private
Parameters
[in,out]thisthis grid connection instance
indexthe local cell index
v_modelthe cell's model
Returns
the index in the interface model

Definition at line 944 of file GridConnection.f90.

945  class(GridConnectionType), intent(inout) :: this !< this grid connection instance
946  integer(I4B) :: index !< the local cell index
947  class(VirtualModelType), pointer :: v_model !< the cell's model
948  integer(I4B) :: iface_idx !< the index in the interface model
949  ! local
950  integer(I4B) :: offset, region_idx
951 
952  offset = this%get_regional_offset(v_model)
953  region_idx = offset + index
954  iface_idx = this%region_to_iface_map(region_idx)
955 
Here is the caller graph for this function:

◆ isperiodic()

logical function gridconnectionmodule::isperiodic ( class(gridconnectiontype), intent(in)  this,
integer(i4b), intent(in)  n,
integer(i4b), intent(in)  m 
)
Parameters
[in]thisthis grid connection instance
[in]nfirst node of the connection
[in]msecond node of the connection
Returns
true when periodic

Definition at line 1268 of file GridConnection.f90.

1269  class(GridConnectionType), intent(in) :: this !< this grid connection instance
1270  integer(I4B), intent(in) :: n !< first node of the connection
1271  integer(I4B), intent(in) :: m !< second node of the connection
1272  logical :: periodic !< true when periodic
1273  ! local
1274  integer(I4B) :: icell
1275 
1276  periodic = .false.
1277  do icell = 1, this%nrOfBoundaryCells
1278  if (.not. this%boundaryCells(icell)%cell%v_model == &
1279  this%connectedCells(icell)%cell%v_model) then
1280  cycle
1281  end if
1282 
1283  ! one way
1284  if (this%boundaryCells(icell)%cell%index == n) then
1285  if (this%connectedCells(icell)%cell%index == m) then
1286  periodic = .true.
1287  return
1288  end if
1289  end if
1290  ! or the other
1291  if (this%boundaryCells(icell)%cell%index == m) then
1292  if (this%connectedCells(icell)%cell%index == n) then
1293  periodic = .true.
1294  return
1295  end if
1296  end if
1297 
1298  end do
1299 

◆ makeprimaryconnections()

subroutine gridconnectionmodule::makeprimaryconnections ( class(gridconnectiontype), intent(inout)  this,
type(sparsematrix), pointer  sparse 
)
Parameters
[in,out]thisthis grid connection instance
sparsethe sparse data structure to hold the connections

Definition at line 628 of file GridConnection.f90.

629  class(GridConnectionType), intent(inout) :: this !< this grid connection instance
630  type(SparseMatrix), pointer :: sparse !< the sparse data structure to hold the connections
631  ! local
632  integer(I4B) :: icell
633  integer(I4B) :: ifaceIdx, ifaceIdxNbr
634 
635  do icell = 1, this%nrOfBoundaryCells
636  ifaceidx = this%getInterfaceIndex(this%boundaryCells(icell)%cell)
637  ifaceidxnbr = this%getInterfaceIndex(this%connectedCells(icell)%cell)
638 
639  ! add diagonals to sparse
640  call sparse%addconnection(ifaceidx, ifaceidx, 1)
641  call sparse%addconnection(ifaceidxnbr, ifaceidxnbr, 1)
642 
643  ! and cross terms
644  call sparse%addconnection(ifaceidx, ifaceidxnbr, 1)
645  call sparse%addconnection(ifaceidxnbr, ifaceidx, 1)
646  end do
647 

◆ maskinternalconnections()

recursive subroutine gridconnectionmodule::maskinternalconnections ( class(gridconnectiontype), intent(inout)  this,
type(cellwithnbrstype), intent(inout)  cell,
type(cellwithnbrstype), intent(inout)  nbrCell,
integer(i4b), intent(in)  level 
)
private
Parameters
[in,out]thisthis grid connection instance
[in,out]cellcell 1 in the connection to mask
[in,out]nbrcellcell 2 in the connection to mask

Definition at line 869 of file GridConnection.f90.

870  class(GridConnectionType), intent(inout) :: this !< this grid connection instance
871  type(CellWithNbrsType), intent(inout) :: cell !< cell 1 in the connection to mask
872  type(CellWithNbrsType), intent(inout) :: nbrCell !< cell 2 in the connection to mask
873  integer(I4B), intent(in) :: level
874  ! local
875  integer(I4B) :: inbr, newLevel
876 
877  ! only set the mask for internal connections, leaving the
878  ! others at 0
879  if (cell%cell%v_model == this%model .and. &
880  nbrcell%cell%v_model == this%model) then
881  ! this will set a mask on both diagonal, and both cross terms
882  call this%setMaskOnConnection(cell, nbrcell, level)
883  call this%setMaskOnConnection(nbrcell, cell, level)
884  end if
885 
886  ! recurse on nbrs-of-nbrs
887  newlevel = level + 1
888  do inbr = 1, nbrcell%nrOfNbrs
889  call this%maskInternalConnections(nbrcell, &
890  nbrcell%neighbors(inbr), &
891  newlevel)
892  end do
893 

◆ registerinterfacecells()

recursive subroutine gridconnectionmodule::registerinterfacecells ( class(gridconnectiontype), intent(inout)  this,
type(cellwithnbrstype cellWithNbrs 
)
private
Parameters
[in,out]thisthis grid connection instance
cellwithnbrsthe cell from where to start registering neighbors

Definition at line 502 of file GridConnection.f90.

503  class(GridConnectionType), intent(inout) :: this !< this grid connection instance
504  type(CellWithNbrsType) :: cellWithNbrs !< the cell from where to start registering neighbors
505  ! local
506  integer(I4B) :: offset, inbr
507  integer(I4B) :: regionIdx ! unique idx in the region (all connected models)
508  integer(I4B) :: ifaceIdx ! unique idx in the interface grid
509 
510  offset = this%get_regional_offset(cellwithnbrs%cell%v_model)
511  regionidx = offset + cellwithnbrs%cell%index
512  ifaceidx = this%getInterfaceIndex(cellwithnbrs%cell)
513  if (ifaceidx == -1) then
514  this%indexCount = this%indexCount + 1
515  ifaceidx = this%indexCount
516  call this%addToGlobalMap(ifaceidx, cellwithnbrs%cell)
517  this%region_to_iface_map(regionidx) = ifaceidx
518  end if
519 
520  ! and also for its neighbors
521  do inbr = 1, cellwithnbrs%nrOfNbrs
522  call this%registerInterfaceCells(cellwithnbrs%neighbors(inbr))
523  end do
524 

◆ setmaskonconnection()

subroutine gridconnectionmodule::setmaskonconnection ( class(gridconnectiontype), intent(inout)  this,
type(cellwithnbrstype), intent(inout)  cell,
type(cellwithnbrstype), intent(inout)  nbrCell,
integer(i4b), intent(in)  level 
)
private
Parameters
[in,out]thisthis grid connection instance
[in,out]cellcell 1 in the connection
[in,out]nbrcellcell 2 in the connection
[in]levelthe level value to set the mask to

Definition at line 899 of file GridConnection.f90.

900  class(GridConnectionType), intent(inout) :: this !< this grid connection instance
901  type(CellWithNbrsType), intent(inout) :: cell !< cell 1 in the connection
902  type(CellWithNbrsType), intent(inout) :: nbrCell !< cell 2 in the connection
903  integer(I4B), intent(in) :: level !< the level value to set the mask to
904  ! local
905  integer(I4B) :: ifaceIdx, ifaceIdxNbr
906  integer(I4B) :: iposdiag, ipos
907  integer(I4B) :: currentLevel
908 
909  ifaceidx = this%getInterfaceIndex(cell%cell)
910  ifaceidxnbr = this%getInterfaceIndex(nbrcell%cell)
911 
912  ! diagonal
913  iposdiag = this%connections%getjaindex(ifaceidx, ifaceidx)
914  currentlevel = this%connections%mask(iposdiag)
915  if (currentlevel == 0 .or. level < currentlevel) then
916  call this%connections%set_mask(iposdiag, level)
917  end if
918  ! cross term
919  ipos = this%connections%getjaindex(ifaceidx, ifaceidxnbr)
920  currentlevel = this%connections%mask(ipos)
921  if (currentlevel == 0 .or. level < currentlevel) then
922  call this%connections%set_mask(ipos, level)
923  end if
924 

◆ sortinterfacegrid()

subroutine gridconnectionmodule::sortinterfacegrid ( class(gridconnectiontype), intent(inout)  this)
private
Parameters
[in,out]thisthis grid connection instance

Definition at line 574 of file GridConnection.f90.

575  use gridsorting, only: quicksortgrid
576  class(GridConnectionType), intent(inout) :: this !< this grid connection instance
577  ! local
578  integer(I4B), dimension(:), allocatable :: newToOldIdx
579  integer(I4B), dimension(:), allocatable :: oldToNewIdx
580  integer(I4B) :: idxOld
581  integer(I4B) :: i
582  type(GlobalCellType), dimension(:), allocatable :: sortedGlobalMap
583  integer(I4B), dimension(:), allocatable :: sortedRegionMap
584  logical(LGP) :: z_only
585 
586  ! only use z coordinate for sorting when no cell2d info is available
587  z_only = (this%icondir == 0)
588 
589  ! sort based on coordinates
590  newtooldidx = (/(i, i=1, size(this%idxToGlobal))/)
591  call quicksortgrid(newtooldidx, size(newtooldidx), this%idxToGlobal, z_only)
592 
593  ! and invert
594  allocate (oldtonewidx(size(newtooldidx)))
595  do i = 1, size(oldtonewidx)
596  oldtonewidx(newtooldidx(i)) = i
597  end do
598 
599  ! reorder global table
600  allocate (sortedglobalmap(size(this%idxToGlobal)))
601  do i = 1, size(newtooldidx)
602  sortedglobalmap(i) = this%idxToGlobal(newtooldidx(i))
603  end do
604  do i = 1, size(newtooldidx)
605  this%idxToGlobal(i) = sortedglobalmap(i)
606  end do
607  deallocate (sortedglobalmap)
608 
609  ! reorder regional lookup table
610  allocate (sortedregionmap(size(this%region_to_iface_map)))
611  do i = 1, size(sortedregionmap)
612  if (this%region_to_iface_map(i) /= -1) then
613  idxold = this%region_to_iface_map(i)
614  sortedregionmap(i) = oldtonewidx(idxold)
615  else
616  sortedregionmap(i) = -1
617  end if
618  end do
619  do i = 1, size(sortedregionmap)
620  this%region_to_iface_map(i) = sortedregionmap(i)
621  end do
622  deallocate (sortedregionmap)
623 
subroutine, public quicksortgrid(array, arraySize, idxToGlobal, z_only)
Definition: GridSorting.f90:15
Here is the call graph for this function:

Variable Documentation

◆ initnrneighbors

integer(i4b), parameter gridconnectionmodule::initnrneighbors = 7
private

Definition at line 31 of file GridConnection.f90.

31  integer(I4B), parameter :: InitNrNeighbors = 7