text:
!     .................................................
!             ____  _       _   ____  _____   _
!            |  _ \| |     |_| |  _ \|  ___| |_|
!            | |_) | |___   _  | |_) | |___   _
!            |  _ /|  _  | | | |  _ /|___  | | |
!            | |   | | | | | | | |    ___| | | |
!            |_|   |_| |_| |_| |_|   |_____| |_|
!     .................................................
!     PhiPsi:     a general-purpose computational
!                 mechanics program written in Fortran.
!     Website:    http://phipsi.top
!     Author:     Shi Fang from Huaiyin Institute of
!                 Technology, HuaiAn, JiangSu, China
!     Contact me: shifang@hyit.edu.cn
!     ------------------------------------------------
!     Please cite the following papers:
!     (1)Shi F, Wang X L, Liu C, Liu H, Wu H A. An
!        XFEM-based method with reduction technique
!        for modeling hydraulic fracture propagation
!        in formations containing frictional natural
!        fractures. Engineering Fracture Mechanics,
!        2017, 173: 64-90.
!     (2)Shi F, Wang X L, Liu C, Liu H, Wu H A. A
!        coupled extended finite element approach
!        for modeling hydraulic fracturing in
!        consideration of proppant. Journal of
!        Natural Gas Science and Engineering, 2016,
!        33: 885-897.
!     (3)Shi F, Wang X L, Liu C, Liu H, Wu H A. An
!        XFEM-based numerical model to calculate
!        conductivity of propped fracture considering
!        proppant transport, embedment and crushing.
!        Journal of Petroleum Science and Engineering,
!        2018, 167: 615-626..

      subroutine Cal_Ele_Num_by_Coors(x,y,OUT_Elem)
C     get element number using coordinates.

      use Global_Model

      implicit none

      integer i
      double precision x,y
      integer,intent(out)::OUT_Elem
      double precision c_x_max,c_x_min,c_y_max,c_y_min

      integer c_count,Potent_Elem(5)
      logical Yes_In, Yes_On
      double precision xpol(5),ypol(5)

      c_count = 0

      !get the potential elements number.
      do i=1,Num_Elem
          c_x_max = maxval(G_X_NODES(i,:))
          c_x_min = minval(G_X_NODES(i,:))
          c_y_max = maxval(G_Y_NODES(i,:))
          c_y_min = minval(G_Y_NODES(i,:))
          if   ((x.le.c_x_max).and.(x.ge.c_x_min).
     &     and. (y.le.c_y_max).and.(y.ge.c_y_min))then
              c_count = c_count +1
              Potent_Elem(c_count) = i
          end if
      end do

      !Looking for the correct element.
      do i =1,c_count
          xpol(1:4) = G_X_NODES(Potent_Elem(i),1:4)
          ypol(1:4) = G_Y_NODES(Potent_Elem(i),1:4)
          xpol(5) = G_X_NODES(Potent_Elem(i),1)
          ypol(5) = G_Y_NODES(Potent_Elem(i),1)
          call Tool_Yes_In_Poly(x,y,xpol,ypol,5,Yes_In)
          call Tool_Yes_On_Poly(x,y,xpol,ypol,5,Yes_On)
          if ((Yes_In.eqv..True.).or.(Yes_On.eqv..True.))then
              OUT_Elem = Potent_Elem(i)
              exit
          end if
      end do

      return
      end SUBROUTINE Cal_Ele_Num_by_Coors

ext: f


text:
      SUBROUTINE WRTMSG (FILEX)
C
      EXTERNAL        LSHIFT,RSHIFT,ANDF,ORF,COMPLF
      INTEGER         FILE,FILEX,TITLE,TTLSAV(32,6),COUNT,LST(50),
     1                FOR(100),RET,EJECT,REW,BLANK,FORMAX,MASK1(5),
     2                MASK2(5),POS,ANDF,ORF,RSHIFT,COMPLF,SYSX
CWKBI
      CHARACTER*1     FORMT(400)
      COMMON /OUTPUT/ TITLE(32,6)
      COMMON /MACHIN/ MACH
      COMMON /SYSTEM/ SYSX(41)
CWKBI
      EQUIVALENCE     (FORMT, FOR)
      EQUIVALENCE  (XLST,LST)
      EQUIVALENCE     (SYSX( 2),MO   ), (SYSX( 9),MAXLIN),
     1                (SYSX(12),COUNT), (SYSX(39),NBPC  ),
     2                (SYSX(40),NBPW ), (SYSX(41),NCPW  )
      DATA    LSTMAX, REW,FORMAX,BLANK/ 50,1,100,4H     /
