Actual source code: pcis.c
petsc-3.13.6 2020-09-29
2: #include <../src/ksp/pc/impls/is/pcis.h>
4: static PetscErrorCode PCISSetUseStiffnessScaling_IS(PC pc, PetscBool use)
5: {
6: PC_IS *pcis = (PC_IS*)pc->data;
9: pcis->use_stiffness_scaling = use;
10: return(0);
11: }
13: /*@
14: PCISSetUseStiffnessScaling - Tells PCIS to construct partition of unity using
15: local matrices' diagonal.
17: Not collective
19: Input Parameters:
20: + pc - the preconditioning context
21: - use - whether or not pcis use matrix diagonal to build partition of unity.
23: Level: intermediate
25: Notes:
27: .seealso: PCBDDC
28: @*/
29: PetscErrorCode PCISSetUseStiffnessScaling(PC pc, PetscBool use)
30: {
36: PetscTryMethod(pc,"PCISSetUseStiffnessScaling_C",(PC,PetscBool),(pc,use));
37: return(0);
38: }
40: static PetscErrorCode PCISSetSubdomainDiagonalScaling_IS(PC pc, Vec scaling_factors)
41: {
43: PC_IS *pcis = (PC_IS*)pc->data;
46: PetscObjectReference((PetscObject)scaling_factors);
47: VecDestroy(&pcis->D);
48: pcis->D = scaling_factors;
49: if (pc->setupcalled) {
50: PetscInt sn;
52: VecGetSize(pcis->D,&sn);
53: if (sn == pcis->n) {
54: VecScatterBegin(pcis->N_to_B,pcis->D,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
55: VecScatterEnd(pcis->N_to_B,pcis->D,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
56: VecDestroy(&pcis->D);
57: VecDuplicate(pcis->vec1_B,&pcis->D);
58: VecCopy(pcis->vec1_B,pcis->D);
59: } else if (sn != pcis->n_B) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Invalid size for scaling vector. Expected %D (or full %D), found %D",pcis->n_B,pcis->n,sn);
60: }
61: return(0);
62: }
64: /*@
65: PCISSetSubdomainDiagonalScaling - Set diagonal scaling for PCIS.
67: Not collective
69: Input Parameters:
70: + pc - the preconditioning context
71: - scaling_factors - scaling factors for the subdomain
73: Level: intermediate
75: Notes:
76: Intended to use with jumping coefficients cases.
78: .seealso: PCBDDC
79: @*/
80: PetscErrorCode PCISSetSubdomainDiagonalScaling(PC pc, Vec scaling_factors)
81: {
87: PetscTryMethod(pc,"PCISSetSubdomainDiagonalScaling_C",(PC,Vec),(pc,scaling_factors));
88: return(0);
89: }
91: static PetscErrorCode PCISSetSubdomainScalingFactor_IS(PC pc, PetscScalar scal)
92: {
93: PC_IS *pcis = (PC_IS*)pc->data;
96: pcis->scaling_factor = scal;
97: if (pcis->D) {
100: VecSet(pcis->D,pcis->scaling_factor);
101: }
102: return(0);
103: }
105: /*@
106: PCISSetSubdomainScalingFactor - Set scaling factor for PCIS.
108: Not collective
110: Input Parameters:
111: + pc - the preconditioning context
112: - scal - scaling factor for the subdomain
114: Level: intermediate
116: Notes:
117: Intended to use with jumping coefficients cases.
119: .seealso: PCBDDC
120: @*/
121: PetscErrorCode PCISSetSubdomainScalingFactor(PC pc, PetscScalar scal)
122: {
127: PetscTryMethod(pc,"PCISSetSubdomainScalingFactor_C",(PC,PetscScalar),(pc,scal));
128: return(0);
129: }
132: /* -------------------------------------------------------------------------- */
133: /*
134: PCISSetUp -
135: */
136: PetscErrorCode PCISSetUp(PC pc, PetscBool computematrices, PetscBool computesolvers)
137: {
138: PC_IS *pcis = (PC_IS*)(pc->data);
139: Mat_IS *matis;
140: MatReuse reuse;
142: PetscBool flg,issbaij;
145: PetscObjectTypeCompare((PetscObject)pc->pmat,MATIS,&flg);
146: if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Requires preconditioning matrix of type MATIS");
147: matis = (Mat_IS*)pc->pmat->data;
148: if (pc->useAmat) {
149: PetscObjectTypeCompare((PetscObject)pc->mat,MATIS,&flg);
150: if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Requires linear system matrix of type MATIS");
151: }
153: /* first time creation, get info on substructuring */
154: if (!pc->setupcalled) {
155: PetscInt n_I;
156: PetscInt *idx_I_local,*idx_B_local,*idx_I_global,*idx_B_global;
157: PetscBT bt;
158: PetscInt i,j;
160: /* get info on mapping */
161: PetscObjectReference((PetscObject)pc->pmat->rmap->mapping);
162: ISLocalToGlobalMappingDestroy(&pcis->mapping);
163: pcis->mapping = pc->pmat->rmap->mapping;
164: ISLocalToGlobalMappingGetSize(pcis->mapping,&pcis->n);
165: ISLocalToGlobalMappingGetInfo(pcis->mapping,&(pcis->n_neigh),&(pcis->neigh),&(pcis->n_shared),&(pcis->shared));
167: /* Identifying interior and interface nodes, in local numbering */
168: PetscBTCreate(pcis->n,&bt);
169: for (i=0;i<pcis->n_neigh;i++)
170: for (j=0;j<pcis->n_shared[i];j++) {
171: PetscBTSet(bt,pcis->shared[i][j]);
172: }
174: /* Creating local and global index sets for interior and inteface nodes. */
175: PetscMalloc1(pcis->n,&idx_I_local);
176: PetscMalloc1(pcis->n,&idx_B_local);
177: for (i=0, pcis->n_B=0, n_I=0; i<pcis->n; i++) {
178: if (!PetscBTLookup(bt,i)) {
179: idx_I_local[n_I] = i;
180: n_I++;
181: } else {
182: idx_B_local[pcis->n_B] = i;
183: pcis->n_B++;
184: }
185: }
187: /* Getting the global numbering */
188: idx_B_global = idx_I_local + n_I; /* Just avoiding allocating extra memory, since we have vacant space */
189: idx_I_global = idx_B_local + pcis->n_B;
190: ISLocalToGlobalMappingApply(pcis->mapping,pcis->n_B,idx_B_local,idx_B_global);
191: ISLocalToGlobalMappingApply(pcis->mapping,n_I,idx_I_local,idx_I_global);
193: /* Creating the index sets */
194: ISCreateGeneral(PETSC_COMM_SELF,pcis->n_B,idx_B_local,PETSC_COPY_VALUES, &pcis->is_B_local);
195: ISCreateGeneral(PetscObjectComm((PetscObject)pc),pcis->n_B,idx_B_global,PETSC_COPY_VALUES,&pcis->is_B_global);
196: ISCreateGeneral(PETSC_COMM_SELF,n_I,idx_I_local,PETSC_COPY_VALUES, &pcis->is_I_local);
197: ISCreateGeneral(PetscObjectComm((PetscObject)pc),n_I,idx_I_global,PETSC_COPY_VALUES,&pcis->is_I_global);
199: /* Freeing memory */
200: PetscFree(idx_B_local);
201: PetscFree(idx_I_local);
202: PetscBTDestroy(&bt);
204: /* Creating work vectors and arrays */
205: VecDuplicate(matis->x,&pcis->vec1_N);
206: VecDuplicate(pcis->vec1_N,&pcis->vec2_N);
207: VecCreate(PETSC_COMM_SELF,&pcis->vec1_D);
208: VecSetSizes(pcis->vec1_D,pcis->n-pcis->n_B,PETSC_DECIDE);
209: VecSetType(pcis->vec1_D,((PetscObject)pcis->vec1_N)->type_name);
210: VecDuplicate(pcis->vec1_D,&pcis->vec2_D);
211: VecDuplicate(pcis->vec1_D,&pcis->vec3_D);
212: VecDuplicate(pcis->vec1_D,&pcis->vec4_D);
213: VecCreate(PETSC_COMM_SELF,&pcis->vec1_B);
214: VecSetSizes(pcis->vec1_B,pcis->n_B,PETSC_DECIDE);
215: VecSetType(pcis->vec1_B,((PetscObject)pcis->vec1_N)->type_name);
216: VecDuplicate(pcis->vec1_B,&pcis->vec2_B);
217: VecDuplicate(pcis->vec1_B,&pcis->vec3_B);
218: MatCreateVecs(pc->pmat,&pcis->vec1_global,0);
219: PetscMalloc1(pcis->n,&pcis->work_N);
220: /* scaling vector */
221: if (!pcis->D) { /* it can happen that the user passed in a scaling vector via PCISSetSubdomainDiagonalScaling */
222: VecDuplicate(pcis->vec1_B,&pcis->D);
223: VecSet(pcis->D,pcis->scaling_factor);
224: }
226: /* Creating the scatter contexts */
227: VecScatterCreate(pcis->vec1_N,pcis->is_I_local,pcis->vec1_D,(IS)0,&pcis->N_to_D);
228: VecScatterCreate(pcis->vec1_global,pcis->is_I_global,pcis->vec1_D,(IS)0,&pcis->global_to_D);
229: VecScatterCreate(pcis->vec1_N,pcis->is_B_local,pcis->vec1_B,(IS)0,&pcis->N_to_B);
230: VecScatterCreate(pcis->vec1_global,pcis->is_B_global,pcis->vec1_B,(IS)0,&pcis->global_to_B);
232: /* map from boundary to local */
233: ISLocalToGlobalMappingCreateIS(pcis->is_B_local,&pcis->BtoNmap);
234: }
236: {
237: PetscInt sn;
239: VecGetSize(pcis->D,&sn);
240: if (sn == pcis->n) {
241: VecScatterBegin(pcis->N_to_B,pcis->D,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
242: VecScatterEnd(pcis->N_to_B,pcis->D,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
243: VecDestroy(&pcis->D);
244: VecDuplicate(pcis->vec1_B,&pcis->D);
245: VecCopy(pcis->vec1_B,pcis->D);
246: } else if (sn != pcis->n_B) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Invalid size for scaling vector. Expected %D (or full %D), found %D",pcis->n_B,pcis->n,sn);
247: }
249: /*
250: Extracting the blocks A_II, A_BI, A_IB and A_BB from A. If the numbering
251: is such that interior nodes come first than the interface ones, we have
253: [ A_II | A_IB ]
254: A = [------+------]
255: [ A_BI | A_BB ]
256: */
257: if (computematrices) {
258: PetscBool amat = (PetscBool)(pc->mat != pc->pmat && pc->useAmat);
259: PetscInt bs,ibs;
261: reuse = MAT_INITIAL_MATRIX;
262: if (pcis->reusesubmatrices && pc->setupcalled) {
263: if (pc->flag == SAME_NONZERO_PATTERN) {
264: reuse = MAT_REUSE_MATRIX;
265: } else {
266: reuse = MAT_INITIAL_MATRIX;
267: }
268: }
269: if (reuse == MAT_INITIAL_MATRIX) {
270: MatDestroy(&pcis->A_II);
271: MatDestroy(&pcis->pA_II);
272: MatDestroy(&pcis->A_IB);
273: MatDestroy(&pcis->A_BI);
274: MatDestroy(&pcis->A_BB);
275: }
277: ISLocalToGlobalMappingGetBlockSize(pcis->mapping,&ibs);
278: MatGetBlockSize(matis->A,&bs);
279: MatCreateSubMatrix(matis->A,pcis->is_I_local,pcis->is_I_local,reuse,&pcis->pA_II);
280: if (amat) {
281: Mat_IS *amatis = (Mat_IS*)pc->mat->data;
282: MatCreateSubMatrix(amatis->A,pcis->is_I_local,pcis->is_I_local,reuse,&pcis->A_II);
283: } else {
284: PetscObjectReference((PetscObject)pcis->pA_II);
285: MatDestroy(&pcis->A_II);
286: pcis->A_II = pcis->pA_II;
287: }
288: MatSetBlockSize(pcis->A_II,bs == ibs ? bs : 1);
289: MatSetBlockSize(pcis->pA_II,bs == ibs ? bs : 1);
290: MatCreateSubMatrix(matis->A,pcis->is_B_local,pcis->is_B_local,reuse,&pcis->A_BB);
291: PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);
292: if (!issbaij) {
293: MatCreateSubMatrix(matis->A,pcis->is_I_local,pcis->is_B_local,reuse,&pcis->A_IB);
294: MatCreateSubMatrix(matis->A,pcis->is_B_local,pcis->is_I_local,reuse,&pcis->A_BI);
295: } else {
296: Mat newmat;
298: MatConvert(matis->A,MATSEQBAIJ,MAT_INITIAL_MATRIX,&newmat);
299: MatCreateSubMatrix(newmat,pcis->is_I_local,pcis->is_B_local,reuse,&pcis->A_IB);
300: MatCreateSubMatrix(newmat,pcis->is_B_local,pcis->is_I_local,reuse,&pcis->A_BI);
301: MatDestroy(&newmat);
302: }
303: MatSetBlockSize(pcis->A_BB,bs == ibs ? bs : 1);
304: }
306: /* Creating scaling vector D */
307: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_is_use_stiffness_scaling",&pcis->use_stiffness_scaling,NULL);
308: if (pcis->use_stiffness_scaling) {
309: PetscScalar *a;
310: PetscInt i,n;
312: if (pcis->A_BB) {
313: MatGetDiagonal(pcis->A_BB,pcis->D);
314: } else {
315: MatGetDiagonal(matis->A,pcis->vec1_N);
316: VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);
317: VecScatterEnd(pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);
318: }
319: VecAbs(pcis->D);
320: VecGetLocalSize(pcis->D,&n);
321: VecGetArray(pcis->D,&a);
322: for (i=0;i<n;i++) if (PetscAbsScalar(a[i])<PETSC_SMALL) a[i] = 1.0;
323: VecRestoreArray(pcis->D,&a);
324: }
325: VecSet(pcis->vec1_global,0.0);
326: VecScatterBegin(pcis->global_to_B,pcis->D,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
327: VecScatterEnd(pcis->global_to_B,pcis->D,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);
328: VecScatterBegin(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
329: VecScatterEnd(pcis->global_to_B,pcis->vec1_global,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
330: VecPointwiseDivide(pcis->D,pcis->D,pcis->vec1_B);
331: /* See historical note 01, at the bottom of this file. */
333: /* Creating the KSP contexts for the local Dirichlet and Neumann problems */
334: if (computesolvers) {
335: PC pc_ctx;
337: pcis->pure_neumann = matis->pure_neumann;
338: /* Dirichlet */
339: KSPCreate(PETSC_COMM_SELF,&pcis->ksp_D);
340: KSPSetErrorIfNotConverged(pcis->ksp_D,pc->erroriffailure);
341: PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_D,(PetscObject)pc,1);
342: KSPSetOperators(pcis->ksp_D,pcis->A_II,pcis->A_II);
343: KSPSetOptionsPrefix(pcis->ksp_D,"is_localD_");
344: KSPGetPC(pcis->ksp_D,&pc_ctx);
345: PCSetType(pc_ctx,PCLU);
346: KSPSetType(pcis->ksp_D,KSPPREONLY);
347: KSPSetFromOptions(pcis->ksp_D);
348: /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
349: KSPSetUp(pcis->ksp_D);
350: /* Neumann */
351: KSPCreate(PETSC_COMM_SELF,&pcis->ksp_N);
352: KSPSetErrorIfNotConverged(pcis->ksp_N,pc->erroriffailure);
353: PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_N,(PetscObject)pc,1);
354: KSPSetOperators(pcis->ksp_N,matis->A,matis->A);
355: KSPSetOptionsPrefix(pcis->ksp_N,"is_localN_");
356: KSPGetPC(pcis->ksp_N,&pc_ctx);
357: PCSetType(pc_ctx,PCLU);
358: KSPSetType(pcis->ksp_N,KSPPREONLY);
359: KSPSetFromOptions(pcis->ksp_N);
360: {
361: PetscBool damp_fixed = PETSC_FALSE,
362: remove_nullspace_fixed = PETSC_FALSE,
363: set_damping_factor_floating = PETSC_FALSE,
364: not_damp_floating = PETSC_FALSE,
365: not_remove_nullspace_floating = PETSC_FALSE;
366: PetscReal fixed_factor,
367: floating_factor;
369: PetscOptionsGetReal(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_damp_fixed",&fixed_factor,&damp_fixed);
370: if (!damp_fixed) fixed_factor = 0.0;
371: PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_damp_fixed",&damp_fixed,NULL);
373: PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_remove_nullspace_fixed",&remove_nullspace_fixed,NULL);
375: PetscOptionsGetReal(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_set_damping_factor_floating",
376: &floating_factor,&set_damping_factor_floating);
377: if (!set_damping_factor_floating) floating_factor = 0.0;
378: PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_set_damping_factor_floating",&set_damping_factor_floating,NULL);
379: if (!set_damping_factor_floating) floating_factor = 1.e-12;
381: PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_not_damp_floating",¬_damp_floating,NULL);
383: PetscOptionsGetBool(((PetscObject)pc_ctx)->options,((PetscObject)pc_ctx)->prefix,"-pc_is_not_remove_nullspace_floating",¬_remove_nullspace_floating,NULL);
385: if (pcis->pure_neumann) { /* floating subdomain */
386: if (!(not_damp_floating)) {
387: PCFactorSetShiftType(pc_ctx,MAT_SHIFT_NONZERO);
388: PCFactorSetShiftAmount(pc_ctx,floating_factor);
389: }
390: if (!(not_remove_nullspace_floating)) {
391: MatNullSpace nullsp;
392: MatNullSpaceCreate(PETSC_COMM_SELF,PETSC_TRUE,0,NULL,&nullsp);
393: MatSetNullSpace(matis->A,nullsp);
394: MatNullSpaceDestroy(&nullsp);
395: }
396: } else { /* fixed subdomain */
397: if (damp_fixed) {
398: PCFactorSetShiftType(pc_ctx,MAT_SHIFT_NONZERO);
399: PCFactorSetShiftAmount(pc_ctx,floating_factor);
400: }
401: if (remove_nullspace_fixed) {
402: MatNullSpace nullsp;
403: MatNullSpaceCreate(PETSC_COMM_SELF,PETSC_TRUE,0,NULL,&nullsp);
404: MatSetNullSpace(matis->A,nullsp);
405: MatNullSpaceDestroy(&nullsp);
406: }
407: }
408: }
409: /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
410: KSPSetUp(pcis->ksp_N);
411: }
412: return(0);
413: }
415: /* -------------------------------------------------------------------------- */
416: /*
417: PCISDestroy -
418: */
419: PetscErrorCode PCISDestroy(PC pc)
420: {
421: PC_IS *pcis = (PC_IS*)(pc->data);
425: ISDestroy(&pcis->is_B_local);
426: ISDestroy(&pcis->is_I_local);
427: ISDestroy(&pcis->is_B_global);
428: ISDestroy(&pcis->is_I_global);
429: MatDestroy(&pcis->A_II);
430: MatDestroy(&pcis->pA_II);
431: MatDestroy(&pcis->A_IB);
432: MatDestroy(&pcis->A_BI);
433: MatDestroy(&pcis->A_BB);
434: VecDestroy(&pcis->D);
435: KSPDestroy(&pcis->ksp_N);
436: KSPDestroy(&pcis->ksp_D);
437: VecDestroy(&pcis->vec1_N);
438: VecDestroy(&pcis->vec2_N);
439: VecDestroy(&pcis->vec1_D);
440: VecDestroy(&pcis->vec2_D);
441: VecDestroy(&pcis->vec3_D);
442: VecDestroy(&pcis->vec4_D);
443: VecDestroy(&pcis->vec1_B);
444: VecDestroy(&pcis->vec2_B);
445: VecDestroy(&pcis->vec3_B);
446: VecDestroy(&pcis->vec1_global);
447: VecScatterDestroy(&pcis->global_to_D);
448: VecScatterDestroy(&pcis->N_to_B);
449: VecScatterDestroy(&pcis->N_to_D);
450: VecScatterDestroy(&pcis->global_to_B);
451: PetscFree(pcis->work_N);
452: if (pcis->n_neigh > -1) {
453: ISLocalToGlobalMappingRestoreInfo(pcis->mapping,&(pcis->n_neigh),&(pcis->neigh),&(pcis->n_shared),&(pcis->shared));
454: }
455: ISLocalToGlobalMappingDestroy(&pcis->mapping);
456: ISLocalToGlobalMappingDestroy(&pcis->BtoNmap);
457: PetscObjectComposeFunction((PetscObject)pc,"PCISSetUseStiffnessScaling_C",NULL);
458: PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainScalingFactor_C",NULL);
459: PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainDiagonalScaling_C",NULL);
460: return(0);
461: }
463: /* -------------------------------------------------------------------------- */
464: /*
465: PCISCreate -
466: */
467: PetscErrorCode PCISCreate(PC pc)
468: {
469: PC_IS *pcis = (PC_IS*)(pc->data);
473: pcis->n_neigh = -1;
474: pcis->scaling_factor = 1.0;
475: pcis->reusesubmatrices = PETSC_TRUE;
476: /* composing functions */
477: PetscObjectComposeFunction((PetscObject)pc,"PCISSetUseStiffnessScaling_C",PCISSetUseStiffnessScaling_IS);
478: PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainScalingFactor_C",PCISSetSubdomainScalingFactor_IS);
479: PetscObjectComposeFunction((PetscObject)pc,"PCISSetSubdomainDiagonalScaling_C",PCISSetSubdomainDiagonalScaling_IS);
480: return(0);
481: }
483: /* -------------------------------------------------------------------------- */
484: /*
485: PCISApplySchur -
487: Input parameters:
488: . pc - preconditioner context
489: . v - vector to which the Schur complement is to be applied (it is NOT modified inside this function, UNLESS vec2_B is null)
491: Output parameters:
492: . vec1_B - result of Schur complement applied to chunk
493: . vec2_B - garbage (used as work space), or null (and v is used as workspace)
494: . vec1_D - garbage (used as work space)
495: . vec2_D - garbage (used as work space)
497: */
498: PetscErrorCode PCISApplySchur(PC pc, Vec v, Vec vec1_B, Vec vec2_B, Vec vec1_D, Vec vec2_D)
499: {
501: PC_IS *pcis = (PC_IS*)(pc->data);
504: if (!vec2_B) vec2_B = v;
506: MatMult(pcis->A_BB,v,vec1_B);
507: MatMult(pcis->A_IB,v,vec1_D);
508: KSPSolve(pcis->ksp_D,vec1_D,vec2_D);
509: KSPCheckSolve(pcis->ksp_D,pc,vec2_D);
510: MatMult(pcis->A_BI,vec2_D,vec2_B);
511: VecAXPY(vec1_B,-1.0,vec2_B);
512: return(0);
513: }
515: /* -------------------------------------------------------------------------- */
516: /*
517: PCISScatterArrayNToVecB - Scatters interface node values from a big array (of all local nodes, interior or interface,
518: including ghosts) into an interface vector, when in SCATTER_FORWARD mode, or vice-versa, when in SCATTER_REVERSE
519: mode.
521: Input parameters:
522: . pc - preconditioner context
523: . array_N - [when in SCATTER_FORWARD mode] Array to be scattered into the vector
524: . v_B - [when in SCATTER_REVERSE mode] Vector to be scattered into the array
526: Output parameter:
527: . array_N - [when in SCATTER_REVERSE mode] Array to receive the scattered vector
528: . v_B - [when in SCATTER_FORWARD mode] Vector to receive the scattered array
530: Notes:
531: The entries in the array that do not correspond to interface nodes remain unaltered.
532: */
533: PetscErrorCode PCISScatterArrayNToVecB(PetscScalar *array_N, Vec v_B, InsertMode imode, ScatterMode smode, PC pc)
534: {
535: PetscInt i;
536: const PetscInt *idex;
538: PetscScalar *array_B;
539: PC_IS *pcis = (PC_IS*)(pc->data);
542: VecGetArray(v_B,&array_B);
543: ISGetIndices(pcis->is_B_local,&idex);
545: if (smode == SCATTER_FORWARD) {
546: if (imode == INSERT_VALUES) {
547: for (i=0; i<pcis->n_B; i++) array_B[i] = array_N[idex[i]];
548: } else { /* ADD_VALUES */
549: for (i=0; i<pcis->n_B; i++) array_B[i] += array_N[idex[i]];
550: }
551: } else { /* SCATTER_REVERSE */
552: if (imode == INSERT_VALUES) {
553: for (i=0; i<pcis->n_B; i++) array_N[idex[i]] = array_B[i];
554: } else { /* ADD_VALUES */
555: for (i=0; i<pcis->n_B; i++) array_N[idex[i]] += array_B[i];
556: }
557: }
558: ISRestoreIndices(pcis->is_B_local,&idex);
559: VecRestoreArray(v_B,&array_B);
560: return(0);
561: }
563: /* -------------------------------------------------------------------------- */
564: /*
565: PCISApplyInvSchur - Solves the Neumann problem related to applying the inverse of the Schur complement.
566: More precisely, solves the problem:
567: [ A_II A_IB ] [ . ] [ 0 ]
568: [ ] [ ] = [ ]
569: [ A_BI A_BB ] [ x ] [ b ]
571: Input parameters:
572: . pc - preconditioner context
573: . b - vector of local interface nodes (including ghosts)
575: Output parameters:
576: . x - vector of local interface nodes (including ghosts); returns the Section 1.5 Writing Application Codes with PETSc of the inverse of the Schur
577: complement to b
578: . vec1_N - vector of local nodes (interior and interface, including ghosts); returns garbage (used as work space)
579: . vec2_N - vector of local nodes (interior and interface, including ghosts); returns garbage (used as work space)
581: */
582: PetscErrorCode PCISApplyInvSchur(PC pc, Vec b, Vec x, Vec vec1_N, Vec vec2_N)
583: {
585: PC_IS *pcis = (PC_IS*)(pc->data);
588: /*
589: Neumann solvers.
590: Applying the inverse of the local Schur complement, i.e, solving a Neumann
591: Problem with zero at the interior nodes of the RHS and extracting the interface
592: part of the solution. inverse Schur complement is applied to b and the result
593: is stored in x.
594: */
595: /* Setting the RHS vec1_N */
596: VecSet(vec1_N,0.0);
597: VecScatterBegin(pcis->N_to_B,b,vec1_N,INSERT_VALUES,SCATTER_REVERSE);
598: VecScatterEnd (pcis->N_to_B,b,vec1_N,INSERT_VALUES,SCATTER_REVERSE);
599: /* Checking for consistency of the RHS */
600: {
601: PetscBool flg = PETSC_FALSE;
602: PetscOptionsGetBool(NULL,NULL,"-pc_is_check_consistency",&flg,NULL);
603: if (flg) {
604: PetscScalar average;
605: PetscViewer viewer;
606: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc),&viewer);
608: VecSum(vec1_N,&average);
609: average = average / ((PetscReal)pcis->n);
610: PetscViewerASCIIPushSynchronized(viewer);
611: if (pcis->pure_neumann) {
612: PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d is floating. Average = % 1.14e\n",PetscGlobalRank,PetscAbsScalar(average));
613: } else {
614: PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d is fixed. Average = % 1.14e\n",PetscGlobalRank,PetscAbsScalar(average));
615: }
616: PetscViewerFlush(viewer);
617: PetscViewerASCIIPopSynchronized(viewer);
618: }
619: }
620: /* Solving the system for vec2_N */
621: KSPSolve(pcis->ksp_N,vec1_N,vec2_N);
622: KSPCheckSolve(pcis->ksp_N,pc,vec2_N);
623: /* Extracting the local interface vector out of the solution */
624: VecScatterBegin(pcis->N_to_B,vec2_N,x,INSERT_VALUES,SCATTER_FORWARD);
625: VecScatterEnd (pcis->N_to_B,vec2_N,x,INSERT_VALUES,SCATTER_FORWARD);
626: return(0);
627: }