← Back to branch summary

~ubuntu-branches/ubuntu/warty/petsc/warty

« back to all changes in this revision

Viewing changes to src/snes/utils/damgsnes.c

  • Committer: Bazaar Package Importer
  • Author(s): Adam C. Powell, IV
  • Date: 2004-06-07 13:41:43 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20040607134143-92p586zrauvie0le
Tags: 2.2.0-2
http://bugs.debian.org/249617
* Upstream patch level 2.
* New PETSC_BOPT_EXTRA option for different BOPT and lib names, with _c++
  symlinks only for plain and single (closes: #249617).
* New DEBIAN_DIST=contrib option to link with hypre, parmetis (closes:
  #249619).
* Combined petsc-c and petsc-fortran substvars into petsc-compilers.
* Extra quote in -dev prerm eliminates "too many arguments" problem.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/snes/utils/damgsnes.c
5
5
 
6
6
 
7
7
/*
 
8
      period of -1 indicates update only on zeroth iteration of SNES
 
9
*/
 
10
#define ShouldUpdate(l,it) (((dmmg[l-1]->updatejacobianperiod == -1) && (it == 0)) || \
 
11
                            ((dmmg[l-1]->updatejacobianperiod >   0) && !(it % dmmg[l-1]->updatejacobianperiod)))
 
12
/*
8
13
   Evaluates the Jacobian on all of the grids. It is used by DMMG to provide the 
9
14
   ComputeJacobian() function that SNESSetJacobian() requires.
10
15
*/
12
17
#define __FUNCT__ "DMMGComputeJacobian_Multigrid"
13
18
int DMMGComputeJacobian_Multigrid(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
14
19
{
15
 
  DMMG       *dmmg = (DMMG*)ptr;
16
 
  int        ierr,i,nlevels = dmmg[0]->nlevels;
17
 
  SLES       sles,lsles;
18
 
  PC         pc;
19
 
  PetscTruth ismg;
20
 
  Vec        W;
 
20
  DMMG         *dmmg = (DMMG*)ptr;
 
21
  int          ierr,i,nlevels = dmmg[0]->nlevels,it;
 
22
  KSP          ksp,lksp;
 
23
  PC           pc;
 
24
  PetscTruth   ismg;
 
25
  Vec          W;
 
26
  MatStructure flg;
21
27
 
22
28
  PetscFunctionBegin;
23
29
  if (!dmmg) SETERRQ(1,"Passing null as user context which should contain DMMG");
 
30
  ierr = SNESGetIterationNumber(snes,&it);CHKERRQ(ierr);
24
31
 
25
32
  /* compute Jacobian on finest grid */
26
 
  ierr = (*DMMGGetFine(dmmg)->computejacobian)(snes,X,J,B,flag,DMMGGetFine(dmmg));CHKERRQ(ierr);
 
33
  if (dmmg[nlevels-1]->updatejacobian && ShouldUpdate(nlevels,it)) {
 
34
    ierr = (*DMMGGetFine(dmmg)->computejacobian)(snes,X,J,B,flag,DMMGGetFine(dmmg));CHKERRQ(ierr);
 
35
  } else {
 
36
    PetscLogInfo(0,"DMMGComputeJacobian_Multigrid:Skipping Jacobian, SNES iteration %d frequence %d level %d\n",it,dmmg[nlevels-1]->updatejacobianperiod,nlevels-1);
 
37
    *flag = SAME_PRECONDITIONER;
 
38
  }
27
39
  ierr = MatSNESMFSetBase(DMMGGetFine(dmmg)->J,X);CHKERRQ(ierr);
28
40
 
29
41
  /* create coarser grid Jacobians for preconditioner if multigrid is the preconditioner */
30
 
  ierr = SNESGetSLES(snes,&sles);CHKERRQ(ierr);
31
 
  ierr = SLESGetPC(sles,&pc);CHKERRQ(ierr);
 
42
  ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr);
 
43
  ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
32
44
  ierr = PetscTypeCompare((PetscObject)pc,PCMG,&ismg);CHKERRQ(ierr);
33
45
  if (ismg) {
34
46
 
35
 
    ierr = MGGetSmoother(pc,nlevels-1,&lsles);CHKERRQ(ierr);
36
 
    ierr = SLESSetOperators(lsles,DMMGGetFine(dmmg)->J,DMMGGetFine(dmmg)->B,*flag);CHKERRQ(ierr);
37
 
 
38
 
    for (i=nlevels-1; i>0; i--) {
39
 
 
40
 
      if (!dmmg[i-1]->w) {
41
 
        ierr = VecDuplicate(dmmg[i-1]->x,&dmmg[i-1]->w);CHKERRQ(ierr);
 
47
    ierr = MGGetSmoother(pc,nlevels-1,&lksp);CHKERRQ(ierr);
 
48
    ierr = KSPSetOperators(lksp,DMMGGetFine(dmmg)->J,DMMGGetFine(dmmg)->B,*flag);CHKERRQ(ierr);
 
49
 
 
50
    if (dmmg[0]->galerkin) {
 
51
      for (i=nlevels-2; i>-1; i--) {
 
52
        PetscTruth JeqB = (PetscTruth)( dmmg[i]->B == dmmg[i]->J);
 
53
        ierr = MatDestroy(dmmg[i]->B);CHKERRQ(ierr);
 
54
        ierr = MatSeqAIJPtAP(dmmg[i+1]->B,dmmg[i+1]->R,&dmmg[i]->B);CHKERRQ(ierr);
 
55
        if (JeqB) dmmg[i]->J = dmmg[i]->B;
 
56
        ierr = MGGetSmoother(pc,i,&lksp);CHKERRQ(ierr);
 
57
        ierr = KSPSetOperators(lksp,dmmg[i]->J,dmmg[i]->B,*flag);CHKERRQ(ierr);
 
58
      }   
 
59
    } else {
 
60
      for (i=nlevels-1; i>0; i--) {
 
61
        if (!dmmg[i-1]->w) {
 
62
          ierr = VecDuplicate(dmmg[i-1]->x,&dmmg[i-1]->w);CHKERRQ(ierr);
 
63
        }
 
64
        W    = dmmg[i-1]->w;
 
65
        /* restrict X to coarser grid */
 
66
        ierr = MatRestrict(dmmg[i]->R,X,W);CHKERRQ(ierr);
 
67
        X    = W;      
 
68
        /* scale to "natural" scaling for that grid */
 
69
        ierr = VecPointwiseMult(dmmg[i]->Rscale,X,X);CHKERRQ(ierr);
 
70
        /* tell the base vector for matrix free multiplies */
 
71
        ierr = MatSNESMFSetBase(dmmg[i-1]->J,X);CHKERRQ(ierr);
 
72
        /* compute Jacobian on coarse grid */
 
73
        if (dmmg[i-1]->updatejacobian && ShouldUpdate(i,it)) {
 
74
          ierr = (*dmmg[i-1]->computejacobian)(snes,X,&dmmg[i-1]->J,&dmmg[i-1]->B,&flg,dmmg[i-1]);CHKERRQ(ierr);
 
75
        } else {
 
76
          PetscLogInfo(0,"DMMGComputeJacobian_Multigrid:Skipping Jacobian, SNES iteration %d frequence %d level %d\n",it,dmmg[i-1]->updatejacobianperiod,i-1);
 
77
          flg = SAME_PRECONDITIONER;
 
78
        }
 
79
        ierr = MGGetSmoother(pc,i-1,&lksp);CHKERRQ(ierr);
 
80
        ierr = KSPSetOperators(lksp,dmmg[i-1]->J,dmmg[i-1]->B,flg);CHKERRQ(ierr);
42
81
      }
43
 
 
44
 
      W    = dmmg[i-1]->w;
45
 
      /* restrict X to coarser grid */
46
 
      ierr = MatRestrict(dmmg[i]->R,X,W);CHKERRQ(ierr);
47
 
      X    = W;      
48
 
 
49
 
      /* scale to "natural" scaling for that grid */
50
 
      ierr = VecPointwiseMult(dmmg[i]->Rscale,X,X);CHKERRQ(ierr);
51
 
 
52
 
      /* tell the base vector for matrix free multiplies */
53
 
      ierr = MatSNESMFSetBase(dmmg[i-1]->J,X);CHKERRQ(ierr);
54
 
 
55
 
      /* compute Jacobian on coarse grid */
56
 
      ierr = (*dmmg[i-1]->computejacobian)(snes,X,&dmmg[i-1]->J,&dmmg[i-1]->B,flag,dmmg[i-1]);CHKERRQ(ierr);
57
 
 
58
 
      ierr = MGGetSmoother(pc,i-1,&lsles);CHKERRQ(ierr);
59
 
      ierr = SLESSetOperators(lsles,dmmg[i-1]->J,dmmg[i-1]->B,*flag);CHKERRQ(ierr);
60
82
    }
61
83
  }
62
84
  PetscFunctionReturn(0);
94
116
  */
95
117
  ierr = DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX);CHKERRQ(ierr);
96
118
  ierr = DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX);CHKERRQ(ierr);
97
 
  ierr = DAFormFunction1(da,localX,F,dmmg->user);
 
119
  ierr = DAFormFunction1(da,localX,F,dmmg->user);CHKERRQ(ierr);
98
120
  ierr = DARestoreLocalVector(da,&localX);CHKERRQ(ierr);
99
121
  PetscFunctionReturn(0); 
100
122
} 
135
157
  */
136
158
  ierr = DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX);CHKERRQ(ierr);
