Actual source code: mp.c

  2: static char help[] = "Model multi-physics solver. Modified from src/snes/examples/tutorials/ex19.c \n\\n";

  4: /* ------------------------------------------------------------------------
  5:     See ex19.c for discussion of the problem 

  7:     Examples of command line options:
  8:       ./mp -dmmg_jacobian_mf_fd_operator
  9:       ./mp -dmcomposite_dense_jacobian #inefficient, but compute entire Jacobian for testing
 10:   ----------------------------------------------------------------------------------------- */
 11:  #include mp.h


 20: int main(int argc,char **argv)
 21: {
 22:   DMMG           *dmmg_comp;          /* multilevel grid structure */
 23:   AppCtx         user;                /* user-defined work context */
 24:   PetscInt       mx,my,its;
 26:   MPI_Comm       comm;
 27:   SNES           snes;
 28:   DA             da1,da2;
 29:   DMComposite    pack;
 30:   PetscTruth     couple = PETSC_FALSE;

 32:   PetscInitialize(&argc,&argv,(char *)0,help);
 33:   PetscLogEventRegister("FormFunc1", 0,&EVENT_FORMFUNCTIONLOCAL1);
 34:   PetscLogEventRegister("FormFunc2", 0,&EVENT_FORMFUNCTIONLOCAL2);
 35:   comm = PETSC_COMM_WORLD;

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

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

 54:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 55:      Setup Physics 2: 
 56:         - Lap(T) + PR*Div([U*T,V*T]) = 0        
 57:         where U and V are given by the given x.u and x.v
 58:         - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 59:   DACreate2d(comm,DA_NONPERIODIC,DA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da2);
 60:   DASetFieldName(da2,0,"temperature");

 62:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 63:     Create the DMComposite object to manage the two grids/physics. 
 64:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 65:   DMCompositeCreate(comm,&pack);
 66:   DMCompositeAddDM(pack,(DM)da1);
 67:   DMCompositeAddDM(pack,(DM)da2);

 69:   PetscOptionsHasName(PETSC_NULL,"-couple",&couple);
 70:   if (couple) {
 71:     DMCompositeSetCoupling(pack,FormCoupleLocations);
 72:   }

 74:   /* Create the solver object and attach the grid/physics info */
 75:   DMMGCreate(comm,1,&user,&dmmg_comp);
 76:   DMMGSetDM(dmmg_comp,(DM)pack);
 77:   DMMGSetISColoringType(dmmg_comp,IS_COLORING_GLOBAL);

 79:   DMMGSetInitialGuess(dmmg_comp,FormInitialGuessComp);
 80:   DMMGSetSNES(dmmg_comp,FormFunctionComp,0);
 81:   DMMGSetFromOptions(dmmg_comp);
 82:   DMMGSetUp(dmmg_comp);

 84:   /* Problem parameters (velocity of lid, prandtl, and grashof numbers) */
 85:   DAGetInfo(da1,PETSC_NULL,&mx,&my,0,0,0,0,0,0,0,0);
 86:   user.lidvelocity = 1.0/(mx*my);
 87:   user.prandtl     = 1.0;
 88:   user.grashof     = 1000.0;
 89:   PetscOptionsGetReal(PETSC_NULL,"-lidvelocity",&user.lidvelocity,PETSC_NULL);
 90:   PetscOptionsGetReal(PETSC_NULL,"-prandtl",&user.prandtl,PETSC_NULL);
 91:   PetscOptionsGetReal(PETSC_NULL,"-grashof",&user.grashof,PETSC_NULL);

 93:   /* Solve the nonlinear system */
 94:   DMMGSolve(dmmg_comp);
 95:   snes = DMMGGetSNES(dmmg_comp);
 96:   SNESGetIterationNumber(snes,&its);
 97:   PetscPrintf(comm,"Composite Physics: Number of Newton iterations = %D\n\n", its);

 99:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100:      Free spaces 
101:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
102:   DMCompositeDestroy(pack);
103:   DADestroy(da1);
104:   DADestroy(da2);
105:   DMMGDestroy(dmmg_comp);
106:   PetscFinalize();
107:   return 0;
108: }

112: /* 
113:    FormInitialGuessComp - 
114:               Forms the initial guess for the composite model
115:               Unwraps the global solution vector and passes its local pieces into the user functions
116:  */
117: PetscErrorCode FormInitialGuessComp(DMMG dmmg,Vec X)
118: {
120:   AppCtx         *user = (AppCtx*)dmmg->user;
121:   DMComposite    dm = (DMComposite)dmmg->dm;
122:   Vec            X1,X2;
123:   Field1         **x1;
124:   Field2         **x2;
125:   DALocalInfo    info1,info2;
126:   DA             da1,da2;

129:   DMCompositeGetEntries(dm,&da1,&da2);
130:   /* Access the subvectors in X */
131:   DMCompositeGetAccess(dm,X,&X1,&X2);
132:   /* Access the arrays inside the subvectors of X */
133:   DAVecGetArray(da1,X1,(void**)&x1);
134:   DAVecGetArray(da2,X2,(void**)&x2);

136:   DAGetLocalInfo(da1,&info1);
137:   DAGetLocalInfo(da2,&info2);

139:   /* Evaluate local user provided function */
140:   FormInitialGuessLocal1(&info1,x1);
141:   FormInitialGuessLocal2(&info2,x2,user);

143:   DAVecRestoreArray(da1,X1,(void**)&x1);
144:   DAVecRestoreArray(da2,X2,(void**)&x2);
145:   DMCompositeRestoreAccess(dm,X,&X1,&X2);
146:   return(0);
147: }

