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