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