Actual source code: damgsnesad.c
1: #define PETSCSNES_DLL
2:
3: #include petscda.h
4: #include petscmg.h
5: #include petscdmmg.h
6: #include src/inline/ilu.h
7: #include src/snes/impls/ls/ls.h
10: EXTERN PetscErrorCode NLFRelax_DAAD(NLF,MatSORType,PetscInt,Vec);
11: EXTERN PetscErrorCode NLFRelax_DAAD4(NLF,MatSORType,PetscInt,Vec);
12: EXTERN PetscErrorCode NLFRelax_DAAD9(NLF,MatSORType,PetscInt,Vec);
13: EXTERN PetscErrorCode NLFRelax_DAADb(NLF,MatSORType,PetscInt,Vec);
15: EXTERN PetscErrorCode DMMGFormFunction(SNES,Vec,Vec,void *);
16: EXTERN PetscErrorCode SNESLSCheckLocalMin_Private(Mat,Vec,Vec,PetscReal,PetscTruth*);
20: /*
21: DMMGComputeJacobianWithAdic - Evaluates the Jacobian via Adic when the user has provided
22: a local function evaluation routine.
23: */
24: PetscErrorCode DMMGComputeJacobianWithAdic(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
25: {
26: DMMG dmmg = (DMMG) ptr;
28: Vec localX;
29: DA da = (DA) dmmg->dm;
32: DAGetLocalVector(da,&localX);
33: DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
34: DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
35: DAComputeJacobian1WithAdic(da,localX,*B,dmmg->user);
36: DARestoreLocalVector(da,&localX);
37: /* Assemble true Jacobian; if it is different */
38: if (*J != *B) {
39: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
40: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
41: }
42: MatSetOption(*B,MAT_NEW_NONZERO_LOCATION_ERR);
43: *flag = SAME_NONZERO_PATTERN;
44: return(0);
45: }
49: /*@
50: SNESDAComputeJacobianWithAdic - This is a universal Jacobian evaluation routine
51: that may be used with SNESSetJacobian() as long as the user context has a DA as
52: its first record and DASetLocalAdicFunction() has been called.
54: Collective on SNES
56: Input Parameters:
57: + snes - the SNES context
58: . X - input vector
59: . J - Jacobian
60: . B - Jacobian used in preconditioner (usally same as J)
61: . flag - indicates if the matrix changed its structure
62: - ptr - optional user-defined context, as set by SNESSetFunction()
64: Level: intermediate
66: .seealso: DASetLocalFunction(), DASetLocalAdicFunction(), SNESSetFunction(), SNESSetJacobian()
68: @*/
69: PetscErrorCode SNESDAComputeJacobianWithAdic(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
70: {
71: DA da = *(DA*) ptr;
73: Vec localX;
76: DAGetLocalVector(da,&localX);
77: DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
78: DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
79: DAComputeJacobian1WithAdic(da,localX,*B,ptr);
80: DARestoreLocalVector(da,&localX);
81: /* Assemble true Jacobian; if it is different */
82: if (*J != *B) {
83: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
84: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
85: }
86: MatSetOption(*B,MAT_NEW_NONZERO_LOCATION_ERR);
87: *flag = SAME_NONZERO_PATTERN;
88: return(0);
89: }
91: #include src/ksp/pc/impls/mg/mgimpl.h
92: /*
93: This is pre-beta FAS code. It's design should not be taken seriously!
95: R is the usual multigrid restriction (e.g. the tranpose of peicewise linear interpolation)
96: Q is either a scaled injection or the usual R
97: */
100: PetscErrorCode DMMGSolveFAS(DMMG *dmmg,PetscInt level)
101: {
103: PetscInt i,j,k;
104: PetscReal norm;
105: PC_MG **mg;
106: PC pc;
109: VecSet(dmmg[level]->r,0.0);
110: for (j=1; j<=level; j++) {
111: if (!dmmg[j]->inject) {
112: DMGetInjection(dmmg[j-1]->dm,dmmg[j]->dm,&dmmg[j]->inject);
113: }
114: }
116: KSPGetPC(dmmg[level]->ksp,&pc);
117: mg = ((PC_MG**)pc->data);
119: for (i=0; i<100; i++) {
121: for (j=level; j>0; j--) {
123: /* Relax residual_fine - F(x_fine) = 0 */
124: for (k=0; k<dmmg[j]->presmooth; k++) {
125: NLFRelax_DAAD(dmmg[j]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[j]->x);
126: }
128: /* R*(residual_fine - F(x_fine)) */
129: DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
130: VecAYPX(dmmg[j]->w,-1.0,dmmg[j]->r);
132: if (j == level || dmmg[j]->monitorall) {
133: /* norm( residual_fine - f(x_fine) ) */
134: VecNorm(dmmg[j]->w,NORM_2,&norm);
135: if (j == level) {
136: if (norm < dmmg[level]->abstol) goto theend;
137: if (i == 0) {
138: dmmg[level]->rrtol = norm*dmmg[level]->rtol;
139: } else {
140: if (norm < dmmg[level]->rrtol) goto theend;
141: }
142: }
143: }
145: if (dmmg[j]->monitorall) {
146: for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm," ");}
147: PetscPrintf(dmmg[j]->comm,"FAS function norm %G\n",norm);
148: }
149: MatRestrict(mg[j]->restrct,dmmg[j]->w,dmmg[j-1]->r);
150:
151: /* F(Q*x_fine) */
152: VecScatterBegin(dmmg[j]->inject,dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD);
153: VecScatterEnd(dmmg[j]->inject,dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD);
154: DMMGFormFunction(0,dmmg[j-1]->x,dmmg[j-1]->w,dmmg[j-1]);
156: /* residual_coarse = F(Q*x_fine) + R*(residual_fine - F(x_fine)) */
157: VecAYPX(dmmg[j-1]->r,1.0,dmmg[j-1]->w);
159: /* save Q*x_fine into b (needed when interpolating compute x back up */
160: VecCopy(dmmg[j-1]->x,dmmg[j-1]->b);
161: }
163: for (j=0; j<dmmg[0]->coarsesmooth; j++) {
164: NLFRelax_DAAD(dmmg[0]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[0]->x);
165: }
166: if (dmmg[0]->monitorall){
167: DMMGFormFunction(0,dmmg[0]->x,dmmg[0]->w,dmmg[0]);
168: VecAXPY(dmmg[0]->w,-1.0,dmmg[0]->r);
169: VecNorm(dmmg[0]->w,NORM_2,&norm);
170: for (k=0; k<level+1; k++) {PetscPrintf(dmmg[0]->comm," ");}
171: PetscPrintf(dmmg[0]->comm,"FAS coarse grid function norm %G\n",norm);
172: }
174: for (j=1; j<=level; j++) {
175: /* x_fine = x_fine + R'*(x_coarse - Q*x_fine) */
176: VecAXPY(dmmg[j-1]->x,-1.0,dmmg[j-1]->b);
177: MatInterpolateAdd(mg[j]->interpolate,dmmg[j-1]->x,dmmg[j]->x,dmmg[j]->x);
179: if (dmmg[j]->monitorall) {
180: /* norm( F(x_fine) - residual_fine ) */
181: DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
182: VecAXPY(dmmg[j]->w,-1.0,dmmg[j]->r);
183: VecNorm(dmmg[j]->w,NORM_2,&norm);
184: for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm," ");}
185: PetscPrintf(dmmg[j]->comm,"FAS function norm %G\n",norm);
186: }
188: /* Relax residual_fine - F(x_fine) = 0 */
189: for (k=0; k<dmmg[j]->postsmooth; k++) {
190: NLFRelax_DAAD(dmmg[j]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[j]->x);
191: }
193: if (dmmg[j]->monitorall) {
194: /* norm( F(x_fine) - residual_fine ) */
195: DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
196: VecAXPY(dmmg[j]->w,-1.0,dmmg[j]->r);
197: VecNorm(dmmg[j]->w,NORM_2,&norm);
198: for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm," ");}
199: PetscPrintf(dmmg[j]->comm,"FAS function norm %G\n",norm);
200: }
201: }
203: if (dmmg[level]->monitor){
204: DMMGFormFunction(0,dmmg[level]->x,dmmg[level]->w,dmmg[level]);
205: VecNorm(dmmg[level]->w,NORM_2,&norm);
206: PetscPrintf(dmmg[level]->comm,"%D FAS function norm %G\n",i+1,norm);
207: }
208: }
209: theend:
210: return(0);
211: }
213: /*
214: This is the point-block version of FAS
215: */
218: PetscErrorCode DMMGSolveFASb(DMMG *dmmg,PetscInt level)
219: {
221: PetscInt i,j,k;
222: PetscReal norm;
223: PC_MG **mg;
224: PC pc;
227: VecSet(dmmg[level]->r,0.0);
228: for (j=1; j<=level; j++) {
229: if (!dmmg[j]->inject) {
230: DMGetInjection(dmmg[j-1]->dm,dmmg[j]->dm,&dmmg[j]->inject);
231: }
232: }
234: KSPGetPC(dmmg[level]->ksp,&pc);
235: mg = ((PC_MG**)pc->data);
237: for (i=0; i<100; i++) {
239: for (j=level; j>0; j--) {
241: /* Relax residual_fine - F(x_fine) = 0 */
242: for (k=0; k<dmmg[j]->presmooth; k++) {
243: NLFRelax_DAADb(dmmg[j]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[j]->x);
244: }
246: /* R*(residual_fine - F(x_fine)) */
247: DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
248: VecAYPX(dmmg[j]->w,-1.0,dmmg[j]->r);
250: if (j == level || dmmg[j]->monitorall) {
251: /* norm( residual_fine - f(x_fine) ) */
252: VecNorm(dmmg[j]->w,NORM_2,&norm);
253: if (j == level) {
254: if (norm < dmmg[level]->abstol) goto theend;
255: if (i == 0) {
256: dmmg[level]->rrtol = norm*dmmg[level]->rtol;
257: } else {
258: if (norm < dmmg[level]->rrtol) goto theend;
259: }
260: }
261: }
263: if (dmmg[j]->monitorall) {
264: for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm," ");}
265: PetscPrintf(dmmg[j]->comm,"FAS function norm %G\n",norm);
266: }
267: MatRestrict(mg[j]->restrct,dmmg[j]->w,dmmg[j-1]->r);
268:
269: /* F(Q*x_fine) */
270: VecScatterBegin(dmmg[j]->inject,dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD);
271: VecScatterEnd(dmmg[j]->inject,dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD);
272: DMMGFormFunction(0,dmmg[j-1]->x,dmmg[j-1]->w,dmmg[j-1]);
274: /* residual_coarse = F(Q*x_fine) + R*(residual_fine - F(x_fine)) */
275: VecAYPX(dmmg[j-1]->r,1.0,dmmg[j-1]->w);
277: /* save Q*x_fine into b (needed when interpolating compute x back up */
278: VecCopy(dmmg[j-1]->x,dmmg[j-1]->b);
279: }
281: for (j=0; j<dmmg[0]->coarsesmooth; j++) {
282: NLFRelax_DAADb(dmmg[0]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[0]->x);
283: }
284: if (dmmg[0]->monitorall){
285: DMMGFormFunction(0,dmmg[0]->x,dmmg[0]->w,dmmg[0]);
286: VecAXPY(dmmg[0]->w,-1.0,dmmg[0]->r);
287: VecNorm(dmmg[0]->w,NORM_2,&norm);
288: for (k=0; k<level+1; k++) {PetscPrintf(dmmg[0]->comm," ");}
289: PetscPrintf(dmmg[0]->comm,"FAS coarse grid function norm %G\n",norm);
290: }
292: for (j=1; j<=level; j++) {
293: /* x_fine = x_fine + R'*(x_coarse - Q*x_fine) */
294: VecAXPY(dmmg[j-1]->x,-1.0,dmmg[j-1]->b);
295: MatInterpolateAdd(mg[j]->interpolate,dmmg[j-1]->x,dmmg[j]->x,dmmg[j]->x);
297: if (dmmg[j]->monitorall) {
298: /* norm( F(x_fine) - residual_fine ) */
299: DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
300: VecAXPY(dmmg[j]->w,-1.0,dmmg[j]->r);
301: VecNorm(dmmg[j]->w,NORM_2,&norm);
302: for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm," ");}
303: PetscPrintf(dmmg[j]->comm,"FAS function norm %G\n",norm);
304: }
306: /* Relax residual_fine - F(x_fine) = 0 */
307: for (k=0; k<dmmg[j]->postsmooth; k++) {
308: NLFRelax_DAADb(dmmg[j]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[j]->x);
309: }
311: if (dmmg[j]->monitorall) {
312: /* norm( F(x_fine) - residual_fine ) */
313: DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
314: VecAXPY(dmmg[j]->w,-1.0,dmmg[j]->r);
315: VecNorm(dmmg[j]->w,NORM_2,&norm);
316: for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm," ");}
317: PetscPrintf(dmmg[j]->comm,"FAS function norm %G\n",norm);
318: }
319: }
321: if (dmmg[level]->monitor){
322: DMMGFormFunction(0,dmmg[level]->x,dmmg[level]->w,dmmg[level]);
323: VecNorm(dmmg[level]->w,NORM_2,&norm);
324: PetscPrintf(dmmg[level]->comm,"%D FAS function norm %G\n",i+1,norm);
325: }
326: }
327: theend:
328: return(0);
329: }
332: #include "adic/ad_utils.h"
337: PetscErrorCode PetscADView(PetscInt N,PetscInt nc,double *ptr,PetscViewer viewer)
338: {
339: PetscInt i,j,nlen = PetscADGetDerivTypeSize();
340: char *cptr = (char*)ptr;
341: double *values;
345: for (i=0; i<N; i++) {
346: PetscPrintf(PETSC_COMM_SELF,"Element %D value %G derivatives: ",i,*(double*)cptr);
347: values = PetscADGetGradArray(cptr);
348: for (j=0; j<nc; j++) {
349: PetscPrintf(PETSC_COMM_SELF,"%G ",*values++);
350: }
351: PetscPrintf(PETSC_COMM_SELF,"\n");
352: cptr += nlen;
353: }
355: return(0);
356: }
360: PetscErrorCode DMMGSolveFAS4(DMMG *dmmg,PetscInt level)
361: {
363: PetscInt i,j,k;
364: PetscReal norm;
365: PetscScalar zero = 0.0,mone = -1.0,one = 1.0;
366: PC_MG **mg;
367: PC pc;
370: VecSet(dmmg[level]->r,zero);
371: for (j=1; j<=level; j++) {
372: if (!dmmg[j]->inject) {
373: DMGetInjection(dmmg[j-1]->dm,dmmg[j]->dm,&dmmg[j]->inject);
374: }
375: }
377: KSPGetPC(dmmg[level]->ksp,&pc);
378: mg = ((PC_MG**)pc->data);
379: for (i=0; i<100; i++) {
381: for (j=level; j>0; j--) {
382: PetscPrintf(PETSC_COMM_WORLD,"I am here");
383: /* Relax residual_fine - F(x_fine) = 0 */
384: for (k=0; k<dmmg[j]->presmooth; k++) {
385: NLFRelax_DAAD4(dmmg[j]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[j]->x);
386: }
388: /* R*(residual_fine - F(x_fine)) */
389: DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
390: VecAYPX(dmmg[j]->w,mone,dmmg[j]->r);
392: if (j == level || dmmg[j]->monitorall) {
393: /* norm( residual_fine - f(x_fine) ) */
394: VecNorm(dmmg[j]->w,NORM_2,&norm);
395: if (j == level) {
396: if (norm < dmmg[level]->abstol) goto theend;
397: if (i == 0) {
398: dmmg[level]->rrtol = norm*dmmg[level]->rtol;
399: } else {
400: if (norm < dmmg[level]->rrtol) goto theend;
401: }
402: }
403: } PetscPrintf(PETSC_COMM_WORLD,"I am here");
405: if (dmmg[j]->monitorall) {
406: for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm," ");}
407: PetscPrintf(dmmg[j]->comm,"FAS function norm %g\n",norm);
408: }
409: MatRestrict(mg[j]->restrct,dmmg[j]->w,dmmg[j-1]->r);
410:
411: /* F(R*x_fine) */
412: VecScatterBegin(dmmg[j]->inject,dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD);
413: VecScatterEnd(dmmg[j]->inject,dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD);
414: DMMGFormFunction(0,dmmg[j-1]->x,dmmg[j-1]->w,dmmg[j-1]);
416: /* residual_coarse = F(R*x_fine) + R*(residual_fine - F(x_fine)) */
417: VecAYPX(dmmg[j-1]->r,one,dmmg[j-1]->w);
419: /* save R*x_fine into b (needed when interpolating compute x back up */
420: VecCopy(dmmg[j-1]->x,dmmg[j-1]->b);
421: }
423: for (j=0; j<dmmg[0]->presmooth; j++) {
424: NLFRelax_DAAD4(dmmg[0]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[0]->x);
425: }
426: if (dmmg[0]->monitorall){
427: DMMGFormFunction(0,dmmg[0]->x,dmmg[0]->w,dmmg[0]);
428: VecAXPY(dmmg[0]->w,mone,dmmg[0]->r);
429: VecNorm(dmmg[0]->w,NORM_2,&norm);
430: for (k=0; k<level+1; k++) {PetscPrintf(dmmg[0]->comm," ");}
431: PetscPrintf(dmmg[0]->comm,"FAS coarse grid function norm %g\n",norm);
432: }
434: for (j=1; j<=level; j++) {
435: /* x_fine = x_fine + R'*(x_coarse - R*x_fine) */
436: VecAXPY(dmmg[j-1]->x,mone,dmmg[j-1]->b);
437: MatInterpolateAdd(mg[j]->interpolate,dmmg[j-1]->x,dmmg[j]->x,dmmg[j]->x);
439: if (dmmg[j]->monitorall) {
440: /* norm( F(x_fine) - residual_fine ) */
441: DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
442: VecAXPY(dmmg[j]->w,mone,dmmg[j]->r);
443: VecNorm(dmmg[j]->w,NORM_2,&norm);
444: for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm," ");}
445: PetscPrintf(dmmg[j]->comm,"FAS function norm %g\n",norm);
446: }
448: /* Relax residual_fine - F(x_fine) = 0 */
449: for (k=0; k<dmmg[j]->postsmooth; k++) {
450: NLFRelax_DAAD4(dmmg[j]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[j]->x);
451: }
453: if (dmmg[j]->monitorall) {
454: /* norm( F(x_fine) - residual_fine ) */
455: DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
456: VecAXPY(dmmg[j]->w,mone,dmmg[j]->r);
457: VecNorm(dmmg[j]->w,NORM_2,&norm);
458: for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm," ");}
459: PetscPrintf(dmmg[j]->comm,"FAS function norm %g\n",norm);
460: }
461: }
463: if (dmmg[level]->monitor){
464: DMMGFormFunction(0,dmmg[level]->x,dmmg[level]->w,dmmg[level]);
465: VecNorm(dmmg[level]->w,NORM_2,&norm);
466: PetscPrintf(dmmg[level]->comm,"%D FAS function norm %g\n",i,norm);
467: }
468: }
469: theend:
470: return(0);
471: }
473: /*
474: This function provide several FAS v_cycle iteration
476: iter: the number of FAS it run
478: */
481: PetscErrorCode DMMGSolveFASn(DMMG *dmmg,PetscInt level,PetscInt iter)
482: {
484: PetscInt i,j,k;
485: PetscReal norm;
486: PC_MG **mg;
487: PC pc;
490: VecSet(dmmg[level]->r,0.0);
491: for (j=1; j<=level; j++) {
492: if (!dmmg[j]->inject) {
493: DMGetInjection(dmmg[j-1]->dm,dmmg[j]->dm,&dmmg[j]->inject);
494: }
495: }
497: KSPGetPC(dmmg[level]->ksp,&pc);
498: mg = ((PC_MG**)pc->data);
500: for (i=0; i<iter; i++) {
502: for (j=level; j>0; j--) {
504: /* Relax residual_fine - F(x_fine) = 0 */
505: for (k=0; k<dmmg[j]->presmooth; k++) {
506: NLFRelax_DAAD(dmmg[j]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[j]->x);
507: }
509: /* R*(residual_fine - F(x_fine)) */
510: DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
511: VecAYPX(dmmg[j]->w,-1.0,dmmg[j]->r);
513: if (j == level || dmmg[j]->monitorall) {
514: /* norm( residual_fine - f(x_fine) ) */
515: VecNorm(dmmg[j]->w,NORM_2,&norm);
516: if (j == level) {
517: if (norm < dmmg[level]->abstol) goto theend;
518: if (i == 0) {
519: dmmg[level]->rrtol = norm*dmmg[level]->rtol;
520: } else {
521: if (norm < dmmg[level]->rrtol) goto theend;
522: }
523: }
524: }
526: if (dmmg[j]->monitorall) {
527: for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm," ");}
528: PetscPrintf(dmmg[j]->comm,"FAS function norm %g\n",norm);
529: }
530: MatRestrict(mg[j]->restrct,dmmg[j]->w,dmmg[j-1]->r);
531:
532: /* F(RI*x_fine) */
533: VecScatterBegin(dmmg[j]->inject,dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD);
534: VecScatterEnd(dmmg[j]->inject,dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD);
535: DMMGFormFunction(0,dmmg[j-1]->x,dmmg[j-1]->w,dmmg[j-1]);
537: /* residual_coarse = F(RI*x_fine) + R*(residual_fine - F(x_fine)) */
538: VecAYPX(dmmg[j-1]->r,1.0,dmmg[j-1]->w);
540: /* save RI*x_fine into b (needed when interpolating compute x back up */
541: VecCopy(dmmg[j-1]->x,dmmg[j-1]->b);
542: }
544: for (j=0; j<dmmg[0]->coarsesmooth; j++) {
545: NLFRelax_DAAD(dmmg[0]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[0]->x);
546: }
547: if (dmmg[0]->monitorall){
548: DMMGFormFunction(0,dmmg[0]->x,dmmg[0]->w,dmmg[0]);
549: VecAXPY(dmmg[0]->w,-1.0,dmmg[0]->r);
550: VecNorm(dmmg[0]->w,NORM_2,&norm);
551: for (k=0; k<level+1; k++) {PetscPrintf(dmmg[0]->comm," ");}
552: PetscPrintf(dmmg[0]->comm,"FAS coarse grid function norm %g\n",norm);
553: }
555: for (j=1; j<=level; j++) {
556: /* x_fine = x_fine + R'*(x_coarse - RI*x_fine) */
557: VecAXPY(dmmg[j-1]->x,-1.0,dmmg[j-1]->b);
558: MatInterpolateAdd(mg[j]->interpolate,dmmg[j-1]->x,dmmg[j]->x,dmmg[j]->x);
560: if (dmmg[j]->monitorall) {
561: /* norm( F(x_fine) - residual_fine ) */
562: DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
563: VecAXPY(dmmg[j]->w,-1.0,dmmg[j]->r);
564: VecNorm(dmmg[j]->w,NORM_2,&norm);
565: for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm," ");}
566: PetscPrintf(dmmg[j]->comm,"FAS function norm %g\n",norm);
567: }
569: /* Relax residual_fine - F(x_fine) = 0 */
570: for (k=0; k<dmmg[j]->postsmooth; k++) {
571: NLFRelax_DAAD(dmmg[j]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[j]->x);
572: }
574: if (dmmg[j]->monitorall) {
575: /* norm( F(x_fine) - residual_fine ) */
576: DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
577: VecAXPY(dmmg[j]->w,-1.0,dmmg[j]->r);
578: VecNorm(dmmg[j]->w,NORM_2,&norm);
579: for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm," ");}
580: PetscPrintf(dmmg[j]->comm,"FAS function norm %g\n",norm);
581: }
582: }
584: if (dmmg[level]->monitor){
585: DMMGFormFunction(0,dmmg[level]->x,dmmg[level]->w,dmmg[level]);
586: VecNorm(dmmg[level]->w,NORM_2,&norm);
587: PetscPrintf(dmmg[level]->comm,"%D FAS function norm %g\n",i+1,norm);
588: }
589: }
590: theend:
591: return(0);
592: }
593: /*
594: This is a simple FAS setup function
595: */
600: PetscErrorCode DMMGSolveFASSetUp(DMMG *dmmg,PetscInt level)
601: {
603: PetscInt j;//,nlevels=dmmg[0]->nlevels-1;
604: //PC pc;
607: VecSet(dmmg[level]->r,0.0);
608: for (j=1; j<=level; j++) {
609: if (!dmmg[j]->inject) {
610: DMGetInjection(dmmg[j-1]->dm,dmmg[j]->dm,&dmmg[j]->inject);
611: }
612: }
613: VecSet(dmmg[level]->r,0.0);
614: dmmg[level]->rrtol = 0.0001*dmmg[level]->rtol;//I want to get more precise solution with FAS
615: return(0);
616: }
619: /*
620: This is function to implement multiplicative FAS
623: Options:
624:
625: -dmmg_fas_cycles 1 : FAS v-cycle
626: 2 : FAS w-cycle
629: */
633: PetscErrorCode DMMGSolveFASMCycle(DMMG *dmmg,PetscInt level,PetscTruth* converged)
634: {
636: PetscInt j,k,cycles=1,nlevels=level;//nlevels=dmmg[0]->nlevels-1;
637: // I need to put nlevels=level in order to get grid sequence correctly
638: PetscReal norm;
639: PC_MG **mg;
640: PC pc;
641:
643:
645: PetscOptionsGetInt(PETSC_NULL,"-dmmg_fas_cycles",&cycles,PETSC_NULL);
646:
647: KSPGetPC(dmmg[level]->ksp,&pc);
648: mg = ((PC_MG**)pc->data);
650: j=level;
652: if(j) {/* not the coarsest level */
653: /* Relax residual_fine - F(x_fine) = 0 */
654: for (k=0; k<dmmg[j]->presmooth; k++) {
655: NLFRelax_DAAD(dmmg[j]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[j]->x);
656: }
657:
658:
660: /* R*(residual_fine - F(x_fine)) */
661: DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
662: VecAYPX(dmmg[j]->w,-1.0,dmmg[j]->r);
664: if (j == nlevels || dmmg[j]->monitorall) {
665: /* norm( residual_fine - f(x_fine) ) */
666: VecNorm(dmmg[j]->w,NORM_2,&norm);
667:
668: if (j == nlevels) {
669: if (norm < dmmg[level]->abstol) {
670: *converged = PETSC_TRUE;
671: goto theend;
672: }
674: if (norm < dmmg[level]->rrtol){
675: *converged = PETSC_TRUE;
676: goto theend;
677:
678: }
679: }
680: }
682: if (dmmg[j]->monitorall) {
683: for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm," ");}
684: PetscPrintf(dmmg[j]->comm,"FAS function norm %g\n",norm);
685: }
686: MatRestrict(mg[j]->restrct,dmmg[j]->w,dmmg[j-1]->r);
687:
688: /* F(RI*x_fine) */
689: VecScatterBegin(dmmg[j]->inject,dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD);
690: VecScatterEnd(dmmg[j]->inject,dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD);
691: DMMGFormFunction(0,dmmg[j-1]->x,dmmg[j-1]->w,dmmg[j-1]);
693: /* residual_coarse = F(RI*x_fine) + R*(residual_fine - F(x_fine)) */
694: VecAYPX(dmmg[j-1]->r,1.0,dmmg[j-1]->w);
696: /* save RI*x_fine into b (needed when interpolating compute x back up */
697: VecCopy(dmmg[j-1]->x,dmmg[j-1]->b);
698:
699:
700: while (cycles--) {
701: DMMGSolveFASMCycle(dmmg,level-1,converged);
702: }
703: }
704: else { /* for the coarsest level */
705: for (k=0; k<dmmg[0]->coarsesmooth; k++) {
706: NLFRelax_DAAD(dmmg[0]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[0]->x);
707: }
708: if (dmmg[0]->monitorall){
709: DMMGFormFunction(0,dmmg[0]->x,dmmg[0]->w,dmmg[0]);
710: VecAXPY(dmmg[0]->w,-1.0,dmmg[0]->r);
711: VecNorm(dmmg[0]->w,NORM_2,&norm);
712: for (k=0; k<level+1; k++) {PetscPrintf(dmmg[0]->comm," ");}
713: PetscPrintf(dmmg[0]->comm,"FAS coarse grid function norm %g\n",norm);
714: }
715: if (j == nlevels || dmmg[j]->monitorall) {
716: /* norm( residual_fine - f(x_fine) ) */
717: VecNorm(dmmg[j]->w,NORM_2,&norm);
718:
719: if (j == nlevels) {
720: if (norm < dmmg[level]->abstol) {
721: *converged = PETSC_TRUE;
722: goto theend;
723: }
724:
725: if (norm < dmmg[level]->rrtol){
726: *converged = PETSC_TRUE;
727: goto theend;
728:
729: }
730: }
731: }
733:
734: }
735: j=level;
736: if(j) { /* not for the coarsest level */
737: /* x_fine = x_fine + R'*(x_coarse - RI*x_fine) */
738: VecAXPY(dmmg[j-1]->x,-1.0,dmmg[j-1]->b);
739: MatInterpolateAdd(mg[j]->interpolate,dmmg[j-1]->x,dmmg[j]->x,dmmg[j]->x);
741: if (dmmg[j]->monitorall) {
742: /* norm( F(x_fine) - residual_fine ) */
743: DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
744: VecAXPY(dmmg[j]->w,-1.0,dmmg[j]->r);
745: VecNorm(dmmg[j]->w,NORM_2,&norm);
746: for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm," ");}
747: PetscPrintf(dmmg[j]->comm,"FAS function norm %g\n",norm);
748: }
750: /* Relax residual_fine - F(x_fine) = 0 */
751: for (k=0; k<dmmg[j]->postsmooth; k++) {
752: NLFRelax_DAAD(dmmg[j]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[j]->x);
753: }
755: if (dmmg[j]->monitorall) {
756: /* norm( F(x_fine) - residual_fine ) */
757: DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
758: VecAXPY(dmmg[j]->w,-1.0,dmmg[j]->r);
759: VecNorm(dmmg[j]->w,NORM_2,&norm);
760: for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm," ");}
761: PetscPrintf(dmmg[j]->comm,"FAS function norm %g\n",norm);
762: }
763:
764:
765: }
766: theend:
767: return(0);
768: }
770: /*
771: This is function to implement multiplicative FAS with block smoother
774: Options:
775:
776: -dmmg_fas_cycles 1 : FAS v-cycle
777: 2 : FAS w-cycle
780: */
785: PetscErrorCode DMMGSolveFASMCycle9(DMMG *dmmg,PetscInt level,PetscTruth* converged)
786: {
788: PetscInt j,k,cycles=1,nlevels=level;//nlevels=dmmg[0]->nlevels-1;
789: // I need to put nlevels=level in order to get grid sequence correctly
790: PetscReal norm;
791: PC_MG **mg;
792: PC pc;
793:
795:
796:
797: PetscOptionsGetInt(PETSC_NULL,"-dmmg_fas_cycles",&cycles,PETSC_NULL);
798:
799: KSPGetPC(dmmg[level]->ksp,&pc);
800: mg = ((PC_MG**)pc->data);
801: // for (j=level; j>0; j--) {
802: j=level;
803: //PetscPrintf(dmmg[level]->comm,"j=%d,nlevels=%d",j,nlevels);
804: if(j) {/* not the coarsest level */
805: /* Relax residual_fine - F(x_fine) = 0 */
806: for (k=0; k<dmmg[j]->presmooth; k++) {
807: NLFRelax_DAAD9(dmmg[j]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[j]->x);
808: }
809:
810:
812: /* R*(residual_fine - F(x_fine)) */
813: DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
814: VecAYPX(dmmg[j]->w,-1.0,dmmg[j]->r);
816: if (j == nlevels || dmmg[j]->monitorall) {
817: /* norm( residual_fine - f(x_fine) ) */
818: VecNorm(dmmg[j]->w,NORM_2,&norm);
819:
820: if (j == nlevels) {
821: if (norm < dmmg[level]->abstol) {
822: *converged = PETSC_TRUE;
823: goto theend;
824: }
825: /* if (i == 0) {
826: dmmg[level]->rrtol = norm*dmmg[level]->rtol;
827: } else {*/
828: if (norm < dmmg[level]->rrtol){
829: *converged = PETSC_TRUE;
830: goto theend;
831:
832: }
833: }
834: }
836: if (dmmg[j]->monitorall) {
837: for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm," ");}
838: PetscPrintf(dmmg[j]->comm,"FAS function norm %g\n",norm);
839: }
840: MatRestrict(mg[j]->restrct,dmmg[j]->w,dmmg[j-1]->r);
841:
842: /* F(RI*x_fine) */
843: VecScatterBegin(dmmg[j]->inject,dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD);
844: VecScatterEnd(dmmg[j]->inject,dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD);
845: DMMGFormFunction(0,dmmg[j-1]->x,dmmg[j-1]->w,dmmg[j-1]);
847: /* residual_coarse = F(RI*x_fine) + R*(residual_fine - F(x_fine)) */
848: VecAYPX(dmmg[j-1]->r,1.0,dmmg[j-1]->w);
850: /* save RI*x_fine into b (needed when interpolating compute x back up */
851: VecCopy(dmmg[j-1]->x,dmmg[j-1]->b);
852:
853:
854: while (cycles--) {
855: DMMGSolveFASMCycle9(dmmg,level-1,converged);
856: }
857: }
858: else { /* for the coarsest level */
859: for (k=0; k<dmmg[0]->coarsesmooth; k++) {
860: NLFRelax_DAAD9(dmmg[0]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[0]->x);
861: }
862: if (dmmg[0]->monitorall){
863: DMMGFormFunction(0,dmmg[0]->x,dmmg[0]->w,dmmg[0]);
864: VecAXPY(dmmg[0]->w,-1.0,dmmg[0]->r);
865: VecNorm(dmmg[0]->w,NORM_2,&norm);
866: for (k=0; k<level+1; k++) {PetscPrintf(dmmg[0]->comm," ");}
867: PetscPrintf(dmmg[0]->comm,"FAS coarse grid function norm %g\n",norm);
868: }
869: if (j == nlevels || dmmg[j]->monitorall) {
870: /* norm( residual_fine - f(x_fine) ) */
871: VecNorm(dmmg[j]->w,NORM_2,&norm);
872:
873: if (j == nlevels) {
874: if (norm < dmmg[level]->abstol) {
875: *converged = PETSC_TRUE;
876: goto theend;
877: }
878: /* if (i == 0) {
879: dmmg[level]->rrtol = norm*dmmg[level]->rtol;
880: } else {*/
881: if (norm < dmmg[level]->rrtol){
882: *converged = PETSC_TRUE;
883: goto theend;
884:
885: }
886: }
887: }
889:
890: }
891: j=level;
892: if(j) { /* not for the coarsest level */
893: /* x_fine = x_fine + R'*(x_coarse - RI*x_fine) */
894: VecAXPY(dmmg[j-1]->x,-1.0,dmmg[j-1]->b);
895: MatInterpolateAdd(mg[j]->interpolate,dmmg[j-1]->x,dmmg[j]->x,dmmg[j]->x);
897: if (dmmg[j]->monitorall) {
898: /* norm( F(x_fine) - residual_fine ) */
899: DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
900: VecAXPY(dmmg[j]->w,-1.0,dmmg[j]->r);
901: VecNorm(dmmg[j]->w,NORM_2,&norm);
902: for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm," ");}
903: PetscPrintf(dmmg[j]->comm,"FAS function norm %g\n",norm);
904: }
906: /* Relax residual_fine - F(x_fine) = 0 */
907: for (k=0; k<dmmg[j]->postsmooth; k++) {
908: NLFRelax_DAAD9(dmmg[j]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[j]->x);
909: }
911: if (dmmg[j]->monitorall) {
912: /* norm( F(x_fine) - residual_fine ) */
913: DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
914: VecAXPY(dmmg[j]->w,-1.0,dmmg[j]->r);
915: VecNorm(dmmg[j]->w,NORM_2,&norm);
916: for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm," ");}
917: PetscPrintf(dmmg[j]->comm,"FAS function norm %g\n",norm);
918: }
919:
920: /* if(j==nlevels){
921: if (dmmg[level]->monitor){
922: DMMGFormFunction(0,dmmg[level]->x,dmmg[level]->w,dmmg[level]);
923: VecNorm(dmmg[level]->w,NORM_2,&norm);
924:
925: PetscPrintf(dmmg[level]->comm," FAS function norm %g\n",norm);
926:
927: }
928: }*/
929: }
930: theend:
931: return(0);
932: }
934: /*
935: This is function to implement full FAS with block smoother(9 points together)
938: Options:
939:
940: -dmmg_fas_cycles 1 : FAS v-cycle
941: 2 : FAS w-cycle
944: */
949: PetscErrorCode DMMGSolveFASFCycle(DMMG *dmmg,PetscInt l,PetscTruth* converged)
950: {
952: PetscInt j;//l = dmmg[0]->nlevels-1;
953: PC_MG **mg;
954: PC pc;
957: KSPGetPC(dmmg[l]->ksp,&pc);
958: mg = ((PC_MG**)pc->data);
959: // restriction all the way down to the coarse level
960: if(l>0) {
961: for(j=l;j>0;j--) {
962:
963: /* R*(residual_fine - F(x_fine)) */
964: DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
965: VecAYPX(dmmg[j]->w,-1.0,dmmg[j]->r);
967: MatRestrict(mg[j]->restrct,dmmg[j]->w,dmmg[j-1]->r);
968:
969: /* F(RI*x_fine) */
970: VecScatterBegin(dmmg[j]->inject,dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD);
971: VecScatterEnd(dmmg[j]->inject,dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD);
972: DMMGFormFunction(0,dmmg[j-1]->x,dmmg[j-1]->w,dmmg[j-1]);
974: /* residual_coarse = F(RI*x_fine) + R*(residual_fine - F(x_fine)) */
975: VecAYPX(dmmg[j-1]->r,1.0,dmmg[j-1]->w);
977: /* save RI*x_fine into b (needed when interpolating compute x back up */
978: VecCopy(dmmg[j-1]->x,dmmg[j-1]->b);
980: }
982: // all the way up to the finest level
983: for (j=0; j<l; j++) {
984: DMMGSolveFASMCycle(dmmg,j,PETSC_NULL);
985: /* x_fine = x_fine + R'*(x_coarse - RI*x_fine) */
986: VecAXPY(dmmg[j]->x,-1.0,dmmg[j]->b);
987: MatInterpolateAdd(mg[j+1]->interpolate,dmmg[j]->x,dmmg[j+1]->x,dmmg[j+1]->x);
989: }
990: }
991: DMMGSolveFASMCycle(dmmg,l,converged);
992: return(0);
993: }
997: /*
998: This is function to implement full FAS with block smoother ( 9 points together)
1001: Options:
1002:
1003: -dmmg_fas_cycles 1 : FAS v-cycle
1004: 2 : FAS w-cycle
1007: */
1010: PetscErrorCode DMMGSolveFASFCycle9(DMMG *dmmg,PetscInt l,PetscTruth* converged)
1011: {
1013: PetscInt j;//l = dmmg[0]->nlevels-1;
1014: PC_MG **mg;
1015: PC pc;
1018: KSPGetPC(dmmg[l]->ksp,&pc);
1019: mg = ((PC_MG**)pc->data);
1020: // restriction all the way down to the coarse level
1021: if(l>0) {
1022: for(j=l;j>0;j--) {
1023:
1024: /* R*(residual_fine - F(x_fine)) */
1025: DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
1026: VecAYPX(dmmg[j]->w,-1.0,dmmg[j]->r);
1028: MatRestrict(mg[j]->restrct,dmmg[j]->w,dmmg[j-1]->r);
1029:
1030: /* F(RI*x_fine) */
1031: VecScatterBegin(dmmg[j]->inject,dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD);
1032: VecScatterEnd(dmmg[j]->inject,dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD);
1033: DMMGFormFunction(0,dmmg[j-1]->x,dmmg[j-1]->w,dmmg[j-1]);
1035: /* residual_coarse = F(RI*x_fine) + R*(residual_fine - F(x_fine)) */
1036: VecAYPX(dmmg[j-1]->r,1.0,dmmg[j-1]->w);
1038: /* save RI*x_fine into b (needed when interpolating compute x back up */
1039: VecCopy(dmmg[j-1]->x,dmmg[j-1]->b);
1041: }
1043: // all the way up to the finest level
1044: for (j=0; j<l; j++) {
1045: DMMGSolveFASMCycle9(dmmg,j,PETSC_NULL);
1046: /* x_fine = x_fine + R'*(x_coarse - RI*x_fine) */
1047: VecAXPY(dmmg[j]->x,-1.0,dmmg[j]->b);
1048: MatInterpolateAdd(mg[j+1]->interpolate,dmmg[j]->x,dmmg[j+1]->x,dmmg[j+1]->x);
1050: }
1051: }
1052: DMMGSolveFASMCycle9(dmmg,l,converged);
1053: return(0);
1054: }
1056: /*
1057: This is function is to solve nonlinear system with FAS
1059: Options:
1061: -dmmg_fas_9: using block smoother
1062:
1063: -dmmg_fas_full: using full FAS
1066: */
1069: PetscErrorCode DMMGSolveFASCycle(DMMG *dmmg,PetscInt level)
1070: {
1072: PetscInt i;
1073: PetscTruth converged = PETSC_FALSE, flg,flgb;
1074: PetscReal norm;
1077: DMMGSolveFASSetUp(dmmg,level);
1078: PetscOptionsHasName(PETSC_NULL,"-dmmg_fas_9",&flgb);
1079:
1080: if(flgb){
1082: PetscOptionsHasName(PETSC_NULL,"-dmmg_fas_full",&flg);
1083: if (flg) {
1084: for(i=0;i<1000;i++){
1085: PetscPrintf(dmmg[level]->comm,"%D ",i+1);
1086: DMMGSolveFASFCycle9(dmmg,level,&converged);
1087: if (dmmg[level]->monitor){
1088: DMMGFormFunction(0,dmmg[level]->x,dmmg[level]->w,dmmg[level]);
1089: VecNorm(dmmg[level]->w,NORM_2,&norm);
1090: PetscPrintf(dmmg[level]->comm," FAS function norm %g\n",norm);
1091: }
1092: if (converged) return(0);
1093: }
1094: }
1095: else{
1096: for(i=0;i<1000;i++){
1097: PetscPrintf(dmmg[level]->comm,"%D ",i+1);
1098: DMMGSolveFASMCycle9(dmmg,level,&converged);
1099: if (dmmg[level]->monitor){
1100: DMMGFormFunction(0,dmmg[level]->x,dmmg[level]->w,dmmg[level]);
1101: VecNorm(dmmg[level]->w,NORM_2,&norm);
1102: PetscPrintf(dmmg[level]->comm," FAS function norm %g\n",norm);
1103: }
1105: if (converged) return(0);
1106: }
1107: }
1108: }
1109: else {
1110: PetscOptionsHasName(PETSC_NULL,"-dmmg_fas_full",&flg);
1111: if (flg) {
1112: for(i=0;i<1000;i++){
1113: PetscPrintf(dmmg[level]->comm,"%D ",i+1);
1114: DMMGSolveFASFCycle(dmmg,level,&converged);
1115: if (dmmg[level]->monitor){
1116: DMMGFormFunction(0,dmmg[level]->x,dmmg[level]->w,dmmg[level]);
1117: VecNorm(dmmg[level]->w,NORM_2,&norm);
1118: PetscPrintf(dmmg[level]->comm," FAS function norm %g\n",norm);
1119: }
1120: if (converged) return(0);
1121: }
1122: }
1123: else{
1124: for(i=0;i<1000;i++){
1125: PetscPrintf(dmmg[level]->comm,"%D ",i+1);
1126: DMMGSolveFASMCycle(dmmg,level,&converged);
1127: if (dmmg[level]->monitor){
1128: DMMGFormFunction(0,dmmg[level]->x,dmmg[level]->w,dmmg[level]);
1129: VecNorm(dmmg[level]->w,NORM_2,&norm);
1130: PetscPrintf(dmmg[level]->comm," FAS function norm %g\n",norm);
1131: }
1133: if (converged) return(0);
1134: }
1136: }
1137: }
1138: return(0);
1139: }
1141: /*
1142: This is function is to implement one FAS iteration
1144: Options:
1146: -dmmg_fas_9: using block smoother
1147:
1148: -dmmg_fas_full: using full FAS
1150: */
1153: PetscErrorCode DMMGSolveFASCyclen(DMMG *dmmg,PetscInt level)
1154: {
1156: PetscTruth converged = PETSC_FALSE, flg,flgb;
1157: PetscReal norm;
1158: // PC_MG **mg;
1159: //PC pc;
1162: DMMGSolveFASSetUp(dmmg,level);
1163: PetscOptionsHasName(PETSC_NULL,"-dmmg_fas_9",&flgb);
1164:
1165: if(flgb){
1167: PetscOptionsHasName(PETSC_NULL,"-dmmg_fas_full",&flg);
1168: if (flg) {
1169:
1170: DMMGSolveFASFCycle9(dmmg,level,&converged);
1171: if (dmmg[level]->monitor){
1172: DMMGFormFunction(0,dmmg[level]->x,dmmg[level]->w,dmmg[level]);
1173: VecNorm(dmmg[level]->w,NORM_2,&norm);
1174: PetscPrintf(dmmg[level]->comm," FAS function norm %g\n",norm);
1175: }
1177: }
1178: else{
1179:
1180: DMMGSolveFASMCycle9(dmmg,level,&converged);
1181: if (dmmg[level]->monitor){
1182: DMMGFormFunction(0,dmmg[level]->x,dmmg[level]->w,dmmg[level]);
1183: VecNorm(dmmg[level]->w,NORM_2,&norm);
1184: PetscPrintf(dmmg[level]->comm," FAS function norm %g\n",norm);
1185: }
1188: }
1189: }
1190: else {
1191: PetscOptionsHasName(PETSC_NULL,"-dmmg_fas_full",&flg);
1192: if (flg) {
1193:
1194: DMMGSolveFASFCycle(dmmg,level,&converged);
1195: if (dmmg[level]->monitor){
1196: DMMGFormFunction(0,dmmg[level]->x,dmmg[level]->w,dmmg[level]);
1197: VecNorm(dmmg[level]->w,NORM_2,&norm);
1198: PetscPrintf(dmmg[level]->comm," FAS function norm %g\n",norm);
1199: }
1201: }
1202: else{
1203:
1204: DMMGSolveFASMCycle(dmmg,level,&converged);
1205: if (dmmg[level]->monitor){
1206: DMMGFormFunction(0,dmmg[level]->x,dmmg[level]->w,dmmg[level]);
1207: VecNorm(dmmg[level]->w,NORM_2,&norm);
1208: PetscPrintf(dmmg[level]->comm," FAS function norm %g\n",norm);
1209: }
1213: }
1214: }
1215:
1217: return(0);
1218: }
1221: /*
1223: This is function to implement Nonlinear CG to accelerate FAS
1225: In order to use this acceleration, the option is
1227: -dmmg_fas_NCG
1231: */
1236: PetscErrorCode DMMGSolveFAS_NCG(DMMG *dmmg, PetscInt level)
1237: {
1238: SNES snes = dmmg[level]->snes;
1239: SNES_LS *neP = (SNES_LS*)snes->data;
1241: PetscInt maxits,i,lits;
1242: PetscTruth lssucceed;
1243: // MatStructure flg = DIFFERENT_NONZERO_PATTERN;
1244: PetscReal fnorm,gnorm,xnorm,ynorm,betaFR,betaPR,beta,betaHS,betaDY;
1245: Vec Y,X,F,G,W,Gradold,Sk;
1246: //KSP ksp;
1251: if (!snes->solve) SETERRQ(PETSC_ERR_ORDER,"SNESSetType() or SNESSetFromOptions() must be called before SNESSolve()");
1253: VecDuplicate(dmmg[level]->x,&Sk);
1254: snes->vec_sol = snes->vec_sol_always = dmmg[level]->x;
1255: if (!snes->setupcalled) {
1256: SNESSetUp(snes);
1257: }
1258: if (snes->conv_hist_reset) snes->conv_hist_len = 0;
1260: snes->nfuncs = 0; snes->linear_its = 0; snes->numFailures = 0;
1263: snes->reason = SNES_CONVERGED_ITERATING;
1265: maxits = snes->max_its; /* maximum number of iterations */
1266: X = snes->vec_sol; /* solution vector */
1267: F = snes->vec_func; /* residual vector */
1268: Y = snes->work[0]; /* work vectors */
1269: G = snes->work[1];
1270: W = snes->work[2];
1272: PetscObjectTakeAccess(snes);
1273: snes->iter = 0;
1274: PetscObjectGrantAccess(snes);
1275:
1276: X = dmmg[level]->x;
1277: VecCopy(X,Y);
1278: VecCopy(X,G);
1280: // to get the residual for the F
1281: SNESComputeFunction(snes,X,F);
1282:
1283: VecNorm(F,NORM_2,&fnorm); /* fnorm <- ||F|| */
1284: if (fnorm != fnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
1285: PetscObjectTakeAccess(snes);
1286: snes->norm = fnorm;
1287: PetscObjectGrantAccess(snes);
1288: SNESLogConvHistory(snes,fnorm,0);
1289: SNESMonitor(snes,0,fnorm);
1291: if (fnorm < snes->abstol) {snes->reason = SNES_CONVERGED_FNORM_ABS; return(0);}
1293: /* set parameter for default relative tolerance convergence test */
1294: snes->ttol = fnorm*snes->rtol;
1295: // dmmg[level]->rrtol= snes->ttol;
1297: // set this to store the old grad
1298: Gradold=snes->vec_sol_update_always;
1299:
1300: // compute the search direction Y
1301: // I need to put Q(x)=x-FAS(x) here
1302: DMMGSolveFASCyclen(dmmg,level);
1303: // F = X - dmmg[level]->x; this is the gradient direction
1304: VecAXPY(Y,-1.0,X);
1305: // copy the gradient to the old
1306: VecCopy(Y,Gradold);
1307: // copy X back
1308: VecCopy(G,X);
1309: VecWAXPY(Sk,-1.0,X,X);
1310: // so far I put X=X_c, F= F(x_c), Gradold= Y=grad(x_c)
1312: // for (i=0; i<maxits; i++) {
1314: for (i=0; i<10000; i++) {
1317: PetscPrintf(PETSC_COMM_WORLD,"iter=%d",i+1);
1318:
1319:
1320: // X=x_c, F=F(x_c),Y search direction; G=F(x_new), W=x_new,
1321: VecNorm(X,NORM_2,&xnorm);
1322: (*neP->LineSearch)(snes,neP->lsP,X,F,G,Y,W,fnorm,&ynorm,&gnorm,&lssucceed);
1323: PetscInfo4(snes,"SNESSolve_LS: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",fnorm,gnorm,ynorm,(int)lssucceed);
1324: if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
1325: PetscPrintf(PETSC_COMM_WORLD,"step=%g,oldnorm=%g,norm=%g ",ynorm,fnorm,gnorm);
1326:
1327: fnorm=gnorm; //copy the new function norm; this will effect the line_search
1328: VecWAXPY(Sk,-1.0,X,W);
1329: //update the new solution
1330: ierr=VecCopy(W,X);
1331: ierr=VecCopy(G,F);
1333:
1334: // compute the new search direction G
1335: // I need to put Q(x)=x-FAS(x) here
1337: DMMGSolveFASCyclen(dmmg,level);
1338: //G = X - dmmg[level]->x; G is the new gradient, Y is old gradient
1339:
1340: VecWAXPY(G,-1.0,X,W);
1341: // copy W back to X
1342: VecCopy(W,X);
1343: VecNorm(G,NORM_2,&gnorm);
1344: VecNorm(Gradold,NORM_2,&ynorm);
1345: betaFR = gnorm*gnorm/(ynorm*ynorm); //FR_beta
1346:
1347: VecWAXPY(W,-1.0,Gradold,G);
1348: VecDot(W,G,&gnorm);
1349: // VecNorm(G,NORM_2,&gnorm);
1350: VecNorm(Gradold,NORM_2,&ynorm);
1351: betaPR = gnorm/(ynorm*ynorm); //PR_beta
1352:
1353: if ( betaPR<-betaFR)
1354: {
1355: beta =- betaFR;
1356: }