C
      N2CPW  = NCPW/2
      N2CPW1 = N2CPW - 1
      NBPC2  = 2*NBPC
      MASK1(1) = RSHIFT(COMPLF(0),NBPC2)
      MASK2(1) = COMPLF(MASK1(1))
      DO 10 I  = 2,N2CPW
      MASK1(I) = ORF(MASK2(1),RSHIFT(MASK1(I-1),NBPC2))
      MASK2(I) = COMPLF(MASK1(I))
   10 CONTINUE
      FILE = FILEX
C
      DO 20 J = 1,6
      DO 20 I = 1,32
      TTLSAV(I,J) = TITLE(I,J)
   20 CONTINUE
C
   30 COUNT = MAXLIN
   40 CALL READ (*500,*30,FILE,N,1,0,NF)
      IF (N) 100,130,110
C
C     A TITLE OR SUBTITLE FOLLOWS.
C
  100 N = -N
      IF (N .LE. 6) CALL FREAD (FILE,TITLE(1,N),32,0)
      IF (N .GT. 6) CALL FREAD (FILE,0,-32,0)
      GO TO 30
C
C     A MESSAGE FOLLOWS...N = NUMBER OF LIST ITEMS.
C
  110 IF (N .LE. LSTMAX) GO TO 120
      CALL FREAD (FILE,0,-N,0)
      GO TO 130
  120 IF (N .NE. 0) CALL FREAD (FILE,LST,N,0)
C
C     READ THE CORRESPONDING FORMAT...NF = SIZE OF THE FORMAT.
C
  130 CALL FREAD (FILE,NF,1,0)
      IF (NF) 140,150,160
  140 COUNT = COUNT - NF
      GO TO 130
  150 COUNT = MAXLIN
      GO TO 130
  160 IF (NF .LE. FORMAX) GO TO 170
      CALL FREAD (FILE,0,-NF,0)
      GO TO 30
  170 CALL FREAD (FILE,FOR,NF,0)
C
C     CONDENSE FOR ARRAY TO ACQUIRE CONTIGUOUS HOLLERITH STRINGS.
C
      IF (NCPW .EQ. 4) GO TO 300
      DO 290 I = 2,NF
      K1 = 1
      POS= 2*I - 1
      J  = (POS+N2CPW1)/N2CPW
      K2 = POS - N2CPW*(J-1)
      ASSIGN 200 TO RET
      GO TO 240
  200 CONTINUE
      K1 = 2
      IF (K2+1 .LE. N2CPW) GO TO 210
      K2 = 1
      J  = J + 1
      GO TO 220
  210 K2 = K2 + 1
  220 CONTINUE
      ASSIGN 230 TO RET
      GO TO 240
  230 CONTINUE
      GO TO 290
  240 IF (K2-K1) 250,260,270
  250 FOR(J) = ORF(ANDF(FOR(J),MASK1(K2)),
     1         LSHIFT(ANDF(FOR(I),MASK2(K1)),(NBPC2*(K1-K2))))
      GO TO 280
  260 FOR(J) = ORF(ANDF(FOR(J),MASK1(K2)),ANDF(FOR(I),MASK2(K1)))
      GO TO 280
  270 FOR(J) = ORF(ANDF(FOR(J),MASK1(K2)),
     1         RSHIFT(ANDF(FOR(I),MASK2(K1)),(NBPC2*(K2-K1))))
      GO TO 280
  280 CONTINUE
      GO TO RET, (200,230)
  290 CONTINUE
  300 CONTINUE
C
C     PRINT THE LINE
C
      IF (EJECT(1) .EQ. 0) GO TO 450
      DO 440 J = 4,6
      DO 410 I = 1,32
      IF (TITLE(I,J) .NE. BLANK) GO TO 420
  410 CONTINUE
      COUNT = COUNT - 1
      GO TO 440
  420 WRITE  (MO,430) (TITLE(I,J),I=1,32)
  430 FORMAT (2X,32A4)
  440 CONTINUE
      WRITE  (MO,430)
      COUNT = COUNT + 1
