Actual source code: pcis.c
petsc-3.3-p7 2013-05-11
2: #include pcis.h
4: EXTERN_C_BEGIN
7: static PetscErrorCode PCISSetSubdomainScalingFactor_IS(PC pc, PetscScalar scal)
8: {
9: PC_IS *pcis = (PC_IS*)pc->data;
12: pcis->scaling_factor = scal;
13: return(0);
14: }
15: EXTERN_C_END
19: /*@
20: PCISSetSubdomainScalingFactor - Set scaling factor for PCIS.
22: Not collective
24: Input Parameters:
25: + pc - the preconditioning context
26: - scal - scaling factor for the subdomain
28: Level: intermediate
30: Notes:
31: Intended to use with jumping coefficients cases.
33: .seealso: PCBDDC
34: @*/
35: PetscErrorCode PCISSetSubdomainScalingFactor(PC pc, PetscScalar scal)
36: {
41: PetscTryMethod(pc,"PCISSetSubdomainScalingFactor_C",(PC,PetscScalar),(pc,scal));
42: return(0);
43: }
46: /* -------------------------------------------------------------------------- */
47: /*
48: PCISSetUp -
49: */
52: PetscErrorCode PCISSetUp(PC pc)
53: {
54: PC_IS *pcis = (PC_IS*)(pc->data);
55: Mat_IS *matis = (Mat_IS*)pc->mat->data;
56: PetscInt i;
57: PetscErrorCode ierr;
58: PetscBool flg;
59: Vec counter;
60:
62: PetscObjectTypeCompare((PetscObject)pc->mat,MATIS,&flg);
63: if (!flg) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG,"Preconditioner type of Neumann Neumman requires matrix of type MATIS");
65: pcis->pure_neumann = matis->pure_neumann;
67: /*
68: Creating the local vector vec1_N, containing the inverse of the number
69: of subdomains to which each local node (either owned or ghost)
70: pertains. To accomplish that, we scatter local vectors of 1's to
71: a global vector (adding the values); scatter the result back to
72: local vectors and finally invert the result.
73: */
74: VecDuplicate(matis->x,&pcis->vec1_N);
75: MatGetVecs(pc->pmat,&counter,0); /* temporary auxiliar vector */
76: VecSet(counter,0.0);
77: VecSet(pcis->vec1_N,1.0);
78: VecScatterBegin(matis->ctx,pcis->vec1_N,counter,ADD_VALUES,SCATTER_REVERSE);
79: VecScatterEnd (matis->ctx,pcis->vec1_N,counter,ADD_VALUES,SCATTER_REVERSE);
80: VecScatterBegin(matis->ctx,counter,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);
81: VecScatterEnd (matis->ctx,counter,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);
82: /*
83: Creating local and global index sets for interior and
84: inteface nodes. Notice that interior nodes have D[i]==1.0.
85: */
86: {
87: PetscInt n_I;
88: PetscInt *idx_I_local,*idx_B_local,*idx_I_global,*idx_B_global;
89: PetscScalar *array;
90: /* Identifying interior and interface nodes, in local numbering */
91: VecGetSize(pcis->vec1_N,&pcis->n);
92: VecGetArray(pcis->vec1_N,&array);
93: PetscMalloc(pcis->n*sizeof(PetscInt),&idx_I_local);
94: PetscMalloc(pcis->n*sizeof(PetscInt),&idx_B_local);
95: for (i=0, pcis->n_B=0, n_I=0; i<pcis->n; i++) {
96: if (array[i] == 1.0) { idx_I_local[n_I] = i; n_I++; }
97: else { idx_B_local[pcis->n_B] = i; pcis->n_B++; }
98: }
99: /* Getting the global numbering */
100: idx_B_global = idx_I_local + n_I; /* Just avoiding allocating extra memory, since we have vacant space */
101: idx_I_global = idx_B_local + pcis->n_B;
102: ISLocalToGlobalMappingApply(matis->mapping,pcis->n_B,idx_B_local,idx_B_global);
103: ISLocalToGlobalMappingApply(matis->mapping,n_I, idx_I_local,idx_I_global);
104: /* Creating the index sets. */
105: ISCreateGeneral(MPI_COMM_SELF,pcis->n_B,idx_B_local,PETSC_COPY_VALUES, &pcis->is_B_local);
106: ISCreateGeneral(MPI_COMM_SELF,pcis->n_B,idx_B_global,PETSC_COPY_VALUES,&pcis->is_B_global);
107: ISCreateGeneral(MPI_COMM_SELF,n_I ,idx_I_local,PETSC_COPY_VALUES, &pcis->is_I_local);
108: ISCreateGeneral(MPI_COMM_SELF,n_I ,idx_I_global,PETSC_COPY_VALUES,&pcis->is_I_global);
109: /* Freeing memory and restoring arrays */
110: PetscFree(idx_B_local);
111: PetscFree(idx_I_local);
112: VecRestoreArray(pcis->vec1_N,&array);
113: }
115: /*
116: Extracting the blocks A_II, A_BI, A_IB and A_BB from A. If the numbering
117: is such that interior nodes come first than the interface ones, we have
119: [ | ]
120: [ A_II | A_IB ]
121: A = [ | ]
122: [-----------+------]
123: [ A_BI | A_BB ]
124: */
126: MatGetSubMatrix(matis->A,pcis->is_I_local,pcis->is_I_local,MAT_INITIAL_MATRIX,&pcis->A_II);
127: MatGetSubMatrix(matis->A,pcis->is_I_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_IB);
128: MatGetSubMatrix(matis->A,pcis->is_B_local,pcis->is_I_local,MAT_INITIAL_MATRIX,&pcis->A_BI);
129: MatGetSubMatrix(matis->A,pcis->is_B_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_BB);
131: /*
132: Creating work vectors and arrays
133: */
134: /* pcis->vec1_N has already been created */
135: VecDuplicate(pcis->vec1_N,&pcis->vec2_N);
136: VecCreateSeq(PETSC_COMM_SELF,pcis->n-pcis->n_B,&pcis->vec1_D);
137: VecDuplicate(pcis->vec1_D,&pcis->vec2_D);
138: VecDuplicate(pcis->vec1_D,&pcis->vec3_D);
139: VecCreateSeq(PETSC_COMM_SELF,pcis->n_B,&pcis->vec1_B);
140: VecDuplicate(pcis->vec1_B,&pcis->vec2_B);
141: VecDuplicate(pcis->vec1_B,&pcis->vec3_B);
142: MatGetVecs(pc->pmat,&pcis->vec1_global,0);
143: PetscMalloc((pcis->n)*sizeof(PetscScalar),&pcis->work_N);
145: /* Creating the scatter contexts */
146: VecScatterCreate(pcis->vec1_global,pcis->is_I_global,pcis->vec1_D,(IS)0,&pcis->global_to_D);
147: VecScatterCreate(pcis->vec1_N,pcis->is_B_local,pcis->vec1_B,(IS)0,&pcis->N_to_B);
148: VecScatterCreate(pcis->vec1_global,pcis->is_B_global,pcis->vec1_B,(IS)0,&pcis->global_to_B);
150: /* Creating scaling "matrix" D */
151: if( !pcis->D ) {
152: VecSet(counter,0.0);
153: VecSet(pcis->vec1_N,pcis->scaling_factor);
154: VecScatterBegin(matis->ctx,pcis->vec1_N,counter,ADD_VALUES,SCATTER_REVERSE);
155: VecScatterEnd (matis->ctx,pcis->vec1_N,counter,ADD_VALUES,SCATTER_REVERSE);
156: VecScatterBegin(matis->ctx,counter,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);
157: VecScatterEnd (matis->ctx,counter,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);
158: VecDuplicate(pcis->vec1_B,&pcis->D);
159: VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);
160: VecScatterEnd (pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);
161: VecReciprocal(pcis->D);
162: VecScale(pcis->D,pcis->scaling_factor);
163: }
165: /* See historical note 01, at the bottom of this file. */
167: /*
168: Creating the KSP contexts for the local Dirichlet and Neumann problems.
169: */
170: {
171: PC pc_ctx;
172: /* Dirichlet */
173: KSPCreate(PETSC_COMM_SELF,&pcis->ksp_D);
174: PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_D,(PetscObject)pc,1);
175: KSPSetOperators(pcis->ksp_D,pcis->A_II,pcis->A_II,SAME_PRECONDITIONER);
176: KSPSetOptionsPrefix(pcis->ksp_D,"is_localD_");
177: KSPGetPC(pcis->ksp_D,&pc_ctx);
178: PCSetType(pc_ctx,PCLU);
179: KSPSetType(pcis->ksp_D,KSPPREONLY);
180: KSPSetFromOptions(pcis->ksp_D);
181: /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
182: KSPSetUp(pcis->ksp_D);
183: /* Neumann */
184: KSPCreate(PETSC_COMM_SELF,&pcis->ksp_N);
185: PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_N,(PetscObject)pc,1);
186: KSPSetOperators(pcis->ksp_N,matis->A,matis->A,SAME_PRECONDITIONER);
187: KSPSetOptionsPrefix(pcis->ksp_N,"is_localN_");
188: KSPGetPC(pcis->ksp_N,&pc_ctx);
189: PCSetType(pc_ctx,PCLU);
190: KSPSetType(pcis->ksp_N,KSPPREONLY);
191: KSPSetFromOptions(pcis->ksp_N);
192: {
193: PetscBool damp_fixed = PETSC_FALSE,
194: remove_nullspace_fixed = PETSC_FALSE,
195: set_damping_factor_floating = PETSC_FALSE,
196: not_damp_floating = PETSC_FALSE,
197: not_remove_nullspace_floating = PETSC_FALSE;
198: PetscReal fixed_factor,
199: floating_factor;
201: PetscOptionsGetReal(((PetscObject)pc_ctx)->prefix,"-pc_is_damp_fixed",&fixed_factor,&damp_fixed);
202: if (!damp_fixed) { fixed_factor = 0.0; }
203: PetscOptionsGetBool(((PetscObject)pc_ctx)->prefix,"-pc_is_damp_fixed",&damp_fixed,PETSC_NULL);
205: PetscOptionsGetBool(((PetscObject)pc_ctx)->prefix,"-pc_is_remove_nullspace_fixed",&remove_nullspace_fixed,PETSC_NULL);
207: PetscOptionsGetReal(((PetscObject)pc_ctx)->prefix,"-pc_is_set_damping_factor_floating",
208: &floating_factor,&set_damping_factor_floating);
209: if (!set_damping_factor_floating) { floating_factor = 0.0; }
210: PetscOptionsGetBool(((PetscObject)pc_ctx)->prefix,"-pc_is_set_damping_factor_floating",&set_damping_factor_floating,PETSC_NULL);
211: if (!set_damping_factor_floating) { floating_factor = 1.e-12; }
213: PetscOptionsGetBool(((PetscObject)pc_ctx)->prefix,"-pc_is_not_damp_floating",¬_damp_floating,PETSC_NULL);
215: PetscOptionsGetBool(((PetscObject)pc_ctx)->prefix,"-pc_is_not_remove_nullspace_floating",¬_remove_nullspace_floating,PETSC_NULL);
217: if (pcis->pure_neumann) { /* floating subdomain */
218: if (!(not_damp_floating)) {
219: PCFactorSetShiftType(pc_ctx,MAT_SHIFT_NONZERO);
220: PCFactorSetShiftAmount(pc_ctx,floating_factor);
221: }
222: if (!(not_remove_nullspace_floating)){
223: MatNullSpace nullsp;
224: MatNullSpaceCreate(PETSC_COMM_SELF,PETSC_TRUE,0,PETSC_NULL,&nullsp);
225: KSPSetNullSpace(pcis->ksp_N,nullsp);
226: MatNullSpaceDestroy(&nullsp);
227: }
228: } else { /* fixed subdomain */
229: if (damp_fixed) {
230: PCFactorSetShiftType(pc_ctx,MAT_SHIFT_NONZERO);
231: PCFactorSetShiftAmount(pc_ctx,floating_factor);
232: }
233: if (remove_nullspace_fixed) {
234: MatNullSpace nullsp;
235: MatNullSpaceCreate(PETSC_COMM_SELF,PETSC_TRUE,0,PETSC_NULL,&nullsp);
236: KSPSetNullSpace(pcis->ksp_N,nullsp);
237: MatNullSpaceDestroy(&nullsp);
238: }
239: }
240: }
241: /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
242: KSPSetUp(pcis->ksp_N);
243: }
245: ISLocalToGlobalMappingGetInfo(((Mat_IS*)(pc->mat->data))->mapping,&(pcis->n_neigh),&(pcis->neigh),&(pcis->n_shared),&(pcis->shared));
246: pcis->ISLocalToGlobalMappingGetInfoWasCalled = PETSC_TRUE;
247: VecDestroy(&counter);
248: return(0);
249: }
251: /* -------------------------------------------------------------------------- */
252: /*
253: PCISDestroy -
254: */
257: PetscErrorCode PCISDestroy(PC pc)
258: {
259: PC_IS *pcis = (PC_IS*)(pc->data);
263: ISDestroy(&pcis->is_B_local);
264: ISDestroy(&pcis->is_I_local);
265: ISDestroy(&pcis->is_B_global);
266: ISDestroy(&pcis->is_I_global);
267: MatDestroy(&pcis->A_II);
268: MatDestroy(&pcis->A_IB);
269: MatDestroy(&pcis->A_BI);
270: MatDestroy(&pcis->A_BB);
271: VecDestroy(&pcis->D);
272: KSPDestroy(&pcis->ksp_N);
273: KSPDestroy(&pcis->ksp_D);
274: VecDestroy(&pcis->vec1_N);
275: VecDestroy(&pcis->vec2_N);
276: VecDestroy(&pcis->vec1_D);
277: VecDestroy(&pcis->vec2_D);
278: VecDestroy(&pcis->vec3_D);
279: VecDestroy(&pcis->vec1_B);
280: VecDestroy(&pcis->vec2_B);
281: VecDestroy(&pcis->vec3_B);
282: VecDestroy(&pcis->vec1_global);
283: VecScatterDestroy(&pcis->global_to_D);
284: VecScatterDestroy(&pcis->N_to_B);
285: VecScatterDestroy(&pcis->global_to_B);
286: PetscFree(pcis->work_N);
287: if (pcis->ISLocalToGlobalMappingGetInfoWasCalled) {
288: ISLocalToGlobalMappingRestoreInfo((ISLocalToGlobalMapping)0,&(pcis->n_neigh),&(pcis->neigh),&(pcis->n_shared),&(pcis->shared));
289: }
290: return(0);
291: }
293: /* -------------------------------------------------------------------------- */
294: /*
295: PCISCreate -
296: */
299: PetscErrorCode PCISCreate(PC pc)
300: {
301: PC_IS *pcis = (PC_IS*)(pc->data);
305: pcis->is_B_local = 0;
306: pcis->is_I_local = 0;
307: pcis->is_B_global = 0;
308: pcis->is_I_global = 0;
309: pcis->A_II = 0;
310: pcis->A_IB = 0;
311: pcis->A_BI = 0;
312: pcis->A_BB = 0;
313: pcis->D = 0;
314: pcis->ksp_N = 0;
315: pcis->ksp_D = 0;
316: pcis->vec1_N = 0;
317: pcis->vec2_N = 0;
318: pcis->vec1_D = 0;
319: pcis->vec2_D = 0;
320: pcis->vec3_D = 0;
321: pcis->vec1_B = 0;
322: pcis->vec2_B = 0;
323: pcis->vec3_B = 0;
324: pcis->vec1_global = 0;
325: pcis->work_N = 0;
326: pcis->global_to_D = 0;
327: pcis->N_to_B = 0;
328: pcis->global_to_B = 0;
329: pcis->ISLocalToGlobalMappingGetInfoWasCalled = PETSC_FALSE;
330: pcis->scaling_factor = 1.0;
331: /* composing functions */
332: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCISSetSubdomainScalingFactor_C","PCISSetSubdomainScalingFactor_IS",
333: PCISSetSubdomainScalingFactor_IS);
334: return(0);
335: }
337: /* -------------------------------------------------------------------------- */
338: /*
339: PCISApplySchur -
341: Input parameters:
342: . pc - preconditioner context
343: . v - vector to which the Schur complement is to be applied (it is NOT modified inside this function, UNLESS vec2_B is null)
345: Output parameters:
346: . vec1_B - result of Schur complement applied to chunk
347: . vec2_B - garbage (used as work space), or null (and v is used as workspace)
348: . vec1_D - garbage (used as work space)
349: . vec2_D - garbage (used as work space)
351: */
354: PetscErrorCode PCISApplySchur(PC pc, Vec v, Vec vec1_B, Vec vec2_B, Vec vec1_D, Vec vec2_D)
355: {
357: PC_IS *pcis = (PC_IS*)(pc->data);
360: if (!vec2_B) { vec2_B = v; }
362: MatMult(pcis->A_BB,v,vec1_B);
363: MatMult(pcis->A_IB,v,vec1_D);
364: KSPSolve(pcis->ksp_D,vec1_D,vec2_D);
365: MatMult(pcis->A_BI,vec2_D,vec2_B);
366: VecAXPY(vec1_B,-1.0,vec2_B);
367: return(0);
368: }
370: /* -------------------------------------------------------------------------- */
371: /*
372: PCISScatterArrayNToVecB - Scatters interface node values from a big array (of all local nodes, interior or interface,
373: including ghosts) into an interface vector, when in SCATTER_FORWARD mode, or vice-versa, when in SCATTER_REVERSE
374: mode.
376: Input parameters:
377: . pc - preconditioner context
378: . array_N - [when in SCATTER_FORWARD mode] Array to be scattered into the vector
379: . v_B - [when in SCATTER_REVERSE mode] Vector to be scattered into the array
381: Output parameter:
382: . array_N - [when in SCATTER_REVERSE mode] Array to receive the scattered vector
383: . v_B - [when in SCATTER_FORWARD mode] Vector to receive the scattered array
385: Notes:
386: The entries in the array that do not correspond to interface nodes remain unaltered.
387: */
390: PetscErrorCode PCISScatterArrayNToVecB (PetscScalar *array_N, Vec v_B, InsertMode imode, ScatterMode smode, PC pc)
391: {
392: PetscInt i;
393: const PetscInt *idex;
395: PetscScalar *array_B;
396: PC_IS *pcis = (PC_IS*)(pc->data);
399: VecGetArray(v_B,&array_B);
400: ISGetIndices(pcis->is_B_local,&idex);
402: if (smode == SCATTER_FORWARD) {
403: if (imode == INSERT_VALUES) {
404: for (i=0; i<pcis->n_B; i++) { array_B[i] = array_N[idex[i]]; }
405: } else { /* ADD_VALUES */
406: for (i=0; i<pcis->n_B; i++) { array_B[i] += array_N[idex[i]]; }
407: }
408: } else { /* SCATTER_REVERSE */
409: if (imode == INSERT_VALUES) {
410: for (i=0; i<pcis->n_B; i++) { array_N[idex[i]] = array_B[i]; }
411: } else { /* ADD_VALUES */
412: for (i=0; i<pcis->n_B; i++) { array_N[idex[i]] += array_B[i]; }
413: }
414: }
415: ISRestoreIndices(pcis->is_B_local,&idex);
416: VecRestoreArray(v_B,&array_B);
417: return(0);
418: }
420: /* -------------------------------------------------------------------------- */
421: /*
422: PCISApplyInvSchur - Solves the Neumann problem related to applying the inverse of the Schur complement.
423: More precisely, solves the problem:
424: [ A_II A_IB ] [ . ] [ 0 ]
425: [ ] [ ] = [ ]
426: [ A_BI A_BB ] [ x ] [ b ]
428: Input parameters:
429: . pc - preconditioner context
430: . b - vector of local interface nodes (including ghosts)
432: Output parameters:
433: . x - vector of local interface nodes (including ghosts); returns the application of the inverse of the Schur
434: complement to b
435: . vec1_N - vector of local nodes (interior and interface, including ghosts); returns garbage (used as work space)
436: . vec2_N - vector of local nodes (interior and interface, including ghosts); returns garbage (used as work space)
438: */
441: PetscErrorCode PCISApplyInvSchur (PC pc, Vec b, Vec x, Vec vec1_N, Vec vec2_N)
442: {
444: PC_IS *pcis = (PC_IS*)(pc->data);
447: /*
448: Neumann solvers.
449: Applying the inverse of the local Schur complement, i.e, solving a Neumann
450: Problem with zero at the interior nodes of the RHS and extracting the interface
451: part of the solution. inverse Schur complement is applied to b and the result
452: is stored in x.
453: */
454: /* Setting the RHS vec1_N */
455: VecSet(vec1_N,0.0);
456: VecScatterBegin(pcis->N_to_B,b,vec1_N,INSERT_VALUES,SCATTER_REVERSE);
457: VecScatterEnd (pcis->N_to_B,b,vec1_N,INSERT_VALUES,SCATTER_REVERSE);
458: /* Checking for consistency of the RHS */
459: {
460: PetscBool flg = PETSC_FALSE;
461: PetscOptionsGetBool(PETSC_NULL,"-pc_is_check_consistency",&flg,PETSC_NULL);
462: if (flg) {
463: PetscScalar average;
464: PetscViewer viewer;
465: PetscViewerASCIIGetStdout(((PetscObject)pc)->comm,&viewer);
467: VecSum(vec1_N,&average);
468: average = average / ((PetscReal)pcis->n);
469: PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
470: if (pcis->pure_neumann) {
471: PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d is floating. Average = % 1.14e\n",PetscGlobalRank,PetscAbsScalar(average));
472: } else {
473: PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d is fixed. Average = % 1.14e\n",PetscGlobalRank,PetscAbsScalar(average));
474: }
475: PetscViewerFlush(viewer);
476: PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
477: }
478: }
479: /* Solving the system for vec2_N */
480: KSPSolve(pcis->ksp_N,vec1_N,vec2_N);
481: /* Extracting the local interface vector out of the solution */
482: VecScatterBegin(pcis->N_to_B,vec2_N,x,INSERT_VALUES,SCATTER_FORWARD);
483: VecScatterEnd (pcis->N_to_B,vec2_N,x,INSERT_VALUES,SCATTER_FORWARD);
484: return(0);
485: }