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