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