integration.f90

Go to the documentation of this file.

module integration
    implicit none
    private
    public integrate
contains
!********************************BEGINNING*************************************
subroutine integrate
    use constants
    use globals
    use physicalproperties
    use conductivity, only: equcond
    use ioutils, only: newunit
    use model, only: modelpu,ngas,nfv,dz,dif,sol,bc,mor
    use inout, only: input,output,print_header
    integer :: i, j, k, l, counter, fi
    integer :: multiplicator
    integer :: itol, itask, istate, iopt
    integer :: MF, ML, MU, LRW, LIW, LENRAT, NNZ, LWM, NEQ
    integer, allocatable :: IWORK(:)

    real(dp) :: tin, tout, keq
    real(dp) :: rtol, atol, jdem

    real(dp), allocatable :: ystate(:), yprime(:) ! vector of state
    real(dp), allocatable :: RWORK(:)
    real(dp), allocatable :: yinit(:)
    real(dp), allocatable :: pp(:) ! partial pressure

    ! model should be general, but physical properties and conductivity are
    ! hardcoded for ngas=4
    ngas=4
    gasname(1)="O2"
    gasname(2)="N2"
    gasname(3)="CO2"
    gasname(4)="CyP"
    allocate(solmodel(ngas),diffmodel(ngas))
    allocate(pp(ngas),sg(ngas),dg(ngas),xg(ngas),sheetsg(ngas),sheetdg(ngas))
    allocate(pbg(ngas),kfoamxg(ngas),kgasxg(ngas))
    allocate(sgmodena(ngas),sginputs(ngas),sgoutputs(ngas))
    allocate(sgtemppos(ngas),sgxl1pos(ngas),sgxl2pos(ngas))
    allocate(dgmodena(ngas),dginputs(ngas),dgoutputs(ngas),dgtemppos(ngas))
    allocate(kgmodena(ngas),kginputs(ngas),kgoutputs(ngas),kgtemppos(ngas))
! -----------------------------------
! load inputs
! -----------------------------------
    call input
! -----------------------------------
! find out physical properties
! -----------------------------------
    dgas= gasdiffusivity(temp)
    call createmodels(ngas)
    eps=1-rhof/rhop
    do i=1,ngas
        if (solmodel(i)==1) then
            sg(i)=solubility(temp,i)
        endif
        if (diffmodel(i)==1) then
            dg(i)=diffusivity(temp,i)
        endif
    enddo
    call print_header
    sg=sg*rg*temp*1100._dp/(1e5*mg)
    sheetsg=sheetsg*rg*temp*1100._dp/(1e5*mg)
! -----------------------------------
! Allocate memory for working arrays
! -----------------------------------
    if (sheet) then
        ncell=nint((dfoam-dsheet)/(dcell+dwall))
    else
        ncell = nint(dfoam/(dcell+dwall))
        print*, ncell
        divsheet=0
    endif
    nfv = divsheet+ncell*(divwall+divcell)
    neq = ngas*nfv
    allocate(ystate(1:neq))
    allocate(yprime(1:neq))
    allocate(dz(nfv))
    allocate(mor(nfv))
    allocate(dif(neq))
    allocate(sol(neq))
    allocate(bc(ngas))
    allocate(yinit(ngas))
! ----------------------------------
! mesh and initial conditions
! ----------------------------------
    k=1
    do i=1,divsheet
        dz(k)=dsheet/divsheet
        mor(k)=3
        do l=1,ngas
            dif(ngas*(k-1)+l)=sheetdg(l)
            sol(ngas*(k-1)+l)=sheetsg(l)
            ystate(ngas*(k-1)+l)=xg(l)*pressure/rg/temp
        enddo
        k=k+1
    enddo
    do i = 1, ncell
        do j=1,divcell
            dz(k)=dcell/divcell
            mor(k)=1
            do l=1,ngas
                dif(ngas*(k-1)+l)=dgas
                sol(ngas*(k-1)+l)=1
                ystate(ngas*(k-1)+l)=xg(l)*pressure/rg/temp
            enddo
            k=k+1
        enddo
        do j=1,divwall
            dz(k)=dwall/divwall
            mor(k)=2
            do l=1,ngas
                dif(ngas*(k-1)+l)=dg(l)
                sol(ngas*(k-1)+l)=sg(l)
                ystate(ngas*(k-1)+l)=xg(l)*pressure/rg/temp
            enddo
            k=k+1
        enddo
    end do
