1! Example program demonstrating projection between particle and finite element spaces 2#include <petsc/finclude/petscdmplex.h> 3#include <petsc/finclude/petscdmswarm.h> 4#include <petsc/finclude/petscksp.h> 5program DMSwarmTestProjection 6 use petscdmplex 7 use petscdmswarm 8 use petscksp 9 implicit none 10 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 24 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)) 31 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)) 37 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)) 60 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)) 67 68! Visualize mesh field 69 PetscCallA(DMSetOutputSequenceNumber(dm, timestep, time, ierr)) 70 PetscCallA(PetscObjectViewFromOptions(PetscObjectCast(rho), PETSC_NULL_OBJECT, '-rho_view', ierr)) 71 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)) 90 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)) 101 102 PetscCallA(PetscFinalize(ierr)) 103end program DMSwarmTestProjection 104 105!/*TEST 106! build: 107! requires: 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*/ 117