GweGweConnection.f90

Go to the documentation of this file.

module gwegweconnectionmodule
  use kindmodule, only: i4b, dp, lgp
  use constantsmodule, only: linelength, lencomponentname, dzero, lenbudtxt
  use csrutilsmodule, only: getcsrindex
  use simmodule, only: ustop
  use memorymanagermodule, only: mem_allocate, mem_deallocate, mem_checkin
  use spatialmodelconnectionmodule
  use numericalmodelmodule
  use gwemodule
  use disconnexchangemodule
  use gwegweexchangemodule
  use gweinterfacemodelmodule
  use sparsemodule, only: sparsematrix
  use connectionsmodule, only: connectionstype
  use cellwithnbrsmodule, only: globalcelltype
  use distvariablemodule
  use simstagesmodule
  use matrixbasemodule
 
  implicit none
  private
 
  public :: castasgwegweconnection
 
  !> Connects a GWE model to other GWE models in space. Derives
  !! from NumericalExchangeType so the solution can use it to
  !! fetch the coefficients for this connection.
  !<
  type, public, extends(spatialmodelconnectiontype) :: gwegweconnectiontype
 
    class(gwemodeltype), pointer :: gwemodel => null() !< the model for which this connection exists
    class(gweexchangetype), pointer :: gweexchange => null() !< the primary exchange, cast to GWE-GWE
    class(gweinterfacemodeltype), pointer :: gweinterfacemodel => null() !< the interface model
    integer(I4B), pointer :: iifaceadvscheme => null() !< the advection scheme at the interface:
                                                       !! 0 = upstream, 1 = central, 2 = TVD
    integer(I4B), pointer :: iifacext3d => null() !< XT3D in the interface CND package: 0 = no, 1 = lhs, 2 = rhs
    integer(I4B), pointer :: exgflowsign => null() !< indicates the flow direction of exgflowja
    real(dp), dimension(:), pointer, contiguous :: exgflowjagwe => null() !< gwe-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 :: conc => null() !< pointer to concentration array
    integer(I4B), dimension(:), pointer, contiguous :: icbound => null() !< store pointer to gwe ibound array
 
    integer(I4B) :: iout = 0 !< the list file for the interface model
 
  contains
 
    procedure, pass(this) :: gwegweconnection_ctor
    generic, public :: construct => gwegweconnection_ctor
 
    procedure :: exg_ar => gwegwecon_ar
    procedure :: exg_df => gwegwecon_df
    procedure :: exg_rp => gwegwecon_rp
    procedure :: exg_ad => gwegwecon_ad
    procedure :: exg_fc => gwegwecon_fc
    procedure :: exg_da => gwegwecon_da
    procedure :: exg_cq => gwegwecon_cq
    procedure :: exg_bd => gwegwecon_bd
    procedure :: exg_ot => gwegwecon_ot
 
    ! overriding 'protected'
    procedure :: validateconnection
 
    ! local stuff
    procedure, private :: allocate_scalars
    procedure, private :: allocate_arrays
    procedure, private :: cfg_dist_vars
    procedure, private :: setgridextent
    procedure, private :: setflowtoexchange
 
  end type gwegweconnectiontype
 