137
159
  ierr = DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX);CHKERRQ(ierr);
138
 
  ierr = DAFormFunction1(da,localX,F,ptr);
 
160
  ierr = DAFormFunction1(da,localX,F,ptr);CHKERRQ(ierr);
139
161
  ierr = DARestoreLocalVector(da,&localX);CHKERRQ(ierr);
140
162
  PetscFunctionReturn(0); 
141
163
} 
330
352
#define __FUNCT__ "DMMGSolveSNES"
331
353
int DMMGSolveSNES(DMMG *dmmg,int level)
332
354
{
333
 
  int  ierr,nlevels = dmmg[0]->nlevels,its;
 
355
  int  ierr,nlevels = dmmg[0]->nlevels;
334
356
 
335
357
  PetscFunctionBegin;
336
358
  dmmg[0]->nlevels = level+1;
337
 
  ierr             = SNESSolve(dmmg[level]->snes,dmmg[level]->x,&its);CHKERRQ(ierr);
 
359
  ierr             = SNESSolve(dmmg[level]->snes,dmmg[level]->x);CHKERRQ(ierr);
338
360
  dmmg[0]->nlevels = nlevels;
339
361
  PetscFunctionReturn(0);
340
362
}
341
363
 
 
364
EXTERN_C_BEGIN
 
