Actual source code: dmproject.c

petsc-master 2020-07-14
Report Typos and Errors

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

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

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

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

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

 41:   for (f = 0; f < Nf; f++) {
 42:     u[f] = 1.;
 43:   }
 44:   return(0);
 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.  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
 50:   a least-squares solution.  It is also assumed that DMLocalToGlobalBegin()/DMLocalToGlobalEnd() with mode = ADD_VALUES is the adjoint of the
 51:   global-to-local map, so that the least-squares solution may be found by the normal equations.

 53:   collective

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

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

 63:   Level: advanced

 65:   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
 66:   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
 67:   closure of any local cell (see DMPlexGetAnchors()/DMPlexSetAnchors()).

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

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

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

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

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

160:   return(0);
161: }

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

166:   Collective on DM

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

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

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

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

203:   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.
204:   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
205:   a subdomain. You can also output a different number of fields than the input, with different discretizations. Last the auxiliary DM, attached to the
206:   auxiliary field vector, which is attached to dm, can also be different. It can have a different topology, number of fields, and discretizations.

208:   Level: intermediate

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

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

236:     DMGetDefaultConstraints(dm, NULL, &cMat);
237:     if (cMat) {
238:       DMGlobalToLocalSolve(dm, localX, X);
239:     }
240:   }
241:   DMRestoreLocalVector(dm, &localX);
242:   if (U != X) {DMRestoreLocalVector(dm, &localU);}
243:   return(0);
244: }