1. Matthew Emmett
  2. libpfasst

Commits

Matthew Emmett  committed d945825

NS: Add echo_energy_hook, which shows difference of fine energy.

  • Participants
  • Parent commits 1b2060e
  • Branches master

Comments (0)

Files changed (5)

File examples/mpi-navier-stokes/Makefile

View file
 #
 
 build/feval.o:     build/pfasst.o build/encap.o build/fftwpp.o
-build/hooks.o:     build/numpy.o build/initial.o
+build/hooks.o:     build/numpy.o build/initial.o build/transfer.o
 build/transfer.o:  build/feval.o
 build/main.o:      build/feval.o build/hooks.o build/transfer.o build/initial.o build/probin.o

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

View file
 
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
+  subroutine energy(q, e)
+    type(carray4), intent(in   ) :: q
+    real(pfdp),    intent(  out) :: e
+
+    integer       :: n, c, i1, i2, i3
+    complex(pfdp) :: z
+
+    n = q%shape(1)
+
+    e = 0
+
+    do c = 1, 3
+       do i1 = 1,  n
+          do  i2 = 1,  n
+             do  i3 = 1,  n
+                z = q%array(i3, i2, i1, c)
+                e = e + real(z * conjg(z))
+             end do
+          end do
+       end do
+    end do
+
+    e = 0.5 * e
+  end subroutine energy
+
+  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
   subroutine project(cptr, n1, n2, n3, ustar, u)
     ! Project ustar to divergence free u
     !

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

View file
     deallocate(qex%array)
   end subroutine echo_error_hook
 
+  subroutine echo_energy_hook(pf, level, state, levelctx)
+    use feval
+    use transfer
+    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), pointer :: q
+    type(feval_t), pointer :: ctx
+
+    type(carray4) :: g
+    integer :: gshape(4)
+
+    real(pfdp) :: e0, e1
+    real(pfdp), save :: energy0 = 0
+
+    if (level%level /= pf%nlevels) return
+
+    call c_f_pointer(levelctx, ctx)
+    call c_f_pointer(level%qend, q)
+    call energy(q, e1)
+
+    gshape(1:3) = q%shape(1:3) / 2
+    gshape(4) = 3
+    call carray4_create(g, gshape)
+    call restrictq(q, g)
+    call energy(g, e0)
+
+    print *, 'step, iter, level, energy: ', state%step, state%iter, level%level, e1, e0, abs(e1-e0), abs(e1-e0-energy0)
+
+    energy0 = e1-e0
+
+    deallocate(g%array)
+  end subroutine echo_energy_hook
+
   subroutine dump(dname, fname, q)
     use pf_mod_ndarray
     character(len=*), intent(in   ) :: dname, fname

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

View file
 
   call pf_add_hook(pf, -1, PF_POST_SWEEP, project_hook)
   ! call pf_add_hook(pf, -1, PF_POST_SWEEP, echo_error_hook)
+  call pf_add_hook(pf, -1, PF_POST_SWEEP, echo_energy_hook)
 
   ! run
   if (pf%rank == 0) then

File examples/mpi-navier-stokes/src/transfer.f90

View file
     !$omp end parallel do
   end subroutine interpolate
 
+  subroutine restrictq(qF, qG)
+    type(carray4), intent(in   ) :: qF
+    type(carray4), intent(inout) :: qG
+
+    integer :: i, j, k, ii, jj, kk, c, nG, nF
+
+    nG = qG%shape(1)
+    nF = qF%shape(1)
+
+    if (nF == nG) then
+       qG%array = qF%array
+       return
+    end if
+
+    !$omp parallel workshare
+    qG%array = 0.0d0
+    !$omp end parallel workshare
+
+    !$omp parallel do private(i,j,k,ii,jj,kk,c)
+    do k = 1, nG
+       if (k <= nG/2) then
+          kk = k
+       else if (k > nG/2+1) then
+          kk = nF - nG + k
+       else
+          cycle
+       end if
+
+       do j = 1, nG
+          if (j <= nG/2) then
+             jj = j
+          else if (j > nG/2+1) then
+             jj = nF - nG + j
+          else
+             cycle
+          end if
+
+          do i = 1, nG
+             if (i <= nG/2) then
+                ii = i
+             else if (i > nG/2+1) then
+                ii = nF - nG + i
+             else
+                cycle
+             end if
+
+             do c = 1, 3
+                qG%array(i, j, k, c) = qF%array(ii, jj, kk, c)
+             end do
+          end do
+       end do
+    end do
+    !$omp end parallel do
+  end subroutine restrictq
+
   subroutine restrict(qFptr, qGptr, levelF, ctxF, levelG, ctxG, t)
     type(c_ptr), intent(in   ), value :: qFptr, qGptr, ctxF, ctxG
     integer,     intent(in   )        :: levelF, levelG