Actual source code: ex13f90aux.F90

      1: #include <petsc/finclude/petscdmda.h>
      2: module ex13f90auxmodule
      3:   use petscdm
      4:   use petscdmda
      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_PETSC_REAL_KIND
     32:     !
     33:     dfdt_vdp(1, :, :, :) = f(2, 1, 1, 1)
     34:     dfdt_vdp(2, :, :, :) = mu*(1.0_PETSC_REAL_KIND - 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:         use petscsys
     49:         PetscReal, intent(in) :: t, dt
     50:         PetscInt, intent(in) :: ib1, ibn, jb1, jbn, kb1, kbn, imax, jmax, kmax, n
     51:         PetscReal, dimension(n, ib1:ibn, jb1:jbn, kb1:kbn), intent(inout) :: f
     52:         PetscReal, dimension(n, imax, jmax, kmax) :: dfdt
     53:       end function dfdt
     54:     end interface
     55:     !--------------------------------------------------------------------------
     56:     !
     57:     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)
     58:   end subroutine forw_euler
     59:   !
     60:   ! The following 4 subroutines handle the mapping of coordinates. I'll explain
     61:   ! this in detail:
     62:   !    PETSc gives you local arrays which are indexed using the global indices.
     63:   ! This is probably handy in some cases, but when you are re-writing an
     64:   ! existing serial code and want to use DMDAs, you have tons of loops going
     65:   ! from 1 to imax etc. that you don't want to change.
     66:   !    These subroutines re-map the arrays so that all the local arrays go from
     67:   ! 1 to the (local) imax.
     68:   !
     69:   subroutine petsc_to_local(da, vec, array, f, dof, stw)
     70:     DM :: da
     71:     Vec, intent(in)     :: vec
     72:     PetscReal, pointer  :: array(:, :, :, :)
     73:     PetscInt, intent(in):: dof, stw
     74:     PetscReal, intent(inout), dimension(:, 1 - stw:, 1 - stw:, 1 - stw:) :: f
     75:     PetscErrorCode                                                :: ierr
     76:     !
     77:     PetscCall(DMDAVecGetArray(da, vec, array, ierr))
     78:     call transform_petsc_us(array, f, stw)
     79:   end subroutine petsc_to_local
     80:   subroutine transform_petsc_us(array, f, stw)
     81:     !Note: this assumed shape-array is what does the "coordinate transformation"
     82:     PetscInt, intent(in)                                   :: stw
     83:     PetscReal, intent(in), dimension(:, 1 - stw:, 1 - stw:, 1 - stw:)  :: array
     84:     PetscReal, intent(inout), dimension(:, 1 - stw:, 1 - stw:, 1 - stw:) :: f
     85:     f(:, :, :, :) = array(:, :, :, :)
     86:   end subroutine transform_petsc_us
     87:   subroutine local_to_petsc(da, vec, array, f, dof, stw)
     88:     DM :: da
     89:     Vec, intent(inout)   :: vec
     90:     PetscReal, pointer   :: array(:, :, :, :)
     91:     PetscInt, intent(in) :: dof, stw
     92:     PetscReal, intent(inout), dimension(:, 1 - stw:, 1 - stw:, 1 - stw:)  :: f
     93:     PetscErrorCode :: ierr
     94:     call transform_us_petsc(array, f, stw)
     95:     PetscCall(DMDAVecRestoreArray(da, vec, array, ierr))
     96:   end subroutine local_to_petsc
     97:   subroutine transform_us_petsc(array, f, stw)
     98:     !Note: this assumed shape-array is what does the "coordinate transformation"
     99:     PetscInt, intent(in) :: stw
    100:     PetscReal, intent(inout), dimension(:, 1 - stw:, 1 - stw:, 1 - stw:) :: array
    101:     PetscReal, intent(in), dimension(:, 1 - stw:, 1 - stw:, 1 - stw:)      :: f
    102:     array(:, :, :, :) = f(:, :, :, :)
    103:   end subroutine transform_us_petsc
    104: end module ex13f90auxmodule