Actual source code: ex23.c

      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: }