177d968b7SBarry Smithmodule ex13f90auxmodule 2c4762a1bSJed Brown implicit none 3c4762a1bSJed Browncontains 4c4762a1bSJed Brown ! 5c4762a1bSJed Brown ! A subroutine which returns the boundary conditions. 6c4762a1bSJed Brown ! 7c4762a1bSJed Brown subroutine get_boundary_cond(b_x,b_y,b_z) 8c4762a1bSJed Brown#include <petsc/finclude/petscdm.h> 9*ce78bad3SBarry Smith#include <petsc/finclude/petscdmda.h> 10c4762a1bSJed Brown use petscdm 11c4762a1bSJed Brown DMBoundaryType,intent(inout) :: b_x,b_y,b_z 12c4762a1bSJed Brown 13c4762a1bSJed Brown ! Here you may set the BC types you want 14c4762a1bSJed Brown b_x = DM_BOUNDARY_GHOSTED 15c4762a1bSJed Brown b_y = DM_BOUNDARY_GHOSTED 16c4762a1bSJed Brown b_z = DM_BOUNDARY_GHOSTED 17c4762a1bSJed Brown 18c4762a1bSJed Brown end subroutine get_boundary_cond 19c4762a1bSJed Brown ! 20c4762a1bSJed Brown ! A function which returns the RHS of the equation we are solving 21c4762a1bSJed Brown ! 22c4762a1bSJed Brown function dfdt_vdp(t,dt,ib1,ibn,jb1,jbn,kb1,kbn,imax,jmax,kmax,n,f) 23c4762a1bSJed Brown ! 24c4762a1bSJed Brown ! Right-hand side for the van der Pol oscillator. Very simple system of two 25c4762a1bSJed Brown ! ODEs. See Iserles, eq (5.2). 26c4762a1bSJed Brown ! 27c4762a1bSJed Brown PetscReal, intent(in) :: t,dt 28c4762a1bSJed Brown PetscInt, intent(in) :: ib1,ibn,jb1,jbn,kb1,kbn,imax,jmax,kmax,n 29c4762a1bSJed Brown PetscReal, dimension(n,ib1:ibn,jb1:jbn,kb1:kbn), intent(inout) :: f 30c4762a1bSJed Brown PetscReal, dimension(n,imax,jmax,kmax) :: dfdt_vdp 31c4762a1bSJed Brown PetscReal, parameter :: mu=1.4, one=1.0 32c4762a1bSJed Brown ! 33c4762a1bSJed Brown dfdt_vdp(1,:,:,:) = f(2,1,1,1) 34c4762a1bSJed Brown dfdt_vdp(2,:,:,:) = mu*(one - f(1,1,1,1)**2)*f(2,1,1,1) - f(1,1,1,1) 35c4762a1bSJed Brown end function dfdt_vdp 36c4762a1bSJed Brown ! 37c4762a1bSJed Brown ! The standard Forward Euler time-stepping method. 38c4762a1bSJed Brown ! 39c4762a1bSJed Brown recursive subroutine forw_euler(t,dt,ib1,ibn,jb1,jbn,kb1,kbn,imax,jmax,kmax,neq,y,dfdt) 40c4762a1bSJed Brown PetscReal, intent(in) :: t,dt 41c4762a1bSJed Brown PetscInt, intent(in) :: ib1,ibn,jb1,jbn,kb1,kbn,imax,jmax,kmax,neq 42c4762a1bSJed Brown PetscReal, dimension(neq,ib1:ibn,jb1:jbn,kb1:kbn), intent(inout) :: y 43c4762a1bSJed Brown ! 44c4762a1bSJed Brown ! Define the right-hand side function 45c4762a1bSJed Brown ! 46c4762a1bSJed Brown interface 47c4762a1bSJed Brown function dfdt(t,dt,ib1,ibn,jb1,jbn,kb1,kbn,imax,jmax,kmax,n,f) 48c4762a1bSJed Brown PetscReal, intent(in) :: t,dt 49c4762a1bSJed Brown PetscInt, intent(in) :: ib1,ibn,jb1,jbn,kb1,kbn,imax,jmax,kmax,n 50c4762a1bSJed Brown PetscReal, dimension(n,ib1:ibn,jb1:jbn,kb1:kbn), intent(inout) :: f 51c4762a1bSJed Brown PetscReal, dimension(n,imax,jmax,kmax) :: dfdt 52c4762a1bSJed Brown end function dfdt 53c4762a1bSJed Brown end interface 54c4762a1bSJed Brown !-------------------------------------------------------------------------- 55c4762a1bSJed Brown ! 56c4762a1bSJed Brown 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) 57c4762a1bSJed Brown end subroutine forw_euler 58c4762a1bSJed Brown ! 59c4762a1bSJed Brown ! The following 4 subroutines handle the mapping of coordinates. I'll explain 60c4762a1bSJed Brown ! this in detail: 61c4762a1bSJed Brown ! PETSc gives you local arrays which are indexed using the global indices. 62c4762a1bSJed Brown ! This is probably handy in some cases, but when you are re-writing an 63c4762a1bSJed Brown ! existing serial code and want to use DMDAs, you have tons of loops going 64c4762a1bSJed Brown ! from 1 to imax etc. that you don't want to change. 65c4762a1bSJed Brown ! These subroutines re-map the arrays so that all the local arrays go from 66c4762a1bSJed Brown ! 1 to the (local) imax. 67c4762a1bSJed Brown ! 68c4762a1bSJed Brown subroutine petsc_to_local(da,vec,array,f,dof,stw) 69c4762a1bSJed Brown use petscdmda 70c4762a1bSJed Brown DM :: da 71c4762a1bSJed Brown Vec,intent(in) :: vec 72c4762a1bSJed Brown PetscReal, pointer :: array(:,:,:,:) 73c4762a1bSJed Brown PetscInt,intent(in) :: dof,stw 74c4762a1bSJed Brown PetscReal, intent(inout), dimension(:,1-stw:,1-stw:,1-stw:) :: f 75c4762a1bSJed Brown PetscErrorCode :: ierr 76c4762a1bSJed Brown ! 77*ce78bad3SBarry Smith PetscCall(DMDAVecGetArray(da,vec,array,ierr)) 78c4762a1bSJed Brown call transform_petsc_us(array,f,stw) 79c4762a1bSJed Brown end subroutine petsc_to_local 80c4762a1bSJed Brown subroutine transform_petsc_us(array,f,stw) 81c4762a1bSJed Brown !Note: this assumed shape-array is what does the "coordinate transformation" 82c4762a1bSJed Brown PetscInt,intent(in) :: stw 83c4762a1bSJed Brown PetscReal, intent(in), dimension(:,1-stw:,1-stw:,1-stw:) :: array 84c4762a1bSJed Brown PetscReal,intent(inout),dimension(:,1-stw:,1-stw:,1-stw:) :: f 85c4762a1bSJed Brown f(:,:,:,:) = array(:,:,:,:) 86c4762a1bSJed Brown end subroutine transform_petsc_us 87c4762a1bSJed Brown subroutine local_to_petsc(da,vec,array,f,dof,stw) 88c4762a1bSJed Brown use petscdmda 89c4762a1bSJed Brown DM :: da 90c4762a1bSJed Brown Vec,intent(inout) :: vec 91c4762a1bSJed Brown PetscReal, pointer :: array(:,:,:,:) 92c4762a1bSJed Brown PetscInt,intent(in) :: dof,stw 93c4762a1bSJed Brown PetscReal,intent(inout),dimension(:,1-stw:,1-stw:,1-stw:) :: f 94c4762a1bSJed Brown PetscErrorCode :: ierr 95c4762a1bSJed Brown call transform_us_petsc(array,f,stw) 96*ce78bad3SBarry Smith PetscCall(DMDAVecRestoreArray(da,vec,array,ierr)) 97c4762a1bSJed Brown end subroutine local_to_petsc 98c4762a1bSJed Brown subroutine transform_us_petsc(array,f,stw) 99c4762a1bSJed Brown !Note: this assumed shape-array is what does the "coordinate transformation" 100c4762a1bSJed Brown PetscInt,intent(in) :: stw 101c4762a1bSJed Brown PetscReal, intent(inout), dimension(:,1-stw:,1-stw:,1-stw:) :: array 102c4762a1bSJed Brown PetscReal, intent(in),dimension(:,1-stw:,1-stw:,1-stw:) :: f 103c4762a1bSJed Brown array(:,:,:,:) = f(:,:,:,:) 104c4762a1bSJed Brown end subroutine transform_us_petsc 10577d968b7SBarry Smithend module ex13f90auxmodule 106