Actual source code: snesnoise.c
petsc-3.14.6 2021-03-30
2: #include <petsc/private/snesimpl.h>
4: PETSC_INTERN PetscErrorCode SNESDiffParameterCreate_More(SNES,Vec,void**);
5: PETSC_INTERN PetscErrorCode SNESDiffParameterCompute_More(SNES,void*,Vec,Vec,PetscReal*,PetscReal*);
6: PETSC_INTERN PetscErrorCode SNESDiffParameterDestroy_More(void*);
8: /* Data used by Jorge's diff parameter computation method */
9: typedef struct {
10: Vec *workv; /* work vectors */
11: FILE *fp; /* output file */
12: PetscInt function_count; /* count of function evaluations for diff param estimation */
13: double fnoise_min; /* minimim allowable noise */
14: double hopt_min; /* minimum allowable hopt */
15: double h_first_try; /* first try for h used in diff parameter estimate */
16: PetscInt fnoise_resets; /* number of times we've reset the noise estimate */
17: PetscInt hopt_resets; /* number of times we've reset the hopt estimate */
18: } DIFFPAR_MORE;
21: PETSC_INTERN PetscErrorCode SNESDefaultMatrixFreeSetParameters2(Mat,double,double,double);
22: PETSC_INTERN PetscErrorCode SNESUnSetMatrixFreeParameter(SNES snes);
23: PETSC_INTERN PetscErrorCode SNESNoise_dnest_(PetscInt*,PetscScalar*,PetscScalar*,PetscScalar*,PetscScalar*,PetscScalar*,PetscInt*,PetscScalar*);
25: static PetscErrorCode JacMatMultCompare(SNES,Vec,Vec,double);
27: PetscErrorCode SNESDiffParameterCreate_More(SNES snes,Vec x,void **outneP)
28: {
29: DIFFPAR_MORE *neP;
30: Vec w;
31: PetscRandom rctx; /* random number generator context */
33: PetscBool flg;
34: char noise_file[PETSC_MAX_PATH_LEN];
37: PetscNewLog(snes,&neP);
39: neP->function_count = 0;
40: neP->fnoise_min = 1.0e-20;
41: neP->hopt_min = 1.0e-8;
42: neP->h_first_try = 1.0e-3;
43: neP->fnoise_resets = 0;
44: neP->hopt_resets = 0;
46: /* Create work vectors */
47: VecDuplicateVecs(x,3,&neP->workv);
48: w = neP->workv[0];
50: /* Set components of vector w to random numbers */
51: PetscRandomCreate(PetscObjectComm((PetscObject)snes),&rctx);
52: PetscRandomSetFromOptions(rctx);
53: VecSetRandom(w,rctx);
54: PetscRandomDestroy(&rctx);
56: /* Open output file */
57: PetscOptionsGetString(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_mf_noise_file",noise_file,sizeof(noise_file),&flg);
58: if (flg) neP->fp = fopen(noise_file,"w");
59: else neP->fp = fopen("noise.out","w");
60: if (!neP->fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_OPEN,"Cannot open file");
61: PetscInfo(snes,"Creating Jorge's differencing parameter context\n");
63: *outneP = neP;
64: return(0);
65: }
67: PetscErrorCode SNESDiffParameterDestroy_More(void *nePv)
68: {
69: DIFFPAR_MORE *neP = (DIFFPAR_MORE*)nePv;
71: int err;
74: /* Destroy work vectors and close output file */
75: VecDestroyVecs(3,&neP->workv);
76: err = fclose(neP->fp);
77: if (err) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SYS,"fclose() failed on file");
78: PetscFree(neP);
79: return(0);
80: }
82: PetscErrorCode SNESDiffParameterCompute_More(SNES snes,void *nePv,Vec x,Vec p,double *fnoise,double *hopt)
83: {
84: DIFFPAR_MORE *neP = (DIFFPAR_MORE*)nePv;
85: Vec w, xp, fvec; /* work vectors to use in computing h */
86: double zero = 0.0, hl, hu, h, fnoise_s, fder2_s;
87: PetscScalar alpha;
88: PetscScalar fval[7], tab[7][7], eps[7], f = -1;
89: double rerrf = -1., fder2;
91: PetscInt iter, k, i, j, info;
92: PetscInt nf = 7; /* number of function evaluations */
93: PetscInt fcount;
94: MPI_Comm comm;
95: FILE *fp;
96: PetscBool noise_test = PETSC_FALSE;
99: PetscObjectGetComm((PetscObject)snes,&comm);
100: /* Call to SNESSetUp() just to set data structures in SNES context */
101: if (!snes->setupcalled) {SNESSetUp(snes);}
103: w = neP->workv[0];
104: xp = neP->workv[1];
105: fvec = neP->workv[2];
106: fp = neP->fp;
108: /* Initialize parameters */
109: hl = zero;
110: hu = zero;
111: h = neP->h_first_try;
112: fnoise_s = zero;
113: fder2_s = zero;
114: fcount = neP->function_count;
116: /* We have 5 tries to attempt to compute a good hopt value */
117: SNESGetIterationNumber(snes,&i);
118: PetscFPrintf(comm,fp,"\n ------- SNES iteration %D ---------\n",i);
119: for (iter=0; iter<5; iter++) {
120: neP->h_first_try = h;
122: /* Compute the nf function values needed to estimate the noise from
123: the difference table */
124: for (k=0; k<nf; k++) {
125: alpha = h * (k+1 - (nf+1)/2);
126: VecWAXPY(xp,alpha,p,x);
127: SNESComputeFunction(snes,xp,fvec);
128: neP->function_count++;
129: VecDot(fvec,w,&fval[k]);
130: }
131: f = fval[(nf+1)/2 - 1];
133: /* Construct the difference table */
134: for (i=0; i<nf; i++) tab[i][0] = fval[i];
136: for (j=0; j<nf-1; j++) {
137: for (i=0; i<nf-j-1; i++) {
138: tab[i][j+1] = tab[i+1][j] - tab[i][j];
139: }
140: }
142: /* Print the difference table */
143: PetscFPrintf(comm,fp,"Difference Table: iter = %D\n",iter);
144: for (i=0; i<nf; i++) {
145: for (j=0; j<nf-i; j++) {
146: PetscFPrintf(comm,fp," %10.2e ",tab[i][j]);
147: }
148: PetscFPrintf(comm,fp,"\n");
149: }
151: /* Call the noise estimator */
152: SNESNoise_dnest_(&nf,fval,&h,fnoise,&fder2,hopt,&info,eps);
154: /* Output statements */
155: rerrf = *fnoise/PetscAbsScalar(f);
156: if (info == 1) {PetscFPrintf(comm,fp,"%s\n","Noise detected");}
157: if (info == 2) {PetscFPrintf(comm,fp,"%s\n","Noise not detected; h is too small");}
158: if (info == 3) {PetscFPrintf(comm,fp,"%s\n","Noise not detected; h is too large");}
159: if (info == 4) {PetscFPrintf(comm,fp,"%s\n","Noise detected, but unreliable hopt");}
160: PetscFPrintf(comm,fp,"Approximate epsfcn %g %g %g %g %g %g\n",(double)eps[0],(double)eps[1],(double)eps[2],(double)eps[3],(double)eps[4],(double)eps[5]);
161: PetscFPrintf(comm,fp,"h = %g, fnoise = %g, fder2 = %g, rerrf = %g, hopt = %g\n\n",(double)h, (double)*fnoise, (double)fder2, (double)rerrf, (double)*hopt);
163: /* Save fnoise and fder2. */
164: if (*fnoise) fnoise_s = *fnoise;
165: if (fder2) fder2_s = fder2;
167: /* Check for noise detection. */
168: if (fnoise_s && fder2_s) {
169: *fnoise = fnoise_s;
170: fder2 = fder2_s;
171: *hopt = 1.68*sqrt(*fnoise/PetscAbsScalar(fder2));
172: goto theend;
173: } else {
175: /* Update hl and hu, and determine new h */
176: if (info == 2 || info == 4) {
177: hl = h;
178: if (hu == zero) h = 100*h;
179: else h = PetscMin(100*h,0.1*hu);
180: } else if (info == 3) {
181: hu = h;
182: h = PetscMax(1.0e-3,sqrt(hl/hu))*hu;
183: }
184: }
185: }
186: theend:
188: if (*fnoise < neP->fnoise_min) {
189: PetscFPrintf(comm,fp,"Resetting fnoise: fnoise1 = %g, fnoise_min = %g\n",(double)*fnoise,(double)neP->fnoise_min);
190: *fnoise = neP->fnoise_min;
191: neP->fnoise_resets++;
192: }
193: if (*hopt < neP->hopt_min) {
194: PetscFPrintf(comm,fp,"Resetting hopt: hopt1 = %g, hopt_min = %g\n",(double)*hopt,(double)neP->hopt_min);
195: *hopt = neP->hopt_min;
196: neP->hopt_resets++;
197: }
199: PetscFPrintf(comm,fp,"Errors in derivative:\n");
200: PetscFPrintf(comm,fp,"f = %g, fnoise = %g, fder2 = %g, hopt = %g\n",(double)f,(double)*fnoise,(double)fder2,(double)*hopt);
202: /* For now, compute h **each** MV Mult!! */
203: /*
204: PetscOptionsHasName(NULL,"-matrix_free_jorge_each_mvp",&flg);
205: if (!flg) {
206: Mat mat;
207: SNESGetJacobian(snes,&mat,NULL,NULL);
208: SNESDefaultMatrixFreeSetParameters2(mat,PETSC_DEFAULT,PETSC_DEFAULT,*hopt);
209: }
210: */
211: fcount = neP->function_count - fcount;
212: PetscInfo5(snes,"fct_now = %D, fct_cum = %D, rerrf=%g, sqrt(noise)=%g, h_more=%g\n",fcount,neP->function_count,(double)rerrf,(double)PetscSqrtReal(*fnoise),(double)*hopt);
214: PetscOptionsGetBool(NULL,NULL,"-noise_test",&noise_test,NULL);
215: if (noise_test) {
216: JacMatMultCompare(snes,x,p,*hopt);
217: }
218: return(0);
219: }
221: PetscErrorCode JacMatMultCompare(SNES snes,Vec x,Vec p,double hopt)
222: {
223: Vec yy1, yy2; /* work vectors */
224: PetscViewer view2; /* viewer */
225: Mat J; /* analytic Jacobian (set as preconditioner matrix) */
226: Mat Jmf; /* matrix-free Jacobian (set as true system matrix) */
227: double h; /* differencing parameter */
228: Vec f;
229: PetscScalar alpha;
230: PetscReal yy1n,yy2n,enorm;
232: PetscInt i;
233: PetscBool printv = PETSC_FALSE;
234: char filename[32];
235: MPI_Comm comm;
238: PetscObjectGetComm((PetscObject)snes,&comm);
239: /* Compute function and analytic Jacobian at x */
240: SNESGetJacobian(snes,&Jmf,&J,NULL,NULL);
241: SNESComputeJacobian(snes,x,Jmf,J);
242: SNESGetFunction(snes,&f,NULL,NULL);
243: SNESComputeFunction(snes,x,f);
245: /* Duplicate work vectors */
246: VecDuplicate(x,&yy2);
247: VecDuplicate(x,&yy1);
249: /* Compute true matrix-vector product */
250: MatMult(J,p,yy1);
251: VecNorm(yy1,NORM_2,&yy1n);
253: /* View product vector if desired */
254: PetscOptionsGetBool(NULL,NULL,"-print_vecs",&printv,NULL);
255: if (printv) {
256: PetscViewerASCIIOpen(comm,"y1.out",&view2);
257: PetscViewerPushFormat(view2,PETSC_VIEWER_ASCII_COMMON);
258: VecView(yy1,view2);
259: PetscViewerPopFormat(view2);
260: PetscViewerDestroy(&view2);
261: }
263: /* Test Jacobian-vector product computation */
264: alpha = -1.0;
265: h = 0.01 * hopt;
266: for (i=0; i<5; i++) {
267: /* Set differencing parameter for matrix-free multiplication */
268: SNESDefaultMatrixFreeSetParameters2(Jmf,PETSC_DEFAULT,PETSC_DEFAULT,h);
270: /* Compute matrix-vector product via differencing approximation */
271: MatMult(Jmf,p,yy2);
272: VecNorm(yy2,NORM_2,&yy2n);
274: /* View product vector if desired */
275: if (printv) {
276: sprintf(filename,"y2.%d.out",(int)i);
277: PetscViewerASCIIOpen(comm,filename,&view2);
278: PetscViewerPushFormat(view2,PETSC_VIEWER_ASCII_COMMON);
279: VecView(yy2,view2);
280: PetscViewerPopFormat(view2);
281: PetscViewerDestroy(&view2);
282: }
284: /* Compute relative error */
285: VecAXPY(yy2,alpha,yy1);
286: VecNorm(yy2,NORM_2,&enorm);
287: enorm = enorm/yy1n;
288: PetscFPrintf(comm,stdout,"h = %g: relative error = %g\n",(double)h,(double)enorm);
289: h *= 10.0;
290: }
291: return(0);
292: }
294: static PetscInt lin_its_total = 0;
296: PetscErrorCode SNESNoiseMonitor(SNES snes,PetscInt its,double fnorm,void *dummy)
297: {
299: PetscInt lin_its;
302: SNESGetLinearSolveIterations(snes,&lin_its);
303: lin_its_total += lin_its;
304: PetscPrintf(PetscObjectComm((PetscObject)snes), "iter = %D, SNES Function norm = %g, lin_its = %D, total_lin_its = %D\n",its,(double)fnorm,lin_its,lin_its_total);
306: SNESUnSetMatrixFreeParameter(snes);
307: return(0);
308: }