Actual source code: nladic.c
petsc-3.3-p7 2013-05-11
2: /*
3: ADIC based nonlinear operator object that can be used with FAS
5: This does not really belong in the matrix directories but since it
6: was cloned off of Mat_DAAD I'm leaving it here until I have a better place
8: */
9: #include <petscsys.h>
10: #include <petscdm.h>
12: EXTERN_C_BEGIN
13: #include <adic/ad_utils.h>
14: EXTERN_C_END
16: #include <petsc-private/daimpl.h>
17: #include <../src/mat/blockinvert.h>
19: struct NLF_DAAD {
20: DM da;
21: void *ctx;
22: Vec residual;
23: int newton_its;
24: };
26: /*
27: Solves the one dimensional equation using Newton's method
28: */
31: PetscErrorCode NLFNewton_DAAD(NLF A,DMDALocalInfo *info,MatStencil *stencil,void *ad_vu,PetscScalar *ad_vustart,int nI,int gI,PetscScalar residual)
32: {
34: PetscInt cnt = A->newton_its;
35: PetscScalar ad_f[2],J,f;
38: ad_vustart[1+2*gI] = 1.0;
39: do {
40: /* compute the function and Jacobian */
41: (*A->da->adicmf_lfi)(info,stencil,ad_vu,ad_f,A->ctx);
42: J = -ad_f[1];
43: f = -ad_f[0] + residual;
44: if (f != f) SETERRQ(PETSC_COMM_SELF,1,"nan");
45: ad_vustart[2*gI] = ad_vustart[2*gI] - f/J;
46: } while (--cnt > 0 && PetscAbsScalar(f) > 1.e-14);
48: ad_vustart[1+2*gI] = 0.0;
49: return(0);
50: }
52: /*
53: Solves the four dimensional equation using Newton's method
54: */
57: PetscErrorCode NLFNewton_DAAD4(NLF A,DMDALocalInfo *info,MatStencil *stencil,void *ad_vu,PetscScalar *ad_vustart,int nI,int gI,PetscScalar *residual)
58: {
60: PetscInt cnt = A->newton_its;
61: PetscScalar ad_f[20], J[16],f[4], res, dd[5];
65: /* This sets the identity as the seed matrix for ADIC */
66: CHKMEMQ;
67: ad_vustart[1+5*gI ] = 1.0;
68: CHKMEMQ;
69: ad_vustart[2+5*gI+5 ] = 1.0;
70: CHKMEMQ;
71: ad_vustart[3+5*gI+10] = 1.0;
72: CHKMEMQ;
73: ad_vustart[4+5*gI+15] = 1.0;
74: CHKMEMQ;
76: do {
77: /* compute the function and Jacobian */
78: CHKMEMQ;
79: (*A->da->adicmf_lfib)(info,stencil,ad_vu,ad_f,A->ctx);
80: CHKMEMQ;
81: /* copy ADIC formated Jacobian into regular C array */
82: J[0] = ad_f[1] ; J[1] = ad_f[2] ; J[2] = ad_f[3] ; J[3] = ad_f[4] ;
83: J[4] = ad_f[6] ; J[5] = ad_f[7] ; J[6] = ad_f[8] ; J[7] = ad_f[9] ;
84: J[8] = ad_f[11]; J[9] = ad_f[12]; J[10]= ad_f[13]; J[11]= ad_f[14];
85: J[12]= ad_f[16]; J[13]= ad_f[17]; J[14]= ad_f[18]; J[15]= ad_f[19];
86: CHKMEMQ;
87: f[0] = -ad_f[0] + residual[0];
88: f[1] = -ad_f[5] + residual[1];
89: f[2] = -ad_f[10] + residual[2];
90: f[3] = -ad_f[15] + residual[3];
92: /* solve Jacobian * dd = ff */
94: /* could use PETSc kernel code to solve system with pivoting */
96: /* could put code in here to compute the solution directly using ADIC data structures instead of copying first */
97: dd[0]=J[0]*(J[5]*(J[10]*J[15]-J[11]*J[14])-J[6]*(J[9]*J[15]-J[11]*J[13])+J[7]*(J[9]*J[14]-J[10]*J[13]))-
98: J[1]*(J[4]*(J[10]*J[15]-J[11]*J[14])-J[6]*(J[8]*J[15]-J[11]*J[12])+J[7]*(J[8]*J[14]-J[10]*J[12]))+
99: J[2]*(J[4]*(J[ 9]*J[15]-J[11]*J[13])-J[5]*(J[8]*J[15]-J[11]*J[12])+J[7]*(J[8]*J[13]-J[ 9]*J[12]))-
100: J[3]*(J[4]*(J[ 9]*J[14]-J[10]*J[13])-J[5]*(J[8]*J[14]-J[10]*J[12])+J[6]*(J[8]*J[13]-J[ 9]*J[12]));
102: dd[1]=(f[0]*(J[5]*(J[10]*J[15]-J[11]*J[14])-J[6]*(J[9]*J[15]-J[11]*J[13])+J[7]*(J[9]*J[14]-J[10]*J[13]))-
103: J[1]*(f[1]*(J[10]*J[15]-J[11]*J[14])-J[6]*(f[2]*J[15]-J[11]*f[ 3])+J[7]*(f[2]*J[14]-J[10]*f[ 3]))+
104: J[2]*(f[1]*(J[ 9]*J[15]-J[11]*J[13])-J[5]*(f[2]*J[15]-J[11]*f[ 3])+J[7]*(f[2]*J[13]-J[ 9]*f[ 3]))-
105: J[3]*(f[1]*(J[ 9]*J[14]-J[10]*J[13])-J[5]*(f[2]*J[14]-J[10]*f[ 3])+J[6]*(f[2]*J[13]-J[ 9]*f[ 3])))/dd[0];
107: dd[2]=(J[0]*(f[1]*(J[10]*J[15]-J[11]*J[14])-J[6]*(f[2]*J[15]-J[11]*f[ 3])+J[7]*(f[2]*J[14]-J[10]*f[ 3]))-
108: f[0]*(J[4]*(J[10]*J[15]-J[11]*J[14])-J[6]*(J[8]*J[15]-J[11]*J[12])+J[7]*(J[8]*J[14]-J[10]*J[12]))+
109: J[2]*(J[4]*(f[ 2]*J[15]-J[11]*f[ 3])-f[1]*(J[8]*J[15]-J[11]*J[12])+J[7]*(J[8]*f[ 3]-f[ 2]*J[12]))-
110: J[3]*(J[4]*(f[ 2]*J[14]-J[10]*f[ 3])-f[2]*(J[8]*J[14]-J[10]*J[12])+J[6]*(J[8]*f[ 3]-f[ 2]*J[12])))/dd[0];
112: dd[3]=(J[0]*(J[5]*(f[ 2]*J[15]-J[11]*f[ 3])-f[1]*(J[9]*J[15]-J[11]*J[13])+J[7]*(J[9]*f[ 3]-f[ 2]*J[13]))-
113: J[1]*(J[4]*(f[ 2]*J[15]-J[11]*f[ 3])-f[1]*(J[8]*J[15]-J[11]*J[12])+J[7]*(J[8]*f[ 3]-f[ 2]*J[12]))+
114: f[0]*(J[4]*(J[ 9]*J[15]-J[11]*J[13])-J[5]*(J[8]*J[15]-J[11]*J[12])+J[7]*(J[8]*J[13]-J[ 9]*J[12]))-
115: J[3]*(J[4]*(J[ 9]*f[ 3]-f[ 2]*J[13])-J[5]*(J[8]*f[ 3]-f[ 2]*J[12])+f[1]*(J[8]*J[13]-J[ 9]*J[12])))/dd[0];
117: dd[4]=(J[0]*(J[5]*(J[10]*f[ 3]-f[ 2]*J[14])-J[6]*(J[9]*f[ 3]-f[ 2]*J[13])+f[1]*(J[9]*J[14]-J[10]*J[13]))-
118: J[1]*(J[4]*(J[10]*f[ 3]-f[ 2]*J[14])-J[6]*(J[8]*f[ 3]-f[ 2]*J[12])+f[1]*(J[8]*J[14]-J[10]*J[12]))+
119: J[2]*(J[4]*(J[ 9]*f[ 3]-f[ 2]*J[13])-J[5]*(J[8]*f[ 3]-f[ 2]*J[12])+f[1]*(J[8]*J[13]-J[ 9]*J[12]))-
120: f[0]*(J[4]*(J[ 9]*J[14]-J[10]*J[13])-J[5]*(J[8]*J[14]-J[10]*J[12])+J[6]*(J[8]*J[13]-J[ 9]*J[12])))/dd[0];
121: CHKMEMQ;
122: /* copy solution back into ADIC data structure */
123: ad_vustart[5*(gI+0)] += dd[1];
124: ad_vustart[5*(gI+1)] += dd[2];
125: ad_vustart[5*(gI+2)] += dd[3];
126: ad_vustart[5*(gI+3)] += dd[4];
127: CHKMEMQ;
128: res = f[0]*f[0];
129: res += f[1]*f[1];
130: res += f[2]*f[2];
131: res += f[3]*f[3];
132: res = sqrt(res);
133:
134: } while (--cnt > 0 && res > 1.e-14);
136: /* zero out this part of the seed matrix that was set initially */
137: ad_vustart[1+5*gI ] = 0.0;
138: ad_vustart[2+5*gI+5 ] = 0.0;
139: ad_vustart[3+5*gI+10] = 0.0;
140: ad_vustart[4+5*gI+15] = 0.0;
141: CHKMEMQ;
142: return(0);
143: }
147: PetscErrorCode NLFNewton_DAAD9(NLF A,DMDALocalInfo *info,MatStencil *stencil,void *ad_vu,PetscScalar *ad_vustart,int nI,int gI,PetscScalar *residual)
148: {
150: PetscInt cnt = A->newton_its;
151: PetscScalar ad_f[100], J[81],f[9], res;
152: PetscInt i,j,ngI[9];
154:
155: /* the order of the nodes */
156: /*
157: (6) (7) (8)
158: i-1,j+1 --- i,j+1 --- i+1,j+1
159: | | |
160: | | |
161: i-1,j --- i,j --- i+1,j
162: |(3) |(4) |(5)
163: | | |
164: i-1,j-1 --- i,j-1--- i+1,j-1
165: (0) (1) (2)
166: */
167:
168: /* the order of the derivative for the center nodes */
169: /*
170: (7) (8) (9)
171: i-1,j+1 --- i,j+1 --- i+1,j+1
172: | | |
173: | | |
174: i-1,j --- i,j --- i+1,j
175: |(4) |(5) |(6)
176: | | |
177: i-1,j-1 --- i,j-1--- i+1,j-1
178: (1) (2) (3)
179: */
180: if( (*stencil).i==0 || (*stencil).i==1||(*stencil).i==(*info).gxs+(*info).gxm-1 || (*stencil).i==(*info).gxs+(*info).gxm-2 || (*stencil).j==0 || (*stencil).j==1 ||(*stencil).j==(*info).gys+(*info).gym-1 || (*stencil).j==(*info).gys+(*info).gym -2) {
182: ad_vustart[1+10*gI] = 1.0;
183:
184: do {
185: /* compute the function and Jacobian */
186: (*A->da->adicmf_lfi)(info,stencil,ad_vu,ad_f,A->ctx);
187: J[0] = -ad_f[1];
188: f[0] = -ad_f[0] + residual[gI];
189: ad_vustart[10*gI] = ad_vustart[10*gI] - f[0]/J[0];
190: } while (--cnt > 0 && PetscAbsScalar(f[0]) > 1.e-14);
192: ad_vustart[1+10*gI] = 0.0;
193: return(0);
196: }
197:
198: ngI[0] = ((*stencil).i -1 - (*info).gxs)*(*info).dof + ((*stencil).j -1 - (*info).gys)*(*info).dof*(*info).gxm + ((*stencil).k - (*info).gzs)*(*info).dof*(*info).gxm*(*info).gym;
199: ngI[1] = ngI[0] + 1;
200: ngI[2] = ngI[1] + 1;
201: ngI[3] = gI - 1;
202: ngI[4] = gI ;
203: ngI[5] = gI + 1;
204: ngI[6] = ((*stencil).i -1 - (*info).gxs)*(*info).dof + ((*stencil).j +1 - (*info).gys)*(*info).dof*(*info).gxm + ((*stencil).k - (*info).gzs)*(*info).dof*(*info).gxm*(*info).gym;
205: ngI[7] = ngI[6] + 1;
206: ngI[8] = ngI[7] + 1;
207:
209: for(j=0 ; j<9; j++){
210: ad_vustart[ngI[j]*10+j+1] = 1.0;
211: }
212:
213: do{
214: /* compute the function and the Jacobian */
215:
216: (*A->da->adicmf_lfi)(info,stencil,ad_vu,ad_f,A->ctx);
217:
218: for(i=0; i<9; i++){
219: for(j=0; j<9; j++){
220: J[i*9+j] = -ad_f[i*10+j+1];
221: }
222: f[i]= -ad_f[i*10] + residual[ngI[i]];
223: }
225: PetscKernel_A_gets_inverse_A_9(J,0.0);
226:
227: res =0 ;
228: for(i=0; i<9; i++){
229: for(j=0;j<9;j++){
230: ad_vustart[10*ngI[i]]= ad_vustart[10*ngI[i]] - J[i*9 +j]*f[j];
231: }
232: res = res + f[i]*f[i];
233: }
234: res = sqrt(res);
235: } while(--cnt>0 && res>1.e-14);
237: for(j=0; j<9; j++){
238: ad_vustart[10*ngI[j]+j+1]=0.0;
239: }
241: return(0);
242: }
244: /*
245: Nonlinear relax on all the equations with an initial guess in x
246: */
247: EXTERN_C_BEGIN
250: PetscErrorCode NLFRelax_DAAD(NLF A,MatSORType flag,int its,Vec xx)
251: {
253: PetscInt j,gtdof,nI,gI;
254: PetscScalar *avu,*av,*ad_vustart,*residual;
255: Vec localxx;
256: DMDALocalInfo info;
257: MatStencil stencil;
258: void* *ad_vu;
261: if (its <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D positive",its);
263: DMGetLocalVector(A->da,&localxx);
264: /* get space for derivative object. */
265: DMDAGetAdicMFArray(A->da,PETSC_TRUE,(void **)&ad_vu,(void**)&ad_vustart,>dof);
266: VecGetArray(A->residual,&residual);
269: /* tell ADIC we will be computing one dimensional Jacobians */
270: PetscADResetIndep();
271: PetscADIncrementTotalGradSize(1);
272: PetscADSetIndepDone();
274: DMDAGetLocalInfo(A->da,&info);
275: while (its--) {
277: /* get initial solution properly ghosted */
278: DMGlobalToLocalBegin(A->da,xx,INSERT_VALUES,localxx);
279: DMGlobalToLocalEnd(A->da,xx,INSERT_VALUES,localxx);
281: /* copy input vector into derivative object */
282: VecGetArray(localxx,&avu);
283: for (j=0; j<gtdof; j++) {
284: ad_vustart[2*j] = avu[j];
285: ad_vustart[2*j+1] = 0.0;
286: }
287:
288: VecRestoreArray(localxx,&avu);
290: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){
291: nI = 0;
292: for (stencil.k = info.zs; stencil.k<info.zs+info.zm; stencil.k++) {
293: for (stencil.j = info.ys; stencil.j<info.ys+info.ym; stencil.j++) {
294: for (stencil.i = info.xs; stencil.i<info.xs+info.xm; stencil.i++) {
295: for (stencil.c = 0; stencil.c<info.dof; stencil.c++) {
296: gI = stencil.c + (stencil.i - info.gxs)*info.dof + (stencil.j - info.gys)*info.dof*info.gxm + (stencil.k - info.gzs)*info.dof*info.gxm*info.gym;
297: NLFNewton_DAAD(A,&info,&stencil,ad_vu,ad_vustart,nI,gI,residual[nI]);
298: nI++;
299: }
300: }
301: }
302: }
303: }
304: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){
305: nI = info.dof*info.xm*info.ym*info.zm - 1;
306: for (stencil.k = info.zs+info.zm-1; stencil.k>=info.zs; stencil.k--) {
307: for (stencil.j = info.ys+info.ym-1; stencil.j>=info.ys; stencil.j--) {
308: for (stencil.i = info.xs+info.xm-1; stencil.i>=info.xs; stencil.i--) {
309: for (stencil.c = info.dof-1; stencil.c>=0; stencil.c--) {
310: gI = stencil.c + (stencil.i - info.gxs)*info.dof + (stencil.j - info.gys)*info.dof*info.gxm + (stencil.k - info.gzs)*info.dof*info.gxm*info.gym;
311: NLFNewton_DAAD(A,&info,&stencil,ad_vu,ad_vustart,nI,gI,residual[nI]);
312: nI--;
313: }
314: }
315: }
316: }
317: }
319: /* copy solution back into ghosted vector from derivative object */
320: VecGetArray(localxx,&av);
321: for (j=0; j<gtdof; j++) {
322: av[j] = ad_vustart[2*j];
323: }
324: VecRestoreArray(localxx,&av);
325: /* stick relaxed solution back into global solution */
326: DMLocalToGlobalBegin(A->da,localxx,INSERT_VALUES,xx);
327: DMLocalToGlobalEnd(A->da,localxx,INSERT_VALUES,xx);
328: }
331: VecRestoreArray(A->residual,&residual);
332: DMRestoreLocalVector(A->da,&localxx);
333: DMDARestoreAdicMFArray(A->da,PETSC_TRUE,(void **)&ad_vu,(void**)&ad_vustart,>dof);
334: return(0);
335: }
336: EXTERN_C_END
338: EXTERN_C_BEGIN
341: PetscErrorCode NLFRelax_DAAD4(NLF A,MatSORType flag,int its,Vec xx)
342: {
344: PetscInt j,gtdof,nI,gI;
345: PetscScalar *avu,*av,*ad_vustart,*residual;
346: Vec localxx;
347: DMDALocalInfo info;
348: MatStencil stencil;
349: void* *ad_vu;
352: if (its <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D positive",its);
353:
354: DMGetLocalVector(A->da,&localxx);
355: /* get space for derivative object. */
356: DMDAGetAdicMFArray4(A->da,PETSC_TRUE,(void **)&ad_vu,(void**)&ad_vustart,>dof);
357: VecGetArray(A->residual,&residual);
360: /* tell ADIC we will be computing four dimensional Jacobians */
361: PetscADResetIndep();
362: PetscADIncrementTotalGradSize(4);
363: PetscADSetIndepDone();
365: DMDAGetLocalInfo(A->da,&info);
366: while (its--) {
368: /* get initial solution properly ghosted */
369: DMGlobalToLocalBegin(A->da,xx,INSERT_VALUES,localxx);
370: DMGlobalToLocalEnd(A->da,xx,INSERT_VALUES,localxx);
372: /* copy input vector into derivative object */
373: VecGetArray(localxx,&avu);
374: for (j=0; j<gtdof; j++) {
375: ad_vustart[5*j ] = avu[j];
376: ad_vustart[5*j+1] = 0.0;
377: ad_vustart[5*j+2] = 0.0;
378: ad_vustart[5*j+3] = 0.0;
379: ad_vustart[5*j+4] = 0.0;
380: }
381:
382: VecRestoreArray(localxx,&avu);
384: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){
385: nI = 0;
386: for (stencil.k = info.zs; stencil.k<info.zs+info.zm; stencil.k++) {
387: for (stencil.j = info.ys; stencil.j<info.ys+info.ym; stencil.j++) {
388: for (stencil.i = info.xs; stencil.i<info.xs+info.xm; stencil.i++) {
389: CHKMEMQ;
390: gI = (stencil.i - info.gxs)*info.dof + (stencil.j - info.gys)*info.dof*info.gxm + (stencil.k - info.gzs)*info.dof*info.gxm*info.gym;
392: NLFNewton_DAAD4(A,&info,&stencil,ad_vu,ad_vustart,nI,gI,residual+nI);
393: nI=nI+4;
394: CHKMEMQ;
395: }
396: }
397: }
398: }
399: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){
400: nI = info.dof*info.xm*info.ym*info.zm - 4;
401:
402: for (stencil.k = info.zs+info.zm-1; stencil.k>=info.zs; stencil.k--) {
403: for (stencil.j = info.ys+info.ym-1; stencil.j>=info.ys; stencil.j--) {
404: for (stencil.i = info.xs+info.xm-1; stencil.i>=info.xs; stencil.i--) {
405: gI = (stencil.i - info.gxs)*info.dof + (stencil.j - info.gys)*info.dof*info.gxm + (stencil.k - info.gzs)*info.dof*info.gxm*info.gym;
406: NLFNewton_DAAD4(A,&info,&stencil,ad_vu,ad_vustart,nI,gI,residual+nI);
407: nI=nI-4;
408: }
409: }
410: }
411: }
412:
413: /* copy solution back into ghosted vector from derivative object */
414: VecGetArray(localxx,&av);
415: for (j=0; j<gtdof; j++) {
416: av[j] = ad_vustart[5*j];
417: }
418: VecRestoreArray(localxx,&av);
419: /* stick relaxed solution back into global solution */
420: DMLocalToGlobalBegin(A->da,localxx,INSERT_VALUES,xx);
421: DMLocalToGlobalEnd(A->da,localxx,INSERT_VALUES,xx);
422: }
425: VecRestoreArray(A->residual,&residual);
426: DMRestoreLocalVector(A->da,&localxx);
427: DMDARestoreAdicMFArray(A->da,PETSC_TRUE,(void **)&ad_vu,(void**)&ad_vustart,>dof);
428: return(0);
429: }
430: EXTERN_C_END
434: PetscErrorCode NLFRelax_DAAD9(NLF A,MatSORType flag,int its,Vec xx)
435: {
437: PetscInt j,gtdof,nI,gI;
438: PetscScalar *avu,*av,*ad_vustart,*residual;
439: Vec localxx;
440: DMDALocalInfo info;
441: MatStencil stencil;
442: void* *ad_vu;
445: if (its <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D positive",its);
446:
447: DMGetLocalVector(A->da,&localxx);
448: /* get space for derivative object. */
449: DMDAGetAdicMFArray9(A->da,PETSC_TRUE,(void **)&ad_vu,(void**)&ad_vustart,>dof);
450: VecGetArray(A->residual,&residual);
453: /* tell ADIC we will be computing nine dimensional Jacobians */
454: PetscADResetIndep();
455: PetscADIncrementTotalGradSize(9);
456: PetscADSetIndepDone();
458: DMDAGetLocalInfo(A->da,&info);
459: while (its--) {
461: /* get initial solution properly ghosted */
462: DMGlobalToLocalBegin(A->da,xx,INSERT_VALUES,localxx);
463: DMGlobalToLocalEnd(A->da,xx,INSERT_VALUES,localxx);
465: /* copy input vector into derivative object */
466: VecGetArray(localxx,&avu);
467: for (j=0; j<gtdof; j++) {
468: ad_vustart[10*j ] = avu[j];
469: ad_vustart[10*j+1] = 0.0;
470: ad_vustart[10*j+2] = 0.0;
471: ad_vustart[10*j+3] = 0.0;
472: ad_vustart[10*j+4] = 0.0;
473: ad_vustart[10*j+5] = 0.0;
474: ad_vustart[10*j+6] = 0.0;
475: ad_vustart[10*j+7] = 0.0;
476: ad_vustart[10*j+8] = 0.0;
477: ad_vustart[10*j+9] = 0.0;
478: }
479:
480: VecRestoreArray(localxx,&avu);
482: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){
483: nI = 0;
484: for (stencil.k = info.zs; stencil.k<info.zs+info.zm; stencil.k++) {
485: for (stencil.j = info.ys; stencil.j<info.ys+info.ym; stencil.j++) {
486: for (stencil.i = info.xs; stencil.i<info.xs+info.xm; stencil.i++) {
487: CHKMEMQ;
488: gI = (stencil.i - info.gxs)*info.dof + (stencil.j - info.gys)*info.dof*info.gxm + (stencil.k - info.gzs)*info.dof*info.gxm*info.gym;
489: NLFNewton_DAAD9(A,&info,&stencil,ad_vu,ad_vustart,nI,gI,residual);
490: nI=nI+1;
491: CHKMEMQ;
492: }
493: }
494: }
495: }
496: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){
497: nI = info.dof*info.xm*info.ym*info.zm - 4;
498:
499: for (stencil.k = info.zs+info.zm-1; stencil.k>=info.zs; stencil.k--) {
500: for (stencil.j = info.ys+info.ym-1; stencil.j>=info.ys; stencil.j--) {
501: for (stencil.i = info.xs+info.xm-1; stencil.i>=info.xs; stencil.i--) {
502: gI = (stencil.i - info.gxs)*info.dof + (stencil.j - info.gys)*info.dof*info.gxm + (stencil.k - info.gzs)*info.dof*info.gxm*info.gym;
503: NLFNewton_DAAD9(A,&info,&stencil,ad_vu,ad_vustart,nI,gI,residual);
504: nI=nI-1;
505: }
506: }
507: }
508: }
509:
510: /* copy solution back into ghosted vector from derivative object */
511: VecGetArray(localxx,&av);
512: for (j=0; j<gtdof; j++) {
513: av[j] = ad_vustart[10*j];
514: }
515: VecRestoreArray(localxx,&av);
516: /* stick relaxed solution back into global solution */
517: DMLocalToGlobalBegin(A->da,localxx,INSERT_VALUES,xx);
518: DMLocalToGlobalEnd(A->da,localxx,INSERT_VALUES,xx);
519: }
522: VecRestoreArray(A->residual,&residual);
523: DMRestoreLocalVector(A->da,&localxx);
524: DMDARestoreAdicMFArray(A->da,PETSC_TRUE,(void **)&ad_vu,(void**)&ad_vustart,>dof);
525: return(0);
526: }
528: /*
529: Point-block nonlinear relax on all the equations with an initial guess in xx using
530: */
531: EXTERN_C_BEGIN
534: PetscErrorCode NLFRelax_DAADb(NLF A,MatSORType flag,int its,Vec xx)
535: {
537: PetscInt i,j,gtdof,nI,gI, bs = A->da->w, bs1 = bs + 1;
538: PetscScalar *avu,*av,*ad_vustart,*residual;
539: Vec localxx;
540: DMDALocalInfo info;
541: MatStencil stencil;
542: void* *ad_vu;
543: PetscErrorCode (*NLFNewton_DAADb)(NLF,DMDALocalInfo*,MatStencil*,void*,PetscScalar*,int,int,PetscScalar*);
546: if (its <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D positive",its);
547: if (bs == 4) {
548: NLFNewton_DAADb = NLFNewton_DAAD4;
549: } else {
550: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Point block nonlinear relaxation currently not for this block size",bs);
551: }
553: DMGetLocalVector(A->da,&localxx);
554: /* get space for derivative object. */
555: DMDAGetAdicMFArrayb(A->da,PETSC_TRUE,(void **)&ad_vu,(void**)&ad_vustart,>dof);
556: VecGetArray(A->residual,&residual);
559: /* tell ADIC we will be computing bs dimensional Jacobians */
560: PetscADResetIndep();
561: PetscADIncrementTotalGradSize(bs);
562: PetscADSetIndepDone();
564: DMDAGetLocalInfo(A->da,&info);
565: while (its--) {
567: /* get initial solution properly ghosted */
568: DMGlobalToLocalBegin(A->da,xx,INSERT_VALUES,localxx);
569: DMGlobalToLocalEnd(A->da,xx,INSERT_VALUES,localxx);
571: /* copy input vector into derivative object */
572: VecGetArray(localxx,&avu);
573: for (j=0; j<gtdof; j++) {
574: ad_vustart[bs1*j] = avu[j];
575: for (i=0; i<bs; i++) {
576: ad_vustart[bs1*j+1+i] = 0.0;
577: }
578: }
579: VecRestoreArray(localxx,&avu);
581: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){
582: nI = 0;
583: for (stencil.k = info.zs; stencil.k<info.zs+info.zm; stencil.k++) {
584: for (stencil.j = info.ys; stencil.j<info.ys+info.ym; stencil.j++) {
585: for (stencil.i = info.xs; stencil.i<info.xs+info.xm; stencil.i++) {
586: gI = (stencil.i - info.gxs)*info.dof + (stencil.j - info.gys)*info.dof*info.gxm + (stencil.k - info.gzs)*info.dof*info.gxm*info.gym;
587: (*NLFNewton_DAADb)(A,&info,&stencil,ad_vu,ad_vustart,nI,gI,residual+nI);
588: nI += bs;
589: }
590: }
591: }
592: }
593: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){
594: nI = info.dof*info.xm*info.ym*info.zm - bs;
595: for (stencil.k = info.zs+info.zm-1; stencil.k>=info.zs; stencil.k--) {
596: for (stencil.j = info.ys+info.ym-1; stencil.j>=info.ys; stencil.j--) {
597: for (stencil.i = info.xs+info.xm-1; stencil.i>=info.xs; stencil.i--) {
598: gI = (stencil.i - info.gxs)*info.dof + (stencil.j - info.gys)*info.dof*info.gxm + (stencil.k - info.gzs)*info.dof*info.gxm*info.gym;
599: (*NLFNewton_DAADb)(A,&info,&stencil,ad_vu,ad_vustart,nI,gI,residual+nI);
600: nI -= bs;
601: }
602: }
603: }
604: }
606: /* copy solution back into ghosted vector from derivative object */
607: VecGetArray(localxx,&av);
608: for (j=0; j<gtdof; j++) {
609: av[j] = ad_vustart[bs1*j];
610: }
611: VecRestoreArray(localxx,&av);
612: /* stick relaxed solution back into global solution */
613: DMLocalToGlobalBegin(A->da,localxx,INSERT_VALUES,xx);
614: DMLocalToGlobalEnd(A->da,localxx,INSERT_VALUES,xx);
615: }
617: VecRestoreArray(A->residual,&residual);
618: DMRestoreLocalVector(A->da,&localxx);
619: DMDARestoreAdicMFArray(A->da,PETSC_TRUE,(void **)&ad_vu,(void**)&ad_vustart,>dof);
620: return(0);
621: }
622: EXTERN_C_END
627: PetscErrorCode NLFDestroy_DAAD(NLF A)
628: {
632: DMDestroy(A->da);
633: PetscFree(A);
634: return(0);
635: }
637: EXTERN_C_BEGIN
640: PetscErrorCode NLFDAADSetDA_DAAD(NLF A,DM da)
641: {
645: PetscObjectReference((PetscObject)da);
646: if (A->da) {DMDestroy(A->da);}
647: A->da = da;
648: return(0);
649: }
650: EXTERN_C_END
652: EXTERN_C_BEGIN
655: PetscErrorCode NLFDAADSetNewtonIterations_DAAD(NLF A,int its)
656: {
658: A->newton_its = its;
659: return(0);
660: }
661: EXTERN_C_END
663: EXTERN_C_BEGIN
666: PetscErrorCode NLFDAADSetResidual_DAAD(NLF A,Vec residual)
667: {
669: A->residual = residual;
670: return(0);
671: }
672: EXTERN_C_END
675: EXTERN_C_BEGIN
678: PetscErrorCode NLFDAADSetCtx_DAAD(NLF A,void *ctx)
679: {
681: A->ctx = ctx;
682: return(0);
683: }
684: EXTERN_C_END
686: EXTERN_C_BEGIN
689: PetscErrorCode NLFCreate_DAAD(NLF *A)
690: {
694: PetscNew(struct NLF_DAAD,A);
695: (*A)->da = 0;
696: (*A)->ctx = 0;
697: (*A)->newton_its = 2;
698: return(0);
699: }
700: EXTERN_C_END