365
extern int NLFCreate_DAAD(NLF*);
 
366
extern int NLFRelax_DAAD(NLF,MatSORType,int,Vec);
 
367
extern int NLFDAADSetDA_DAAD(NLF,DA);
 
368
extern int NLFDAADSetCtx_DAAD(NLF,void*);
 
369
extern int NLFDAADSetResidual_DAAD(NLF,Vec);
 
370
extern int NLFDAADSetNewtonIterations_DAAD(NLF,int);
 
371
EXTERN_C_END
 
372
 
 
373
#if defined(PETSC_HAVE_ADIC) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE)
 
374
#include "src/ksp/pc/impls/mg/mgimpl.h"                    /*I "petscmg.h" I*/
 
375
/*
 
376
          This is pre-beta FAS code. It's design should not be taken seriously!
 
377
*/
 
378
#undef __FUNCT__  
 
379
#define __FUNCT__ "DMMGSolveFAS"
 
380
int DMMGSolveFAS(DMMG *dmmg,int level)
 
381
{
 
382
  int         ierr,i,j,k;
 
383
  PetscReal   norm;
 
384
  PetscScalar zero = 0.0,mone = -1.0,one = 1.0;
 
385
  MG          *mg;
 
386
  PC          pc;
 
387
  KSP         ksp;
 
388
 
 
389
  PetscFunctionBegin;
 
390
  ierr = VecSet(&zero,dmmg[level]->r);CHKERRQ(ierr);
 
391
  for (j=1; j<=level; j++) {
 
392
    if (!dmmg[j]->inject) {
 
393
      ierr = DMGetInjection(dmmg[j-1]->dm,dmmg[j]->dm,&dmmg[j]->inject);CHKERRQ(ierr);
 
394
    }
 
395
  }
 
396
 
 
397
  ierr = SNESGetKSP(dmmg[level]->snes,&ksp);CHKERRQ(ierr);
 
398
  ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
 
399
  mg   = ((MG*)pc->data);
 
400
 
 
401
  for (i=0; i<100; i++) {
 
402
 
 
403
    for (j=level; j>0; j--) {
 
404
 
 
405
      /* Relax residual_fine - F(x_fine) = 0 */
 
406
      for (k=0; k<dmmg[j]->presmooth; k++) {
 
407
        ierr = NLFRelax_DAAD(dmmg[j]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[j]->x);CHKERRQ(ierr);
 
408
      }
 
409
 
 
410
      /* R*(residual_fine - F(x_fine)) */
 
411
      ierr = DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);CHKERRQ(ierr);
 
412
      ierr = VecAYPX(&mone,dmmg[j]->r,dmmg[j]->w);CHKERRQ(ierr);
 
413
 
 
414
      if (j == level || dmmg[j]->monitorall) {
 
415
        /* norm( residual_fine - f(x_fine) ) */
 
416
        ierr = VecNorm(dmmg[j]->w,NORM_2,&norm);CHKERRQ(ierr);
 
417
        if (j == level) {
 
418
          if (norm < dmmg[level]->atol) goto theend; 
 
419
          if (i == 0) {
 
420
            dmmg[level]->rrtol = norm*dmmg[level]->rtol;
 
421
          } else {
 
422
            if (norm < dmmg[level]->rrtol) goto theend;
 
423
          }
 
424
        }
 
425
      }
 
426
 
 
427
      if (dmmg[j]->monitorall) {
 
428
        for (k=0; k<level-j+1; k++) {ierr = PetscPrintf(dmmg[j]->comm,"  ");CHKERRQ(ierr);}
 
429
        ierr = PetscPrintf(dmmg[j]->comm,"FAS function norm %g\n",norm);CHKERRQ(ierr);
 
430
      }
 
431
      ierr = MatRestrict(mg[j]->restrct,dmmg[j]->w,dmmg[j-1]->r);CHKERRQ(ierr); 
 
432
      
 
433
      /* F(R*x_fine) */
 
434
      ierr = VecScatterBegin(dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD,dmmg[j]->inject);CHKERRQ(ierr);
 
435
      ierr = VecScatterEnd(dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD,dmmg[j]->inject);CHKERRQ(ierr);
 
436
      ierr = DMMGFormFunction(0,dmmg[j-1]->x,dmmg[j-1]->w,dmmg[j-1]);CHKERRQ(ierr);
 
437
 
 
438
      /* residual_coarse = F(R*x_fine) + R*(residual_fine - F(x_fine)) */
 
439
      ierr = VecAYPX(&one,dmmg[j-1]->w,dmmg[j-1]->r);CHKERRQ(ierr);
 
440
 
 
441
      /* save R*x_fine into b (needed when interpolating compute x back up */
 
442
      ierr = VecCopy(dmmg[j-1]->x,dmmg[j-1]->b);CHKERRQ(ierr);
 
443
    }
 
444
 
 
445
    for (j=0; j<dmmg[0]->presmooth; j++) {
 
446
      ierr = NLFRelax_DAAD(dmmg[0]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[0]->x);CHKERRQ(ierr);
 
447
    }
 
448
    if (dmmg[0]->monitorall){ 
 
449
      ierr = DMMGFormFunction(0,dmmg[0]->x,dmmg[0]->w,dmmg[0]);CHKERRQ(ierr);
 
450
      ierr = VecAXPY(&mone,dmmg[0]->r,dmmg[0]->w);CHKERRQ(ierr);
 
451
      ierr = VecNorm(dmmg[0]->w,NORM_2,&norm);CHKERRQ(ierr);
 
452
      for (k=0; k<level+1; k++) {ierr = PetscPrintf(dmmg[0]->comm,"  ");CHKERRQ(ierr);}
 
453
      ierr = PetscPrintf(dmmg[0]->comm,"FAS coarse grid function norm %g\n",norm);CHKERRQ(ierr);
 
454
    }
 
455
 
 
456
    for (j=1; j<=level; j++) {
 
457
      /* x_fine = x_fine + R'*(x_coarse - R*x_fine) */
 
458
      ierr = VecAXPY(&mone,dmmg[j-1]->b,dmmg[j-1]->x);CHKERRQ(ierr);
 
459
      ierr = MatInterpolateAdd(mg[j]->restrct,dmmg[j-1]->x,dmmg[j]->x,dmmg[j]->x);CHKERRQ(ierr);
 
460
 
 
461
      if (dmmg[j]->monitorall) {
 
462
        /* norm( F(x_fine) - residual_fine ) */
 
463
        ierr = DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);CHKERRQ(ierr);
 
464
        ierr = VecAXPY(&mone,dmmg[j]->r,dmmg[j]->w);CHKERRQ(ierr);
 
465
        ierr = VecNorm(dmmg[j]->w,NORM_2,&norm);CHKERRQ(ierr);
 
466
        for (k=0; k<level-j+1; k++) {ierr = PetscPrintf(dmmg[j]->comm,"  ");CHKERRQ(ierr);}
 
467
        ierr = PetscPrintf(dmmg[j]->comm,"FAS function norm %g\n",norm);CHKERRQ(ierr);
 
468
      }
 
469
 
 
470
      /* Relax residual_fine - F(x_fine)  = 0 */
 
471
      for (k=0; k<dmmg[j]->postsmooth; k++) {
 
472
        ierr = NLFRelax_DAAD(dmmg[j]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[j]->x);CHKERRQ(ierr);
 
473
      }
 
474
 
 
475
      if (dmmg[j]->monitorall) {
 
476
        /* norm( F(x_fine) - residual_fine ) */
 
477
        ierr = DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);CHKERRQ(ierr);
 
478
        ierr = VecAXPY(&mone,dmmg[j]->r,dmmg[j]->w);CHKERRQ(ierr);
 
479
        ierr = VecNorm(dmmg[j]->w,NORM_2,&norm);CHKERRQ(ierr);
 
480
        for (k=0; k<level-j+1; k++) {ierr = PetscPrintf(dmmg[j]->comm,"  ");CHKERRQ(ierr);}
 
481
        ierr = PetscPrintf(dmmg[j]->comm,"FAS function norm %g\n",norm);CHKERRQ(ierr);
 
482
      }
 
