1module ex13f90aux 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 use petscdm 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 PetscReal, intent(in) :: t,dt 48 PetscInt, intent(in) :: ib1,ibn,jb1,jbn,kb1,kbn,imax,jmax,kmax,n 49 PetscReal, dimension(n,ib1:ibn,jb1:jbn,kb1:kbn), intent(inout) :: f 50 PetscReal, dimension(n,imax,jmax,kmax) :: dfdt 51 end function dfdt 52 end interface 53 !-------------------------------------------------------------------------- 54 ! 55 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) 56 end subroutine forw_euler 57 ! 58 ! The following 4 subroutines handle the mapping of coordinates. I'll explain 59 ! this in detail: 60 ! PETSc gives you local arrays which are indexed using the global indices. 61 ! This is probably handy in some cases, but when you are re-writing an 62 ! existing serial code and want to use DMDAs, you have tons of loops going 63 ! from 1 to imax etc. that you don't want to change. 64 ! These subroutines re-map the arrays so that all the local arrays go from 65 ! 1 to the (local) imax. 66 ! 67 subroutine petsc_to_local(da,vec,array,f,dof,stw) 68 use petscdmda 69 DM :: da 70 Vec,intent(in) :: vec 71 PetscReal, pointer :: array(:,:,:,:) 72 PetscInt,intent(in) :: dof,stw 73 PetscReal, intent(inout), dimension(:,1-stw:,1-stw:,1-stw:) :: f 74 PetscErrorCode :: ierr 75 ! 76 call DMDAVecGetArrayF90(da,vec,array,ierr);CHKERRQ(ierr); 77 call transform_petsc_us(array,f,stw) 78 end subroutine petsc_to_local 79 subroutine transform_petsc_us(array,f,stw) 80 !Note: this assumed shape-array is what does the "coordinate transformation" 81 PetscInt,intent(in) :: stw 82 PetscReal, intent(in), dimension(:,1-stw:,1-stw:,1-stw:) :: array 83 PetscReal,intent(inout),dimension(:,1-stw:,1-stw:,1-stw:) :: f 84 f(:,:,:,:) = array(:,:,:,:) 85 end subroutine transform_petsc_us 86 subroutine local_to_petsc(da,vec,array,f,dof,stw) 87 use petscdmda 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 call DMDAVecRestoreArrayF90(da,vec,array,ierr);CHKERRQ(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 104end module ex13f90aux 105