1: /*$Id: ex23.c,v 1.13 2001/08/22 20:35:33 bsmith Exp $*/

  3: static char help[] = "Solves PDE problem from ex22.cnn";

 5:  #include petscda.h
 6:  #include petscpf.h
 7:  #include petscsnes.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;

 30: extern int FormFunction(SNES,Vec,Vec,void*);
 31: extern int FormFunctionLocali(DALocalInfo*,MatStencil*,PetscScalar*,PetscScalar*,void*);

 33: #undef __FUNCT__
 35: int main(int argc,char **argv)
 36: {
 37:   int     ierr;
 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: int FormFunction(SNES snes,Vec U,Vec FU,void* dummy)
 90: {
 91:   DMMG    dmmg = (DMMG)dummy;
 92:   int     ierr,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,(void**)&u);
105:   DAVecGetArray(da,FU,(void**)&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,(void**)&u);
116:   DAVecRestoreArray(da,FU,(void**)&fu);
117:   DARestoreLocalVector(da,&vu);
118:   PetscLogFlops(9*N);
119:   return(0);
120: }

122: int FormFunctionLocali(DALocalInfo *info,MatStencil *pt,PetscScalar *u,PetscScalar *fu,void* dummy)
123: {
124:   int          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: }