Actual source code: dmproject.c

petsc-3.13.4 2020-08-01
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;
11: }
12: projectConstraintsCtx;

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

22:   MatShellGetContext(CtC,&ctx);
23:   dm   = ctx->dm;
25:   DMGetLocalVector(dm,&local);
26:   DMGlobalToLocalBegin(dm,x,INSERT_VALUES,local);
27:   DMGlobalToLocalEnd(dm,x,INSERT_VALUES,local);
29:   VecSet(y,0.);
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

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;
79:   projectConstraintsCtx ctx;

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

94:       PetscMalloc1(numValues,&ones);
95:       for (i = 0; i < numValues; i++) {
96:         ones[i] = 1.;
97:       }
98:       for (c = cStart; c < cEnd; c++) {
100:       }
101:       PetscFree(ones);
102:     }
103:   }
104:   else {

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:       }
122:       PetscFree2(func, ctx);
123:     }
125:   }
126:   ctx.dm   = dm;
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.);
148:   KSPSolve(ksp,global,y);
149:   DMRestoreGlobalVector(dm,&global);
150:   /* clean up */
151:   KSPDestroy(&ksp);
152:   MatDestroy(&CtC);
153:   if (isPlex) {
155:   }
156:   else {
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: }