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