| |
| c--------------------------------------------------------------------- |
| c--------------------------------------------------------------------- |
| |
| subroutine initialize |
| |
| c--------------------------------------------------------------------- |
| c--------------------------------------------------------------------- |
| |
| c--------------------------------------------------------------------- |
| c This subroutine initializes the field variable u using |
| c tri-linear transfinite interpolation of the boundary values |
| c--------------------------------------------------------------------- |
| |
| include 'header.h' |
| |
| integer c, i, j, k, m, ii, jj, kk, ix, iy, iz |
| double precision xi, eta, zeta, Pface(5,3,2), Pxi, Peta, |
| > Pzeta, temp(5) |
| |
| |
| c--------------------------------------------------------------------- |
| c Later (in compute_rhs) we compute 1/u for every element. A few of |
| c the corner elements are not used, but it convenient (and faster) |
| c to compute the whole thing with a simple loop. Make sure those |
| c values are nonzero by initializing the whole thing here. |
| c--------------------------------------------------------------------- |
| do c = 1, ncells |
| do kk = -1, IMAX |
| do jj = -1, IMAX |
| do ii = -1, IMAX |
| u(ii, jj, kk, 1, c) = 1.0 |
| u(ii, jj, kk, 2, c) = 0.0 |
| u(ii, jj, kk, 3, c) = 0.0 |
| u(ii, jj, kk, 4, c) = 0.0 |
| u(ii, jj, kk, 5, c) = 1.0 |
| end do |
| end do |
| end do |
| end do |
| |
| c--------------------------------------------------------------------- |
| c first store the "interpolated" values everywhere on the grid |
| c--------------------------------------------------------------------- |
| do c=1, ncells |
| kk = 0 |
| do k = cell_low(3,c), cell_high(3,c) |
| zeta = dble(k) * dnzm1 |
| jj = 0 |
| do j = cell_low(2,c), cell_high(2,c) |
| eta = dble(j) * dnym1 |
| ii = 0 |
| do i = cell_low(1,c), cell_high(1,c) |
| xi = dble(i) * dnxm1 |
| |
| do ix = 1, 2 |
| call exact_solution(dble(ix-1), eta, zeta, |
| > Pface(1,1,ix)) |
| end do |
| |
| do iy = 1, 2 |
| call exact_solution(xi, dble(iy-1) , zeta, |
| > Pface(1,2,iy)) |
| end do |
| |
| do iz = 1, 2 |
| call exact_solution(xi, eta, dble(iz-1), |
| > Pface(1,3,iz)) |
| end do |
| |
| do m = 1, 5 |
| Pxi = xi * Pface(m,1,2) + |
| > (1.0d0-xi) * Pface(m,1,1) |
| Peta = eta * Pface(m,2,2) + |
| > (1.0d0-eta) * Pface(m,2,1) |
| Pzeta = zeta * Pface(m,3,2) + |
| > (1.0d0-zeta) * Pface(m,3,1) |
| |
| u(ii,jj,kk,m,c) = Pxi + Peta + Pzeta - |
| > Pxi*Peta - Pxi*Pzeta - Peta*Pzeta + |
| > Pxi*Peta*Pzeta |
| |
| end do |
| ii = ii + 1 |
| end do |
| jj = jj + 1 |
| end do |
| kk = kk+1 |
| end do |
| end do |
| |
| c--------------------------------------------------------------------- |
| c now store the exact values on the boundaries |
| c--------------------------------------------------------------------- |
| |
| c--------------------------------------------------------------------- |
| c west face |
| c--------------------------------------------------------------------- |
| c = slice(1,1) |
| ii = 0 |
| xi = 0.0d0 |
| kk = 0 |
| do k = cell_low(3,c), cell_high(3,c) |
| zeta = dble(k) * dnzm1 |
| jj = 0 |
| do j = cell_low(2,c), cell_high(2,c) |
| eta = dble(j) * dnym1 |
| call exact_solution(xi, eta, zeta, temp) |
| do m = 1, 5 |
| u(ii,jj,kk,m,c) = temp(m) |
| end do |
| jj = jj + 1 |
| end do |
| kk = kk + 1 |
| end do |
| |
| c--------------------------------------------------------------------- |
| c east face |
| c--------------------------------------------------------------------- |
| c = slice(1,ncells) |
| ii = cell_size(1,c)-1 |
| xi = 1.0d0 |
| kk = 0 |
| do k = cell_low(3,c), cell_high(3,c) |
| zeta = dble(k) * dnzm1 |
| jj = 0 |
| do j = cell_low(2,c), cell_high(2,c) |
| eta = dble(j) * dnym1 |
| call exact_solution(xi, eta, zeta, temp) |
| do m = 1, 5 |
| u(ii,jj,kk,m,c) = temp(m) |
| end do |
| jj = jj + 1 |
| end do |
| kk = kk + 1 |
| end do |
| |
| c--------------------------------------------------------------------- |
| c south face |
| c--------------------------------------------------------------------- |
| c = slice(2,1) |
| jj = 0 |
| eta = 0.0d0 |
| kk = 0 |
| do k = cell_low(3,c), cell_high(3,c) |
| zeta = dble(k) * dnzm1 |
| ii = 0 |
| do i = cell_low(1,c), cell_high(1,c) |
| xi = dble(i) * dnxm1 |
| call exact_solution(xi, eta, zeta, temp) |
| do m = 1, 5 |
| u(ii,jj,kk,m,c) = temp(m) |
| end do |
| ii = ii + 1 |
| end do |
| kk = kk + 1 |
| end do |
| |
| |
| c--------------------------------------------------------------------- |
| c north face |
| c--------------------------------------------------------------------- |
| c = slice(2,ncells) |
| jj = cell_size(2,c)-1 |
| eta = 1.0d0 |
| kk = 0 |
| do k = cell_low(3,c), cell_high(3,c) |
| zeta = dble(k) * dnzm1 |
| ii = 0 |
| do i = cell_low(1,c), cell_high(1,c) |
| xi = dble(i) * dnxm1 |
| call exact_solution(xi, eta, zeta, temp) |
| do m = 1, 5 |
| u(ii,jj,kk,m,c) = temp(m) |
| end do |
| ii = ii + 1 |
| end do |
| kk = kk + 1 |
| end do |
| |
| c--------------------------------------------------------------------- |
| c bottom face |
| c--------------------------------------------------------------------- |
| c = slice(3,1) |
| kk = 0 |
| zeta = 0.0d0 |
| jj = 0 |
| do j = cell_low(2,c), cell_high(2,c) |
| eta = dble(j) * dnym1 |
| ii = 0 |
| do i =cell_low(1,c), cell_high(1,c) |
| xi = dble(i) *dnxm1 |
| call exact_solution(xi, eta, zeta, temp) |
| do m = 1, 5 |
| u(ii,jj,kk,m,c) = temp(m) |
| end do |
| ii = ii + 1 |
| end do |
| jj = jj + 1 |
| end do |
| |
| c--------------------------------------------------------------------- |
| c top face |
| c--------------------------------------------------------------------- |
| c = slice(3,ncells) |
| kk = cell_size(3,c)-1 |
| zeta = 1.0d0 |
| jj = 0 |
| do j = cell_low(2,c), cell_high(2,c) |
| eta = dble(j) * dnym1 |
| ii = 0 |
| do i =cell_low(1,c), cell_high(1,c) |
| xi = dble(i) * dnxm1 |
| call exact_solution(xi, eta, zeta, temp) |
| do m = 1, 5 |
| u(ii,jj,kk,m,c) = temp(m) |
| end do |
| ii = ii + 1 |
| end do |
| jj = jj + 1 |
| end do |
| |
| return |
| end |
| |
| |
| subroutine lhsinit |
| |
| include 'header.h' |
| |
| integer i, j, k, d, c, n |
| |
| c--------------------------------------------------------------------- |
| c loop over all cells |
| c--------------------------------------------------------------------- |
| do c = 1, ncells |
| |
| c--------------------------------------------------------------------- |
| c first, initialize the start and end arrays |
| c--------------------------------------------------------------------- |
| do d = 1, 3 |
| if (cell_coord(d,c) .eq. 1) then |
| start(d,c) = 1 |
| else |
| start(d,c) = 0 |
| endif |
| if (cell_coord(d,c) .eq. ncells) then |
| end(d,c) = 1 |
| else |
| end(d,c) = 0 |
| endif |
| end do |
| |
| c--------------------------------------------------------------------- |
| c zap the whole left hand side for starters |
| c--------------------------------------------------------------------- |
| do n = 1, 15 |
| do k = 0, cell_size(3,c)-1 |
| do j = 0, cell_size(2,c)-1 |
| do i = 0, cell_size(1,c)-1 |
| lhs(i,j,k,n,c) = 0.0d0 |
| end do |
| end do |
| end do |
| end do |
| |
| c--------------------------------------------------------------------- |
| c next, set all diagonal values to 1. This is overkill, but convenient |
| c--------------------------------------------------------------------- |
| do n = 1, 3 |
| do k = 0, cell_size(3,c)-1 |
| do j = 0, cell_size(2,c)-1 |
| do i = 0, cell_size(1,c)-1 |
| lhs(i,j,k,5*n-2,c) = 1.0d0 |
| end do |
| end do |
| end do |
| end do |
| |
| end do |
| |
| return |
| end |
| |
| |
| |