483
    }
 
484
 
 
485
    if (dmmg[level]->monitor){
 
486
      ierr = DMMGFormFunction(0,dmmg[level]->x,dmmg[level]->w,dmmg[level]);CHKERRQ(ierr);
 
487
      ierr = VecNorm(dmmg[level]->w,NORM_2,&norm);CHKERRQ(ierr);
 
488
      ierr = PetscPrintf(dmmg[level]->comm,"%d FAS function norm %g\n",i,norm);CHKERRQ(ierr);
 
489
    }
 
490
  }
 
491
  theend:
 
492
  PetscFunctionReturn(0);
 
493
}
 
494
#endif
 
495
 
342
496
/* ===========================================================================================================*/
343
497
 
344
498
#undef __FUNCT__  
354
508
.   function - the function that defines the nonlinear system
355
509
-   jacobian - optional function to compute Jacobian
356
510
 
357
 
    Options Database:
 
511
    Options Database Keys:
358
512
+    -dmmg_snes_monitor
359
513
.    -dmmg_jacobian_fd
360
514
.    -dmmg_jacobian_ad
361
515
.    -dmmg_jacobian_mf_fd_operator
362
516
.    -dmmg_jacobian_mf_fd
363
517
.    -dmmg_jacobian_mf_ad_operator
364
 