C
  450 IF(N.EQ.0 .AND. (MACH.EQ.5 .OR. MACH.EQ.12) )GO TO 470
      IF (MACH .EQ. 5 .OR. MACH .EQ. 12 ) GO TO 460
      CALL FORWRT ( FORMT, LST, N )
      GO TO 40
  460 WRITE (MO,FOR,ERR=465) (LST(J),J=1,N)
  465 CONTINUE
      GO TO 40
  470 WRITE (MO,FOR)
      GO TO 40
C
C     END OF MESSAGE FILE
C
  500 CALL CLOSE (FILE,REW)
      DO 510 J = 1,6
      DO 510 I = 1,32
      TITLE(I,J) = TTLSAV(I,J)
  510 CONTINUE
      RETURN
      END


text:
module grid_generic

  use core_lib, only : sp, dp, hid_t, warn, error
  use mpi_hdf5_io, only : mp_create_group
  use mpi_core, only : main_process

  use grid_io, only : write_grid_3d, write_grid_4d, write_grid_5d
  use grid_geometry, only : geo
  use grid_physics, only : n_photons, last_photon_id, specific_energy_sum, specific_energy_sum_nu, specific_energy, density, density_original
  use settings, only : output_n_photons, output_specific_energy, output_density, output_density_diff, physics_io_type, compute_isrf

  use dust_main, only: d

  implicit none
  save

  private
  public :: grid_reset_energy
  public :: output_grid

contains

  subroutine grid_reset_energy()
    implicit none
    if(allocated(n_photons)) n_photons = 0
    if(allocated(last_photon_id)) last_photon_id = 0
    if(allocated(specific_energy_sum)) specific_energy_sum = 0._dp
    if(allocated(specific_energy_sum_nu)) specific_energy_sum_nu = 0._dp
  end subroutine grid_reset_energy

  subroutine output_grid(group, iter, n_iter)

    implicit none

    integer(hid_t),intent(in) :: group
    integer,intent(in) :: iter, n_iter
    integer :: n_cells, n_dust, n_isrf_lam
    integer :: i,j,k
    real, dimension(d(1)%n_nu) :: energy_frequency_bins

    if(main_process()) write(*,'(" [output_grid] outputting grid arrays for iteration")')

    ! NUMBER OF PHOTONS IN EACH CELL

    if(trim(output_n_photons)=='all' .or. (trim(output_n_photons)=='last'.and.iter==n_iter)) then
       if(allocated(n_photons)) then
          call write_grid_3d(group, 'n_photons', n_photons, geo)
       else
          call warn("output_grid","n_photons array is not allocated")
       end if
    end if

    ! ENERGY/PATH LENGTHS

    if(trim(output_specific_energy)=='all' .or. (trim(output_specific_energy)=='last'.and.iter==n_iter)) then
       if(allocated(specific_energy)) then
          select case(physics_io_type)
          case(sp)
             call write_grid_4d(group, 'specific_energy', real(specific_energy, sp), geo)
          case(dp)
             call write_grid_4d(group, 'specific_energy', real(specific_energy, dp), geo)
          case default
             call error("output_grid","unexpected value of physics_io_type (should be sp or dp)")
          end select
       else
          call warn("output_grid","specific_energy array is not allocated")
       end if
    end if

    ! WRITE THE ISRF
    if (compute_isrf) then

       n_cells = size(specific_energy_sum_nu, 1)
       n_dust = size(specific_energy_sum_nu, 2)
       n_isrf_lam = size(specific_energy_sum_nu,3)

       do i=1,d(1)%n_nu
          energy_frequency_bins(i) = d(1)%nu(i)
       end do

       if(trim(output_specific_energy)=='all' .or. (trim(output_specific_energy)=='last'.and.iter==n_iter)) then
          !if(allocated(energy_frequency_bins)) then
          call write_grid_3d(group, 'ISRF_frequency_bins',energy_frequency_bins, geo)
       else
          call warn("output_grid","energy_frequency bins [ISRF wavelengths] array is not allocated")
          !end if
       end if

       if(trim(output_specific_energy)=='all' .or. (trim(output_specific_energy)=='last'.and.iter==n_iter)) then
          if(allocated(specific_energy_sum_nu)) then

             select case(physics_io_type)

             case(sp)
                call write_grid_5d(group, 'specific_energy_sum_nu', real(specific_energy_sum_nu, sp), geo)
             case(dp)
                call write_grid_5d(group, 'specific_energy_sum_nu', real(specific_energy_sum_nu, dp), geo)
             case default
                call error("output_grid","unexpected value of physics_io_type (should be sp or dp)")
             end select
          else
             call warn("output_grid","specific_energy_sum_nu array is not allocated")
          end if
       end if

    endif

    ! DENSITY

    if(trim(output_density)=='all' .or. (trim(output_density)=='last'.and.iter==n_iter)) then
       if(allocated(density)) then
          select case(physics_io_type)
          case(sp)
             call write_grid_4d(group, 'density', real(density, sp), geo)
          case(dp)
             call write_grid_4d(group, 'density', real(density, dp), geo)
          case default
             call error("output_grid","unexpected value of physics_io_type (should be sp or dp)")
          end select
       else
          call warn("output_grid","density array is not allocated")
       end if
    end if

    ! DENSITY DIFFERENCE

    if(trim(output_density_diff)=='all' .or. (trim(output_density_diff)=='last'.and.iter==n_iter)) then
       if(allocated(density).and.allocated(density_original)) then
          select case(physics_io_type)
          case(sp)
             call write_grid_4d(group, 'density_diff', real(density - density_original, sp), geo)
          case(dp)
             call write_grid_4d(group, 'density_diff', real(density - density_original, dp), geo)
          case default
             call error("output_grid","unexpected value of physics_io_type (should be sp or dp)")
          end select
       else
          if(.not.allocated(density)) call warn("output_grid","density array is not allocated")
          if(.not.allocated(density_original)) call warn("output_grid","density_original array is not allocated")
       end if
    end if

  end subroutine output_grid

