Actual source code: ex13f90.F90

petsc-3.7.3 2016-08-01
Report Typos and Errors
  1: program main
  2: !
  3: ! This example intends to show how DMDA is used to solve a PDE on a decomposed
  4: ! domain. The equation we are solving is not a PDE, but a toy example: van der
  5: ! Pol's 2-variable ODE duplicated onto a 3D grid:
  6: ! dx/dt = y
  7: ! dy/dt = mu(1-x**2)y - x
  8: !
  9: ! So we are solving the same equation on all grid points, with no spatial
 10: ! dependencies. Still we tell PETSc to communicate (stencil width >0) so we
 11: ! have communication between different parts of the domain.
 12: !
 13: ! The example is structured so that one can replace the RHS function and
 14: ! the forw_euler routine with a suitable RHS and a suitable time-integration
 15: ! scheme and make little or no modifications to the DMDA parts. In particular,
 16: ! the "inner" parts of the RHS and time-integration do not "know about" the
 17: ! decomposed domain.
 18: !
 19: !     See:     http://dx.doi.org/10.6084/m9.figshare.1368581
 20: !
 21: !     Contributed by Aasmund Ervik (asmunder at pvv.org)
 22: !
 23:   use ex13f90aux
 24:   implicit none
 25: #include <petsc/finclude/petscsys.h>
 26: #include <petsc/finclude/petscvec.h>
 27: #include <petsc/finclude/petscdmda.h>
 28: #include <petsc/finclude/petscvec.h90>
 29: #include <petsc/finclude/petscdmda.h90>
 30:   PetscErrorCode   ierr
 31:   PetscMPIInt      rank,size
 32:   MPI_Comm         comm
 33:   Vec              Lvec,coords
 34:   DM               SolScal,CoordDM
 35:   DMBoundaryType b_x,b_y,b_z
 36:   PetscReal, pointer :: array(:,:,:,:)
 37:   PetscReal :: t,tend,dt,xmin,xmax,ymin,ymax,zmin,zmax,&
 38:           xgmin,xgmax,ygmin,ygmax,zgmin,zgmax
 39:   PetscReal, allocatable :: f(:,:,:,:), grid(:,:,:,:)
 40:   PetscInt :: i,j,k,igmax,jgmax,kgmax,ib1,ibn,jb1,jbn,kb1,kbn,&
 41:              imax,jmax,kmax,itime,maxstep,nscreen,dof,stw,ndim

 43:   ! Fire up PETSc:
 44:   call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 45:   comm = PETSC_COMM_WORLD
 46:   call MPI_Comm_rank(comm,rank,ierr)
 47:   call MPI_Comm_size(comm,size,ierr)
 48:   if (rank == 0) then
 49:     write(*,*) 'Hi! We are solving van der Pol using ',size,' processes.'
 50:     write(*,*) ' '
 51:     write(*,*) '  t     x1         x2'
 52:   endif

 54:   ! Set up the global grid
 55:   igmax = 50
 56:   jgmax = 50
 57:   kgmax = 50
 58:   xgmin = 0.0
 59:   ygmin = 0.0
 60:   zgmin = 0.0
 61:   xgmax = 1.0
 62:   ygmax = 1.0
 63:   zgmax = 1.0
 64:   stw = 1 ! stencil width
 65:   dof = 2 ! number of variables in this DA
 66:   ndim = 3 ! 3D code

 68:   ! Get the BCs and create the DMDA
 69:   call get_boundary_cond(b_x,b_y,b_z)
 70:   call DMDACreate3d(comm,b_x,b_y,b_z,DMDA_STENCIL_STAR,igmax,jgmax,kgmax,&
 71:                     PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,dof,stw,&
 72:                     PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,&
 73:                     SolScal,ierr)
 74:   ! Set global coordinates, get a global and a local work vector
 75:   call DMDASetUniformCoordinates(SolScal,xgmin,xgmax,ygmin,ygmax,zgmin,zgmax,&
 76:                                  ierr)
 77:   call DMCreateLocalVector(SolScal,Lvec,ierr)
 78: 
 79:   ! Get ib1,imax,ibn etc. of the local grid.
 80:   ! Our convention is:
 81:   ! the first local ghost cell is ib1
 82:   ! the first local       cell is 1
 83:   ! the last  local       cell is imax
 84:   ! the last  local ghost cell is ibn.
 85:   !
 86:   ! i,j,k must be in this call, but are not used
 87:   call DMDAGetCorners(SolScal,i,j,k,imax,jmax,kmax,ierr)
 88:   ib1=1-stw
 89:   jb1=1-stw
 90:   kb1=1-stw
 91:   ibn=imax+stw
 92:   jbn=jmax+stw
 93:   kbn=kmax+stw
 94:   allocate(f(dof,ib1:ibn,jb1:jbn,kb1:kbn))
 95:   allocate(grid(ndim,ib1:ibn,jb1:jbn,kb1:kbn))

 97:   ! Get xmin,xmax etc. for the local grid
 98:   ! The "coords" local vector here is borrowed, so we shall not destroy it.
 99:   call DMGetCoordinatesLocal(SolScal,coords,ierr)
