Actual source code: ex13f90aux.F90

petsc-3.11.3 2019-06-26
Report Typos and Errors
```  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
```