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