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:   Vec            u_f;
 13:   DMPolytopeType ct;
 14:   PetscInt       dim, Nc = 1, faces[3];
 15:   PetscInt       Np = 10, field = 0, zero = 0, bs, cStart;
 16:   PetscReal      energy_0 = 0, energy_1 = 0;
 17:   PetscReal      lo[3], hi[3], h[3];
 18:   PetscBool      removePoints = PETSC_TRUE;
 19:   PetscReal     *wq, *coords;
 20:   PetscDataType  dtype;

 22:   PetscFunctionBeginUser;
 23:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 24:   /* Create a mesh */
 25:   PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
 26:   PetscCall(DMSetType(dm, DMPLEX));
 27:   PetscCall(DMSetFromOptions(dm));
 28:   PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));

 30:   PetscCall(DMGetDimension(dm, &dim));
 31:   bs = dim;
 32:   PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-dm_plex_box_faces", faces, &bs, NULL));
 33:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-np", &Np, NULL));
 34:   PetscCall(DMGetBoundingBox(dm, lo, hi));
 35:   for (PetscInt i = 0; i < dim; ++i) {
 36:     h[i] = (hi[i] - lo[i]) / faces[i];
 37:     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]));
 38:   }
 39:   // Create FE space
 40:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
 41:   PetscCall(DMPlexGetCellType(dm, cStart, &ct));
 42:   PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, Nc, ct, NULL, PETSC_DECIDE, &fe));
 43:   PetscCall(PetscFESetFromOptions(fe));
 44:   PetscCall(PetscObjectSetName((PetscObject)fe, "fe"));
 45:   PetscCall(DMSetField(dm, field, NULL, (PetscObject)fe));
 46:   PetscCall(DMCreateDS(dm));
 47:   PetscCall(PetscFEDestroy(&fe));
 48:   PetscCall(DMCreateGlobalVector(dm, &u_f));
 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 (PetscInt p = 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"));

 73:   // Project between particles and continuum field
 74:   const char *fieldnames[1] = {"w_q"};
 75:   Vec         fields[1]     = {u_f};
 76:   PetscCall(DMSwarmProjectFields(sw, 1, fieldnames, fields, SCATTER_FORWARD));
 77:   PetscCall(DMSwarmProjectFields(sw, 1, fieldnames, fields, SCATTER_REVERSE));

 79:   // Compute energy
 80:   PetscCall(DMSwarmGetField(sw, "w_q", &bs, &dtype, (void **)&wq));
 81:   PetscCall(DMSwarmGetField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords));
 82:   for (PetscInt p = 0; p < Np; ++p) energy_1 += wq[p] * (PetscSqr(coords[p * 2 + 0]) + PetscSqr(coords[p * 2 + 1]));
 83:   PetscCall(DMSwarmRestoreField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords));
 84:   PetscCall(DMSwarmRestoreField(sw, "w_q", &bs, &dtype, (void **)&wq));
 85:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "Energy = %20.12e error = %20.12e\n", (double)energy_0, (double)((energy_1 - energy_0) / energy_0)));

 87:   // Cleanup
 88:   PetscCall(VecDestroy(&u_f));
 89:   PetscCall(DMDestroy(&sw));
 90:   PetscCall(DMDestroy(&dm));
 91:   PetscCall(PetscFinalize());
 92:   return 0;
 93: }

 95: /*TEST

 97:   build:
 98:     requires: !complex

100:   test:
101:     suffix: 0
102:     requires: double triangle
103:     args: -dm_plex_simplex 0 -dm_plex_box_faces 4,2 -dm_plex_box_lower -2.0,0.0 -dm_plex_box_upper 2.0,2.0 \
104:            -np 50 -petscspace_degree 2 \
105:            -ptof_ksp_type cg -ptof_pc_type ilu -ptof_ksp_rtol 1.e-14 \
106:            -ftop_ksp_type lsqr -ftop_pc_type none -ftop_ksp_rtol 1.e-14 \
107:            -dm_view -swarm_view
108:     filter: grep -v DM_ | grep -v atomic

110:   test:
111:     suffix: bjacobi
112:     requires: double triangle
113:     args: -dm_plex_simplex 0 -dm_plex_box_faces 4,2 -dm_plex_box_lower -2.0,0.0 -dm_plex_box_upper 2.0,2.0 \
114:           -np 50 -petscspace_degree 2 -dm_plex_hash_location \
115:           -ptof_ksp_type cg -ptof_pc_type ilu -ptof_ksp_rtol 1.e-14 \
116:           -ftop_ksp_type lsqr -ftop_pc_type bjacobi -ftop_sub_pc_type lu -ftop_sub_pc_factor_shift_type nonzero \
117:           -dm_view -swarm_view -ftop_ksp_rtol 1.e-14
118:     filter: grep -v DM_ | grep -v atomic

120: TEST*/