Implement case class.

src/CMakeLists.txt

@@ -1,6 +1,7 @@
 set(SRC
   allocator.f90
   backend.f90
+  case.f90
   common.f90
   field.f90
   mesh.f90

src/case.f90

@@ -0,0 +1,119 @@
+module m_case
+  use m_allocator, only: allocator_t, field_t
+  use m_common, only: dp, DIR_X
+  use m_solver, only: solver_t
+  use m_time_integrator, only: time_intg_t
+
+  implicit none
+
+  type :: case_t
+    !! Base case type
+    !!
+    !! Each specific case extends this base class and overwrites its procedures
+    class(base_backend_t), pointer :: backend
+    class(solver_t), pointer :: solver
+    !class(time_intg_t) :: time_integrator
+    type(allocator_t) :: host_allocator
+  contains
+    procedure :: set_ICs
+    procedure :: impose_BCs
+    procedure :: post_transeq
+    procedure :: postprocess
+    procedure :: run
+  end type case_t
+
+  interface case_t
+    module procedure init
+  end interface case_t 
+
+contains
+
+  function init(solver) result(case)
+    implicit none
+
+    class(solver_t), target, intent(inout) :: solver
+    type(case_t) :: case
+
+    case%solver => solver
+
+    call case%set_ICs()
+
+  end function init
+
+  subroutine set_ICs(self)
+    implicit none
+
+    class(case_t) :: self
+
+  end subroutine set_ICs
+
+  subroutine impose_BCs(self)
+    implicit none
+
+    class(case_t) :: self
+
+  end subroutine impose_BCs
+
+  subroutine post_transeq(self, du, dv, dw)
+    !! Use this for forcings; e.g. constant pressure gradient
+    implicit none
+
+    class(case_t) :: self
+    class(field_t), intent(inout) :: du, dv, dw
+
+  end subroutine post_transeq
+
+  subroutine postprocess(self, t)
+    implicit none
+
+    class(case_t) :: self
+    real(dp), intent(in) :: t
+
+    ! for example call enstrophy from vector_calculus
+  end subroutine postprocess
+
+  subroutine run(self)
+    implicit none
+
+    class(case_t) :: self
+
+    class(field_t), pointer :: du, dv, dw
+
+    real(dp) :: t
+    integer :: i, j
+
+    do i = 1, self%solver%n_iters
+      do j = 1, self%solver%time_integrator%nstage
+        du => self%solver%backend%allocator%get_block(DIR_X)
+        dv => self%solver%backend%allocator%get_block(DIR_X)
+        dw => self%solver%backend%allocator%get_block(DIR_X)
+
+        call self%solver%transeq(du, dv, dw, self%solver%u, self%solver%v, self%solver%w)
+
+        call self%post_transeq(du, dv, dw)
+
+        ! time integration
+        call self%solver%time_integrator%step( &
+          self%solver%u, self%solver%v, self%solver%w, du, dv, dw, self%solver%dt &
+          )
+
+        call self%solver%backend%allocator%release_block(du)
+        call self%solver%backend%allocator%release_block(dv)
+        call self%solver%backend%allocator%release_block(dw)
+
+        ! impose boundary conditions
+        call self%impose_BCs()
+
+        ! pressure
+        call self%solver%pressure_correction()
+      end do
+
+      if (mod(i, self%solver%n_output) == 0) then
+        t = i*self%solver%dt
+        call self%postprocess(t)
+      end if
+    end do
+
+  end subroutine run
+
+end module m_case

src/solver.f90

@@ -58,6 +58,7 @@ module m_solver
     procedure(poisson_solver), pointer :: poisson => null()
   contains
     procedure :: transeq
+    procedure :: pressure_correction
     procedure :: divergence_v2p
     procedure :: gradient_p2v
     procedure :: curl
@@ -406,6 +407,42 @@ subroutine poisson_cg(self, pressure, div_u)
 
   end subroutine poisson_cg
 
+  subroutine pressure_correction(self)
+    implicit none
+
+    class(solver_t) :: self
+
+    class(field_t), pointer :: div_u, pressure, dpdx, dpdy, dpdz
+
+    div_u => self%backend%allocator%get_block(DIR_Z)
+
+    call self%divergence_v2p(div_u, self%u, self%v, self%w)
+
+    pressure => self%backend%allocator%get_block(DIR_Z, CELL)
+
+    call self%poisson(pressure, div_u)
+
+    call self%backend%allocator%release_block(div_u)
+
+    dpdx => self%backend%allocator%get_block(DIR_X)
+    dpdy => self%backend%allocator%get_block(DIR_X)
+    dpdz => self%backend%allocator%get_block(DIR_X)
+
+    call self%gradient_p2v(dpdx, dpdy, dpdz, pressure)
+
+    call self%backend%allocator%release_block(pressure)
+
+    ! velocity correction
+    call self%backend%vecadd(-1._dp, dpdx, 1._dp, self%u)
+    call self%backend%vecadd(-1._dp, dpdy, 1._dp, self%v)
+    call self%backend%vecadd(-1._dp, dpdz, 1._dp, self%w)
+
+    call self%backend%allocator%release_block(dpdx)
+    call self%backend%allocator%release_block(dpdy)
+    call self%backend%allocator%release_block(dpdz)
+
+  end subroutine pressure_correction
+
   subroutine output(self, t)
     implicit none
 
@@ -489,33 +526,7 @@ subroutine run(self)
         call self%backend%allocator%release_block(dv)
         call self%backend%allocator%release_block(dw)
 
-        ! pressure
-        div_u => self%backend%allocator%get_block(DIR_Z)
-
-        call self%divergence_v2p(div_u, self%u, self%v, self%w)
-
-        pressure => self%backend%allocator%get_block(DIR_Z, CELL)
-
-        call self%poisson(pressure, div_u)
-
-        call self%backend%allocator%release_block(div_u)
-
-        dpdx => self%backend%allocator%get_block(DIR_X)
-        dpdy => self%backend%allocator%get_block(DIR_X)
-        dpdz => self%backend%allocator%get_block(DIR_X)
-
-        call self%gradient_p2v(dpdx, dpdy, dpdz, pressure)
-
-        call self%backend%allocator%release_block(pressure)
-
-        ! velocity correction
-        call self%backend%vecadd(-1._dp, dpdx, 1._dp, self%u)
-        call self%backend%vecadd(-1._dp, dpdy, 1._dp, self%v)
-        call self%backend%vecadd(-1._dp, dpdz, 1._dp, self%w)
-
-        call self%backend%allocator%release_block(dpdx)
-        call self%backend%allocator%release_block(dpdy)
-        call self%backend%allocator%release_block(dpdz)
+        call self%pressure_correction()
       end do
 
       if (mod(i, self%n_output) == 0) then