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