Actual source code: ex31.c

  2: static char help[] = "\n\n";

  4: /*T
  5:    Concepts: SNES^solving a system of nonlinear equations (parallel multicomponent example);
  6:    Concepts: DA^using distributed arrays;
  7:    Concepts: multicomponent
  8:    Processors: n
  9: T*/

 11: /* 
 12:    Include "petscda.h" so that we can use distributed arrays (DAs).
 13:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
 14:    file automatically includes:
 15:      petsc.h       - base PETSc routines   petscvec.h - vectors
 16:      petscsys.h    - system routines       petscmat.h - matrices
 17:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 18:      petscviewer.h - viewers               petscpc.h  - preconditioners
 19:      petscksp.h   - linear solvers 
 20: */
 21:  #include petscsnes.h
 22:  #include petscda.h
 23:  #include petscdmmg.h

 25: /* 
 26:    User-defined routines and data structures
 27: */
 28: typedef struct {
 29:   PetscScalar h,uh;
 30: } Field;


 37: int main(int argc,char **argv)
 38: {
 39:   DMMG           *dmmg;               /* multilevel grid structure */
 41:   MPI_Comm       comm;
 42:   DA             da;

 44:   PetscInitialize(&argc,&argv,(char *)0,help);
 45:   comm = PETSC_COMM_WORLD;


 48:   PreLoadBegin(PETSC_TRUE,"SetUp");
 49:     DMMGCreate(comm,2,0,&dmmg);

 51:     /*
 52:       Create distributed array multigrid object (DMMG) to manage parallel grid and vectors
 53:       for principal unknowns (x) and governing residuals (f)

 55:       The problem is actually not periodic; but we declare it as periodic so that we have
 56:       two ghost points at each end where we can put "ghost" boundary conditions.
 57:     */
 58:     DACreate1d(comm,DA_XPERIODIC,-6,2,2,0,&da);
 59:     DMMGSetDM(dmmg,(DM)da);
 60:     DADestroy(da);

 62:     /* 
 63:      Problem parameters 
 64:     */
 65:     DASetFieldName(DMMGGetDA(dmmg),0,"h");
 66:     DASetFieldName(DMMGGetDA(dmmg),1,"uh");

 68:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 69:        Create user context, set problem data, create vector data structures.
 70:        Also, compute the initial guess.
 71:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 73:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 74:        Create nonlinear solver context

 76:        Process adiC(36): FormFunctionLocal 
 77:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 78:     DMMGSetSNESLocal(dmmg,FormFunctionLocal,0,ad_FormFunctionLocal,admf_FormFunctionLocal);

 80:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 81:        Solve the nonlinear system
 82:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 83:     DMMGSetInitialGuessLocal(dmmg,(PetscErrorCode (*)(DMMG,void*)) FormInitialGuessLocal);

 85:   PreLoadStage("Solve");
 86:     DMMGSolve(dmmg);
 87:     DMMGSetInitialGuess(dmmg,PETSC_NULL);

 89:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 90:        Free work space.  All PETSc objects should be destroyed when they
 91:        are no longer needed.
 92:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 94:     DMMGDestroy(dmmg);
 95:   PreLoadEnd();

 97:   PetscFinalize();
 98:   return 0;
 99: }

101: /* ------------------------------------------------------------------- */


106: /* 
107:    FormInitialGuessLocal - Forms initial approximation.

109:    Input Parameters:
110:    user - user-defined application context
111:    X - vector

113:    Output Parameter:
114:    X - vector
115:  */
116: PetscErrorCode FormInitialGuessLocal(DALocalInfo* info,Field x[])
117: {
118:   PetscInt       i;
120:   PetscReal      dhx,hx;

123:   dhx = (PetscReal)(info->mx-1);
124:   hx = 1.0/dhx;

126:   for (i=info->xs; i<info->xs+info->xm; i++) {
127:     x[i].h   = 22.3;
128:     x[i].uh  = 0.0;
129:   }
130:   return(0);
131: }
132: 
133: PetscErrorCode FormFunctionLocal(DALocalInfo *info,Field *x,Field *f,void *ptr)
134:  {
136:   PetscInt       i;
137:   PetscReal      hx,dhx;

140:   /* 
141:      Define mesh intervals ratios for uniform grid.
142:      [Note: FD formulae below are normalized by multiplying through by
143:      local volume element to obtain coefficients O(1) in two dimensions.]
144:   */
145:   dhx = (PetscReal)(info->mx-1);
146:   hx = 1.0/dhx;

148:   /* 
149:      Put the ghost boundary conditions at the left and right end; we have the available space because
150:      we declared the problem with periodic boundary conditions and two sets of ghost points.

152:      This is actually a staggered grid with the h's living on cell vertices and uh on cell centers,
153:      the "extra" x[mx-1].uh is set to zero.
154:        
155:   */
156:   if (info->xs = 0) {
157:     x[-2].uh  = x[-1].uh = 0.0;
158:     x[-2].h   = x[-1].h  = x[0].h;
159:   } else if (info->xs+>info->xm == mx) {
160:     x[mx+1].uh = x[mx].uh = x[mx-1].uh  = 0.0;
161:     x[mx+1].h = x[mx].h  = x[mx-1].h;
162:   }

164:   for (i=info->xs; i<info->xs+info->xm; i++) {
165:     f[i].h   = x[i].h;
166:     f[i].uh  = x[i].uh;
167:   }
168:   return(0);
169: }