Commits

Anonymous committed cff82c1

NS: Various fixes/updates.

Comments (0)

Files changed (6)

examples/mpi-navier-stokes/Makefile

 # dependencies
 #
 
-build/solutions.o: build/pfasst.o
 build/feval.o:     build/pfasst.o build/encap.o build/fftwpp.o
-build/hooks.o:     build/numpy.o build/solutions.o
+build/hooks.o:     build/numpy.o build/initial.o
+build/transfer.o:  build/feval.o
 build/main.o:      build/feval.o build/hooks.o build/transfer.o build/initial.o

examples/mpi-navier-stokes/npy2bov.py

+
+import sys
+import numpy as np
+
+for fname in sys.argv[1:]:
+    z = np.load(fname)
+    u = np.real(np.fft.ifftn(z[:,:,:,0]))
+    v = np.real(np.fft.ifftn(z[:,:,:,1]))
+    w = np.real(np.fft.ifftn(z[:,:,:,2]))
+
+    basename = fname.split('.')[0]
+    dat = basename + '.dat'
+    bov = basename + '.bov'
+
+    n = u.shape[0]
+
+    print 'converting:', fname
+
+    b = np.empty((3,) + u.shape)
+    b[0] = u
+    b[1] = v
+    b[2] = w
+    with open(dat, 'wb') as f:
+        f.write(b)
+
+    with open(bov, 'w') as f:
+        f.write('DATA_FILE: %s\n' % dat)
+        f.write('DATA_SIZE: %d %d %d\n' % (n, n, n))
+        f.write('DATA_COMPONENTS: 3\n')
+
+

examples/mpi-navier-stokes/src/feval.f90

 
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-  subroutine divergence(cptr, n1, n2, n3, u, div)
-    implicit none
-    type(c_ptr),               intent(in), value :: cptr
-    integer(c_int),            intent(in), value :: n1, n2, n3
-    complex(c_double_complex), intent(in)        :: u(n3, n2, n1, 3)
-    complex(c_double_complex), intent(out)       :: div(n3, n2, n1)
-
-    type(feval_t), pointer    :: fptr
-    integer                   :: i1, i2, i3
-    real(c_double)            :: kk(n1)
-
-    call c_f_pointer(cptr, fptr)
-    kk = fptr%k
-
-    !$omp parallel do private(i1, i2, i3)
-    do i1 = 1,  n1
-       do  i2 = 1,  n2
-          do  i3 = 1,  n3
-
-             ! phi = div(ustar)
-             div(i3, i2, i1) = &
-                  kk(i1) * (0.0d0,1.0d0) * u(i3, i2, i1, 1) + &
-                  kk(i2) * (0.0d0,1.0d0) * u(i3, i2, i1, 2) + &
-                  kk(i3) * (0.0d0,1.0d0) * u(i3, i2, i1, 3)
-
-          end do
-       end do
-    end do
-    !$omp end parallel do
-
-  end subroutine divergence
-
-  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
   subroutine project(cptr, n1, n2, n3, ustar, u)
     ! Project ustar to divergence free u
     !
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
   subroutine copy_for_convolve(y, ctx)
-    type(carray4), intent(in) :: y
-    type(feval_t),    intent(in) :: ctx
+    type(carray4), intent(in   ) :: y
+    type(feval_t), intent(in   ) :: ctx
 
-    complex(c_double), dimension(:,:,:), pointer :: &
+    complex(c_double_complex), dimension(:,:,:), pointer :: &
          u1, v1, w1, v2, w2, w3, uu, uv, uw, vv, vw, ww
 
     call c_f_pointer(ctx%u1, u1, [ ctx%n, ctx%n, ctx%n ])
     call c_f_pointer(ctx%vw, vw, [ ctx%n, ctx%n, ctx%n ])
     call c_f_pointer(ctx%ww, ww, [ ctx%n, ctx%n, ctx%n ])
 
