Actual source code: dmproject.c
petsc-3.13.6 2020-09-29
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: }