X-Git-Url: http://info.iut-bm.univ-fcomte.fr/pub/gitweb/simgrid.git/blobdiff_plain/d12142c9dc382512752775d0cca043bb1aef3f86..4c74bb7b6f2398da81ce462cbdfd9c5a77ffa683:/examples/smpi/NAS/BT/initialize.f diff --git a/examples/smpi/NAS/BT/initialize.f b/examples/smpi/NAS/BT/initialize.f new file mode 100644 index 0000000000..274cdb1899 --- /dev/null +++ b/examples/smpi/NAS/BT/initialize.f @@ -0,0 +1,308 @@ +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, KMAX + do jj = -1, JMAX + do ii = -1, IMAX + do m = 1, 5 + u(m, ii, jj, kk, c) = 1.0 + end do + end do + end do + end do + end do +c--------------------------------------------------------------------- + + + +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)) + enddo + + do iy = 1, 2 + call exact_solution(xi, dble(iy-1) , zeta, + > Pface(1,2,iy)) + enddo + + do iz = 1, 2 + call exact_solution(xi, eta, dble(iz-1), + > Pface(1,3,iz)) + enddo + + 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(m,ii,jj,kk,c) = Pxi + Peta + Pzeta - + > Pxi*Peta - Pxi*Pzeta - Peta*Pzeta + + > Pxi*Peta*Pzeta + + enddo + ii = ii + 1 + enddo + jj = jj + 1 + enddo + kk = kk+1 + enddo + enddo + +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(m,ii,jj,kk,c) = temp(m) + enddo + jj = jj + 1 + enddo + kk = kk + 1 + enddo + +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(m,ii,jj,kk,c) = temp(m) + enddo + jj = jj + 1 + enddo + kk = kk + 1 + enddo + +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(m,ii,jj,kk,c) = temp(m) + enddo + ii = ii + 1 + enddo + kk = kk + 1 + enddo + + +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(m,ii,jj,kk,c) = temp(m) + enddo + ii = ii + 1 + enddo + kk = kk + 1 + enddo + +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(m,ii,jj,kk,c) = temp(m) + enddo + ii = ii + 1 + enddo + jj = jj + 1 + enddo + +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(m,ii,jj,kk,c) = temp(m) + enddo + ii = ii + 1 + enddo + jj = jj + 1 + enddo + + return + end + + +c--------------------------------------------------------------------- +c--------------------------------------------------------------------- + + subroutine lhsinit + +c--------------------------------------------------------------------- +c--------------------------------------------------------------------- + + include 'header.h' + + integer i, j, k, d, c, m, 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 + enddo + +c--------------------------------------------------------------------- +c zero the whole left hand side for starters +c--------------------------------------------------------------------- + do k = 0, cell_size(3,c)-1 + do j = 0, cell_size(2,c)-1 + do i = 0, cell_size(1,c)-1 + do m = 1,5 + do n = 1, 5 + lhsc(m,n,i,j,k,c) = 0.0d0 + enddo + enddo + enddo + enddo + enddo + + enddo + + return + end + + +c--------------------------------------------------------------------- +c--------------------------------------------------------------------- + + subroutine lhsabinit(lhsa, lhsb, size) + implicit none + + integer size + double precision lhsa(5, 5, -1:size), lhsb(5, 5, -1:size) + + integer i, m, n + +c--------------------------------------------------------------------- +c next, set all diagonal values to 1. This is overkill, but convenient +c--------------------------------------------------------------------- + do i = 0, size + do m = 1, 5 + do n = 1, 5 + lhsa(m,n,i) = 0.0d0 + lhsb(m,n,i) = 0.0d0 + enddo + lhsb(m,m,i) = 1.0d0 + enddo + enddo + + return + end + + +