Actual source code: dmproject.c

      1: #include <petsc/private/dmimpl.h>
      2: #include <petscdm.h>
      3: #include <petscdmda.h>
      4: #include <petscdmplex.h>
      5: #include <petscdmswarm.h>
      6: #include <petscksp.h>
      7: #include <petscblaslapack.h>

      9: #include <petsc/private/dmswarmimpl.h>
     10: #include "../src/dm/impls/swarm/data_bucket.h" // For DataBucket internals

     12: typedef struct _projectConstraintsCtx {
     13:   DM  dm;
     14:   Vec mask;
     15: } projectConstraintsCtx;

     17: static PetscErrorCode MatMult_GlobalToLocalNormal(Mat CtC, Vec x, Vec y)
     18: {
     19:   DM                     dm;
     20:   Vec                    local, mask;
     21:   projectConstraintsCtx *ctx;

     23:   PetscFunctionBegin;
     24:   PetscCall(MatShellGetContext(CtC, &ctx));
     25:   dm   = ctx->dm;
     26:   mask = ctx->mask;
     27:   PetscCall(DMGetLocalVector(dm, &local));
     28:   PetscCall(DMGlobalToLocalBegin(dm, x, INSERT_VALUES, local));
     29:   PetscCall(DMGlobalToLocalEnd(dm, x, INSERT_VALUES, local));
     30:   if (mask) PetscCall(VecPointwiseMult(local, mask, local));
     31:   PetscCall(VecSet(y, 0.));
     32:   PetscCall(DMLocalToGlobalBegin(dm, local, ADD_VALUES, y));
     33:   PetscCall(DMLocalToGlobalEnd(dm, local, ADD_VALUES, y));
     34:   PetscCall(DMRestoreLocalVector(dm, &local));
     35:   PetscFunctionReturn(PETSC_SUCCESS);
     36: }

     38: static PetscErrorCode DMGlobalToLocalSolve_project1(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx)
     39: {
     40:   PetscInt f;

     42:   PetscFunctionBegin;
     43:   for (f = 0; f < Nf; f++) u[f] = 1.;
     44:   PetscFunctionReturn(PETSC_SUCCESS);
     45: }

     47: /*@
     48:   DMGlobalToLocalSolve - Solve for the global vector that is mapped to a given local vector by `DMGlobalToLocalBegin()`/`DMGlobalToLocalEnd()` with mode
     49:   `INSERT_VALUES`.

     51:   Collective

     53:   Input Parameters:
     54: + dm - The `DM` object
     55: . x  - The local vector
     56: - y  - The global vector: the input value of globalVec is used as an initial guess

     58:   Output Parameter:
     59: . y - The least-squares solution

     61:   Level: advanced

     63:   Note:
     64:   It is assumed that the sum of all the local vector sizes is greater than or equal to the global vector size, so the solution is
     65:   a least-squares solution.  It is also assumed that `DMLocalToGlobalBegin()`/`DMLocalToGlobalEnd()` with mode `ADD_VALUES` is the adjoint of the
     66:   global-to-local map, so that the least-squares solution may be found by the normal equations.

     68:   If the `DM` is of type `DMPLEX`, then `y` is the solution of $ L^T * D * L * y = L^T * D * x $, where $D$ is a diagonal mask that is 1 for every point in
     69:   the union of the closures of the local cells and 0 otherwise.  This difference is only relevant if there are anchor points that are not in the
     70:   closure of any local cell (see `DMPlexGetAnchors()`/`DMPlexSetAnchors()`).

     72:   What is L?

     74:   If this solves for a global vector from a local vector why is not called `DMLocalToGlobalSolve()`?

     76: .seealso: [](ch_dmbase), `DM`, `DMGlobalToLocalBegin()`, `DMGlobalToLocalEnd()`, `DMLocalToGlobalBegin()`, `DMLocalToGlobalEnd()`, `DMPlexGetAnchors()`, `DMPlexSetAnchors()`
     77: @*/
     78: PetscErrorCode DMGlobalToLocalSolve(DM dm, Vec x, Vec y)
     79: {
     80:   Mat                   CtC;
     81:   PetscInt              n, N, cStart, cEnd, c;
     82:   PetscBool             isPlex;
     83:   KSP                   ksp;
     84:   PC                    pc;
     85:   Vec                   global, mask = NULL;
     86:   projectConstraintsCtx ctx;

     88:   PetscFunctionBegin;
     89:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isPlex));
     90:   if (isPlex) {
     91:     /* mark points in the closure */
     92:     PetscCall(DMCreateLocalVector(dm, &mask));
     93:     PetscCall(VecSet(mask, 0.0));
     94:     PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
     95:     if (cEnd > cStart) {
     96:       PetscScalar *ones;
     97:       PetscInt     numValues, i;

     99:       PetscCall(DMPlexVecGetClosure(dm, NULL, mask, cStart, &numValues, NULL));
    100:       PetscCall(PetscMalloc1(numValues, &ones));
    101:       for (i = 0; i < numValues; i++) ones[i] = 1.;
    102:       for (c = cStart; c < cEnd; c++) PetscCall(DMPlexVecSetClosure(dm, NULL, mask, c, ones, INSERT_VALUES));
    103:       PetscCall(PetscFree(ones));
    104:     }
    105:   } else {
    106:     PetscBool hasMask;

    108:     PetscCall(DMHasNamedLocalVector(dm, "_DMGlobalToLocalSolve_mask", &hasMask));
    109:     if (!hasMask) {
    110:       PetscErrorCode (**func)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
    111:       void   **ctx;
    112:       PetscInt Nf, f;

    114:       PetscCall(DMGetNumFields(dm, &Nf));
    115:       PetscCall(PetscMalloc2(Nf, &func, Nf, &ctx));
    116:       for (f = 0; f < Nf; ++f) {
    117:         func[f] = DMGlobalToLocalSolve_project1;
    118:         ctx[f]  = NULL;
    119:       }
    120:       PetscCall(DMGetNamedLocalVector(dm, "_DMGlobalToLocalSolve_mask", &mask));
    121:       PetscCall(DMProjectFunctionLocal(dm, 0.0, func, ctx, INSERT_ALL_VALUES, mask));
    122:       PetscCall(DMRestoreNamedLocalVector(dm, "_DMGlobalToLocalSolve_mask", &mask));
    123:       PetscCall(PetscFree2(func, ctx));
    124:     }
    125:     PetscCall(DMGetNamedLocalVector(dm, "_DMGlobalToLocalSolve_mask", &mask));
    126:   }
    127:   ctx.dm   = dm;
    128:   ctx.mask = mask;
    129:   PetscCall(VecGetSize(y, &N));
    130:   PetscCall(VecGetLocalSize(y, &n));
    131:   PetscCall(MatCreate(PetscObjectComm((PetscObject)dm), &CtC));
    132:   PetscCall(MatSetSizes(CtC, n, n, N, N));
    133:   PetscCall(MatSetType(CtC, MATSHELL));
    134:   PetscCall(MatSetUp(CtC));
    135:   PetscCall(MatShellSetContext(CtC, &ctx));
    136:   PetscCall(MatShellSetOperation(CtC, MATOP_MULT, (void (*)(void))MatMult_GlobalToLocalNormal));
    137:   PetscCall(KSPCreate(PetscObjectComm((PetscObject)dm), &ksp));
    138:   PetscCall(KSPSetOperators(ksp, CtC, CtC));
    139:   PetscCall(KSPSetType(ksp, KSPCG));
    140:   PetscCall(KSPGetPC(ksp, &pc));
    141:   PetscCall(PCSetType(pc, PCNONE));
    142:   PetscCall(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
    143:   PetscCall(KSPSetUp(ksp));
    144:   PetscCall(DMGetGlobalVector(dm, &global));
    145:   PetscCall(VecSet(global, 0.));
    146:   if (mask) PetscCall(VecPointwiseMult(x, mask, x));
    147:   PetscCall(DMLocalToGlobalBegin(dm, x, ADD_VALUES, global));
    148:   PetscCall(DMLocalToGlobalEnd(dm, x, ADD_VALUES, global));
    149:   PetscCall(KSPSolve(ksp, global, y));
    150:   PetscCall(DMRestoreGlobalVector(dm, &global));
    151:   /* clean up */
    152:   PetscCall(KSPDestroy(&ksp));
    153:   PetscCall(MatDestroy(&CtC));
    154:   if (isPlex) {
    155:     PetscCall(VecDestroy(&mask));
    156:   } else {
    157:     PetscCall(DMRestoreNamedLocalVector(dm, "_DMGlobalToLocalSolve_mask", &mask));
    158:   }
    159:   PetscFunctionReturn(PETSC_SUCCESS);
    160: }

    162: /*@C
    163:   DMProjectField - This projects a given function of the input fields into the function space provided by a `DM`, putting the coefficients in a global vector.

    165:   Collective

    167:   Input Parameters:
    168: + dm    - The `DM`
    169: . time  - The time
    170: . U     - The input field vector
    171: . funcs - The functions to evaluate, one per field
    172: - mode  - The insertion mode for values

    174:   Output Parameter:
    175: . X - The output vector

    177:   Calling sequence of `funcs`:
    178: + dim          - The spatial dimension
    179: . Nf           - The number of input fields
    180: . NfAux        - The number of input auxiliary fields
    181: . uOff         - The offset of each field in `u`
    182: . uOff_x       - The offset of each field in `u_x`
    183: . u            - The field values at this point in space
    184: . u_t          - The field time derivative at this point in space (or `NULL`)
    185: . u_x          - The field derivatives at this point in space
    186: . aOff         - The offset of each auxiliary field in `u`
    187: . aOff_x       - The offset of each auxiliary field in `u_x`
    188: . a            - The auxiliary field values at this point in space
    189: . a_t          - The auxiliary field time derivative at this point in space (or `NULL`)
    190: . a_x          - The auxiliary field derivatives at this point in space
    191: . t            - The current time
    192: . x            - The coordinates of this point
    193: . numConstants - The number of constants
    194: . constants    - The value of each constant
    195: - f            - The value of the function at this point in space

    197:   Level: advanced

    199:   Note:
    200:   There are three different `DM`s that potentially interact in this function. The output `dm`, specifies the layout of the values calculates by the function.
    201:   The input `DM`, attached to `U`, may be different. For example, you can input the solution over the full domain, but output over a piece of the boundary, or
    202:   a subdomain. You can also output a different number of fields than the input, with different discretizations. Last the auxiliary `DM`, attached to the
    203:   auxiliary field vector, which is attached to `dm`, can also be different. It can have a different topology, number of fields, and discretizations.

    205: .seealso: [](ch_dmbase), `DM`, `DMProjectFieldLocal()`, `DMProjectFieldLabelLocal()`, `DMProjectFunction()`, `DMComputeL2Diff()`
    206: @*/
    207: PetscErrorCode DMProjectField(DM dm, PetscReal time, Vec U, void (**funcs)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[]), InsertMode mode, Vec X)
    208: {
    209:   Vec localX, localU;
    210:   DM  dmIn;

    212:   PetscFunctionBegin;
    214:   PetscCall(DMGetLocalVector(dm, &localX));
    215:   /* We currently check whether locU == locX to see if we need to apply BC */
    216:   if (U != X) {
    217:     PetscCall(VecGetDM(U, &dmIn));
    218:     PetscCall(DMGetLocalVector(dmIn, &localU));
    219:   } else {
    220:     dmIn   = dm;
    221:     localU = localX;
    222:   }
    223:   PetscCall(DMGlobalToLocalBegin(dmIn, U, INSERT_VALUES, localU));
    224:   PetscCall(DMGlobalToLocalEnd(dmIn, U, INSERT_VALUES, localU));
    225:   PetscCall(DMProjectFieldLocal(dm, time, localU, funcs, mode, localX));
    226:   PetscCall(DMLocalToGlobalBegin(dm, localX, mode, X));
    227:   PetscCall(DMLocalToGlobalEnd(dm, localX, mode, X));
    228:   if (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES) {
    229:     Mat cMat;

    231:     PetscCall(DMGetDefaultConstraints(dm, NULL, &cMat, NULL));
    232:     if (cMat) PetscCall(DMGlobalToLocalSolve(dm, localX, X));
    233:   }
    234:   PetscCall(DMRestoreLocalVector(dm, &localX));
    235:   if (U != X) PetscCall(DMRestoreLocalVector(dmIn, &localU));
    236:   PetscFunctionReturn(PETSC_SUCCESS);
    237: }

    239: /********************* Adaptive Interpolation **************************/

    241: /* See the discussion of Adaptive Interpolation in manual/high_level_mg.rst */
    242: PetscErrorCode DMAdaptInterpolator(DM dmc, DM dmf, Mat In, KSP smoother, Mat MF, Mat MC, Mat *InAdapt, void *user)
    243: {
    244:   Mat                globalA, AF;
    245:   Vec                tmp;
    246:   const PetscScalar *af, *ac;
    247:   PetscScalar       *A, *b, *x, *workscalar;
    248:   PetscReal         *w, *sing, *workreal, rcond = PETSC_SMALL;
    249:   PetscBLASInt       M, N, one = 1, irank, lwrk, info;
    250:   PetscInt           debug = 0, rStart, rEnd, r, maxcols = 0, k, Nc, ldac, ldaf;
    251:   PetscBool          allocVc = PETSC_FALSE;

    253:   PetscFunctionBegin;
    254:   PetscCall(PetscLogEventBegin(DM_AdaptInterpolator, dmc, dmf, 0, 0));
    255:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-dm_interpolator_adapt_debug", &debug, NULL));
    256:   PetscCall(MatGetSize(MF, NULL, &Nc));
    257:   PetscCall(MatDuplicate(In, MAT_SHARE_NONZERO_PATTERN, InAdapt));
    258:   PetscCall(MatGetOwnershipRange(In, &rStart, &rEnd));
    259: #if 0
    260:   PetscCall(MatGetMaxRowLen(In, &maxcols));
    261: #else
    262:   for (r = rStart; r < rEnd; ++r) {
    263:     PetscInt ncols;

    265:     PetscCall(MatGetRow(In, r, &ncols, NULL, NULL));
    266:     maxcols = PetscMax(maxcols, ncols);
    267:     PetscCall(MatRestoreRow(In, r, &ncols, NULL, NULL));
    268:   }
    269: #endif
    270:   if (Nc < maxcols) PetscCall(PetscPrintf(PETSC_COMM_SELF, "The number of input vectors %" PetscInt_FMT " < %" PetscInt_FMT " the maximum number of column entries\n", Nc, maxcols));
    271:   for (k = 0; k < Nc && debug; ++k) {
    272:     char        name[PETSC_MAX_PATH_LEN];
    273:     const char *prefix;
    274:     Vec         vc, vf;

    276:     PetscCall(PetscObjectGetOptionsPrefix((PetscObject)smoother, &prefix));

    278:     if (MC) {
    279:       PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "%sCoarse Vector %" PetscInt_FMT, prefix ? prefix : NULL, k));
    280:       PetscCall(MatDenseGetColumnVecRead(MC, k, &vc));
    281:       PetscCall(PetscObjectSetName((PetscObject)vc, name));
    282:       PetscCall(VecViewFromOptions(vc, NULL, "-dm_adapt_interp_view_coarse"));
    283:       PetscCall(MatDenseRestoreColumnVecRead(MC, k, &vc));
    284:     }
    285:     PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "%sFine Vector %" PetscInt_FMT, prefix ? prefix : NULL, k));
    286:     PetscCall(MatDenseGetColumnVecRead(MF, k, &vf));
    287:     PetscCall(PetscObjectSetName((PetscObject)vf, name));
    288:     PetscCall(VecViewFromOptions(vf, NULL, "-dm_adapt_interp_view_fine"));
    289:     PetscCall(MatDenseRestoreColumnVecRead(MF, k, &vf));
    290:   }
    291:   PetscCall(PetscBLASIntCast(3 * PetscMin(Nc, maxcols) + PetscMax(2 * PetscMin(Nc, maxcols), PetscMax(Nc, maxcols)), &lwrk));
    292:   PetscCall(PetscMalloc7(Nc * maxcols, &A, PetscMax(Nc, maxcols), &b, Nc, &w, maxcols, &x, maxcols, &sing, lwrk, &workscalar, 5 * PetscMin(Nc, maxcols), &workreal));
    293:   /* w_k = \frac{\HC{v_k} B_l v_k}{\HC{v_k} A_l v_k} or the inverse Rayleigh quotient, which we calculate using \frac{\HC{v_k} v_k}{\HC{v_k} B^{-1}_l A_l v_k} */
    294:   PetscCall(KSPGetOperators(smoother, &globalA, NULL));

    296:   PetscCall(MatMatMult(globalA, MF, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &AF));
    297:   for (k = 0; k < Nc; ++k) {
    298:     PetscScalar vnorm, vAnorm;
    299:     Vec         vf;

    301:     w[k] = 1.0;
    302:     PetscCall(MatDenseGetColumnVecRead(MF, k, &vf));
    303:     PetscCall(MatDenseGetColumnVecRead(AF, k, &tmp));
    304:     PetscCall(VecDot(vf, vf, &vnorm));
    305: #if 0
    306:     PetscCall(DMGetGlobalVector(dmf, &tmp2));
    307:     PetscCall(KSPSolve(smoother, tmp, tmp2));
    308:     PetscCall(VecDot(vf, tmp2, &vAnorm));
    309:     PetscCall(DMRestoreGlobalVector(dmf, &tmp2));
    310: #else
    311:     PetscCall(VecDot(vf, tmp, &vAnorm));
    312: #endif
    313:     w[k] = PetscRealPart(vnorm) / PetscRealPart(vAnorm);
    314:     PetscCall(MatDenseRestoreColumnVecRead(MF, k, &vf));
    315:     PetscCall(MatDenseRestoreColumnVecRead(AF, k, &tmp));
    316:   }
    317:   PetscCall(MatDestroy(&AF));
    318:   if (!MC) {
    319:     allocVc = PETSC_TRUE;
    320:     PetscCall(MatTransposeMatMult(In, MF, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &MC));
    321:   }
    322:   /* Solve a LS system for each fine row
    323:      MATT: Can we generalize to the case where Nc for the fine space
    324:      is different for Nc for the coarse? */
    325:   PetscCall(MatDenseGetArrayRead(MF, &af));
    326:   PetscCall(MatDenseGetLDA(MF, &ldaf));
    327:   PetscCall(MatDenseGetArrayRead(MC, &ac));
    328:   PetscCall(MatDenseGetLDA(MC, &ldac));
    329:   for (r = rStart; r < rEnd; ++r) {
    330:     PetscInt           ncols, c;
    331:     const PetscInt    *cols;
    332:     const PetscScalar *vals;

    334:     PetscCall(MatGetRow(In, r, &ncols, &cols, &vals));
    335:     for (k = 0; k < Nc; ++k) {
    336:       /* Need to fit lowest mode exactly */
    337:       const PetscReal wk = ((ncols == 1) && (k > 0)) ? 0.0 : PetscSqrtReal(w[k]);

    339:       /* b_k = \sqrt{w_k} f^{F,k}_r */
    340:       b[k] = wk * af[r - rStart + k * ldaf];
    341:       /* A_{kc} = \sqrt{w_k} f^{C,k}_c */
    342:       /* TODO Must pull out VecScatter from In, scatter in vc[k] values up front, and access them indirectly just as in MatMult() */
    343:       for (c = 0; c < ncols; ++c) {
    344:         /* This is element (k, c) of A */
    345:         A[c * Nc + k] = wk * ac[cols[c] - rStart + k * ldac];
    346:       }
    347:     }
    348:     PetscCall(PetscBLASIntCast(Nc, &M));
    349:     PetscCall(PetscBLASIntCast(ncols, &N));
    350:     if (debug) {
    351: #if defined(PETSC_USE_COMPLEX)
    352:       PetscScalar *tmp;
    353:       PetscInt     j;

    355:       PetscCall(DMGetWorkArray(dmc, Nc, MPIU_SCALAR, (void *)&tmp));
    356:       for (j = 0; j < Nc; ++j) tmp[j] = w[j];
    357:       PetscCall(DMPrintCellMatrix(r, "Interpolator Row LS weights", Nc, 1, tmp));
    358:       PetscCall(DMPrintCellMatrix(r, "Interpolator Row LS matrix", Nc, ncols, A));
    359:       for (j = 0; j < Nc; ++j) tmp[j] = b[j];
    360:       PetscCall(DMPrintCellMatrix(r, "Interpolator Row LS rhs", Nc, 1, tmp));
    361:       PetscCall(DMRestoreWorkArray(dmc, Nc, MPIU_SCALAR, (void *)&tmp));
    362: #else
    363:       PetscCall(DMPrintCellMatrix(r, "Interpolator Row LS weights", Nc, 1, w));
    364:       PetscCall(DMPrintCellMatrix(r, "Interpolator Row LS matrix", Nc, ncols, A));
    365:       PetscCall(DMPrintCellMatrix(r, "Interpolator Row LS rhs", Nc, 1, b));
    366: #endif
    367:     }
    368: #if defined(PETSC_USE_COMPLEX)
    369:     /* ZGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, RWORK, INFO) */
    370:     PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &one, A, &M, b, M > N ? &M : &N, sing, &rcond, &irank, workscalar, &lwrk, workreal, &info));
    371: #else
    372:     /* DGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, INFO) */
    373:     PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&M, &N, &one, A, &M, b, M > N ? &M : &N, sing, &rcond, &irank, workscalar, &lwrk, &info));
    374: #endif
    375:     PetscCheck(info >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to GELSS");
    376:     PetscCheck(info <= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "SVD failed to converge");
    377:     if (debug) {
    378:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "rank %" PetscBLASInt_FMT " rcond %g\n", irank, (double)rcond));
    379: #if defined(PETSC_USE_COMPLEX)
    380:       {
    381:         PetscScalar *tmp;
    382:         PetscInt     j;

    384:         PetscCall(DMGetWorkArray(dmc, Nc, MPIU_SCALAR, (void *)&tmp));
    385:         for (j = 0; j < PetscMin(Nc, ncols); ++j) tmp[j] = sing[j];
    386:         PetscCall(DMPrintCellMatrix(r, "Interpolator Row LS singular values", PetscMin(Nc, ncols), 1, tmp));
    387:         PetscCall(DMRestoreWorkArray(dmc, Nc, MPIU_SCALAR, (void *)&tmp));
    388:       }
    389: #else
    390:       PetscCall(DMPrintCellMatrix(r, "Interpolator Row LS singular values", PetscMin(Nc, ncols), 1, sing));
    391: #endif
    392:       PetscCall(DMPrintCellMatrix(r, "Interpolator Row LS old P", ncols, 1, vals));
    393:       PetscCall(DMPrintCellMatrix(r, "Interpolator Row LS sol", ncols, 1, b));
    394:     }
    395:     PetscCall(MatSetValues(*InAdapt, 1, &r, ncols, cols, b, INSERT_VALUES));
    396:     PetscCall(MatRestoreRow(In, r, &ncols, &cols, &vals));
    397:   }
    398:   PetscCall(MatDenseRestoreArrayRead(MF, &af));
    399:   PetscCall(MatDenseRestoreArrayRead(MC, &ac));
    400:   PetscCall(PetscFree7(A, b, w, x, sing, workscalar, workreal));
    401:   if (allocVc) PetscCall(MatDestroy(&MC));
    402:   PetscCall(MatAssemblyBegin(*InAdapt, MAT_FINAL_ASSEMBLY));
    403:   PetscCall(MatAssemblyEnd(*InAdapt, MAT_FINAL_ASSEMBLY));
    404:   PetscCall(PetscLogEventEnd(DM_AdaptInterpolator, dmc, dmf, 0, 0));
    405:   PetscFunctionReturn(PETSC_SUCCESS);
    406: }

    408: PetscErrorCode DMCheckInterpolator(DM dmf, Mat In, Mat MC, Mat MF, PetscReal tol)
    409: {
    410:   Vec       tmp;
    411:   PetscReal norminf, norm2, maxnorminf = 0.0, maxnorm2 = 0.0;
    412:   PetscInt  k, Nc;

    414:   PetscFunctionBegin;
    415:   PetscCall(DMGetGlobalVector(dmf, &tmp));
    416:   PetscCall(MatViewFromOptions(In, NULL, "-dm_interpolator_adapt_error"));
    417:   PetscCall(MatGetSize(MF, NULL, &Nc));
    418:   for (k = 0; k < Nc; ++k) {
    419:     Vec vc, vf;

    421:     PetscCall(MatDenseGetColumnVecRead(MC, k, &vc));
    422:     PetscCall(MatDenseGetColumnVecRead(MF, k, &vf));
    423:     PetscCall(MatMult(In, vc, tmp));
    424:     PetscCall(VecAXPY(tmp, -1.0, vf));
    425:     PetscCall(VecViewFromOptions(vc, NULL, "-dm_interpolator_adapt_error"));
    426:     PetscCall(VecViewFromOptions(vf, NULL, "-dm_interpolator_adapt_error"));
    427:     PetscCall(VecViewFromOptions(tmp, NULL, "-dm_interpolator_adapt_error"));
    428:     PetscCall(VecNorm(tmp, NORM_INFINITY, &norminf));
    429:     PetscCall(VecNorm(tmp, NORM_2, &norm2));
    430:     maxnorminf = PetscMax(maxnorminf, norminf);
    431:     maxnorm2   = PetscMax(maxnorm2, norm2);
    432:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dmf), "Coarse vec %" PetscInt_FMT " ||vf - P vc||_\\infty %g, ||vf - P vc||_2 %g\n", k, (double)norminf, (double)norm2));
    433:     PetscCall(MatDenseRestoreColumnVecRead(MC, k, &vc));
    434:     PetscCall(MatDenseRestoreColumnVecRead(MF, k, &vf));
    435:   }
    436:   PetscCall(DMRestoreGlobalVector(dmf, &tmp));
    437:   PetscCheck(maxnorm2 <= tol, PetscObjectComm((PetscObject)dmf), PETSC_ERR_ARG_WRONG, "max_k ||vf_k - P vc_k||_2 %g > tol %g", (double)maxnorm2, (double)tol);
    438:   PetscFunctionReturn(PETSC_SUCCESS);
    439: }

    441: // Project particles to field
    442: //   M_f u_f = M_p u_p
    443: //   u_f = M^{-1}_f M_p u_p
    444: static PetscErrorCode DMSwarmProjectField_Conservative_PLEX(DM sw, DM dm, Vec u_p, Vec u_f)
    445: {
    446:   KSP         ksp;
    447:   Mat         M_f, M_p; // TODO Should cache these
    448:   Vec         rhs;
    449:   const char *prefix;

    451:   PetscFunctionBegin;
    452:   PetscCall(DMCreateMassMatrix(dm, dm, &M_f));
    453:   PetscCall(DMCreateMassMatrix(sw, dm, &M_p));
    454:   PetscCall(DMGetGlobalVector(dm, &rhs));
    455:   PetscCall(MatMultTranspose(M_p, u_p, rhs));

    457:   PetscCall(KSPCreate(PetscObjectComm((PetscObject)sw), &ksp));
    458:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)sw, &prefix));
    459:   PetscCall(KSPSetOptionsPrefix(ksp, prefix));
    460:   PetscCall(KSPAppendOptionsPrefix(ksp, "ptof_"));
    461:   PetscCall(KSPSetFromOptions(ksp));

    463:   PetscCall(KSPSetOperators(ksp, M_f, M_f));
    464:   PetscCall(KSPSolve(ksp, rhs, u_f));

    466:   PetscCall(DMRestoreGlobalVector(dm, &rhs));
    467:   PetscCall(KSPDestroy(&ksp));
    468:   PetscCall(MatDestroy(&M_f));
    469:   PetscCall(MatDestroy(&M_p));
    470:   PetscFunctionReturn(PETSC_SUCCESS);
    471: }

    473: // Project field to particles
    474: //   M_p u_p = M_f u_f
    475: //   u_p = M^+_p M_f u_f
    476: static PetscErrorCode DMSwarmProjectParticles_Conservative_PLEX(DM sw, DM dm, Vec u_p, Vec u_f)
    477: {
    478:   KSP         ksp;
    479:   PC          pc;
    480:   Mat         M_f, M_p, PM_p;
    481:   Vec         rhs;
    482:   PetscBool   isBjacobi;
    483:   const char *prefix;

    485:   PetscFunctionBegin;
    486:   PetscCall(DMCreateMassMatrix(dm, dm, &M_f));
    487:   PetscCall(DMCreateMassMatrix(sw, dm, &M_p));
    488:   PetscCall(DMGetGlobalVector(dm, &rhs));
    489:   PetscCall(MatMultTranspose(M_f, u_f, rhs));

    491:   PetscCall(KSPCreate(PetscObjectComm((PetscObject)sw), &ksp));
    492:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)sw, &prefix));
    493:   PetscCall(KSPSetOptionsPrefix(ksp, prefix));
    494:   PetscCall(KSPAppendOptionsPrefix(ksp, "ftop_"));
    495:   PetscCall(KSPSetFromOptions(ksp));

    497:   PetscCall(KSPGetPC(ksp, &pc));
    498:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCBJACOBI, &isBjacobi));
    499:   if (isBjacobi) {
    500:     PetscCall(DMSwarmCreateMassMatrixSquare(sw, dm, &PM_p));
    501:   } else {
    502:     PM_p = M_p;
    503:     PetscCall(PetscObjectReference((PetscObject)PM_p));
    504:   }
    505:   PetscCall(KSPSetOperators(ksp, M_p, PM_p));
    506:   PetscCall(KSPSolveTranspose(ksp, rhs, u_p));

    508:   PetscCall(DMRestoreGlobalVector(dm, &rhs));
    509:   PetscCall(KSPDestroy(&ksp));
    510:   PetscCall(MatDestroy(&M_f));
    511:   PetscCall(MatDestroy(&M_p));
    512:   PetscCall(MatDestroy(&PM_p));
    513:   PetscFunctionReturn(PETSC_SUCCESS);
    514: }

    516: static PetscErrorCode DMSwarmProjectFields_Plex_Internal(DM sw, DM dm, PetscInt Nf, const char *fieldnames[], Vec vec, ScatterMode mode)
    517: {
    518:   PetscDS  ds;
    519:   Vec      u;
    520:   PetscInt f = 0, bs, *Nc;

    522:   PetscFunctionBegin;
    523:   PetscCall(DMGetDS(dm, &ds));
    524:   PetscCall(PetscDSGetComponents(ds, &Nc));
    525:   PetscCall(PetscCitationsRegister(SwarmProjCitation, &SwarmProjcite));
    526:   PetscCheck(Nf == 1, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "Currently supported only for a single field");
    527:   PetscCall(DMSwarmVectorDefineField(sw, fieldnames[f]));
    528:   PetscCall(DMSwarmCreateGlobalVectorFromField(sw, fieldnames[f], &u));
    529:   PetscCall(VecGetBlockSize(u, &bs));
    530:   PetscCheck(Nc[f] == bs, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "Field %" PetscInt_FMT " components %" PetscInt_FMT " != %" PetscInt_FMT " blocksize for swarm field %s", f, Nc[f], bs, fieldnames[f]);
    531:   if (mode == SCATTER_FORWARD) {
    532:     PetscCall(DMSwarmProjectField_Conservative_PLEX(sw, dm, u, vec));
    533:   } else {
    534:     PetscCall(DMSwarmProjectParticles_Conservative_PLEX(sw, dm, u, vec));
    535:   }
    536:   PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, fieldnames[0], &u));
    537:   PetscFunctionReturn(PETSC_SUCCESS);
    538: }

    540: static PetscErrorCode DMSwarmProjectField_ApproxQ1_DA_2D(DM swarm, PetscReal *swarm_field, DM dm, Vec v_field)
    541: {
    542:   Vec                v_field_l, denom_l, coor_l, denom;
    543:   PetscScalar       *_field_l, *_denom_l;
    544:   PetscInt           k, p, e, npoints, nel, npe;
    545:   PetscInt          *mpfield_cell;
    546:   PetscReal         *mpfield_coor;
    547:   const PetscInt    *element_list;
    548:   const PetscInt    *element;
    549:   PetscScalar        xi_p[2], Ni[4];
    550:   const PetscScalar *_coor;

    552:   PetscFunctionBegin;
    553:   PetscCall(VecZeroEntries(v_field));

    555:   PetscCall(DMGetLocalVector(dm, &v_field_l));
    556:   PetscCall(DMGetGlobalVector(dm, &denom));
    557:   PetscCall(DMGetLocalVector(dm, &denom_l));
    558:   PetscCall(VecZeroEntries(v_field_l));
    559:   PetscCall(VecZeroEntries(denom));
    560:   PetscCall(VecZeroEntries(denom_l));

    562:   PetscCall(VecGetArray(v_field_l, &_field_l));
    563:   PetscCall(VecGetArray(denom_l, &_denom_l));

    565:   PetscCall(DMGetCoordinatesLocal(dm, &coor_l));
    566:   PetscCall(VecGetArrayRead(coor_l, &_coor));

    568:   PetscCall(DMDAGetElements(dm, &nel, &npe, &element_list));
    569:   PetscCall(DMSwarmGetLocalSize(swarm, &npoints));
    570:   PetscCall(DMSwarmGetField(swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&mpfield_coor));
    571:   PetscCall(DMSwarmGetField(swarm, DMSwarmPICField_cellid, NULL, NULL, (void **)&mpfield_cell));

    573:   for (p = 0; p < npoints; p++) {
    574:     PetscReal         *coor_p;
    575:     const PetscScalar *x0;
    576:     const PetscScalar *x2;
    577:     PetscScalar        dx[2];

    579:     e       = mpfield_cell[p];
    580:     coor_p  = &mpfield_coor[2 * p];
    581:     element = &element_list[npe * e];

    583:     /* compute local coordinates: (xp-x0)/dx = (xip+1)/2 */
    584:     x0 = &_coor[2 * element[0]];
    585:     x2 = &_coor[2 * element[2]];

    587:     dx[0] = x2[0] - x0[0];
    588:     dx[1] = x2[1] - x0[1];

    590:     xi_p[0] = 2.0 * (coor_p[0] - x0[0]) / dx[0] - 1.0;
    591:     xi_p[1] = 2.0 * (coor_p[1] - x0[1]) / dx[1] - 1.0;

    593:     /* evaluate basis functions */
    594:     Ni[0] = 0.25 * (1.0 - xi_p[0]) * (1.0 - xi_p[1]);
    595:     Ni[1] = 0.25 * (1.0 + xi_p[0]) * (1.0 - xi_p[1]);
    596:     Ni[2] = 0.25 * (1.0 + xi_p[0]) * (1.0 + xi_p[1]);
    597:     Ni[3] = 0.25 * (1.0 - xi_p[0]) * (1.0 + xi_p[1]);

    599:     for (k = 0; k < npe; k++) {
    600:       _field_l[element[k]] += Ni[k] * swarm_field[p];
    601:       _denom_l[element[k]] += Ni[k];
    602:     }
    603:   }

    605:   PetscCall(DMSwarmRestoreField(swarm, DMSwarmPICField_cellid, NULL, NULL, (void **)&mpfield_cell));
    606:   PetscCall(DMSwarmRestoreField(swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&mpfield_coor));
    607:   PetscCall(DMDARestoreElements(dm, &nel, &npe, &element_list));
    608:   PetscCall(VecRestoreArrayRead(coor_l, &_coor));
    609:   PetscCall(VecRestoreArray(v_field_l, &_field_l));
    610:   PetscCall(VecRestoreArray(denom_l, &_denom_l));

    612:   PetscCall(DMLocalToGlobalBegin(dm, v_field_l, ADD_VALUES, v_field));
    613:   PetscCall(DMLocalToGlobalEnd(dm, v_field_l, ADD_VALUES, v_field));
    614:   PetscCall(DMLocalToGlobalBegin(dm, denom_l, ADD_VALUES, denom));
    615:   PetscCall(DMLocalToGlobalEnd(dm, denom_l, ADD_VALUES, denom));

    617:   PetscCall(VecPointwiseDivide(v_field, v_field, denom));

    619:   PetscCall(DMRestoreLocalVector(dm, &v_field_l));
    620:   PetscCall(DMRestoreLocalVector(dm, &denom_l));
    621:   PetscCall(DMRestoreGlobalVector(dm, &denom));
    622:   PetscFunctionReturn(PETSC_SUCCESS);
    623: }

    625: static PetscErrorCode DMSwarmProjectFields_DA_Internal(DM swarm, DM celldm, PetscInt nfields, DMSwarmDataField dfield[], Vec vecs[], ScatterMode mode)
    626: {
    627:   PetscInt        f, dim;
    628:   DMDAElementType etype;

    630:   PetscFunctionBegin;
    631:   PetscCall(DMDAGetElementType(celldm, &etype));
    632:   PetscCheck(etype != DMDA_ELEMENT_P1, PetscObjectComm((PetscObject)swarm), PETSC_ERR_SUP, "Only Q1 DMDA supported");
    633:   PetscCheck(mode == SCATTER_FORWARD, PetscObjectComm((PetscObject)swarm), PETSC_ERR_SUP, "Mapping the continuum to particles is not currently supported for DMDA");

    635:   PetscCall(DMGetDimension(swarm, &dim));
    636:   switch (dim) {
    637:   case 2:
    638:     for (f = 0; f < nfields; f++) {
    639:       PetscReal *swarm_field;

    641:       PetscCall(DMSwarmDataFieldGetEntries(dfield[f], (void **)&swarm_field));
    642:       PetscCall(DMSwarmProjectField_ApproxQ1_DA_2D(swarm, swarm_field, celldm, vecs[f]));
    643:     }
    644:     break;
    645:   case 3:
    646:     SETERRQ(PetscObjectComm((PetscObject)swarm), PETSC_ERR_SUP, "No support for 3D");
    647:   default:
    648:     break;
    649:   }
    650:   PetscFunctionReturn(PETSC_SUCCESS);
    651: }

    653: /*@C
    654:   DMSwarmProjectFields - Project a set of swarm fields onto another `DM`

    656:   Collective

    658:   Input Parameters:
    659: + sw         - the `DMSWARM`
    660: . dm         - the `DM`, or `NULL` to use the cell `DM`
    661: . nfields    - the number of swarm fields to project
    662: . fieldnames - the textual names of the swarm fields to project
    663: . fields     - an array of `Vec`'s of length nfields
    664: - mode       - if `SCATTER_FORWARD` then map particles to the continuum, and if `SCATTER_REVERSE` map the continuum to particles

    666:   Level: beginner

    668:   Notes:
    669:   Currently, there are two available projection methods. The first is conservative projection, used for a `DMPLEX` cell `DM`.
    670:   The second is the averaging which is used for a `DMDA` cell `DM`

    672:   $$
    673:   \phi_i = \sum_{p=0}^{np} N_i(x_p) \phi_p dJ / \sum_{p=0}^{np} N_i(x_p) dJ
    674:   $$

    676:   where $\phi_p $ is the swarm field at point $p$, $N_i()$ is the cell `DM` basis function at vertex $i$, $dJ$ is the determinant of the cell Jacobian and
    677:   $\phi_i$ is the projected vertex value of the field $\phi$.

    679:   The user is responsible for destroying both the array and the individual `Vec` objects.

    681:   For the `DMPLEX` case, there is only a single vector, so the field layout in the `DMPLEX` must match the requested fields from the `DMSwarm`.

    683:   For averaging projection, nly swarm fields registered with data type of `PETSC_REAL` can be projected onto the cell `DM`, and only swarm fields of block size = 1 can currently be projected.

    685: .seealso: [](ch_dmbase), `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSetCellDM()`, `DMSwarmType`
    686: @*/
    687: PetscErrorCode DMSwarmProjectFields(DM sw, DM dm, PetscInt nfields, const char *fieldnames[], Vec fields[], ScatterMode mode)
    688: {
    689:   DM_Swarm         *swarm = (DM_Swarm *)sw->data;
    690:   DMSwarmDataField *gfield;
    691:   PetscBool         isDA, isPlex;
    692:   MPI_Comm          comm;

    694:   PetscFunctionBegin;
    695:   DMSWARMPICVALID(sw);
    696:   PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
    697:   if (!dm) PetscCall(DMSwarmGetCellDM(sw, &dm));
    698:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMDA, &isDA));
    699:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isPlex));
    700:   PetscCall(PetscMalloc1(nfields, &gfield));
    701:   for (PetscInt f = 0; f < nfields; ++f) PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldnames[f], &gfield[f]));

    703:   if (isDA) {
    704:     for (PetscInt f = 0; f < nfields; f++) {
    705:       PetscCheck(gfield[f]->petsc_type == PETSC_REAL, comm, PETSC_ERR_SUP, "Projection only valid for fields using a data type = PETSC_REAL");
    706:       PetscCheck(gfield[f]->bs == 1, comm, PETSC_ERR_SUP, "Projection only valid for fields with block size = 1");
    707:     }
    708:     PetscCall(DMSwarmProjectFields_DA_Internal(sw, dm, nfields, gfield, fields, mode));
    709:   } else if (isPlex) {
    710:     PetscInt Nf;

    712:     PetscCall(DMGetNumFields(dm, &Nf));
    713:     PetscCheck(Nf == nfields, comm, PETSC_ERR_ARG_WRONG, "Number of DM fields %" PetscInt_FMT " != %" PetscInt_FMT " number of requested Swarm fields", Nf, nfields);
    714:     PetscCall(DMSwarmProjectFields_Plex_Internal(sw, dm, nfields, fieldnames, fields[0], mode));
    715:   } else SETERRQ(PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "Only supported for cell DMs of type DMDA and DMPLEX");

    717:   PetscCall(PetscFree(gfield));
    718:   PetscFunctionReturn(PETSC_SUCCESS);
    719: }