stag/mkgrid.f90 at master · sscollis/stag · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
program mkgrid

! Simple two-dimensional grid generator
!
! S. Scott Collis
! Rice University
! MEMS, MS 321
! Houston, TX 77005
! 713-348-3617
! collis@rice.edu

implicit none

integer :: nx, ny, i, j
real, allocatable :: x(:), y(:), xi(:), eta(:)

!real :: xmin=-32, xmax=32, ymin=0, ymax=32
!real :: xmin=-20, xmax=20, ymin=0, ymax=30
real :: xmin=-20, xmax=20, ymin=0, ymax=40
real :: dx, dy, dxi, deta
real :: xs, ys, aa, bb

write(*,"('Enter number of cells (nx, ny) ==> ',$)")
read(*,*) nx, ny
write(*,"('Enter Xs, Ys (Xs = Ys = 0 makes a uniform mesh) ==> ',$)")
read(*,*) xs, ys

allocate( x(nx+2), y(ny+2), xi(nx+2), eta(ny+2) )

dx = (xmax - xmin) / real(nx)
dy = (ymax - ymin) / real(ny)

if (xs.eq.0) then
  do i = 1, nx+1
    x(i) = xmin + dx * (i-1)
  end do
  x(nx+2) = x(nx+1)
else
  dxi = 1.0 / float(nx)
  write(*,"('Algebraic grid in x')")
  aa = (xmax-xmin) * xs / ( (xmax-xmin) - 2.0 * xs )
  bb = 1.0 + aa / (xmax-xmin)

  do i = 1, nx+1
    xi(i) = -0.5 + (i-1) * dxi
    x(i)  = 4.0 * ( aa * xi(i) / (bb - abs(xi(i))) )
  end do
  x(nx+2) = x(nx+1)
  do i = 1, nx+1
    write(11,"(2(1pe13.6,1x))") xi(i), x(i)
  end do
end if

if (ys.eq.0) then
  do j = 1, ny+1
    y(j) = ymin + dy * (j-1)
  end do
  y(ny+2) = y(ny+1)
else
  deta = 1.0 / float(ny)
  write(*,"('Algebraic grid in y')")
  aa = ymax * ys / ( ymax - 2.0 * ys )
  bb = 1.0 + aa / ymax
  do j = 1, ny+1
    eta(j) = (j-1) * deta
    y(j)   = aa * eta(j) / (bb - eta(j))
  end do
  y(ny+2) = y(ny+1)
  do j = 1, ny+1
    write(12,"(2(1pe13.6,1x))") eta(j), y(j)
  end do
end if

open(10,file='grid.in',form='formatted',status='unknown')
write(10,*) nx+2, ny+2
write(10,*) x
write(10,*) y
close(10)

stop
end program mkgrid