Restart simulation fixes

src/case.f90

@@ -206,6 +206,17 @@ subroutine case_init_common(C)
     call json_get(C%params, 'case.numerics.polynomial_order', lx)
     lx = lx + 1 ! add 1 to get poly order
     call C%fluid%init(C%msh, lx, C%params, C%usr, C%material_properties)
+    C%fluid%chkp%tlag => C%tlag
+    C%fluid%chkp%dtlag => C%dtlag
+    select type(f => C%fluid)
+    type is(fluid_pnpn_t)
+      f%chkp%abx1 => f%abx1
+      f%chkp%abx2 => f%abx2
+      f%chkp%aby1 => f%aby1
+      f%chkp%aby2 => f%aby2
+      f%chkp%abz1 => f%abz1
+      f%chkp%abz2 => f%abz2
+    end select
 
 
     !
@@ -232,6 +243,9 @@ subroutine case_init_common(C)
        call C%scalar%init(C%msh, C%fluid%c_Xh, C%fluid%gs_Xh, C%params, C%usr,&
                           C%material_properties)
        call C%fluid%chkp%add_scalar(C%scalar%s)
+       C%fluid%chkp%abs1 => C%scalar%abx1
+       C%fluid%chkp%abs2 => C%scalar%abx2
+       C%fluid%chkp%slag => C%scalar%slag
     end if
 
     !
@@ -381,11 +395,11 @@ subroutine case_init_common(C)
     ! Save checkpoints (if nothing specified, default to saving at end of sim)
     !
     call json_get_or_default(C%params, 'case.output_checkpoints',&
-                             logical_val, .false.)
+                             logical_val, .true.)
     if (logical_val) then
        C%f_chkp = chkp_output_t(C%fluid%chkp, path=output_directory)
-       call json_get(C%params, 'case.checkpoint_control', string_val)
-       call json_get(C%params, 'case.checkpoint_value', real_val)
+       call json_get_or_default(C%params, 'case.checkpoint_control', string_val,"simulationtime")
+       call json_get_or_default(C%params, 'case.checkpoint_value', real_val,1e10_rp)
        call C%s%add(C%f_chkp, real_val, string_val)
     end if
 

src/common/checkpoint.f90

@@ -50,14 +50,30 @@ module checkpoint
      type(field_t), pointer :: w => null()
      type(field_t), pointer :: p => null()
 
+
      !
      ! Optional payload
      !
      type(field_series_t), pointer :: ulag => null()
      type(field_series_t), pointer :: vlag => null()
      type(field_series_t), pointer :: wlag => null()
+
+     real(kind=rp), pointer :: tlag(:) => null()
+     real(kind=rp), pointer :: dtlag(:) => null()
+
+     !> for pnpn
+     type(field_t), pointer :: abx1 => null()
+     type(field_t), pointer :: abx2 => null()
+     type(field_t), pointer :: aby1 => null()
+     type(field_t), pointer :: aby2 => null()
+     type(field_t), pointer :: abz1 => null()
+     type(field_t), pointer :: abz2 => null()
 
      type(field_t), pointer :: s => null()
+     type(field_series_t), pointer :: slag => null()
+
+     type(field_t), pointer :: abs1 => null()
+     type(field_t), pointer :: abs2 => null()
 
      real(kind=dp) :: t         !< Restart time (valid after load)
      type(mesh_t) :: previous_mesh
@@ -125,13 +141,15 @@ subroutine chkp_sync_host(this)
 
     if (NEKO_BCKND_DEVICE .eq. 1) then
        associate(u=>this%u, v=>this%v, w=>this%w, &
-            ulag=>this%ulag, vlag=>this%vlag, wlag=>this%wlag)
+            ulag=>this%ulag, vlag=>this%vlag, wlag=>this%wlag, &
+            p=>this%p)
 
          if (associated(this%u) .and. associated(this%v) .and. &
-              associated(this%w)) then
+              associated(this%w) .and. associated(this%p)) then
             call device_memcpy(u%x, u%x_d, u%dof%size(), DEVICE_TO_HOST)
             call device_memcpy(v%x, v%x_d, v%dof%size(), DEVICE_TO_HOST)
             call device_memcpy(w%x, w%x_d, w%dof%size(), DEVICE_TO_HOST)