-    -dmmg_jacobian_mf_ad
 
518
.    -dmmg_jacobian_mf_ad
 
519
-    -dmmg_jacobian_period <p> - Indicates how often in the SNES solve the Jacobian is recomputed (on all levels)
 
520
                                 as suggested by Florin Dobrian if p is -1 then Jacobian is computed only on first
 
521
                                 SNES iteration (i.e. -1 is equivalent to infinity) 
365
522
 
366
523
    Level: advanced
367
524
 
368
 
.seealso DMMGCreate(), DMMGDestroy, DMMGSetSLES(), DMMGSetSNESLocal()
 
525
.seealso DMMGCreate(), DMMGDestroy, DMMGSetKSP(), DMMGSetSNESLocal()
369
526
 
370
527
@*/
371
528
int DMMGSetSNES(DMMG *dmmg,int (*function)(SNES,Vec,Vec,void*),int (*jacobian)(SNES,Vec,Mat*,Mat*,MatStructure*,void*))
372
529
{
373
 
  int         ierr,i,nlevels = dmmg[0]->nlevels;
 
530
  int         ierr,size,i,nlevels = dmmg[0]->nlevels,period = 1;
374
531
  PetscTruth  snesmonitor,mffdoperator,mffd,fdjacobian;
375
532
#if defined(PETSC_HAVE_ADIC) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE)
376
533
  PetscTruth  mfadoperator,mfad,adjacobian;
377
534
#endif
378
 
  SLES        sles;
 
535
  KSP        ksp;
379
536
  PetscViewer ascii;
380
537
  MPI_Comm    comm;
381
538
 
436
593
    }
