solution2.f90

Go to the documentation of this file.

! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +    

SUBROUTINE solution2 (SOLVER, ifail)

!-------------------------------------------------------!
!                                                       !
!     This subroutine is prepared to solve single       !    
!     transport equation in standardised form adopted   !
!     by the ETS.                                       !
!                                                       !
!-------------------------------------------------------!  
!                                                       !
!     Source:       provided by M.Tokar on 17.02.2009   !
!     Developers:   M.Tokar, D.Kalupin                  !
!     Kontacts:     M.Tokar@fz-juelich.de               !
!                   D.Kalupin@fz-juelich.de             !
!                                                       !
!     Comments:     equation is solved in               !
!                   INTEGRAL form                       !
!                                                       !
!-------------------------------------------------------!  




      USE type_transport_solver_numerics

      USE ids_schemas

      IMPLICIT NONE

      INTEGER, INTENT (INOUT)         :: ifail

! +++ Input/Output with ETS:
      TYPE (solver_io), INTENT (INOUT) :: solver               !contains all I/O quantities to numerics part


! +++ Internal input/output parameters:
      INTEGER   :: idim,    ndim                              !equation index and total number of equations to be solved

      INTEGER   :: irho,    nrho                              !radius index, number of radial points

      INTEGER   :: flag                                       !flag for equation: 0 - interpretative (not solved), 1 - predictive (solved)

      REAL (IDS_REAL) :: rho(solver%nrho)                           !radii

      REAL (IDS_REAL) :: amix                                       !fraction of new sollution mixed

      REAL (IDS_REAL) :: y(solver%nrho), ym(solver%nrho)            !function at the current amd previous time steps
      REAL (IDS_REAL) :: dy(solver%nrho)                            !derivative of function

      REAL (IDS_REAL) :: a(solver%nrho), b(solver%nrho)             !coefficients
      REAL (IDS_REAL) :: c(solver%nrho), d(solver%nrho)             !coefficients 
      REAL (IDS_REAL) :: e(solver%nrho), f(solver%nrho)             !coefficients 
      REAL (IDS_REAL) :: g(solver%nrho), h                          !coefficients 
      REAL (IDS_REAL) :: dd(solver%nrho), de(solver%nrho)           !derivatives of coefficients 

      REAL (IDS_REAL) :: v(2), u(2), w(2)                           !boundary conditions 






! + + + + + + + + + INTERFACE PART  + + + + + + + + + + +    

! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! +++ Set up local variables (input, obtained from ETS):
!     Control parameters:
  ndim           = solver%NDIM
  nrho           = solver%NRHO
  amix           = solver%AMIX



! +++ Solution of equations starting from 1 to NDIM:
  equation_loop1: DO idim = 1, ndim

     flag         = solver%EQ_FLAG(idim)

     IF (flag.EQ.0) goto 20                                !equation is not solved



