solution_interface.f90

Go to the documentation of this file.

! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! + + + + + + INTERFACE TO NUMERICAL SOLVER + + + + + + +    

      SUBROUTINE solution_interface (SOLVER, ifail)

! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
!     This interface is prepared to test the ETS        +
!     and will be replaced by KEPLER                    +
!                                                       +
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
!     Source:       --                                  +
!     Developers:   D.Kalupin                           +
!     Kontacts:     D.Kalupin@fz-juelich.de             +
!                                                       +
!     Comments:     should be replaced by KEPLER        +
!                                                       +
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  


      USE ets_plasma
      USE type_solver
      USE itm_types
      USE ets_math

      IMPLICIT NONE


      TYPE (numerics), INTENT (INOUT) :: solver          
      INTEGER,         INTENT (INOUT) :: ifail
      INTEGER                         :: i
      REAL (R8)                       :: fun(solver%nrho) 
      REAL (R8)                       :: dfun(solver%nrho)
      REAL (R8)                       :: intfun(solver%nrho)
      LOGICAL, parameter              :: debug = .true.



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

      choose_numerical_solver: IF (solver%TYPE.EQ.1) THEN 
!       call numerical solver extracted from RITM, with differential form of equation
         CALL solution1(solver, ifail)

      ELSE IF (solver%TYPE.EQ.2) THEN
!       call numerical solver extracted from RITM, with integral form of equation
         CALL solution2(solver, ifail)

      ELSE IF (solver%TYPE.EQ.3) THEN                        
!       call matrix numerical solver extracted from COREDIF 
         CALL solution3(solver, ifail) 

      ELSE IF (solver%TYPE.EQ.4) THEN
!       call the improved solver 3 (progonka/block thomas) from Roman Stankiewicz
         CALL solution4(solver, ifail)

      ELSE IF (solver%TYPE.EQ.6) THEN
!       call solver from Guido Huysmans
         CALL solution6(solver, ifail)

      ELSE IF (solver%TYPE.EQ.7 .OR. solver%TYPE.EQ.8) THEN
!       call numerical solver extracted from RITM, with differential form of equation
         CALL solution7(solver, ifail)

      ELSE IF (solver%TYPE.EQ.9) THEN
!       call solver from Guido Huysmans
!        CALL SOLUTION9 (SOLVER, ifail)
         
      ELSE IF (solver%TYPE.EQ.10) THEN
!       call solver from CRONOS
         CALL solution10(solver, ifail) 

      ELSE IF (solver%TYPE.EQ.11) THEN
!       call solver from FESB (Anna Šušnjara)
         CALL solutionfem(solver, ifail)

      ELSE 
         WRITE(*,*) 'Invalid SOLVER%TYPE = ',solver%TYPE
      END IF choose_numerical_solver

      IF (ifail.GT.2) ifail = 0




      check_derivative: DO i = 1,solver%NDIM
         derivative_recalculate: IF (solver%DERIVATIVE_FLAG(i).EQ.1) THEN
            fun(:)                      = solver%Y(i,:)
            CALL   deriv_fun(solver%NRHO, solver%RHO, fun, dfun)
            solver%DY(i,:)              = dfun(:)

         ELSE IF (solver%DERIVATIVE_FLAG(i).EQ.2) THEN
            fun(:)                      = solver%C(i,:) &
                                         *((solver%F(i,:) + solver%B(i,:)*solver%YM(i,:)/solver%H) &
                                         -(solver%G(i,:) + solver%A(i,:)/solver%H)*solver%Y(i,:))

            CALL integr_fun(solver%NRHO, solver%RHO, fun, intfun)

            solver%DY(i,:)              = (solver%E(i,:)*solver%Y(i,:)-intfun(:)) / solver%D(i,:)

            check_inner_boundary: IF (solver%V(i,1).NE.0.0_r8) THEN
               solver%DY(i,1)           = (solver%W(i,1) - solver%U(i,1)*solver%Y(i,1)) / solver%V(i,1)
            END IF check_inner_boundary
            check_outer_boundary: IF (solver%V(i,2).NE.0.0_r8) THEN
               solver%DY(i,solver%NRHO) = (solver%W(i,2) - solver%U(i,2)*solver%Y(i,solver%NRHO)) / solver%V(i,2)
            END IF check_outer_boundary

         END IF derivative_recalculate

      END DO check_derivative