437
594
    
438
595
    if (!dmmg[i]->B) {
439
 
      ierr = DMGetMatrix(dmmg[i]->dm,MATMPIAIJ,&dmmg[i]->B);CHKERRQ(ierr);
 
596
      ierr = MPI_Comm_size(dmmg[i]->comm,&size);CHKERRQ(ierr);
 
597
      ierr = DMGetMatrix(dmmg[i]->dm,MATAIJ,&dmmg[i]->B);CHKERRQ(ierr);
440
598
    } 
441
599
    if (!dmmg[i]->J) {
442
600
      dmmg[i]->J = dmmg[i]->B;
443
601
    }
444
602
 
445
 
    ierr = SNESGetSLES(dmmg[i]->snes,&sles);CHKERRQ(ierr);
446
 
    ierr = DMMGSetUpLevel(dmmg,sles,i+1);CHKERRQ(ierr);
 
603
    ierr = SNESGetKSP(dmmg[i]->snes,&ksp);CHKERRQ(ierr);
 
604
    ierr = DMMGSetUpLevel(dmmg,ksp,i+1);CHKERRQ(ierr);
447
605
    
448
606
    /*
449
607
       if the number of levels is > 1 then we want the coarse solve in the grid sequencing to use LU
451
609
    */
452
610
    if (nlevels > 1 && i == 0) {
453
611
      PC         pc;
454
 
      SLES       csles;
 
612
      KSP        cksp;
455
613
      PetscTruth flg1,flg2,flg3;
456
614
 
457
 
      ierr = SLESGetPC(sles,&pc);CHKERRQ(ierr);
458
 
      ierr = MGGetCoarseSolve(pc,&csles);CHKERRQ(ierr);
459
 
      ierr = SLESGetPC(csles,&pc);CHKERRQ(ierr);
 
615
      ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
 
616
      ierr = MGGetCoarseSolve(pc,&cksp);CHKERRQ(ierr);
 
617
      ierr = KSPGetPC(cksp,&pc);CHKERRQ(ierr);
460
618
      ierr = PetscTypeCompare((PetscObject)pc,PCILU,&flg1);CHKERRQ(ierr);
461
619
      ierr = PetscTypeCompare((PetscObject)pc,PCSOR,&flg2);CHKERRQ(ierr);
462
620
      ierr = PetscTypeCompare((PetscObject)pc,PETSC_NULL,&flg3);CHKERRQ(ierr);
502
660
    ierr = DMGetInterpolationScale(dmmg[i-1]->dm,dmmg[i]->dm,dmmg[i]->R,&dmmg[i]->Rscale);CHKERRQ(ierr);
503
661
  }
504
662
 
 
663
  ierr = PetscOptionsGetInt(PETSC_NULL,"-dmmg_jacobian_period",&period,PETSC_NULL);CHKERRQ(ierr);
 
664
  for (i=0; i<nlevels; i++) {
 
665
    dmmg[i]->updatejacobian       = PETSC_TRUE;
 
666
    dmmg[i]->updatejacobianperiod = period;
 
667
  }
 
668
 
 
669
#if defined(PETSC_HAVE_ADIC) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE)
 
