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