2: static char help[] = "Solves PDE problem from ex22.c\n\n";

 4:  #include petscda.h
 5:  #include petscpf.h
 6:  #include petscsnes.h
 7:  #include petscdmmg.h

  9: /*

 11:        In this example the PDE is 
 12:                              Uxx + U^2 = 2, 
 13:                             u(0) = .25
 14:                             u(1) = 0

 16:        The exact solution for u is given by u(x) = x*x - 1.25*x + .25

 18:        Use the usual centered finite differences.


 21:        See ex22.c for a design optimization problem

 23: */

 25: typedef struct {
 26:   PetscViewer  u_viewer;
 27:   PetscViewer  fu_viewer;
 28: } UserCtx;


 35: int main(int argc,char **argv)
 36: {
 38:   UserCtx        user;
 39:   DA             da;
 40:   DMMG           *dmmg;

 42:   PetscInitialize(&argc,&argv,PETSC_NULL,help);

 44:   /* Hardwire several options; can be changed at command line */
 45:   PetscOptionsSetValue("-dmmg_grid_sequence",PETSC_NULL);
 46:   PetscOptionsSetValue("-ksp_type","fgmres");
 47:   PetscOptionsSetValue("-ksp_max_it","5");
 48:   PetscOptionsSetValue("-pc_mg_type","full");
 49:   PetscOptionsSetValue("-mg_coarse_ksp_type","gmres");
 50:   PetscOptionsSetValue("-mg_levels_ksp_type","gmres");
 51:   PetscOptionsSetValue("-mg_coarse_ksp_max_it","3");
 52:   PetscOptionsSetValue("-mg_levels_ksp_max_it","3");
 53:   /* PetscOptionsSetValue("-snes_mf_type","wp"); */
 54:   /* PetscOptionsSetValue("-snes_mf_compute_norma","no"); */
 55:   /* PetscOptionsSetValue("-snes_mf_compute_normu","no"); */
 56:   PetscOptionsSetValue("-snes_ls","basic");
 57:   /*  PetscOptionsSetValue("-dmmg_jacobian_mf_fd",0); */
 58:   /*  PetscOptionsSetValue("-snes_ls","basicnonorms"); */
 59:   PetscOptionsInsert(&argc,&argv,PETSC_NULL);
 60: 
 61:   /* Create a global vector from a da arrays */
 62:   DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,-5,1,1,PETSC_NULL,&da);

 64:   /* create graphics windows */
 65:   PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"u - state variables",-1,-1,-1,-1,&user.u_viewer);
 66:   PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"fu - discretization of function",-1,-1,-1,-1,&user.fu_viewer);

 68:   /* create nonlinear multi-level solver */
 69:   DMMGCreate(PETSC_COMM_WORLD,2,&user,&dmmg);
 70:   DMMGSetDM(dmmg,(DM)da);
 71:   DMMGSetSNES(dmmg,FormFunction,PETSC_NULL);
 72:   DMMGSetSNESLocali(dmmg,FormFunctionLocali,0,0);
 73:   DMMGSolve(dmmg);
 74:   DMMGDestroy(dmmg);

 76:   DADestroy(da);
 77:   PetscViewerDestroy(user.u_viewer);
 78:   PetscViewerDestroy(user.fu_viewer);

 80:   PetscFinalize();
 81:   return 0;
 82: }
 83: 
 84: /*
 85:      This local function acts on the ghosted version of U (accessed via DAGetLocalVector())
 86:      BUT the global, nonghosted version of FU

 88: */
 89: PetscErrorCode FormFunction(SNES snes,Vec U,Vec FU,void* dummy)
 90: {
 91:   DMMG           dmmg = (DMMG)dummy;
 93:   PetscInt       xs,xm,i,N;
 94:   PetscScalar    *u,*fu,d,h;
 95:   Vec            vu;
 96:   DA             da = (DA) dmmg->dm;

 99:   DAGetLocalVector(da,&vu);
100:   DAGlobalToLocalBegin(da,U,INSERT_VALUES,vu);
101:   DAGlobalToLocalEnd(da,U,INSERT_VALUES,vu);

103:   DAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
104:   DAGetInfo(da,0,&N,0,0,0,0,0,0,0,0,0);
105:   DAVecGetArray(da,vu,&u);
106:   DAVecGetArray(da,FU,&fu);
107:   d    = N-1.0;
108:   h    = 1.0/d;

110:   for (i=xs; i<xs+xm; i++) {
111:     if      (i == 0)   fu[0]   = 2.0*d*(u[0] - .25) + h*u[0]*u[0];
112:     else if (i == N-1) fu[N-1] = 2.0*d*u[N-1] + h*u[N-1]*u[N-1];
113:     else               fu[i]   = -(d*(u[i+1] - 2.0*u[i] + u[i-1]) - 2.0*h) + h*u[i]*u[i];
114:   }

116:   DAVecRestoreArray(da,vu,&u);
117:   DAVecRestoreArray(da,FU,&fu);
118:   DARestoreLocalVector(da,&vu);
119:   PetscLogFlops(9*N);
120:   return(0);
121: }

123: PetscErrorCode FormFunctionLocali(DALocalInfo *info,MatStencil *pt,PetscScalar *u,PetscScalar *fu,void* dummy)
124: {
125:   PetscInt     i,N = info->mx;
126:   PetscScalar  d,h;

129:   d    = N-1.0;
130:   h    = 1.0/d;

132:   i = pt->i;
133:   if      (i == 0)   *fu = 2.0*d*(u[0] - .25) + h*u[0]*u[0];
134:   else if (i == N-1) *fu = 2.0*d*u[N-1] + h*u[N-1]*u[N-1];
135:   else               *fu = -(d*(u[i+1] - 2.0*u[i] + u[i-1]) - 2.0*h) + h*u[i]*u[i];

137:   return(0);
138: }