Actual source code: swarm_ex1.c

      1: #include "petscis.h"
      2: #include "petscsys.h"
      3: #include "petscsystypes.h"
      4: #include "petscvec.h"
      5: static char help[] = "Tests DMSwarm\n\n";

      7: #include <petscdm.h>
      8: #include <petscdmda.h>
      9: #include <petscdmswarm.h>

     11: PETSC_EXTERN PetscErrorCode DMSwarmCollect_General(DM, PetscErrorCode (*)(DM, void *, PetscInt *, PetscInt **), size_t, void *, PetscInt *);
     12: PETSC_EXTERN PetscErrorCode DMSwarmCollect_DMDABoundingBox(DM, PetscInt *);

     14: PetscErrorCode ex1_1(void)
     15: {
     16:   DM          dms;
     17:   Vec         x;
     18:   PetscMPIInt rank, size;
     19:   PetscInt    p;

     21:   PetscFunctionBegin;
     22:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
     23:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
     24:   PetscCheck(!(size > 1) || !(size != 4), PETSC_COMM_WORLD, PETSC_ERR_SUP, "Must be run with 4 MPI ranks");

     26:   PetscCall(DMCreate(PETSC_COMM_WORLD, &dms));
     27:   PetscCall(DMSetType(dms, DMSWARM));
     28:   PetscCall(PetscObjectSetName((PetscObject)dms, "Particles"));

     30:   PetscCall(DMSwarmInitializeFieldRegister(dms));
     31:   PetscCall(DMSwarmRegisterPetscDatatypeField(dms, "viscosity", 1, PETSC_REAL));
     32:   PetscCall(DMSwarmRegisterPetscDatatypeField(dms, "strain", 1, PETSC_REAL));
     33:   PetscCall(DMSwarmFinalizeFieldRegister(dms));
     34:   PetscCall(DMSwarmSetLocalSizes(dms, 5 + rank, 4));
     35:   PetscCall(DMView(dms, PETSC_VIEWER_STDOUT_WORLD));

     37:   {
     38:     PetscReal *array;
     39:     PetscCall(DMSwarmGetField(dms, "viscosity", NULL, NULL, (void **)&array));
     40:     for (p = 0; p < 5 + rank; p++) array[p] = 11.1 + p * 0.1 + rank * 100.0;
     41:     PetscCall(DMSwarmRestoreField(dms, "viscosity", NULL, NULL, (void **)&array));
     42:   }

     44:   {
     45:     PetscReal *array;
     46:     PetscCall(DMSwarmGetField(dms, "strain", NULL, NULL, (void **)&array));
     47:     for (p = 0; p < 5 + rank; p++) array[p] = 2.0e-2 + p * 0.001 + rank * 1.0;
     48:     PetscCall(DMSwarmRestoreField(dms, "strain", NULL, NULL, (void **)&array));
     49:   }

     51:   PetscCall(DMSwarmCreateGlobalVectorFromField(dms, "viscosity", &x));
     52:   PetscCall(DMSwarmDestroyGlobalVectorFromField(dms, "viscosity", &x));

     54:   PetscCall(DMSwarmVectorDefineField(dms, "strain"));
     55:   PetscCall(DMCreateGlobalVector(dms, &x));
     56:   PetscCall(VecDestroy(&x));

     58:   {
     59:     PetscInt *rankval;
     60:     PetscInt  npoints[2], npoints_orig[2];

     62:     PetscCall(DMSwarmGetLocalSize(dms, &npoints_orig[0]));
     63:     PetscCall(DMSwarmGetSize(dms, &npoints_orig[1]));
     64:     PetscCall(DMSwarmGetField(dms, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
     65:     if ((rank == 0) && (size > 1)) {
     66:       rankval[0] = 1;
     67:       rankval[3] = 1;
     68:     }
     69:     if (rank == 3) rankval[2] = 1;
     70:     PetscCall(DMSwarmRestoreField(dms, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
     71:     PetscCall(DMSwarmMigrate(dms, PETSC_TRUE));
     72:     PetscCall(DMSwarmGetLocalSize(dms, &npoints[0]));
     73:     PetscCall(DMSwarmGetSize(dms, &npoints[1]));
     74:     PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
     75:     PetscCall(PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD, "rank[%d] before(%" PetscInt_FMT ",%" PetscInt_FMT ") after(%" PetscInt_FMT ",%" PetscInt_FMT ")\n", rank, npoints_orig[0], npoints_orig[1], npoints[0], npoints[1]));
     76:     PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
     77:     PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));
     78:   }
     79:   {
     80:     PetscCall(DMSwarmCreateGlobalVectorFromField(dms, "viscosity", &x));
     81:     PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
     82:     PetscCall(DMSwarmDestroyGlobalVectorFromField(dms, "viscosity", &x));
     83:   }
     84:   {
     85:     PetscCall(DMSwarmCreateGlobalVectorFromField(dms, "strain", &x));
     86:     PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
     87:     PetscCall(DMSwarmDestroyGlobalVectorFromField(dms, "strain", &x));
     88:   }

     90:   PetscCall(DMDestroy(&dms));
     91:   PetscFunctionReturn(PETSC_SUCCESS);
     92: }

     94: PetscErrorCode ex1_2(void)
     95: {
     96:   DM          dms;
     97:   Vec         x;
     98:   PetscMPIInt rank, size;
     99:   PetscInt    p;

    101:   PetscFunctionBegin;
    102:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
    103:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
    104:   PetscCall(DMCreate(PETSC_COMM_WORLD, &dms));
    105:   PetscCall(DMSetType(dms, DMSWARM));
    106:   PetscCall(PetscObjectSetName((PetscObject)dms, "Particles"));
    107:   PetscCall(DMSwarmInitializeFieldRegister(dms));

    109:   PetscCall(DMSwarmRegisterPetscDatatypeField(dms, "viscosity", 1, PETSC_REAL));
    110:   PetscCall(DMSwarmRegisterPetscDatatypeField(dms, "strain", 1, PETSC_REAL));
    111:   PetscCall(DMSwarmFinalizeFieldRegister(dms));
    112:   PetscCall(DMSwarmSetLocalSizes(dms, 5 + rank, 4));
    113:   PetscCall(DMView(dms, PETSC_VIEWER_STDOUT_WORLD));
    114:   {
    115:     PetscReal *array;
    116:     PetscCall(DMSwarmGetField(dms, "viscosity", NULL, NULL, (void **)&array));
    117:     for (p = 0; p < 5 + rank; p++) array[p] = 11.1 + p * 0.1 + rank * 100.0;
    118:     PetscCall(DMSwarmRestoreField(dms, "viscosity", NULL, NULL, (void **)&array));
    119:   }
    120:   {
    121:     PetscReal *array;
    122:     PetscCall(DMSwarmGetField(dms, "strain", NULL, NULL, (void **)&array));
    123:     for (p = 0; p < 5 + rank; p++) array[p] = 2.0e-2 + p * 0.001 + rank * 1.0;
    124:     PetscCall(DMSwarmRestoreField(dms, "strain", NULL, NULL, (void **)&array));
    125:   }
    126:   {
    127:     PetscInt *rankval;
    128:     PetscInt  npoints[2], npoints_orig[2];

    130:     PetscCall(DMSwarmGetLocalSize(dms, &npoints_orig[0]));
    131:     PetscCall(DMSwarmGetSize(dms, &npoints_orig[1]));
    132:     PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
    133:     PetscCall(PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD, "rank[%d] before(%" PetscInt_FMT ",%" PetscInt_FMT ")\n", rank, npoints_orig[0], npoints_orig[1]));
    134:     PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
    135:     PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));

    137:     PetscCall(DMSwarmGetField(dms, "DMSwarm_rank", NULL, NULL, (void **)&rankval));

    139:     if (rank == 1) rankval[0] = -1;
    140:     if (rank == 2) rankval[1] = -1;
    141:     if (rank == 3) {
    142:       rankval[3] = -1;
    143:       rankval[4] = -1;
    144:     }
    145:     PetscCall(DMSwarmRestoreField(dms, "DMSwarm_rank", NULL, NULL, (void **)&rankval));
    146:     PetscCall(DMSwarmCollectViewCreate(dms));
    147:     PetscCall(DMSwarmGetLocalSize(dms, &npoints[0]));
    148:     PetscCall(DMSwarmGetSize(dms, &npoints[1]));
    149:     PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
    150:     PetscCall(PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD, "rank[%d] after(%" PetscInt_FMT ",%" PetscInt_FMT ")\n", rank, npoints[0], npoints[1]));
    151:     PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
    152:     PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));

    154:     PetscCall(DMSwarmCreateGlobalVectorFromField(dms, "viscosity", &x));
    155:     PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
    156:     PetscCall(DMSwarmDestroyGlobalVectorFromField(dms, "viscosity", &x));

    158:     PetscCall(DMSwarmCollectViewDestroy(dms));
    159:     PetscCall(DMSwarmGetLocalSize(dms, &npoints[0]));
    160:     PetscCall(DMSwarmGetSize(dms, &npoints[1]));
    161:     PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
    162:     PetscCall(PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD, "rank[%d] after_v(%" PetscInt_FMT ",%" PetscInt_FMT ")\n", rank, npoints[0], npoints[1]));
    163:     PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
    164:     PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));

    166:     PetscCall(DMSwarmCreateGlobalVectorFromField(dms, "viscosity", &x));
    167:     PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
    168:     PetscCall(DMSwarmDestroyGlobalVectorFromField(dms, "viscosity", &x));
    169:   }
    170:   PetscCall(DMDestroy(&dms));
    171:   PetscFunctionReturn(PETSC_SUCCESS);
    172: }

    174: /*
    175:  splot "c-rank0.gp","c-rank1.gp","c-rank2.gp","c-rank3.gp"
    176: */
    177: PetscErrorCode ex1_3(void)
    178: {
    179:   DM          dms;
    180:   PetscMPIInt rank, size;
    181:   PetscInt    is, js, ni, nj, overlap;
    182:   DM          dmcell;

    184:   PetscFunctionBegin;
    185:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
    186:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
    187:   overlap = 2;
    188:   PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 13, 13, PETSC_DECIDE, PETSC_DECIDE, 1, overlap, NULL, NULL, &dmcell));
    189:   PetscCall(DMSetFromOptions(dmcell));
    190:   PetscCall(DMSetUp(dmcell));
    191:   PetscCall(DMDASetUniformCoordinates(dmcell, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
    192:   PetscCall(DMDAGetCorners(dmcell, &is, &js, NULL, &ni, &nj, NULL));
    193:   PetscCall(DMCreate(PETSC_COMM_WORLD, &dms));
    194:   PetscCall(DMSetType(dms, DMSWARM));
    195:   PetscCall(PetscObjectSetName((PetscObject)dms, "Particles"));
    196:   PetscCall(DMSwarmSetCellDM(dms, dmcell));

    198:   /* load in data types */
    199:   PetscCall(DMSwarmInitializeFieldRegister(dms));
    200:   PetscCall(DMSwarmRegisterPetscDatatypeField(dms, "viscosity", 1, PETSC_REAL));
    201:   PetscCall(DMSwarmRegisterPetscDatatypeField(dms, "coorx", 1, PETSC_REAL));
    202:   PetscCall(DMSwarmRegisterPetscDatatypeField(dms, "coory", 1, PETSC_REAL));
    203:   PetscCall(DMSwarmFinalizeFieldRegister(dms));
    204:   PetscCall(DMSwarmSetLocalSizes(dms, ni * nj * 4, 4));
    205:   PetscCall(DMView(dms, PETSC_VIEWER_STDOUT_WORLD));

    207:   /* set values within the swarm */
    208:   {
    209:     PetscReal   *array_x, *array_y;
    210:     PetscInt     npoints, i, j, cnt;
    211:     DMDACoor2d **LA_coor;
    212:     Vec          coor;
    213:     DM           dmcellcdm;

    215:     PetscCall(DMGetCoordinateDM(dmcell, &dmcellcdm));
    216:     PetscCall(DMGetCoordinates(dmcell, &coor));
    217:     PetscCall(DMDAVecGetArray(dmcellcdm, coor, &LA_coor));
    218:     PetscCall(DMSwarmGetLocalSize(dms, &npoints));
    219:     PetscCall(DMSwarmGetField(dms, "coorx", NULL, NULL, (void **)&array_x));
    220:     PetscCall(DMSwarmGetField(dms, "coory", NULL, NULL, (void **)&array_y));
    221:     cnt = 0;
    222:     for (j = js; j < js + nj; j++) {
    223:       for (i = is; i < is + ni; i++) {
    224:         PetscReal xp, yp;
    225:         xp                   = PetscRealPart(LA_coor[j][i].x);
    226:         yp                   = PetscRealPart(LA_coor[j][i].y);
    227:         array_x[4 * cnt + 0] = xp - 0.05;
    228:         if (array_x[4 * cnt + 0] < -1.0) array_x[4 * cnt + 0] = -1.0 + 1.0e-12;
    229:         array_x[4 * cnt + 1] = xp + 0.05;
    230:         if (array_x[4 * cnt + 1] > 1.0) array_x[4 * cnt + 1] = 1.0 - 1.0e-12;
    231:         array_x[4 * cnt + 2] = xp - 0.05;
    232:         if (array_x[4 * cnt + 2] < -1.0) array_x[4 * cnt + 2] = -1.0 + 1.0e-12;
    233:         array_x[4 * cnt + 3] = xp + 0.05;
    234:         if (array_x[4 * cnt + 3] > 1.0) array_x[4 * cnt + 3] = 1.0 - 1.0e-12;

    236:         array_y[4 * cnt + 0] = yp - 0.05;
    237:         if (array_y[4 * cnt + 0] < -1.0) array_y[4 * cnt + 0] = -1.0 + 1.0e-12;
    238:         array_y[4 * cnt + 1] = yp - 0.05;
    239:         if (array_y[4 * cnt + 1] < -1.0) array_y[4 * cnt + 1] = -1.0 + 1.0e-12;
    240:         array_y[4 * cnt + 2] = yp + 0.05;
    241:         if (array_y[4 * cnt + 2] > 1.0) array_y[4 * cnt + 2] = 1.0 - 1.0e-12;
    242:         array_y[4 * cnt + 3] = yp + 0.05;
    243:         if (array_y[4 * cnt + 3] > 1.0) array_y[4 * cnt + 3] = 1.0 - 1.0e-12;
    244:         cnt++;
    245:       }
    246:     }
    247:     PetscCall(DMSwarmRestoreField(dms, "coory", NULL, NULL, (void **)&array_y));
    248:     PetscCall(DMSwarmRestoreField(dms, "coorx", NULL, NULL, (void **)&array_x));
    249:     PetscCall(DMDAVecRestoreArray(dmcellcdm, coor, &LA_coor));
    250:   }
    251:   {
    252:     PetscInt npoints[2], npoints_orig[2], ng;

    254:     PetscCall(DMSwarmGetLocalSize(dms, &npoints_orig[0]));
    255:     PetscCall(DMSwarmGetSize(dms, &npoints_orig[1]));
    256:     PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
    257:     PetscCall(PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD, "rank[%d] before(%" PetscInt_FMT ",%" PetscInt_FMT ")\n", rank, npoints_orig[0], npoints_orig[1]));
    258:     PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
    259:     PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));
    260:     PetscCall(DMSwarmCollect_DMDABoundingBox(dms, &ng));

    262:     PetscCall(DMSwarmGetLocalSize(dms, &npoints[0]));
    263:     PetscCall(DMSwarmGetSize(dms, &npoints[1]));
    264:     PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
    265:     PetscCall(PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD, "rank[%d] after(%" PetscInt_FMT ",%" PetscInt_FMT ")\n", rank, npoints[0], npoints[1]));
    266:     PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
    267:     PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));
    268:   }
    269:   {
    270:     PetscReal *array_x, *array_y;
    271:     PetscInt   npoints, p;
    272:     FILE      *fp = NULL;
    273:     char       name[PETSC_MAX_PATH_LEN];

    275:     PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "c-rank%d.gp", rank));
    276:     fp = fopen(name, "w");
    277:     PetscCheck(fp, PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open file %s", name);
    278:     PetscCall(DMSwarmGetLocalSize(dms, &npoints));
    279:     PetscCall(DMSwarmGetField(dms, "coorx", NULL, NULL, (void **)&array_x));
    280:     PetscCall(DMSwarmGetField(dms, "coory", NULL, NULL, (void **)&array_y));
    281:     for (p = 0; p < npoints; p++) fprintf(fp, "%+1.4e %+1.4e %1.4e\n", array_x[p], array_y[p], (double)rank);
    282:     PetscCall(DMSwarmRestoreField(dms, "coory", NULL, NULL, (void **)&array_y));
    283:     PetscCall(DMSwarmRestoreField(dms, "coorx", NULL, NULL, (void **)&array_x));
    284:     fclose(fp);
    285:   }
    286:   PetscCall(DMDestroy(&dmcell));
    287:   PetscCall(DMDestroy(&dms));
    288:   PetscFunctionReturn(PETSC_SUCCESS);
    289: }

    291: typedef struct {
    292:   PetscReal cx[2];
    293:   PetscReal radius;
    294: } CollectZoneCtx;

    296: PetscErrorCode collect_zone(DM dm, PetscCtx ctx, PetscInt *nfound, PetscInt **foundlist)
    297: {
    298:   CollectZoneCtx *zone = (CollectZoneCtx *)ctx;
    299:   PetscInt        p, npoints;
    300:   PetscReal      *array_x, *array_y, r2;
    301:   PetscInt        p2collect, *plist;
    302:   PetscMPIInt     rank;

    304:   PetscFunctionBegin;
    305:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
    306:   PetscCall(DMSwarmGetLocalSize(dm, &npoints));
    307:   PetscCall(DMSwarmGetField(dm, "coorx", NULL, NULL, (void **)&array_x));
    308:   PetscCall(DMSwarmGetField(dm, "coory", NULL, NULL, (void **)&array_y));

    310:   r2        = zone->radius * zone->radius;
    311:   p2collect = 0;
    312:   for (p = 0; p < npoints; p++) {
    313:     PetscReal sep2;

    315:     sep2 = (array_x[p] - zone->cx[0]) * (array_x[p] - zone->cx[0]);
    316:     sep2 += (array_y[p] - zone->cx[1]) * (array_y[p] - zone->cx[1]);
    317:     if (sep2 < r2) p2collect++;
    318:   }

    320:   PetscCall(PetscMalloc1(p2collect + 1, &plist));
    321:   p2collect = 0;
    322:   for (p = 0; p < npoints; p++) {
    323:     PetscReal sep2;

    325:     sep2 = (array_x[p] - zone->cx[0]) * (array_x[p] - zone->cx[0]);
    326:     sep2 += (array_y[p] - zone->cx[1]) * (array_y[p] - zone->cx[1]);
    327:     if (sep2 < r2) {
    328:       plist[p2collect] = p;
    329:       p2collect++;
    330:     }
    331:   }
    332:   PetscCall(DMSwarmRestoreField(dm, "coory", NULL, NULL, (void **)&array_y));
    333:   PetscCall(DMSwarmRestoreField(dm, "coorx", NULL, NULL, (void **)&array_x));

    335:   *nfound    = p2collect;
    336:   *foundlist = plist;
    337:   PetscFunctionReturn(PETSC_SUCCESS);
    338: }

    340: PetscErrorCode ex1_4(void)
    341: {
    342:   DM              dms;
    343:   PetscMPIInt     rank, size;
    344:   PetscInt        is, js, ni, nj, overlap, nn;
    345:   DM              dmcell;
    346:   CollectZoneCtx *zone;
    347:   PetscReal       dx;

    349:   PetscFunctionBegin;
    350:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
    351:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
    352:   nn      = 101;
    353:   dx      = 2.0 / (PetscReal)(nn - 1);
    354:   overlap = 0;
    355:   PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, nn, nn, PETSC_DECIDE, PETSC_DECIDE, 1, overlap, NULL, NULL, &dmcell));
    356:   PetscCall(DMSetFromOptions(dmcell));
    357:   PetscCall(DMSetUp(dmcell));
    358:   PetscCall(DMDASetUniformCoordinates(dmcell, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
    359:   PetscCall(DMDAGetCorners(dmcell, &is, &js, NULL, &ni, &nj, NULL));
    360:   PetscCall(DMCreate(PETSC_COMM_WORLD, &dms));
    361:   PetscCall(DMSetType(dms, DMSWARM));
    362:   PetscCall(PetscObjectSetName((PetscObject)dms, "Particles"));

    364:   /* load in data types */
    365:   PetscCall(DMSwarmInitializeFieldRegister(dms));
    366:   PetscCall(DMSwarmRegisterPetscDatatypeField(dms, "viscosity", 1, PETSC_REAL));
    367:   PetscCall(DMSwarmRegisterPetscDatatypeField(dms, "coorx", 1, PETSC_REAL));
    368:   PetscCall(DMSwarmRegisterPetscDatatypeField(dms, "coory", 1, PETSC_REAL));
    369:   PetscCall(DMSwarmFinalizeFieldRegister(dms));
    370:   PetscCall(DMSwarmSetLocalSizes(dms, ni * nj * 4, 4));
    371:   PetscCall(DMView(dms, PETSC_VIEWER_STDOUT_WORLD));

    373:   /* set values within the swarm */
    374:   {
    375:     PetscReal   *array_x, *array_y;
    376:     PetscInt     npoints, i, j, cnt;
    377:     DMDACoor2d **LA_coor;
    378:     Vec          coor;
    379:     DM           dmcellcdm;

    381:     PetscCall(DMGetCoordinateDM(dmcell, &dmcellcdm));
    382:     PetscCall(DMGetCoordinates(dmcell, &coor));
    383:     PetscCall(DMDAVecGetArray(dmcellcdm, coor, &LA_coor));
    384:     PetscCall(DMSwarmGetLocalSize(dms, &npoints));
    385:     PetscCall(DMSwarmGetField(dms, "coorx", NULL, NULL, (void **)&array_x));
    386:     PetscCall(DMSwarmGetField(dms, "coory", NULL, NULL, (void **)&array_y));
    387:     cnt = 0;
    388:     for (j = js; j < js + nj; j++) {
    389:       for (i = is; i < is + ni; i++) {
    390:         PetscReal xp, yp;

    392:         xp                   = PetscRealPart(LA_coor[j][i].x);
    393:         yp                   = PetscRealPart(LA_coor[j][i].y);
    394:         array_x[4 * cnt + 0] = xp - dx * 0.1; /*if (array_x[4*cnt+0] < -1.0) array_x[4*cnt+0] = -1.0+1.0e-12;*/
    395:         array_x[4 * cnt + 1] = xp + dx * 0.1; /*if (array_x[4*cnt+1] > 1.0)  { array_x[4*cnt+1] =  1.0-1.0e-12; }*/
    396:         array_x[4 * cnt + 2] = xp - dx * 0.1; /*if (array_x[4*cnt+2] < -1.0) array_x[4*cnt+2] = -1.0+1.0e-12;*/
    397:         array_x[4 * cnt + 3] = xp + dx * 0.1; /*if (array_x[4*cnt+3] > 1.0)  { array_x[4*cnt+3] =  1.0-1.0e-12; }*/
    398:         array_y[4 * cnt + 0] = yp - dx * 0.1; /*if (array_y[4*cnt+0] < -1.0) array_y[4*cnt+0] = -1.0+1.0e-12;*/
    399:         array_y[4 * cnt + 1] = yp - dx * 0.1; /*if (array_y[4*cnt+1] < -1.0) array_y[4*cnt+1] = -1.0+1.0e-12;*/
    400:         array_y[4 * cnt + 2] = yp + dx * 0.1; /*if (array_y[4*cnt+2] > 1.0)  { array_y[4*cnt+2] =  1.0-1.0e-12; }*/
    401:         array_y[4 * cnt + 3] = yp + dx * 0.1; /*if (array_y[4*cnt+3] > 1.0)  { array_y[4*cnt+3] =  1.0-1.0e-12; }*/
    402:         cnt++;
    403:       }
    404:     }
    405:     PetscCall(DMSwarmRestoreField(dms, "coory", NULL, NULL, (void **)&array_y));
    406:     PetscCall(DMSwarmRestoreField(dms, "coorx", NULL, NULL, (void **)&array_x));
    407:     PetscCall(DMDAVecRestoreArray(dmcellcdm, coor, &LA_coor));
    408:   }
    409:   PetscCall(PetscMalloc1(1, &zone));
    410:   if (size == 4) {
    411:     if (rank == 0) {
    412:       zone->cx[0]  = 0.5;
    413:       zone->cx[1]  = 0.5;
    414:       zone->radius = 0.3;
    415:     }
    416:     if (rank == 1) {
    417:       zone->cx[0]  = -0.5;
    418:       zone->cx[1]  = 0.5;
    419:       zone->radius = 0.25;
    420:     }
    421:     if (rank == 2) {
    422:       zone->cx[0]  = 0.5;
    423:       zone->cx[1]  = -0.5;
    424:       zone->radius = 0.2;
    425:     }
    426:     if (rank == 3) {
    427:       zone->cx[0]  = -0.5;
    428:       zone->cx[1]  = -0.5;
    429:       zone->radius = 0.1;
    430:     }
    431:   } else {
    432:     if (rank == 0) {
    433:       zone->cx[0]  = 0.5;
    434:       zone->cx[1]  = 0.5;
    435:       zone->radius = 0.8;
    436:     } else {
    437:       zone->cx[0]  = 10.0;
    438:       zone->cx[1]  = 10.0;
    439:       zone->radius = 0.0;
    440:     }
    441:   }
    442:   {
    443:     PetscInt npoints[2], npoints_orig[2], ng;

    445:     PetscCall(DMSwarmGetLocalSize(dms, &npoints_orig[0]));
    446:     PetscCall(DMSwarmGetSize(dms, &npoints_orig[1]));
    447:     PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
    448:     PetscCall(PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD, "rank[%d] before(%" PetscInt_FMT ",%" PetscInt_FMT ")\n", rank, npoints_orig[0], npoints_orig[1]));
    449:     PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
    450:     PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));
    451:     PetscCall(DMSwarmCollect_General(dms, collect_zone, sizeof(CollectZoneCtx), zone, &ng));
    452:     PetscCall(DMSwarmGetLocalSize(dms, &npoints[0]));
    453:     PetscCall(DMSwarmGetSize(dms, &npoints[1]));
    454:     PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
    455:     PetscCall(PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD, "rank[%d] after(%" PetscInt_FMT ",%" PetscInt_FMT ")\n", rank, npoints[0], npoints[1]));
    456:     PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
    457:     PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));
    458:   }
    459:   {
    460:     PetscReal *array_x, *array_y;
    461:     PetscInt   npoints, p;
    462:     FILE      *fp = NULL;
    463:     char       name[PETSC_MAX_PATH_LEN];

    465:     PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "c-rank%d.gp", rank));
    466:     fp = fopen(name, "w");
    467:     PetscCheck(fp, PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open file %s", name);
    468:     PetscCall(DMSwarmGetLocalSize(dms, &npoints));
    469:     PetscCall(DMSwarmGetField(dms, "coorx", NULL, NULL, (void **)&array_x));
    470:     PetscCall(DMSwarmGetField(dms, "coory", NULL, NULL, (void **)&array_y));
    471:     for (p = 0; p < npoints; p++) fprintf(fp, "%+1.4e %+1.4e %1.4e\n", array_x[p], array_y[p], (double)rank);
    472:     PetscCall(DMSwarmRestoreField(dms, "coory", NULL, NULL, (void **)&array_y));
    473:     PetscCall(DMSwarmRestoreField(dms, "coorx", NULL, NULL, (void **)&array_x));
    474:     fclose(fp);
    475:   }
    476:   PetscCall(DMDestroy(&dmcell));
    477:   PetscCall(DMDestroy(&dms));
    478:   PetscCall(PetscFree(zone));
    479:   PetscFunctionReturn(PETSC_SUCCESS);
    480: }

    482: int main(int argc, char **argv)
    483: {
    484:   PetscInt test_mode = 4;

    486:   PetscFunctionBeginUser;
    487:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    488:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-test_mode", &test_mode, NULL));
    489:   if (test_mode == 1) PetscCall(ex1_1());
    490:   else if (test_mode == 2) PetscCall(ex1_2());
    491:   else if (test_mode == 3) PetscCall(ex1_3());
    492:   else {
    493:     PetscCheck(test_mode == 4, PETSC_COMM_SELF, PETSC_ERR_USER, "Unknown test_mode value, should be 1,2,3,4");
    494:     PetscCall(ex1_4());
    495:   }
    496:   PetscCall(PetscFinalize());
    497:   return 0;
    498: }

    500: /*TEST

    502:    build:
    503:       requires: !complex double

    505:    test:
    506:       args: -test_mode 1
    507:       filter: grep -v atomic
    508:       filter_output: grep -v atomic

    510:    test:
    511:       suffix: 2
    512:       args: -test_mode 2
    513:       filter: grep -v atomic
    514:       filter_output: grep -v atomic

    516:    test:
    517:       suffix: 3
    518:       args: -test_mode 3
    519:       filter: grep -v atomic
    520:       filter_output: grep -v atomic
    521:       TODO: broken

    523:    test:
    524:       suffix: 4
    525:       args: -test_mode 4
    526:       filter: grep -v atomic
    527:       filter_output: grep -v atomic

    529:    test:
    530:       suffix: 5
    531:       nsize: 4
    532:       args: -test_mode 1
    533:       filter: grep -v atomic
    534:       filter_output: grep -v atomic

    536:    test:
    537:       suffix: 6
    538:       nsize: 4
    539:       args: -test_mode 2
    540:       filter: grep -v atomic
    541:       filter_output: grep -v atomic

    543:    test:
    544:       suffix: 7
    545:       nsize: 4
    546:       args: -test_mode 3
    547:       filter: grep -v atomic
    548:       filter_output: grep -v atomic
    549:       TODO: broken

    551:    test:
    552:       suffix: 8
    553:       nsize: 4
    554:       args: -test_mode 4
    555:       filter: grep -v atomic
    556:       filter_output: grep -v atomic

    558: TEST*/