1: ! Example program demonstrating projection between particle and finite element spaces
  2: program DMSwarmTestProjection
  3: #include <petsc/finclude/petscdmplex.h>
  4: #include <petsc/finclude/petscdmswarm.h>
  5: #include <petsc/finclude/petscksp.h>
  6:   use petscdmplex
  7:   use petscdmswarm
  8:   use petscksp
  9:   implicit none

 11:   DM ::          dm, sw
 12:   PetscFE ::     fe
 13:   KSP ::         ksp
 14:   Mat ::         M_p, M
 15:   Vec ::         f, rho, rhs
 16:   PetscInt ::    dim, Nc = 1, degree = 1, timestep = 0
 17:   PetscInt ::    Np = 100, p, field = 0, zero = 0, bs
 18:   PetscReal ::   time = 0.0, norm
 19:   PetscBool ::   removePoints = PETSC_TRUE
 20:   PetscDataType :: dtype
 21:   PetscScalar, pointer :: coords(:)
 22:   PetscScalar, pointer :: wq(:)
 23:   PetscErrorCode :: ierr

 25:   PetscCallA(PetscInitialize(PETSC_NULL_CHARACTER, ierr))
 26:   PetscCallA(DMCreate(PETSC_COMM_WORLD, dm, ierr))
 27:   PetscCallA(DMSetType(dm, DMPLEX, ierr))
 28:   PetscCallA(DMSetFromOptions(dm, ierr))
 29:   PetscCallA(DMGetDimension(dm, dim, ierr))
 30:   PetscCallA(DMViewFromOptions(dm, PETSC_NULL_OBJECT, '-dm_view', ierr))

 32: !     Create finite element space
 33:   PetscCallA(PetscFECreateLagrange(PETSC_COMM_SELF, dim, Nc, PETSC_FALSE, degree, PETSC_DETERMINE, fe, ierr))
 34:   PetscCallA(DMSetField(dm, field, PETSC_NULL_DMLABEL, PetscObjectCast(fe), ierr))
 35:   PetscCallA(DMCreateDS(dm, ierr))
 36:   PetscCallA(PetscFEDestroy(fe, ierr))

 38: !     Create particle swarm
 39:   PetscCallA(DMCreate(PETSC_COMM_WORLD, sw, ierr))
 40:   PetscCallA(DMSetType(sw, DMSWARM, ierr))
 41:   PetscCallA(DMSetDimension(sw, dim, ierr))
 42:   PetscCallA(DMSwarmSetType(sw, DMSWARM_PIC, ierr))
 43:   PetscCallA(DMSwarmSetCellDM(sw, dm, ierr))
 44:   PetscCallA(DMSwarmRegisterPetscDatatypeField(sw, 'w_q', Nc, PETSC_SCALAR, ierr))
 45:   PetscCallA(DMSwarmFinalizeFieldRegister(sw, ierr))
 46:   PetscCallA(DMSwarmSetLocalSizes(sw, Np, zero, ierr))
 47:   PetscCallA(DMSetFromOptions(sw, ierr))
 48:   PetscCallA(DMSwarmGetField(sw, 'w_q', bs, dtype, wq, ierr))
 49:   PetscCallA(DMSwarmGetField(sw, 'DMSwarmPIC_coor', bs, dtype, coords, ierr))
 50:   do p = 1, Np
 51:     coords(p*2 - 1) = -cos(dble(p)/dble(Np + 1)*PETSC_PI)
 52:     coords(p*2 - 0) = sin(dble(p)/dble(Np + 1)*PETSC_PI)
 53:     wq(p) = 1.0
 54:   end do
 55:   PetscCallA(DMSwarmRestoreField(sw, 'DMSwarmPIC_coor', bs, dtype, coords, ierr))
 56:   PetscCallA(DMSwarmRestoreField(sw, 'w_q', bs, dtype, wq, ierr))
 57:   PetscCallA(DMSwarmMigrate(sw, removePoints, ierr))
 58:   PetscCallA(DMSwarmVectorDefineField(sw, 'w_q', ierr))
 59:   PetscCallA(DMViewFromOptions(sw, PETSC_NULL_OBJECT, '-swarm_view', ierr))

 61: !     Project particles to field
 62: !       This gives M f = \int_\Omega \phi f, which looks like a rhs for a PDE
 63:   PetscCallA(DMCreateMassMatrix(sw, dm, M_p, ierr))
 64:   PetscCallA(DMCreateGlobalVector(dm, rho, ierr))
 65:   PetscCallA(DMSwarmCreateGlobalVectorFromField(sw, 'w_q', f, ierr))
 66:   PetscCallA(MatMultTranspose(M_p, f, rho, ierr))

 68: !     Visualize mesh field
 69:   PetscCallA(DMSetOutputSequenceNumber(dm, timestep, time, ierr))
 70:   PetscCallA(PetscObjectViewFromOptions(PetscObjectCast(rho), PETSC_NULL_OBJECT, '-rho_view', ierr))

 72: !     Project field to particles
 73: !       This gives f_p = M_p^+ M f
 74:   PetscCallA(DMCreateMassMatrix(dm, dm, M, ierr))
 75:   PetscCallA(DMCreateGlobalVector(dm, rhs, ierr))
 76:   if (.false.) then
 77:     PetscCallA(MatMult(M, rho, rhs, ierr)) ! this is what you would do for and FE solve
 78:   else
 79:     PetscCallA(VecCopy(rho, rhs, ierr)) ! Identity: M^1 M rho
 80:   end if
 81:   PetscCallA(KSPCreate(PETSC_COMM_WORLD, ksp, ierr))
 82:   PetscCallA(KSPSetOptionsPrefix(ksp, 'ftop_', ierr))
 83:   PetscCallA(KSPSetFromOptions(ksp, ierr))
 84:   PetscCallA(KSPSetOperators(ksp, M_p, M_p, ierr))
 85:   PetscCallA(KSPSolveTranspose(ksp, rhs, f, ierr))
 86:   PetscCallA(KSPDestroy(ksp, ierr))
 87:   PetscCallA(VecDestroy(rhs, ierr))
 88:   PetscCallA(MatDestroy(M_p, ierr))
 89:   PetscCallA(MatDestroy(M, ierr))

 91: !     Visualize particle field
 92:   PetscCallA(DMSetOutputSequenceNumber(sw, timestep, time, ierr))
 93:   PetscCallA(PetscObjectViewFromOptions(PetscObjectCast(f), PETSC_NULL_OBJECT, '-weights_view', ierr))
 94:   PetscCallA(VecNorm(f, NORM_1, norm, ierr))
 95:   print *, 'Total number density = ', norm
 96: !     Cleanup
 97:   PetscCallA(DMSwarmDestroyGlobalVectorFromField(sw, 'w_q', f, ierr))
 98:   PetscCallA(VecDestroy(rho, ierr))
 99:   PetscCallA(DMDestroy(sw, ierr))
100:   PetscCallA(DMDestroy(dm, ierr))

102:   PetscCallA(PetscFinalize(ierr))
103: end program DMSwarmTestProjection

105: !/*TEST
106: !  build:
107: !    requires: defined(PETSC_USING_F90FREEFORM) double !complex
108: !
109: !  test:
110: !    suffix: 0
111: !    requires: double
112: !    args: -dm_plex_simplex 0 -dm_plex_box_faces 40,20 -dm_plex_box_lower -2.0,0.0 -dm_plex_box_upper 2.0,2.0 -ftop_ksp_type lsqr -ftop_pc_type none -dm_view -swarm_view
113: !    filter: grep -v DM_ | grep -v atomic
114: !    filter_output: grep -v atomic
115: !
116: !TEST*/