Actual source code: dmproject.c

petsc-3.11.3 2019-06-26
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, cEndInterior, 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:     DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
 89:     DMPlexGetHybridBounds(dm,&cEndInterior,NULL,NULL,NULL);
 90:     cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
 91:     if (cEnd > cStart) {
 92:       PetscScalar *ones;
 93:       PetscInt numValues, i;

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

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

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

162:   return(0);
163: }

165: /*@C
166:   DMProjectField - This projects the given function of the fields into the function space provided.

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:   Level: developer

180: .seealso: DMProjectFunction(), DMComputeL2Diff()
181: @*/
182: PetscErrorCode DMProjectField(DM dm, PetscReal time, Vec U,
183:                               void (**funcs)(PetscInt, PetscInt, PetscInt,
184:                                              const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
185:                                              const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
186:                                              PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
187:                               InsertMode mode, Vec X)
188: {
189:   Vec            localX, localU;

194:   DMGetLocalVector(dm, &localX);
195:   /* We currently check whether locU == locX to see if we need to apply BC */
196:   if (U != X) {DMGetLocalVector(dm, &localU);}
197:   else        {localU = localX;}
198:   DMGlobalToLocalBegin(dm, U, INSERT_VALUES, localU);
199:   DMGlobalToLocalEnd(dm, U, INSERT_VALUES, localU);
200:   DMProjectFieldLocal(dm, time, localU, funcs, mode, localX);
201:   DMLocalToGlobalBegin(dm, localX, mode, X);
202:   DMLocalToGlobalEnd(dm, localX, mode, X);
203:   if (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES) {
204:     Mat cMat;

206:     DMGetDefaultConstraints(dm, NULL, &cMat);
207:     if (cMat) {
208:       DMGlobalToLocalSolve(dm, localX, X);
209:     }
210:   }
211:   DMRestoreLocalVector(dm, &localX);
212:   if (U != X) {DMRestoreLocalVector(dm, &localU);}
213:   return(0);
214: }