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: DMBoundaryType,intent(inout) :: b_x,b_y,b_z
10: 11: ! Here you may set the BC types you want
12: b_x = DM_BOUNDARY_GHOSTED
13: b_y = DM_BOUNDARY_GHOSTED
14: b_z = DM_BOUNDARY_GHOSTED
15: 16: end subroutine get_boundary_cond
17: !
18: ! A function which returns the RHS of the equation we are solving
19: !
20: function dfdt_vdp(t,dt,ib1,ibn,jb1,jbn,kb1,kbn,imax,jmax,kmax,n,f)
21: !
22: ! Right-hand side for the van der Pol oscillator. Very simple system of two
23: ! ODEs. See Iserles, eq (5.2).
24: !
25: PetscReal, intent(in) :: t,dt
26: PetscInt, intent(in) :: ib1,ibn,jb1,jbn,kb1,kbn,imax,jmax,kmax,n
27: PetscReal, dimension(n,ib1:ibn,jb1:jbn,kb1:kbn), intent(inout) :: f
28: PetscReal, dimension(n,imax,jmax,kmax) :: dfdt_vdp
29: PetscReal, parameter :: mu=1.4, one=1.0
30: !
31: dfdt_vdp(1,:,:,:) = f(2,1,1,1)
32: dfdt_vdp(2,:,:,:) = mu*(one - f(1,1,1,1)**2)*f(2,1,1,1) - f(1,1,1,1)
33: end function dfdt_vdp
34: !
35: ! The standard Forward Euler time-stepping method.
36: !
37: recursive subroutine forw_euler(t,dt,ib1,ibn,jb1,jbn,kb1,kbn,&
38: 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: #include <petsc/finclude/petscsys.h>
69: #include <petsc/finclude/petscvec.h>
70: #include <petsc/finclude/petscdmda.h>
71: #include <petsc/finclude/petscvec.h90>
72: #include <petsc/finclude/petscdmda.h90>
73: DM :: da
74: Vec,intent(in) :: vec
75: PetscReal, pointer :: array(:,:,:,:)
76: PetscInt,intent(in) :: dof,stw
77: PetscReal, intent(inout), dimension(:,1-stw:,1-stw:,1-stw:) :: f
78: PetscErrorCode :: ierr
79: !
80: call DMDAVecGetArrayF90(da,vec,array,ierr)
81: call transform_petsc_us(array,f,stw)
82: end subroutine petsc_to_local
83: subroutine transform_petsc_us(array,f,stw)
84: !Note: this assumed shape-array is what does the "coordinate transformation"
85: PetscInt,intent(in) :: stw
86: PetscReal, intent(in), dimension(:,1-stw:,1-stw:,1-stw:) :: array
87: PetscReal,intent(inout),dimension(:,1-stw:,1-stw:,1-stw:) :: f
88: f(:,:,:,:) = array(:,:,:,:)
89: end subroutine transform_petsc_us
90: subroutine local_to_petsc(da,vec,array,f,dof,stw)
91: #include <petsc/finclude/petscsys.h>
92: #include <petsc/finclude/petscvec.h>
93: #include <petsc/finclude/petscdmda.h>
94: #include <petsc/finclude/petscvec.h90>
95: #include <petsc/finclude/petscdmda.h90>
96: DM :: da
97: Vec,intent(inout) :: vec
98: PetscReal, pointer :: array(:,:,:,:)
99: PetscInt,intent(in) :: dof,stw
100: PetscReal,intent(inout),dimension(:,1-stw:,1-stw:,1-stw:) :: f
101: PetscErrorCode :: ierr
102: call transform_us_petsc(array,f,stw)
103: call DMDAVecRestoreArrayF90(da,vec,array,ierr)
104: end subroutine local_to_petsc
105: subroutine transform_us_petsc(array,f,stw)
106: !Note: this assumed shape-array is what does the "coordinate transformation"
107: PetscInt,intent(in) :: stw
108: PetscReal, intent(inout), dimension(:,1-stw:,1-stw:,1-stw:) :: array
109: PetscReal, intent(in),dimension(:,1-stw:,1-stw:,1-stw:) :: f
110: array(:,:,:,:) = f(:,:,:,:)
111: end subroutine transform_us_petsc
112: end module ex13f90aux