diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 8f352757..36d7a86d 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,6 +1,7 @@ set(SRC allocator.f90 backend.f90 + case.f90 common.f90 field.f90 mesh.f90 diff --git a/src/case.f90 b/src/case.f90 new file mode 100644 index 00000000..f0efba67 --- /dev/null +++ b/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 diff --git a/src/solver.f90 b/src/solver.f90 index 4f9547b3..b6820a13 100644 --- a/src/solver.f90 +++ b/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