670
  { 
 
671
    PetscTruth flg;
 
672
    ierr = PetscOptionsHasName(PETSC_NULL,"-dmmg_fas",&flg);CHKERRQ(ierr);
 
673
    if (flg) {
 
674
      int newton_its;
 
675
      ierr = PetscOptionsHasName(0,"-fas_view",&flg);CHKERRQ(ierr);
 
676
      for (i=0; i<nlevels; i++) {
 
677
        ierr = NLFCreate_DAAD(&dmmg[i]->nlf);CHKERRQ(ierr);
 
678
        ierr = NLFDAADSetDA_DAAD(dmmg[i]->nlf,(DA)dmmg[i]->dm);CHKERRQ(ierr);
 
679
        ierr = NLFDAADSetCtx_DAAD(dmmg[i]->nlf,dmmg[i]->user);CHKERRQ(ierr);
 
680
        ierr = NLFDAADSetResidual_DAAD(dmmg[i]->nlf,dmmg[i]->r);CHKERRQ(ierr);
 
681
        ierr = VecDuplicate(dmmg[i]->b,&dmmg[i]->w);CHKERRQ(ierr);
 
682
 
 
683
        dmmg[i]->monitor    = PETSC_FALSE;
 
684
        ierr = PetscOptionsHasName(0,"-dmmg_fas_monitor",&dmmg[i]->monitor);CHKERRQ(ierr);
 
685
        dmmg[i]->monitorall = PETSC_FALSE;
 
686
        ierr = PetscOptionsHasName(0,"-dmmg_fas_monitor_all",&dmmg[i]->monitorall);CHKERRQ(ierr);
 
687
        dmmg[i]->presmooth  = 2;
 
688
        ierr = PetscOptionsGetInt(0,"-dmmg_fas_presmooth",&dmmg[i]->presmooth,0);CHKERRQ(ierr);
 
689
        dmmg[i]->postsmooth = 2;
 
690
        ierr = PetscOptionsGetInt(0,"-dmmg_fas_postsmooth",&dmmg[i]->postsmooth,0);CHKERRQ(ierr);
 
691
        dmmg[i]->coarsesmooth = 2;
 
692
        ierr = PetscOptionsGetInt(0,"-dmmg_fas_coarsesmooth",&dmmg[i]->coarsesmooth,0);CHKERRQ(ierr);
 
693
 
 
694
        dmmg[i]->rtol = 1.e-8;
 
695
        ierr = PetscOptionsGetReal(0,"-dmmg_fas_rtol",&dmmg[i]->rtol,0);CHKERRQ(ierr);
 
696
        dmmg[i]->atol = 1.e-50;
 
697
        ierr = PetscOptionsGetReal(0,"-dmmg_fas_atol",&dmmg[i]->atol,0);CHKERRQ(ierr);
 
698
 
 
699
        newton_its = 2;
 
700
        ierr = PetscOptionsGetInt(0,"-dmmg_fas_newton_its",&newton_its,0);CHKERRQ(ierr);
 
701
        ierr = NLFDAADSetNewtonIterations_DAAD(dmmg[i]->nlf,newton_its);CHKERRQ(ierr);
 
702
 
 
703
        if (flg) {
 
704
          if (i == 0) {
 
705
            ierr = PetscPrintf(dmmg[i]->comm,"FAS Solver Parameters\n");CHKERRQ(ierr);
 
706
            ierr = PetscPrintf(dmmg[i]->comm,"  rtol %g atol %g\n",dmmg[i]->rtol,dmmg[i]->atol);CHKERRQ(ierr);
 
707
            ierr = PetscPrintf(dmmg[i]->comm,"             coarsesmooths %d\n",dmmg[i]->coarsesmooth);CHKERRQ(ierr);
 
708
            ierr = PetscPrintf(dmmg[i]->comm,"             Newton iterations %d\n",newton_its);CHKERRQ(ierr);
 
709
          } else {
 
710
            ierr = PetscPrintf(dmmg[i]->comm,"  level %d   presmooths    %d\n",i,dmmg[i]->presmooth);CHKERRQ(ierr);
 
711
            ierr = PetscPrintf(dmmg[i]->comm,"             postsmooths   %d\n",dmmg[i]->postsmooth);CHKERRQ(ierr);
 
712
            ierr = PetscPrintf(dmmg[i]->comm,"             Newton iterations %d\n",newton_its);CHKERRQ(ierr);
 
713
          }
 
714
        }
 
715
        dmmg[i]->solve = DMMGSolveFAS;
 
716
      }
 
717
    }
 
718
  }
 
719
#endif
 
720
   