100:   ! We need a new DM for coordinate stuff since PETSc supports unstructured grid
101:   call DMGetCoordinateDM(SolScal,CoordDM,ierr)
102:   ! petsc_to_local and local_to_petsc are convenience functions, see
103:   ! ex13f90aux.F90.
104:   call petsc_to_local(CoordDM,coords,array,grid,ndim,stw)
105:   xmin=grid(1,1,1,1)
106:   ymin=grid(2,1,1,1)
107:   zmin=grid(3,1,1,1)
108:   xmax=grid(1,imax,jmax,kmax)
109:   ymax=grid(2,imax,jmax,kmax)
110:   zmax=grid(3,imax,jmax,kmax)
111:   call local_to_petsc(CoordDM,coords,array,grid,ndim,stw)
112: 
113:   ! Note that we never use xmin,xmax in this example, but the preceding way of
114:   ! getting the local xmin,xmax etc. from PETSc for a structured uniform grid
115:   ! is not documented in any other examples I could find.

117:   ! Set up the time-stepping
118:   t = 0.0
119:   tend = 100.0
120:   dt = 1e-3
121:   maxstep=ceiling((tend-t)/dt)
122:   ! Write output every second (of simulation-time)
123:   nscreen = int(1.0/dt)+1

125:   ! Set initial condition
126:   call DMDAVecGetArrayF90(SolScal,Lvec,array,ierr)
127:   array(0,:,:,:) = 0.5
128:   array(1,:,:,:) = 0.5
129:   call DMDAVecRestoreArrayF90(SolScal,Lvec,array,ierr)
130: 
131:   ! Initial set-up finished.
132:   ! Time loop
133:   maxstep = 5
134:   do itime=1,maxstep

136:     ! Communicate such that everyone has the correct values in ghost cells
137:     call DMLocalToLocalBegin(SolScal,Lvec,INSERT_VALUES,Lvec,ierr)
138:     call DMLocalToLocalEnd(SolScal,Lvec,INSERT_VALUES,Lvec,ierr)
139: 
140:     ! Get the old solution from the PETSc data structures
141:     call petsc_to_local(SolScal,Lvec,array,f,dof,stw)
142: 
143:     ! Do the time step
144:     call forw_euler(t,dt,ib1,ibn,jb1,jbn,kb1,kbn,imax,jmax,kmax,dof,f,dfdt_vdp)
145:     t=t+dt

147:     ! Write result to screen (if main process and it's time to)
148:     if (rank == 0 .and. mod(itime,nscreen) == 0) then
149:       write(*,101) t,f(1,1,1,1),f(2,1,1,1)
150:     endif
151: 
152:     ! Put our new solution in the PETSc data structures
153:     call local_to_petsc(SolScal,Lvec,array,f,dof,stw)
154:   end do
155: 
156:   ! Deallocate and finalize
157:   call DMRestoreLocalVector(SolScal,Lvec,ierr)
158:   call DMDestroy(SolScal,ierr)
159:   deallocate(f)
160:   deallocate(grid)
161:   call PetscFinalize(ierr)

163:   ! Format for writing output to screen
164: 101 format(F5.1,2F11.6)

166: end program main