-!xxx    !$omp parallel workshare
+    !$omp parallel workshare
     u1 = y%array(:,:,:,1)
     v1 = y%array(:,:,:,2)
     w1 = y%array(:,:,:,3)
     uw = y%array(:,:,:,1)
     vw = y%array(:,:,:,2)
     ww = y%array(:,:,:,3)
-!xxx    !$omp end parallel workshare
+    !$omp end parallel workshare
   end subroutine copy_for_convolve
 
   subroutine eval_f1(yptr, t, level, ctxptr, f1ptr)
     real(pfdp),  intent(in   )        :: t
     integer,     intent(in   )        :: level
 
-    complex(c_double), dimension(:,:,:), pointer :: uu, uv, uw, vv, vw, ww
+    complex(c_double_complex), dimension(:,:,:), pointer :: uu, uv, uw, vv, vw, ww
+    type(carray4)                                        :: u
 
     type(feval_t), pointer :: ctx
     type(carray4), pointer :: y, f1
     call c_f_pointer(f1ptr, f1)
     call c_f_pointer(ctxptr, ctx)
 
-    call copy_for_convolve(y, ctx)
+    call carray4_create(u, y%shape)
+    call project(ctxptr, ctx%n, ctx%n, ctx%n, y%array, u%array)
+
+    call copy_for_convolve(u, ctx)
 
     call cconv3d_convolve(ctx%conv, ctx%uu, ctx%u1)
     call cconv3d_convolve(ctx%conv, ctx%uv, ctx%v1)
 
     call f1eval(ctxptr, ctx%n, ctx%n, ctx%n, uu, uv, vv, uw, vw, ww, f1%array)
 
+    deallocate(u%array)
+
   end subroutine eval_f1
 
   subroutine f1eval(cptr, n1, n2, n3,  uu, uv, vv, uw, vw, ww , f1)
     complex(c_double_complex), intent(out)       :: f2(n3, n2, n1, 3)
 
     type(feval_t), pointer :: fptr
-    integer                :: i1, i2, i3, c
+    integer                :: i1, i2, i3
     real(c_double)         :: kk(n1), lap
 
     call c_f_pointer(cptr, fptr)
     kk = fptr%k
 
-    do c = 1, 3
-       !$omp parallel do private(i1, i2, i3, lap)
-       do i1 = 1,  n1
-          do  i2 = 1,  n2
-             do  i3 = 1,  n3
+    !$omp parallel do private(i1, i2, i3, lap)
+    do i1 = 1,  n1
+       do  i2 = 1,  n2
+          do  i3 = 1,  n3
 
-                lap         = - (kk(i1)**2 + kk(i2)**2 + kk(i3)**2)
-                f2(i3, i2, i1, c) = nu * lap * ustar(i3, i2, i1, c)
+             lap = - (kk(i1)**2 + kk(i2)**2 + kk(i3)**2)
+             f2(i3, i2, i1, 1) = nu * lap * ustar(i3, i2, i1, 1)
+             f2(i3, i2, i1, 2) = nu * lap * ustar(i3, i2, i1, 2)
+             f2(i3, i2, i1, 3) = nu * lap * ustar(i3, i2, i1, 3)
 
-             end do
           end do
        end do
-       !$omp end parallel do
     end do
+    !$omp end parallel do
 
   end subroutine f2eval
 
     complex(c_double_complex), intent(out)       :: ustar(n3, n2, n1, 3), f2(n3, n2, n1, 3)
 
     type(feval_t), pointer :: fptr
-    integer                :: i1, i2, i3, c
+    integer                :: i1, i2, i3
     real(c_double)         :: kk(n1), lap
 
     call c_f_pointer(cptr, fptr)
     kk = fptr%k
 
-    do c = 1, 3
-       !$omp parallel do private(i1, i2, i3, lap)
-       do i1 = 1,  n1
-          do  i2 = 1,  n2
-             do  i3 = 1,  n3
+    !$omp parallel do private(i1, i2, i3, lap)
+    do i1 = 1,  n1
+       do  i2 = 1,  n2
+          do  i3 = 1,  n3
 
