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