end module grid_generic

max_stars_repo_name: dnarayanan/hyperion

max_issues_repo_name: dnarayanan/hyperion

max_forks_repo_name: dnarayanan/hyperion


text:
  subroutine func(argc, argv, output)

    implicit none

    integer, intent(in) :: argc
    type (string_array) :: argv
    double precision, intent(out) :: output

max_stars_repo_name: jensteff/Swift-Prov

max_issues_repo_name: jensteff/Swift-Prov

max_forks_repo_name: jensteff/Swift-Prov


text:
module multigrid_2d

  use shared_data

  implicit none

  private

  public :: solve_mg_2d

  contains

  subroutine solve_mg_2d(mgphi, f, tol) ! outside facing interface

    real(num), dimension(:,:), allocatable, intent(inout) :: mgphi
    real(num), dimension(1:nx,1:ny), intent(in) :: f
    real(num), intent(in) :: tol
    ! add an argument to select type of MG solver when more than one exists

    call mg_2level(mgphi, f, tol)

    print *, 'called solve_mg_2d'
  end subroutine solve_mg_2d

  ! different solvers

  subroutine mg_2level(mgphi, f, tol)

    real(num), dimension(:,:), allocatable, intent(inout) :: mgphi
    real(num), dimension(1:nx,1:ny), intent(in) :: f
    real(num), intent(in) :: tol

    integer :: nsteps
    integer :: ncoarse = 3
    integer :: nfine = 50
    integer :: c

    real(num), dimension(1:nx,1:ny) :: residual
    real(num) :: L2, L2_old
    real(num) :: Lap
    real(num) :: start
    real(num) :: finish

    ! coarse stuff
    real(num), dimension(1:nx/2,1:ny/2) :: residual_c
    real(num), dimension(:,:), allocatable :: epsi
    integer :: ixc, iyc
    if (.not.allocated(epsi)) allocate(epsi(0:nx/2+1,0:ny/2+1))

    print *, 'called mg_2level'
    call cpu_time (start)

    ! main cycle
    nsteps = 0
    L2_old = 1e6_num
    do
      nsteps = nsteps + 1
      ! ncoarse iterations of GS relaxation on finest grid
      do c = 1, ncoarse

        call bcs_fine(mgphi)

        ! redblack, odd
        do iy = 1, ny
        do ix = 1, nx
          if (modulo(ix+iy,2) == 1) then
            mgphi(ix,iy) = 0.25_num * ( &
              & mgphi(ix+1,iy) + mgphi(ix-1,iy) + mgphi(ix,iy+1) + mgphi(ix,iy-1) &
              - dx**2 * f(ix,iy) )
          endif
        end do
        end do
        ! even iteration
        do iy = 1, ny
        do ix = 1, nx
          if (modulo(ix+iy,2) == 0) then
            mgphi(ix,iy) = 0.25_num * ( &
              & mgphi(ix+1,iy) + mgphi(ix-1,iy) + mgphi(ix,iy+1) + mgphi(ix,iy-1) &
              - dx**2 * f(ix,iy) )
          endif
        end do
        end do
      enddo

      ! see if converged to tolerance and exit if satisfied

      call bcs_fine(mgphi)

      L2 = 0.0_num
      do iy = 1, ny
      do ix = 1, nx
        Lap = (mgphi(ix+1,iy) - 2.0_num*mgphi(ix,iy) + mgphi(ix-1,iy)) / dx**2 + &
            & (mgphi(ix,iy+1) - 2.0_num*mgphi(ix,iy) + mgphi(ix,iy-1)) / dy**2
        residual(ix,iy) = Lap - f(ix,iy)
      enddo
      enddo
      L2 = sqrt(sum(abs(residual)**2)/real(nx*ny,num))

      if (abs(L2-L2_old) <= tol) exit
      L2_old = L2