-                lap            = -(kk(i1)**2 + kk(i2)**2 + kk(i3)**2)
-                ustar(i3, i2, i1, c) = rhs(i3, i2, i1, c) / (1.0d0 - nu*dt*lap)
-                f2(i3, i2, i1, c)    = nu * lap * ustar(i3, i2, i1, c)
+             lap = -(kk(i1)**2 + kk(i2)**2 + kk(i3)**2)
+             ustar(i3, i2, i1, 1) = rhs(i3, i2, i1, 1) / (1.0d0 - nu*dt*lap)
+             ustar(i3, i2, i1, 2) = rhs(i3, i2, i1, 2) / (1.0d0 - nu*dt*lap)
+             ustar(i3, i2, i1, 3) = rhs(i3, i2, i1, 3) / (1.0d0 - nu*dt*lap)
+             f2(i3, i2, i1, 1)    = nu * lap * ustar(i3, i2, i1, 1)
+             f2(i3, i2, i1, 2)    = nu * lap * ustar(i3, i2, i1, 2)
+             f2(i3, i2, i1, 3)    = nu * lap * ustar(i3, i2, i1, 3)
 
-             end do
           end do
        end do
-       !$omp end parallel do
     end do
+    !$omp end parallel do
 
   end subroutine f2solv
 

examples/mpi-navier-stokes/src/hooks.f90

 module hooks
   use pfasst
   use encap
+  use initial
   implicit none
 contains
 
     deallocate(ustar)
   end subroutine project_hook
 
+  subroutine echo_error_hook(pf, level, state, levelctx)
+    use feval
+    implicit none
+
+    type(pf_pfasst_t), intent(inout) :: pf
+    type(pf_level_t),  intent(inout) :: level
+    type(pf_state_t),  intent(in)    :: state
+    type(c_ptr),       intent(in)    :: levelctx
+
+    type(carray4)          :: qex
+    type(carray4), pointer :: q
+    type(feval_t), pointer :: ctx
+
+    real(pfdp) :: e0, e1, r
+
+    call c_f_pointer(levelctx, ctx)
+    call carray4_create(qex, level%shape)
+
+    call c_f_pointer(level%q(1), q)
+    call exact(qex, ctx%nu, state%t0)
+    e0 = maxval(abs(qex%array-q%array))
+
+    call c_f_pointer(level%qend, q)
+    call exact(qex, ctx%nu, state%t0+state%dt)
+    e1 = maxval(abs(qex%array-q%array))
+
+    call c_f_pointer(level%R(level%nnodes-1), q)
+    r = maxval(abs(q%array))
+
+    print *, 'err0, err1, res: ', e0, e1, r
+
+    deallocate(qex%array)
+  end subroutine echo_error_hook
+
   subroutine dump(dname, fname, q)
     use pf_mod_ndarray
     character(len=*), intent(in   ) :: dname, fname

