1 #include <../src/snes/impls/vi/rs/virsimpl.h> /*I "petscsnes.h" I*/
2 #include <petsc/private/dmimpl.h>
3 #include <petsc/private/vecimpl.h>
4
5 /*@
6 SNESVIGetInactiveSet - Gets the global indices for the inactive set variables (these correspond to the degrees of freedom the linear
7 system is solved on)
8
9 Input Parameter:
10 . snes - the `SNES` context
11
12 Output Parameter:
13 . inact - inactive set index set
14
15 Level: advanced
16
17 .seealso: [](ch_snes), `SNES`, `SNESVINEWTONRSLS`
18 @*/
SNESVIGetInactiveSet(SNES snes,IS * inact)19 PetscErrorCode SNESVIGetInactiveSet(SNES snes, IS *inact)
20 {
21 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;
22
23 PetscFunctionBegin;
24 *inact = vi->IS_inact;
25 PetscFunctionReturn(PETSC_SUCCESS);
26 }
27
28 /*
29 Provides a wrapper to a DM to allow it to be used to generated the interpolation/restriction from the DM for the smaller matrices and vectors
30 defined by the reduced space method.
31
32 Simple calls the regular DM interpolation and restricts it to operation on the variables not associated with active constraints.
33 */
34 typedef struct {
35 PetscInt n; /* size of vectors in the reduced DM space */
36 IS inactive;
37
38 PetscErrorCode (*createinterpolation)(DM, DM, Mat *, Vec *); /* DM's original routines */
39 PetscErrorCode (*coarsen)(DM, MPI_Comm, DM *);
40 PetscErrorCode (*createglobalvector)(DM, Vec *);
41 PetscErrorCode (*createinjection)(DM, DM, Mat *);
42 PetscErrorCode (*hascreateinjection)(DM, PetscBool *);
43
44 DM dm; /* when destroying this object we need to reset the above function into the base DM */
45 } DM_SNESVI;
46
47 /*
48 DMCreateGlobalVector_SNESVI - Creates global vector of the size of the reduced space
49 */
DMCreateGlobalVector_SNESVI(DM dm,Vec * vec)50 static PetscErrorCode DMCreateGlobalVector_SNESVI(DM dm, Vec *vec)
51 {
52 PetscContainer isnes;
53 DM_SNESVI *dmsnesvi;
54
55 PetscFunctionBegin;
56 PetscCall(PetscObjectQuery((PetscObject)dm, "VI", (PetscObject *)&isnes));
57 PetscCheck(isnes, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Composed SNES is missing");
58 PetscCall(PetscContainerGetPointer(isnes, &dmsnesvi));
59 PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)dm), dmsnesvi->n, PETSC_DETERMINE, vec));
60 PetscCall(VecSetDM(*vec, dm));
61 PetscFunctionReturn(PETSC_SUCCESS);
62 }
63
64 /*
65 DMCreateInterpolation_SNESVI - Modifieds the interpolation obtained from the DM by removing all rows and columns associated with active constraints.
66 */
DMCreateInterpolation_SNESVI(DM dm1,DM dm2,Mat * mat,Vec * vec)67 static PetscErrorCode DMCreateInterpolation_SNESVI(DM dm1, DM dm2, Mat *mat, Vec *vec)
68 {
69 PetscContainer isnes;
70 DM_SNESVI *dmsnesvi1, *dmsnesvi2;
71 Mat interp;
72
73 PetscFunctionBegin;
74 PetscCall(PetscObjectQuery((PetscObject)dm1, "VI", (PetscObject *)&isnes));
75 PetscCheck(isnes, PetscObjectComm((PetscObject)dm1), PETSC_ERR_PLIB, "Composed VI data structure is missing");
76 PetscCall(PetscContainerGetPointer(isnes, &dmsnesvi1));
77 PetscCall(PetscObjectQuery((PetscObject)dm2, "VI", (PetscObject *)&isnes));
78 PetscCheck(isnes, PetscObjectComm((PetscObject)dm2), PETSC_ERR_PLIB, "Composed VI data structure is missing");
79 PetscCall(PetscContainerGetPointer(isnes, &dmsnesvi2));
80
81 PetscCall((*dmsnesvi1->createinterpolation)(dm1, dm2, &interp, NULL));
82 PetscCall(MatCreateSubMatrix(interp, dmsnesvi2->inactive, dmsnesvi1->inactive, MAT_INITIAL_MATRIX, mat));
83 PetscCall(MatDestroy(&interp));
84 *vec = NULL;
85 PetscFunctionReturn(PETSC_SUCCESS);
86 }
87
88 /*
89 DMCoarsen_SNESVI - Computes the regular coarsened DM then computes additional information about its inactive set
90 */
DMCoarsen_SNESVI(DM dm1,MPI_Comm comm,DM * dm2)91 static PetscErrorCode DMCoarsen_SNESVI(DM dm1, MPI_Comm comm, DM *dm2)
92 {
93 PetscContainer isnes;
94 DM_SNESVI *dmsnesvi1;
95 Vec finemarked, coarsemarked;
96 IS inactive;
97 Mat inject;
98 const PetscInt *index;
99 PetscInt n, k, cnt = 0, rstart, *coarseindex;
100 PetscScalar *marked;
101
102 PetscFunctionBegin;
103 PetscCall(PetscObjectQuery((PetscObject)dm1, "VI", (PetscObject *)&isnes));
104 PetscCheck(isnes, PetscObjectComm((PetscObject)dm1), PETSC_ERR_PLIB, "Composed VI data structure is missing");
105 PetscCall(PetscContainerGetPointer(isnes, &dmsnesvi1));
106
107 /* get the original coarsen */
108 PetscCall((*dmsnesvi1->coarsen)(dm1, comm, dm2));
109
110 /* not sure why this extra reference is needed, but without the dm2 disappears too early */
111 /* Updating the KSPCreateVecs() to avoid using DMGetGlobalVector() when matrix is available removes the need for this reference? */
112 /* PetscCall(PetscObjectReference((PetscObject)*dm2));*/
113
114 /* need to set back global vectors in order to use the original injection */
115 PetscCall(DMClearGlobalVectors(dm1));
116
117 dm1->ops->createglobalvector = dmsnesvi1->createglobalvector;
118
119 PetscCall(DMCreateGlobalVector(dm1, &finemarked));
120 PetscCall(DMCreateGlobalVector(*dm2, &coarsemarked));
121
122 /*
123 fill finemarked with locations of inactive points
124 */
125 PetscCall(ISGetIndices(dmsnesvi1->inactive, &index));
126 PetscCall(ISGetLocalSize(dmsnesvi1->inactive, &n));
127 PetscCall(VecSet(finemarked, 0.0));
128 for (k = 0; k < n; k++) PetscCall(VecSetValue(finemarked, index[k], 1.0, INSERT_VALUES));
129 PetscCall(VecAssemblyBegin(finemarked));
130 PetscCall(VecAssemblyEnd(finemarked));
131
132 PetscCall(DMCreateInjection(*dm2, dm1, &inject));
133 PetscCall(MatRestrict(inject, finemarked, coarsemarked));
134 PetscCall(MatDestroy(&inject));
135
136 /*
137 create index set list of coarse inactive points from coarsemarked
138 */
139 PetscCall(VecGetLocalSize(coarsemarked, &n));
140 PetscCall(VecGetOwnershipRange(coarsemarked, &rstart, NULL));
141 PetscCall(VecGetArray(coarsemarked, &marked));
142 for (k = 0; k < n; k++) {
143 if (marked[k] != 0.0) cnt++;
144 }
145 PetscCall(PetscMalloc1(cnt, &coarseindex));
146 cnt = 0;
147 for (k = 0; k < n; k++) {
148 if (marked[k] != 0.0) coarseindex[cnt++] = k + rstart;
149 }
150 PetscCall(VecRestoreArray(coarsemarked, &marked));
151 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)coarsemarked), cnt, coarseindex, PETSC_OWN_POINTER, &inactive));
152
153 PetscCall(DMClearGlobalVectors(dm1));
154
155 dm1->ops->createglobalvector = DMCreateGlobalVector_SNESVI;
156
157 PetscCall(DMSetVI(*dm2, inactive));
158
159 PetscCall(VecDestroy(&finemarked));
160 PetscCall(VecDestroy(&coarsemarked));
161 PetscCall(ISDestroy(&inactive));
162 PetscFunctionReturn(PETSC_SUCCESS);
163 }
164
DMDestroy_SNESVI(PetscCtxRt ctx)165 static PetscErrorCode DMDestroy_SNESVI(PetscCtxRt ctx)
166 {
167 DM_SNESVI *dmsnesvi = *(DM_SNESVI **)ctx;
168
169 PetscFunctionBegin;
170 /* reset the base methods in the DM object that were changed when the DM_SNESVI was reset */
171 dmsnesvi->dm->ops->createinterpolation = dmsnesvi->createinterpolation;
172 dmsnesvi->dm->ops->coarsen = dmsnesvi->coarsen;
173 dmsnesvi->dm->ops->createglobalvector = dmsnesvi->createglobalvector;
174 dmsnesvi->dm->ops->createinjection = dmsnesvi->createinjection;
175 dmsnesvi->dm->ops->hascreateinjection = dmsnesvi->hascreateinjection;
176 /* need to clear out this vectors because some of them may not have a reference to the DM
177 but they are counted as having references to the DM in DMDestroy() */
178 PetscCall(DMClearGlobalVectors(dmsnesvi->dm));
179
180 PetscCall(ISDestroy(&dmsnesvi->inactive));
181 PetscCall(PetscFree(dmsnesvi));
182 PetscFunctionReturn(PETSC_SUCCESS);
183 }
184
185 /*@
186 DMSetVI - Marks a `DM` as associated with a VI problem. This causes the interpolation/restriction operators to
187 be restricted to only those variables NOT associated with active constraints.
188
189 Logically Collective
190
191 Input Parameters:
192 + dm - the `DM` object
193 - inactive - an `IS` indicating which points are currently not active
194
195 Level: intermediate
196
197 .seealso: [](ch_snes), `SNES`, `SNESVINEWTONRSLS`, `SNESVIGetInactiveSet()`
198 @*/
DMSetVI(DM dm,IS inactive)199 PetscErrorCode DMSetVI(DM dm, IS inactive)
200 {
201 PetscContainer isnes;
202 DM_SNESVI *dmsnesvi;
203
204 PetscFunctionBegin;
205 if (!dm) PetscFunctionReturn(PETSC_SUCCESS);
206
207 PetscCall(PetscObjectReference((PetscObject)inactive));
208
209 PetscCall(PetscObjectQuery((PetscObject)dm, "VI", (PetscObject *)&isnes));
210 if (!isnes) {
211 PetscCall(PetscNew(&dmsnesvi));
212 PetscCall(PetscObjectContainerCompose((PetscObject)dm, "VI", dmsnesvi, DMDestroy_SNESVI));
213
214 dmsnesvi->createinterpolation = dm->ops->createinterpolation;
215 dm->ops->createinterpolation = DMCreateInterpolation_SNESVI;
216 dmsnesvi->coarsen = dm->ops->coarsen;
217 dm->ops->coarsen = DMCoarsen_SNESVI;
218 dmsnesvi->createglobalvector = dm->ops->createglobalvector;
219 dm->ops->createglobalvector = DMCreateGlobalVector_SNESVI;
220 dmsnesvi->createinjection = dm->ops->createinjection;
221 dm->ops->createinjection = NULL;
222 dmsnesvi->hascreateinjection = dm->ops->hascreateinjection;
223 dm->ops->hascreateinjection = NULL;
224 } else {
225 PetscCall(PetscContainerGetPointer(isnes, &dmsnesvi));
226 PetscCall(ISDestroy(&dmsnesvi->inactive));
227 }
228 PetscCall(DMClearGlobalVectors(dm));
229 PetscCall(ISGetLocalSize(inactive, &dmsnesvi->n));
230
231 dmsnesvi->inactive = inactive;
232 dmsnesvi->dm = dm;
233 PetscFunctionReturn(PETSC_SUCCESS);
234 }
235
236 /*
237 DMDestroyVI - Frees the DM_SNESVI object contained in the DM
238 - also resets the function pointers in the DM for createinterpolation() etc to use the original DM
239 */
DMDestroyVI(DM dm)240 PetscErrorCode DMDestroyVI(DM dm)
241 {
242 PetscFunctionBegin;
243 if (!dm) PetscFunctionReturn(PETSC_SUCCESS);
244 PetscCall(PetscObjectCompose((PetscObject)dm, "VI", (PetscObject)NULL));
245 PetscFunctionReturn(PETSC_SUCCESS);
246 }
247
248 /* Create active and inactive set vectors. The local size of this vector is set and PETSc computes the global size */
SNESCreateSubVectors_VINEWTONRSLS(SNES snes,PetscInt n,Vec * newv)249 static PetscErrorCode SNESCreateSubVectors_VINEWTONRSLS(SNES snes, PetscInt n, Vec *newv)
250 {
251 Vec v;
252
253 PetscFunctionBegin;
254 PetscCall(VecCreate(PetscObjectComm((PetscObject)snes), &v));
255 PetscCall(VecSetSizes(v, n, PETSC_DECIDE));
256 PetscCall(VecSetType(v, VECSTANDARD));
257 *newv = v;
258 PetscFunctionReturn(PETSC_SUCCESS);
259 }
260
261 /* Resets the snes PC and KSP when the active set sizes change */
SNESVIResetPCandKSP(SNES snes,Mat Amat,Mat Pmat)262 static PetscErrorCode SNESVIResetPCandKSP(SNES snes, Mat Amat, Mat Pmat)
263 {
264 KSP snesksp;
265
266 PetscFunctionBegin;
267 PetscCall(SNESGetKSP(snes, &snesksp));
268 PetscCall(KSPReset(snesksp));
269 PetscCall(KSPResetFromOptions(snesksp));
270
271 /*
272 KSP kspnew;
273 PC pcnew;
274 MatSolverType stype;
275
276 PetscCall(KSPCreate(PetscObjectComm((PetscObject)snes),&kspnew));
277 kspnew->pc_side = snesksp->pc_side;
278 kspnew->rtol = snesksp->rtol;
279 kspnew->abstol = snesksp->abstol;
280 kspnew->max_it = snesksp->max_it;
281 PetscCall(KSPSetType(kspnew,((PetscObject)snesksp)->type_name));
282 PetscCall(KSPGetPC(kspnew,&pcnew));
283 PetscCall(PCSetType(kspnew->pc,((PetscObject)snesksp->pc)->type_name));
284 PetscCall(PCSetOperators(kspnew->pc,Amat,Pmat));
285 PetscCall(PCFactorGetMatSolverType(snesksp->pc,&stype));
286 PetscCall(PCFactorSetMatSolverType(kspnew->pc,stype));
287 PetscCall(KSPDestroy(&snesksp));
288 snes->ksp = kspnew;
289 PetscCall(KSPSetFromOptions(kspnew));*/
290 PetscFunctionReturn(PETSC_SUCCESS);
291 }
292
293 /* Variational Inequality solver using reduce space method. No semismooth algorithm is
294 implemented in this algorithm. It basically identifies the active constraints and does
295 a linear solve on the other variables (those not associated with the active constraints). */
SNESSolve_VINEWTONRSLS(SNES snes)296 static PetscErrorCode SNESSolve_VINEWTONRSLS(SNES snes)
297 {
298 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;
299 PetscInt maxits, i, lits;
300 SNESLineSearchReason lsreason;
301 PetscReal fnorm, gnorm, xnorm = 0, ynorm;
302 Vec Y, X, F;
303 KSPConvergedReason kspreason;
304 KSP ksp;
305 PC pc;
306
307 PetscFunctionBegin;
308 /* Multigrid must use Galerkin for coarse grids with active set/reduced space methods; cannot rediscretize on coarser grids*/
309 PetscCall(SNESGetKSP(snes, &ksp));
310 PetscCall(KSPGetPC(ksp, &pc));
311 PetscCall(PCMGSetGalerkin(pc, PC_MG_GALERKIN_BOTH));
312
313 snes->numFailures = 0;
314 snes->numLinearSolveFailures = 0;
315 snes->reason = SNES_CONVERGED_ITERATING;
316
317 maxits = snes->max_its; /* maximum number of iterations */
318 X = snes->vec_sol; /* solution vector */
319 F = snes->vec_func; /* residual vector */
320 Y = snes->work[0]; /* work vectors */
321
322 PetscCall(SNESLineSearchSetVIFunctions(snes->linesearch, SNESVIProjectOntoBounds, SNESVIComputeInactiveSetFnorm, SNESVIComputeInactiveSetFtY));
323 PetscCall(SNESLineSearchSetVecs(snes->linesearch, X, NULL, NULL, NULL, NULL));
324 PetscCall(SNESLineSearchSetUp(snes->linesearch));
325
326 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
327 snes->iter = 0;
328 snes->norm = 0.0;
329 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
330
331 PetscCall(SNESVIProjectOntoBounds(snes, X));
332 PetscCall(SNESComputeFunction(snes, X, F));
333 PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm));
334 PetscCall(VecNorm(X, NORM_2, &xnorm)); /* xnorm <- ||x|| */
335 SNESCheckFunctionDomainError(snes, fnorm);
336 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
337 snes->norm = fnorm;
338 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
339 PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0));
340
341 /* test convergence */
342 PetscCall(SNESConverged(snes, 0, 0.0, 0.0, fnorm));
343 PetscCall(SNESMonitor(snes, 0, fnorm));
344 if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
345
346 for (i = 0; i < maxits; i++) {
347 IS IS_act; /* _act -> active set _inact -> inactive set */
348 IS IS_redact; /* redundant active set */
349 VecScatter scat_act, scat_inact;
350 PetscInt nis_act, nis_inact;
351 Vec Y_act, Y_inact, F_inact;
352 Mat jac_inact_inact, prejac_inact_inact;
353 PetscBool isequal;
354
355 /* Call general purpose update function */
356 PetscTryTypeMethod(snes, update, snes->iter);
357 PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
358 SNESCheckJacobianDomainError(snes);
359
360 /* Create active and inactive index sets */
361
362 /*original
363 PetscCall(SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&vi->IS_inact));
364 */
365 PetscCall(SNESVIGetActiveSetIS(snes, X, F, &IS_act));
366
367 if (vi->checkredundancy) {
368 PetscCall((*vi->checkredundancy)(snes, IS_act, &IS_redact, vi->ctxP));
369 if (IS_redact) {
370 PetscCall(ISSort(IS_redact));
371 PetscCall(ISComplement(IS_redact, X->map->rstart, X->map->rend, &vi->IS_inact));
372 PetscCall(ISDestroy(&IS_redact));
373 } else {
374 PetscCall(ISComplement(IS_act, X->map->rstart, X->map->rend, &vi->IS_inact));
375 }
376 } else {
377 PetscCall(ISComplement(IS_act, X->map->rstart, X->map->rend, &vi->IS_inact));
378 }
379
380 /* Create inactive set submatrix */
381 PetscCall(MatCreateSubMatrix(snes->jacobian, vi->IS_inact, vi->IS_inact, MAT_INITIAL_MATRIX, &jac_inact_inact));
382
383 if (0) { /* Dead code (temporary developer hack) */
384 IS keptrows;
385 PetscCall(MatFindNonzeroRows(jac_inact_inact, &keptrows));
386 if (keptrows) {
387 PetscInt cnt, *nrows, k;
388 const PetscInt *krows, *inact;
389 PetscInt rstart;
390
391 PetscCall(MatGetOwnershipRange(jac_inact_inact, &rstart, NULL));
392 PetscCall(MatDestroy(&jac_inact_inact));
393 PetscCall(ISDestroy(&IS_act));
394
395 PetscCall(ISGetLocalSize(keptrows, &cnt));
396 PetscCall(ISGetIndices(keptrows, &krows));
397 PetscCall(ISGetIndices(vi->IS_inact, &inact));
398 PetscCall(PetscMalloc1(cnt, &nrows));
399 for (k = 0; k < cnt; k++) nrows[k] = inact[krows[k] - rstart];
400 PetscCall(ISRestoreIndices(keptrows, &krows));
401 PetscCall(ISRestoreIndices(vi->IS_inact, &inact));
402 PetscCall(ISDestroy(&keptrows));
403 PetscCall(ISDestroy(&vi->IS_inact));
404
405 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)snes), cnt, nrows, PETSC_OWN_POINTER, &vi->IS_inact));
406 PetscCall(ISComplement(vi->IS_inact, F->map->rstart, F->map->rend, &IS_act));
407 PetscCall(MatCreateSubMatrix(snes->jacobian, vi->IS_inact, vi->IS_inact, MAT_INITIAL_MATRIX, &jac_inact_inact));
408 }
409 }
410 PetscCall(DMSetVI(snes->dm, vi->IS_inact));
411 /* remove later */
412
413 /*
414 PetscCall(VecView(vi->xu,PETSC_VIEWER_BINARY_(((PetscObject)vi->xu)->comm)));
415 PetscCall(VecView(vi->xl,PETSC_VIEWER_BINARY_(((PetscObject)vi->xl)->comm)));
416 PetscCall(VecView(X,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)X))));
417 PetscCall(VecView(F,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)F))));
418 PetscCall(ISView(vi->IS_inact,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)vi->IS_inact))));
419 */
420
421 /* Get sizes of active and inactive sets */
422 PetscCall(ISGetLocalSize(IS_act, &nis_act));
423 PetscCall(ISGetLocalSize(vi->IS_inact, &nis_inact));
424
425 /* Create active and inactive set vectors */
426 PetscCall(SNESCreateSubVectors_VINEWTONRSLS(snes, nis_inact, &F_inact));
427 PetscCall(SNESCreateSubVectors_VINEWTONRSLS(snes, nis_act, &Y_act));
428 PetscCall(SNESCreateSubVectors_VINEWTONRSLS(snes, nis_inact, &Y_inact));
429
430 /* Create scatter contexts */
431 PetscCall(VecScatterCreate(Y, IS_act, Y_act, NULL, &scat_act));
432 PetscCall(VecScatterCreate(Y, vi->IS_inact, Y_inact, NULL, &scat_inact));
433
434 /* Do a vec scatter to active and inactive set vectors */
435 PetscCall(VecScatterBegin(scat_inact, F, F_inact, INSERT_VALUES, SCATTER_FORWARD));
436 PetscCall(VecScatterEnd(scat_inact, F, F_inact, INSERT_VALUES, SCATTER_FORWARD));
437
438 PetscCall(VecScatterBegin(scat_act, Y, Y_act, INSERT_VALUES, SCATTER_FORWARD));
439 PetscCall(VecScatterEnd(scat_act, Y, Y_act, INSERT_VALUES, SCATTER_FORWARD));
440
441 PetscCall(VecScatterBegin(scat_inact, Y, Y_inact, INSERT_VALUES, SCATTER_FORWARD));
442 PetscCall(VecScatterEnd(scat_inact, Y, Y_inact, INSERT_VALUES, SCATTER_FORWARD));
443
444 /* Active set direction = 0 */
445 PetscCall(VecSet(Y_act, 0));
446 if (snes->jacobian != snes->jacobian_pre) PetscCall(MatCreateSubMatrix(snes->jacobian_pre, vi->IS_inact, vi->IS_inact, MAT_INITIAL_MATRIX, &prejac_inact_inact));
447 else prejac_inact_inact = jac_inact_inact;
448
449 PetscCall(ISEqual(vi->IS_inact_prev, vi->IS_inact, &isequal));
450 if (!isequal) {
451 PetscCall(SNESVIResetPCandKSP(snes, jac_inact_inact, prejac_inact_inact));
452 PetscCall(PCFieldSplitRestrictIS(pc, vi->IS_inact));
453 }
454
455 /* PetscCall(ISView(vi->IS_inact,0)); */
456 /* PetscCall(ISView(IS_act,0));*/
457 /* ierr = MatView(snes->jacobian_pre,0); */
458
459 PetscCall(KSPSetOperators(snes->ksp, jac_inact_inact, prejac_inact_inact));
460 PetscCall(KSPSetUp(snes->ksp));
461 {
462 PC pc;
463 PetscBool flg;
464 PetscCall(KSPGetPC(snes->ksp, &pc));
465 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &flg));
466 if (flg) {
467 KSP *subksps;
468 PetscCall(PCFieldSplitGetSubKSP(pc, NULL, &subksps));
469 PetscCall(KSPGetPC(subksps[0], &pc));
470 PetscCall(PetscFree(subksps));
471 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCBJACOBI, &flg));
472 if (flg) {
473 PetscInt n, N = 101 * 101, j, cnts[3] = {0, 0, 0};
474 const PetscInt *ii;
475
476 PetscCall(ISGetSize(vi->IS_inact, &n));
477 PetscCall(ISGetIndices(vi->IS_inact, &ii));
478 for (j = 0; j < n; j++) {
479 if (ii[j] < N) cnts[0]++;
480 else if (ii[j] < 2 * N) cnts[1]++;
481 else if (ii[j] < 3 * N) cnts[2]++;
482 }
483 PetscCall(ISRestoreIndices(vi->IS_inact, &ii));
484
485 PetscCall(PCBJacobiSetTotalBlocks(pc, 3, cnts));
486 }
487 }
488 }
489
490 PetscCall(KSPSolve(snes->ksp, F_inact, Y_inact));
491 PetscCall(VecScatterBegin(scat_act, Y_act, Y, INSERT_VALUES, SCATTER_REVERSE));
492 PetscCall(VecScatterEnd(scat_act, Y_act, Y, INSERT_VALUES, SCATTER_REVERSE));
493 PetscCall(VecScatterBegin(scat_inact, Y_inact, Y, INSERT_VALUES, SCATTER_REVERSE));
494 PetscCall(VecScatterEnd(scat_inact, Y_inact, Y, INSERT_VALUES, SCATTER_REVERSE));
495
496 PetscCall(VecDestroy(&F_inact));
497 PetscCall(VecDestroy(&Y_act));
498 PetscCall(VecDestroy(&Y_inact));
499 PetscCall(VecScatterDestroy(&scat_act));
500 PetscCall(VecScatterDestroy(&scat_inact));
501 PetscCall(ISDestroy(&IS_act));
502 if (!isequal) {
503 PetscCall(ISDestroy(&vi->IS_inact_prev));
504 PetscCall(ISDuplicate(vi->IS_inact, &vi->IS_inact_prev));
505 }
506 PetscCall(ISDestroy(&vi->IS_inact));
507 PetscCall(MatDestroy(&jac_inact_inact));
508 if (snes->jacobian != snes->jacobian_pre) PetscCall(MatDestroy(&prejac_inact_inact));
509
510 PetscCall(KSPGetConvergedReason(snes->ksp, &kspreason));
511 if (kspreason < 0) {
512 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
513 PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", number linear solve failures %" PetscInt_FMT " greater than current SNES allowed, stopping solve\n", snes->iter, snes->numLinearSolveFailures));
514 snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
515 break;
516 }
517 }
518
519 PetscCall(KSPGetIterationNumber(snes->ksp, &lits));
520 snes->linear_its += lits;
521 PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits));
522 /*
523 if (snes->ops->precheck) {
524 PetscBool changed_y = PETSC_FALSE;
525 PetscUseTypeMethod(snes,precheck ,X,Y,snes->precheck,&changed_y);
526 }
527
528 if (PetscLogPrintInfo) PetscCall(SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W));
529 */
530 /* Compute a (scaled) negative update in the line search routine:
531 Y <- X - lambda*Y
532 and evaluate G = function(Y) (depends on the line search).
533 */
534 PetscCall(VecCopy(Y, snes->vec_sol_update));
535 ynorm = 1;
536 gnorm = fnorm;
537 PetscCall(SNESLineSearchApply(snes->linesearch, X, F, &gnorm, Y));
538 PetscCall(DMDestroyVI(snes->dm));
539 if (snes->reason) break;
540 PetscCall(SNESLineSearchGetReason(snes->linesearch, &lsreason));
541 if (lsreason) {
542 if (snes->stol * xnorm > ynorm) {
543 snes->reason = SNES_CONVERGED_SNORM_RELATIVE;
544 break;
545 } else if (lsreason == SNES_LINESEARCH_FAILED_FUNCTION_DOMAIN) {
546 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
547 break;
548 } else if (lsreason == SNES_LINESEARCH_FAILED_NANORINF) {
549 snes->reason = SNES_DIVERGED_FUNCTION_NANORINF;
550 break;
551 } else if (lsreason == SNES_LINESEARCH_FAILED_OBJECTIVE_DOMAIN) {
552 snes->reason = SNES_DIVERGED_OBJECTIVE_DOMAIN;
553 break;
554 } else if (lsreason == SNES_LINESEARCH_FAILED_JACOBIAN_DOMAIN) {
555 snes->reason = SNES_DIVERGED_JACOBIAN_DOMAIN;
556 break;
557 } else if (++snes->numFailures >= snes->maxFailures) {
558 PetscBool ismin;
559
560 snes->reason = SNES_DIVERGED_LINE_SEARCH;
561 PetscCall(SNESVICheckLocalMin_Private(snes, snes->jacobian, F, X, gnorm, &ismin));
562 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
563 break;
564 }
565 }
566 PetscCall(SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm));
567 PetscCall(PetscInfo(snes, "fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n", (double)fnorm, (double)gnorm, (double)ynorm, (int)lsreason));
568 /* Update function and solution vectors */
569 fnorm = gnorm;
570 /* Monitor convergence */
571 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
572 snes->iter = i + 1;
573 snes->norm = fnorm;
574 snes->xnorm = xnorm;
575 snes->ynorm = ynorm;
576 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
577 PetscCall(SNESLogConvergenceHistory(snes, snes->norm, lits));
578 /* Test for convergence, xnorm = || X || */
579 if (snes->ops->converged != SNESConvergedSkip) PetscCall(VecNorm(X, NORM_2, &xnorm));
580 PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm));
581 PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
582 if (snes->reason) break;
583 }
584 /* make sure that the VI information attached to the DM is removed if the for loop above was broken early due to some exceptional conditional */
585 PetscCall(DMDestroyVI(snes->dm));
586 PetscFunctionReturn(PETSC_SUCCESS);
587 }
588
589 /*@C
590 SNESVISetRedundancyCheck - Provide a function to check for any redundancy in the VI active set
591
592 Logically Collective
593
594 Input Parameters:
595 + snes - the `SNESVINEWTONRSLS` context
596 . func - the function to check of redundancies
597 - ctx - optional context used by the function
598
599 Level: advanced
600
601 Note:
602 Sometimes the inactive set will result in a singular sub-Jacobian problem that needs to be solved, this allows the user,
603 when they know more about their specific problem to provide a function that removes the redundancy that results in the singular linear system
604
605 .seealso: [](ch_snes), `SNES`, `SNESVINEWTONRSLS`, `SNESVIGetInactiveSet()`, `DMSetVI()`
606 @*/
SNESVISetRedundancyCheck(SNES snes,PetscErrorCode (* func)(SNES,IS,IS *,void *),PetscCtx ctx)607 PetscErrorCode SNESVISetRedundancyCheck(SNES snes, PetscErrorCode (*func)(SNES, IS, IS *, void *), PetscCtx ctx)
608 {
609 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;
610
611 PetscFunctionBegin;
612 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
613 vi->checkredundancy = func;
614 vi->ctxP = ctx;
615 PetscFunctionReturn(PETSC_SUCCESS);
616 }
617
618 #if defined(PETSC_HAVE_MATLAB)
619 #include <engine.h>
620 #include <mex.h>
621 typedef struct {
622 char *funcname;
623 mxArray *ctx;
624 } SNESMatlabContext;
625
SNESVIRedundancyCheck_Matlab(SNES snes,IS is_act,IS * is_redact,PetscCtx ctx)626 PetscErrorCode SNESVIRedundancyCheck_Matlab(SNES snes, IS is_act, IS *is_redact, PetscCtx ctx)
627 {
628 SNESMatlabContext *sctx = (SNESMatlabContext *)ctx;
629 int nlhs = 1, nrhs = 5;
630 mxArray *plhs[1], *prhs[5];
631 long long int l1 = 0, l2 = 0, ls = 0;
632 PetscInt *indices = NULL;
633
634 PetscFunctionBegin;
635 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
636 PetscValidHeaderSpecific(is_act, IS_CLASSID, 2);
637 PetscAssertPointer(is_redact, 3);
638 PetscCheckSameComm(snes, 1, is_act, 2);
639
640 /* Create IS for reduced active set of size 0, its size and indices will
641 bet set by the MATLAB function */
642 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)snes), 0, indices, PETSC_OWN_POINTER, is_redact));
643 /* call MATLAB function in ctx */
644 PetscCall(PetscArraycpy(&ls, &snes, 1));
645 PetscCall(PetscArraycpy(&l1, &is_act, 1));
646 PetscCall(PetscArraycpy(&l2, is_redact, 1));
647 prhs[0] = mxCreateDoubleScalar((double)ls);
648 prhs[1] = mxCreateDoubleScalar((double)l1);
649 prhs[2] = mxCreateDoubleScalar((double)l2);
650 prhs[3] = mxCreateString(sctx->funcname);
651 prhs[4] = sctx->ctx;
652 PetscCall(mexCallMATLAB(nlhs, plhs, nrhs, prhs, "PetscSNESVIRedundancyCheckInternal"));
653 PetscCall(mxGetScalar(plhs[0]));
654 mxDestroyArray(prhs[0]);
655 mxDestroyArray(prhs[1]);
656 mxDestroyArray(prhs[2]);
657 mxDestroyArray(prhs[3]);
658 mxDestroyArray(plhs[0]);
659 PetscFunctionReturn(PETSC_SUCCESS);
660 }
661
SNESVISetRedundancyCheckMatlab(SNES snes,const char * func,mxArray * ctx)662 PetscErrorCode SNESVISetRedundancyCheckMatlab(SNES snes, const char *func, mxArray *ctx)
663 {
664 SNESMatlabContext *sctx;
665
666 PetscFunctionBegin;
667 /* currently sctx is memory bleed */
668 PetscCall(PetscNew(&sctx));
669 PetscCall(PetscStrallocpy(func, &sctx->funcname));
670 sctx->ctx = mxDuplicateArray(ctx);
671 PetscCall(SNESVISetRedundancyCheck(snes, SNESVIRedundancyCheck_Matlab, sctx));
672 PetscFunctionReturn(PETSC_SUCCESS);
673 }
674
675 #endif
676
SNESSetUp_VINEWTONRSLS(SNES snes)677 static PetscErrorCode SNESSetUp_VINEWTONRSLS(SNES snes)
678 {
679 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;
680 PetscInt *indices;
681 PetscInt i, n, rstart, rend;
682 SNESLineSearch linesearch;
683
684 PetscFunctionBegin;
685 PetscCall(SNESSetUp_VI(snes));
686
687 /* Set up previous active index set for the first snes solve
688 vi->IS_inact_prev = 0,1,2,....N */
689
690 PetscCall(VecGetOwnershipRange(snes->work[0], &rstart, &rend));
691 PetscCall(VecGetLocalSize(snes->work[0], &n));
692 PetscCall(PetscMalloc1(n, &indices));
693 for (i = 0; i < n; i++) indices[i] = rstart + i;
694 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)snes), n, indices, PETSC_OWN_POINTER, &vi->IS_inact_prev));
695
696 /* set the line search functions */
697 if (!snes->linesearch) {
698 PetscCall(SNESGetLineSearch(snes, &linesearch));
699 PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT));
700 }
701 PetscFunctionReturn(PETSC_SUCCESS);
702 }
703
SNESReset_VINEWTONRSLS(SNES snes)704 static PetscErrorCode SNESReset_VINEWTONRSLS(SNES snes)
705 {
706 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;
707
708 PetscFunctionBegin;
709 PetscCall(SNESReset_VI(snes));
710 PetscCall(ISDestroy(&vi->IS_inact_prev));
711 PetscFunctionReturn(PETSC_SUCCESS);
712 }
713
714 /*MC
715 SNESVINEWTONRSLS - Reduced space active set solvers for variational inequalities based on Newton's method
716
717 Options Database Keys:
718 + -snes_type <vinewtonssls,vinewtonrsls> - a semi-smooth solver or a reduced space active set method
719 - -snes_vi_monitor - prints the number of active constraints at each iteration.
720
721 Level: beginner
722
723 Note:
724 At each set of this methods the algorithm produces an inactive set of variables that are constrained to their current values
725 (because changing these values would result in those variables no longer satisfying the inequality constraints)
726 and produces a step direction by solving the linear system arising from the Jacobian with the inactive variables removed. In other
727 words on a reduced space of the solution space. Based on the Newton update it then adjusts the inactive set for the next iteration.
728
729 See {cite}`benson2006flexible`
730
731 .seealso: [](ch_snes), `SNESVISetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESVINEWTONSSLS`, `SNESNEWTONTR`, `SNESLineSearchSetType()`, `SNESLineSearchSetPostCheck()`, `SNESLineSearchSetPreCheck()`, `SNESVIGetInactiveSet()`, `DMSetVI()`, `SNESVISetRedundancyCheck()`
732 M*/
SNESCreate_VINEWTONRSLS(SNES snes)733 PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONRSLS(SNES snes)
734 {
735 SNES_VINEWTONRSLS *vi;
736 SNESLineSearch linesearch;
737
738 PetscFunctionBegin;
739 snes->ops->reset = SNESReset_VINEWTONRSLS;
740 snes->ops->setup = SNESSetUp_VINEWTONRSLS;
741 snes->ops->solve = SNESSolve_VINEWTONRSLS;
742 snes->ops->destroy = SNESDestroy_VI;
743 snes->ops->setfromoptions = SNESSetFromOptions_VI;
744 snes->ops->view = NULL;
745 snes->ops->converged = SNESConvergedDefault_VI;
746
747 snes->usesksp = PETSC_TRUE;
748 snes->usesnpc = PETSC_FALSE;
749
750 PetscCall(SNESGetLineSearch(snes, &linesearch));
751 if (!((PetscObject)linesearch)->type_name) PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT));
752 PetscCall(SNESLineSearchBTSetAlpha(linesearch, 0.0));
753
754 snes->alwayscomputesfinalresidual = PETSC_TRUE;
755
756 PetscCall(SNESParametersInitialize(snes));
757
758 PetscCall(PetscNew(&vi));
759 snes->data = (void *)vi;
760 vi->checkredundancy = NULL;
761
762 PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESVISetVariableBounds_C", SNESVISetVariableBounds_VI));
763 PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESVISetComputeVariableBounds_C", SNESVISetComputeVariableBounds_VI));
764 PetscFunctionReturn(PETSC_SUCCESS);
765 }
766