151: /* 
152:    FormFunctionComp  - Unwraps the input vector and passes its local ghosted pieces into the user function
153: */
154: PetscErrorCode FormFunctionComp(SNES snes,Vec X,Vec F,void *ctx)
155: {
157:   DMMG           dmmg = (DMMG)ctx;
158:   AppCtx         *user = (AppCtx*)dmmg->user;
159:   DMComposite    dm = (DMComposite)dmmg->dm;
160:   DALocalInfo    info1,info2;
161:   DA             da1,da2;
162:   Field1         **x1,**f1;
163:   Field2         **x2,**f2;
164:   Vec            X1,X2,F1,F2;


168:   DMCompositeGetEntries(dm,&da1,&da2);
169:   DAGetLocalInfo(da1,&info1);
170:   DAGetLocalInfo(da2,&info2);

172:   /* Get local vectors to hold ghosted parts of X */
173:   DMCompositeGetLocalVectors(dm,&X1,&X2);
174:   DMCompositeScatter(dm,X,X1,X2);

176:   /* Access the arrays inside the subvectors of X */
177:   DAVecGetArray(da1,X1,(void**)&x1);
178:   DAVecGetArray(da2,X2,(void**)&x2);

180:   /* Access the subvectors in F. 
181:      These are not ghosted so directly access the memory locations in F */
182:   DMCompositeGetAccess(dm,F,&F1,&F2);

184:   /* Access the arrays inside the subvectors of F */
185:   DAVecGetArray(da1,F1,(void**)&f1);
186:   DAVecGetArray(da2,F2,(void**)&f2);

188:   /* Evaluate local user provided function */
189:   FormFunctionLocal1(&info1,x1,x2,f1,(void**)user);
190:   FormFunctionLocal2(&info2,x1,x2,f2,(void**)user);

192:   DAVecRestoreArray(da1,F1,(void**)&f1);
193:   DAVecRestoreArray(da2,F2,(void**)&f2);
194:   DMCompositeRestoreAccess(dm,F,&F1,&F2);
195:   DAVecRestoreArray(da1,X1,(void**)&x1);
196:   DAVecRestoreArray(da2,X2,(void**)&x2);
197:   DMCompositeRestoreLocalVectors(dm,&X1,&X2);
198:   return(0);
199: }

203: /* 
204:    Computes the coupling between DA1 and DA2. This determines the location of each coupling between DA1 and DA2.
205: */
206: PetscErrorCode FormCoupleLocations(DMComposite dmcomposite,Mat A,PetscInt *dnz,PetscInt *onz,PetscInt __rstart,PetscInt __nrows,PetscInt __start,PetscInt __end)
207: {
208:   PetscInt       i,j,cols[2],istart,jstart,in,jn,row,col,M;
210:   DA             da1,da2;

213:    DMCompositeGetEntries(dmcomposite,&da1,&da2);
214:    DAGetInfo(da1,0,&M,0,0,0,0,0,0,0,0,0);
215:   DAGetCorners(da1,&istart,&jstart,PETSC_NULL,&in,&jn,PETSC_NULL);

217:   /* coupling from physics 1 to physics 2 */
218:   row = __rstart + 2;  /* global location of first omega on this process */
219:   col = __rstart + 3*in*jn;  /* global location of first temp on this process */
220:   for (j=jstart; j<jstart+jn; j++) {
221:     for (i=istart; i<istart+in; i++) {

223:       /* each omega is coupled to the temp to the left and right */
224:       if (i == 0) {
225:         cols[0] = col + 1;
226:         MatPreallocateLocation(A,row,1,cols,dnz,onz);
227:       } else if (i == M-1) {
228:         cols[0] = col - 1;
229:         MatPreallocateLocation(A,row,1,cols,dnz,onz);
230:       } else {
231:         cols[0] = col - 1;
232:         cols[1] = col + 1;
233:         MatPreallocateLocation(A,row,2,cols,dnz,onz);
234:       }
235:       row += 3;
236:       col += 1;
237:     }
238:   }

240:   /* coupling from physics 2 to physics 1 */
241:   col = __rstart;  /* global location of first u on this process */
242:   row = __rstart + 3*in*jn;  /* global location of first temp on this process */
243:   for (j=jstart; j<jstart+jn; j++) {
244:     for (i=istart; i<istart+in; i++) {

246:       /* temp is coupled to both u and v at each point */
247:       cols[0] = col;
248:       cols[1] = col + 1;
249:       MatPreallocateLocation(A,row,2,cols,dnz,onz);
250:       row += 1;
251:       col += 3;
252:     }
253:   }

255:   return(0);
256: }