generate_chunk_kernel.f90

!Crown Copyright 2014 AWE.
!
! This file is part of TeaLeaf.
!
! TeaLeaf is free software: you can redistribute it and/or modify it under 
! the terms of the GNU General Public License as published by the 
! Free Software Foundation, either version 3 of the License, or (at your option) 
! any later version.
!
! TeaLeaf is distributed in the hope that it will be useful, but 
! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 
! FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 
! details.
!
! You should have received a copy of the GNU General Public License along with 
! TeaLeaf. If not, see http://www.gnu.org/licenses/.

!>  @brief Fortran mesh chunk generator
!>  @author David Beckingsale, Wayne Gaudin
!>  @details Generates the field data on a mesh chunk based on the user specified
!>  input for the states.
!>
!>  Note that state one is always used as the background state, which is then
!>  overwritten by further state definitions.

MODULE generate_chunk_kernel_module

CONTAINS

SUBROUTINE generate_chunk_kernel(x_min,x_max,y_min,y_max, &
                                 vertexx,                 &
                                 vertexy,                 &
                                 cellx,                   &
                                 celly,                   &
                                 density,                 &
                                 energy0,                 &
                                 u0,                      &
                                 number_of_states,        &
                                 state_density,           &
                                 state_energy,            &
                                 state_xmin,              &
                                 state_xmax,              &
                                 state_ymin,              &
                                 state_ymax,              &
                                 state_radius,            &
                                 state_geometry,          &
                                 g_rect,                  &
                                 g_circ,                  &
                                 g_point)

  IMPLICIT NONE

  INTEGER      :: x_min,x_max,y_min,y_max
  REAL(KIND=8), DIMENSION(x_min-2:x_max+3) :: vertexx
  REAL(KIND=8), DIMENSION(y_min-2:y_max+3) :: vertexy
  REAL(KIND=8), DIMENSION(y_min-2:y_max+3) :: vertexdy
  REAL(KIND=8), DIMENSION(x_min-2:x_max+2) :: cellx
  REAL(KIND=8), DIMENSION(y_min-2:y_max+2) :: celly
  REAL(KIND=8), DIMENSION(x_min-2:x_max+2,y_min-2:y_max+2) :: density,energy0
  REAL(KIND=8), DIMENSION(x_min-2:x_max+2,y_min-2:y_max+2) :: u0
  INTEGER      :: number_of_states
  REAL(KIND=8), DIMENSION(number_of_states) :: state_density
  REAL(KIND=8), DIMENSION(number_of_states) :: state_energy
  REAL(KIND=8), DIMENSION(number_of_states) :: state_xmin
  REAL(KIND=8), DIMENSION(number_of_states) :: state_xmax
  REAL(KIND=8), DIMENSION(number_of_states) :: state_ymin
  REAL(KIND=8), DIMENSION(number_of_states) :: state_ymax
  REAL(KIND=8), DIMENSION(number_of_states) :: state_radius
  INTEGER     , DIMENSION(number_of_states) :: state_geometry
  INTEGER      :: g_rect
  INTEGER      :: g_circ
  INTEGER      :: g_point

  REAL(KIND=8) :: radius,x_cent,y_cent
  INTEGER      :: state

  INTEGER      :: j,k,jt,kt

  ! State 1 is always the background state

!$OMP PARALLEL SHARED(x_cent,y_cent)
!$OMP DO
  DO k=y_min-2,y_max+2
    DO j=x_min-2,x_max+2
      energy0(j,k)=state_energy(1)
    ENDDO
  ENDDO
!$OMP END DO
!$OMP DO
  DO k=y_min-2,y_max+2
    DO j=x_min-2,x_max+2
      density(j,k)=state_density(1)
    ENDDO
  ENDDO
!$OMP END DO

  DO state=2,number_of_states

    x_cent=state_xmin(state)
    y_cent=state_ymin(state)

!$OMP DO PRIVATE(radius,jt,kt)
    DO k=y_min-2,y_max+2
      DO j=x_min-2,x_max+2
        IF(state_geometry(state).EQ.g_rect ) THEN
          IF(vertexx(j+1).GE.state_xmin(state).AND.vertexx(j).LT.state_xmax(state)) THEN
            IF(vertexy(k+1).GE.state_ymin(state).AND.vertexy(k).LT.state_ymax(state)) THEN
              energy0(j,k)=state_energy(state)
              density(j,k)=state_density(state)
            ENDIF
          ENDIF
        ELSEIF(state_geometry(state).EQ.g_circ ) THEN
          radius=SQRT((cellx(j)-x_cent)*(cellx(j)-x_cent)+(celly(k)-y_cent)*(celly(k)-y_cent))
          IF(radius.LE.state_radius(state))THEN
            energy0(j,k)=state_energy(state)
            density(j,k)=state_density(state)
          ENDIF
        ELSEIF(state_geometry(state).EQ.g_point) THEN
          IF(vertexx(j).EQ.x_cent .AND. vertexy(k).EQ.y_cent) THEN
            energy0(j,k)=state_energy(state)
            density(j,k)=state_density(state)
          ENDIF
        ENDIF
      ENDDO
    ENDDO
!$OMP END DO

  ENDDO

!$OMP DO 
  DO k=y_min-1, y_max+1
    DO j=x_min-1, x_max+1
      u0(j,k) =  energy0(j,k) * density(j,k)
    ENDDO
  ENDDO
!$OMP END DO

!$OMP END PARALLEL

END SUBROUTINE generate_chunk_kernel

END MODULE generate_chunk_kernel_module