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 = 2, 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 DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, PETSC_FALSE, PETSC_NULL_INTEGER, PETSC_NULL_REAL, PETSC_NULL_REAL, PETSC_NULL_INTEGER, PETSC_TRUE, dm, ierr);CHKERRA(ierr) 34 call DMViewFromOptions(dm, PETSC_NULL_VEC, '-dm_view', ierr);CHKERRA(ierr) 35 36! Create finite element space 37 call PetscFECreateLagrange(PETSC_COMM_SELF, dim, Nc, PETSC_FALSE, degree, PETSC_DETERMINE, fe, ierr);CHKERRA(ierr) 38 call DMSetField(dm, field, PETSC_NULL_DMLABEL, fe, ierr);CHKERRA(ierr) 39 call DMCreateDS(dm, ierr);CHKERRA(ierr) 40 call PetscFEDestroy(fe, ierr);CHKERRA(ierr) 41 42! Create particle swarm 43 call DMCreate(PETSC_COMM_WORLD, sw, ierr);CHKERRA(ierr) 44 call DMSetType(sw, DMSWARM, ierr);CHKERRA(ierr) 45 call DMSetDimension(sw, dim, ierr);CHKERRA(ierr) 46 call DMSwarmSetType(sw, DMSWARM_PIC, ierr);CHKERRA(ierr) 47 call DMSwarmSetCellDM(sw, dm, ierr);CHKERRA(ierr) 48 call DMSwarmRegisterPetscDatatypeField(sw, 'w_q', Nc, PETSC_SCALAR, ierr);CHKERRA(ierr) 49 call DMSwarmFinalizeFieldRegister(sw, ierr);CHKERRA(ierr) 50 call DMSwarmSetLocalSizes(sw, Np, zero, ierr);CHKERRA(ierr) 51 call DMSetFromOptions(sw, ierr);CHKERRA(ierr) 52 call DMSwarmGetField(sw, 'w_q', bs, dtype, wq, ierr);CHKERRA(ierr) 53 call DMSwarmGetField(sw, 'DMSwarmPIC_coor', bs, dtype, coords, ierr);CHKERRA(ierr) 54 do p = 1, Np 55 coords(p*2-1) = -cos(dble(p)/dble(Np+1) * PETSC_PI) 56 coords(p*2-0) = sin(dble(p)/dble(Np+1) * PETSC_PI) 57 wq(p) = 1.0 58 end do 59 call DMSwarmRestoreField(sw, 'DMSwarmPIC_coor', bs, dtype, coords, ierr);CHKERRA(ierr) 60 call DMSwarmRestoreField(sw, 'w_q', bs, dtype, wq, ierr);CHKERRA(ierr) 61 call DMSwarmMigrate(sw, removePoints, ierr);CHKERRA(ierr) 62 call DMViewFromOptions(sw, PETSC_NULL_VEC, '-swarm_view', ierr);CHKERRA(ierr) 63 64! Project particles to field 65! This gives M f = \int_\Omega \phi f, which looks like a rhs for a PDE 66 call DMCreateMassMatrix(sw, dm, M_p, ierr);CHKERRA(ierr) 67 call DMCreateGlobalVector(dm, rho, ierr);CHKERRA(ierr) 68 call DMSwarmCreateGlobalVectorFromField(sw, 'w_q', f, ierr);CHKERRA(ierr) 69 call MatMultTranspose(M_p, f, rho, ierr);CHKERRA(ierr) 70 71! Visualize mesh field 72 call DMSetOutputSequenceNumber(dm, timestep, time, ierr);CHKERRA(ierr) 73 call PetscObjectViewFromOptions(rho, PETSC_NULL_VEC, '-rho_view', ierr);CHKERRA(ierr) 74 75! Project field to particles 76! This gives f_p = M_p^+ M f 77 call DMCreateMassMatrix(dm, dm, M, ierr);CHKERRA(ierr) 78 call DMCreateGlobalVector(dm, rhs, ierr);CHKERRA(ierr) 79 if (.false.) then 80 call MatMult(M, rho, rhs, ierr);CHKERRA(ierr) ! this is what you would do for and FE solve 81 else 82 call VecCopy(rho, rhs, ierr);CHKERRA(ierr) ! Indentity: M^1 M rho 83 end if 84 call KSPCreate(PETSC_COMM_WORLD, ksp, ierr);CHKERRA(ierr) 85 call KSPSetOptionsPrefix(ksp, 'ftop_', ierr);CHKERRA(ierr) 86 call KSPSetFromOptions(ksp, ierr);CHKERRA(ierr) 87 call KSPSetOperators(ksp, M_p, M_p, ierr);CHKERRA(ierr) 88 call KSPSolveTranspose(ksp, rhs, f, ierr);CHKERRA(ierr) 89 call KSPDestroy(ksp, ierr);CHKERRA(ierr) 90 call VecDestroy(rhs, ierr);CHKERRA(ierr) 91 call MatDestroy(M_p, ierr);CHKERRA(ierr) 92 call MatDestroy(M, ierr);CHKERRA(ierr) 93 94! Visualize particle field 95 call DMSetOutputSequenceNumber(sw, timestep, time, ierr);CHKERRA(ierr) 96 call PetscObjectViewFromOptions(f, PETSC_NULL_VEC, '-weights_view', ierr);CHKERRA(ierr) 97 call VecNorm(f,NORM_1,norm,ierr);CHKERRA(ierr) 98 print *, 'Total number density = ', norm 99! Cleanup 100 call DMSwarmDestroyGlobalVectorFromField(sw, 'w_q', f, ierr);CHKERRA(ierr) 101 call MatDestroy(M_p, ierr);CHKERRA(ierr) 102 call VecDestroy(rho, ierr);CHKERRA(ierr) 103 call DMDestroy(sw, ierr);CHKERRA(ierr) 104 call DMDestroy(dm, ierr);CHKERRA(ierr) 105 106 call PetscFinalize(ierr) 107 end program DMSwarmTestProjection 108 109!/*TEST 110! build: 111! requires: define(PETSC_USING_F90FREEFORM) double !complex 112! 113! test: 114! suffix: 0 115! requires: double 116! args: -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 117! filter: grep -v DM_ | grep -v atomic 118! filter_output: grep -v atomic 119! 120!TEST*/ 121