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