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