+            call device_memcpy(p%x, p%x_d, p%dof%size(), DEVICE_TO_HOST)
          end if
 
          if (associated(this%ulag) .and. associated(this%vlag) .and. &
@@ -150,12 +168,33 @@ subroutine chkp_sync_host(this)
                                w%dof%size(), DEVICE_TO_HOST)
             call device_memcpy(wlag%lf(2)%x, wlag%lf(2)%x_d, &
                                w%dof%size(), DEVICE_TO_HOST)
+            call device_memcpy(this%abx1%x, this%abx1%x_d, &
+                               w%dof%size(), DEVICE_TO_HOST)
+            call device_memcpy(this%abx2%x, this%abx2%x_d, &
+                               w%dof%size(), DEVICE_TO_HOST)
+            call device_memcpy(this%aby1%x, this%aby1%x_d, &
+                               w%dof%size(), DEVICE_TO_HOST)
+            call device_memcpy(this%aby2%x, this%aby2%x_d, &
+                               w%dof%size(), DEVICE_TO_HOST)
+            call device_memcpy(this%abz1%x, this%abz1%x_d, &
+                               w%dof%size(), DEVICE_TO_HOST)
+            call device_memcpy(this%abz2%x, this%abz2%x_d, &
+                               w%dof%size(), DEVICE_TO_HOST)
          end if
          if (associated(this%s)) then
             call device_memcpy(this%s%x, this%s%x_d, &
                                this%s%dof%size(), DEVICE_TO_HOST)
+            call device_memcpy(this%slag%lf(1)%x, this%slag%lf(1)%x_d, &
+                               this%s%dof%size(), DEVICE_TO_HOST)
+            call device_memcpy(this%slag%lf(2)%x, this%slag%lf(2)%x_d, &
+                               this%s%dof%size(), DEVICE_TO_HOST)
+            call device_memcpy(this%abs1%x, this%abs1%x_d, &
+                               w%dof%size(), DEVICE_TO_HOST)
+            call device_memcpy(this%abs2%x, this%abs2%x_d, &
+                               w%dof%size(), DEVICE_TO_HOST)
          end if
        end associate
+    call device_sync()
     end if
 
   end subroutine chkp_sync_host
@@ -166,13 +205,15 @@ subroutine chkp_sync_device(this)
 
     if (NEKO_BCKND_DEVICE .eq. 1) then
        associate(u=>this%u, v=>this%v, w=>this%w, &
-            ulag=>this%ulag, vlag=>this%vlag, wlag=>this%wlag)
+            ulag=>this%ulag, vlag=>this%vlag, wlag=>this%wlag,&
+            p=>this%p)
 
          if (associated(this%u) .and. associated(this%v) .and. &
               associated(this%w)) then
             call device_memcpy(u%x, u%x_d, u%dof%size(), HOST_TO_DEVICE)
             call device_memcpy(v%x, v%x_d, v%dof%size(), HOST_TO_DEVICE)
             call device_memcpy(w%x, w%x_d, w%dof%size(), HOST_TO_DEVICE)
