Actual source code: ex13f90aux.F90

petsc-3.6.4 2016-04-12
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:     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, intent(inout)                           :: 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, intent(inout)                   :: 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