contains
 
  !> @brief Basic construction of the connection
  !<
  subroutine gwegweconnection_ctor(this, model, gweEx)
    ! -- modules
    use inputoutputmodule, only: openfile
    ! -- dummy
    class(gwegweconnectiontype) :: this !< the connection
    class(numericalmodeltype), pointer :: model !< the model owning this connection,
                                                !! this must be a GweModelType
    class(disconnexchangetype), pointer :: gweEx !< the GWE-GWE 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%gweModel => castasgwemodel(objptr)
    objptr => gweex
    this%gweExchange => castasgweexchange(objptr)
    !
    if (gweex%v_model1%is_local .and. gweex%v_model2%is_local) then
      this%owns_exchange = associated(model, gweex%model1)
    else
      this%owns_exchange = .true.
    end if
    !
    if (gweex%v_model1 == model) then
      write (name, '(a,i0)') 'GWECON1_', gweex%id
    else
      write (name, '(a,i0)') 'GWECON2_', gweex%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 GWE-GWE connection for model ', &
        trim(this%gweModel%name), 'from exchange ', &
        trim(gweex%name)
    end if
    !
    ! -- First call base constructor
    call this%SpatialModelConnectionType%spatialConnection_ctor(model, &
                                                                gweex, &
                                                                name)
    !
    call this%allocate_scalars()
    this%typename = 'GWE-GWE'
    this%iIfaceAdvScheme = 0
    this%iIfaceXt3d = 0
    this%exgflowSign = 1
    !
    allocate (this%gweInterfaceModel)
    this%interface_model => this%gweInterfaceModel
    !
  end subroutine gwegweconnection_ctor
 
  !> @brief Allocate scalar variables for this connection
  !<
  subroutine allocate_scalars(this)
    ! -- dummy
    class(gwegweconnectiontype) :: 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)
    !
  end subroutine allocate_scalars
 
  !> @brief define the GWE-GWE connection
  !<
  subroutine gwegwecon_df(this)
    ! -- dummy
    class(gwegweconnectiontype) :: this !< the connection
    ! -- local
    character(len=LENCOMPONENTNAME) :: imName
 
    ! -- Determine advection scheme (the GWE-GWE exchange
    !    has been read at this point)
    this%iIfaceAdvScheme = this%gweExchange%iAdvScheme
    !
    ! -- Determine xt3d setting on interface
    this%iIfaceXt3d = this%gweExchange%ixt3d
 
    ! -- Turn off when off in the owning model
    if (this%gweModel%incnd > 0) then
      this%iIfaceXt3d = this%gweModel%cnd%ixt3d
    end if
 
    ! -- 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%prim_exchange%v_model1 == this%owner) then
      write (imname, '(a,i0)') 'GWEIM1_', this%gweExchange%id
    else
      write (imname, '(a,i0)') 'GWEIM2_', this%gweExchange%id
    end if
    call this%gweInterfaceModel%gweifmod_cr(imname, &
                                            this%iout, &
                                            this%ig_builder)
    call this%gweInterfaceModel%set_idsoln(this%gweModel%idsoln)
    this%gweInterfaceModel%iAdvScheme = this%iIfaceAdvScheme
    this%gweInterfaceModel%ixt3d = this%iIfaceXt3d
    call this%gweInterfaceModel%model_df()
 
    call this%cfg_dist_vars()
 
    call this%allocate_arrays()
    call this%gweInterfaceModel%allocate_fmi()
 
    ! -- Connect X, RHS, IBOUND, and flowja
    call this%spatialcon_setmodelptrs()
 
    ! -- Connect pointers (used by BUY)
    this%conc => this%gweInterfaceModel%x
    this%icbound => this%gweInterfaceModel%ibound
 
    ! -- Add connections from the interface model to solution matrix
    call this%spatialcon_connect()
 
  end subroutine gwegwecon_df
 
  !> @brief Configure distributed variables for this interface model
  !<
  subroutine cfg_dist_vars(this)
    ! -- dummy
    class(gwegweconnectiontype) :: this !< the connection
    !
    call this%cfg_dv('X', '', sync_nds, &
                     (/stg_bfr_con_ar, stg_bfr_exg_ad, stg_bfr_exg_cf/))
    call this%cfg_dv('IBOUND', '', sync_nds, (/stg_bfr_con_ar/))
    call this%cfg_dv('TOP', 'DIS', sync_nds, (/stg_bfr_con_ar/))
    call this%cfg_dv('BOT', 'DIS', sync_nds, (/stg_bfr_con_ar/))
    call this%cfg_dv('AREA', 'DIS', sync_nds, (/stg_bfr_con_ar/))
    !
    if (this%gweInterfaceModel%cnd%idisp > 0) then
      call this%cfg_dv('ALH', 'CND', sync_nds, (/stg_bfr_con_ar/))
      call this%cfg_dv('ALV', 'CND', sync_nds, (/stg_bfr_con_ar/))
      call this%cfg_dv('ATH1', 'CND', sync_nds, (/stg_bfr_con_ar/))
      call this%cfg_dv('ATH2', 'CND', sync_nds, (/stg_bfr_con_ar/))
      call this%cfg_dv('ATV', 'CND', sync_nds, (/stg_bfr_con_ar/))
      call this%cfg_dv('KTW', 'CND', sync_nds, (/stg_bfr_con_ar/))
      call this%cfg_dv('KTS', 'CND', sync_nds, (/stg_bfr_con_ar/))
    end if
    call this%cfg_dv('GWFHEAD', 'FMI', sync_nds, (/stg_bfr_exg_ad/))
    call this%cfg_dv('GWFSAT', 'FMI', sync_nds, (/stg_bfr_exg_ad/))
    call this%cfg_dv('GWFSPDIS', 'FMI', sync_nds, (/stg_bfr_exg_ad/))
    call this%cfg_dv('GWFFLOWJA', 'FMI', sync_con, (/stg_bfr_exg_ad/))
    call this%cfg_dv('GWFFLOWJA', 'FMI', sync_exg, (/stg_bfr_exg_ad/), &
                     exg_var_name='GWFSIMVALS')
    ! -- Fill porosity from est packages, needed for cnd
    if (this%gweModel%incnd > 0 .and. this%gweModel%inest > 0) then
      call this%cfg_dv('POROSITY', 'EST', sync_nds, (/stg_aft_con_ar/))
    end if
    !
  end subroutine cfg_dist_vars
 
  !> @brief Allocate array variables for this connection
  !<
  subroutine allocate_arrays(this)
    class(gwegweconnectiontype) :: this !< the connection
 
    call mem_allocate(this%exgflowjaGwe, this%ig_builder%nrOfBoundaryCells, &
                      'EXGFLOWJAGWT', this%memoryPath)
 
  end subroutine allocate_arrays
 
  !> @brief Set required extent of the interface grid from
  !< the configuration
  subroutine setgridextent(this)
    ! -- dummy
    class(gwegweconnectiontype) :: this !< the connection
    ! -- local
    logical(LGP) :: hasAdv, hasCnd
    !
    hasadv = this%gweModel%inadv > 0
    hascnd = this%gweModel%incnd > 0
    !
    if (hasadv) then
      if (this%iIfaceAdvScheme == 2) then
        this%exg_stencil_depth = 2
        if (this%gweModel%adv%iadvwt == 2) then
          this%int_stencil_depth = 2
        end if
      end if
    end if
    !
    if (hascnd) then
      if (this%iIfaceXt3d > 0) then
        this%exg_stencil_depth = 2
        if (this%gweModel%cnd%ixt3d > 0) then
          this%int_stencil_depth = 2
        end if
      end if
    end if
    !
  end subroutine setgridextent
 
  !> @brief allocate and read/set the connection's data structures
  !<
  subroutine gwegwecon_ar(this)
    class(gwegweconnectiontype) :: this !< the connection
 
    ! 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()
 
    ! ... and now the interface model
    call this%gweInterfaceModel%model_ar()
 
    ! -- Set a pointer in the interface model to the gwecommon data
    if (this%gweModel%inest > 0) then
      this%gweInterfaceModel%gwecommon%gwecpw = this%gweModel%gwecommon%gwecpw
      this%gweInterfaceModel%gwecommon%gwerhow = this%gweModel%gwecommon%gwerhow
    end if
 
    !-- Set the equation scaling factor in the interface model to that of
    !   underlying GWE model
    if (this%gweModel%incnd > 0) then
      this%gweInterfaceModel%ieqnsclfac = this%gweModel%cnd%eqnsclfac
    end if
 
    ! -- AR the movers and obs through the exchange
    if (this%owns_exchange) then
      !cdl implement this when MVT is ready
      !cdl if (this%gweExchange%inmvt > 0) then
      !cdl   call this%gweExchange%mvt%mvt_ar()
      !cdl end if
      if (this%gweExchange%inobs > 0) then
        call this%gweExchange%obs%obs_ar()
      end if
    end if
 
  end subroutine gwegwecon_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(gwegweconnectiontype) :: this !< this connection
 
    ! -- Base validation, the spatial/geometry part
    call this%SpatialModelConnectionType%validateConnection()
 
    ! -- We cannot validate this (yet) in parallel mode
    if (.not. this%gweExchange%v_model1%is_local) return
    if (.not. this%gweExchange%v_model2%is_local) return
 
    ! -- GWE related matters
    if ((this%gweExchange%gwemodel1%inadv > 0 .and. &
         this%gweExchange%gwemodel2%inadv == 0) .or. &
        (this%gweExchange%gwemodel2%inadv > 0 .and. &
         this%gweExchange%gwemodel1%inadv == 0)) then
      write (errmsg, '(1x,a,a,a)') 'Cannot connect GWE models in exchange ', &
        trim(this%gweExchange%name), ' because one model is configured with ADV &
        &and the other one is not'
      call store_error(errmsg)
    end if
    !
    if ((this%gweExchange%gwemodel1%incnd > 0 .and. &
         this%gweExchange%gwemodel2%incnd == 0) .or. &
        (this%gweExchange%gwemodel2%incnd > 0 .and. &
         this%gweExchange%gwemodel1%incnd == 0)) then
      write (errmsg, '(1x,a,a,a)') 'Cannot connect GWE models in exchange ', &
        trim(this%gweExchange%name), ' because one model is configured with CND &
        &and the other one is not'
      call store_error(errmsg)
    end if
    !
    ! Abort on errors
    if (count_errors() > 0) then
      write (errmsg, '(a)') 'Errors occurred while processing exchange(s)'
      call ustop()
    end if
    !
  end subroutine validateconnection
 
  subroutine gwegwecon_rp(this)
    ! -- dummy
    class(gwegweconnectiontype) :: this !< the connection
    !
    ! Call exchange rp routines
    if (this%owns_exchange) then
      call this%gweExchange%exg_rp()
    end if
    !
  end subroutine gwegwecon_rp
 
  !> @brief Advance this connection
  !<
  subroutine gwegwecon_ad(this)
    !
    class(gwegweconnectiontype) :: this !< this connection
    !
    ! -- Recalculate conduction ellipse
    if (this%gweInterfaceModel%incnd > 0) call this%gweInterfaceModel%cnd%cnd_ad()
    !
    if (this%owns_exchange) then
      call this%gweExchange%exg_ad()
    end if
    !
  end subroutine gwegwecon_ad
 
  subroutine gwegwecon_fc(this, kiter, matrix_sln, rhs_sln, inwtflag)
    ! -- dummy
    class(gwegweconnectiontype) :: this !< the connection
    integer(I4B), intent(in) :: kiter !< the iteration counter
    class(matrixbasetype), pointer :: matrix_sln !< the system matrix
    real(DP), dimension(:), intent(inout) :: rhs_sln !< global right-hand-side
    integer(I4B), optional, intent(in) :: inwtflag !< newton-raphson flag
    !
    call this%SpatialModelConnectionType%spatialcon_fc( &
      kiter, matrix_sln, rhs_sln, inwtflag)
    !
    ! FC the movers through the exchange
    if (this%owns_exchange) then
      if (this%gweExchange%inmvt > 0) then
        call this%gweExchange%mvt%mvt_fc(this%gweExchange%gwemodel1%x, &
                                         this%gweExchange%gwemodel2%x)
      end if
    end if
    !
  end subroutine gwegwecon_fc
 
  subroutine gwegwecon_cq(this, icnvg, isuppress_output, isolnid)
    ! -- dummy
    class(gwegweconnectiontype) :: 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%gweInterfaceModel%model_cq(icnvg, isuppress_output)
    call this%setFlowToExchange()
    !
  end subroutine gwegwecon_cq
 
  !> @brief Set the flows (flowja from interface model) to the
  !< simvals in the exchange, leaving the budget calcution in there
  subroutine setflowtoexchange(this)
    ! -- modules
    use indexmapmodule
    ! -- dummy
    class(gwegweconnectiontype) :: this !< this connection
    ! -- local
    integer(I4B) :: i
    class(gweexchangetype), pointer :: gweEx
    type(indexmapsgntype), pointer :: map
    !
    if (this%owns_exchange) then
      gweex => this%gweExchange
      map => this%interface_map%exchange_maps(this%interface_map%prim_exg_idx)
      !
      ! -- Use (half of) the exchange map in reverse:
      do i = 1, size(map%src_idx)
        if (map%sign(i) < 0) cycle ! simvals is defined from exg%m1 => exg%m2
        gweex%simvals(map%src_idx(i)) = &
          this%gweInterfaceModel%flowja(map%tgt_idx(i))
      end do
    end if
    !
  end subroutine setflowtoexchange
 
  subroutine gwegwecon_bd(this, icnvg, isuppress_output, isolnid)
    ! -- modules
    use budgetmodule, only: rate_accumulator
    ! -- dummy
    class(gwegweconnectiontype) :: 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%owns_exchange) then
      call this%gweExchange%exg_bd(icnvg, isuppress_output, isolnid)
    end if
    !
  end subroutine gwegwecon_bd
 
  subroutine gwegwecon_ot(this)
    ! -- dummy
    class(gwegweconnectiontype) :: this !< the connection
    !
    ! -- Call exg_ot() here as it handles all output processing
    !    based on gweExchange%simvals(:), which was correctly
    !    filled from gwegwecon
    if (this%owns_exchange) then
      call this%gweExchange%exg_ot()
    end if
    !
  end subroutine gwegwecon_ot
 
  subroutine gwegwecon_da(this)
    ! -- dummy
    class(gwegweconnectiontype) :: 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%exgflowjaGwe)
    !
    ! -- Interface model
    call this%gweInterfaceModel%model_da()
    deallocate (this%gweInterfaceModel)
    !
    ! -- Dealloc base
    call this%spatialcon_da()
    !
    inquire (this%iout, opened=isopen)
    if (isopen) then
      close (this%iout)
    end if
    !
    ! -- We need to deallocate the exchange we own:
    if (this%owns_exchange) then
      call this%gweExchange%exg_da()
    end if
    !
  end subroutine gwegwecon_da
 
  !> @brief Cast to GweGweConnectionType
  !<
  function castasgwegweconnection(obj) result(res)
    implicit none
    ! -- dummy
    class(*), pointer, intent(inout) :: obj !< object to be cast
    ! -- return
    class(gwegweconnectiontype), pointer :: res !< the GweGweConnection
    !
    res => null()
    if (.not. associated(obj)) return
    !
    select type (obj)
    class is (gwegweconnectiontype)
      res => obj
    end select
    !
  end function castasgwegweconnection
 
end module