505
721
  PetscFunctionReturn(0);
506
722
}
507
723
 
519
735
 
520
736
    Level: advanced
521
737
 
522
 
.seealso DMMGCreate(), DMMGDestroy, DMMGSetSLES()
 
738
.seealso DMMGCreate(), DMMGDestroy, DMMGSetKSP()
523
739
 
524
740
@*/
525
741
int DMMGSetInitialGuess(DMMG *dmmg,int (*guess)(SNES,Vec,void*))
533
749
  PetscFunctionReturn(0);
534
750
}
535
751
 
536
 
 
537
752
/*M
538
753
    DMMGSetSNESLocal - Sets the local user function that defines the nonlinear set of equations
539
754
    that will use the grid hierarchy and (optionally) its derivative.
553
768
-   admf_function - the name of the function with an ad_ prefix. This is ignored if ADIC is
554
769
                  not installed
555
770
 
 
771
    Options Database Keys:
 
772
+    -dmmg_snes_monitor
 
773
.    -dmmg_jacobian_fd
 
774
.    -dmmg_jacobian_ad
 
775
.    -dmmg_jacobian_mf_fd_operator
 
776
.    -dmmg_jacobian_mf_fd
 
777
.    -dmmg_jacobian_mf_ad_operator
 
778
.    -dmmg_jacobian_mf_ad
 
779
-    -dmmg_jacobian_period <p> - Indicates how often in the SNES solve the Jacobian is recomputed (on all levels)
 
780
                                 as suggested by Florin Dobrian if p is -1 then Jacobian is computed only on first
 
781
                                 SNES iteration (i.e. -1 is equivalent to infinity) 
 
782
 
 
783
 
556
784
    Level: intermediate
557
785
 
558
786
    Notes: 
563
791
    If ADIC/ADIFOR have not been installed and the Jacobian is not provided, this routine
564
792
    uses finite differencing to approximate the Jacobian.
565
793
 
566
 
.seealso DMMGCreate(), DMMGDestroy, DMMGSetSLES(), DMMGSetSNES()
 
794
.seealso DMMGCreate(), DMMGDestroy, DMMGSetKSP(), DMMGSetSNES()
567
795
 
568
796
M*/
569
797
 
605
833
  ierr = DAGetScatter((DA)dmmg->dm,0,&gtol,0);CHKERRQ(ierr);
606
834
  ierr = VecScatterBegin(u,U,INSERT_VALUES,SCATTER_FORWARD_LOCAL,gtol);CHKERRQ(ierr);
607
835
  ierr = VecScatterEnd(u,U,INSERT_VALUES,SCATTER_FORWARD_LOCAL,gtol);CHKERRQ(ierr);
608
 
 
609
836
  ierr = DAFormFunctioni1((DA)dmmg->dm,i,U,r,dmmg->user);CHKERRQ(ierr);
610
837
  PetscFunctionReturn(0);
611
838
}
631
858
  int ierr,i,nlevels = dmmg[0]->nlevels;
632
859
 
633
860
  PetscFunctionBegin;
634
 
CHKMEMQ;
635
861
  for (i=0; i<nlevels; i++) {
636
 
CHKMEMQ;
637
862
    ierr = DASetLocalFunctioni((DA)dmmg[i]->dm,functioni);CHKERRQ(ierr);
638
 
CHKMEMQ;
639
863
    ierr = DASetLocalAdicFunctioni((DA)dmmg[i]->dm,adi);CHKERRQ(ierr);
640
 
CHKMEMQ;
641
864
    ierr = DASetLocalAdicMFFunctioni((DA)dmmg[i]->dm,adimf);CHKERRQ(ierr);
642
 
CHKMEMQ;
643
865
    ierr = MatSNESMFSetFunctioni(dmmg[i]->J,DMMGFunctioni);CHKERRQ(ierr);
644
 
CHKMEMQ;
645
866
    ierr = MatSNESMFSetFunctioniBase(dmmg[i]->J,DMMGFunctioniBase);CHKERRQ(ierr);    
646
 
CHKMEMQ;
647
867
    ierr = DACreateLocalVector((DA)dmmg[i]->dm,&dmmg[i]->lwork1);CHKERRQ(ierr);
648
868
  }
649
869
  PetscFunctionReturn(0);