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_normu","no"); */
 55:   PetscOptionsSetValue("-snes_ls","basic");
 56:   /*  PetscOptionsSetValue("-dmmg_jacobian_mf_fd",0); */
 57:   /*  PetscOptionsSetValue("-snes_ls","basicnonorms"); */
 58:   PetscOptionsInsert(&argc,&argv,PETSC_NULL);
 59: 
 60:   /* Create a global vector from a da arrays */
 61:   DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,-5,1,1,PETSC_NULL,&da);

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

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

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

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

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

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

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

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

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

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

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

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

136:   return(0);
137: }