! ----------------------------------
! boundary conditions
! ----------------------------------
    bc=pbg/rg/temp
! ----------------------------------
! initialize integration
! ----------------------------------
    lenrat = 2        ! usually for real(dp)
    nnz = neq**2      ! nonzero elements  - MV CHECK
    lwm = 2*nnz + 2*neq + (nnz+10*neq)/lenrat     ! MITER = 2
    liw = 31 + neq + nnz +100
    lrw = 20 + (2 + 1._dp/lenrat)*nnz + (11 + 9._dp/lenrat)*neq
    mf = 222    ! 222 stiff, internally generated jacobi matrix
                ! 10 'normal' - chage the allocation to smaller fields
                ! 210 structure obtained NEQ+1 calls, implicit Adams,
    allocate(rwork(1:lrw))
    allocate(iwork(1:liw))
    open(10,file='dcmp_progress.out', status='unknown')
    counter = 1
    itol = 1
    rtol = 1.0e-10_dp
    atol = 1.0e-6_dp
    itask = 1
    istate = 1
    iopt = 0
    multiplicator = 100
! ----------------------------------
! Integration loop
! ----------------------------------
    open (newunit(fi),file='keq_time.out')
    write(fi,'(10A23)') '#time', 'eq_conductivity', 'total_pressure', 'p_O2', &
        'p_N2', 'p_CO2', 'p_CP'
    call output(0, 0.0_dp, ystate, neq, pp)
    call equcond(keq,ystate,ngas,nfv,mor,eps,dcell,fstrut,temp_cond)
    write(fi,'(10es23.15)') tbeg/(3600*24),keq*1.0e3_dp,sum(pp),pp
    do i = 1, nroutputs*multiplicator       ! stabilizing multiplicator
        tin  = dble(i-1)*(tend-tbeg)/dble(nroutputs*multiplicator)+tbeg
        tout = dble(i  )*(tend-tbeg)/dble(nroutputs*multiplicator)+tbeg
100     continue    ! try to make another run for the initial step simulation
        call dlsodes(modelpu, neq, ystate, tin, tout, itol, rtol, atol, itask,&
            istate, iopt, rwork, lrw, iwork, liw, jdem, mf)
        ! evaluating the integration
        if (istate.lt.0) then
            write(*,*) 'Something is wrong, look for ISTATE =', istate
            if (istate.eq.-1) then  ! not enough steps to reach tout
                istate = 3
                iopt = 1   ! start to change something
                rwork(5:8)=0.0e0_dp
                iwork(5) = 0
                iwork(6) = counter*1000
                iwork(7) = 0
                counter = counter + 2
                write(*,*) 'MAXSTEP', iwork(6)
                write(10,*) 'MAXSTEP', iwork(6)
                goto 100
            endif
            ! stop
        elseif (counter>3) then
            counter=counter-1
            iwork(6) = counter*1000
            write(*,*) 'MAXSTEP', iwork(6)
            write(10,*) 'MAXSTEP', iwork(6)
        endif
        ! some output
        if (mod(i,multiplicator).eq.0) then
            write(*,'(2x,A,1x,f6.1,1x,A)') 'time:', tout/(3600*24),'days'
            write(10,'(2x,A,1x,es9.3,1x,A)') 'time:', tout/(3600*24),'days'
            call output(i/multiplicator, tout, ystate, neq, pp)
            call equcond(keq,ystate,ngas,nfv,mor,eps,dcell,fstrut,temp_cond)
            write(fi,'(10es23.15)') tout/(3600*24),keq*1e3,sum(pp),pp
        endif
    enddo
    close(10)
    close(fi)
    call destroymodels(ngas)
end subroutine integrate
!***********************************END****************************************
end module integration