1: module ex13f90aux
2: implicit none
3: contains
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);
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);
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
104: end module ex13f90aux