Actual source code: ex23.c

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

  4: #include <petscdmda.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:   DM             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:   DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,-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:   DMMGSetFromOptions(dmmg);
 72:   DMMGSetSNESLocali(dmmg,FormFunctionLocali,0,0);
 73:   DMMGSolve(dmmg);
 74:   DMMGDestroy(dmmg);

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

 80:   PetscFinalize();
 81:   return 0;
 82: }

 86: /*
 87:      This local function acts on the ghosted version of U (accessed via DMGetLocalVector())
 88:      BUT the global, nonghosted version of FU

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

101:   DMGetLocalVector(da,&vu);
102:   DMGlobalToLocalBegin(da,U,INSERT_VALUES,vu);
103:   DMGlobalToLocalEnd(da,U,INSERT_VALUES,vu);

105:   DMDAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
106:   DMDAGetInfo(da,0,&N,0,0,0,0,0,0,0,0,0,0,0);
107:   DMDAVecGetArray(da,vu,&u);
108:   DMDAVecGetArray(da,FU,&fu);
109:   d    = N-1.0;
110:   h    = 1.0/d;

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

118:   DMDAVecRestoreArray(da,vu,&u);
119:   DMDAVecRestoreArray(da,FU,&fu);
120:   DMRestoreLocalVector(da,&vu);
121:   PetscLogFlops(9.0*N);
122:   return(0);
123: }

127: PetscErrorCode FormFunctionLocali(DMDALocalInfo *info,MatStencil *pt,PetscScalar *u,PetscScalar *fu,void* dummy)
128: {
129:   PetscInt     i,N = info->mx;
130:   PetscScalar  d,h;

133:   d    = N-1.0;
134:   h    = 1.0/d;

136:   i = pt->i;
137:   if      (i == 0)   *fu = 2.0*d*(u[0] - .25) + h*u[0]*u[0];
138:   else if (i == N-1) *fu = 2.0*d*u[N-1] + h*u[N-1]*u[N-1];
139:   else               *fu = -(d*(u[i+1] - 2.0*u[i] + u[i-1]) - 2.0*h) + h*u[i]*u[i];

141:   return(0);
142: }