+            call device_memcpy(p%x, p%x_d, p%dof%size(), HOST_TO_DEVICE)
          end if
 
          if (associated(this%ulag) .and. associated(this%vlag) .and. &

src/fluid/bcknd/device/pnpn_res_device.F90

@@ -255,6 +255,7 @@ subroutine pnpn_prs_res_device_compute(p, p_res, u, v, w, u_e, v_e, w_e, &
     call curl(ta1, ta2, ta3, u_e, v_e, w_e, work1, work2, c_Xh)
     call curl(wa1, wa2, wa3, ta1, ta2, ta3, work1, work2, c_Xh)
 
+
 #ifdef HAVE_HIP
     call pnpn_prs_res_part1_hip(ta1%x_d, ta2%x_d, ta3%x_d, &
          wa1%x_d, wa2%x_d, wa3%x_d, f_x%x_d, f_y%x_d, f_z%x_d, &
@@ -270,6 +271,7 @@ subroutine pnpn_prs_res_device_compute(p, p_res, u, v, w, u_e, v_e, w_e, &
          c_Xh%B_d, c_Xh%h1_d, mu, rho, n) 
 #endif
      c_Xh%ifh2 = .false.
+    call device_cfill(c_Xh%h1_d,1.0_rp / rho,n)
 
     call gs_Xh%op(ta1, GS_OP_ADD) 
     call gs_Xh%op(ta2, GS_OP_ADD) 

src/fluid/fluid_pnpn.f90

@@ -44,6 +44,7 @@ module fluid_pnpn
   use fluid_aux    
   use time_scheme_controller
   use projection
+  use device
   use logger
   use advection
   use profiler
@@ -107,6 +108,7 @@ module fluid_pnpn
      procedure, pass(this) :: init => fluid_pnpn_init
      procedure, pass(this) :: free => fluid_pnpn_free
      procedure, pass(this) :: step => fluid_pnpn_step
+     procedure, pass(this) :: restart => fluid_pnpn_restart
   end type fluid_pnpn_t
 
 contains
@@ -158,10 +160,15 @@ subroutine fluid_pnpn_init(this, msh, lx, params, user, material_properties)
       call this%abx1%init(dm_Xh, "abx1")
       call this%aby1%init(dm_Xh, "aby1")
       call this%abz1%init(dm_Xh, "abz1")
-
       call this%abx2%init(dm_Xh, "abx2")
       call this%aby2%init(dm_Xh, "aby2")
       call this%abz2%init(dm_Xh, "abz2")
+      this%abx1 = 0.0_rp
+      this%aby1 = 0.0_rp
+      this%abz1 = 0.0_rp
+      this%abx2 = 0.0_rp
+      this%aby2 = 0.0_rp
+      this%abz2 = 0.0_rp
 
       call this%du%init(dm_Xh, 'du')
       call this%dv%init(dm_Xh, 'dv')
@@ -271,6 +278,139 @@ subroutine fluid_pnpn_init(this, msh, lx, params, user, material_properties)
 
   end subroutine fluid_pnpn_init
 
+  subroutine fluid_pnpn_restart(this,dtlag, tlag)
+    class(fluid_pnpn_t), target, intent(inout) :: this
+    real(kind=rp) :: dtlag(10), tlag(10)
+    type(field_t) :: u_temp, v_temp, w_temp
+    integer :: i, n
+
+    n = this%u%dof%size()
+    ! Make sure that continuity is maintained (important for interpolation) 
+    call col2(this%u%x,this%c_Xh%mult,this%u%dof%size())
+    call col2(this%v%x,this%c_Xh%mult,this%u%dof%size())
+    call col2(this%w%x,this%c_Xh%mult,this%u%dof%size())
+    call col2(this%p%x,this%c_Xh%mult,this%u%dof%size())
+    do i = 1, this%ulag%size()
+       call col2(this%ulag%lf(i)%x,this%c_Xh%mult,this%u%dof%size())
+       call col2(this%vlag%lf(i)%x,this%c_Xh%mult,this%u%dof%size())
+       call col2(this%wlag%lf(i)%x,this%c_Xh%mult,this%u%dof%size())
+    end do
+
+    call col2(this%abx1%x,this%c_Xh%mult,this%abx1%dof%size())
+    call col2(this%abx2%x,this%c_Xh%mult,this%abx2%dof%size())
+    call col2(this%aby1%x,this%c_Xh%mult,this%aby1%dof%size())
+    call col2(this%aby2%x,this%c_Xh%mult,this%aby2%dof%size())
+    call col2(this%abz1%x,this%c_Xh%mult,this%abz1%dof%size())
+    call col2(this%abz2%x,this%c_Xh%mult,this%abz2%dof%size())
+
+
+    if (NEKO_BCKND_DEVICE .eq. 1) then
+       associate(u=>this%u, v=>this%v, w=>this%w, &
+            ulag=>this%ulag, vlag=>this%vlag, wlag=>this%wlag,&
+            p=>this%p)
+       call device_memcpy(u%x, u%x_d, u%dof%size(), HOST_TO_DEVICE)
+       call device_memcpy(v%x, v%x_d, v%dof%size(), HOST_TO_DEVICE)
+       call device_memcpy(w%x, w%x_d, w%dof%size(), HOST_TO_DEVICE)
+       call device_memcpy(p%x, p%x_d, p%dof%size(), HOST_TO_DEVICE)
+       call device_memcpy(ulag%lf(1)%x, ulag%lf(1)%x_d, &
+                          u%dof%size(), HOST_TO_DEVICE)
+       call device_memcpy(ulag%lf(2)%x, ulag%lf(2)%x_d, &
+                          u%dof%size(), HOST_TO_DEVICE)
+
+       call device_memcpy(vlag%lf(1)%x, vlag%lf(1)%x_d, &
+                          v%dof%size(), HOST_TO_DEVICE)
+       call device_memcpy(vlag%lf(2)%x, vlag%lf(2)%x_d, &
+                          v%dof%size(), HOST_TO_DEVICE)
+
+       call device_memcpy(wlag%lf(1)%x, wlag%lf(1)%x_d, &
+                          w%dof%size(), HOST_TO_DEVICE)
+       call device_memcpy(wlag%lf(2)%x, wlag%lf(2)%x_d, &
+                          w%dof%size(), HOST_TO_DEVICE)
+       call device_memcpy(this%abx1%x, this%abx1%x_d, &
+                          w%dof%size(), HOST_TO_DEVICE)
+       call device_memcpy(this%abx2%x, this%abx2%x_d, &
+                          w%dof%size(), HOST_TO_DEVICE)
+       call device_memcpy(this%aby1%x, this%aby1%x_d, &
+                          w%dof%size(), HOST_TO_DEVICE)
+       call device_memcpy(this%aby2%x, this%aby2%x_d, &
+                          w%dof%size(), HOST_TO_DEVICE)
+       call device_memcpy(this%abz1%x, this%abz1%x_d, &
+                          w%dof%size(), HOST_TO_DEVICE)
+       call device_memcpy(this%abz2%x, this%abz2%x_d, &
+                          w%dof%size(), HOST_TO_DEVICE)
+       end associate
+   end if
+
+
+
+    call this%gs_Xh%op(this%u,GS_OP_ADD)
+    call this%gs_Xh%op(this%v,GS_OP_ADD)
+    call this%gs_Xh%op(this%w,GS_OP_ADD)
+    call this%gs_Xh%op(this%p,GS_OP_ADD)
+
+    do i = 1, this%ulag%size()
+       call this%gs_Xh%op(this%ulag%lf(i),GS_OP_ADD)
+       call this%gs_Xh%op(this%vlag%lf(i),GS_OP_ADD)
+       call this%gs_Xh%op(this%wlag%lf(i),GS_OP_ADD)
+    end do
+    call this%gs_Xh%op(this%abx1,GS_OP_ADD)
+    call this%gs_Xh%op(this%abx2,GS_OP_ADD)
+    call this%gs_Xh%op(this%aby1,GS_OP_ADD)
+    call this%gs_Xh%op(this%aby2,GS_OP_ADD)
+    call this%gs_Xh%op(this%abz1,GS_OP_ADD)
+    call this%gs_Xh%op(this%abz2,GS_OP_ADD)
+
+    !! If we would decide to only restart from lagged fields instead of asving abx1, aby1 etc.
+    !! Observe that one also needs to recompute the focing at the old time steps
+    !u_temp = this%ulag%lf(2)
+    !v_temp = this%vlag%lf(2)
+    !w_temp = this%wlag%lf(2)
+    !! Compute the source terms
+    !call this%source_term%compute(tlag(2), -1)
+    !  
+    !! Pre-multiply the source terms with the mass matrix.
+    !if (NEKO_BCKND_DEVICE .eq. 1) then
+    !   call device_opcolv(this%f_x%x_d, this%f_y%x_d, this%f_z%x_d, this%c_Xh%B_d, this%msh%gdim, n)
+    !else
+    !   call opcolv(this%f_x%x, this%f_y%x, this%f_z%x, this%c_Xh%B, this%msh%gdim, n)
+    !end if
+
+    !! Add the advection operators to the right-hand-side.
+    !call this%adv%compute(u_temp, v_temp, w_temp, &
+    !                      this%f_x%x, this%f_y%x, this%f_z%x, &
+    !                      this%Xh, this%c_Xh, this%dm_Xh%size())
+    !this%abx2 = this%f_x
+    !this%aby2 = this%f_y
+    !this%abz2 = this%f_z
+    !
+    !u_temp = this%ulag%lf(1)
+    !v_temp = this%vlag%lf(1)
+    !w_temp = this%wlag%lf(1)
+    !call this%source_term%compute(tlag(1), 0)
+
+    !! Pre-multiply the source terms with the mass matrix.
+    !if (NEKO_BCKND_DEVICE .eq. 1) then
+    !   call device_opcolv(this%f_x%x_d, this%f_y%x_d, this%f_z%x_d, this%c_Xh%B_d, this%msh%gdim, n)
+    !else
+    !   call opcolv(this%f_x%x, this%f_y%x, this%f_z%x, this%c_Xh%B, this%msh%gdim, n)
+    !end if
+
+    !! Pre-multiply the source terms with the mass matrix.
+    !if (NEKO_BCKND_DEVICE .eq. 1) then
+    !   call device_opcolv(this%f_x%x_d, this%f_y%x_d, this%f_z%x_d, this%c_Xh%B_d, this%msh%gdim, n)
+    !else
+    !   call opcolv(this%f_x%x, this%f_y%x, this%f_z%x, this%c_Xh%B, this%msh%gdim, n)
+    !end if
+
+    !call this%adv%compute(u_temp, v_temp, w_temp, &
+    !                      this%f_x%x, this%f_y%x, this%f_z%x, &
+    !                      this%Xh, this%c_Xh, this%dm_Xh%size())
+    !this%abx1 = this%f_x
+    !this%aby1 = this%f_y
+    !this%abz1 = this%f_z
+
+  end subroutine fluid_pnpn_restart
+
   subroutine fluid_pnpn_free(this)
     class(fluid_pnpn_t), intent(inout) :: this
 
@@ -384,7 +524,6 @@ subroutine fluid_pnpn_step(this, t, tstep, dt, ext_bdf)
       call this%scratch%request_field(u_e, temp_indices(1))
       call this%scratch%request_field(v_e, temp_indices(2))
       call this%scratch%request_field(w_e, temp_indices(3))
-
       call sumab%compute_fluid(u_e, v_e, w_e, u, v, w, &
            ulag, vlag, wlag, ext_bdf%advection_coeffs, ext_bdf%nadv)
 

src/fluid/fluid_scheme.f90

@@ -137,6 +137,7 @@ module fluid_scheme
      procedure(fluid_scheme_init_intrf), pass(this), deferred :: init
      procedure(fluid_scheme_free_intrf), pass(this), deferred :: free
      procedure(fluid_scheme_step_intrf), pass(this), deferred :: step
+     procedure(fluid_scheme_restart_intrf), pass(this), deferred :: restart
      generic :: scheme_init => fluid_scheme_init_all, fluid_scheme_init_uvw
   end type fluid_scheme_t
 
@@ -180,6 +181,18 @@ subroutine fluid_scheme_step_intrf(this, t, tstep, dt, ext_bdf)
      end subroutine fluid_scheme_step_intrf
   end interface
 
+  !> Abstract interface to restart a fluid scheme
+  abstract interface
+     subroutine fluid_scheme_restart_intrf(this, dtlag, tlag)
+       import fluid_scheme_t
+       import rp
+       class(fluid_scheme_t), target, intent(inout) :: this
+       real(kind=rp) :: dtlag(10), tlag(10)
+
+     end subroutine fluid_scheme_restart_intrf
+  end interface
+
+
 contains
 
   !> Initialize common data for the current scheme