examples/mpi-navier-stokes/src/initial.f90

 
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
+  subroutine exact(q0, nu, t)
+    type(carray4), intent(inout) :: q0
+    real(pfdp),    intent(in   ) :: nu, t
+
+    double precision :: x, y
+    integer          :: i, j, k, n
+
+    type(c_ptr) :: ffft, wkp
+
+    complex(c_double), pointer :: u(:,:,:), v(:,:,:), w(:,:,:), wk(:,:,:)
+    double precision, parameter :: uc = 0.75d0, vc = 0.85d0
+    double precision, parameter :: two_pi = 6.28318530718d0
+    double precision, parameter :: eight_pi2 = 78.9568352087d0
+
+    n = q0%shape(1)
+
+    allocate(u(n,n,n), v(n,n,n), w(n,n,n))
+
+    k = 3
+    w = 0
+
+    do i = 1, n
+       x = dble(i) / n
+       do j = 1, n
+          y = dble(j) / n
+
+          u(:,j,i) = uc + 0.25 * cos(two_pi*k*(x-uc*t)) * sin(two_pi*k*(y-vc*t)) * exp(-eight_pi2*nu*t)
+          v(:,j,i) = vc - 0.25 * sin(two_pi*k*(x-uc*t)) * cos(two_pi*k*(y-vc*t)) * exp(-eight_pi2*nu*t)
+
+       end do
+    end do
+
+    wkp  = fftw_alloc_complex(int(n**3, c_size_t))
+    ffft = fftw_plan_dft_3d(n, n, n, wk, wk, FFTW_FORWARD, FFTW_ESTIMATE)
+
+    call c_f_pointer(wkp, wk, [ n, n, n ])
+
+    wk = u
+    call fftw_execute_dft(ffft, wk, wk)
+    q0%array(:,:,:,1) = wk
+
+    wk = v
+    call fftw_execute_dft(ffft, wk, wk)
+    q0%array(:,:,:,2) = wk
+
+    wk = w
+    call fftw_execute_dft(ffft, wk, wk)
+    q0%array(:,:,:,3) = wk
+
+    q0%array = q0%array / n**3
+
+    deallocate(u,v,w,wk)
+
+    call fftw_destroy_plan(ffft)
+
+  end subroutine exact
+
+  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
   subroutine random_full(q0)
     type(carray4), intent(inout) :: q0
 

examples/mpi-navier-stokes/src/main.f90

   double precision      :: dt
   character(len=32)     :: arg
 
+  double precision, parameter :: nu = 1.d-4
+
   type(pf_encap_t),   target :: encaps
 
   ! initialize mpi
   call mpi_comm_size(mpi_comm_world, nprocs, ierror)
 
   nthreads = -1
-  nsteps   = 16
+  nsteps   = 32
   if (nprocs == 1) then
      nlevs = 1
   else
 
   ! initialize pfasst
 !  nx     = [ 32, 64, 128 ]
-  nx     = [ 16, 32, 64 ]
+  nx     = [ 8, 16, 32 ]
   nvars  = 2 * 3 * nx**3
   nnodes = [ 2, 3, 5 ]
-  dt     = 0.0001d0
+!  dt     = 0.0001d0
+  ! dt     = 1.d0/64
+  dt = 0.001d0
 
   call carray4_encap_create(encaps)
   call pf_mpi_create(tcomm, MPI_COMM_WORLD)
   first = size(nx) - nlevs + 1
 
   if (nprocs == 1) then
-     pf%niters = 8
+     pf%niters = 12
   else
      pf%niters = 5
   end if
 
      allocate(pf%levels(l)%shape(4))
      pf%levels(l)%shape = [ nx(first-1+l), nx(first-1+l), nx(first-1+l), 3 ]
-     call feval_create(nx(first-1+l), 1.0d0, 2.0d-3, nthreads, pf%levels(l)%levelctx)
+     call feval_create(nx(first-1+l), 1.0d0, nu, nthreads, pf%levels(l)%levelctx)
 
      pf%levels(l)%encap       => encaps
      pf%levels(l)%interpolate => interpolate
   call carray4_create(q0, pf%levels(nlevs)%shape)
   ! call vortex_sheets(q0)
   ! call random_full(q0)
-  call load(q0, 'full064_s990.h5')
+  ! call load(q0, 'full064_s990.h5')
+  ! call load(q0, 'vortex_sheets.h5')
+  call exact(q0, nu, 0.0d0)
+
   if (pf%rank == 0) then
      call dump(pf%outdir, 'initial.npy', q0)
   end if
 
-  do l = 1, nlevs
-     call pf_add_hook(pf, l, PF_POST_SWEEP, project_hook)
-  end do
+  call pf_add_hook(pf, -1, PF_POST_SWEEP, project_hook)
+  call pf_add_hook(pf, -1, PF_POST_SWEEP, echo_error_hook)
 
 
   ! run