NumericalModel.f90

Go to the documentation of this file.

module numericalmodelmodule
 
  use kindmodule, only: dp, i4b
  use constantsmodule, only: linelength, lenbudtxt, lenpakloc
  use basemodelmodule, only: basemodeltype
  use sparsemodule, only: sparsematrix
  use timearrayseriesmanagermodule, only: timearrayseriesmanagertype
  use listmodule, only: listtype
  use matrixbasemodule
  use vectorbasemodule
 
  implicit none
  private
  public :: numericalmodeltype, addnumericalmodeltolist, &
            getnumericalmodelfromlist
 
  type, extends(basemodeltype) :: numericalmodeltype
    integer(I4B), pointer :: neq => null() !< number of equations
    integer(I4B), pointer :: moffset => null() !< offset of this model in the solution
    integer(I4B), pointer :: icnvg => null() !< convergence flag
    integer(I4B), dimension(:), pointer, contiguous :: ia => null() !< csr row pointer
    integer(I4B), dimension(:), pointer, contiguous :: ja => null() !< csr columns
    real(dp), dimension(:), pointer, contiguous :: x => null() !< dependent variable (head, conc, etc)
    real(dp), dimension(:), pointer, contiguous :: rhs => null() !< right-hand side vector
    real(dp), dimension(:), pointer, contiguous :: cond => null() !< conductance matrix
    integer(I4B), dimension(:), pointer, contiguous :: idxglo => null() !< pointer to position in solution matrix
    real(dp), dimension(:), pointer, contiguous :: xold => null() !< dependent variable for previous timestep
 
  contains
    ! Overridden methods
    procedure :: model_df
    procedure :: model_ar
    procedure :: model_fp
    procedure :: model_da
    ! Specific methods
    procedure :: model_ac
    procedure :: model_mc
    procedure :: model_rp
    procedure :: model_ad
    procedure :: model_reset
    procedure :: model_solve
    procedure :: model_cf
    procedure :: model_fc
    procedure :: model_ptcchk
    procedure :: model_ptc
    procedure :: model_nr
    procedure :: model_cc
    procedure :: model_nur
    procedure :: model_cq
    procedure :: model_bd
    procedure :: model_bdcalc
    procedure :: model_bdsave
    procedure :: model_ot
    procedure :: model_bdentry
    ! Utility methods
    procedure :: allocate_scalars
    procedure :: allocate_arrays
    procedure :: set_moffset
    procedure :: set_idsoln
    procedure :: set_xptr
    procedure :: set_rhsptr
    procedure :: set_iboundptr
    procedure :: get_mrange
    procedure :: get_mcellid
    procedure :: get_mnodeu
    procedure :: get_iasym
    procedure :: get_idv_scale
  end type numericalmodeltype
 