! skip for test purposes
!GOTO 9999 ! naughty

      ! restrict the residue down to a grid with half the resolution in each
      ! dimension
      iy=  1
      do iyc = 1, ny/2
      ix = 1
      do ixc = 1, nx/2
        residual_c(ixc,iyc) = 0.25_num * (residual(ix,iy) + residual(ix+1,iy) &
          & + residual(ix,iy+1) + residual(ix+1,iy+1) )
        ix = ix + 2
      enddo
      iy = iy + 2
      enddo

      ! perform smoothing on the coarser grid for Lap * epsi = residual_c

      do c = 1, nfine

        call bcs_coarse(epsi)

        ! redblack, odd
        do iy = 1, ny/2
        do ix = 1, nx/2
          if (modulo(ix+iy,2) == 1) then
            epsi(ix,iy) = 0.25_num * ( &
              & epsi(ix+1,iy) + epsi(ix-1,iy) + epsi(ix,iy+1) + epsi(ix,iy-1) &
              - (2.0_num*dx)**2 * residual_c(ix,iy) )
          endif
        end do
        end do
        ! even iteration
        do iy = 1, ny/2
        do ix = 1, nx/2
          if (modulo(ix+iy,2) == 0) then
            epsi(ix,iy) = 0.25_num * ( &
              & epsi(ix+1,iy) + epsi(ix-1,iy) + epsi(ix,iy+1) + epsi(ix,iy-1) &
              - (2.0_num*dx)**2 * residual_c(ix,iy) )
          endif
        end do
        end do

      enddo

      ! prolongate via direct injection back into fine residual

      iy=  1
      do iyc = 1, ny/2
      ix = 1
      do ixc = 1, nx/2
        mgphi(ix,iy) = mgphi(ix,iy) - epsi(ixc,iyc)
        mgphi(ix+1,iy) = mgphi(ix+1,iy) - epsi(ixc,iyc)
        mgphi(ix,iy+1) = mgphi(ix,iy+1) - epsi(ixc,iyc)
        mgphi(ix+1,iy+1) = mgphi(ix+1,iy+1) - epsi(ixc,iyc)
        ix = ix + 2
      enddo
      iy = iy + 2
      enddo
9999 CONTINUE ! goto is naughty

    enddo ! exit main clycle

    call cpu_time (finish)

    print *, 'nsteps', nsteps
    print *,' L2', L2
    print '(" ****** cpu_time: ",f20.3," seconds.")',finish-start

  end subroutine mg_2level

  subroutine bcs_fine(mgphi)

    real(num), dimension(:,:), allocatable, intent(inout) :: mgphi

    mgphi(0,:) = mgphi(nx,:)
    mgphi(nx+1,:) = mgphi(1,:)
    mgphi(:,0) = mgphi(:,ny)
    mgphi(:,ny+1) = mgphi(:,1)

  end subroutine bcs_fine

  subroutine bcs_coarse(epsi)

    real(num), dimension(:,:), allocatable, intent(inout) :: epsi

    epsi(0,:) = epsi(nx/2,:)
    epsi(nx/2+1,:) = epsi(1,:)
    epsi(:,0) = epsi(:,ny/2)
    epsi(:,ny/2+1) = epsi(:,1)

  end subroutine bcs_coarse

end module multigrid_2d

max_stars_repo_name: JOThurgood/SimpleCFD

max_issues_repo_name: JOThurgood/SimpleCFD

max_forks_repo_name: JOThurgood/SimpleCFD