Actual source code: p1.c

  2: static char help[] = "Model single-physics solver. Modified from ex19.c \n\\n";

  4: /* ------------------------------------------------------------------------

  6:     See ex19.c for discussion of the problem 

  8:   ------------------------------------------------------------------------- */
 9:  #include mp.h


 17: int main(int argc,char **argv)
 18: {
 19:   DMMG           *dmmg;               /* multilevel grid structure */
 20:   AppCtx         user;                /* user-defined work context */
 21:   PetscInt       mx,my,its;
 23:   MPI_Comm       comm;
 24:   SNES           snes;
 25:   DA             da1;

 27:   PetscInitialize(&argc,&argv,(char *)0,help);
 28:   PetscLogEventRegister("FormFunc1", 0,&EVENT_FORMFUNCTIONLOCAL1);
 29:   comm = PETSC_COMM_WORLD;

 31:   /* Problem parameters (velocity of lid, prandtl, and grashof numbers) */
 32:   PetscOptionsGetReal(PETSC_NULL,"-lidvelocity",&user.lidvelocity,PETSC_NULL);
 33:   PetscOptionsGetReal(PETSC_NULL,"-prandtl",&user.prandtl,PETSC_NULL);
 34:   PetscOptionsGetReal(PETSC_NULL,"-grashof",&user.grashof,PETSC_NULL);

 36:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 37:      Create user context, set problem data, create vector data structures.
 38:      Also, compute the initial guess.
 39:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 41:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 42:      Setup Physics 1: 
 43:         - Lap(U) - Grad_y(Omega) = 0
 44:         - Lap(V) + Grad_x(Omega) = 0
 45:         - Lap(Omega) + Div([U*Omega,V*Omega]) - GR*Grad_x(T) = 0
 46:         where T is given by the given x.temp
 47:         - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 48:   DACreate2d(comm,DA_NONPERIODIC,DA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,3,1,0,0,&da1);
 49:   DASetFieldName(da1,0,"x-velocity");
 50:   DASetFieldName(da1,1,"y-velocity");
 51:   DASetFieldName(da1,2,"Omega");

 53:   /* Create the solver object and attach the grid/physics info */
 54:   DMMGCreate(comm,1,&user,&dmmg);
 55:   DMMGSetDM(dmmg,(DM)da1);
 56:   DMMGSetISColoringType(dmmg,IS_COLORING_GLOBAL);

 58:   DMMGSetInitialGuess(dmmg,FormInitialGuess);
 59:   DMMGSetSNES(dmmg,FormFunction,0);
 60:   DMMGSetFromOptions(dmmg);

 62:   DAGetInfo(da1,PETSC_NULL,&mx,&my,0,0,0,0,0,0,0,0);
 63:   user.lidvelocity = 1.0/(mx*my);
 64:   user.prandtl     = 1.0;
 65:   user.grashof     = 1.0;

 67:   /* Solve the nonlinear system */
 68:   DMMGSolve(dmmg);
 69:   snes = DMMGGetSNES(dmmg);
 70:   SNESGetIterationNumber(snes,&its);
 71:   PetscPrintf(comm,"Physics 1: Number of Newton iterations = %D\n\n", its);

 73:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 74:      Free spaces 
 75:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 76:   DADestroy(da1);
 77:   DMMGDestroy(dmmg);
 78:   PetscFinalize();
 79:   return 0;
 80: }

 84: PetscErrorCode FormInitialGuess(DMMG dmmg,Vec X)
 85: {
 87:   DA             da1 = (DA)dmmg->dm;
 88:   Field1         **x1;
 89:   DALocalInfo    info1;


 93:   /* Access the array inside of X */
 94:   DAVecGetArray(da1,X,(void**)&x1);

 96:   DAGetLocalInfo(da1,&info1);

 98:   /* Evaluate local user provided function */
 99:   FormInitialGuessLocal1(&info1,x1);

101:   DAVecRestoreArray(da1,X,(void**)&x1);
102:   return(0);
103: }

107: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ctx)
108: {
110:   DMMG           dmmg = (DMMG)ctx;
111:   AppCtx         *user = (AppCtx*)dmmg->user;
112:   DA             da1 = (DA)dmmg->dm;
113:   DALocalInfo    info1;
114:   Field1         **x1,**f1;
115:   Vec            X1;


119:   DAGetLocalInfo(da1,&info1);

121:   /* Get local vectors to hold ghosted parts of X */
122:   DAGetLocalVector(da1,&X1);
123:   DAGlobalToLocalBegin(da1,X,INSERT_VALUES,X1);
124:   DAGlobalToLocalEnd(da1,X,INSERT_VALUES,X1);

126:   /* Access the arrays inside X1 */
127:   DAVecGetArray(da1,X1,(void**)&x1);

129:   /* Access the subvectors in F. 
130:      These are not ghosted so directly access the memory locations in F */
131:   DAVecGetArray(da1,F,(void**)&f1);

133:   /* Evaluate local user provided function */
134:   FormFunctionLocal1(&info1,x1,0,f1,(void**)user);


137:   DAVecRestoreArray(da1,F,(void**)&f1);
138:   DAVecRestoreArray(da1,X1,(void**)&x1);
139:   DARestoreLocalVector(da1,&X1);
140:   return(0);
141: }