Actual source code: ex1.c

  1: static char help[] = "Example program demonstrating projection between particle and finite element spaces\n\n";

  3: #include "petscdmplex.h"
  4: #include "petscds.h"
  5: #include "petscdmswarm.h"
  6: #include "petscksp.h"

  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;

 26:   PetscInitialize(&argc, &argv, NULL, help);
 27:   /* Create a mesh */
 28:   DMCreate(PETSC_COMM_WORLD, &dm);
 29:   DMSetType(dm, DMPLEX);
 30:   DMSetFromOptions(dm);
 31:   DMViewFromOptions(dm, NULL, "-dm_view");

 33:   DMGetDimension(dm, &dim);
 34:   i = dim;
 35:   PetscOptionsGetIntArray(NULL, NULL, "-dm_plex_box_faces", faces, &i, NULL);
 36:   PetscOptionsGetInt(NULL, NULL, "-np", &Np, NULL);
 37:   DMGetBoundingBox(dm, lo, hi);
 38:   for (i = 0; i < dim; i++) {
 39:     h[i] = (hi[i] - lo[i]) / faces[i];
 40:     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:   }

 43:   PetscFECreateDefault(PETSC_COMM_SELF, dim, Nc, PETSC_FALSE, "", PETSC_DECIDE, &fe);
 44:   PetscFESetFromOptions(fe);
 45:   PetscObjectSetName((PetscObject)fe, "fe");
 46:   DMSetField(dm, field, NULL, (PetscObject)fe);
 47:   DMCreateDS(dm);
 48:   PetscFEDestroy(&fe);
 49:   /* Create particle swarm */
 50:   DMCreate(PETSC_COMM_SELF, &sw);
 51:   DMSetType(sw, DMSWARM);
 52:   DMSetDimension(sw, dim);
 53:   DMSwarmSetType(sw, DMSWARM_PIC);
 54:   DMSwarmSetCellDM(sw, dm);
 55:   DMSwarmRegisterPetscDatatypeField(sw, "w_q", Nc, PETSC_SCALAR);
 56:   DMSwarmFinalizeFieldRegister(sw);
 57:   DMSwarmSetLocalSizes(sw, Np, zero);
 58:   DMSetFromOptions(sw);
 59:   DMSwarmGetField(sw, "w_q", &bs, &dtype, (void **)&wq);
 60:   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:   DMSwarmRestoreField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords);
 68:   DMSwarmRestoreField(sw, "w_q", &bs, &dtype, (void **)&wq);
 69:   DMSwarmMigrate(sw, removePoints);
 70:   PetscObjectSetName((PetscObject)sw, "Particle Grid");
 71:   DMViewFromOptions(sw, NULL, "-swarm_view");

 73:   /* Project particles to field */
 74:   /* This gives M f = \int_\Omega \phi f, which looks like a rhs for a PDE */
 75:   DMCreateMassMatrix(sw, dm, &M_p);
 76:   DMCreateGlobalVector(dm, &rho);
 77:   PetscObjectSetName((PetscObject)rho, "rho");

 79:   DMSwarmCreateGlobalVectorFromField(sw, "w_q", &f);
 80:   PetscObjectSetName((PetscObject)f, "weights");
 81:   MatMultTranspose(M_p, f, rho);

 83:   /* Visualize mesh field */
 84:   DMSetOutputSequenceNumber(dm, timestep, time);
 85:   VecViewFromOptions(rho, NULL, "-rho_view");

 87:   /* Project field to particles */
 88:   /*   This gives f_p = M_p^+ M f */
 89:   DMCreateGlobalVector(dm, &rhs);
 90:   VecCopy(rho, rhs); /* Identity: M^1 M rho */

 92:   KSPCreate(PETSC_COMM_WORLD, &ksp);
 93:   KSPSetOptionsPrefix(ksp, "ftop_");
 94:   KSPSetFromOptions(ksp);
 95:   KSPGetPC(ksp, &pc);
 96:   PetscObjectTypeCompare((PetscObject)pc, PCBJACOBI, &is_bjac);
 97:   if (is_bjac) {
 98:     DMSwarmCreateMassMatrixSquare(sw, dm, &PM_p);
 99:     KSPSetOperators(ksp, M_p, PM_p);
100:   } else {
101:     KSPSetOperators(ksp, M_p, M_p);
102:   }
103:   KSPSolveTranspose(ksp, rhs, f);
104:   KSPDestroy(&ksp);
105:   VecDestroy(&rhs);

107:   /* Visualize particle field */
108:   DMSetOutputSequenceNumber(sw, timestep, time);
109:   VecViewFromOptions(f, NULL, "-weights_view");
110:   VecNorm(f, NORM_1, &norm);
111:   DMSwarmDestroyGlobalVectorFromField(sw, "w_q", &f);

113:   /* compute energy */
114:   DMSwarmGetField(sw, "w_q", &bs, &dtype, (void **)&wq);
115:   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:   DMSwarmRestoreField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords);
118:   DMSwarmRestoreField(sw, "w_q", &bs, &dtype, (void **)&wq);
119:   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:   MatDestroy(&M_p);
122:   MatDestroy(&PM_p);
123:   VecDestroy(&rho);
124:   DMDestroy(&sw);
125:   DMDestroy(&dm);
126:   PetscFinalize();
127:   return 0;
128: }

130: /*TEST

132:   build:
133:     requires: !complex

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

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

147: TEST*/