1: module ex13f90auxmodule
  2: #include <petsc/finclude/petscdm.h>
  3: #include <petsc/finclude/petscdmda.h>
  4:   use petscdm
  5:   implicit none
  6: contains
  7:   !
  8:   ! A subroutine which returns the boundary conditions.
  9:   !
 10:   subroutine get_boundary_cond(b_x, b_y, b_z)
 11:     DMBoundaryType, intent(inout) :: b_x, b_y, b_z

 13:     ! Here you may set the BC types you want
 14:     b_x = DM_BOUNDARY_GHOSTED
 15:     b_y = DM_BOUNDARY_GHOSTED
 16:     b_z = DM_BOUNDARY_GHOSTED

 18:   end subroutine get_boundary_cond
 19:   !
 20:   ! A function which returns the RHS of the equation we are solving
 21:   !
 22:   function dfdt_vdp(t, dt, ib1, ibn, jb1, jbn, kb1, kbn, imax, jmax, kmax, n, f)
 23:     !
 24:     ! Right-hand side for the van der Pol oscillator.  Very simple system of two
 25:     ! ODEs.  See Iserles, eq (5.2).
 26:     !
 27:     PetscReal, intent(in) :: t, dt
 28:     PetscInt, intent(in) :: ib1, ibn, jb1, jbn, kb1, kbn, imax, jmax, kmax, n
 29:     PetscReal, dimension(n, ib1:ibn, jb1:jbn, kb1:kbn), intent(inout) :: f
 30:     PetscReal, dimension(n, imax, jmax, kmax) :: dfdt_vdp
 31:     PetscReal, parameter :: mu = 1.4, one = 1.0
 32:     !
 33:     dfdt_vdp(1, :, :, :) = f(2, 1, 1, 1)
 34:     dfdt_vdp(2, :, :, :) = mu*(one - f(1, 1, 1, 1)**2)*f(2, 1, 1, 1) - f(1, 1, 1, 1)
 35:   end function dfdt_vdp
 36:   !
 37:   ! The standard Forward Euler time-stepping method.
 38:   !
 39:   recursive subroutine forw_euler(t, dt, ib1, ibn, jb1, jbn, kb1, kbn, imax, jmax, kmax, neq, y, dfdt)
 40:     PetscReal, intent(in) :: t, dt
 41:     PetscInt, intent(in) :: ib1, ibn, jb1, jbn, kb1, kbn, imax, jmax, kmax, neq
 42:     PetscReal, dimension(neq, ib1:ibn, jb1:jbn, kb1:kbn), intent(inout) :: y
 43:     !
 44:     ! Define the right-hand side function
 45:     !
 46:     interface
 47:       function dfdt(t, dt, ib1, ibn, jb1, jbn, kb1, kbn, imax, jmax, kmax, n, f)
 48: #include <petsc/finclude/petscsys.h>
 49:         use petscsys
 50:         PetscReal, intent(in) :: t, dt
 51:         PetscInt, intent(in) :: ib1, ibn, jb1, jbn, kb1, kbn, imax, jmax, kmax, n
 52:         PetscReal, dimension(n, ib1:ibn, jb1:jbn, kb1:kbn), intent(inout) :: f
 53:         PetscReal, dimension(n, imax, jmax, kmax) :: dfdt
 54:       end function dfdt
 55:     end interface
 56:     !--------------------------------------------------------------------------
 57:     !
 58:     y(:, 1:imax, 1:jmax, 1:kmax) = y(:, 1:imax, 1:jmax, 1:kmax) + dt*dfdt(t, dt, ib1, ibn, jb1, jbn, kb1, kbn, imax, jmax, kmax, neq, y)
 59:   end subroutine forw_euler
 60:   !
 61:   ! The following 4 subroutines handle the mapping of coordinates. I'll explain
 62:   ! this in detail:
 63:   !    PETSc gives you local arrays which are indexed using the global indices.
 64:   ! This is probably handy in some cases, but when you are re-writing an
 65:   ! existing serial code and want to use DMDAs, you have tons of loops going
 66:   ! from 1 to imax etc. that you don't want to change.
 67:   !    These subroutines re-map the arrays so that all the local arrays go from
 68:   ! 1 to the (local) imax.
 69:   !
 70:   subroutine petsc_to_local(da, vec, array, f, dof, stw)
 71:     use petscdmda
 72:     DM                                                            :: da
 73:     Vec, intent(in)                                                :: vec
 74:     PetscReal, pointer                                            :: array(:, :, :, :)
 75:     PetscInt, intent(in)                                           :: dof, stw
 76:     PetscReal, intent(inout), dimension(:, 1 - stw:, 1 - stw:, 1 - stw:) :: f
 77:     PetscErrorCode                                                :: ierr
 78:     !
 79:     PetscCall(DMDAVecGetArray(da, vec, array, ierr))
 80:     call transform_petsc_us(array, f, stw)
 81:   end subroutine petsc_to_local
 82:   subroutine transform_petsc_us(array, f, stw)
 83:     !Note: this assumed shape-array is what does the "coordinate transformation"
 84:     PetscInt, intent(in)                                   :: stw
 85:     PetscReal, intent(in), dimension(:, 1 - stw:, 1 - stw:, 1 - stw:)  :: array
 86:     PetscReal, intent(inout), dimension(:, 1 - stw:, 1 - stw:, 1 - stw:) :: f
 87:     f(:, :, :, :) = array(:, :, :, :)
 88:   end subroutine transform_petsc_us
 89:   subroutine local_to_petsc(da, vec, array, f, dof, stw)
 90:     use petscdmda
 91:     DM                                                    :: da
 92:     Vec, intent(inout)                                     :: vec
 93:     PetscReal, pointer                                    :: array(:, :, :, :)
 94:     PetscInt, intent(in)                                    :: dof, stw
 95:     PetscReal, intent(inout), dimension(:, 1 - stw:, 1 - stw:, 1 - stw:)  :: f
 96:     PetscErrorCode                                        :: ierr
 97:     call transform_us_petsc(array, f, stw)
 98:     PetscCall(DMDAVecRestoreArray(da, vec, array, ierr))
 99:   end subroutine local_to_petsc
100:   subroutine transform_us_petsc(array, f, stw)
101:     !Note: this assumed shape-array is what does the "coordinate transformation"
102:     PetscInt, intent(in)                                     :: stw
103:     PetscReal, intent(inout), dimension(:, 1 - stw:, 1 - stw:, 1 - stw:) :: array
104:     PetscReal, intent(in), dimension(:, 1 - stw:, 1 - stw:, 1 - stw:)      :: f
105:     array(:, :, :, :) = f(:, :, :, :)
106:   end subroutine transform_us_petsc
107: end module ex13f90auxmodule