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 petscdt 9 use petscksp 10 use petscsys 11 implicit none 12 13 DM :: dm, sw 14 PetscFE :: fe 15 KSP :: ksp 16 Mat :: M_p, M 17 Vec :: f, rho, rhs 18 PetscInt :: dim, Nc = 1, degree = 1, timestep = 0 19 PetscInt :: Np = 100, p, field = 0, zero = 0, bs 20 PetscReal :: time = 0.0, norm 21 PetscBool :: removePoints = PETSC_TRUE 22 PetscDataType :: dtype 23 PetscScalar, pointer :: coords(:) 24 PetscScalar, pointer :: wq(:) 25 PetscErrorCode :: ierr 26 27 call PetscInitialize(PETSC_NULL_CHARACTER, ierr) 28 if (ierr .ne. 0) then 29 print*,'Unable to initialize PETSc' 30 stop 31 endif 32 33 call DMCreate(PETSC_COMM_WORLD, dm, ierr);CHKERRA(ierr) 34 call DMSetType(dm, DMPLEX, ierr);CHKERRA(ierr) 35 call DMSetFromOptions(dm, ierr);CHKERRA(ierr) 36 call DMGetDimension(dm, dim, ierr);CHKERRA(ierr) 37 call DMViewFromOptions(dm, PETSC_NULL_VEC, '-dm_view', ierr);CHKERRA(ierr) 38 39! Create finite element space 40 call PetscFECreateLagrange(PETSC_COMM_SELF, dim, Nc, PETSC_FALSE, degree, PETSC_DETERMINE, fe, ierr);CHKERRA(ierr) 41 call DMSetField(dm, field, PETSC_NULL_DMLABEL, fe, ierr);CHKERRA(ierr) 42 call DMCreateDS(dm, ierr);CHKERRA(ierr) 43 call PetscFEDestroy(fe, ierr);CHKERRA(ierr) 44 45! Create particle swarm 46 call DMCreate(PETSC_COMM_WORLD, sw, ierr);CHKERRA(ierr) 47 call DMSetType(sw, DMSWARM, ierr);CHKERRA(ierr) 48 call DMSetDimension(sw, dim, ierr);CHKERRA(ierr) 49 call DMSwarmSetType(sw, DMSWARM_PIC, ierr);CHKERRA(ierr) 50 call DMSwarmSetCellDM(sw, dm, ierr);CHKERRA(ierr) 51 call DMSwarmRegisterPetscDatatypeField(sw, 'w_q', Nc, PETSC_SCALAR, ierr);CHKERRA(ierr) 52 call DMSwarmFinalizeFieldRegister(sw, ierr);CHKERRA(ierr) 53 call DMSwarmSetLocalSizes(sw, Np, zero, ierr);CHKERRA(ierr) 54 call DMSetFromOptions(sw, ierr);CHKERRA(ierr) 55 call DMSwarmGetField(sw, 'w_q', bs, dtype, wq, ierr);CHKERRA(ierr) 56 call DMSwarmGetField(sw, 'DMSwarmPIC_coor', bs, dtype, coords, ierr);CHKERRA(ierr) 57 do p = 1, Np 58 coords(p*2-1) = -cos(dble(p)/dble(Np+1) * PETSC_PI) 59 coords(p*2-0) = sin(dble(p)/dble(Np+1) * PETSC_PI) 60 wq(p) = 1.0 61 end do 62 call DMSwarmRestoreField(sw, 'DMSwarmPIC_coor', bs, dtype, coords, ierr);CHKERRA(ierr) 63 call DMSwarmRestoreField(sw, 'w_q', bs, dtype, wq, ierr);CHKERRA(ierr) 64 call DMSwarmMigrate(sw, removePoints, ierr);CHKERRA(ierr) 65 call DMViewFromOptions(sw, PETSC_NULL_VEC, '-swarm_view', ierr);CHKERRA(ierr) 66 67! Project particles to field 68! This gives M f = \int_\Omega \phi f, which looks like a rhs for a PDE 69 call DMCreateMassMatrix(sw, dm, M_p, ierr);CHKERRA(ierr) 70 call DMCreateGlobalVector(dm, rho, ierr);CHKERRA(ierr) 71 call DMSwarmCreateGlobalVectorFromField(sw, 'w_q', f, ierr);CHKERRA(ierr) 72 call MatMultTranspose(M_p, f, rho, ierr);CHKERRA(ierr) 73 74! Visualize mesh field 75 call DMSetOutputSequenceNumber(dm, timestep, time, ierr);CHKERRA(ierr) 76 call PetscObjectViewFromOptions(rho, PETSC_NULL_VEC, '-rho_view', ierr);CHKERRA(ierr) 77 78! Project field to particles 79! This gives f_p = M_p^+ M f 80 call DMCreateMassMatrix(dm, dm, M, ierr);CHKERRA(ierr) 81 call DMCreateGlobalVector(dm, rhs, ierr);CHKERRA(ierr) 82 if (.false.) then 83 call MatMult(M, rho, rhs, ierr);CHKERRA(ierr) ! this is what you would do for and FE solve 84 else 85 call VecCopy(rho, rhs, ierr);CHKERRA(ierr) ! Identity: M^1 M rho 86 end if 87 call KSPCreate(PETSC_COMM_WORLD, ksp, ierr);CHKERRA(ierr) 88 call KSPSetOptionsPrefix(ksp, 'ftop_', ierr);CHKERRA(ierr) 89 call KSPSetFromOptions(ksp, ierr);CHKERRA(ierr) 90 call KSPSetOperators(ksp, M_p, M_p, ierr);CHKERRA(ierr) 91 call KSPSolveTranspose(ksp, rhs, f, ierr);CHKERRA(ierr) 92 call KSPDestroy(ksp, ierr);CHKERRA(ierr) 93 call VecDestroy(rhs, ierr);CHKERRA(ierr) 94 call MatDestroy(M_p, ierr);CHKERRA(ierr) 95 call MatDestroy(M, ierr);CHKERRA(ierr) 96 97! Visualize particle field 98 call DMSetOutputSequenceNumber(sw, timestep, time, ierr);CHKERRA(ierr) 99 call PetscObjectViewFromOptions(f, PETSC_NULL_VEC, '-weights_view', ierr);CHKERRA(ierr) 100 call VecNorm(f,NORM_1,norm,ierr);CHKERRA(ierr) 101 print *, 'Total number density = ', norm 102! Cleanup 103 call DMSwarmDestroyGlobalVectorFromField(sw, 'w_q', f, ierr);CHKERRA(ierr) 104 call MatDestroy(M_p, ierr);CHKERRA(ierr) 105 call VecDestroy(rho, ierr);CHKERRA(ierr) 106 call DMDestroy(sw, ierr);CHKERRA(ierr) 107 call DMDestroy(dm, ierr);CHKERRA(ierr) 108 109 call PetscFinalize(ierr) 110 end program DMSwarmTestProjection 111 112!/*TEST 113! build: 114! requires: defined(PETSC_USING_F90FREEFORM) double !complex 115! 116! test: 117! suffix: 0 118! requires: double 119! 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 120! filter: grep -v DM_ | grep -v atomic 121! filter_output: grep -v atomic 122! 123!TEST*/ 124