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