contains
 
  subroutine model_df(this)
    class(numericalmodeltype) :: this
  end subroutine model_df
 
  subroutine model_ac(this, sparse)
    class(numericalmodeltype) :: this
    type(sparsematrix), intent(inout) :: sparse
  end subroutine model_ac
 
  subroutine model_mc(this, matrix_sln)
    class(numericalmodeltype) :: this
    class(matrixbasetype), pointer :: matrix_sln
  end subroutine model_mc
 
  subroutine model_ar(this)
    class(numericalmodeltype) :: this
  end subroutine model_ar
 
  subroutine model_rp(this)
    class(numericalmodeltype) :: this
  end subroutine model_rp
 
  subroutine model_ad(this)
    class(numericalmodeltype) :: this
  end subroutine model_ad
 
  subroutine model_reset(this)
    use bndmodule, only: bndtype, getbndfromlist
    class(numericalmodeltype) :: this
    ! local
    class(bndtype), pointer :: packobj
    integer(I4B) :: ip
 
    do ip = 1, this%bndlist%Count()
      packobj => getbndfromlist(this%bndlist, ip)
      call packobj%bnd_reset()
    end do
 
  end subroutine model_reset
 
  subroutine model_solve(this)
    class(numericalmodeltype) :: this
  end subroutine model_solve
 
  subroutine model_cf(this, kiter)
    class(numericalmodeltype) :: this
    integer(I4B), intent(in) :: kiter
  end subroutine model_cf
 
  subroutine model_fc(this, kiter, matrix_sln, inwtflag)
    class(numericalmodeltype) :: this
    integer(I4B), intent(in) :: kiter
    class(matrixbasetype), pointer :: matrix_sln
    integer(I4B), intent(in) :: inwtflag
  end subroutine model_fc
 
  subroutine model_ptcchk(this, iptc)
    class(numericalmodeltype) :: this
    integer(I4B), intent(inout) :: iptc
    iptc = 0
  end subroutine model_ptcchk
 
  subroutine model_ptc(this, vec_residual, iptc, ptcf)
    class(numericalmodeltype) :: this
    class(vectorbasetype), pointer :: vec_residual
    integer(I4B), intent(inout) :: iptc
    real(DP), intent(inout) :: ptcf
  end subroutine model_ptc
 
  subroutine model_nr(this, kiter, matrix, inwtflag)
    class(numericalmodeltype) :: this
    integer(I4B), intent(in) :: kiter
    class(matrixbasetype), pointer :: matrix
    integer(I4B), intent(in) :: inwtflag
  end subroutine model_nr
 
  subroutine model_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
    class(numericalmodeltype) :: this
    integer(I4B), intent(in) :: innertot
    integer(I4B), intent(in) :: kiter
    integer(I4B), intent(in) :: iend
    integer(I4B), intent(in) :: icnvgmod
    character(len=LENPAKLOC), intent(inout) :: cpak
    integer(I4B), intent(inout) :: ipak
    real(DP), intent(inout) :: dpak
  end subroutine model_cc
 
  subroutine model_nur(this, neqmod, x, xtemp, dx, inewtonur, dxmax, locmax)
    class(numericalmodeltype) :: this
    integer(I4B), intent(in) :: neqmod
    real(DP), dimension(neqmod), intent(inout) :: x
    real(DP), dimension(neqmod), intent(in) :: xtemp
    real(DP), dimension(neqmod), intent(inout) :: dx
    integer(I4B), intent(inout) :: inewtonur
    real(DP), intent(inout) :: dxmax
    integer(I4B), intent(inout) :: locmax
  end subroutine model_nur
 
  subroutine model_cq(this, icnvg, isuppress_output)
    class(numericalmodeltype) :: this
    integer(I4B), intent(in) :: icnvg
    integer(I4B), intent(in) :: isuppress_output
  end subroutine model_cq
 
  subroutine model_bd(this, icnvg, isuppress_output)
    class(numericalmodeltype) :: this
    integer(I4B), intent(in) :: icnvg
    integer(I4B), intent(in) :: isuppress_output
  end subroutine model_bd
 
  subroutine model_bdcalc(this, icnvg)
    class(numericalmodeltype) :: this
    integer(I4B), intent(in) :: icnvg
  end subroutine model_bdcalc
 
  subroutine model_bdsave(this, icnvg)
    class(numericalmodeltype) :: this
    integer(I4B), intent(in) :: icnvg
  end subroutine model_bdsave
 
  subroutine model_ot(this)
    class(numericalmodeltype) :: this
  end subroutine model_ot
 
  subroutine model_bdentry(this, budterm, budtxt, rowlabel)
    class(numericalmodeltype) :: this
    real(DP), dimension(:, :), intent(in) :: budterm
    character(len=LENBUDTXT), dimension(:), intent(in) :: budtxt
    character(len=*), intent(in) :: rowlabel
  end subroutine model_bdentry
 
  subroutine model_fp(this)
    class(numericalmodeltype) :: this
  end subroutine model_fp
 
  subroutine model_da(this)
    ! modules
    use memorymanagermodule, only: mem_deallocate
    class(numericalmodeltype) :: this
 
    ! deallocate scalars
    call mem_deallocate(this%neq)
    call mem_deallocate(this%nja)
    call mem_deallocate(this%icnvg)
    call mem_deallocate(this%moffset)
    deallocate (this%filename)
 
    ! deallocate arrays
    call mem_deallocate(this%xold)
    call mem_deallocate(this%flowja, 'FLOWJA', this%memoryPath)
    call mem_deallocate(this%idxglo)
 
    ! deallocate derived types
    call this%bndlist%Clear()
    deallocate (this%bndlist)
 
    ! nullify pointers
    call mem_deallocate(this%x, 'X', this%memoryPath)
    call mem_deallocate(this%rhs, 'RHS', this%memoryPath)
    call mem_deallocate(this%ibound, 'IBOUND', this%memoryPath)
 
    ! base deallocation
    call this%BaseModelType%model_da()
 
  end subroutine model_da
 
  subroutine set_moffset(this, moffset)
    class(numericalmodeltype) :: this
    integer(I4B), intent(in) :: moffset
    this%moffset = moffset
  end subroutine set_moffset
 
  subroutine get_mrange(this, mstart, mend)
    class(numericalmodeltype) :: this
    integer(I4B), intent(inout) :: mstart
    integer(I4B), intent(inout) :: mend
    mstart = this%moffset + 1
    mend = mstart + this%neq - 1
  end subroutine get_mrange
 
  subroutine set_idsoln(this, id)
    class(numericalmodeltype) :: this
    integer(I4B), intent(in) :: id
    this%idsoln = id
  end subroutine set_idsoln
 
  subroutine allocate_scalars(this, modelname)
    use memorymanagermodule, only: mem_allocate
    class(numericalmodeltype) :: this
    character(len=*), intent(in) :: modelname
 
    ! allocate base members
    call this%BaseModelType%allocate_scalars(modelname)
 
    ! allocate members from this type
    call mem_allocate(this%neq, 'NEQ', this%memoryPath)
    call mem_allocate(this%nja, 'NJA', this%memoryPath)
    call mem_allocate(this%icnvg, 'ICNVG', this%memoryPath)
    call mem_allocate(this%moffset, 'MOFFSET', this%memoryPath)
    allocate (this%filename)
    allocate (this%bndlist)
 
    this%filename = ''
    this%neq = 0
    this%nja = 0
    this%icnvg = 0
    this%moffset = 0
  end subroutine allocate_scalars
 
  subroutine allocate_arrays(this)
    use constantsmodule, only: dzero
    use memorymanagermodule, only: mem_allocate
    class(numericalmodeltype) :: this
    integer(I4B) :: i
 
    call mem_allocate(this%xold, this%neq, 'XOLD', this%memoryPath)
    call mem_allocate(this%flowja, this%nja, 'FLOWJA', this%memoryPath)
    call mem_allocate(this%idxglo, this%nja, 'IDXGLO', this%memoryPath)
 
    do i = 1, size(this%flowja)
      this%flowja(i) = dzero
    end do
  end subroutine allocate_arrays
 
  subroutine set_xptr(this, xsln, sln_offset, varNameTgt, memPathTgt)
    use memorymanagermodule, only: mem_checkin
    ! dummy
    class(numericalmodeltype) :: this
    real(DP), dimension(:), pointer, contiguous, intent(in) :: xsln
    integer(I4B) :: sln_offset
    character(len=*), intent(in) :: varNameTgt
    character(len=*), intent(in) :: memPathTgt
    ! local
    integer(I4B) :: offset
 
    offset = this%moffset - sln_offset
    this%x => xsln(offset + 1:offset + this%neq)
    call mem_checkin(this%x, 'X', this%memoryPath, varnametgt, mempathtgt)
  end subroutine set_xptr
 
  subroutine set_rhsptr(this, rhssln, sln_offset, varNameTgt, memPathTgt)
    use memorymanagermodule, only: mem_checkin
    ! dummy
    class(numericalmodeltype) :: this
    real(DP), dimension(:), pointer, contiguous, intent(in) :: rhssln
    integer(I4B) :: sln_offset
    character(len=*), intent(in) :: varNameTgt
    character(len=*), intent(in) :: memPathTgt
    ! local
    integer(I4B) :: offset
 
    offset = this%moffset - sln_offset
    this%rhs => rhssln(offset + 1:offset + this%neq)
    call mem_checkin(this%rhs, 'RHS', this%memoryPath, varnametgt, mempathtgt)
  end subroutine set_rhsptr
 
  subroutine set_iboundptr(this, iboundsln, sln_offset, varNameTgt, memPathTgt)
    use memorymanagermodule, only: mem_checkin
    ! dummy
    class(numericalmodeltype) :: this
    integer(I4B), dimension(:), pointer, contiguous, intent(in) :: iboundsln
    integer(I4B) :: sln_offset
    character(len=*), intent(in) :: varNameTgt
    character(len=*), intent(in) :: memPathTgt
    ! local
    integer(I4B) :: offset
 
    offset = this%moffset - sln_offset
    this%ibound => iboundsln(offset + 1:offset + this%neq)
    call mem_checkin(this%ibound, 'IBOUND', this%memoryPath, varnametgt, &
                     mempathtgt)
  end subroutine set_iboundptr
 
  subroutine get_mcellid(this, node, mcellid)
    use bndmodule, only: bndtype, getbndfromlist
    ! dummy
    class(numericalmodeltype) :: this
    integer(I4B), intent(in) :: node
    character(len=*), intent(inout) :: mcellid
    ! local
    character(len=20) :: cellid
    integer(I4B) :: ip, ipaknode, istart, istop
    class(bndtype), pointer :: packobj
 
    if (node < 1) then
      cellid = ''
    else if (node <= this%dis%nodes) then
      call this%dis%noder_to_string(node, cellid)
    else
      cellid = '***ERROR***'
      ipaknode = node - this%dis%nodes
      istart = 1
      do ip = 1, this%bndlist%Count()
        packobj => getbndfromlist(this%bndlist, ip)
        if (packobj%npakeq == 0) cycle
        istop = istart + packobj%npakeq - 1
        if (istart <= ipaknode .and. ipaknode <= istop) then
          write (cellid, '(a, a, a, i0, a, i0, a)') '(', &
            trim(packobj%filtyp), '_', &
            packobj%ibcnum, '-', ipaknode - packobj%ioffset, ')'
          exit
        end if
        istart = istop + 1
      end do
    end if
    write (mcellid, '(i0, a, a, a, a)') this%id, '_', this%macronym, '-', &
      trim(adjustl(cellid))
  end subroutine get_mcellid
 
  subroutine get_mnodeu(this, node, nodeu)
    use bndmodule, only: bndtype, getbndfromlist
    class(numericalmodeltype) :: this
    integer(I4B), intent(in) :: node
    integer(I4B), intent(inout) :: nodeu
 
    if (node <= this%dis%nodes) then
      nodeu = this%dis%get_nodeuser(node)
    else
      nodeu = -(node - this%dis%nodes)
    end if
  end subroutine get_mnodeu
 
  function get_iasym(this) result(iasym)
    class(numericalmodeltype) :: this
    integer(I4B) :: iasym
    iasym = 0
  end function get_iasym
 
  function get_idv_scale(this) result(idv_scale)
    class(numericalmodeltype) :: this
    integer(I4B) :: idv_scale
    idv_scale = 0
  end function get_idv_scale
 
  function castasnumericalmodelclass(obj) result(res)
    implicit none
    class(*), pointer, intent(inout) :: obj
    class(numericalmodeltype), pointer :: res
 
    res => null()
    if (.not. associated(obj)) return
 
    select type (obj)
    class is (numericalmodeltype)
      res => obj
    end select
  end function castasnumericalmodelclass
 
  subroutine addnumericalmodeltolist(list, model)
    implicit none
    ! dummy
    type(listtype), intent(inout) :: list
    class(numericalmodeltype), pointer, intent(inout) :: model
    ! local
    class(*), pointer :: obj
 
    obj => model
    call list%Add(obj)
  end subroutine addnumericalmodeltolist
 
  function getnumericalmodelfromlist(list, idx) result(res)
    implicit none
    ! dummy
    type(listtype), intent(inout) :: list
    integer(I4B), intent(in) :: idx
    class(numericalmodeltype), pointer :: res
    ! local
    class(*), pointer :: obj
 
    obj => list%GetItem(idx)
    res => castasnumericalmodelclass(obj)
  end function getnumericalmodelfromlist
 
end module numericalmodelmodule