xref: /petsc/src/ksp/pc/impls/fieldsplit/fieldsplit.c (revision 0971522ceecba37df062696def42d3ac8bbed20b)
1 /*
2 
3 */
4 #include "src/ksp/pc/pcimpl.h"     /*I "petscpc.h" I*/
5 
6 typedef struct _PC_FieldSplitLink *PC_FieldSplitLink;
7 struct _PC_FieldSplitLink {
8   PC                pc;
9   Vec               x,y;
10   PetscInt          nfields;
11   PetscInt          *fields;
12   PC_FieldSplitLink next;
13 };
14 
15 typedef struct {
16   PetscInt          bs;
17   PetscInt          nsplits;
18   PC_FieldSplitLink head;
19   Vec               *x,*y;
20 } PC_FieldSplit;
21 
22 #undef __FUNCT__
23 #define __FUNCT__ "PCView_FieldSplit"
24 static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer)
25 {
26   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
27   PetscErrorCode    ierr;
28   PetscTruth        iascii;
29   PetscInt          i,j;
30   PC_FieldSplitLink link = jac->head;
31 
32   PetscFunctionBegin;
33   ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr);
34   if (iascii) {
35     ierr = PetscViewerASCIIPrintf(viewer,"  FieldSplit: total splits = %D",jac->nsplits);CHKERRQ(ierr);
36     ierr = PetscViewerASCIIPrintf(viewer,"  Solver info for each split is in the following PC objects:\n");CHKERRQ(ierr);
37     ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
38     for (i=0; i<jac->nsplits; i++) {
39       ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);CHKERRQ(ierr);
40       for (j=0; j<link->nfields; j++) {
41         ierr = PetscViewerASCIIPrintf(viewer,"%D \n",link->fields[j]);CHKERRQ(ierr);
42       }
43       ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr);
44       ierr = PCView(link->pc,viewer);CHKERRQ(ierr);
45       link = link->next;
46     }
47     ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
48   } else {
49     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCFieldSplit",((PetscObject)viewer)->type_name);
50   }
51   PetscFunctionReturn(0);
52 }
53 
54 #undef __FUNCT__
55 #define __FUNCT__ "PCSetUp_FieldSplit"
56 static PetscErrorCode PCSetUp_FieldSplit(PC pc)
57 {
58   PC_FieldSplit     *jac  = (PC_FieldSplit*)pc->data;
59   PetscErrorCode    ierr;
60   PC_FieldSplitLink link = jac->head;
61 
62   PetscFunctionBegin;
63   while (link) {
64     ierr = PCSetUp(link->pc);CHKERRQ(ierr);
65     link = link->next;
66   }
67   PetscFunctionReturn(0);
68 }
69 
70 #undef __FUNCT__
71 #define __FUNCT__ "PCApply_FieldSplit"
72 static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y)
73 {
74   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
75   PetscErrorCode    ierr;
76   PC_FieldSplitLink link = jac->head;
77 
78   PetscFunctionBegin;
79   ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr);
80   while (link) {
81     ierr = PCApply(link->pc,link->x,link->y);CHKERRQ(ierr);
82     link = link->next;
83   }
84   ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr);
85   PetscFunctionReturn(0);
86 }
87 
88 #undef __FUNCT__
89 #define __FUNCT__ "PCDestroy_FieldSplit"
90 static PetscErrorCode PCDestroy_FieldSplit(PC pc)
91 {
92   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
93   PetscErrorCode    ierr;
94   PC_FieldSplitLink link = jac->head;
95 
96   PetscFunctionBegin;
97   while (link) {
98     ierr = PCDestroy(link->pc);CHKERRQ(ierr);
99     ierr = PetscFree2(link,link->fields);CHKERRQ(ierr);
100   }
101   if (jac->x) {ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr);}
102   ierr = PetscFree(jac);CHKERRQ(ierr);
103   PetscFunctionReturn(0);
104 }
105 
106 #undef __FUNCT__
107 #define __FUNCT__ "PCSetFromOptions_FieldSplit"
108 static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc)
109 {
110   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
111   PetscErrorCode    ierr;
112   PC_FieldSplitLink link = jac->head;
113 
114   PetscFunctionBegin;
115   if (!link) { /* user never set fields, so set them ourselves */
116     link = jac->head;
117   }
118   while (link) {
119     ierr = PCApply(link->pc,link->x,link->y);CHKERRQ(ierr);
120     link = link->next;
121   }
122   PetscFunctionReturn(0);
123 }
124 
125 /*------------------------------------------------------------------------------------*/
126 
127 EXTERN_C_BEGIN
128 #undef __FUNCT__
129 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit"
130 PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,PetscInt n,PetscInt *fields)
131 {
132   PC_FieldSplit     *jac;
133   PetscErrorCode    ierr;
134   PetscInt          *myfields;
135   PC_FieldSplitLink link,next = jac->head;
136 
137   PetscFunctionBegin;
138   if (n <= 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Negative number of fields requested");
139   ierr = PetscMalloc2(1,struct _PC_FieldSplitLink,&link,n,PetscInt,&myfields);CHKERRQ(ierr);
140   ierr = PetscMemcpy(myfields,fields,n*sizeof(PetscInt));CHKERRQ(ierr);
141   link->fields  = myfields;
142   link->nfields = n;
143 
144   if (!next) {
145     jac->head = link;
146   } else {
147     while (next->next) {
148       next = next->next;
149     }
150     next->next = link;
151   }
152   jac->nsplits++;
153   PetscFunctionReturn(0);
154 }
155 EXTERN_C_END
156 
157 
158 EXTERN_C_BEGIN
159 #undef __FUNCT__
160 #define __FUNCT__ "PCFieldSplitGetSubPC_FieldSplit"
161 PetscErrorCode PCFieldSplitGetSubPC_FieldSplit(PC pc,PetscInt *n,PC **subpc)
162 {
163   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
164   PetscErrorCode    ierr;
165   PetscInt          cnt = 0;
166   PC_FieldSplitLink link;
167 
168   PetscFunctionBegin;
169   ierr = PetscMalloc(jac->nsplits*sizeof(PC*),subpc);CHKERRQ(ierr);
170   while (link) {
171     (*subpc)[cnt++] = link->pc;
172     link = link->next;
173   }
174   if (cnt != jac->nsplits) SETERRQ2(PETSC_ERR_PLIB,"Corrupt PCFIELDSPLIT object: number splits in linked list %D in object %D",cnt,jac->nsplits);
175   *n = jac->nsplits;
176   PetscFunctionReturn(0);
177 }
178 EXTERN_C_END
179 
180 #undef __FUNCT__
181 #define __FUNCT__ "PCFieldSplitSetFields"
182 /*@
183     PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner
184 
185     Collective on PC
186 
187     Input Parameters:
188 +   pc  - the preconditioner context
189 .   n - the number of fields in this split
190 .   fields - the fields in this split
191 
192     Level: intermediate
193 
194 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT
195 
196 @*/
197 PetscErrorCode PCFieldSplitSetFields(PC pc,PetscInt n, PetscInt *fields)
198 {
199   PetscErrorCode ierr,(*f)(PC,PetscInt,PetscInt *);
200 
201   PetscFunctionBegin;
202   PetscValidHeaderSpecific(pc,PC_COOKIE,1);
203   ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitSetFields_C",(void (**)(void))&f);CHKERRQ(ierr);
204   if (f) {
205     ierr = (*f)(pc,n,fields);CHKERRQ(ierr);
206   }
207   PetscFunctionReturn(0);
208 }
209 
210 #undef __FUNCT__
211 #define __FUNCT__ "PCFieldSplitGetSubPC"
212 /*@C
213    PCFieldSplitGetSubPC - Gets the PC contexts for all splits
214 
215    Collective on PC
216 
217    Input Parameter:
218 .  pc - the preconditioner context
219 
220    Output Parameters:
221 +  n - the number of split
222 -  pc - the array of PC contexts
223 
224    Note:
225    After PCFieldSplitGetSubPC() the array of PCs IS to be freed
226 
227    You must call KSPSetUp() before calling PCFieldSplitGetSubPC().
228 
229    Level: advanced
230 
231 .keywords: PC,
232 
233 .seealso: PCFIELDSPLIT
234 @*/
235 PetscErrorCode PCFieldSplitGetSubPC(PC pc,PetscInt *n,PC *subpc[])
236 {
237   PetscErrorCode ierr,(*f)(PC,PetscInt*,PC **);
238 
239   PetscFunctionBegin;
240   PetscValidHeaderSpecific(pc,PC_COOKIE,1);
241   PetscValidIntPointer(n,2);
242   ierr = PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitGetSubPC_C",(void (**)(void))&f);CHKERRQ(ierr);
243   if (f) {
244     ierr = (*f)(pc,n,subpc);CHKERRQ(ierr);
245   } else {
246     SETERRQ(PETSC_ERR_ARG_WRONG,"Cannot get subpc for this type of PC");
247   }
248 
249  PetscFunctionReturn(0);
250 }
251 
252 /* -------------------------------------------------------------------------------------*/
253 /*MC
254    PCFieldSplit - Preconditioner created by combining seperate preconditioners for individual
255                   fields or groups of fields
256 
257 
258      To set options on the solvers for each block append -sub_ to all the PC
259         options database keys. For example, -sub_pc_type ilu -sub_pc_ilu_levels 1
260 
261      To set the options on the solvers seperate for each block call PCFieldSplitGetSubPC()
262          and set the options directly on the resulting PC object
263 
264    Level: intermediate
265 
266    Concepts: physics based preconditioners
267 
268 .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
269            PCFieldSplitGetSubPC(), PCFieldSplitSetFields()
270 M*/
271 
272 EXTERN_C_BEGIN
273 #undef __FUNCT__
274 #define __FUNCT__ "PCCreate_FieldSplit"
275 PetscErrorCode PCCreate_FieldSplit(PC pc)
276 {
277   PetscErrorCode ierr;
278   PC_FieldSplit  *jac;
279 
280   PetscFunctionBegin;
281   ierr = PetscNew(PC_FieldSplit,&jac);CHKERRQ(ierr);
282   PetscLogObjectMemory(pc,sizeof(PC_FieldSplit));
283   ierr = PetscMemzero(jac,sizeof(PC_FieldSplit));CHKERRQ(ierr);
284   jac->bs      = -1;
285   jac->nsplits = 0;
286   pc->data     = (void*)jac;
287 
288   pc->ops->apply             = PCApply_FieldSplit;
289   pc->ops->setup             = PCSetUp_FieldSplit;
290   pc->ops->destroy           = PCDestroy_FieldSplit;
291   pc->ops->setfromoptions    = PCSetFromOptions_FieldSplit;
292   pc->ops->view              = PCView_FieldSplit;
293   pc->ops->applyrichardson   = 0;
294 
295   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubPC_C","PCFieldSplitGetSubPC_FieldSplit",
296                     PCFieldSplitGetSubPC_FieldSplit);CHKERRQ(ierr);
297   ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","PCFieldSplitSetFields_FieldSplit",
298                     PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr);
299   PetscFunctionReturn(0);
300 }
301 EXTERN_C_END
302 
303 
304