1 static char help[] = "Example program demonstrating projection between particle and finite element spaces\n\n"; 2 3 #include "petscdmplex.h" 4 #include "petscds.h" 5 #include "petscdmswarm.h" 6 #include "petscksp.h" 7 8 int main(int argc, char **argv) 9 { 10 DM dm, sw; 11 PetscFE fe; 12 KSP ksp; 13 PC pc; 14 Mat M_p, PM_p = NULL; 15 Vec f, rho, rhs; 16 PetscInt dim, Nc = 1, timestep = 0, i, faces[3]; 17 PetscInt Np = 10, p, field = 0, zero = 0, bs; 18 PetscReal time = 0.0, norm, energy_0, energy_1; 19 PetscReal lo[3], hi[3], h[3]; 20 PetscBool removePoints = PETSC_TRUE; 21 PetscReal *wq, *coords; 22 PetscDataType dtype; 23 PetscBool is_bjac; 24 25 PetscFunctionBeginUser; 26 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 27 /* Create a mesh */ 28 PetscCall(DMCreate(PETSC_COMM_WORLD, &dm)); 29 PetscCall(DMSetType(dm, DMPLEX)); 30 PetscCall(DMSetFromOptions(dm)); 31 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 32 33 PetscCall(DMGetDimension(dm, &dim)); 34 i = dim; 35 PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-dm_plex_box_faces", faces, &i, NULL)); 36 PetscCall(PetscOptionsGetInt(NULL, NULL, "-np", &Np, NULL)); 37 PetscCall(DMGetBoundingBox(dm, lo, hi)); 38 for (i = 0; i < dim; i++) { 39 h[i] = (hi[i] - lo[i]) / faces[i]; 40 PetscCall(PetscPrintf(PETSC_COMM_SELF, " lo = %g hi = %g n = %" PetscInt_FMT " h = %g\n", (double)lo[i], (double)hi[i], faces[i], (double)h[i])); 41 } 42 43 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, Nc, PETSC_FALSE, "", PETSC_DECIDE, &fe)); 44 PetscCall(PetscFESetFromOptions(fe)); 45 PetscCall(PetscObjectSetName((PetscObject)fe, "fe")); 46 PetscCall(DMSetField(dm, field, NULL, (PetscObject)fe)); 47 PetscCall(DMCreateDS(dm)); 48 PetscCall(PetscFEDestroy(&fe)); 49 /* Create particle swarm */ 50 PetscCall(DMCreate(PETSC_COMM_SELF, &sw)); 51 PetscCall(DMSetType(sw, DMSWARM)); 52 PetscCall(DMSetDimension(sw, dim)); 53 PetscCall(DMSwarmSetType(sw, DMSWARM_PIC)); 54 PetscCall(DMSwarmSetCellDM(sw, dm)); 55 PetscCall(DMSwarmRegisterPetscDatatypeField(sw, "w_q", Nc, PETSC_SCALAR)); 56 PetscCall(DMSwarmFinalizeFieldRegister(sw)); 57 PetscCall(DMSwarmSetLocalSizes(sw, Np, zero)); 58 PetscCall(DMSetFromOptions(sw)); 59 PetscCall(DMSwarmGetField(sw, "w_q", &bs, &dtype, (void **)&wq)); 60 PetscCall(DMSwarmGetField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords)); 61 for (p = 0, energy_0 = 0; p < Np; p++) { 62 coords[p * 2 + 0] = -PetscCosReal((PetscReal)(p + 1) / (PetscReal)(Np + 1) * PETSC_PI); 63 coords[p * 2 + 1] = PetscSinReal((PetscReal)(p + 1) / (PetscReal)(Np + 1) * PETSC_PI); 64 wq[p] = 1.0; 65 energy_0 += wq[p] * (PetscSqr(coords[p * 2 + 0]) + PetscSqr(coords[p * 2 + 1])); 66 } 67 PetscCall(DMSwarmRestoreField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords)); 68 PetscCall(DMSwarmRestoreField(sw, "w_q", &bs, &dtype, (void **)&wq)); 69 PetscCall(DMSwarmMigrate(sw, removePoints)); 70 PetscCall(PetscObjectSetName((PetscObject)sw, "Particle Grid")); 71 PetscCall(DMViewFromOptions(sw, NULL, "-swarm_view")); 72 73 /* Project particles to field */ 74 /* This gives M f = \int_\Omega \phi f, which looks like a rhs for a PDE */ 75 PetscCall(DMCreateMassMatrix(sw, dm, &M_p)); 76 PetscCall(DMCreateGlobalVector(dm, &rho)); 77 PetscCall(PetscObjectSetName((PetscObject)rho, "rho")); 78 79 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f)); 80 PetscCall(PetscObjectSetName((PetscObject)f, "weights")); 81 PetscCall(MatMultTranspose(M_p, f, rho)); 82 83 /* Visualize mesh field */ 84 PetscCall(DMSetOutputSequenceNumber(dm, timestep, time)); 85 PetscCall(VecViewFromOptions(rho, NULL, "-rho_view")); 86 87 /* Project field to particles */ 88 /* This gives f_p = M_p^+ M f */ 89 PetscCall(DMCreateGlobalVector(dm, &rhs)); 90 PetscCall(VecCopy(rho, rhs)); /* Identity: M^1 M rho */ 91 92 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp)); 93 PetscCall(KSPSetOptionsPrefix(ksp, "ftop_")); 94 PetscCall(KSPSetFromOptions(ksp)); 95 PetscCall(KSPGetPC(ksp, &pc)); 96 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCBJACOBI, &is_bjac)); 97 if (is_bjac) { 98 PetscCall(DMSwarmCreateMassMatrixSquare(sw, dm, &PM_p)); 99 PetscCall(KSPSetOperators(ksp, M_p, PM_p)); 100 } else { 101 PetscCall(KSPSetOperators(ksp, M_p, M_p)); 102 } 103 PetscCall(KSPSolveTranspose(ksp, rhs, f)); 104 PetscCall(KSPDestroy(&ksp)); 105 PetscCall(VecDestroy(&rhs)); 106 107 /* Visualize particle field */ 108 PetscCall(DMSetOutputSequenceNumber(sw, timestep, time)); 109 PetscCall(VecViewFromOptions(f, NULL, "-weights_view")); 110 PetscCall(VecNorm(f, NORM_1, &norm)); 111 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f)); 112 113 /* compute energy */ 114 PetscCall(DMSwarmGetField(sw, "w_q", &bs, &dtype, (void **)&wq)); 115 PetscCall(DMSwarmGetField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords)); 116 for (p = 0, energy_1 = 0; p < Np; p++) energy_1 += wq[p] * (PetscSqr(coords[p * 2 + 0]) + PetscSqr(coords[p * 2 + 1])); 117 PetscCall(DMSwarmRestoreField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords)); 118 PetscCall(DMSwarmRestoreField(sw, "w_q", &bs, &dtype, (void **)&wq)); 119 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Total number = %20.12e. energy = %20.12e error = %20.12e\n", (double)norm, (double)energy_0, (double)((energy_1 - energy_0) / energy_0))); 120 /* Cleanup */ 121 PetscCall(MatDestroy(&M_p)); 122 PetscCall(MatDestroy(&PM_p)); 123 PetscCall(VecDestroy(&rho)); 124 PetscCall(DMDestroy(&sw)); 125 PetscCall(DMDestroy(&dm)); 126 PetscCall(PetscFinalize()); 127 return 0; 128 } 129 130 /*TEST 131 132 build: 133 requires: !complex 134 135 test: 136 suffix: 0 137 requires: double triangle 138 args: -dm_plex_simplex 0 -dm_plex_box_faces 4,2 -np 50 -dm_plex_box_lower -2.0,0.0 -dm_plex_box_upper 2.0,2.0 -petscspace_degree 2 -ftop_ksp_type lsqr -ftop_pc_type none -dm_view -swarm_view -ftop_ksp_rtol 1.e-14 139 filter: grep -v DM_ | grep -v atomic 140 141 test: 142 suffix: bjacobi 143 requires: double triangle 144 args: -dm_plex_simplex 0 -dm_plex_box_faces 4,2 -np 50 -dm_plex_box_lower -2.0,0.0 -dm_plex_box_upper 2.0,2.0 -petscspace_degree 2 -dm_plex_hash_location -ftop_ksp_type lsqr -ftop_pc_type bjacobi -ftop_sub_pc_type lu -ftop_sub_pc_factor_shift_type nonzero -dm_view -swarm_view -ftop_ksp_rtol 1.e-14 145 filter: grep -v DM_ | grep -v atomic 146 147 TEST*/ 148