1module ex13f90auxmodule 2#include <petsc/finclude/petscdm.h> 3#include <petsc/finclude/petscdmda.h> 4 use petscdm 5 implicit none 6contains 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 12 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 17 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 107end module ex13f90auxmodule 108