1#include <petsc/finclude/petscdm.h> 2#include <petsc/finclude/petscdmda.h> 3#include <petsc/finclude/petscsys.h> 4module ex13f90auxmodule 5 use petscdm 6 implicit none 7contains 8 ! 9 ! A subroutine which returns the boundary conditions. 10 ! 11 subroutine get_boundary_cond(b_x, b_y, b_z) 12 DMBoundaryType, intent(inout) :: b_x, b_y, b_z 13 14 ! Here you may set the BC types you want 15 b_x = DM_BOUNDARY_GHOSTED 16 b_y = DM_BOUNDARY_GHOSTED 17 b_z = DM_BOUNDARY_GHOSTED 18 19 end subroutine get_boundary_cond 20 ! 21 ! A function which returns the RHS of the equation we are solving 22 ! 23 function dfdt_vdp(t, dt, ib1, ibn, jb1, jbn, kb1, kbn, imax, jmax, kmax, n, f) 24 ! 25 ! Right-hand side for the van der Pol oscillator. Very simple system of two 26 ! ODEs. See Iserles, eq (5.2). 27 ! 28 PetscReal, intent(in) :: t, dt 29 PetscInt, intent(in) :: ib1, ibn, jb1, jbn, kb1, kbn, imax, jmax, kmax, n 30 PetscReal, dimension(n, ib1:ibn, jb1:jbn, kb1:kbn), intent(inout) :: f 31 PetscReal, dimension(n, imax, jmax, kmax) :: dfdt_vdp 32 PetscReal, parameter :: mu = 1.4, one = 1.0 33 ! 34 dfdt_vdp(1, :, :, :) = f(2, 1, 1, 1) 35 dfdt_vdp(2, :, :, :) = mu*(one - f(1, 1, 1, 1)**2)*f(2, 1, 1, 1) - f(1, 1, 1, 1) 36 end function dfdt_vdp 37 ! 38 ! The standard Forward Euler time-stepping method. 39 ! 40 recursive subroutine forw_euler(t, dt, ib1, ibn, jb1, jbn, kb1, kbn, imax, jmax, kmax, neq, y, dfdt) 41 PetscReal, intent(in) :: t, dt 42 PetscInt, intent(in) :: ib1, ibn, jb1, jbn, kb1, kbn, imax, jmax, kmax, neq 43 PetscReal, dimension(neq, ib1:ibn, jb1:jbn, kb1:kbn), intent(inout) :: y 44 ! 45 ! Define the right-hand side function 46 ! 47 interface 48 function dfdt(t, dt, ib1, ibn, jb1, jbn, kb1, kbn, imax, jmax, kmax, n, f) 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