Actual source code: dmproject.c

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

      8: typedef struct _projectConstraintsCtx
      9: {
     10:   DM  dm;
     11:   Vec mask;
     12: }
     13: projectConstraintsCtx;

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

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

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

     42:   for (f = 0; f < Nf; f++) {
     43:     u[f] = 1.;
     44:   }
     45:   return(0);
     46: }

     48: /*@
     49:   DMGlobalToLocalSolve - Solve for the global vector that is mapped to a given local vector by DMGlobalToLocalBegin()/DMGlobalToLocalEnd() with mode
     50:   = INSERT_VALUES.  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
     51:   a least-squares solution.  It is also assumed that DMLocalToGlobalBegin()/DMLocalToGlobalEnd() with mode = ADD_VALUES is the adjoint of the
     52:   global-to-local map, so that the least-squares solution may be found by the normal equations.

     54:   collective

     56:   Input Parameters:
     57: + dm - The DM object
     58: . x - The local vector
     59: - y - The global vector: the input value of globalVec is used as an initial guess

     61:   Output Parameters:
     62: . y - The least-squares solution

     64:   Level: advanced

     66:   Note: If the DM is of type DMPLEX, then y is the solution of L' * D * L * y = L' * D * x, where D is a diagonal mask that is 1 for every point in
     67:   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
     68:   closure of any local cell (see DMPlexGetAnchors()/DMPlexSetAnchors()).

     70: .seealso: DMGlobalToLocalBegin(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMLocalToGlobalEnd(), DMPlexGetAnchors(), DMPlexSetAnchors()
     71: @*/
     72: PetscErrorCode DMGlobalToLocalSolve(DM dm, Vec x, Vec y)
     73: {
     74:   Mat                   CtC;
     75:   PetscInt              n, N, cStart, cEnd, c;
     76:   PetscBool             isPlex;
     77:   KSP                   ksp;
     78:   PC                    pc;
     79:   Vec                   global, mask=NULL;
     80:   projectConstraintsCtx ctx;

     84:   PetscObjectTypeCompare((PetscObject)dm,DMPLEX,&isPlex);
     85:   if (isPlex) {
     86:     /* mark points in the closure */
     87:     DMCreateLocalVector(dm,&mask);
     88:     VecSet(mask,0.0);
     89:     DMPlexGetSimplexOrBoxCells(dm,0,&cStart,&cEnd);
     90:     if (cEnd > cStart) {
     91:       PetscScalar *ones;
     92:       PetscInt numValues, i;

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

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

    161:   return(0);
    162: }

    164: /*@C
    165:   DMProjectField - This projects the given function of the input fields into the function space provided, putting the coefficients in a global vector.

    167:   Collective on DM

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

    176:   Output Parameter:
    177: . X       - The output vector

    179:    Calling sequence of func:
    180: $    func(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    181: $         const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    182: $         const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    183: $         PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[]);

    185: +  dim          - The spatial dimension
    186: .  Nf           - The number of input fields
    187: .  NfAux        - The number of input auxiliary fields
    188: .  uOff         - The offset of each field in u[]
    189: .  uOff_x       - The offset of each field in u_x[]
    190: .  u            - The field values at this point in space
    191: .  u_t          - The field time derivative at this point in space (or NULL)
    192: .  u_x          - The field derivatives at this point in space
    193: .  aOff         - The offset of each auxiliary field in u[]
    194: .  aOff_x       - The offset of each auxiliary field in u_x[]
    195: .  a            - The auxiliary field values at this point in space
    196: .  a_t          - The auxiliary field time derivative at this point in space (or NULL)
    197: .  a_x          - The auxiliary field derivatives at this point in space
    198: .  t            - The current time
    199: .  x            - The coordinates of this point
    200: .  numConstants - The number of constants
    201: .  constants    - The value of each constant
    202: -  f            - The value of the function at this point in space

    204:   Note: There are three different DMs that potentially interact in this function. The output DM, dm, specifies the layout of the values calculates by funcs.
    205:   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
    206:   a subdomain. You can also output a different number of fields than the input, with different discretizations. Last the auxiliary DM, attached to the
    207:   auxiliary field vector, which is attached to dm, can also be different. It can have a different topology, number of fields, and discretizations.

    209:   Level: intermediate

    211: .seealso: DMProjectFieldLocal(), DMProjectFieldLabelLocal(), DMProjectFunction(), DMComputeL2Diff()
    212: @*/
    213: PetscErrorCode DMProjectField(DM dm, PetscReal time, Vec U,
    214:                               void (**funcs)(PetscInt, PetscInt, PetscInt,
    215:                                              const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
    216:                                              const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
    217:                                              PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
    218:                               InsertMode mode, Vec X)
    219: {
    220:   Vec            localX, localU;
    221:   DM             dmIn;

    226:   DMGetLocalVector(dm, &localX);
    227:   /* We currently check whether locU == locX to see if we need to apply BC */
    228:   if (U != X) {
    229:     VecGetDM(U, &dmIn);
    230:     DMGetLocalVector(dmIn, &localU);
    231:   } else {
    232:     dmIn   = dm;
    233:     localU = localX;
    234:   }
    235:   DMGlobalToLocalBegin(dmIn, U, INSERT_VALUES, localU);
    236:   DMGlobalToLocalEnd(dmIn, U, INSERT_VALUES, localU);
    237:   DMProjectFieldLocal(dm, time, localU, funcs, mode, localX);
    238:   DMLocalToGlobalBegin(dm, localX, mode, X);
    239:   DMLocalToGlobalEnd(dm, localX, mode, X);
    240:   if (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES) {
    241:     Mat cMat;

    243:     DMGetDefaultConstraints(dm, NULL, &cMat);
    244:     if (cMat) {
    245:       DMGlobalToLocalSolve(dm, localX, X);
    246:     }
    247:   }
    248:   DMRestoreLocalVector(dm, &localX);
    249:   if (U != X) {DMRestoreLocalVector(dmIn, &localU);}
    250:   return(0);
    251: }

    253: /********************* Adaptive Interpolation **************************/

    255: /* See the discussion of Adaptive Interpolation in manual/high_level_mg.rst */
    256: PetscErrorCode DMAdaptInterpolator(DM dmc, DM dmf, Mat In, KSP smoother, PetscInt Nc, Vec vf[], Vec vc[], Mat *InAdapt, void *user)
    257: {
    258:   Mat            globalA;
    259:   Vec            tmp, tmp2;
    260:   PetscScalar   *A, *b, *x, *workscalar;
    261:   PetscReal     *w, *sing, *workreal, rcond = PETSC_SMALL;
    262:   PetscBLASInt   M, N, one = 1, irank, lwrk, info;
    263:   PetscInt       debug = 0, rStart, rEnd, r, maxcols = 0, k;
    264:   PetscBool      allocVc = PETSC_FALSE;

    268:   PetscLogEventBegin(DM_AdaptInterpolator,dmc,dmf,0,0);
    269:   PetscOptionsGetInt(NULL, NULL, "-dm_interpolator_adapt_debug", &debug, NULL);
    270:   MatDuplicate(In, MAT_SHARE_NONZERO_PATTERN, InAdapt);
    271:   MatGetOwnershipRange(In, &rStart, &rEnd);
    272:   #if 0
    273:   MatGetMaxRowLen(In, &maxcols);
    274:   #else
    275:   for (r = rStart; r < rEnd; ++r) {
    276:     PetscInt           ncols;
    277:     const PetscInt    *cols;
    278:     const PetscScalar *vals;

    280:     MatGetRow(In, r, &ncols, &cols, &vals);
    281:     maxcols = PetscMax(maxcols, ncols);
    282:     MatRestoreRow(In, r, &ncols, &cols, &vals);
    283:   }
    284:   #endif
    285:   if (Nc < maxcols) PetscPrintf(PETSC_COMM_SELF, "The number of input vectors %D < %D the maximum number of column entries\n", Nc, maxcols);
    286:   for (k = 0; k < Nc; ++k) {
    287:     char        name[PETSC_MAX_PATH_LEN];
    288:     const char *prefix;

    290:     PetscObjectGetOptionsPrefix((PetscObject) smoother, &prefix);
    291:     PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "%sCoarse Vector %D", prefix ? prefix : NULL, k);
    292:     PetscObjectSetName((PetscObject) vc[k], name);
    293:     VecViewFromOptions(vc[k], NULL, "-dm_adapt_interp_view_coarse");
    294:     PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "%sFine Vector %D", prefix ? prefix : NULL, k);
    295:     PetscObjectSetName((PetscObject) vf[k], name);
    296:     VecViewFromOptions(vf[k], NULL, "-dm_adapt_interp_view_fine");
    297:   }
    298:   PetscBLASIntCast(3*PetscMin(Nc, maxcols) + PetscMax(2*PetscMin(Nc, maxcols), PetscMax(Nc, maxcols)), &lwrk);
    299:   PetscMalloc7(Nc*maxcols, &A, PetscMax(Nc, maxcols), &b, Nc, &w, maxcols, &x, maxcols, &sing, lwrk, &workscalar, 5*PetscMin(Nc, maxcols), &workreal);
    300:   /* 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} */
    301:   KSPGetOperators(smoother, &globalA, NULL);
    302:   DMGetGlobalVector(dmf, &tmp);
    303:   DMGetGlobalVector(dmf, &tmp2);
    304:   for (k = 0; k < Nc; ++k) {
    305:     PetscScalar vnorm, vAnorm;
    306:     PetscBool   canMult = PETSC_FALSE;
    307:     const char *type;

    309:     w[k] = 1.0;
    310:     PetscObjectGetType((PetscObject) globalA, &type);
    311:     if (type) {MatAssembled(globalA, &canMult);}
    312:     if (type && canMult) {
    313:       VecDot(vf[k], vf[k], &vnorm);
    314:       MatMult(globalA, vf[k], tmp);
    315: #if 0
    316:       KSPSolve(smoother, tmp, tmp2);
    317:       VecDot(vf[k], tmp2, &vAnorm);
    318: #else
    319:       VecDot(vf[k], tmp, &vAnorm);
    320: #endif
    321:       w[k] = PetscRealPart(vnorm) / PetscRealPart(vAnorm);
    322:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "System matrix is not assembled.");
    323:   }
    324:   DMRestoreGlobalVector(dmf, &tmp);
    325:   DMRestoreGlobalVector(dmf, &tmp2);
    326:   if (!vc) {
    327:     allocVc = PETSC_TRUE;
    328:     PetscMalloc1(Nc, &vc);
    329:     for (k = 0; k < Nc; ++k) {
    330:       DMGetGlobalVector(dmc, &vc[k]);
    331:       MatMultTranspose(In, vf[k], vc[k]);
    332:     }
    333:   }
    334:   /* Solve a LS system for each fine row */
    335:   for (r = rStart; r < rEnd; ++r) {
    336:     PetscInt           ncols, c;
    337:     const PetscInt    *cols;
    338:     const PetscScalar *vals, *a;

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

    345:       /* b_k = \sqrt{w_k} f^{F,k}_r */
    346:       VecGetArrayRead(vf[k], &a);
    347:       b[k] = wk * a[r-rStart];
    348:       VecRestoreArrayRead(vf[k], &a);
    349:       /* A_{kc} = \sqrt{w_k} f^{C,k}_c */
    350:       /* TODO Must pull out VecScatter from In, scatter in vc[k] values up front, and access them indirectly just as in MatMult() */
    351:       VecGetArrayRead(vc[k], &a);
    352:       for (c = 0; c < ncols; ++c) {
    353:         /* This is element (k, c) of A */
    354:         A[c*Nc+k] = wk * a[cols[c]-rStart];
    355:       }
    356:       VecRestoreArrayRead(vc[k], &a);
    357:     }
    358:     PetscBLASIntCast(Nc,    &M);
    359:     PetscBLASIntCast(ncols, &N);
    360:     if (debug) {
    361: #if defined(PETSC_USE_COMPLEX)
    362:       PetscScalar *tmp;
    363:       PetscInt     j;

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

    394:         DMGetWorkArray(dmc, Nc, MPIU_SCALAR, (void *) &tmp);
    395:         for (j = 0; j < PetscMin(Nc, ncols); ++j) tmp[j] = sing[j];
    396:         DMPrintCellMatrix(r, "Interpolator Row LS singular values", PetscMin(Nc, ncols), 1, tmp);
    397:         DMRestoreWorkArray(dmc, Nc, MPIU_SCALAR, (void *) &tmp);
    398:       }
    399: #else
    400:       DMPrintCellMatrix(r, "Interpolator Row LS singular values", PetscMin(Nc, ncols), 1, sing);
    401: #endif
    402:       DMPrintCellMatrix(r, "Interpolator Row LS old P", ncols, 1, vals);
    403:       DMPrintCellMatrix(r, "Interpolator Row LS sol", ncols, 1, b);
    404:     }
    405:     MatSetValues(*InAdapt, 1, &r, ncols, cols, b, INSERT_VALUES);
    406:     MatRestoreRow(In, r, &ncols, &cols, &vals);
    407:   }
    408:   PetscFree7(A, b, w, x, sing, workscalar, workreal);
    409:   if (allocVc) {
    410:     for (k = 0; k < Nc; ++k) {DMRestoreGlobalVector(dmc, &vc[k]);}
    411:     PetscFree(vc);
    412:   }
    413:   MatAssemblyBegin(*InAdapt, MAT_FINAL_ASSEMBLY);
    414:   MatAssemblyEnd(*InAdapt, MAT_FINAL_ASSEMBLY);
    415:   PetscLogEventEnd(DM_AdaptInterpolator,dmc,dmf,0,0);
    416:   return(0);
    417: }

    419: PetscErrorCode DMCheckInterpolator(DM dmf, Mat In, PetscInt Nc, Vec vc[], Vec vf[], PetscReal tol)
    420: {
    421:   Vec            tmp;
    422:   PetscReal      norminf, norm2, maxnorminf = 0.0, maxnorm2 = 0.0;
    423:   PetscInt       k;

    427:   DMGetGlobalVector(dmf, &tmp);
    428:   MatViewFromOptions(In, NULL, "-dm_interpolator_adapt_error");
    429:   for (k = 0; k < Nc; ++k) {
    430:     MatMult(In, vc[k], tmp);
    431:     VecAXPY(tmp, -1.0, vf[k]);
    432:     VecViewFromOptions(vc[k], NULL, "-dm_interpolator_adapt_error");
    433:     VecViewFromOptions(vf[k], NULL, "-dm_interpolator_adapt_error");
    434:     VecViewFromOptions(tmp, NULL, "-dm_interpolator_adapt_error");
    435:     VecNorm(tmp, NORM_INFINITY, &norminf);
    436:     VecNorm(tmp, NORM_2, &norm2);
    437:     maxnorminf = PetscMax(maxnorminf, norminf);
    438:     maxnorm2   = PetscMax(maxnorm2,   norm2);
    439:     PetscPrintf(PetscObjectComm((PetscObject) dmf), "Coarse vec %D ||vf - P vc||_\\infty %g, ||vf - P vc||_2 %g\n", k, norminf, norm2);
    440:   }
    441:   DMRestoreGlobalVector(dmf, &tmp);
    442:   if (maxnorm2 > tol) SETERRQ2(PetscObjectComm((PetscObject) dmf), PETSC_ERR_ARG_WRONG, "max_k ||vf_k - P vc_k||_2 %g > tol %g", maxnorm2, tol);
    443:   return(0);
    444: }