From 469bc4b225de5fcbc9a453174a7e60cfd61c41e2 Mon Sep 17 00:00:00 2001 From: mjreno Date: Wed, 22 Jun 2022 12:26:57 -0400 Subject: [PATCH] style(connections): apply fprettify formatting (#960) * working toward consistent code formatting --- src/Model/Connection/CellWithNbrs.f90 | 31 +- src/Model/Connection/ConnectionBuilder.f90 | 106 +- src/Model/Connection/CsrUtils.f90 | 16 +- src/Model/Connection/GridConnection.f90 | 506 +++++----- src/Model/Connection/GridSorting.f90 | 126 +-- src/Model/Connection/GwfGwfConnection.f90 | 252 ++--- src/Model/Connection/GwfInterfaceModel.f90 | 83 +- src/Model/Connection/GwtGwtConnection.f90 | 928 +++++++++--------- src/Model/Connection/GwtInterfaceModel.f90 | 369 ++++--- .../Connection/SpatialModelConnection.f90 | 267 ++--- 10 files changed, 1359 insertions(+), 1325 deletions(-) diff --git a/src/Model/Connection/CellWithNbrs.f90 b/src/Model/Connection/CellWithNbrs.f90 index 20186bfa4c9..f9b4212c3b4 100644 --- a/src/Model/Connection/CellWithNbrs.f90 +++ b/src/Model/Connection/CellWithNbrs.f90 @@ -7,23 +7,24 @@ module CellWithNbrsModule integer(I4B), parameter :: defaultCapacity = 6 !> Data structure to hold a global cell identifier, - !! using a pointer to the model and its local cell + !! using a pointer to the model and its local cell !< index type, public :: GlobalCellType integer(I4B) :: index !< the index on the model grid class(NumericalModelType), pointer :: model => null() !< the model end type - + ! a global cell with neighbors type, public :: CellWithNbrsType type(GlobalCellType) :: cell integer(I4B) :: nrOfNbrs = 0 - type(CellWithNbrsType), dimension(:), pointer, contiguous :: neighbors => null() + type(CellWithNbrsType), dimension(:), pointer, & + contiguous :: neighbors => null() contains procedure :: addNbrCell end type - contains +contains subroutine addNbrCell(this, index, modelToAdd) class(CellWithNbrsType) :: this @@ -31,35 +32,35 @@ subroutine addNbrCell(this, index, modelToAdd) class(NumericalModelType), pointer :: modelToAdd ! local integer(I4B) :: nbrCnt, currentSize, i - type(CellWithNbrsType), dimension(:), pointer, contiguous :: newNeighbors - type(CellWithNbrsType), dimension(:), pointer, contiguous :: oldNeighbors + type(CellWithNbrsType), dimension(:), pointer, contiguous :: newNeighbors + type(CellWithNbrsType), dimension(:), pointer, contiguous :: oldNeighbors if (.not. associated(this%neighbors)) then - allocate(this%neighbors(defaultCapacity)) + allocate (this%neighbors(defaultCapacity)) this%nrOfNbrs = 0 end if - + nbrCnt = this%nrOfNbrs currentSize = size(this%neighbors) if (nbrCnt + 1 > currentSize) then ! inflate oldNeighbors => this%neighbors - allocate(newNeighbors(currentSize + defaultCapacity)) - do i=1, currentSize + allocate (newNeighbors(currentSize + defaultCapacity)) + do i = 1, currentSize newNeighbors(i) = oldNeighbors(i) end do this%neighbors => newNeighbors ! clean up - deallocate(oldNeighbors) - nullify(oldNeighbors) + deallocate (oldNeighbors) + nullify (oldNeighbors) end if - + this%neighbors(nbrCnt + 1)%cell%index = index this%neighbors(nbrCnt + 1)%cell%model => modelToAdd this%nrOfNbrs = nbrCnt + 1 - + end subroutine addNbrCell -end module \ No newline at end of file +end module diff --git a/src/Model/Connection/ConnectionBuilder.f90 b/src/Model/Connection/ConnectionBuilder.f90 index 66f5e2f08ba..fbbb972121d 100644 --- a/src/Model/Connection/ConnectionBuilder.f90 +++ b/src/Model/Connection/ConnectionBuilder.f90 @@ -6,31 +6,31 @@ module ConnectionBuilderModule use BaseSolutionModule, only: BaseSolutionType use NumericalSolutionModule, only: NumericalSolutionType use BaseExchangeModule, only: BaseExchangeType, GetBaseExchangeFromList - use DisConnExchangeModule, only: DisConnExchangeType, & - GetDisConnExchangeFromList + use DisConnExchangeModule, only: DisConnExchangeType, & + GetDisConnExchangeFromList use NumericalModelModule, only: NumericalModelType use SpatialModelConnectionModule, only: SpatialModelConnectionType, & - CastAsSpatialModelConnectionClass, & - GetSpatialModelConnectionFromList, & - AddSpatialModelConnectionToList - - implicit none + CastAsSpatialModelConnectionClass, & + GetSpatialModelConnectionFromList, & + AddSpatialModelConnectionToList + + implicit none private - + type, public :: ConnectionBuilderType contains procedure, pass(this) :: processSolution procedure, private, pass(this) :: processExchanges procedure, private, pass(this) :: setConnectionsToSolution procedure, private, pass(this) :: assignExchangesToConnections - end type ConnectionBuilderType - - contains + end type ConnectionBuilderType + +contains !> @brief Process the exchanges in the solution into model connections !! !! This routine processes all exchanges in a solution and, - !! when required, creates model connections of the proper + !! when required, creates model connections of the proper !! type (GWF-GWF, GWT-GWT, ...) for a subset. It removes this !! subset of exchanges from the solution and replaces them with the !! created connections. @@ -54,10 +54,10 @@ subroutine processSolution(this, solution) call this%processExchanges(numSol%exchangelist, newConnections) if (newConnections%Count() == 0) then return - end if + end if - write(iout,'(1x,a,i0,a,a)') 'Created ', newConnections%Count(), & - ' model connections for solution ', trim(solution%name) + write (iout, '(1x,a,i0,a,a)') 'Created ', newConnections%Count(), & + ' model connections for solution ', trim(solution%name) ! set the global exchanges from this solution to ! the model connections @@ -73,7 +73,7 @@ end subroutine processSolution !> @brief Create connections from exchanges !! - !! If the configuration demands it, this will create connections, + !! If the configuration demands it, this will create connections, !! for the exchanges (one connection per exchange) add them to !! the global list, and return them as @param newConnections !< @@ -96,10 +96,11 @@ subroutine processExchanges(this, exchanges, newConnections) ! Force use of the interface model dev_always_ifmod = .false. if (IDEVELOPMODE == 1) then - call get_environment_variable('DEV_ALWAYS_USE_IFMOD', value=envvar, status=status) + call get_environment_variable('DEV_ALWAYS_USE_IFMOD', & + value=envvar, status=status) if (status == 0 .and. envvar == '1') then dev_always_ifmod = .true. - write(*,'(a,/)') "### Experimental: forcing interface model ###" + write (*, '(a,/)') "### Experimental: forcing interface model ###" end if end if @@ -109,17 +110,17 @@ subroutine processExchanges(this, exchanges, newConnections) ! if it is not DisConnExchangeType, we can skip it continue end if - + ! for now, if we have XT3D on the interface, we use a connection, - ! (this will be more generic in the future) - if (conEx%use_interface_model() .or. conEx%dev_ifmod_on & + ! (this will be more generic in the future) + if (conEx%use_interface_model() .or. conEx%dev_ifmod_on & .or. dev_always_ifmod) then ! we should not get period connections here isPeriodic = associated(conEx%model1, conEx%model2) if (isPeriodic) then - write(*,*) 'Error (which should never happen): interface model '// & - 'does not support periodic boundary condition' + write (*, *) 'Error (which should never happen): interface model '// & + 'does not support periodic boundary condition' call ustop() end if @@ -128,7 +129,7 @@ subroutine processExchanges(this, exchanges, newConnections) call AddSpatialModelConnectionToList(baseconnectionlist, modelConnection) call AddSpatialModelConnectionToList(newConnections, modelConnection) - ! and for model 2, unless periodic + ! and for model 2, unless periodic modelConnection => createModelConnection(conEx%model2, conEx) call AddSpatialModelConnectionToList(baseconnectionlist, modelConnection) call AddSpatialModelConnectionToList(newConnections, modelConnection) @@ -142,7 +143,7 @@ subroutine processExchanges(this, exchanges, newConnections) exit end if end do - + end if end do @@ -158,40 +159,40 @@ function createModelConnection(model, exchange) result(connection) use GwfGwfConnectionModule, only: GwfGwfConnectionType use GwtGwtConnectionModule, only: GwtGwtConnectionType use GwfModule, only: GwfModelType - - class(NumericalModelType), pointer , intent(in) :: model !< the model for which the connection will be created + + class(NumericalModelType), pointer, intent(in) :: model !< the model for which the connection will be created class(DisConnExchangeType), pointer, intent(in) :: exchange !< the type of connection class(SpatialModelConnectionType), pointer :: connection !< the created connection - + ! different concrete connection types: class(GwfGwfConnectionType), pointer :: flowConnection => null() class(GwtGwtConnectionType), pointer :: transportConnection => null() - + connection => null() - + ! select on type of connection to create - select case(exchange%typename) - case('GWF-GWF') - allocate(GwfGwfConnectionType :: flowConnection) - call flowConnection%construct(model, exchange) - connection => flowConnection - flowConnection => null() - case('GWT-GWT') - allocate(GwtGwtConnectionType :: transportConnection) - call transportConnection%construct(model, exchange) - connection => transportConnection - transportConnection => null() - case default - write(*,*) 'Error (which should never happen): undefined exchangetype found' - call ustop() - end select - + select case (exchange%typename) + case ('GWF-GWF') + allocate (GwfGwfConnectionType :: flowConnection) + call flowConnection%construct(model, exchange) + connection => flowConnection + flowConnection => null() + case ('GWT-GWT') + allocate (GwtGwtConnectionType :: transportConnection) + call transportConnection%construct(model, exchange) + connection => transportConnection + transportConnection => null() + case default + write (*, *) 'Error (which should never happen): '// & + 'undefined exchangetype found' + call ustop() + end select + end function createModelConnection - - + !> @brief Set connections to the solution !! - !! This adds the connections to the solution and removes + !! This adds the connections to the solution and removes !! those exchanges which are replaced by a connection !< subroutine setConnectionsToSolution(this, connections, solution) @@ -211,7 +212,7 @@ subroutine setConnectionsToSolution(this, connections, solution) ! will this exchange be replaced by a connection? keepExchange = .true. do iconn = 1, connections%Count() - conn => GetSpatialModelConnectionFromList(connections,iconn) + conn => GetSpatialModelConnectionFromList(connections, iconn) exPtr2 => conn%primaryExchange if (associated(exPtr2, exPtr)) then ! if so, don't add it to the list @@ -284,8 +285,7 @@ subroutine assignExchangesToConnections(this, exchanges, connections) ! clean call keepList%Clear(destroy=.false.) - + end subroutine assignExchangesToConnections - - + end module ConnectionBuilderModule diff --git a/src/Model/Connection/CsrUtils.f90 b/src/Model/Connection/CsrUtils.f90 index 0bd0847cce5..729e4926c08 100644 --- a/src/Model/Connection/CsrUtils.f90 +++ b/src/Model/Connection/CsrUtils.f90 @@ -2,11 +2,11 @@ module CsrUtilsModule implicit none private - + public :: getCSRIndex - + contains - + !> @brief Return index for element i,j in CSR storage, !< returns -1 when not there function getCSRIndex(i, j, ia, ja) result(csrIndex) @@ -18,15 +18,15 @@ function getCSRIndex(i, j, ia, ja) result(csrIndex) integer(I4B) :: csrIndex !< the CSR ndex of element i,j ! local integer(I4B) :: idx - + csrIndex = -1 - do idx = ia(i), ia(i+1)-1 + do idx = ia(i), ia(i + 1) - 1 if (ja(idx) == j) then csrIndex = idx return end if end do - - end function -end module \ No newline at end of file + end function + +end module diff --git a/src/Model/Connection/GridConnection.f90 b/src/Model/Connection/GridConnection.f90 index 70a915783d0..eb37d36e8c9 100644 --- a/src/Model/Connection/GridConnection.f90 +++ b/src/Model/Connection/GridConnection.f90 @@ -9,17 +9,17 @@ module GridConnectionModule use GwfDisuModule use DisConnExchangeModule use CellWithNbrsModule - use ConnectionsModule + use ConnectionsModule use SparseModule, only: sparsematrix implicit none private - + ! Initial nr of neighbors for sparse matrix allocation integer(I4B), parameter :: InitNrNeighbors = 7 - + !> This class is used to construct the connections object for !! the interface model's spatial discretization/grid. - !! + !! !! It works as follows: !! !! 1: construct basic instance, allocate data structures @@ -48,7 +48,7 @@ module GridConnectionModule integer(I4B) :: exchangeStencilDepth !< stencil size at the interface class(NumericalModelType), pointer :: model => null() !< the model for which this grid connection exists - + integer(I4B), pointer :: nrOfBoundaryCells => null() !< nr of boundary cells with connection to another model type(CellWithNbrsType), dimension(:), pointer :: boundaryCells => null() !< cells on our side of the primary connections type(CellWithNbrsType), dimension(:), pointer :: connectedCells => null() !< cells on the neighbors side of the primary connection @@ -56,20 +56,19 @@ module GridConnectionModule !! the required depth integer, dimension(:), pointer :: primConnections => null() !< table mapping the index in the boundaryCells/connectedCells - !< arrays into a connection index for e.g. access to flowja - + integer(I4B), pointer :: nrOfCells => null() !< the total number of cells in the interface type(GlobalCellType), dimension(:), pointer :: idxToGlobal => null() !< a map from interface index to global coordinate integer(I4B), dimension(:), pointer, contiguous :: idxToGlobalIdx => null() !< a (flat) map from interface index to global index, !! stored in mem. mgr. so can be used for debugging - + integer(I4B), dimension(:), pointer :: regionalToInterfaceIdxMap => null() !< (sparse) mapping from regional index to interface ixd type(ListType) :: regionalModels !< the models participating in the interface integer(I4B), dimension(:), pointer :: regionalModelOffset => null() !< the new offset to compactify the range of indices integer(I4B), pointer :: indexCount => null() !< counts the number of cells in the interface type(ConnectionsType), pointer :: connections => null() !< sparse matrix with the connections integer(I4B), dimension(:), pointer :: connectionMask => null() !< to mask out connections from the amat coefficient calculation - + contains ! public procedure, pass(this) :: construct @@ -79,13 +78,13 @@ module GridConnectionModule procedure, pass(this) :: findModelNeighbors procedure, pass(this) :: extendConnection generic :: getInterfaceIndex => getInterfaceIndexByCell, & - getInterfaceIndexByIndexModel + getInterfaceIndexByIndexModel procedure, pass(this) :: getDiscretization - + ! 'protected' procedure, pass(this) :: isPeriodic - + ! private routines procedure, private, pass(this) :: buildConnections procedure, private, pass(this) :: addNeighbors @@ -109,8 +108,8 @@ module GridConnectionModule procedure, private, pass(this) :: setMaskOnConnection procedure, private, pass(this) :: createLookupTable end type - - contains + +contains !> @brief Construct the GridConnection and allocate !! the data structures for the primary connections @@ -126,21 +125,21 @@ subroutine construct(this, model, nrOfPrimaries, connectionName) this%memoryPath = create_mem_path(connectionName, 'GC') call this%allocateScalars() - - allocate(this%boundaryCells(nrOfPrimaries)) - allocate(this%connectedCells(nrOfPrimaries)) - allocate(this%primConnections(nrOfPrimaries)) - allocate(this%idxToGlobal(2*nrOfPrimaries)) - + + allocate (this%boundaryCells(nrOfPrimaries)) + allocate (this%connectedCells(nrOfPrimaries)) + allocate (this%primConnections(nrOfPrimaries)) + allocate (this%idxToGlobal(2 * nrOfPrimaries)) + call this%addToRegionalModels(model) - + this%nrOfBoundaryCells = 0 this%internalStencilDepth = 1 this%exchangeStencilDepth = 1 end subroutine construct - + !> @brief Connect neighboring cells at the interface by !! storing them in the boundary cell and connected cell !! arrays @@ -151,13 +150,14 @@ subroutine connectCell(this, idx1, model1, idx2, model2) class(NumericalModelType), pointer :: model1 !< model of cell 1 integer(I4B) :: idx2 !< local index cell 2 class(NumericalModelType), pointer :: model2 !< model of cell 2 - - this%nrOfBoundaryCells = this%nrOfBoundaryCells + 1 + + this%nrOfBoundaryCells = this%nrOfBoundaryCells + 1 if (this%nrOfBoundaryCells > size(this%boundaryCells)) then - write(*,*) 'Error: nr of cell connections exceeds capacity in grid connection, terminating...' + write (*, *) 'Error: nr of cell connections exceeds '// & + 'capacity in grid connection, terminating...' call ustop() end if - + if (associated(model1, this%model)) then this%boundaryCells(this%nrOfBoundaryCells)%cell%index = idx1 this%boundaryCells(this%nrOfBoundaryCells)%cell%model => this%model @@ -171,12 +171,11 @@ subroutine connectCell(this, idx1, model1, idx2, model2) this%connectedCells(this%nrOfBoundaryCells)%cell%index = idx1 this%connectedCells(this%nrOfBoundaryCells)%cell%model => model1 else - write(*,*) 'Error: unable to connect cells outside the model' + write (*, *) 'Error: unable to connect cells outside the model' call ustop() end if - - end subroutine connectCell + end subroutine connectCell !> @brief Create the tree structure with all model nbrs, nbrs-of-nbrs, !< etc. for this model up to the specified depth @@ -190,10 +189,11 @@ subroutine findModelNeighbors(this, globalExchanges, depth) end subroutine findModelNeighbors - !> @brief Add neighbors and nbrs-of-nbrs to the model tree - !< - recursive subroutine addModelNeighbors(this, model, globalExchanges, depth, mask) + !< + recursive subroutine addModelNeighbors(this, model, & + globalExchanges, & + depth, mask) class(GridConnectionType), intent(inout) :: this !< this grid connection class(NumericalModelType), pointer, intent(inout) :: model !< the model to add neighbors for type(ListType), intent(inout) :: globalExchanges !< list with all exchanges @@ -203,7 +203,7 @@ recursive subroutine addModelNeighbors(this, model, globalExchanges, depth, mask integer(I4B) :: i, n class(DisConnExchangeType), pointer :: connEx class(NumericalModelType), pointer :: neighborModel - class(NumericalModelType), pointer :: modelMask + class(NumericalModelType), pointer :: modelMask type(ListType) :: nbrModels class(*), pointer :: objPtr procedure(isEqualIface), pointer :: areEqualMethod @@ -214,7 +214,7 @@ recursive subroutine addModelNeighbors(this, model, globalExchanges, depth, mask modelMask => mask end if - ! first find all direct neighbors of the model and add them, + ! first find all direct neighbors of the model and add them, ! avoiding duplicates do i = 1, globalExchanges%Count() @@ -228,7 +228,7 @@ recursive subroutine addModelNeighbors(this, model, globalExchanges, depth, mask ! check if there is a neighbor, and it is not masked ! (to prevent back-and-forth connections) - if (associated(neighborModel) .and. .not. & + if (associated(neighborModel) .and. .not. & associated(neighborModel, modelMask)) then ! add to neighbors @@ -249,11 +249,11 @@ recursive subroutine addModelNeighbors(this, model, globalExchanges, depth, mask end if end do - + ! now recurse on the neighbors up to the specified depth depth = depth - 1 if (depth == 0) return - + do n = 1, nbrModels%Count() neighborModel => GetNumericalModelFromList(nbrModels, n) call this%addModelNeighbors(neighborModel, globalExchanges, depth, model) @@ -264,14 +264,13 @@ recursive subroutine addModelNeighbors(this, model, globalExchanges, depth, mask end subroutine addModelNeighbors - !> @brief Add a model to a list of all regional models !< subroutine addToRegionalModels(this, modelToAdd) class(GridConnectionType), intent(inout) :: this !< this grid connection class(NumericalModelType), pointer :: modelToAdd !< the model to add to the region ! local - class(*), pointer :: mPtr + class(*), pointer :: mPtr procedure(isEqualIface), pointer :: areEqualMethod mPtr => modelToAdd @@ -279,10 +278,10 @@ subroutine addToRegionalModels(this, modelToAdd) if (.not. this%regionalModels%ContainsObject(mPtr, areEqualMethod)) then call AddNumericalModelToList(this%regionalModels, modelToAdd) end if - + end subroutine addToRegionalModels - !> @brief Extend the connection topology to deal with + !> @brief Extend the connection topology to deal with !! higher levels of connectivity (neighbors-of-neighbors, etc.) !! !! The following steps are taken: @@ -291,35 +290,35 @@ end subroutine addToRegionalModels !! 3. Allocate a (sparse) mapping table for the region !! 4. Build connection object for the interface grid, and the mask !< - subroutine extendConnection(this) + subroutine extendConnection(this) class(GridConnectionType), intent(inout) :: this !< this grid connection - ! local + ! local integer(I4B) :: remoteDepth, localDepth integer(I4B) :: icell integer(I4B) :: imod, regionSize, offset class(NumericalModelType), pointer :: numModel - + ! we need (stencildepth-1) extra cells for the interior remoteDepth = this%exchangeStencilDepth - localDepth = 2*this%internalStencilDepth - 1 + localDepth = 2 * this%internalStencilDepth - 1 if (localDepth < remoteDepth) then localDepth = remoteDepth end if - - ! first add the neighbors for the interior + + ! first add the neighbors for the interior ! (possibly extending into other models) do icell = 1, this%nrOfBoundaryCells - call this%addNeighbors(this%boundaryCells(icell), localDepth, & + call this%addNeighbors(this%boundaryCells(icell), localDepth, & this%connectedCells(icell)%cell, .true.) end do ! and for the exterior do icell = 1, this%nrOfBoundaryCells - call this%addNeighbors(this%connectedCells(icell), remoteDepth, & + call this%addNeighbors(this%connectedCells(icell), remoteDepth, & this%boundaryCells(icell)%cell, .false.) end do - + ! set up mapping for the region (models participating in interface model grid) - allocate(this%regionalModelOffset(this%regionalModels%Count())) + allocate (this%regionalModelOffset(this%regionalModels%Count())) regionSize = 0 offset = 0 do imod = 1, this%regionalModels%Count() @@ -329,11 +328,11 @@ subroutine extendConnection(this) offset = offset + numModel%dis%nodes end do ! init to -1, meaning 'interface index was not assigned yet' - allocate(this%regionalToInterfaceIdxMap(regionSize)) + allocate (this%regionalToInterfaceIdxMap(regionSize)) this%regionalToInterfaceIdxMap = -1 - + call this%buildConnections() - + end subroutine extendConnection !> @brief Builds a sparse matrix holding all cell connections, @@ -343,12 +342,12 @@ subroutine buildConnections(this) ! local integer(I4B) :: icell, iconn integer(I4B), dimension(:), allocatable :: nnz - type(SparseMatrix), pointer :: sparse - integer(I4B) :: ierror + type(SparseMatrix), pointer :: sparse + integer(I4B) :: ierror type(ConnectionsType), pointer :: conn - + ! Recursively generate interface cell indices, fill map to global cells, - ! and add to region lookup table + ! and add to region lookup table this%indexCount = 0 do icell = 1, this%nrOfBoundaryCells call this%registerInterfaceCells(this%boundaryCells(icell)) @@ -363,20 +362,20 @@ subroutine buildConnections(this) ! sort interface indexes such that 'n > m' means 'n below m' call this%sortInterfaceGrid() - + ! allocate a map from interface index to global coordinates - call mem_allocate(this%idxToGlobalIdx, this%nrOfCells, & + call mem_allocate(this%idxToGlobalIdx, this%nrOfCells, & 'IDXTOGLOBALIDX', this%memoryPath) - + ! create sparse data structure, to temporarily hold connections - allocate(sparse) - allocate(nnz(this%nrOfCells)) - nnz = InitNrNeighbors+1 + allocate (sparse) + allocate (nnz(this%nrOfCells)) + nnz = InitNrNeighbors + 1 call sparse%init(this%nrOfCells, this%nrOfCells, nnz) - - ! now (recursively) add connections to sparse, start with + + ! now (recursively) add connections to sparse, start with ! the primary connections (n-m from the exchange files) - call this%makePrimaryConnections(sparse) + call this%makePrimaryConnections(sparse) ! then into own domain do icell = 1, this%nrOfBoundaryCells call this%connectNeighborCells(this%boundaryCells(icell), sparse) @@ -385,93 +384,94 @@ subroutine buildConnections(this) do icell = 1, this%nrOfBoundaryCells call this%connectNeighborCells(this%connectedCells(icell), sparse) end do - - ! create connections object - allocate(this%connections) - conn => this%connections + + ! create connections object + allocate (this%connections) + conn => this%connections call conn%allocate_scalars(this%memoryPath) conn%nodes = this%nrOfCells conn%nja = sparse%nnz - conn%njas = (conn%nja - conn%nodes) / 2 + conn%njas = (conn%nja - conn%nodes) / 2 call conn%allocate_arrays() do iconn = 1, conn%njas conn%anglex(iconn) = -999. end do ! fill connection from sparse - call sparse%filliaja(conn%ia, conn%ja, ierror) + call sparse%filliaja(conn%ia, conn%ja, ierror) if (ierror /= 0) then - write(*,*) 'Error filling ia/ja in GridConnection: terminating...' + write (*, *) 'Error filling ia/ja in GridConnection: terminating...' call ustop() - end if + end if call fillisym(conn%nodes, conn%nja, conn%ia, conn%ja, conn%isym) - call filljas(conn%nodes, conn%nja, conn%ia, conn%ja, conn%isym, conn%jas) + call filljas(conn%nodes, conn%nja, conn%ia, conn%ja, conn%isym, conn%jas) call sparse%destroy() - - ! fill connection data (ihc, cl1, cl2, etc.) using data + + ! fill connection data (ihc, cl1, cl2, etc.) using data ! from models and exchanges - call this%fillConnectionDataInternal() + call this%fillConnectionDataInternal() call this%fillConnectionDataFromExchanges() - + ! set the masks on connections call this%createConnectionMask() ! create lookup table(s) call this%createLookupTable() - - end subroutine buildConnections - + end subroutine buildConnections + !< @brief Routine for finding neighbors-of-neighbors, recursively !< recursive subroutine addNeighbors(this, cellNbrs, depth, mask, interior) use SimModule, only: ustop class(GridConnectionType), intent(inout) :: this !< this grid connection - type(CellWithNbrsType), intent(inout) :: cellNbrs !< cell to add to + type(CellWithNbrsType), intent(inout) :: cellNbrs !< cell to add to integer(I4B), intent(inout) :: depth !< current depth (typically decreases in recursion) type(GlobalCellType), optional :: mask !< mask to excluded back-and-forth connection between cells logical(LGP) :: interior !< when true, we are adding from the exchange back into the model ! local integer(I4B) :: nbrIdx, ipos, inbr type(ConnectionsType), pointer :: conn - integer(I4B) :: newDepth - + integer(I4B) :: newDepth + ! if depth == 1, then we are not adding neighbors but use ! the boundary and connected cell only if (depth < 2) then return end if newDepth = depth - 1 - + conn => cellNbrs%cell%model%dis%con - + ! find neighbors local to this cell by looping through grid connections - do ipos=conn%ia(cellNbrs%cell%index) + 1, conn%ia(cellNbrs%cell%index+1) - 1 + do ipos = conn%ia(cellNbrs%cell%index) + 1, & + conn%ia(cellNbrs%cell%index + 1) - 1 nbrIdx = conn%ja(ipos) call this%addNeighborCell(cellNbrs, nbrIdx, cellNbrs%cell%model, mask) end do - + ! add remote nbr using the data from the exchanges call this%addRemoteNeighbors(cellNbrs, mask) - + ! now find nbr-of-nbr - do inbr=1, cellNbrs%nrOfNbrs + do inbr = 1, cellNbrs%nrOfNbrs ! are we leaving the model through another exchange? if (interior .and. associated(cellNbrs%cell%model, this%model)) then - if (.not. associated(cellNbrs%neighbors(inbr)%cell%model, this%model)) then + if (.not. associated(cellNbrs%neighbors(inbr)%cell%model, & + this%model)) then ! decrement by 1, because the connection we are crossing is not ! calculated by this interface newDepth = newDepth - 1 - end if + end if end if ! and add neigbors with the new depth - call this%addNeighbors(cellNbrs%neighbors(inbr), newDepth, & + call this%addNeighbors(cellNbrs%neighbors(inbr), newDepth, & cellNbrs%cell, interior) end do - + end subroutine addNeighbors - + !> @brief Add cell neighbors across models using the stored exchange !! data structures subroutine addRemoteNeighbors(this, cellNbrs, mask) @@ -481,9 +481,9 @@ subroutine addRemoteNeighbors(this, cellNbrs, mask) ! local integer(I4B) :: ix, iexg type(DisConnExchangeType), pointer :: connEx - - ! loop over all exchanges - do ix = 1, this%exchanges%Count() + + ! loop over all exchanges + do ix = 1, this%exchanges%Count() connEx => GetDisConnExchangeFromList(this%exchanges, ix) ! loop over n-m links in the exchange @@ -491,7 +491,7 @@ subroutine addRemoteNeighbors(this, cellNbrs, mask) do iexg = 1, connEx%nexg if (connEx%nodem1(iexg) == cellNbrs%cell%index) then ! we have a link, now add foreign neighbor - call this%addNeighborCell(cellNbrs, connEx%nodem2(iexg), & + call this%addNeighborCell(cellNbrs, connEx%nodem2(iexg), & connEx%model2, mask) end if end do @@ -501,27 +501,26 @@ subroutine addRemoteNeighbors(this, cellNbrs, mask) do iexg = 1, connEx%nexg if (connEx%nodem2(iexg) == cellNbrs%cell%index) then ! we have a link, now add foreign neighbor - call this%addNeighborCell(cellNbrs, connEx%nodem1(iexg), & + call this%addNeighborCell(cellNbrs, connEx%nodem1(iexg), & connEx%model1, mask) end if end do end if - + end do - + end subroutine addRemoteNeighbors - !> @brief Add neighboring cell to tree structure !< - subroutine addNeighborCell(this, cellNbrs, newNbrIdx, nbrModel, mask) + subroutine addNeighborCell(this, cellNbrs, newNbrIdx, nbrModel, mask) class(GridConnectionType), intent(in) :: this !< this grid connection instance type(CellWithNbrsType), intent(inout) :: cellNbrs !< the root cell which to add to integer(I4B), intent(in) :: newNbrIdx !< the neigboring cell's index class(NumericalModelType), pointer :: nbrModel !< the model where the new neighbor lives - type(GlobalCellType), optional :: mask !< don't add connections to this cell (optional) + type(GlobalCellType), optional :: mask !< don't add connections to this cell (optional) ! local - + if (present(mask)) then if (newNbrIdx == mask%index .and. associated(nbrModel, mask%model)) then return @@ -530,7 +529,7 @@ subroutine addNeighborCell(this, cellNbrs, newNbrIdx, nbrModel, mask) call cellNbrs%addNbrCell(newNbrIdx, nbrModel) end subroutine addNeighborCell - + !> @brief Recursively set interface cell indexes and !< add to the region-to-interface loopup table recursive subroutine registerInterfaceCells(this, cellWithNbrs) @@ -540,7 +539,7 @@ recursive subroutine registerInterfaceCells(this, cellWithNbrs) integer(I4B) :: offset, inbr integer(I4B) :: regionIdx ! unique idx in the region (all connected models) integer(I4B) :: ifaceIdx ! unique idx in the interface grid - + offset = this%getRegionalModelOffset(cellWithNbrs%cell%model) regionIdx = offset + cellWithNbrs%cell%index ifaceIdx = this%getInterfaceIndex(cellWithNbrs%cell) @@ -550,12 +549,12 @@ recursive subroutine registerInterfaceCells(this, cellWithNbrs) call this%addToGlobalMap(ifaceIdx, cellWithNbrs%cell) this%regionalToInterfaceIdxMap(regionIdx) = ifaceIdx end if - + ! and also for its neighbors do inbr = 1, cellWithNbrs%nrOfNbrs call this%registerInterfaceCells(cellWithNbrs%neighbors(inbr)) end do - + end subroutine registerInterfaceCells !> @brief Add entry to lookup table, inflating when necessary @@ -571,13 +570,13 @@ subroutine addToGlobalMap(this, ifaceIdx, cell) ! inflate? currentSize = size(this%idxToGlobal) if (ifaceIdx > currentSize) then - newSize = nint(1.5*currentSize) - allocate(tempMap(newSize)) + newSize = nint(1.5 * currentSize) + allocate (tempMap(newSize)) do i = 1, currentSize tempMap(i) = this%idxToGlobal(i) end do - - deallocate(this%idxToGlobal) + + deallocate (this%idxToGlobal) this%idxToGlobal => tempMap end if @@ -593,12 +592,12 @@ subroutine compressGlobalMap(this) type(GlobalCellType), dimension(:), pointer :: tempMap if (size(this%idxToGlobal) > this%nrOfCells) then - allocate(tempMap(this%nrOfCells)) + allocate (tempMap(this%nrOfCells)) tempMap(1:this%nrOfCells) = this%idxToGlobal(1:this%nrOfCells) - deallocate(this%idxToGlobal) - allocate(this%idxToGlobal(this%nrOfCells)) + deallocate (this%idxToGlobal) + allocate (this%idxToGlobal(this%nrOfCells)) this%idxToGlobal(1:this%nrOfCells) = tempMap(1:this%nrOfCells) - deallocate(tempMap) + deallocate (tempMap) end if end subroutine compressGlobalMap @@ -617,28 +616,28 @@ subroutine sortInterfaceGrid(this) integer(I4B), dimension(:), allocatable :: sortedRegionMap ! sort based on coordinates - newToOldIdx = (/ (i, i=1, size(this%idxToGlobal)) /) + newToOldIdx = (/(i, i=1, size(this%idxToGlobal))/) call quickSortGrid(newToOldIdx, size(newToOldIdx), this%idxToGlobal) - + ! and invert - allocate(oldToNewIdx(size(newToOldIdx))) - do i=1, size(oldToNewIdx) + allocate (oldToNewIdx(size(newToOldIdx))) + do i = 1, size(oldToNewIdx) oldToNewIdx(newToOldIdx(i)) = i end do ! reorder global table - allocate(sortedGlobalMap(size(this%idxToGlobal))) - do i=1, size(newToOldIdx) + allocate (sortedGlobalMap(size(this%idxToGlobal))) + do i = 1, size(newToOldIdx) sortedGlobalMap(i) = this%idxToGlobal(newToOldIdx(i)) end do - do i=1, size(newToOldIdx) + do i = 1, size(newToOldIdx) this%idxToGlobal(i) = sortedGlobalMap(i) end do - deallocate(sortedGlobalMap) + deallocate (sortedGlobalMap) - ! reorder regional lookup table - allocate(sortedRegionMap(size(this%regionalToInterfaceIdxMap))) - do i=1, size(sortedRegionMap) + ! reorder regional lookup table + allocate (sortedRegionMap(size(this%regionalToInterfaceIdxMap))) + do i = 1, size(sortedRegionMap) if (this%regionalToInterfaceIdxMap(i) /= -1) then idxOld = this%regionalToInterfaceIdxMap(i) sortedRegionMap(i) = oldToNewIdx(idxOld) @@ -646,13 +645,13 @@ subroutine sortInterfaceGrid(this) sortedRegionMap(i) = -1 end if end do - do i=1, size(sortedRegionMap) + do i = 1, size(sortedRegionMap) this%regionalToInterfaceIdxMap(i) = sortedRegionMap(i) end do - deallocate(sortedRegionMap) - + deallocate (sortedRegionMap) + end subroutine sortInterfaceGrid - + !> @brief Add primary connections to the sparse data structure !< subroutine makePrimaryConnections(this, sparse) @@ -661,22 +660,22 @@ subroutine makePrimaryConnections(this, sparse) ! local integer(I4B) :: icell integer(I4B) :: ifaceIdx, ifaceIdxNbr - - do icell = 1, this%nrOfBoundaryCells + + do icell = 1, this%nrOfBoundaryCells ifaceIdx = this%getInterfaceIndex(this%boundaryCells(icell)%cell) ifaceIdxNbr = this%getInterfaceIndex(this%connectedCells(icell)%cell) - + ! add diagonals to sparse call sparse%addconnection(ifaceIdx, ifaceIdx, 1) call sparse%addconnection(ifaceIdxNbr, ifaceIdxNbr, 1) - + ! and cross terms call sparse%addconnection(ifaceIdx, ifaceIdxNbr, 1) call sparse%addconnection(ifaceIdxNbr, ifaceIdx, 1) end do - + end subroutine makePrimaryConnections - + !> @brief Recursively add higher order connections (from !! cells neighoring the primarily connected cells) to the !< sparse data structure @@ -687,38 +686,38 @@ recursive subroutine connectNeighborCells(this, cell, sparse) ! local integer(I4B) :: ifaceIdx, ifaceIdxNbr ! unique idx in the interface grid integer(I4B) :: inbr - + ifaceIdx = this%getInterfaceIndex(cell%cell) do inbr = 1, cell%nrOfNbrs ifaceIdxNbr = this%getInterfaceIndex(cell%neighbors(inbr)%cell) - + call sparse%addconnection(ifaceIdxNbr, ifaceIdxNbr, 1) call sparse%addconnection(ifaceIdx, ifaceIdxNbr, 1) call sparse%addconnection(ifaceIdxNbr, ifaceIdx, 1) - + ! recurse call this%connectNeighborCells(cell%neighbors(inbr), sparse) end do - + end subroutine connectNeighborCells - + !> @brief Fill connection data (ihc, cl1, ...) for !< connections between cells within the same model. subroutine fillConnectionDataInternal(this) use ConstantsModule, only: DPI, DTWOPI class(GridConnectionType), intent(inout) :: this !< this grid connection instance ! local - type(ConnectionsType), pointer :: conn, connOrig + type(ConnectionsType), pointer :: conn, connOrig integer(I4B) :: n, m, ipos, isym, iposOrig, isymOrig type(GlobalCellType), pointer :: ncell, mcell - + conn => this%connections - + do n = 1, conn%nodes - do ipos=conn%ia(n)+1, conn%ia(n+1)-1 + do ipos = conn%ia(n) + 1, conn%ia(n + 1) - 1 m = conn%ja(ipos) if (n > m) cycle - + isym = conn%jas(ipos) ncell => this%idxToGlobal(n) mcell => this%idxToGlobal(m) @@ -727,15 +726,15 @@ subroutine fillConnectionDataInternal(this) connOrig => ncell%model%dis%con iposOrig = connOrig%getjaindex(ncell%index, mcell%index) if (iposOrig == 0) then - ! periodic boundary conditions can add connections between cells in + ! periodic boundary conditions can add connections between cells in ! the same model, but they are dealt with through the exchange data if (this%isPeriodic(ncell%index, mcell%index)) cycle - + ! this should not be possible - write(*,*) 'Error: cannot find cell connection in model grid' - call ustop() + write (*, *) 'Error: cannot find cell connection in model grid' + call ustop() end if - + isymOrig = connOrig%jas(iposOrig) conn%hwva(isym) = connOrig%hwva(isymOrig) conn%ihc(isym) = connOrig%ihc(isymOrig) @@ -752,12 +751,12 @@ subroutine fillConnectionDataInternal(this) end do end do end subroutine fillConnectionDataInternal - - !> @brief Fill connection data (ihc, cl1, ...) for + + !> @brief Fill connection data (ihc, cl1, ...) for !< all exchanges subroutine fillConnectionDataFromExchanges(this) use ConstantsModule, only: DPI, DTWOPI, DPIO180 - use ArrayHandlersModule, only: ifind + use ArrayHandlersModule, only: ifind class(GridConnectionType), intent(inout) :: this !< this grid connection instance ! local integer(I4B) :: inx, iexg, ivalAngldegx @@ -765,12 +764,12 @@ subroutine fillConnectionDataFromExchanges(this) integer(I4B) :: nOffset, mOffset, nIfaceIdx, mIfaceIdx class(DisConnExchangeType), pointer :: connEx type(ConnectionsType), pointer :: conn - + conn => this%connections - + do inx = 1, this%exchanges%Count() - connEx => GetDisConnExchangeFromList(this%exchanges, inx) - + connEx => GetDisConnExchangeFromList(this%exchanges, inx) + ivalAngldegx = -1 if (connEx%naux > 0) then ivalAngldegx = ifind(connEx%auxname, 'ANGLDEGX') @@ -778,50 +777,51 @@ subroutine fillConnectionDataFromExchanges(this) conn%ianglex = 1 end if end if - + nOffset = this%getRegionalModelOffset(connEx%model1) mOffset = this%getRegionalModelOffset(connEx%model2) do iexg = 1, connEx%nexg nIfaceIdx = this%regionalToInterfaceIdxMap(noffset + connEx%nodem1(iexg)) mIfaceIdx = this%regionalToInterfaceIdxMap(moffset + connEx%nodem2(iexg)) - ! not all nodes from the exchanges are part of the interface grid + ! not all nodes from the exchanges are part of the interface grid ! (think of exchanges between neigboring models, and their neighbors) if (nIFaceIdx == -1 .or. mIFaceIdx == -1) then cycle end if - - ipos = conn%getjaindex(nIfaceIdx, mIfaceIdx) - ! (see prev. remark) sometimes the cells are in the interface grid, + + ipos = conn%getjaindex(nIfaceIdx, mIfaceIdx) + ! (see prev. remark) sometimes the cells are in the interface grid, ! but the connection isn't. This can happen for leaf nodes of the grid. if (ipos == 0) then ! no match, safely cycle cycle - end if + end if isym = conn%jas(ipos) - + ! note: cl1 equals L_nm: the length from cell n to the shared ! face with cell m (and cl2 analogously for L_mn) if (nIfaceIdx < mIfaceIdx) then conn%cl1(isym) = connEx%cl1(iexg) conn%cl2(isym) = connEx%cl2(iexg) if (ivalAngldegx > 0) then - conn%anglex(isym) = connEx%auxvar(ivalAngldegx,iexg) * DPIO180 + conn%anglex(isym) = connEx%auxvar(ivalAngldegx, iexg) * DPIO180 end if else conn%cl1(isym) = connEx%cl2(iexg) conn%cl2(isym) = connEx%cl1(iexg) if (ivalAngldegx > 0) then - conn%anglex(isym) = mod(connEx%auxvar(ivalAngldegx,iexg) + 180.0_DP, 360.0_DP) * DPIO180 + conn%anglex(isym) = mod(connEx%auxvar(ivalAngldegx, iexg) + & + 180.0_DP, 360.0_DP) * DPIO180 end if end if conn%hwva(isym) = connEx%hwva(iexg) conn%ihc(isym) = connEx%ihc(iexg) - - end do + + end do end do - - end subroutine fillConnectionDataFromExchanges - + + end subroutine fillConnectionDataFromExchanges + !> @brief Create the connection masks !! !! The level indicates the nr of connections away from @@ -834,55 +834,59 @@ subroutine createConnectionMask(this) integer(I4B) :: icell, inbr, n, ipos integer(I4B) :: level, newMask type(CellWithNbrsType), pointer :: cell, nbrCell - + ! set all masks to zero to begin with do ipos = 1, this%connections%nja - call this%connections%set_mask(ipos, 0) + call this%connections%set_mask(ipos, 0) end do - + ! remote connections remain masked ! now set mask for exchange connections (level == 1) level = 1 - do icell = 1, this%nrOfBoundaryCells - call this%setMaskOnConnection(this%boundaryCells(icell), this%connectedCells(icell), level) + do icell = 1, this%nrOfBoundaryCells + call this%setMaskOnConnection(this%boundaryCells(icell), & + this%connectedCells(icell), level) ! for cross-boundary connections, we need to apply the mask to both n-m and m-n, ! because if the upper triangular one is disabled, its transposed (lower triangular) ! counter part is skipped in the NPF calculation as well. - call this%setMaskOnConnection(this%connectedCells(icell), this%boundaryCells(icell), level) + call this%setMaskOnConnection(this%connectedCells(icell), & + this%boundaryCells(icell), level) end do - + ! now extend mask recursively into the internal domain (level > 1) do icell = 1, this%nrOfBoundaryCells - cell => this%boundaryCells(icell) + cell => this%boundaryCells(icell) do inbr = 1, cell%nrOfNbrs nbrCell => this%boundaryCells(icell)%neighbors(inbr) level = 2 ! this is incremented within the recursion - call this%maskInternalConnections(this%boundaryCells(icell), this%boundaryCells(icell)%neighbors(inbr), level) - end do + call this%maskInternalConnections(this%boundaryCells(icell), & + this%boundaryCells(icell)% & + neighbors(inbr), level) + end do end do - + ! set normalized mask: ! =1 for links with connectivity <= interior stencil depth ! =0 otherwise - do n = 1, this%connections%nodes + do n = 1, this%connections%nodes ! set diagonals to zero call this%connections%set_mask(this%connections%ia(n), 0) - + do ipos = this%connections%ia(n) + 1, this%connections%ia(n + 1) - 1 newMask = 0 if (this%connections%mask(ipos) > 0) then if (this%connections%mask(ipos) < this%internalStencilDepth + 1) then newMask = 1 end if - end if + end if ! set mask on off-diag call this%connections%set_mask(ipos, newMask) - end do + end do end do - + end subroutine createConnectionMask - !> @brief Create lookup tables for efficient access + !> @brief Create lookup tables for efficient access !< (this needs the connections object to be available) subroutine createLookupTable(this) use CsrUtilsModule, only: getCSRIndex @@ -893,15 +897,15 @@ subroutine createLookupTable(this) do i = 1, this%nrOfBoundaryCells n1 = this%getInterfaceIndexByIndexModel(this%boundaryCells(i)%cell%index, & this%boundaryCells(i)%cell%model) - n2 = this%getInterfaceIndexByIndexModel(this%connectedCells(i)%cell%index,& + n2 = this%getInterfaceIndexByIndexModel(this%connectedCells(i)%cell%index, & this%connectedCells(i)%cell%model) - + ipos = getCSRIndex(n1, n2, this%connections%ia, this%connections%ja) this%primConnections(i) = ipos - end do + end do end subroutine createLookupTable - + !> @brief Recursively mask connections, increasing the level as we go !< recursive subroutine maskInternalConnections(this, cell, nbrCell, level) @@ -911,24 +915,26 @@ recursive subroutine maskInternalConnections(this, cell, nbrCell, level) integer(I4B), intent(in) :: level ! local integer(I4B) :: inbr, newLevel - + ! only set the mask for internal connections, leaving the ! others at 0 - if (associated(cell%cell%model, this%model) .and. & + if (associated(cell%cell%model, this%model) .and. & associated(nbrCell%cell%model, this%model)) then ! this will set a mask on both diagonal, and both cross terms call this%setMaskOnConnection(cell, nbrCell, level) call this%setMaskOnConnection(nbrCell, cell, level) end if - + ! recurse on nbrs-of-nbrs newLevel = level + 1 - do inbr = 1, nbrCell%nrOfNbrs - call this%maskInternalConnections(nbrCell, nbrCell%neighbors(inbr), newLevel) + do inbr = 1, nbrCell%nrOfNbrs + call this%maskInternalConnections(nbrCell, & + nbrCell%neighbors(inbr), & + newLevel) end do - + end subroutine maskInternalConnections - + !> @brief Set a mask on the connection from a cell to its neighbor, !! (and not the transposed!) not overwriting the current level !< of a connection when it is smaller @@ -941,10 +947,10 @@ subroutine setMaskOnConnection(this, cell, nbrCell, level) integer(I4B) :: ifaceIdx, ifaceIdxNbr integer(I4B) :: iposdiag, ipos integer(I4B) :: currentLevel - + ifaceIdx = this%getInterfaceIndex(cell%cell) ifaceIdxNbr = this%getInterfaceIndex(nbrCell%cell) - + ! diagonal iposdiag = this%connections%getjaindex(ifaceIdx, ifaceIdx) currentLevel = this%connections%mask(iposdiag) @@ -957,9 +963,9 @@ subroutine setMaskOnConnection(this, cell, nbrCell, level) if (currentLevel == 0 .or. level < currentLevel) then call this%connections%set_mask(ipos, level) end if - + end subroutine setMaskOnConnection - + !> @brief Get interface index from global cell !< function getInterfaceIndexByCell(this, cell) result(ifaceIdx) @@ -968,7 +974,7 @@ function getInterfaceIndexByCell(this, cell) result(ifaceIdx) integer(I4B) :: ifaceIdx !< the index in the interface model ! local integer(I4B) :: offset, regionIdx - + offset = this%getRegionalModelOffset(cell%model) regionIdx = offset + cell%index ifaceIdx = this%regionalToInterfaceIdxMap(regionIdx) @@ -983,12 +989,12 @@ function getInterfaceIndexByIndexModel(this, index, model) result(ifaceIdx) integer(I4B) :: ifaceIdx !< the index in the interface model ! local integer(I4B) :: offset, regionIdx - + offset = this%getRegionalModelOffset(model) regionIdx = offset + index ifaceIdx = this%regionalToInterfaceIdxMap(regionIdx) end function getInterfaceIndexByIndexModel - + !> @brief Get the offset for a regional model !< function getRegionalModelOffset(this, model) result(offset) @@ -1000,40 +1006,40 @@ function getRegionalModelOffset(this, model) result(offset) class(NumericalModelType), pointer :: modelInList offset = 0 do im = 1, this%regionalModels%Count() - modelInList => GetNumericalModelFromList(this%regionalModels, im) - if (associated(model, modelInList)) then - offset = this%regionalModelOffset(im) - return - end if + modelInList => GetNumericalModelFromList(this%regionalModels, im) + if (associated(model, modelInList)) then + offset = this%regionalModelOffset(im) + return + end if end do - + end function getRegionalModelOffset - + !> @brief Allocate scalar data !< subroutine allocateScalars(this) use MemoryManagerModule, only: mem_allocate class(GridConnectionType) :: this !< this grid connection instance - + call mem_allocate(this%nrOfBoundaryCells, 'NRBNDCELLS', this%memoryPath) call mem_allocate(this%indexCount, 'IDXCOUNT', this%memoryPath) call mem_allocate(this%nrOfCells, 'NRCELLS', this%memoryPath) - + end subroutine allocateScalars !> @brief Sets the discretization (DISU) after all !! preprocessing by this grid connection has been done, !< this comes after disu_cr - subroutine getDiscretization(this, disu) - use ConnectionsModule + subroutine getDiscretization(this, disu) + use ConnectionsModule use SparseModule, only: sparsematrix class(GridConnectionType) :: this !< the grid connection - class(GwfDisuType), pointer :: disu !< the target disu object + class(GwfDisuType), pointer :: disu !< the target disu object ! local integer(I4B) :: icell, nrOfCells, idx type(NumericalModelType), pointer :: model real(DP) :: x, y, xglo, yglo - + ! the following is similar to dis_df nrOfCells = this%nrOfCells disu%nodes = nrOfCells @@ -1041,23 +1047,23 @@ subroutine getDiscretization(this, disu) disu%nja = this%connections%nja call disu%allocate_arrays() - ! these are otherwise allocated in dis%read_dimensions + ! these are otherwise allocated in dis%read_dimensions call disu%allocate_arrays_mem() - + ! fill data do icell = 1, nrOfCells idx = this%idxToGlobal(icell)%index model => this%idxToGlobal(icell)%model - + disu%top(icell) = model%dis%top(idx) disu%bot(icell) = model%dis%bot(idx) disu%area(icell) = model%dis%area(idx) end do - + ! grid connections follow from GridConnection: disu%con => this%connections - disu%njas = disu%con%njas - + disu%njas = disu%con%njas + ! copy cell x,y do icell = 1, nrOfCells idx = this%idxToGlobal(icell)%index @@ -1067,8 +1073,8 @@ subroutine getDiscretization(this, disu) ! we are merging grids with possibly (likely) different origins, ! transform: call model%dis%transform_xy(x, y, xglo, yglo) - disu%cellxy(1,icell) = xglo - disu%cellxy(2,icell) = yglo + disu%cellxy(1, icell) = xglo + disu%cellxy(2, icell) = yglo end do ! if vertices will be needed, it will look like this: @@ -1079,30 +1085,29 @@ subroutine getDiscretization(this, disu) ! 4. get vertex data per cell, add functions to base ! 5. add vertex (x,y) to list and connectivity to sparse ! 6. generate ia/ja from sparse - + end subroutine getDiscretization - !> @brief Deallocate grid connection resources !< subroutine destroy(this) - use MemoryManagerModule, only: mem_deallocate + use MemoryManagerModule, only: mem_deallocate class(GridConnectionType) :: this !< this grid connection instance - + call mem_deallocate(this%nrOfBoundaryCells) call mem_deallocate(this%indexCount) call mem_deallocate(this%nrOfCells) ! arrays - deallocate(this%idxToGlobal) - deallocate(this%boundaryCells) - deallocate(this%connectedCells) - deallocate(this%primConnections) + deallocate (this%idxToGlobal) + deallocate (this%boundaryCells) + deallocate (this%connectedCells) + deallocate (this%primConnections) call mem_deallocate(this%idxToGlobalIdx) - + end subroutine destroy - + !> @brief Test if the connection between nodes within !< the same model is periodic function isPeriodic(this, n, m) result(periodic) @@ -1112,11 +1117,12 @@ function isPeriodic(this, n, m) result(periodic) logical :: periodic !< true when periodic ! local integer(I4B) :: icell - - periodic = .false. + + periodic = .false. do icell = 1, this%nrOfBoundaryCells - if (.not. associated(this%boundaryCells(icell)%cell%model, this%connectedCells(icell)%cell%model)) cycle - + if (.not. associated(this%boundaryCells(icell)%cell%model, & + this%connectedCells(icell)%cell%model)) cycle + ! one way if (this%boundaryCells(icell)%cell%index == n) then if (this%connectedCells(icell)%cell%index == m) then @@ -1131,9 +1137,9 @@ function isPeriodic(this, n, m) result(periodic) return end if end if - + end do - + end function - + end module GridConnectionModule diff --git a/src/Model/Connection/GridSorting.f90 b/src/Model/Connection/GridSorting.f90 index 4c10b28bf84..732816b4e6c 100644 --- a/src/Model/Connection/GridSorting.f90 +++ b/src/Model/Connection/GridSorting.f90 @@ -1,4 +1,4 @@ -module GridSorting +module GridSorting use KindModule, only: I4B, DP, LGP use ConstantsModule, only: DHALF use CellWithNbrsModule, only: GlobalCellType @@ -11,74 +11,76 @@ module GridSorting contains ! Sort an array of integers subroutine quickSortGrid(array, arraySize, idxToGlobal) - integer, intent(inout), dimension(:) :: array - integer, intent(in) :: arraySize - type(GlobalCellType), dimension(:), pointer :: idxToGlobal + integer, intent(inout), dimension(:) :: array + integer, intent(in) :: arraySize + type(GlobalCellType), dimension(:), pointer :: idxToGlobal + ! local + integer :: QSORT_THRESHOLD = 8 + include "qsort_inline.inc" + + contains + subroutine init() + end subroutine init + + ! Compare two grid cells, this doesn't work as + ! smooth for staggered discretizations though... + function lessThan(n, m) result(isLess) + integer(I4B), intent(in) :: n + integer(I4B), intent(in) :: m + logical(LGP) :: isLess ! local - integer :: QSORT_THRESHOLD = 8 - include "qsort_inline.inc" - - contains - subroutine init() - end subroutine init + type(GlobalCellType), pointer :: gcn, gcm + real(DP) :: xnloc, ynloc, xmloc, ymloc + real(DP) :: xn, yn, zn, xm, ym, zm + + ! get coordinates + gcn => idxToGlobal(array(n)) + gcm => idxToGlobal(array(m)) + + ! convert coordinates + call gcn%model%dis%get_cellxy(gcn%index, xnloc, ynloc) + call gcn%model%dis%transform_xy(xnloc, ynloc, xn, yn) + zn = DHALF * (gcn%model%dis%top(gcn%index) + & + gcn%model%dis%bot(gcn%index)) + + call gcm%model%dis%get_cellxy(gcm%index, xmloc, ymloc) + call gcm%model%dis%transform_xy(xmloc, ymloc, xm, ym) + zm = DHALF * (gcm%model%dis%top(gcm%index) + & + gcm%model%dis%bot(gcm%index)) + + ! compare + if (.not. is_same(zn, zm, 10 * epsilon(zn))) then + isLess = zn > zm + else if (.not. is_same(yn, ym, 10 * epsilon(yn))) then + isLess = yn > ym + else if (.not. is_same(xn, xm, 10 * epsilon(xn))) then + isLess = xn < xm + else + isLess = .false. + end if - ! Compare two grid cells, this doesn't work as - ! smooth for staggered discretizations though... - function lessThan(n, m) result(isLess) - integer(I4B), intent(in) :: n - integer(I4B), intent(in) :: m - logical(LGP) :: isLess - ! local - type(GlobalCellType), pointer :: gcn, gcm - real(DP) :: xnloc, ynloc, xmloc, ymloc - real(DP) :: xn, yn, zn, xm, ym, zm - - ! get coordinates - gcn => idxToGlobal(array(n)) - gcm => idxToGlobal(array(m)) - - ! convert coordinates - call gcn%model%dis%get_cellxy(gcn%index, xnloc, ynloc) - call gcn%model%dis%transform_xy(xnloc, ynloc, xn, yn) - zn = DHALF*(gcn%model%dis%top(gcn%index) + gcn%model%dis%bot(gcn%index)) - - call gcm%model%dis%get_cellxy(gcm%index, xmloc, ymloc) - call gcm%model%dis%transform_xy(xmloc, ymloc, xm, ym) - zm = DHALF*(gcm%model%dis%top(gcm%index) + gcm%model%dis%bot(gcm%index)) - - ! compare - if (.not. is_same(zn, zm, 10*epsilon(zn))) then - isLess = zn > zm - else if (.not. is_same(yn, ym, 10*epsilon(yn))) then - isLess = yn > ym - else if (.not. is_same(xn, xm, 10*epsilon(xn))) then - isLess = xn < xm - else - isLess = .false. - end if + end function lessThan - end function lessThan + ! swap indices + subroutine swap(a, b) + integer, intent(in) :: a, b + integer :: hold - ! swap indices - subroutine swap(a,b) - integer, intent(in) :: a,b - integer :: hold + hold = array(a) + array(a) = array(b) + array(b) = hold - hold=array(a) - array(a)=array(b) - array(b)=hold + end subroutine swap - end subroutine swap - - ! circular shift-right by one - subroutine rshift(left,right) - integer, intent(in) :: left, right - integer :: hold + ! circular shift-right by one + subroutine rshift(left, right) + integer, intent(in) :: left, right + integer :: hold - hold=array(right) - array(left+1:right)=array(left:right-1) - array(left)=hold + hold = array(right) + array(left + 1:right) = array(left:right - 1) + array(left) = hold - end subroutine rshift + end subroutine rshift end subroutine quickSortGrid end module GridSorting diff --git a/src/Model/Connection/GwfGwfConnection.f90 b/src/Model/Connection/GwfGwfConnection.f90 index 90a0c70db6f..caf3d17c15b 100644 --- a/src/Model/Connection/GwfGwfConnection.f90 +++ b/src/Model/Connection/GwfGwfConnection.f90 @@ -1,30 +1,30 @@ module GwfGwfConnectionModule use KindModule, only: I4B, DP, LGP - use ConstantsModule, only: DZERO, DONE, DEM6, LENCOMPONENTNAME, LINELENGTH + use ConstantsModule, only: DZERO, DONE, DEM6, LENCOMPONENTNAME, LINELENGTH use CsrUtilsModule, only: getCSRIndex - use SparseModule, only:sparsematrix + use SparseModule, only: sparsematrix use MemoryManagerModule, only: mem_allocate, mem_deallocate use SimModule, only: ustop - use SpatialModelConnectionModule + use SpatialModelConnectionModule use GwfInterfaceModelModule use NumericalModelModule use GwfModule, only: GwfModelType, CastAsGwfModel use DisConnExchangeModule - use GwfGwfExchangeModule, only: GwfExchangeType, GetGwfExchangeFromList, & + use GwfGwfExchangeModule, only: GwfExchangeType, GetGwfExchangeFromList, & CastAsGwfExchange use GwfNpfModule, only: GwfNpfType, hcond, vcond use GwfBuyModule, only: GwfBuyType use BaseDisModule, only: DisBaseType use ConnectionsModule, only: ConnectionsType use CellWithNbrsModule, only: GlobalCellType - + implicit none private public :: CastAsGwfGwfConnection - !> Connecting a GWF model to other models in space, implements - !! NumericalExchangeType so the solution can used this object to determine + !> Connecting a GWF model to other models in space, implements + !! NumericalExchangeType so the solution can used this object to determine !! the coefficients for the coupling between two adjacent models. !< type, public, extends(SpatialModelConnectionType) :: GwfGwfConnectionType @@ -39,13 +39,13 @@ module GwfGwfConnectionModule integer(I4B) :: iout = 0 !< the list file for the interface model real(DP), dimension(:), pointer, contiguous :: exgflowja => null() !< flowja through exchange faces - - contains + + contains procedure, pass(this) :: gwfGwfConnection_ctor generic, public :: construct => gwfGwfConnection_ctor - + ! overriding NumericalExchangeType - procedure :: exg_df => gwfgwfcon_df + procedure :: exg_df => gwfgwfcon_df procedure :: exg_ar => gwfgwfcon_ar procedure :: exg_rp => gwfgwfcon_rp procedure :: exg_ad => gwfgwfcon_ad @@ -58,7 +58,7 @@ module GwfGwfConnectionModule ! overriding 'protected' procedure, pass(this) :: validateConnection - + ! local stuff procedure, pass(this), private :: allocateScalars procedure, pass(this), private :: allocate_arrays @@ -68,18 +68,18 @@ module GwfGwfConnectionModule procedure, pass(this), private :: setFlowToExchange procedure, pass(this), private :: saveExchangeFlows procedure, pass(this), private :: setNpfEdgeProps - + end type GwfGwfConnectionType contains - + !> @brief Basic construction of the connection !< subroutine gwfGwfConnection_ctor(this, model, gwfEx) use NumericalModelModule, only: NumericalModelType use InputOutputModule, only: openfile class(GwfGwfConnectionType) :: this !< the connection - class(NumericalModelType), pointer :: model !< the model owning this connection, + class(NumericalModelType), pointer :: model !< the model owning this connection, !! this must of course be a GwfModelType class(DisConnExchangeType), pointer :: gwfEx !< the exchange the interface model is created for ! local @@ -94,58 +94,60 @@ subroutine gwfGwfConnection_ctor(this, model, gwfEx) this%gwfExchange => CastAsGwfExchange(objPtr) this%exchangeIsOwned = associated(gwfEx%model1, model) - + if (this%exchangeIsOwned) then - write(name,'(a,i0)') 'GWFCON1_', gwfEx%id + write (name, '(a,i0)') 'GWFCON1_', gwfEx%id else - write(name,'(a,i0)') 'GWFCON2_', gwfEx%id + write (name, '(a,i0)') 'GWFCON2_', gwfEx%id end if ! .lst file for interface model if (write_ifmodel_listfile) then fname = trim(name)//'.im.lst' call openfile(this%iout, 0, fname, 'LIST', filstat_opt='REPLACE') - write(this%iout, '(4a)') 'Creating GWF-GWF connection for model ', & - trim(this%gwfModel%name), ' from exchange ', & - trim(gwfEx%name) + write (this%iout, '(4a)') 'Creating GWF-GWF connection for model ', & + trim(this%gwfModel%name), ' from exchange ', & + trim(gwfEx%name) end if - + ! first call base constructor - call this%SpatialModelConnectionType%spatialConnection_ctor(model, gwfEx, name) - + call this%SpatialModelConnectionType%spatialConnection_ctor(model, & + gwfEx, & + name) + call this%allocateScalars() - + this%typename = 'GWF-GWF' this%iXt3dOnExchange = 0 - - allocate(this%gwfInterfaceModel) + + allocate (this%gwfInterfaceModel) this%interfaceModel => this%gwfInterfaceModel - + end subroutine gwfGwfConnection_ctor - + !> @brief Define the connection - !! - !! This sets up the GridConnection (for creating the - !! interface grid), creates and defines the interface + !! + !! This sets up the GridConnection (for creating the + !! interface grid), creates and defines the interface !< model subroutine gwfgwfcon_df(this) - class(GwfGwfConnectionType) :: this !< this connection + class(GwfGwfConnectionType) :: this !< this connection ! local - character(len=LENCOMPONENTNAME) :: imName !< the interface model's name + character(len=LENCOMPONENTNAME) :: imName !< the interface model's name ! determine the required size of the interface grid call this%setGridExtent() - ! this sets up the GridConnection + ! this sets up the GridConnection call this%spatialcon_df() - + ! Now grid conn is defined, we create the interface model ! here, and the remainder of this routine is define. ! we basically follow the logic that is present in sln_df() - if(this%exchangeIsOwned) then - write(imName,'(a,i0)') 'GWFIM1_', this%gwfExchange%id + if (this%exchangeIsOwned) then + write (imName, '(a,i0)') 'GWFIM1_', this%gwfExchange%id else - write(imName,'(a,i0)') 'GWFIM2_', this%gwfExchange%id + write (imName, '(a,i0)') 'GWFIM2_', this%gwfExchange%id end if call this%gwfInterfaceModel%gwfifm_cr(imName, this%iout, this%gridConnection) @@ -160,7 +162,7 @@ subroutine gwfgwfcon_df(this) call this%spatialcon_connect() call this%allocate_arrays() - + end subroutine gwfgwfcon_df !> @brief Set the required size of the interface grid from @@ -169,7 +171,7 @@ subroutine setGridExtent(this) class(GwfGwfConnectionType) :: this !< the connection ! local - this%iXt3dOnExchange = this%gwfExchange%ixt3d + this%iXt3dOnExchange = this%gwfExchange%ixt3d if (this%iXt3dOnExchange > 0) then this%exchangeStencilDepth = 2 if (this%gwfModel%npf%ixt3d > 0) then @@ -178,7 +180,7 @@ subroutine setGridExtent(this) end if end subroutine setGridExtent - + !> @brief allocation of scalars in the connection !< subroutine allocateScalars(this) @@ -198,26 +200,26 @@ subroutine allocate_arrays(this) ! local integer(I4B) :: i - call mem_allocate(this%exgflowja, this%gridConnection%nrOfBoundaryCells, & + call mem_allocate(this%exgflowja, this%gridConnection%nrOfBoundaryCells, & 'EXGFLOWJA', this%memoryPath) do i = 1, size(this%exgflowja) this%exgflowja(i) = 0.0_DP end do end subroutine allocate_arrays - + !> @brief Allocate and read the connection !< subroutine gwfgwfcon_ar(this) - use GridConnectionModule, only: GridConnectionType + use GridConnectionModule, only: GridConnectionType class(GwfGwfConnectionType) :: this !< this connection - ! local + ! local ! check if we can construct an interface model ! NB: only makes sense after the models' allocate&read have been ! called, which is why we do it here call this%validateConnection() - + ! allocate and read base call this%spatialcon_ar() @@ -240,7 +242,7 @@ end subroutine gwfgwfcon_ar !< subroutine gwfgwfcon_rp(this) class(GwfGwfConnectionType) :: this !< this connection - + ! Call exchange rp routines if (this%exchangeIsOwned) then call this%gwfExchange%exg_rp() @@ -259,7 +261,7 @@ subroutine gwfgwfcon_ad(this) ! this triggers the BUY density calculation if (this%gwfInterfaceModel%inbuy > 0) call this%gwfInterfaceModel%buy%buy_ad() - + if (this%exchangeIsOwned) then call this%gwfExchange%exg_ad() end if @@ -274,7 +276,7 @@ subroutine gwfgwfcon_cf(this, kiter) integer(I4B), intent(in) :: kiter !< the iteration counter ! local integer(I4B) :: i - + ! reset interface system do i = 1, this%nja this%amat(i) = 0.0_DP @@ -282,16 +284,16 @@ subroutine gwfgwfcon_cf(this, kiter) do i = 1, this%neq this%rhs(i) = 0.0_DP end do - + ! copy model data into interface model ! (when kiter == 1, this is already done in _ad) if (kiter > 1) call this%syncInterfaceModel() ! calculate (wetting/drying, saturation) call this%gwfInterfaceModel%model_cf(kiter) - + end subroutine gwfgwfcon_cf - + !> @brief Synchronize the interface model !! Fills interface model data from the !! contributing GWF models, at the iteration @@ -301,20 +303,20 @@ subroutine syncInterfaceModel(this) ! local integer(I4B) :: icell, idx class(NumericalModelType), pointer :: model - + ! copy head values - do icell = 1, this%gridConnection%nrOfCells + do icell = 1, this%gridConnection%nrOfCells idx = this%gridConnection%idxToGlobal(icell)%index model => this%gridConnection%idxToGlobal(icell)%model - + this%x(icell) = model%x(idx) this%gwfInterfaceModel%ibound(icell) = model%ibound(idx) this%gwfInterfaceModel%xold(icell) = model%xold(idx) end do - + end subroutine syncInterfaceModel - - !> @brief Write the calculated coefficients into the global + + !> @brief Write the calculated coefficients into the global !< system matrix and the rhs subroutine gwfgwfcon_fc(this, kiter, iasln, amatsln, rhssln, inwtflag) class(GwfGwfConnectionType) :: this !< this connection @@ -325,25 +327,28 @@ subroutine gwfgwfcon_fc(this, kiter, iasln, amatsln, rhssln, inwtflag) integer(I4B), optional, intent(in) :: inwtflag !< newton-raphson flag ! local integer(I4B) :: n, ipos, nglo - + ! fill (and add to...) coefficients for interface call this%gwfInterfaceModel%model_fc(kiter, this%amat, this%nja, inwtflag) - + ! map back to solution matrix do n = 1, this%neq ! we cannot check with the mask here, because cross-terms are not ! necessarily from primary connections. But, we only need the coefficients ! for our own model (i.e. fluxes into cells belonging to this%owner): - if (.not. associated(this%gridConnection%idxToGlobal(n)%model, this%owner)) then + if (.not. associated(this%gridConnection%idxToGlobal(n)%model, & + this%owner)) then ! only add connections for own model to global matrix cycle end if - - nglo = this%gridConnection%idxToGlobal(n)%index + this%gridConnection%idxToGlobal(n)%model%moffset + + nglo = this%gridConnection%idxToGlobal(n)%index + & + this%gridConnection%idxToGlobal(n)%model%moffset rhssln(nglo) = rhssln(nglo) + this%rhs(n) - - do ipos = this%ia(n), this%ia(n+1) - 1 - amatsln(this%mapIdxToSln(ipos)) = amatsln(this%mapIdxToSln(ipos)) + this%amat(ipos) + + do ipos = this%ia(n), this%ia(n + 1) - 1 + amatsln(this%mapIdxToSln(ipos)) = amatsln(this%mapIdxToSln(ipos)) + & + this%amat(ipos) end do end do @@ -358,7 +363,7 @@ subroutine gwfgwfcon_fc(this, kiter, iasln, amatsln, rhssln, inwtflag) end subroutine gwfgwfcon_fc !> @brief Validate this connection - !! This is called before proceeding to construct + !! This is called before proceeding to construct !! the interface model !< subroutine validateConnection(this) @@ -366,14 +371,14 @@ subroutine validateConnection(this) use SimModule, only: count_errors class(GwfGwfConnectionType) :: this !< this connection ! local - + ! base validation (geometry/spatial) call this%SpatialModelConnectionType%validateConnection() call this%validateGwfExchange() ! abort on errors - if(count_errors() > 0) then - write(errmsg, '(1x,a)') 'Errors occurred while processing exchange(s)' + if (count_errors() > 0) then + write (errmsg, '(1x,a)') 'Errors occurred while processing exchange(s)' call ustop() end if @@ -389,7 +394,7 @@ subroutine validateGwfExchange(this) use SimVariablesModule, only: errmsg use SimModule, only: store_error use GwfNpfModule, only: GwfNpfType - class(GwfGwfConnectionType) :: this !< this connection + class(GwfGwfConnectionType) :: this !< this connection ! local class(GwfExchangeType), pointer :: gwfEx class(*), pointer :: modelPtr @@ -406,18 +411,18 @@ subroutine validateGwfExchange(this) ! GNC not allowed if (gwfEx%ingnc /= 0) then - write(errmsg, '(1x,2a)') 'Ghost node correction not supported '// & - 'with interface model for exchange', & - trim(gwfEx%name) + write (errmsg, '(1x,2a)') 'Ghost node correction not supported '// & + 'with interface model for exchange', & + trim(gwfEx%name) call store_error(errmsg) end if - if ((gwfModel1%inbuy > 0 .and. gwfModel2%inbuy == 0) .or. & + if ((gwfModel1%inbuy > 0 .and. gwfModel2%inbuy == 0) .or. & (gwfModel1%inbuy == 0 .and. gwfModel2%inbuy > 0)) then - write(errmsg, '(1x,2a)') 'Buoyancy package should be enabled/disabled '// & - 'simultaneously in models connected with the '// & - 'interface model for exchange ', & - trim(gwfEx%name) + write (errmsg, '(1x,2a)') 'Buoyancy package should be enabled/disabled '// & + 'simultaneously in models connected with the '// & + 'interface model for exchange ', & + trim(gwfEx%name) call store_error(errmsg) end if @@ -425,10 +430,10 @@ subroutine validateGwfExchange(this) if (gwfModel1%inbuy > 0 .and. gwfModel2%inbuy > 0) then ! does not work with XT3D if (this%iXt3dOnExchange > 0) then - write(errmsg, '(1x,2a)') 'Connecting models with BUY package not '// & - 'allowed with XT3D enabled on exchange ', & - trim(gwfEx%name) - call store_error(errmsg) + write (errmsg, '(1x,2a)') 'Connecting models with BUY package not '// & + 'allowed with XT3D enabled on exchange ', & + trim(gwfEx%name) + call store_error(errmsg) end if ! check compatibility of buoyancy @@ -445,12 +450,12 @@ subroutine validateGwfExchange(this) end if if (.not. compatible) then - write(errmsg, '(1x,6a)') 'Buoyancy packages in model ', & - trim(gwfEx%model1%name), ' and ', & - trim(gwfEx%model2%name), & - ' should be equivalent to construct an '// & - ' interface model for exchange ', & - trim(gwfEx%name) + write (errmsg, '(1x,6a)') 'Buoyancy packages in model ', & + trim(gwfEx%model1%name), ' and ', & + trim(gwfEx%model2%name), & + ' should be equivalent to construct an '// & + ' interface model for exchange ', & + trim(gwfEx%name) call store_error(errmsg) end if @@ -460,7 +465,7 @@ end subroutine validateGwfExchange !> @brief Deallocate all resources !< - subroutine gwfgwfcon_da(this) + subroutine gwfgwfcon_da(this) use KindModule, only: LGP class(GwfGwfConnectionType) :: this !< this connection ! local @@ -471,22 +476,22 @@ subroutine gwfgwfcon_da(this) ! arrays call mem_deallocate(this%exgflowja) - + call this%gwfInterfaceModel%model_da() - deallocate(this%gwfInterfaceModel) - + deallocate (this%gwfInterfaceModel) + call this%spatialcon_da() - inquire(this%iout, opened=isOpen) + inquire (this%iout, opened=isOpen) if (isOpen) then - close(this%iout) + close (this%iout) end if ! we need to deallocate the baseexchange we own: if (this%exchangeIsOwned) then call this%gwfExchange%exg_da() end if - + end subroutine gwfgwfcon_da !> @brief Calculate intra-cell flows @@ -498,8 +503,8 @@ subroutine gwfgwfcon_cq(this, icnvg, isuppress_output, isolnid) integer(I4B), intent(inout) :: icnvg !< convergence flag integer(I4B), intent(in) :: isuppress_output !< suppress output when =1 integer(I4B), intent(in) :: isolnid !< solution id - - call this%gwfInterfaceModel%model_cq(icnvg, isuppress_output) + + call this%gwfInterfaceModel%model_cq(icnvg, isuppress_output) call this%setFlowToExchange() @@ -509,7 +514,7 @@ subroutine gwfgwfcon_cq(this, icnvg, isuppress_output, isolnid) ! simvals? ! if needed, we add the edge properties to the model's NPF ! package for its spdis calculation: - if (this%gwfModel%npf%icalcspdis == 1) then + if (this%gwfModel%npf%icalcspdis == 1) then call this%setNpfEdgeProps() end if @@ -521,10 +526,9 @@ subroutine gwfgwfcon_cq(this, icnvg, isuppress_output, isolnid) call this%gwfExchange%gwf_gwf_add_to_flowja() end if - end subroutine gwfgwfcon_cq - !> @brief Set the flows (flowja from interface model) to the + !> @brief Set the flows (flowja from interface model) to the !< simvals in the exchange, leaving the budget calcution in there subroutine setFlowToExchange(this) class(GwfGwfConnectionType) :: this !< this connection @@ -534,16 +538,19 @@ subroutine setFlowToExchange(this) class(GwfExchangeType), pointer :: gwfEx gwfEx => this%gwfExchange - if (this%exchangeIsOwned) then + if (this%exchangeIsOwned) then do i = 1, gwfEx%nexg gwfEx%simvals(i) = DZERO - if (gwfEx%gwfmodel1%ibound(gwfEx%nodem1(i)) /= 0 .and. & + if (gwfEx%gwfmodel1%ibound(gwfEx%nodem1(i)) /= 0 .and. & gwfEx%gwfmodel2%ibound(gwfEx%nodem2(i)) /= 0) then - nIface = this%gridConnection%getInterfaceIndex(gwfEx%nodem1(i), gwfEx%model1) - mIface = this%gridConnection%getInterfaceIndex(gwfEx%nodem2(i), gwfEx%model2) - ipos = getCSRIndex(nIface, mIface, this%gwfInterfaceModel%ia, this%gwfInterfaceModel%ja) + nIface = this%gridConnection%getInterfaceIndex(gwfEx%nodem1(i), & + gwfEx%model1) + mIface = this%gridConnection%getInterfaceIndex(gwfEx%nodem2(i), & + gwfEx%model2) + ipos = getCSRIndex(nIface, mIface, this%gwfInterfaceModel%ia, & + this%gwfInterfaceModel%ja) gwfEx%simvals(i) = this%gwfInterfaceModel%flowja(ipos) end if @@ -563,11 +570,12 @@ subroutine saveExchangeFlows(this) do i = 1, this%gridConnection%nrOfBoundaryCells boundaryCell = this%gridConnection%boundaryCells(i)%cell connectedCell = this%gridConnection%connectedCells(i)%cell - n = this%gridConnection%getInterfaceIndex(boundaryCell%index, & + n = this%gridConnection%getInterfaceIndex(boundaryCell%index, & boundaryCell%model) - m = this%gridConnection%getInterfaceIndex(connectedCell%index, & + m = this%gridConnection%getInterfaceIndex(connectedCell%index, & connectedCell%model) - ipos = getCSRIndex(n, m, this%gwfInterfaceModel%ia, this%gwfInterfaceModel%ja) + ipos = getCSRIndex(n, m, this%gwfInterfaceModel%ia, & + this%gwfInterfaceModel%ja) this%exgflowja(i) = this%gwfInterfaceModel%flowja(ipos) end do @@ -617,26 +625,26 @@ subroutine setNpfEdgeProps(this) nLoc = toGlobal(n)%index - do ipos = imCon%ia(n)+1, imCon%ia(n+1) - 1 + do ipos = imCon%ia(n) + 1, imCon%ia(n + 1) - 1 if (imCon%mask(ipos) < 1) then ! skip this connection, it's masked so not determined by us cycle end if m = imCon%ja(ipos) - mLoc = toGlobal(m)%index + mLoc = toGlobal(m)%index if (.not. associated(toGlobal(m)%model, this%owner)) then ! boundary connection, set edge properties isym = imCon%jas(ipos) ihc = imCon%ihc(isym) - area = imCon%hwva(isym) + area = imCon%hwva(isym) satThick = imNpf%calcSatThickness(n, m, ihc) rrate = this%gwfInterfaceModel%flowja(ipos) - call imDis%connection_normal(n, m, ihc, nx, ny, nz, ipos) + call imDis%connection_normal(n, m, ihc, nx, ny, nz, ipos) call imDis%connection_vector(n, m, nozee, imNpf%sat(n), imNpf%sat(m), & - ihc, cx, cy, cz, conLen) + ihc, cx, cy, cz, conLen) if (ihc == 0) then ! check if n is below m @@ -650,8 +658,8 @@ subroutine setNpfEdgeProps(this) cl = imCon%cl2(isym) end if dist = conLen * cl / (imCon%cl1(isym) + imCon%cl2(isym)) - call this%gwfModel%npf%set_edge_properties(nLoc, ihc, rrate, area, & - nx, ny, dist) + call this%gwfModel%npf%set_edge_properties(nLoc, ihc, rrate, area, & + nx, ny, dist) else ! internal, need to set flowja for n-m ! TODO_MJR: should we mask the flowja calculation in the model? @@ -679,7 +687,7 @@ subroutine gwfgwfcon_bd(this, icnvg, isuppress_output, isolnid) if (this%exchangeIsOwned) then call this%gwfExchange%exg_bd(icnvg, isuppress_output, isolnid) end if - + end subroutine gwfgwfcon_bd !> @brief Write output for exchange (and calls @@ -687,7 +695,7 @@ end subroutine gwfgwfcon_bd subroutine gwfgwfcon_ot(this) class(GwfGwfConnectionType) :: this !< this connection ! local - + ! Call exg_ot() here as it handles all output processing ! based on gwfExchange%simvals(:), which was correctly ! filled from gwfgwfcon @@ -699,14 +707,14 @@ end subroutine gwfgwfcon_ot !> @brief Cast to GwfGwfConnectionType !< - function CastAsGwfGwfConnection(obj) result (res) + function CastAsGwfGwfConnection(obj) result(res) implicit none class(*), pointer, intent(inout) :: obj !< object to be cast class(GwfGwfConnectionType), pointer :: res !< the GwfGwfConnection - + res => null() if (.not. associated(obj)) return - + select type (obj) class is (GwfGwfConnectionType) res => obj diff --git a/src/Model/Connection/GwfInterfaceModel.f90 b/src/Model/Connection/GwfInterfaceModel.f90 index 1fa930292dd..5c87d62f200 100644 --- a/src/Model/Connection/GwfInterfaceModel.f90 +++ b/src/Model/Connection/GwfInterfaceModel.f90 @@ -1,6 +1,6 @@ module GwfInterfaceModelModule use KindModule, only: I4B, DP - use ConstantsModule, only: DZERO + use ConstantsModule, only: DZERO use MemoryManagerModule, only: mem_allocate use MemoryHelperModule, only: create_mem_path use NumericalModelModule, only: NumericalModelType, GetNumericalModelFromList @@ -19,11 +19,11 @@ module GwfInterfaceModelModule private !> The GWF Interface Model is a utility to calculate the solution's - !! exchange coefficients from the interface between a GWF model and - !! its GWF neighbors. The interface model itself will not be part - !! of the solution, it is not being solved. + !! exchange coefficients from the interface between a GWF model and + !! its GWF neighbors. The interface model itself will not be part + !! of the solution, it is not being solved. !! Patching (a part of the) discretizations of two GWF models in a - !! general way, e.g. DIS+DIS with refinement, requires the resulting + !! general way, e.g. DIS+DIS with refinement, requires the resulting !< discretization to be of type DISU. type, public, extends(GwfModelType) :: GwfInterfaceModelType class(GridConnectionType), pointer :: gridConnection => null() !< The grid connection class will provide the interface grid @@ -33,17 +33,17 @@ module GwfInterfaceModelModule procedure, pass(this) :: gwfifm_cr procedure :: model_df => gwfifm_df procedure :: model_ar => gwfifm_ar - procedure :: model_da => gwfifm_da + procedure :: model_da => gwfifm_da ! private procedure, private, pass(this) :: setNpfOptions procedure, private, pass(this) :: setNpfGridData procedure, private, pass(this) :: setBuyData end type - + contains - - !> @brief set up the interface model, analogously to what + + !> @brief set up the interface model, analogously to what !< happens in gwf_cr subroutine gwfifm_cr(this, name, iout, gridConn) class(GwfInterfaceModelType) :: this !< the GWF interface model @@ -52,32 +52,32 @@ subroutine gwfifm_cr(this, name, iout, gridConn) class(GridConnectionType), pointer, intent(in) :: gridConn !< the grid connection for creating a DISU ! local class(*), pointer :: modPtr - + this%memoryPath = create_mem_path(name) call this%allocate_scalars(name) - - this%iout = iout + + this%iout = iout this%gridConnection => gridConn modPtr => this%gridConnection%model this%owner => CastAsGwfModel(modPtr) - + this%innpf = huge(1_I4B) this%inewton = this%owner%inewton this%inewtonur = this%owner%inewtonur - + if (this%owner%inbuy > 0) then this%inbuy = huge(1_I4B) end if - + ! create discretization and packages call disu_cr(this%dis, this%name, -1, this%iout) call npf_cr(this%npf, this%name, this%innpf, this%iout) call xt3d_cr(this%xt3d, this%name, this%innpf, this%iout) call buy_cr(this%buy, this%name, this%inbuy, this%iout) - + end subroutine gwfifm_cr - + !> @brief Define, mostly DISU and the NPF package !< for this interface model subroutine gwfifm_df(this) @@ -106,39 +106,38 @@ subroutine gwfifm_df(this) call this%buy%buy_df(this%dis, buyData) call buyData%destruct() end if - + this%neq = this%dis%nodes this%nja = this%dis%nja - this%ia => this%dis%con%ia - this%ja => this%dis%con%ja - + this%ia => this%dis%con%ia + this%ja => this%dis%con%ja + call this%allocate_arrays() - + end subroutine gwfifm_df - + !> @brief allocate and read the packages !< subroutine gwfifm_ar(this) class(GwfInterfaceModelType) :: this !< the GWF interface model ! local type(GwfNpfGridDataType) :: npfGridData - + call npfGridData%construct(this%dis%nodes) call this%setNpfGridData(npfGridData) call this%npf%npf_ar(this%ic, this%ibound, this%x, npfGridData) call npfGridData%destroy() if (this%inbuy > 0) call this%buy%buy_ar(this%npf, this%ibound) - + end subroutine gwfifm_ar - !> @brief Clean up !< subroutine gwfifm_da(this) - use MemoryManagerModule, only: mem_deallocate + use MemoryManagerModule, only: mem_deallocate class(GwfInterfaceModelType) :: this !< the GWF interface model - + ! -- Internal flow packages deallocate call this%dis%dis_da() call this%npf%npf_da() @@ -146,9 +145,9 @@ subroutine gwfifm_da(this) call this%buy%buy_da() ! ! -- Internal package objects - deallocate(this%dis) - deallocate(this%npf) - deallocate(this%xt3d) + deallocate (this%dis) + deallocate (this%npf) + deallocate (this%xt3d) ! ! -- Scalars call mem_deallocate(this%inic) @@ -166,9 +165,9 @@ subroutine gwfifm_da(this) ! ! -- NumericalModelType call this%NumericalModelType%model_da() - + end subroutine - + !> @brief Copy NPF options from the model owning !! the interface to the data structure !< @@ -190,8 +189,8 @@ subroutine setNpfOptions(this, npfOptions) end subroutine setNpfOptions - !> @brief Loop over the interface grid and fill the structure - !! with NPF grid data, copied from the models that participate + !> @brief Loop over the interface grid and fill the structure + !! with NPF grid data, copied from the models that participate !! in this interface !< subroutine setNpfGridData(this, npfGridData) @@ -205,7 +204,7 @@ subroutine setNpfGridData(this, npfGridData) ! TODO_MJR: deal with inhomogeneity, for now, we assume ! that we can just take the owning model's settings... npfGridData%ik22 = this%owner%npf%ik22 - npfGridData%ik33 = this%owner%npf%ik33 + npfGridData%ik33 = this%owner%npf%ik33 npfGridData%iwetdry = this%owner%npf%iwetdry npfGridData%iangle1 = this%owner%npf%iangle1 npfGridData%iangle2 = this%owner%npf%iangle2 @@ -215,7 +214,7 @@ subroutine setNpfGridData(this, npfGridData) npfGridData%iangle2 = 1 npfGridData%iangle3 = 1 end if - + do icell = 1, this%gridConnection%nrOfCells idx = this%gridConnection%idxToGlobal(icell)%index modelPtr => this%gridConnection%idxToGlobal(icell)%model @@ -227,9 +226,9 @@ subroutine setNpfGridData(this, npfGridData) npfGridData%k33(icell) = gwfModel%npf%k33(idx) ! the K rotation angles, or default (0.0) - if (npfGridData%iangle1 == 1) then + if (npfGridData%iangle1 == 1) then if (gwfModel%npf%iangle1 == 1) then - npfGridData%angle1(icell) = gwfModel%npf%angle1(idx) + npfGridData%angle1(icell) = gwfModel%npf%angle1(idx) else npfGridData%angle1(icell) = DZERO end if @@ -262,8 +261,8 @@ subroutine setNpfGridData(this, npfGridData) end subroutine setNpfGridData - !> @brief Sets the BUY input data from the models that - !! make up this interface. We adopt everything from the + !> @brief Sets the BUY input data from the models that + !! make up this interface. We adopt everything from the !! owning model, but during validation it should be !< checked that the models are compatible. subroutine setBuyData(this, buyData) @@ -284,5 +283,5 @@ subroutine setBuyData(this, buyData) end do end subroutine setBuyData - + end module GwfInterfaceModelModule diff --git a/src/Model/Connection/GwtGwtConnection.f90 b/src/Model/Connection/GwtGwtConnection.f90 index fb9cef32314..b1f8eb781df 100644 --- a/src/Model/Connection/GwtGwtConnection.f90 +++ b/src/Model/Connection/GwtGwtConnection.f90 @@ -1,6 +1,6 @@ module GwtGwtConnectionModule use KindModule, only: I4B, DP, LGP - use ConstantsModule, only: LINELENGTH, LENCOMPONENTNAME, DZERO, LENBUDTXT + use ConstantsModule, only: LINELENGTH, LENCOMPONENTNAME, DZERO, LENBUDTXT use CsrUtilsModule, only: getCSRIndex use SimModule, only: ustop use MemoryManagerModule, only: mem_allocate, mem_deallocate @@ -37,15 +37,15 @@ module GwtGwtConnectionModule integer(I4B), pointer :: exgflowSign => null() !< indicates the flow direction of exgflowja real(DP), dimension(:), pointer, contiguous :: exgflowjaGwt => null() !< gwt-flowja at the interface (this is a subset of the GWT !! interface model flowja's) - + real(DP), dimension(:), pointer, contiguous :: gwfflowja => null() !< gwfflowja for the interface model real(DP), dimension(:), pointer, contiguous :: gwfsat => null() !< gwfsat for the interface model real(DP), dimension(:), pointer, contiguous :: gwfhead => null() !< gwfhead for the interface model - real(DP), dimension(:,:), pointer, contiguous :: gwfspdis => null() !< gwfspdis for the interface model + real(DP), dimension(:, :), pointer, contiguous :: gwfspdis => null() !< gwfspdis for the interface model real(DP), dimension(:), pointer, contiguous :: conc => null() !< pointer to concentration array integer(I4B), dimension(:), pointer, contiguous :: icbound => null() !< store pointer to gwt ibound array - + integer(I4B) :: iout = 0 !< the list file for the interface model contains @@ -81,444 +81,451 @@ module GwtGwtConnectionModule !> @brief Basic construction of the connection !< -subroutine gwtGwtConnection_ctor(this, model, gwtEx) - use InputOutputModule, only: openfile - class(GwtGwtConnectionType) :: this !< the connection - class(NumericalModelType), pointer :: model !< the model owning this connection, + subroutine gwtGwtConnection_ctor(this, model, gwtEx) + use InputOutputModule, only: openfile + class(GwtGwtConnectionType) :: this !< the connection + class(NumericalModelType), pointer :: model !< the model owning this connection, !! this must be a GwtModelType - class(DisConnExchangeType), pointer :: gwtEx !< the GWT-GWT exchange the interface model is created for - ! local - character(len=LINELENGTH) :: fname - character(len=LENCOMPONENTNAME) :: name - class(*), pointer :: objPtr - logical(LGP) :: write_ifmodel_listfile = .false. - - objPtr => model - this%gwtModel => CastAsGwtModel(objPtr) - objPtr => gwtEx - this%gwtExchange => CastAsGwtExchange(objPtr) - - this%exchangeIsOwned = associated(model, gwtEx%model1) - - if (this%exchangeIsOwned) then - write(name,'(a,i0)') 'GWTCON1_', gwtEx%id - else - write(name,'(a,i0)') 'GWTCON2_', gwtEx%id - end if - - ! .lst file for interface model - if (write_ifmodel_listfile) then - fname = trim(name)//'.im.lst' - call openfile(this%iout, 0, fname, 'LIST', filstat_opt='REPLACE') - write(this%iout, '(4a)') 'Creating GWT-GWT connection for model ', & - trim(this%gwtModel%name), 'from exchange ', & - trim(gwtEx%name) - end if - - ! first call base constructor - call this%SpatialModelConnectionType%spatialConnection_ctor(model, gwtEx, name) - - call this%allocate_scalars() - this%typename = 'GWT-GWT' - this%iIfaceAdvScheme = 0 - this%iIfaceXt3d = 1 - this%exgflowSign = 1 - - allocate(this%gwtInterfaceModel) - this%interfaceModel => this%gwtInterfaceModel - -end subroutine gwtGwtConnection_ctor + class(DisConnExchangeType), pointer :: gwtEx !< the GWT-GWT exchange the interface model is created for + ! local + character(len=LINELENGTH) :: fname + character(len=LENCOMPONENTNAME) :: name + class(*), pointer :: objPtr + logical(LGP) :: write_ifmodel_listfile = .false. + + objPtr => model + this%gwtModel => CastAsGwtModel(objPtr) + objPtr => gwtEx + this%gwtExchange => CastAsGwtExchange(objPtr) + + this%exchangeIsOwned = associated(model, gwtEx%model1) + + if (this%exchangeIsOwned) then + write (name, '(a,i0)') 'GWTCON1_', gwtEx%id + else + write (name, '(a,i0)') 'GWTCON2_', gwtEx%id + end if + + ! .lst file for interface model + if (write_ifmodel_listfile) then + fname = trim(name)//'.im.lst' + call openfile(this%iout, 0, fname, 'LIST', filstat_opt='REPLACE') + write (this%iout, '(4a)') 'Creating GWT-GWT connection for model ', & + trim(this%gwtModel%name), 'from exchange ', & + trim(gwtEx%name) + end if + + ! first call base constructor + call this%SpatialModelConnectionType%spatialConnection_ctor(model, & + gwtEx, & + name) + + call this%allocate_scalars() + this%typename = 'GWT-GWT' + this%iIfaceAdvScheme = 0 + this%iIfaceXt3d = 1 + this%exgflowSign = 1 + + allocate (this%gwtInterfaceModel) + this%interfaceModel => this%gwtInterfaceModel + + end subroutine gwtGwtConnection_ctor !> @brief Allocate scalar variables for this connection !< -subroutine allocate_scalars(this) - class(GwtGwtConnectionType) :: this !< the connection + subroutine allocate_scalars(this) + class(GwtGwtConnectionType) :: this !< the connection - call mem_allocate(this%iIfaceAdvScheme, 'IADVSCHEME', this%memoryPath) - call mem_allocate(this%iIfaceXt3d, 'IXT3D', this%memoryPath) - call mem_allocate(this%exgflowSign, 'EXGFLOWSIGN', this%memoryPath) + call mem_allocate(this%iIfaceAdvScheme, 'IADVSCHEME', this%memoryPath) + call mem_allocate(this%iIfaceXt3d, 'IXT3D', this%memoryPath) + call mem_allocate(this%exgflowSign, 'EXGFLOWSIGN', this%memoryPath) -end subroutine allocate_scalars + end subroutine allocate_scalars !> @brief Allocate array variables for this connection !< -subroutine allocate_arrays(this) - class(GwtGwtConnectionType) :: this !< the connection - ! local - integer(I4B) :: i + subroutine allocate_arrays(this) + class(GwtGwtConnectionType) :: this !< the connection + ! local + integer(I4B) :: i - call mem_allocate(this%gwfflowja, this%interfaceModel%nja, 'GWFFLOWJA', & - this%memoryPath) - call mem_allocate(this%gwfsat, this%neq, 'GWFSAT', this%memoryPath) - call mem_allocate(this%gwfhead, this%neq, 'GWFHEAD', this%memoryPath) - call mem_allocate(this%gwfspdis, 3, this%neq, 'GWFSPDIS', this%memoryPath) + call mem_allocate(this%gwfflowja, this%interfaceModel%nja, 'GWFFLOWJA', & + this%memoryPath) + call mem_allocate(this%gwfsat, this%neq, 'GWFSAT', this%memoryPath) + call mem_allocate(this%gwfhead, this%neq, 'GWFHEAD', this%memoryPath) + call mem_allocate(this%gwfspdis, 3, this%neq, 'GWFSPDIS', this%memoryPath) - call mem_allocate(this%exgflowjaGwt, this%gridConnection%nrOfBoundaryCells, & - 'EXGFLOWJAGWT', this%memoryPath) + call mem_allocate(this%exgflowjaGwt, this%gridConnection%nrOfBoundaryCells, & + 'EXGFLOWJAGWT', this%memoryPath) - do i = 1, size(this%gwfflowja) - this%gwfflowja = 0.0_DP - end do + do i = 1, size(this%gwfflowja) + this%gwfflowja = 0.0_DP + end do - do i = 1, this%neq - this%gwfsat = 0.0_DP - end do + do i = 1, this%neq + this%gwfsat = 0.0_DP + end do -end subroutine allocate_arrays + end subroutine allocate_arrays !> @brief define the GWT-GWT connection !< -subroutine gwtgwtcon_df(this) - class(GwtGwtConnectionType) :: this !< the connection - ! local - character(len=LENCOMPONENTNAME) :: imName - - ! determine advection scheme (the GWT-GWT exchange - ! has been read at this point) - this%iIfaceAdvScheme = this%gwtExchange%iAdvScheme - - ! determine xt3d setting on interface - this%iIfaceXt3d = this%gwtExchange%ixt3d - - ! determine the required size of the interface model grid - call this%setGridExtent() - - ! now set up the GridConnection - call this%spatialcon_df() - - ! we have to 'catch up' and create the interface model - ! here, then the remainder of this routine will be define - if (this%exchangeIsOwned) then - write(imName,'(a,i0)') 'GWTIM1_', this%gwtExchange%id - else - write(imName,'(a,i0)') 'GWTIM2_', this%gwtExchange%id - end if - call this%gwtInterfaceModel%gwtifmod_cr(imName, this%iout, this%gridConnection) - this%gwtInterfaceModel%iAdvScheme = this%iIfaceAdvScheme - this%gwtInterfaceModel%ixt3d = this%iIfaceXt3d - call this%gwtInterfaceModel%model_df() - - call this%allocate_arrays() - - ! connect X, RHS, IBOUND, and flowja - call this%spatialcon_setmodelptrs() - - this%gwtInterfaceModel%fmi%gwfflowja => this%gwfflowja - this%gwtInterfaceModel%fmi%gwfsat => this%gwfsat - this%gwtInterfaceModel%fmi%gwfhead => this%gwfhead - this%gwtInterfaceModel%fmi%gwfspdis => this%gwfspdis - - ! connect pointers (used by BUY) - this%conc => this%gwtInterfaceModel%x - this%icbound => this%gwtInterfaceModel%ibound - - ! add connections from the interface model to solution matrix - call this%spatialcon_connect() - -end subroutine gwtgwtcon_df + subroutine gwtgwtcon_df(this) + class(GwtGwtConnectionType) :: this !< the connection + ! local + character(len=LENCOMPONENTNAME) :: imName + + ! determine advection scheme (the GWT-GWT exchange + ! has been read at this point) + this%iIfaceAdvScheme = this%gwtExchange%iAdvScheme + + ! determine xt3d setting on interface + this%iIfaceXt3d = this%gwtExchange%ixt3d + + ! determine the required size of the interface model grid + call this%setGridExtent() + + ! now set up the GridConnection + call this%spatialcon_df() + + ! we have to 'catch up' and create the interface model + ! here, then the remainder of this routine will be define + if (this%exchangeIsOwned) then + write (imName, '(a,i0)') 'GWTIM1_', this%gwtExchange%id + else + write (imName, '(a,i0)') 'GWTIM2_', this%gwtExchange%id + end if + call this%gwtInterfaceModel%gwtifmod_cr(imName, & + this%iout, & + this%gridConnection) + this%gwtInterfaceModel%iAdvScheme = this%iIfaceAdvScheme + this%gwtInterfaceModel%ixt3d = this%iIfaceXt3d + call this%gwtInterfaceModel%model_df() + + call this%allocate_arrays() + + ! connect X, RHS, IBOUND, and flowja + call this%spatialcon_setmodelptrs() + + this%gwtInterfaceModel%fmi%gwfflowja => this%gwfflowja + this%gwtInterfaceModel%fmi%gwfsat => this%gwfsat + this%gwtInterfaceModel%fmi%gwfhead => this%gwfhead + this%gwtInterfaceModel%fmi%gwfspdis => this%gwfspdis + + ! connect pointers (used by BUY) + this%conc => this%gwtInterfaceModel%x + this%icbound => this%gwtInterfaceModel%ibound + + ! add connections from the interface model to solution matrix + call this%spatialcon_connect() + + end subroutine gwtgwtcon_df !> @brief Set required extent of the interface grid from !< the configuration -subroutine setGridExtent(this) - class(GwtGwtConnectionType) :: this !< the connection - ! local - logical(LGP) :: hasAdv, hasDsp - - hasAdv = this%gwtModel%inadv > 0 - hasDsp = this%gwtModel%indsp > 0 - - if (hasAdv) then - if (this%iIfaceAdvScheme == 2) then - this%exchangeStencilDepth = 2 - if (this%gwtModel%adv%iadvwt == 2) then - this%internalStencilDepth = 2 + subroutine setGridExtent(this) + class(GwtGwtConnectionType) :: this !< the connection + ! local + logical(LGP) :: hasAdv, hasDsp + + hasAdv = this%gwtModel%inadv > 0 + hasDsp = this%gwtModel%indsp > 0 + + if (hasAdv) then + if (this%iIfaceAdvScheme == 2) then + this%exchangeStencilDepth = 2 + if (this%gwtModel%adv%iadvwt == 2) then + this%internalStencilDepth = 2 + end if end if end if - end if - if (hasDsp) then - if (this%iIfaceXt3d > 0) then - this%exchangeStencilDepth = 2 - if (this%gwtModel%dsp%ixt3d > 0) then - this%internalStencilDepth = 2 + if (hasDsp) then + if (this%iIfaceXt3d > 0) then + this%exchangeStencilDepth = 2 + if (this%gwtModel%dsp%ixt3d > 0) then + this%internalStencilDepth = 2 + end if end if end if - end if -end subroutine setGridExtent + end subroutine setGridExtent !> @brief allocate and read/set the connection's data structures !< -subroutine gwtgwtcon_ar(this) - class(GwtGwtConnectionType) :: this !< the connection - ! local - integer(I4B) :: i, idx - class(GwtModelType), pointer :: gwtModel - class(*), pointer :: modelPtr - - ! check if we can construct an interface model - ! NB: only makes sense after the models' allocate&read have been - ! called, which is why we do it here - call this%validateConnection() - - ! fill porosity from mst packages, needed for dsp - if (this%gwtModel%inmst > 0) then - do i = 1, this%neq - modelPtr => this%gridConnection%idxToGlobal(i)%model - gwtModel => CastAsGwtModel(modelPtr) - idx = this%gridConnection%idxToGlobal(i)%index - this%gwtInterfaceModel%porosity(i) = gwtModel%mst%porosity(idx) - end do - end if - - ! allocate and read base - call this%spatialcon_ar() - - ! ... and now the interface model - call this%gwtInterfaceModel%model_ar() - - ! AR the movers and obs through the exchange - if (this%exchangeIsOwned) then - !cdl implement this when MVT is ready - !cdl if (this%gwtExchange%inmvt > 0) then - !cdl call this%gwtExchange%mvt%mvt_ar() - !cdl end if - if (this%gwtExchange%inobs > 0) then - call this%gwtExchange%obs%obs_ar() + subroutine gwtgwtcon_ar(this) + class(GwtGwtConnectionType) :: this !< the connection + ! local + integer(I4B) :: i, idx + class(GwtModelType), pointer :: gwtModel + class(*), pointer :: modelPtr + + ! check if we can construct an interface model + ! NB: only makes sense after the models' allocate&read have been + ! called, which is why we do it here + call this%validateConnection() + + ! fill porosity from mst packages, needed for dsp + if (this%gwtModel%inmst > 0) then + do i = 1, this%neq + modelPtr => this%gridConnection%idxToGlobal(i)%model + gwtModel => CastAsGwtModel(modelPtr) + idx = this%gridConnection%idxToGlobal(i)%index + this%gwtInterfaceModel%porosity(i) = gwtModel%mst%porosity(idx) + end do end if - end if -end subroutine gwtgwtcon_ar + ! allocate and read base + call this%spatialcon_ar() + + ! ... and now the interface model + call this%gwtInterfaceModel%model_ar() -!> @brief validate this connection prior to constructing + ! AR the movers and obs through the exchange + if (this%exchangeIsOwned) then + !cdl implement this when MVT is ready + !cdl if (this%gwtExchange%inmvt > 0) then + !cdl call this%gwtExchange%mvt%mvt_ar() + !cdl end if + if (this%gwtExchange%inobs > 0) then + call this%gwtExchange%obs%obs_ar() + end if + end if + + end subroutine gwtgwtcon_ar + +!> @brief validate this connection prior to constructing !< the interface model -subroutine validateConnection(this) - use SimVariablesModule, only: errmsg - use SimModule, only: count_errors, store_error - class(GwtGwtConnectionType) :: this !< this connection - - ! base validation, the spatial/geometry part - call this%SpatialModelConnectionType%validateConnection() - - ! GWT related matters - if ((this%gwtExchange%gwtmodel1%inadv > 0 .and. this%gwtExchange%gwtmodel2%inadv == 0) .or. & - (this%gwtExchange%gwtmodel2%inadv > 0 .and. this%gwtExchange%gwtmodel1%inadv == 0)) then - write(errmsg, '(1x,a,a,a)') 'Cannot connect GWT models in exchange ', & - trim(this%gwtExchange%name), ' because one model is configured with ADV & - &and the other one is not' - call store_error(errmsg) - end if - - if ((this%gwtExchange%gwtmodel1%indsp > 0 .and. this%gwtExchange%gwtmodel2%indsp == 0) .or. & - (this%gwtExchange%gwtmodel2%indsp > 0 .and. this%gwtExchange%gwtmodel1%indsp == 0)) then - write(errmsg, '(1x,a,a,a)') 'Cannot connect GWT models in exchange ', & - trim(this%gwtExchange%name), ' because one model is configured with DSP & - &and the other one is not' - call store_error(errmsg) - end if - - ! abort on errors - if(count_errors() > 0) then - write(errmsg, '(1x,a)') 'Errors occurred while processing exchange(s)' - call ustop() - end if - -end subroutine validateConnection + subroutine validateConnection(this) + use SimVariablesModule, only: errmsg + use SimModule, only: count_errors, store_error + class(GwtGwtConnectionType) :: this !< this connection + ! base validation, the spatial/geometry part + call this%SpatialModelConnectionType%validateConnection() + + ! GWT related matters + if ((this%gwtExchange%gwtmodel1%inadv > 0 .and. & + this%gwtExchange%gwtmodel2%inadv == 0) .or. & + (this%gwtExchange%gwtmodel2%inadv > 0 .and. & + this%gwtExchange%gwtmodel1%inadv == 0)) then + write (errmsg, '(1x,a,a,a)') 'Cannot connect GWT models in exchange ', & + trim(this%gwtExchange%name), ' because one model is configured with ADV & + &and the other one is not' + call store_error(errmsg) + end if + + if ((this%gwtExchange%gwtmodel1%indsp > 0 .and. & + this%gwtExchange%gwtmodel2%indsp == 0) .or. & + (this%gwtExchange%gwtmodel2%indsp > 0 .and. & + this%gwtExchange%gwtmodel1%indsp == 0)) then + write (errmsg, '(1x,a,a,a)') 'Cannot connect GWT models in exchange ', & + trim(this%gwtExchange%name), ' because one model is configured with DSP & + &and the other one is not' + call store_error(errmsg) + end if + + ! abort on errors + if (count_errors() > 0) then + write (errmsg, '(1x,a)') 'Errors occurred while processing exchange(s)' + call ustop() + end if + + end subroutine validateConnection !> @brief add connections to the global system for !< this connection -subroutine gwtgwtcon_ac(this, sparse) - class(GwtGwtConnectionType) :: this !< this connection - type(sparsematrix), intent(inout) :: sparse !< sparse matrix to store the connections - ! local - integer(I4B) :: ic, iglo, jglo - type(GlobalCellType) :: boundaryCell, connectedCell - - ! connections to other models - do ic = 1, this%gridConnection%nrOfBoundaryCells - boundaryCell = this%gridConnection%boundaryCells(ic)%cell - connectedCell = this%gridConnection%connectedCells(ic)%cell - iglo = boundaryCell%index + boundaryCell%model%moffset - jglo = connectedCell%index + connectedCell%model%moffset - call sparse%addconnection(iglo, jglo, 1) - call sparse%addconnection(jglo, iglo, 1) - end do - - ! and internal connections - call this%spatialcon_ac(sparse) - -end subroutine gwtgwtcon_ac - -subroutine gwtgwtcon_rp(this) - class(GwtGwtConnectionType) :: this !< the connection - - ! Call exchange rp routines - if (this%exchangeIsOwned) then - call this%gwtExchange%exg_rp() - end if - -end subroutine gwtgwtcon_rp + subroutine gwtgwtcon_ac(this, sparse) + class(GwtGwtConnectionType) :: this !< this connection + type(sparsematrix), intent(inout) :: sparse !< sparse matrix to store the connections + ! local + integer(I4B) :: ic, iglo, jglo + type(GlobalCellType) :: boundaryCell, connectedCell + + ! connections to other models + do ic = 1, this%gridConnection%nrOfBoundaryCells + boundaryCell = this%gridConnection%boundaryCells(ic)%cell + connectedCell = this%gridConnection%connectedCells(ic)%cell + iglo = boundaryCell%index + boundaryCell%model%moffset + jglo = connectedCell%index + connectedCell%model%moffset + call sparse%addconnection(iglo, jglo, 1) + call sparse%addconnection(jglo, iglo, 1) + end do + ! and internal connections + call this%spatialcon_ac(sparse) + + end subroutine gwtgwtcon_ac + + subroutine gwtgwtcon_rp(this) + class(GwtGwtConnectionType) :: this !< the connection + + ! Call exchange rp routines + if (this%exchangeIsOwned) then + call this%gwtExchange%exg_rp() + end if + + end subroutine gwtgwtcon_rp !> @brief Advance this connection !< -subroutine gwtgwtcon_ad(this) - class(GwtGwtConnectionType) :: this !< this connection - - ! copy model data into interface model - call this%syncInterfaceModel() + subroutine gwtgwtcon_ad(this) + class(GwtGwtConnectionType) :: this !< this connection - ! recalculate dispersion ellipse - if (this%gwtInterfaceModel%indsp > 0) call this%gwtInterfaceModel%dsp%dsp_ad() + ! copy model data into interface model + call this%syncInterfaceModel() - if (this%exchangeIsOwned) then - call this%gwtExchange%exg_ad() - end if + ! recalculate dispersion ellipse + if (this%gwtInterfaceModel%indsp > 0) call this%gwtInterfaceModel%dsp%dsp_ad() -end subroutine gwtgwtcon_ad + if (this%exchangeIsOwned) then + call this%gwtExchange%exg_ad() + end if + end subroutine gwtgwtcon_ad -subroutine gwtgwtcon_cf(this, kiter) - class(GwtGwtConnectionType) :: this !< the connection - integer(I4B), intent(in) :: kiter !< the iteration counter - ! local - integer(I4B) :: i + subroutine gwtgwtcon_cf(this, kiter) + class(GwtGwtConnectionType) :: this !< the connection + integer(I4B), intent(in) :: kiter !< the iteration counter + ! local + integer(I4B) :: i - ! copy model data into interface model - ! (when kiter == 1, this is already done in _ad) - if (kiter > 1) call this%syncInterfaceModel() + ! copy model data into interface model + ! (when kiter == 1, this is already done in _ad) + if (kiter > 1) call this%syncInterfaceModel() - ! reset interface system - do i = 1, this%nja - this%amat(i) = 0.0_DP - end do - do i = 1, this%neq - this%rhs(i) = 0.0_DP - end do + ! reset interface system + do i = 1, this%nja + this%amat(i) = 0.0_DP + end do + do i = 1, this%neq + this%rhs(i) = 0.0_DP + end do - call this%gwtInterfaceModel%model_cf(kiter) - -end subroutine gwtgwtcon_cf + call this%gwtInterfaceModel%model_cf(kiter) + end subroutine gwtgwtcon_cf !> @brief called during advance (*_ad), to copy the data !! from the models into the connection's placeholder arrays !< -subroutine syncInterfaceModel(this) - class(GwtGwtConnectionType) :: this !< the connection - ! local - integer(I4B) :: i, n, m, ipos, iposLoc, idx - type(ConnectionsType), pointer :: imCon !< interface model connections - type(GlobalCellType), dimension(:), pointer :: toGlobal !< map interface index to global cell - type(GlobalCellType), pointer :: boundaryCell, connectedCell - class(GwtModelType), pointer :: gwtModel - class(*), pointer :: modelPtr - - ! for readability - imCon => this%gwtInterfaceModel%dis%con - toGlobal => this%gridConnection%idxToGlobal - - ! loop over connections in interface - do n = 1, this%neq - do ipos = imCon%ia(n) + 1, imCon%ia(n+1) - 1 - m = imCon%ja(ipos) - if (associated(toGlobal(n)%model, toGlobal(m)%model)) then - ! internal connection for a model, copy from its flowja - iposLoc = getCSRIndex(toGlobal(n)%index, toGlobal(m)%index, & - toGlobal(n)%model%ia, toGlobal(n)%model%ja) - modelPtr => toGlobal(n)%model - gwtModel => CastAsGwtModel(modelPtr) - this%gwfflowja(ipos) = gwtModel%fmi%gwfflowja(iposLoc) - end if + subroutine syncInterfaceModel(this) + class(GwtGwtConnectionType) :: this !< the connection + ! local + integer(I4B) :: i, n, m, ipos, iposLoc, idx + type(ConnectionsType), pointer :: imCon !< interface model connections + type(GlobalCellType), dimension(:), pointer :: toGlobal !< map interface index to global cell + type(GlobalCellType), pointer :: boundaryCell, connectedCell + class(GwtModelType), pointer :: gwtModel + class(*), pointer :: modelPtr + + ! for readability + imCon => this%gwtInterfaceModel%dis%con + toGlobal => this%gridConnection%idxToGlobal + + ! loop over connections in interface + do n = 1, this%neq + do ipos = imCon%ia(n) + 1, imCon%ia(n + 1) - 1 + m = imCon%ja(ipos) + if (associated(toGlobal(n)%model, toGlobal(m)%model)) then + ! internal connection for a model, copy from its flowja + iposLoc = getCSRIndex(toGlobal(n)%index, toGlobal(m)%index, & + toGlobal(n)%model%ia, toGlobal(n)%model%ja) + modelPtr => toGlobal(n)%model + gwtModel => CastAsGwtModel(modelPtr) + this%gwfflowja(ipos) = gwtModel%fmi%gwfflowja(iposLoc) + end if + end do end do - end do - - ! the flowja for exchange cells - do i = 1, this%gridConnection%nrOfBoundaryCells - boundaryCell => this%gridConnection%boundaryCells(i)%cell - connectedCell => this%gridConnection%connectedCells(i)%cell - n = this%gridConnection%getInterfaceIndex(boundaryCell%index, & - boundaryCell%model) - m = this%gridConnection%getInterfaceIndex(connectedCell%index, & - connectedCell%model) - ipos = getCSRIndex(n, m, imCon%ia, imCon%ja) - this%gwfflowja(ipos) = this%exgflowja(i) * this%exgflowSign - ipos = getCSRIndex(m, n, imCon%ia, imCon%ja) - this%gwfflowja(ipos) = -this%exgflowja(i) * this%exgflowSign - end do - - ! copy concentrations - do i = 1, this%gridConnection%nrOfCells - idx = this%gridConnection%idxToGlobal(i)%index - this%x(i) = this%gridConnection%idxToGlobal(i)%model%x(idx) - this%gwtInterfaceModel%xold(i) = this%gridConnection%idxToGlobal(i)%model%xold(idx) - end do - - ! copy fmi - do i = 1, this%gridConnection%nrOfCells - idx = this%gridConnection%idxToGlobal(i)%index - modelPtr => this%gridConnection%idxToGlobal(i)%model - gwtModel => CastAsGwtModel(modelPtr) - - this%gwfsat(i) = gwtModel%fmi%gwfsat(idx) - this%gwfhead(i) = gwtModel%fmi%gwfhead(idx) - this%gwfspdis(1, i) = gwtModel%fmi%gwfspdis(1, idx) - this%gwfspdis(2, i) = gwtModel%fmi%gwfspdis(2, idx) - this%gwfspdis(3, i) = gwtModel%fmi%gwfspdis(3, idx) - end do - -end subroutine syncInterfaceModel - - -subroutine gwtgwtcon_fc(this, kiter, iasln, amatsln, rhssln, inwtflag) - class(GwtGwtConnectionType) :: this !< the connection - integer(I4B), intent(in) :: kiter !< the iteration counter - integer(I4B), dimension(:), intent(in) :: iasln !< global system's IA array - real(DP), dimension(:), intent(inout) :: amatsln !< global system matrix coefficients - real(DP), dimension(:), intent(inout) ::rhssln !< global right-hand-side - integer(I4B), optional, intent(in) :: inwtflag !< newton-raphson flag - ! local - integer(I4B) :: n, nglo, ipos - - call this%gwtInterfaceModel%model_fc(kiter, this%amat, this%nja, inwtflag) - - ! map back to solution matrix - do n = 1, this%neq - ! We only need the coefficients for our own model - ! (i.e. rows in the matrix that belong to this%owner): - if (.not. associated(this%gridConnection%idxToGlobal(n)%model, this%owner)) then - cycle - end if - - nglo = this%gridConnection%idxToGlobal(n)%index + this%gridConnection%idxToGlobal(n)%model%moffset - rhssln(nglo) = rhssln(nglo) + this%rhs(n) - - do ipos = this%ia(n), this%ia(n+1) - 1 - amatsln(this%mapIdxToSln(ipos)) = amatsln(this%mapIdxToSln(ipos)) + this%amat(ipos) + + ! the flowja for exchange cells + do i = 1, this%gridConnection%nrOfBoundaryCells + boundaryCell => this%gridConnection%boundaryCells(i)%cell + connectedCell => this%gridConnection%connectedCells(i)%cell + n = this%gridConnection%getInterfaceIndex(boundaryCell%index, & + boundaryCell%model) + m = this%gridConnection%getInterfaceIndex(connectedCell%index, & + connectedCell%model) + ipos = getCSRIndex(n, m, imCon%ia, imCon%ja) + this%gwfflowja(ipos) = this%exgflowja(i) * this%exgflowSign + ipos = getCSRIndex(m, n, imCon%ia, imCon%ja) + this%gwfflowja(ipos) = -this%exgflowja(i) * this%exgflowSign end do - end do - ! FC the movers through the exchange; we can call - ! exg_fc() directly because it only handles mover terms (unlike in GwfExchange%exg_fc) - if (this%exchangeIsOwned) then - call this%gwtExchange%exg_fc(kiter, iasln, amatsln, rhssln, inwtflag) - end if + ! copy concentrations + do i = 1, this%gridConnection%nrOfCells + idx = this%gridConnection%idxToGlobal(i)%index + this%x(i) = this%gridConnection%idxToGlobal(i)%model%x(idx) + this%gwtInterfaceModel%xold(i) = & + this%gridConnection%idxToGlobal(i)%model%xold(idx) + end do -end subroutine gwtgwtcon_fc + ! copy fmi + do i = 1, this%gridConnection%nrOfCells + idx = this%gridConnection%idxToGlobal(i)%index + modelPtr => this%gridConnection%idxToGlobal(i)%model + gwtModel => CastAsGwtModel(modelPtr) + + this%gwfsat(i) = gwtModel%fmi%gwfsat(idx) + this%gwfhead(i) = gwtModel%fmi%gwfhead(idx) + this%gwfspdis(1, i) = gwtModel%fmi%gwfspdis(1, idx) + this%gwfspdis(2, i) = gwtModel%fmi%gwfspdis(2, idx) + this%gwfspdis(3, i) = gwtModel%fmi%gwfspdis(3, idx) + end do -subroutine gwtgwtcon_cq(this, icnvg, isuppress_output, isolnid) - class(GwtGwtConnectionType) :: this !< the connection - integer(I4B), intent(inout) :: icnvg !< convergence flag - integer(I4B), intent(in) :: isuppress_output !< suppress output when =1 - integer(I4B), intent(in) :: isolnid !< solution id + end subroutine syncInterfaceModel - call this%gwtInterfaceModel%model_cq(icnvg, isuppress_output) - call this%setFlowToExchange() + subroutine gwtgwtcon_fc(this, kiter, iasln, amatsln, rhssln, inwtflag) + class(GwtGwtConnectionType) :: this !< the connection + integer(I4B), intent(in) :: kiter !< the iteration counter + integer(I4B), dimension(:), intent(in) :: iasln !< global system's IA array + real(DP), dimension(:), intent(inout) :: amatsln !< global system matrix coefficients + real(DP), dimension(:), intent(inout) ::rhssln !< global right-hand-side + integer(I4B), optional, intent(in) :: inwtflag !< newton-raphson flag + ! local + integer(I4B) :: n, nglo, ipos -end subroutine gwtgwtcon_cq + call this%gwtInterfaceModel%model_fc(kiter, this%amat, this%nja, inwtflag) - !> @brief Set the flows (flowja from interface model) to the + ! map back to solution matrix + do n = 1, this%neq + ! We only need the coefficients for our own model + ! (i.e. rows in the matrix that belong to this%owner): + if (.not. associated(this%gridConnection%idxToGlobal(n)%model, & + this%owner)) then + cycle + end if + + nglo = this%gridConnection%idxToGlobal(n)%index + & + this%gridConnection%idxToGlobal(n)%model%moffset + rhssln(nglo) = rhssln(nglo) + this%rhs(n) + + do ipos = this%ia(n), this%ia(n + 1) - 1 + amatsln(this%mapIdxToSln(ipos)) = amatsln(this%mapIdxToSln(ipos)) + & + this%amat(ipos) + end do + end do + + ! FC the movers through the exchange; we can call + ! exg_fc() directly because it only handles mover terms (unlike in GwfExchange%exg_fc) + if (this%exchangeIsOwned) then + call this%gwtExchange%exg_fc(kiter, iasln, amatsln, rhssln, inwtflag) + end if + + end subroutine gwtgwtcon_fc + + subroutine gwtgwtcon_cq(this, icnvg, isuppress_output, isolnid) + class(GwtGwtConnectionType) :: this !< the connection + integer(I4B), intent(inout) :: icnvg !< convergence flag + integer(I4B), intent(in) :: isuppress_output !< suppress output when =1 + integer(I4B), intent(in) :: isolnid !< solution id + + call this%gwtInterfaceModel%model_cq(icnvg, isuppress_output) + call this%setFlowToExchange() + + end subroutine gwtgwtcon_cq + + !> @brief Set the flows (flowja from interface model) to the !< simvals in the exchange, leaving the budget calcution in there subroutine setFlowToExchange(this) class(GwtGwtConnectionType) :: this !< this connection @@ -528,16 +535,19 @@ subroutine setFlowToExchange(this) class(GwtExchangeType), pointer :: gwtEx gwtEx => this%gwtExchange - if (this%exchangeIsOwned) then + if (this%exchangeIsOwned) then do i = 1, gwtEx%nexg gwtEx%simvals(i) = DZERO - if (gwtEx%gwtmodel1%ibound(gwtEx%nodem1(i)) /= 0 .and. & + if (gwtEx%gwtmodel1%ibound(gwtEx%nodem1(i)) /= 0 .and. & gwtEx%gwtmodel2%ibound(gwtEx%nodem2(i)) /= 0) then - nIface = this%gridConnection%getInterfaceIndex(gwtEx%nodem1(i), gwtEx%model1) - mIface = this%gridConnection%getInterfaceIndex(gwtEx%nodem2(i), gwtEx%model2) - ipos = getCSRIndex(nIface, mIface, this%gwtInterfaceModel%ia, this%gwtInterfaceModel%ja) + nIface = this%gridConnection%getInterfaceIndex(gwtEx%nodem1(i), & + gwtEx%model1) + mIface = this%gridConnection%getInterfaceIndex(gwtEx%nodem2(i), & + gwtEx%model2) + ipos = getCSRIndex(nIface, mIface, this%gwtInterfaceModel%ia, & + this%gwtInterfaceModel%ja) gwtEx%simvals(i) = this%gwtInterfaceModel%flowja(ipos) end if @@ -546,84 +556,84 @@ subroutine setFlowToExchange(this) end subroutine setFlowToExchange -subroutine gwtgwtcon_bd(this, icnvg, isuppress_output, isolnid) - use BudgetModule, only: rate_accumulator - class(GwtGwtConnectionType) :: this !< the connection - integer(I4B), intent(inout) :: icnvg !< convergence flag - integer(I4B), intent(in) :: isuppress_output !< suppress output when =1 - integer(I4B), intent(in) :: isolnid !< solution id - - ! call exchange budget routine, also calls bd - ! for movers. - if (this%exchangeIsOwned) then - call this%gwtExchange%exg_bd(icnvg, isuppress_output, isolnid) - end if - -end subroutine gwtgwtcon_bd - -subroutine gwtgwtcon_ot(this) - class(GwtGwtConnectionType) :: this !< the connection - - ! Call exg_ot() here as it handles all output processing - ! based on gwtExchange%simvals(:), which was correctly - ! filled from gwtgwtcon - if (this%exchangeIsOwned) then - call this%gwtExchange%exg_ot() - end if - -end subroutine gwtgwtcon_ot - -subroutine gwtgwtcon_da(this) - class(GwtGwtConnectionType) :: this !< the connection - ! local - logical(LGP) :: isOpen - - ! scalars - call mem_deallocate(this%iIfaceAdvScheme) - call mem_deallocate(this%iIfaceXt3d) - call mem_deallocate(this%exgflowSign) - - ! arrays - call mem_deallocate(this%gwfflowja) - call mem_deallocate(this%gwfsat) - call mem_deallocate(this%gwfhead) - call mem_deallocate(this%gwfspdis) - call mem_deallocate(this%exgflowjaGwt) - - ! interface model - call this%gwtInterfaceModel%model_da() - deallocate(this%gwtInterfaceModel) - - ! dealloc base - call this%spatialcon_da() - - inquire(this%iout, opened=isOpen) + subroutine gwtgwtcon_bd(this, icnvg, isuppress_output, isolnid) + use BudgetModule, only: rate_accumulator + class(GwtGwtConnectionType) :: this !< the connection + integer(I4B), intent(inout) :: icnvg !< convergence flag + integer(I4B), intent(in) :: isuppress_output !< suppress output when =1 + integer(I4B), intent(in) :: isolnid !< solution id + + ! call exchange budget routine, also calls bd + ! for movers. + if (this%exchangeIsOwned) then + call this%gwtExchange%exg_bd(icnvg, isuppress_output, isolnid) + end if + + end subroutine gwtgwtcon_bd + + subroutine gwtgwtcon_ot(this) + class(GwtGwtConnectionType) :: this !< the connection + + ! Call exg_ot() here as it handles all output processing + ! based on gwtExchange%simvals(:), which was correctly + ! filled from gwtgwtcon + if (this%exchangeIsOwned) then + call this%gwtExchange%exg_ot() + end if + + end subroutine gwtgwtcon_ot + + subroutine gwtgwtcon_da(this) + class(GwtGwtConnectionType) :: this !< the connection + ! local + logical(LGP) :: isOpen + + ! scalars + call mem_deallocate(this%iIfaceAdvScheme) + call mem_deallocate(this%iIfaceXt3d) + call mem_deallocate(this%exgflowSign) + + ! arrays + call mem_deallocate(this%gwfflowja) + call mem_deallocate(this%gwfsat) + call mem_deallocate(this%gwfhead) + call mem_deallocate(this%gwfspdis) + call mem_deallocate(this%exgflowjaGwt) + + ! interface model + call this%gwtInterfaceModel%model_da() + deallocate (this%gwtInterfaceModel) + + ! dealloc base + call this%spatialcon_da() + + inquire (this%iout, opened=isOpen) if (isOpen) then - close(this%iout) + close (this%iout) end if - ! we need to deallocate the exchange we own: - if (this%exchangeIsOwned) then - call this%gwtExchange%exg_da() - end if + ! we need to deallocate the exchange we own: + if (this%exchangeIsOwned) then + call this%gwtExchange%exg_da() + end if -end subroutine gwtgwtcon_da + end subroutine gwtgwtcon_da !> @brief Cast to GwtGwtConnectionType !< -function CastAsGwtGwtConnection(obj) result (res) - implicit none - class(*), pointer, intent(inout) :: obj !< object to be cast - class(GwtGwtConnectionType), pointer :: res !< the GwtGwtConnection - - res => null() - if (.not. associated(obj)) return - - select type (obj) - class is (GwtGwtConnectionType) - res => obj - end select - return -end function CastAsGwtGwtConnection - -end module \ No newline at end of file + function CastAsGwtGwtConnection(obj) result(res) + implicit none + class(*), pointer, intent(inout) :: obj !< object to be cast + class(GwtGwtConnectionType), pointer :: res !< the GwtGwtConnection + + res => null() + if (.not. associated(obj)) return + + select type (obj) + class is (GwtGwtConnectionType) + res => obj + end select + return + end function CastAsGwtGwtConnection + +end module diff --git a/src/Model/Connection/GwtInterfaceModel.f90 b/src/Model/Connection/GwtInterfaceModel.f90 index 0e6f3423cb7..f9b2bbb013a 100644 --- a/src/Model/Connection/GwtInterfaceModel.f90 +++ b/src/Model/Connection/GwtInterfaceModel.f90 @@ -1,5 +1,5 @@ module GwtInterfaceModelModule - use KindModule, only: I4B, DP + use KindModule, only: I4B, DP use MemoryManagerModule, only: mem_allocate, mem_deallocate use MemoryHelperModule, only: create_mem_path use NumericalModelModule, only: NumericalModelType @@ -18,9 +18,9 @@ module GwtInterfaceModelModule private !> The GWT Interface Model is a utility to calculate the solution's - !! exchange coefficients from the interface between a GWT model and - !! its GWT neighbors. The interface model itself will not be part - !! of the solution, it is not being solved. + !! exchange coefficients from the interface between a GWT model and + !! its GWT neighbors. The interface model itself will not be part + !! of the solution, it is not being solved. type, public, extends(GwtModelType) :: GwtInterfaceModelType integer(i4B), pointer :: iAdvScheme => null() !< the advection scheme: 0 = up, 1 = central, 2 = tvd @@ -43,197 +43,194 @@ module GwtInterfaceModelModule contains -!> @brief Create the interface model, analogously to what +!> @brief Create the interface model, analogously to what !< happens in gwt_cr -subroutine gwtifmod_cr(this, name, iout, gridConn) - class(GwtInterfaceModelType) :: this !< the GWT interface model - character(len=*), intent(in) :: name !< the interface model's name - integer(I4B), intent(in) :: iout !< the output unit - class(GridConnectionType), pointer, intent(in) :: gridConn !< the grid connection data for creating a DISU - ! local - class(*), pointer :: modelPtr - integer(I4B), target :: inobs - integer(I4B) :: adv_unit, dsp_unit - - this%memoryPath = create_mem_path(name) - call this%allocate_scalars(name) - - ! defaults - this%iAdvScheme = 0 - this%ixt3d = 0 - - this%iout = iout - this%gridConnection => gridConn - modelPtr => gridConn%model - this%owner => CastAsGwtModel(modelPtr) - - inobs = 0 - adv_unit = 0 - dsp_unit = 0 - if (this%owner%inadv > 0) then - this%inadv = huge(1_I4B) - adv_unit = huge(1_I4B) - end if - if (this%owner%indsp > 0) then - this%indsp = huge(1_I4B) - dsp_unit = huge(1_I4B) - end if - - ! create dis and packages - call disu_cr(this%dis, this%name, -1, this%iout) - call fmi_cr(this%fmi, this%name, 0, this%iout) - call adv_cr(this%adv, this%name, adv_unit, this%iout, this%fmi) - call dsp_cr(this%dsp, this%name, dsp_unit, this%iout, this%fmi) - call gwt_obs_cr(this%obs, inobs) - -end subroutine gwtifmod_cr - -subroutine allocate_scalars(this, modelname) - class(GwtInterfaceModelType) :: this !< the GWT interface model - character(len=*), intent(in) :: modelname !< the model name - - call this%GwtModelType%allocate_scalars(modelname) - - call mem_allocate(this%iAdvScheme, 'ADVSCHEME', this%memoryPath) - call mem_allocate(this%ixt3d, 'IXT3D', this%memoryPath) - -end subroutine allocate_scalars + subroutine gwtifmod_cr(this, name, iout, gridConn) + class(GwtInterfaceModelType) :: this !< the GWT interface model + character(len=*), intent(in) :: name !< the interface model's name + integer(I4B), intent(in) :: iout !< the output unit + class(GridConnectionType), pointer, intent(in) :: gridConn !< the grid connection data for creating a DISU + ! local + class(*), pointer :: modelPtr + integer(I4B), target :: inobs + integer(I4B) :: adv_unit, dsp_unit + + this%memoryPath = create_mem_path(name) + call this%allocate_scalars(name) + + ! defaults + this%iAdvScheme = 0 + this%ixt3d = 0 + + this%iout = iout + this%gridConnection => gridConn + modelPtr => gridConn%model + this%owner => CastAsGwtModel(modelPtr) + + inobs = 0 + adv_unit = 0 + dsp_unit = 0 + if (this%owner%inadv > 0) then + this%inadv = huge(1_I4B) + adv_unit = huge(1_I4B) + end if + if (this%owner%indsp > 0) then + this%indsp = huge(1_I4B) + dsp_unit = huge(1_I4B) + end if -!> @brief Define the GWT interface model -!< -subroutine gwtifmod_df(this) - class(GwtInterfaceModelType) :: this !< the GWT interface model - ! local - class(*), pointer :: disPtr - type(GwtAdvOptionsType) :: adv_options - type(GwtDspOptionsType) :: dsp_options - integer(I4B) :: i - - this%moffset = 0 - adv_options%iAdvScheme = this%iAdvScheme - dsp_options%ixt3d = this%ixt3d - - ! define DISU - disPtr => this%dis - call this%gridConnection%getDiscretization(CastAsDisuType(disPtr)) - call this%fmi%fmi_df(this%dis, 0) - - if (this%inadv > 0) then - call this%adv%adv_df(adv_options) - end if - if (this%indsp > 0) then - call this%dsp%dsp_df(this%dis, dsp_options) - end if - - ! assign or point model members to dis members - this%neq = this%dis%nodes - this%nja = this%dis%nja - this%ia => this%dis%con%ia - this%ja => this%dis%con%ja - ! - ! allocate model arrays, now that neq and nja are assigned - call this%allocate_arrays() - call mem_allocate(this%porosity, this%neq, 'POROSITY', this%memoryPath) - - do i = 1, size(this%flowja) - this%flowja = 0.0_DP - end do - do i = 1, this%neq - this%porosity = 0.0_DP - end do - -end subroutine gwtifmod_df - - -!> @brief Override allocate and read the GWT interface model and its -!! packages so that we can create stuff from memory instead of input -!< files -subroutine gwtifmod_ar(this) - class(GwtInterfaceModelType) :: this !< the GWT interface model - ! local - type(GwtDspGridDataType) :: dspGridData - - call this%fmi%fmi_ar(this%ibound) - if (this%inadv > 0) then - call this%adv%adv_ar(this%dis, this%ibound) - end if - if (this%indsp > 0) then - this%dsp%idiffc = this%owner%dsp%idiffc - this%dsp%idisp = this%owner%dsp%idisp - call dspGridData%construct(this%neq) - call this%setDspGridData(dspGridData) - call this%dsp%dsp_ar(this%ibound, this%porosity, dspGridData) - end if - -end subroutine gwtifmod_ar + ! create dis and packages + call disu_cr(this%dis, this%name, -1, this%iout) + call fmi_cr(this%fmi, this%name, 0, this%iout) + call adv_cr(this%adv, this%name, adv_unit, this%iout, this%fmi) + call dsp_cr(this%dsp, this%name, dsp_unit, this%iout, this%fmi) + call gwt_obs_cr(this%obs, inobs) + end subroutine gwtifmod_cr -!> @brief set dsp grid data from models + subroutine allocate_scalars(this, modelname) + class(GwtInterfaceModelType) :: this !< the GWT interface model + character(len=*), intent(in) :: modelname !< the model name + + call this%GwtModelType%allocate_scalars(modelname) + + call mem_allocate(this%iAdvScheme, 'ADVSCHEME', this%memoryPath) + call mem_allocate(this%ixt3d, 'IXT3D', this%memoryPath) + + end subroutine allocate_scalars + +!> @brief Define the GWT interface model !< -subroutine setDspGridData(this, gridData) - class(GwtInterfaceModelType) :: this !< the GWT interface model - type(GwtDspGridDataType) :: gridData !< the dsp grid data to be set - ! local - integer(I4B) :: i, idx - class(GwtModelType), pointer :: gwtModel - class(*), pointer :: modelPtr - - do i = 1, this%neq - modelPtr => this%gridConnection%idxToGlobal(i)%model - gwtModel => CastAsGwtModel(modelPtr) - idx = this%gridConnection%idxToGlobal(i)%index - - if (this%dsp%idiffc > 0) then - gridData%diffc(i) = gwtModel%dsp%diffc(idx) + subroutine gwtifmod_df(this) + class(GwtInterfaceModelType) :: this !< the GWT interface model + ! local + class(*), pointer :: disPtr + type(GwtAdvOptionsType) :: adv_options + type(GwtDspOptionsType) :: dsp_options + integer(I4B) :: i + + this%moffset = 0 + adv_options%iAdvScheme = this%iAdvScheme + dsp_options%ixt3d = this%ixt3d + + ! define DISU + disPtr => this%dis + call this%gridConnection%getDiscretization(CastAsDisuType(disPtr)) + call this%fmi%fmi_df(this%dis, 0) + + if (this%inadv > 0) then + call this%adv%adv_df(adv_options) + end if + if (this%indsp > 0) then + call this%dsp%dsp_df(this%dis, dsp_options) + end if + + ! assign or point model members to dis members + this%neq = this%dis%nodes + this%nja = this%dis%nja + this%ia => this%dis%con%ia + this%ja => this%dis%con%ja + ! + ! allocate model arrays, now that neq and nja are assigned + call this%allocate_arrays() + call mem_allocate(this%porosity, this%neq, 'POROSITY', this%memoryPath) + + do i = 1, size(this%flowja) + this%flowja = 0.0_DP + end do + do i = 1, this%neq + this%porosity = 0.0_DP + end do + + end subroutine gwtifmod_df + +!> @brief Override allocate and read the GWT interface model and its +!! packages so that we can create stuff from memory instead of input +!< files + subroutine gwtifmod_ar(this) + class(GwtInterfaceModelType) :: this !< the GWT interface model + ! local + type(GwtDspGridDataType) :: dspGridData + + call this%fmi%fmi_ar(this%ibound) + if (this%inadv > 0) then + call this%adv%adv_ar(this%dis, this%ibound) end if - if (this%dsp%idisp > 0) then - gridData%alh(i) = gwtModel%dsp%alh(idx) - gridData%alv(i) = gwtModel%dsp%alv(idx) - gridData%ath1(i) = gwtModel%dsp%ath1(idx) - gridData%ath2(i) = gwtModel%dsp%ath2(idx) - gridData%atv(i) = gwtModel%dsp%atv(idx) + if (this%indsp > 0) then + this%dsp%idiffc = this%owner%dsp%idiffc + this%dsp%idisp = this%owner%dsp%idisp + call dspGridData%construct(this%neq) + call this%setDspGridData(dspGridData) + call this%dsp%dsp_ar(this%ibound, this%porosity, dspGridData) end if - end do + end subroutine gwtifmod_ar -end subroutine setDspGridData +!> @brief set dsp grid data from models +!< + subroutine setDspGridData(this, gridData) + class(GwtInterfaceModelType) :: this !< the GWT interface model + type(GwtDspGridDataType) :: gridData !< the dsp grid data to be set + ! local + integer(I4B) :: i, idx + class(GwtModelType), pointer :: gwtModel + class(*), pointer :: modelPtr + + do i = 1, this%neq + modelPtr => this%gridConnection%idxToGlobal(i)%model + gwtModel => CastAsGwtModel(modelPtr) + idx = this%gridConnection%idxToGlobal(i)%index + + if (this%dsp%idiffc > 0) then + gridData%diffc(i) = gwtModel%dsp%diffc(idx) + end if + if (this%dsp%idisp > 0) then + gridData%alh(i) = gwtModel%dsp%alh(idx) + gridData%alv(i) = gwtModel%dsp%alv(idx) + gridData%ath1(i) = gwtModel%dsp%ath1(idx) + gridData%ath2(i) = gwtModel%dsp%ath2(idx) + gridData%atv(i) = gwtModel%dsp%atv(idx) + end if + + end do + + end subroutine setDspGridData !> @brief Clean up resources !< -subroutine gwtifmod_da(this) - class(GwtInterfaceModelType) :: this !< the GWT interface model - - ! this - call mem_deallocate(this%iAdvScheme) - call mem_deallocate(this%ixt3d) - call mem_deallocate(this%porosity) - - ! gwt packages - call this%dis%dis_da() - call this%fmi%fmi_da() - call this%adv%adv_da() - call this%dsp%dsp_da() - - deallocate(this%dis) - deallocate(this%fmi) - deallocate(this%adv) - deallocate(this%dsp) - - ! gwt scalars - call mem_deallocate(this%inic) - call mem_deallocate(this%infmi) - call mem_deallocate(this%inadv) - call mem_deallocate(this%indsp) - call mem_deallocate(this%inssm) - call mem_deallocate(this%inmst) - call mem_deallocate(this%inmvt) - call mem_deallocate(this%inoc) - call mem_deallocate(this%inobs) - - ! base - call this%NumericalModelType%model_da() - -end subroutine gwtifmod_da - - -end module GwtInterfaceModelModule \ No newline at end of file + subroutine gwtifmod_da(this) + class(GwtInterfaceModelType) :: this !< the GWT interface model + + ! this + call mem_deallocate(this%iAdvScheme) + call mem_deallocate(this%ixt3d) + call mem_deallocate(this%porosity) + + ! gwt packages + call this%dis%dis_da() + call this%fmi%fmi_da() + call this%adv%adv_da() + call this%dsp%dsp_da() + + deallocate (this%dis) + deallocate (this%fmi) + deallocate (this%adv) + deallocate (this%dsp) + + ! gwt scalars + call mem_deallocate(this%inic) + call mem_deallocate(this%infmi) + call mem_deallocate(this%inadv) + call mem_deallocate(this%indsp) + call mem_deallocate(this%inssm) + call mem_deallocate(this%inmst) + call mem_deallocate(this%inmvt) + call mem_deallocate(this%inoc) + call mem_deallocate(this%inobs) + + ! base + call this%NumericalModelType%model_da() + + end subroutine gwtifmod_da + +end module GwtInterfaceModelModule diff --git a/src/Model/Connection/SpatialModelConnection.f90 b/src/Model/Connection/SpatialModelConnection.f90 index 19f1ba3c486..5e766bc9ced 100644 --- a/src/Model/Connection/SpatialModelConnection.f90 +++ b/src/Model/Connection/SpatialModelConnection.f90 @@ -1,6 +1,6 @@ module SpatialModelConnectionModule use KindModule, only: I4B, DP, LGP - use SparseModule, only:sparsematrix + use SparseModule, only: sparsematrix use ConnectionsModule, only: ConnectionsType use CsrUtilsModule, only: getCSRIndex use SimModule, only: ustop @@ -11,18 +11,18 @@ module SpatialModelConnectionModule use MemoryHelperModule, only: create_mem_path use GridConnectionModule, only: GridConnectionType use ListModule, only: ListType - + implicit none private public :: CastAsSpatialModelConnectionClass public :: AddSpatialModelConnectionToList public :: GetSpatialModelConnectionFromList - !> Class to manage spatial connection of a model to one - !! or more models of the same type. Spatial connection here - !! means that the model domains (spatial discretization) are + !> Class to manage spatial connection of a model to one + !! or more models of the same type. Spatial connection here + !! means that the model domains (spatial discretization) are !! adjacent and connected via DisConnExchangeType object(s). - !! The connection itself is a Numerical Exchange as well, + !! The connection itself is a Numerical Exchange as well, !! and part of a Numerical Solution providing the amat and rhs !< values for the exchange. type, public, extends(NumericalExchangeType) :: SpatialModelConnectionType @@ -37,8 +37,7 @@ module SpatialModelConnectionModule !! default = 1, xt3d = 2, ... integer(I4B), pointer :: exchangeStencilDepth => null() !< size of the computational stencil at the interface !! default = 1, xt3d = 2, ... - - + ! The following variables are equivalent to those in Numerical Solution: integer(I4B), pointer :: neq => null() !< nr. of equations in matrix system integer(I4B), pointer :: nja => null() !< nr. of nonzero matrix elements @@ -48,26 +47,26 @@ module SpatialModelConnectionModule real(DP), dimension(:), pointer, contiguous :: rhs => null() !< rhs of interface system real(DP), dimension(:), pointer, contiguous :: x => null() !< dependent variable of interface system integer(I4B), dimension(:), pointer, contiguous :: active => null() !< cell status (c.f. ibound) of interface system - + ! these are not in the memory manager class(GridConnectionType), pointer :: gridConnection => null() !< facility to build the interface grid connection structure integer(I4B), dimension(:), pointer :: mapIdxToSln => null() !< mapping between interface matrix and the solution matrix - + contains - + ! public procedure, pass(this) :: spatialConnection_ctor generic :: construct => spatialConnection_ctor ! partly overriding NumericalExchangeType: - procedure :: exg_df => spatialcon_df + procedure :: exg_df => spatialcon_df procedure :: exg_ar => spatialcon_ar procedure :: exg_ac => spatialcon_ac procedure :: exg_mc => spatialcon_mc procedure :: exg_da => spatialcon_da - + ! protected - procedure, pass(this) :: spatialcon_df + procedure, pass(this) :: spatialcon_df procedure, pass(this) :: spatialcon_ar procedure, pass(this) :: spatialcon_ac procedure, pass(this) :: spatialcon_da @@ -78,16 +77,16 @@ module SpatialModelConnectionModule ! private procedure, private, pass(this) :: setupGridConnection procedure, private, pass(this) :: setExchangeConnections - procedure, private, pass(this) :: getNrOfConnections + procedure, private, pass(this) :: getNrOfConnections procedure, private, pass(this) :: allocateScalars - procedure, private, pass(this) :: allocateArrays - procedure, private, pass(this) :: createCoefficientMatrix + procedure, private, pass(this) :: allocateArrays + procedure, private, pass(this) :: createCoefficientMatrix procedure, private, pass(this) :: maskOwnerConnections - + end type SpatialModelConnectionType contains ! module procedures - + !> @brief Construct the spatial connection base !! !! This constructor is typically called from a derived class. @@ -95,17 +94,17 @@ module SpatialModelConnectionModule subroutine spatialConnection_ctor(this, model, exchange, name) class(SpatialModelConnectionType) :: this !< the connection class(NumericalModelType), intent(in), pointer :: model !< the model that owns the connection - class(DisConnExchangeType), intent(in), pointer :: exchange !< the primary exchange from which + class(DisConnExchangeType), intent(in), pointer :: exchange !< the primary exchange from which !! the connection is created character(len=*), intent(in) :: name !< the connection name (for memory management mostly) - + this%name = name this%memoryPath = create_mem_path(this%name) this%owner => model this%primaryExchange => exchange - allocate(this%gridConnection) + allocate (this%gridConnection) call this%allocateScalars() this%internalStencilDepth = 1 @@ -114,28 +113,29 @@ subroutine spatialConnection_ctor(this, model, exchange, name) ! this should be set in derived ctor this%interfaceModel => null() - + end subroutine spatialConnection_ctor - - + !> @brief Define this connection, mostly sets up the grid !< connection, allocates arrays, and links x,rhs, and ibound subroutine spatialcon_df(this) class(SpatialModelConnectionType) :: this !< this connection - + ! create the grid connection data structure this%nrOfConnections = this%getNrOfConnections() - call this%gridConnection%construct(this%owner, this%nrOfConnections, this%name) + call this%gridConnection%construct(this%owner, & + this%nrOfConnections, & + this%name) this%gridConnection%internalStencilDepth = this%internalStencilDepth this%gridConnection%exchangeStencilDepth = this%exchangeStencilDepth call this%setupGridConnection() - + this%neq = this%gridConnection%nrOfCells call this%allocateArrays() - + end subroutine spatialcon_df - !> @brief Allocate the connection, + !> @brief Allocate the connection, !< subroutine spatialcon_ar(this) class(SpatialModelConnectionType) :: this !< this connection @@ -143,12 +143,12 @@ subroutine spatialcon_ar(this) integer(I4B) :: icell, idx, localIdx class(GridConnectionType), pointer :: gc class(NumericalModelType), pointer :: model - + ! init x and ibound with model data gc => this%gridConnection do icell = 1, gc%nrOfCells idx = gc%idxToGlobal(icell)%index - model => gc%idxToGlobal(icell)%model + model => gc%idxToGlobal(icell)%model this%interfaceModel%x(icell) = model%x(idx) this%interfaceModel%ibound(icell) = model%ibound(idx) end do @@ -156,28 +156,34 @@ subroutine spatialcon_ar(this) ! fill mapping to global index (which can be ! done now because moffset is set in sln_df) do localIdx = 1, gc%nrOfCells - gc%idxToGlobalIdx(localIdx) = gc%idxToGlobal(localIdx)%index + & + gc%idxToGlobalIdx(localIdx) = gc%idxToGlobal(localIdx)%index + & gc%idxToGlobal(localIdx)%model%moffset end do end subroutine spatialcon_ar - !> @brief set model pointers to connection + !> @brief set model pointers to connection !< subroutine spatialcon_setmodelptrs(this) class(SpatialModelConnectionType) :: this !< this connection - + ! point x, ibound, and rhs to connection this%interfaceModel%x => this%x - call mem_checkin(this%interfaceModel%x, 'X', this%interfaceModel%memoryPath, 'X', this%memoryPath) + call mem_checkin(this%interfaceModel%x, 'X', & + this%interfaceModel%memoryPath, 'X', & + this%memoryPath) this%interfaceModel%rhs => this%rhs - call mem_checkin(this%interfaceModel%rhs, 'RHS', this%interfaceModel%memoryPath, 'RHS', this%memoryPath) + call mem_checkin(this%interfaceModel%rhs, 'RHS', & + this%interfaceModel%memoryPath, 'RHS', & + this%memoryPath) this%interfaceModel%ibound => this%active - call mem_checkin(this%interfaceModel%ibound, 'IBOUND', this%interfaceModel%memoryPath, 'IBOUND', this%memoryPath) + call mem_checkin(this%interfaceModel%ibound, 'IBOUND', & + this%interfaceModel%memoryPath, 'IBOUND', & + this%memoryPath) end subroutine spatialcon_setmodelptrs - !> @brief map interface model connections to our sparse matrix, + !> @brief map interface model connections to our sparse matrix, !< analogously to what happens in sln_connect. subroutine spatialcon_connect(this) class(SpatialModelConnectionType) :: this !< this connection @@ -186,21 +192,20 @@ subroutine spatialcon_connect(this) call sparse%init(this%neq, this%neq, 7) call this%interfaceModel%model_ac(sparse) - + ! create amat from sparse call this%createCoefficientMatrix(sparse) call sparse%destroy() - + ! map connections call this%interfaceModel%model_mc(this%ia, this%ja) call this%maskOwnerConnections() end subroutine spatialcon_connect - !> @brief Mask the owner's connections !! - !! Determine which connections are handled by the interface model + !! Determine which connections are handled by the interface model !! (using the connections object in its discretization) and !< set their mask to zero for the owning model. subroutine maskOwnerConnections(this) @@ -209,77 +214,82 @@ subroutine maskOwnerConnections(this) ! local integer(I4B) :: ipos, n, m, nloc, mloc, csrIdx type(ConnectionsType), pointer :: conn - + ! set the mask on connections that are calculated by the interface model conn => this%interfaceModel%dis%con do n = 1, conn%nodes ! only for connections internal to the owning model - if (.not. associated(this%gridConnection%idxToGlobal(n)%model, this%owner)) then + if (.not. associated(this%gridConnection%idxToGlobal(n)%model, & + this%owner)) then cycle end if nloc = this%gridConnection%idxToGlobal(n)%index - + do ipos = conn%ia(n) + 1, conn%ia(n + 1) - 1 m = conn%ja(ipos) - if (.not. associated(this%gridConnection%idxToGlobal(m)%model, this%owner)) then - cycle + if (.not. associated(this%gridConnection%idxToGlobal(m)%model, & + this%owner)) then + cycle end if mloc = this%gridConnection%idxToGlobal(m)%index - + if (conn%mask(ipos) > 0) then ! calculated by interface model, set local model's mask to zero - csrIdx = getCSRIndex(nloc, mloc, this%owner%ia, this%owner%ja) + csrIdx = getCSRIndex(nloc, mloc, this%owner%ia, this%owner%ja) if (csrIdx == -1) then - ! this can only happen with periodic boundary conditions, + ! this can only happen with periodic boundary conditions, ! then there is no need to set the mask if (this%gridConnection%isPeriodic(nloc, mloc)) cycle - - write(*,*) 'Error: cannot find cell connection in global system' + + write (*, *) 'Error: cannot find cell connection in global system' call ustop() - end if + end if if (this%owner%dis%con%mask(csrIdx) > 0) then call this%owner%dis%con%set_mask(csrIdx, 0) else ! edge case, someone will be calculating this connection ! so we ignore it here (TODO_MJR: add name) - write(*,*) 'Debug: overlap detected, ignoring connection ', & - nloc, ':', mloc, ' for model ', trim(this%owner%name), & - ' in Exchange ???' + write (*, *) 'Debug: overlap detected, ignoring connection ', & + nloc, ':', mloc, ' for model ', trim(this%owner%name), & + ' in Exchange ???' call conn%set_mask(ipos, 0) end if end if end do end do - + end subroutine maskOwnerConnections !> @brief Add connections, handled by the interface model, !< to the global system's sparse - subroutine spatialcon_ac(this, sparse) + subroutine spatialcon_ac(this, sparse) class(SpatialModelConnectionType) :: this !< this connection type(sparsematrix), intent(inout) :: sparse !< sparse matrix to store the connections ! local integer(I4B) :: n, m, ipos integer(I4B) :: nglo, mglo - + do n = 1, this%neq - if (.not. associated(this%gridConnection%idxToGlobal(n)%model, this%owner)) then + if (.not. associated(this%gridConnection%idxToGlobal(n)%model, & + this%owner)) then ! only add connections for own model to global matrix cycle end if - nglo = this%gridConnection%idxToGlobal(n)%index + this%gridConnection%idxToGlobal(n)%model%moffset - do ipos = this%ia(n) + 1, this%ia(n+1) - 1 + nglo = this%gridConnection%idxToGlobal(n)%index + & + this%gridConnection%idxToGlobal(n)%model%moffset + do ipos = this%ia(n) + 1, this%ia(n + 1) - 1 m = this%ja(ipos) - mglo = this%gridConnection%idxToGlobal(m)%index + this%gridConnection%idxToGlobal(m)%model%moffset - + mglo = this%gridConnection%idxToGlobal(m)%index + & + this%gridConnection%idxToGlobal(m)%model%moffset + call sparse%addconnection(nglo, mglo, 1) end do end do - + end subroutine spatialcon_ac - !> @brief Creates the mapping from the local system + !> @brief Creates the mapping from the local system !< matrix to the global one subroutine spatialcon_mc(this, iasln, jasln) use SimModule, only: ustop @@ -289,33 +299,36 @@ subroutine spatialcon_mc(this, iasln, jasln) ! local integer(I4B) :: m, n, mglo, nglo, ipos, csrIdx logical(LGP) :: isOwnedConnection - - allocate(this%mapIdxToSln(this%nja)) - - do n = 1, this%neq - isOwnedConnection = associated(this%gridConnection%idxToGlobal(n)%model, this%owner) - do ipos = this%ia(n), this%ia(n+1)-1 - m = this%ja(ipos) - nglo = this%gridConnection%idxToGlobal(n)%index + this%gridConnection%idxToGlobal(n)%model%moffset - mglo = this%gridConnection%idxToGlobal(m)%index + this%gridConnection%idxToGlobal(m)%model%moffset + + allocate (this%mapIdxToSln(this%nja)) + + do n = 1, this%neq + isOwnedConnection = associated(this%gridConnection%idxToGlobal(n)%model, & + this%owner) + do ipos = this%ia(n), this%ia(n + 1) - 1 + m = this%ja(ipos) + nglo = this%gridConnection%idxToGlobal(n)%index + & + this%gridConnection%idxToGlobal(n)%model%moffset + mglo = this%gridConnection%idxToGlobal(m)%index + & + this%gridConnection%idxToGlobal(m)%model%moffset csrIdx = getCSRIndex(nglo, mglo, iasln, jasln) if (csrIdx == -1 .and. isOwnedConnection) then ! this should not be possible - write(*,*) 'Error: cannot find cell connection in global system' + write (*, *) 'Error: cannot find cell connection in global system' call ustop() end if - + this%mapIdxToSln(ipos) = csrIdx end do end do - + end subroutine spatialcon_mc - + !> @brief Deallocation !< subroutine spatialcon_da(this) class(SpatialModelConnectionType) :: this !< this connection - + call mem_deallocate(this%neq) call mem_deallocate(this%nja) call mem_deallocate(this%internalStencilDepth) @@ -325,41 +338,41 @@ subroutine spatialcon_da(this) call mem_deallocate(this%ia) call mem_deallocate(this%ja) call mem_deallocate(this%amat) - + call mem_deallocate(this%x) call mem_deallocate(this%rhs) call mem_deallocate(this%active) - + call this%gridConnection%destroy() - deallocate(this%gridConnection) - deallocate(this%mapIdxToSln) - + deallocate (this%gridConnection) + deallocate (this%mapIdxToSln) + end subroutine spatialcon_da - + !> @brief Set up the grid connection !! !! This works in three steps: !! 1. set the primary connections - !! 2. create the topology of connected models, finding + !! 2. create the topology of connected models, finding !! neighbors of neighboring models when required !! 3. extend the interface grid, using that information !< subroutine setupGridConnection(this) class(SpatialModelConnectionType) :: this !< this connection ! local - + ! set boundary cells call this%setExchangeConnections() - + ! create topology of models - call this%gridConnection%findModelNeighbors(this%globalExchanges, & + call this%gridConnection%findModelNeighbors(this%globalExchanges, & this%exchangeStencilDepth) - - ! now scan for nbr-of-nbrs and create final data structures + + ! now scan for nbr-of-nbrs and create final data structures call this%gridConnection%extendConnection() - + end subroutine setupGridConnection - + !> @brief Set the primary connections from the exchange data !< subroutine setExchangeConnections(this) @@ -367,31 +380,30 @@ subroutine setExchangeConnections(this) ! local integer(I4B) :: iconn type(DisConnExchangeType), pointer :: connEx - + ! set boundary cells connEx => this%primaryExchange - do iconn=1, connEx%nexg + do iconn = 1, connEx%nexg call this%gridConnection%connectCell(connEx%nodem1(iconn), connEx%model1, & connEx%nodem2(iconn), connEx%model2) end do - + end subroutine setExchangeConnections - !> @brief Allocation of scalars !< subroutine allocateScalars(this) use MemoryManagerModule, only: mem_allocate class(SpatialModelConnectionType) :: this !< this connection - + call mem_allocate(this%neq, 'NEQ', this%memoryPath) call mem_allocate(this%nja, 'NJA', this%memoryPath) call mem_allocate(this%internalStencilDepth, 'INTSTDEPTH', this%memoryPath) call mem_allocate(this%exchangeStencilDepth, 'EXGSTDEPTH', this%memoryPath) call mem_allocate(this%nrOfConnections, 'NROFCONNS', this%memoryPath) - + end subroutine allocateScalars - + !> @brief Allocation of arrays !< subroutine allocateArrays(this) @@ -400,31 +412,31 @@ subroutine allocateArrays(this) class(SpatialModelConnectionType) :: this !< this connection ! local integer(I4B) :: i - + call mem_allocate(this%x, this%neq, 'X', this%memoryPath) call mem_allocate(this%rhs, this%neq, 'RHS', this%memoryPath) call mem_allocate(this%active, this%neq, 'IACTIVE', this%memoryPath) - + ! c.f. NumericalSolution do i = 1, this%neq this%x(i) = DZERO this%active(i) = 1 ! default is active this%rhs(i) = DZERO - enddo - + end do + end subroutine allocateArrays - + !> @brief Returns total nr. of primary connections !< function getNrOfConnections(this) result(nrConns) class(SpatialModelConnectionType) :: this !< this connection - integer(I4B) :: nrConns + integer(I4B) :: nrConns !local - + nrConns = this%primaryExchange%nexg - + end function getNrOfConnections - + !> @brief Create connection's matrix (ia,ja,amat) from sparse !< subroutine createCoefficientMatrix(this, sparse) @@ -433,7 +445,7 @@ subroutine createCoefficientMatrix(this, sparse) type(sparsematrix), intent(inout) :: sparse !< the sparse matrix with the cell connections ! local integer(I4B) :: ierror - + this%nja = sparse%nnz call mem_allocate(this%ia, this%neq + 1, 'IA', this%memoryPath) call mem_allocate(this%ja, this%nja, 'JA', this%memoryPath) @@ -443,10 +455,10 @@ subroutine createCoefficientMatrix(this, sparse) call sparse%filliaja(this%ia, this%ja, ierror) if (ierror /= 0) then - write(*,*) 'Error: cannot fill ia/ja for model connection' + write (*, *) 'Error: cannot fill ia/ja for model connection' call ustop() end if - + end subroutine createCoefficientMatrix !> @brief Validate this connection @@ -459,31 +471,30 @@ subroutine validateConnection(this) class(DisConnExchangeType), pointer :: conEx => null() conEx => this%primaryExchange - if (conEx%ixt3d > 0) then + if (conEx%ixt3d > 0) then ! if XT3D, we need these angles: if (conEx%model1%dis%con%ianglex == 0) then - write(errmsg, '(1x,a,a,a,a,a)') 'XT3D configured on the exchange ', & - trim(conEx%name), ' but the discretization in model ', & - trim(conEx%model1%name), ' has no ANGLDEGX specified' + write (errmsg, '(1x,a,a,a,a,a)') 'XT3D configured on the exchange ', & + trim(conEx%name), ' but the discretization in model ', & + trim(conEx%model1%name), ' has no ANGLDEGX specified' call store_error(errmsg) end if if (conEx%model2%dis%con%ianglex == 0) then - write(errmsg, '(1x,a,a,a,a,a)') 'XT3D configured on the exchange ', & - trim(conEx%name), ' but the discretization in model ', & - trim(conEx%model2%name), ' has no ANGLDEGX specified' + write (errmsg, '(1x,a,a,a,a,a)') 'XT3D configured on the exchange ', & + trim(conEx%name), ' but the discretization in model ', & + trim(conEx%model2%name), ' has no ANGLDEGX specified' call store_error(errmsg) end if end if end subroutine validateConnection - !> @brief Cast to SpatialModelConnectionType !< - function CastAsSpatialModelConnectionClass(obj) result (res) + function CastAsSpatialModelConnectionClass(obj) result(res) implicit none class(*), pointer, intent(inout) :: obj !< object to be cast - class(SpatialModelConnectionType), pointer :: res !< the instance of SpatialModelConnectionType + class(SpatialModelConnectionType), pointer :: res !< the instance of SpatialModelConnectionType ! res => null() if (.not. associated(obj)) return @@ -500,7 +511,7 @@ end function CastAsSpatialModelConnectionClass subroutine AddSpatialModelConnectionToList(list, conn) implicit none ! -- dummy - type(ListType), intent(inout) :: list !< the list + type(ListType), intent(inout) :: list !< the list class(SpatialModelConnectionType), pointer, intent(in) :: conn !< the connection ! -- local class(*), pointer :: obj @@ -517,13 +528,13 @@ function GetSpatialModelConnectionFromList(list, idx) result(res) type(ListType), intent(inout) :: list !< the list integer(I4B), intent(in) :: idx !< the index of the connection class(SpatialModelConnectionType), pointer :: res !< the returned connection - + ! local class(*), pointer :: obj obj => list%GetItem(idx) res => CastAsSpatialModelConnectionClass(obj) ! return - end function GetSpatialModelConnectionFromList - + end function GetSpatialModelConnectionFromList + end module SpatialModelConnectionModule