!   CHECK errors on solution and its derivative:

      IF (debug) THEN
         i = solver%NRHO/2
         WRITE(*,*) 'Solution maximum fractional error (mid radius) = ',  &
              maxval(abs(  &
                        solver%c(1:solver%ndim,:)*solver%a(1:solver%ndim,:)*solver%y(1:solver%ndim,:) &
                      - solver%c(1:solver%ndim,:)*solver%b(1:solver%ndim,:)*solver%ym(1:solver%ndim,:) &
                      - solver%h*solver%d(1:solver%ndim,:)*solver%dy(1:solver%ndim,:) &
                      + solver%h*solver%e(1:solver%ndim,:)*solver%y(1:solver%ndim,:) &
                      - solver%h*solver%c(1:solver%ndim,:)*solver%f(1:solver%ndim,:) &
                      + solver%h*solver%c(1:solver%ndim,:)*solver%g(1:solver%ndim,:)*solver%y(1:solver%ndim,:) )) &
                      /  &
               max( tiny(1.0_r8), &
               maxval(abs( solver%c(1:solver%ndim,:)*solver%a(1:solver%ndim,:)*solver%y(1:solver%ndim,:) )), &
               maxval(abs( solver%c(1:solver%ndim,:)*solver%b(1:solver%ndim,:)*solver%ym(1:solver%ndim,:) )), &
               maxval(abs( solver%h*solver%d(1:solver%ndim,:)*solver%dy(1:solver%ndim,:) )), &
               maxval(abs( solver%h*solver%e(1:solver%ndim,:)*solver%y(1:solver%ndim,:) )), &
               maxval(abs( solver%h*solver%c(1:solver%ndim,:)*solver%f(1:solver%ndim,:) )), &
               maxval(abs( solver%h*solver%c(1:solver%ndim,:)*solver%g(1:solver%ndim,:)*solver%y(1:solver%ndim,:) ))  &
               )


         WRITE(*,*) 'Solution maximum fractional error (Boundary 1) = ',  &
              maxval(abs(  &
                        solver%v(1:solver%ndim,1)*solver%dy(1:solver%ndim,1) &
                      + solver%u(1:solver%ndim,1)*solver%y(1:solver%ndim,1) &
                      - solver%w(1:solver%ndim,1))) &
                      /  &
               max( tiny(1.0_r8), &
               maxval(abs( solver%v(1:solver%ndim,1)*solver%dy(1:solver%ndim,1) )), &
               maxval(abs( solver%u(1:solver%ndim,1)*solver%y(1:solver%ndim,1) )), &
               maxval(abs( solver%w(1:solver%ndim,1) )))


         WRITE(*,*) 'Solution maximum fractional error (Boundary 2) = ',  &
              maxval(abs(  &
                        solver%v(1:solver%ndim,2)*solver%dy(1:solver%ndim,solver%nrho) &
                      + solver%u(1:solver%ndim,2)*solver%y(1:solver%ndim,solver%nrho) &
                      - solver%w(1:solver%ndim,2))) &
                      /  &
               max( tiny(1.0_r8), &
               maxval(abs( solver%v(1:solver%ndim,2)*solver%dy(1:solver%ndim,solver%nrho) )), &
               maxval(abs( solver%u(1:solver%ndim,2)*solver%y(1:solver%ndim,solver%nrho) )), &
               maxval(abs( solver%w(1:solver%ndim,2) )))

      END IF

      RETURN

      END SUBROUTINE solution_interface

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