! +++ Numerical coefficients obtained from physics part  in form:
!
!     (A*Y-B*Y(t-1))/H + 1/C * (-D*Y' + E*Y) = F - G*Y

     h            = solver%H

     rho_loop1: DO irho=1,nrho
        rho(irho)  = solver%RHO(irho)

        y(irho)    = solver%Y(idim,irho)
        dy(irho)   = solver%DY(idim,irho)
        ym(irho)   = solver%YM(idim,irho)

        a(irho)    = solver%A(idim,irho) 
        b(irho)    = solver%B(idim,irho)
        c(irho)    = solver%C(idim,irho)
        d(irho)    = solver%D(idim,irho) 
        e(irho)    = solver%E(idim,irho)
        f(irho)    = solver%F(idim,irho)
        g(irho)    = solver%G(idim,irho) 
     END DO rho_loop1


! +++ Boundary conditions for numerical solver in form:
!
!     V*Y' + U*Y = W 

! +++ On axis:

     v(1)         = solver%V(idim,1)
     u(1)         = solver%U(idim,1)
     w(1)         = solver%W(idim,1) 

! +++ At the edge:

     v(2)         = solver%V(idim,2) 
     u(2)         = solver%U(idim,2) 
     w(2)         = solver%W(idim,2) 



! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! +++  Solution of transport equation:      

     CALL solver_ritm_integral        &
          (nrho,rho,amix,  y,ym,dy,   a,b,c,d,e,f,g,h,  u,v,w)
!              radii   b.c.       coeff.      solution   spline


        rho_loop3: DO irho = 1,nrho
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! +++ Return solution to ETS:      
          solver%Y(idim,irho)  = y(irho)
          solver%DY(idim,irho) = dy(irho)

        END DO rho_loop3



20   CONTINUE



  END DO equation_loop1



  RETURN



END SUBROUTINE solution2


! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  





! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  

! + + + + + + + + + NUMERICAL PART  + + + + + + + + + + +    
      SUBROUTINE solver_ritm_integral (n,r,amix,z,zm,dz,       &
          aets,bets,cets,dets,eets,fets,gets,tau,uets,vets,wets)

      USE ids_schemas

      IMPLICIT NONE

      INTEGER         i,n

      REAL (IDS_REAL)         tau,xi,p,amix,error,errornew
!DPC 2009-02-28 Changed h(n) to h(n+1) to match definition in solut3
      REAL (IDS_REAL)         h(n+1),r(n),z(n),zm(n)
      REAL (IDS_REAL)         dz(n),dzm(n)
      REAL (IDS_REAL)         aets(n),bets(n),cets(n)
      REAL (IDS_REAL)         dets(n),eets(n),gets(n),fets(n),gamets(n)
      REAL (IDS_REAL)         uets(2),vets(2),wets(2)
      REAL (IDS_REAL)         cs(n),csn(n)
      REAL (IDS_REAL)         al(n),be(n),so(n),nu(n),dif(n),vel(n),dfl(n)
      REAL (IDS_REAL)         om(n),dom(n),s(n),ints(n)
      REAL (IDS_REAL)         a(n),b(n),f(n),u(2),v(2),w(2)
      REAL (IDS_REAL)         y(n),d(n)
      REAL (IDS_REAL)         dd(n)

! +++ ETS - standard equation in cylindrical geometry
      DO i=2,n
         al(i)     = aets(i)*cets(i)/r(i)
         be(i)     = bets(i)*cets(i)/r(i)
         so(i)     = fets(i)*cets(i)/r(i)
         nu(i)     = gets(i)*cets(i)/r(i)
         dif(i)    = dets(i)/r(i)
         vel(i)    = eets(i)/r(i)
         gamets(i) = 0.d0
         dfl(i)    = gamets(i)/r(i)
         h(i-1)    = r(i)-r(i-1) 
      END DO
      h(n)=0.0_ids_real                        !DPC 2009-02-28

      FORALL(i=1:n) cs(i)=0.d0

      CALL derivn2(n,r,zm,dzm)

! +++ parameters for integral representation
 10   xi=1.d-1

      DO i=2,n
         p     = (1.d0-xi)*al(i)+tau*nu(i)
         om(i) = xi*al(i)+(p+dabs(p))/2.d0
         s(i)  = so(i)+((p-dabs(p))/2.d0*z(i)-(om(i)-be(i))*zm(i))/tau
      END DO

      CALL axis(n,h,r,om)

      CALL axis(n,h,r,s)


! +++ coefficients for differential equation
      CALL derivn2(n,r,om,dom)
      CALL integral(n,h,r,s,ints)

      DO i=2,n
         a(i) = 3.d0/r(i)-dom(i)/om(i)-vel(i)/dif(i)
         b(i) = om(i)/tau/dif(i)+2.d0/r(i)*(3.d0/r(i)-a(i))
         f(i) = om(i)/r(i)/dif(i)*                             &
                (ints(i)/r(i)+dif(i)*dzm(i)-vel(i)*zm(i)-dfl(i))
      END DO


! +++ boundary conditions
      IF(wets(1).NE.0.d0) THEN
         write(6,*) 'wrong boundary condition: w_1 is not 0' 
         stop
      END IF
    
      IF(uets(1).EQ.0.d0) THEN                            !dz/dr(0)=0
         u(1)=b(2)
         v(1)=-1.d0/r(2)-a(2)
         w(1)=f(2)
      END IF
    
      IF(vets(1).EQ.0.d0) THEN                            !z(0)=0
         u(1) = b(2)+2.d0/r(2)**2
         v(1) = -a(2)-2.d0/r(2)
         w(1) = f(2)
      END IF

      w(2) = vets(2)*f(n)+                          &
             om(n)/r(n)*(wets(2)-uets(2)*zm(n)-vets(2)*dzm(n))
      v(2) = uets(2)+vets(2)*vel(n)/dif(n)
      u(2) = vets(2)*om(n)/tau/dif(n)+2.d0/r(n)*v(2)


! +++ solution: cs=z-zm 
      CALL solut3(n,h,a,b,f,u,v,w,y,d)

      CALL axis(n,h,r,y)

      IF(vets(1).EQ.0.d0) y(1)=0.d0

!      CALL DERIVN2 (n,r,y,d)

      DO i=1,n
         csn(i)   = (r(i)*d(i)+2.d0*y(i))/om(i)
         cs(i)    = cs(i)*(1.d0-amix)+csn(i)*amix
         z(i)     = zm(i)+cs(i)
      END DO

      CALL derivn2(n,r,z,dz)

! +++ error in solution
      error=0.d0

      DO i=1,n
         errornew  =dsqrt(dabs(1.d0-csn(i)/cs(i)))
         IF(cs(i).NE.0.d0.AND.error.LT.errornew) error=errornew
      END DO

      IF (error.GT.1.d-5*amix) goto 10


      RETURN


      END SUBROUTINE solver_ritm_integral
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  





! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
      SUBROUTINE solut3 (n,h,a,b,f,u,v,w,y,dy)
!-------------------------------------------------------!
!                                                       !
!     This subroutine provides solution for the equation!
!                                                       !
!     d(dy/dx)/dx + a(x)*dy/dx = b(x)*y - f(x)          !
!                                                       !
!-------------------------------------------------------!
      USE ids_schemas

      IMPLICIT NONE

      INTEGER   i,j,k,l,m,n

      REAL (IDS_REAL)   h(n+1),a(n),b(n),f(n),y(n),dy(n),lam(10001,2)
      REAL (IDS_REAL)   u(2),v(2),w(2)
      REAL (IDS_REAL)   y0(n),y1(10001,2,2)
      REAL (IDS_REAL)   et(2,2),mu(2),det
      REAL (IDS_REAL)   c(n),d(n),g(2),x
      REAL (IDS_REAL)   cdy(n,2)

      DO i=2,n
         y0(i)        = f(i)/b(i)
         det          = a(i)**2/4.d0+b(i)

         DO k=1,2
            lam(i,k)  = -a(i)/2.d0-(-1)**k*dsqrt(det)


            FORALL(l=1:2)                                   &
            y1(i,k,l) = dexp(dmin1(200.d0,(-1)**l*lam(i,k)  &
                        *h(i-2+l)/2.d0))
         END DO

      END DO

     
      c(3)            = v(1)/(v(1)-h(3)*u(1))
      d(3)            = h(3)*w(1)/(h(3)*u(1)-v(1))


      DO i=3,n-1

         DO m=1,2
            l         = 3-m
            j         = i+(-1)**m

            DO k=1,2 
              et(k,m) = -(-1)**m*(lam(j,m)-lam(i,k)        &
                        -(lam(j,l)-lam(i,k))*y1(j,l,l)     &
                        /y1(j,m,l))/(lam(j,1)-lam(j,2))
            END DO

            mu(m)     = -(-1)**m*(y0(i)-y0(j))/y1(i,l,m)   &
                        *(lam(j,m)-lam(j,l)*y1(j,l,l)      &
                        /y1(j,m,l))/(lam(j,1)-lam(j,2))
         END DO

         det          = et(1,2)*et(2,1)-et(1,1)*et(2,2)    &
                        *y1(i,l,1)/y1(i,2,1)               &
                        *y1(i,2,2)/y1(i,l,2)

         dy(i)        = 0.0d0

         DO m=1,2
            l         = 3-m
            g(m)      = y1(i+(-1)**m,l,l)/y1(i,l,m)*       &
                        (et(m,l)-y1(i,l,l)/y1(i,m,l)       &
                        *et(l,l))/det
!        derivative:
            dy(i)     = dy(i)+lam(i,m)*                    &
                        (mu(m)*y1(i,l,l)/y1(i,m,l)*        &
                         et(l,l)-mu(l)*et(l,m))/det
            cdy(i,m)  = (et(m,l)*lam(i,l)-                 &
                         y1(i,l,l)/y1(i,m,l)*et(l,l)*lam(i,m))/det
         END DO

         x            = (mu(1)*(y1(i,2,2)/y1(i,1,2)        &
                        *et(2,2)-et(1,2))+mu(2)*(y1(i,1,1) &
                        /y1(i,2,1)*et(1,1)-et(2,1)))/det-  &
                        y0(i-1)*g(1)-y0(i+1)*g(2)+y0(i)
         c(i+1)       = g(2)/(1.d0-g(1)*c(i))
         d(i+1)       = (g(1)*d(i)+x)/(1.d0-g(1)*c(i))

      END DO

      y(n)            = (h(n-1)*w(2)+v(2)*d(n))            &
                        /(h(n-1)*u(2)+v(2)-v(2)*c(n))

      DO i=n,3,-1
         y(i-1)       = c(i)*y(i)+d(i)
      END DO


!     derivative:
      DO i=3,n-1
         DO m=1,2
            l         = 3-m
            j         = i+(-1)**m
            dy(i)     = dy(i)+(y(j)-y0(j))*y1(j,l,l)/y1(i,l,m)*cdy(i,m)
         END DO
      END DO

      dy(2)=(y(3)-y(3))/h(2)
      dy(n)=(y(n)-y(n-1))/h(n-1)


      RETURN


      END SUBROUTINE solut3
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  






! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
     SUBROUTINE derivn2(N,X,Y,DY1)
!-------------------------------------------------------!
!                                                       !
!    This subroutine calculates first and               !
!    second derivatives, DY1 and DY2,                   !
!    of function Y respect to argument X                !
!                                                       !
!-------------------------------------------------------!
      USE ids_schemas

      IMPLICIT NONE

      INTEGER      n, i
      REAL (IDS_REAL)    x(n), y(n), dy1(n)
      REAL (IDS_REAL)    h(n), dy2(n)

      SAVE




      DO i=2,n-1
        dy1(i)  = (y(i+1)-y(i-1))/(x(i+1)-x(i-1)) 
      END DO


      dy1(1)    = (y(2)-y(1))/(x(2)-x(1))
      dy1(n)    = (y(n)-y(n-1))/(x(n)-x(n-1))


      RETURN


      END SUBROUTINE derivn2
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  






! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
     SUBROUTINE derivn21(N,X,Y,DY1)
!-------------------------------------------------------!
!                                                       !
!    This subroutine calculates first and               !
!    second derivatives, DY1 and DY2,                   !
!    of function Y respect to argument X                !
!                                                       !
!-------------------------------------------------------!
      USE ids_schemas

      IMPLICIT NONE

      INTEGER      n, i
      REAL (IDS_REAL)    x(n), y(n), dy1(n)
      REAL (IDS_REAL)    h(n), dy2(n)

      SAVE


      DO i=1,n-1
        h(i)    = x(i+1)-x(i)
      END DO


      DO i=2,n-1
        dy1(i)  = ((y(i+1)-y(i))*h(i-1)/h(i)                 &
                  +(y(i)-y(i-1))*h(i)/h(i-1))/(h(i)+h(i-1))
        dy2(i)  = 2.d0*((y(i-1)-y(i))/h(i-1)                 &
                  +(y(i+1)-y(i))/h(i))/(h(i)+h(i-1))
      END DO


      dy1(1)    = dy1(2)-dy2(2)*h(1)
      dy1(n)    = dy1(n-1)+dy2(n-1)*h(n-1)


      RETURN


      END SUBROUTINE derivn21
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  






! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
      SUBROUTINE integral (n, h, r, f, int)
!-------------------------------------------------------!
!                                                       !
!     This subroutine calculates integral               !
!     from 0 to r_n of f*r byinterpolating f            !
!     with a local quadratic spline                     !
!                                                       !
!-------------------------------------------------------!

      USE ids_schemas

      IMPLICIT NONE

      INTEGER       i, n

      REAL (IDS_REAL)       h(n), r(n), f(n), int(n)
      REAL (IDS_REAL)       den
      REAL (IDS_REAL)       a(n), b(n)


      DO i=2,n-1
         den     = h(i-1)*(h(i-1)+h(i))*h(i)
         a(i)    = ((f(i+1)-f(i))*h(i-1)**2-(f(i-1)-f(i))*h(i)**2)/den
         b(i)    = ((f(i+1)-f(i))*h(i-1)+(f(i-1)-f(i))*h(i))/den
      END DO


      int(1)     = 0.d0

      int(2)     = f(2)*r(2)**2/2.d0-a(2)*r(2)**3/6.d0+b(2)*r(2)**4/12.d0


      DO i=2,n-2 
        int(i+1) = int(i)+(r(i)*f(i)+r(i+1)*f(i+1))*h(i)/2.d0+           &
                   (f(i)+r(i)*a(i)-f(i+1)-r(i+1)*a(i+1))*h(i)**2/8.d0+   &
                   (a(i)+r(i)*b(i)+a(i+1)+r(i+1)*b(i+1))*h(i)**3/24.d0+  &
                   (b(i)-b(i+1))*h(i)**4/64.d0
      END DO


      int(n)     = int(n-1)+(f(n-1)-a(n-1)*r(n-1)+b(n-1)*r(n-1)**2)*     &
                   (r(n)**2-r(n-1)**2)/2.d0+(a(n-1)-2.d0*b(n-1)*r(n-1))* &
                   (r(n)**3-r(n-1)**3)/3.d0+b(n-1)*(r(n)**4-r(n-1)**4)/4.d0


      RETURN


      END SUBROUTINE integral
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  






! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
      SUBROUTINE axis(n, h, r, f)
!-------------------------------------------------------!
!                                                       !
!     This subroutine finds                             !
!     f(r_1=0) from f(r_2), f(r_3) and f(r_4)           !
!                                                       !
!-------------------------------------------------------!  

      USE ids_schemas

      IMPLICIT NONE

      INTEGER    n
      REAL (IDS_REAL)    h(n), r(n), f(n)

      f(1)     = ((f(2)*r(4)/h(2)+f(4)*r(2)/h(3))*r(3)        &
                 -f(3)*(r(2)/h(2)+r(4)/h(3))*r(2)*r(4)/r(3))  &
                 /(r(4)-r(2))


      RETURN


      END SUBROUTINE axis
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +