1#include <petsc/finclude/petscdmda.h> 2program main 3! 4! This example intends to show how DMDA is used to solve a PDE on a decomposed 5! domain. The equation we are solving is not a PDE, but a toy example: van der 6! Pol's 2-variable ODE duplicated onto a 3D grid: 7! dx/dt = y 8! dy/dt = mu(1-x**2)y - x 9! 10! So we are solving the same equation on all grid points, with no spatial 11! dependencies. Still we tell PETSc to communicate (stencil width >0) so we 12! have communication between different parts of the domain. 13! 14! The example is structured so that one can replace the RHS function and 15! the forw_euler routine with a suitable RHS and a suitable time-integration 16! scheme and make little or no modifications to the DMDA parts. In particular, 17! the "inner" parts of the RHS and time-integration do not "know about" the 18! decomposed domain. 19! 20! See: http://dx.doi.org/10.6084/m9.figshare.1368581 21! 22! Contributed by Aasmund Ervik (asmunder at pvv.org) 23! 24 use ex13f90auxmodule 25 26 use petscdmda 27 28 PetscErrorCode ierr 29 PetscMPIInt rank, size 30 MPIU_Comm comm 31 Vec Lvec, coords 32 DM SolScal, CoordDM 33 DMBoundaryType b_x, b_y, b_z 34 PetscReal, pointer :: array(:, :, :, :) 35 PetscReal :: t, tend, dt, xmin, xmax, ymin, ymax, zmin, zmax, xgmin, xgmax, ygmin, ygmax, zgmin, zgmax 36 PetscReal, allocatable :: f(:, :, :, :), grid(:, :, :, :) 37 PetscInt :: i, j, k, igmax, jgmax, kgmax, ib1, ibn, jb1, jbn, kb1, kbn, imax, jmax, kmax, itime, maxstep, nscreen, dof, stw, ndim 38 39 ! Fire up PETSc: 40 PetscCallA(PetscInitialize(ierr)) 41 42 comm = PETSC_COMM_WORLD 43 PetscCallMPIA(MPI_Comm_rank(comm, rank, ierr)) 44 PetscCallMPIA(MPI_Comm_size(comm, size, ierr)) 45 if (rank == 0) then 46 write (*, *) 'Hi! We are solving van der Pol using ', size, ' processes.' 47 write (*, *) ' ' 48 write (*, *) ' t x1 x2' 49 end if 50 51 ! Set up the global grid 52 igmax = 50 53 jgmax = 50 54 kgmax = 50 55 xgmin = 0.0 56 ygmin = 0.0 57 zgmin = 0.0 58 xgmax = 1.0 59 ygmax = 1.0 60 zgmax = 1.0 61 stw = 1 ! stencil width 62 dof = 2 ! number of variables in this DA 63 ndim = 3 ! 3D code 64 65 ! Get the BCs and create the DMDA 66 call get_boundary_cond(b_x, b_y, b_z) 67 PetscCallA(DMDACreate3d(comm, b_x, b_y, b_z, DMDA_STENCIL_STAR, igmax, jgmax, kgmax, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, dof, stw, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, SolScal, ierr)) 68 PetscCallA(DMSetFromOptions(SolScal, ierr)) 69 PetscCallA(DMSetUp(SolScal, ierr)) 70 71 ! Set global coordinates, get a global and a local work vector 72 PetscCallA(DMDASetUniformCoordinates(SolScal, xgmin, xgmax, ygmin, ygmax, zgmin, zgmax, ierr)) 73 PetscCallA(DMCreateLocalVector(SolScal, Lvec, ierr)) 74 75 ! Get ib1,imax,ibn etc. of the local grid. 76 ! Our convention is: 77 ! the first local ghost cell is ib1 78 ! the first local cell is 1 79 ! the last local cell is imax 80 ! the last local ghost cell is ibn. 81 ! 82 ! i,j,k must be in this call, but are not used 83 PetscCallA(DMDAGetCorners(SolScal, i, j, k, imax, jmax, kmax, ierr)) 84 ib1 = 1 - stw 85 jb1 = 1 - stw 86 kb1 = 1 - stw 87 ibn = imax + stw 88 jbn = jmax + stw 89 kbn = kmax + stw 90 allocate (f(dof, ib1:ibn, jb1:jbn, kb1:kbn)) 91 allocate (grid(ndim, ib1:ibn, jb1:jbn, kb1:kbn)) 92 93 ! Get xmin,xmax etc. for the local grid 94 ! The "coords" local vector here is borrowed, so we shall not destroy it. 95 PetscCallA(DMGetCoordinatesLocal(SolScal, coords, ierr)) 96 ! We need a new DM for coordinate stuff since PETSc supports unstructured grid 97 PetscCallA(DMGetCoordinateDM(SolScal, CoordDM, ierr)) 98 ! petsc_to_local and local_to_petsc are convenience functions, see 99 ! ex13f90auxmodule.F90. 100 call petsc_to_local(CoordDM, coords, array, grid, ndim, stw) 101 xmin = grid(1, 1, 1, 1) 102 ymin = grid(2, 1, 1, 1) 103 zmin = grid(3, 1, 1, 1) 104 xmax = grid(1, imax, jmax, kmax) 105 ymax = grid(2, imax, jmax, kmax) 106 zmax = grid(3, imax, jmax, kmax) 107 call local_to_petsc(CoordDM, coords, array, grid, ndim, stw) 108 109 ! Note that we never use xmin,xmax in this example, but the preceding way of 110 ! getting the local xmin,xmax etc. from PETSc for a structured uniform grid 111 ! is not documented in any other examples I could find. 112 113 ! Set up the time-stepping 114 t = 0.0 115 tend = 100.0 116 dt = 1e-3 117 maxstep = ceiling((tend - t)/dt) 118 ! Write output every second (of simulation-time) 119 nscreen = int(1.0/dt) + 1 120 121 ! Set initial condition 122 PetscCallA(DMDAVecGetArray(SolScal, Lvec, array, ierr)) 123 array(0, :, :, :) = 0.5 124 array(1, :, :, :) = 0.5 125 PetscCallA(DMDAVecRestoreArray(SolScal, Lvec, array, ierr)) 126 127 ! Initial set-up finished. 128 ! Time loop 129 maxstep = 5 130 do itime = 1, maxstep 131 132 ! Communicate such that everyone has the correct values in ghost cells 133 PetscCallA(DMLocalToLocalBegin(SolScal, Lvec, INSERT_VALUES, Lvec, ierr)) 134 PetscCallA(DMLocalToLocalEnd(SolScal, Lvec, INSERT_VALUES, Lvec, ierr)) 135 136 ! Get the old solution from the PETSc data structures 137 call petsc_to_local(SolScal, Lvec, array, f, dof, stw) 138 139 ! Do the time step 140 call forw_euler(t, dt, ib1, ibn, jb1, jbn, kb1, kbn, imax, jmax, kmax, dof, f, dfdt_vdp) 141 t = t + dt 142 143 ! Write result to screen (if main process and it's time to) 144 if (rank == 0 .and. mod(itime, nscreen) == 0) then 145 write (*, 101) t, f(1, 1, 1, 1), f(2, 1, 1, 1) 146 end if 147 148 ! Put our new solution in the PETSc data structures 149 call local_to_petsc(SolScal, Lvec, array, f, dof, stw) 150 end do 151 152 ! Deallocate and finalize 153 PetscCallA(DMRestoreLocalVector(SolScal, Lvec, ierr)) 154 PetscCallA(DMDestroy(SolScal, ierr)) 155 deallocate (f) 156 deallocate (grid) 157 PetscCallA(PetscFinalize(ierr)) 158 159 ! Format for writing output to screen 160101 format(F5.1, 2F11.6) 161 162end program main 163 164!/*TEST 165! 166! build: 167! requires: !complex 168! depends: ex13f90aux.F90 169! 170! test: 171! nsize: 4 172! 173!TEST*/ 174