program mwe_periodic_3d

    use petscsys
    use petscdmplex
#include <petsc/finclude/petscsys.h>
#include <petsc/finclude/petscdmplex.h>

    implicit none

    PetscErrorCode :: ierr
    PetscInt, parameter :: ndim = 3
    PetscInt, parameter :: nvar = 1
    PetscInt, dimension(ndim) :: faces = (/16, 16, 16/)
    PetscReal, dimension(ndim) :: lower = (/-5.0, -5.0, -5.0/)
    PetscReal, dimension(ndim) :: upper = (/5.0, 5.0, 5.0/)
    DMBoundaryType, dimension(ndim) :: periodicity = (/DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC/)
    DM :: dm, dmghost
    PetscViewer :: vtkviewer
    PetscFV :: fvm
    PetscDS :: prob
    Vec :: sol
    PetscScalar, dimension(1) :: context

    call PetscInitialize(PETSC_NULL_CHARACTER,ierr);CHKERRA(ierr)

    call DMPlexCreateBoxMesh(PETSC_COMM_WORLD,ndim,PETSC_FALSE,faces,lower,upper,periodicity,PETSC_TRUE,dm,ierr);CHKERRA(ierr)
    call DMLocalizeCoordinates(dm,ierr);CHKERRA(ierr)
    call DMSetFromOptions(dm,ierr);CHKERRA(ierr)
    call DMPlexConstructGhostCells(dm,PETSC_NULL_CHARACTER,PETSC_NULL_INTEGER,dmghost,ierr);CHKERRA(ierr)
    if (dmghost /= PETSC_NULL_DM) then
        call DMDestroy(dm,ierr);CHKERRA(ierr)
        dm = dmghost
    end if

    call PetscViewerCreate(PETSC_COMM_WORLD,vtkviewer,ierr);CHKERRA(ierr)
    call PetscViewerSetType(vtkviewer,PETSCVIEWERVTK,ierr);CHKERRA(ierr)
    call PetscViewerFileSetName(vtkviewer,"initmesh.vtu",ierr);CHKERRA(ierr)
    call DMView(dm,vtkviewer,ierr);CHKERRA(ierr)
    call PetscViewerDestroy(vtkviewer,ierr);CHKERRA(ierr)

    call PetscFVCreate(PETSC_COMM_WORLD,fvm,ierr);CHKERRA(ierr)
    call PetscFVSetFromOptions(fvm,ierr);CHKERRA(ierr)
    call PetscFVSetNumComponents(fvm,nvar,ierr);CHKERRA(ierr)
    call PetscFVSetSpatialDimension(fvm,ndim,ierr);CHKERRA(ierr)
    call PetscObjectSetName(fvm,"FV solver",ierr);CHKERRA(ierr)
    call PetscFVSetComponentName(fvm,0,"Density",ierr);CHKERRA(ierr)

    call DMAddField(dm,PETSC_NULL_DMLABEL,fvm,ierr);CHKERRA(ierr)
    call DMCreateDS(dm,ierr);CHKERRA(ierr)
    call DMGetDS(dm,prob,ierr);CHKERRA(ierr)
    call PetscDSSetFromOptions(prob,ierr);CHKERRA(ierr)

    call DMCreateGlobalVector(dm,sol,ierr);CHKERRA(ierr)
    call VecZeroEntries(sol,ierr);CHKERRA(ierr)
    call PetscObjectSetNAme(sol,"Solution",ierr);CHKERRA(ierr)
    call DMProjectFunction(dm,0.d0,vortex,context,INSERT_ALL_VALUES,sol,ierr);CHKERRA(ierr)

    call PetscViewerCreate(PETSC_COMM_WORLD,vtkviewer,ierr);CHKERRA(ierr)
    call PetscViewerSetType(vtkviewer,PETSCVIEWERVTK,ierr);CHKERRA(ierr)
    call PetscViewerFileSetName(vtkviewer,"solution.vtu",ierr);CHKERRA(ierr)
    call VecView(sol,vtkviewer,ierr);CHKERRA(ierr)
    call PetscViewerDestroy(vtkviewer,ierr);CHKERRA(ierr)

    call PetscFinalize(ierr);CHKERRA(ierr)

    contains

        subroutine vortex(dim,time,x,Nf,u,ctx,ierr)

            PetscInt                         :: dim, Nf
            PetscReal                        :: time
            PetscReal, dimension(0:dim-1)    :: x
            PetscScalar, dimension(1)        :: u
            PetscScalar, dimension(:)        :: ctx
            PetscErrorCode, intent(out)      :: ierr
            PetscReal                        :: norm_vel, density, &
            &                                   pi, radius, beta, &
            &                                   gamma

            gamma    = 1.4d0
            beta     = 5.0d0
            pi       = 4.d0*atan(1.d0)
            radius   = x(0)**2 + x(1)**2
            density  = 1.0d0*(1.d0 - ((gamma-1.d0)*beta*beta)/(8.d0*gamma*pi*pi) * exp(1-radius))**(1.d0/(gamma-1.d0))

            ! Fill in the vector with the right conservative variables
            u(1) = density        ! rho

        end subroutine vortex

end program mwe_periodic_3d