1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
|
static char help[] = "Model multi-physics solver. Modified from ex19.c \n\\n";
/*T
Concepts: SNES^solving a system of nonlinear equations (parallel multicomponent example);
Concepts: DA^using distributed arrays;
Concepts: multicomponent
Processors: n
T*/
/* ------------------------------------------------------------------------
See ex19.c for discussion of the problem
------------------------------------------------------------------------- */
#include "mp.h"
extern PetscErrorCode FormInitialGuessComp(DMMG,Vec);
extern PetscErrorCode FormFunctionComp(SNES,Vec,Vec,void*);
#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **argv)
{
DMMG *dmmg_comp; /* multilevel grid structure */
AppCtx user; /* user-defined work context */
PetscInt mx,my,its,dof,nlevels=1;
PetscErrorCode ierr;
MPI_Comm comm;
SNES snes;
DA da1,da2;
DMComposite pack;
PetscInitialize(&argc,&argv,(char *)0,help);
comm = PETSC_COMM_WORLD;
ierr = PetscOptionsGetInt(PETSC_NULL,"-nlevels",&nlevels,PETSC_NULL);CHKERRQ(ierr);
/* Problem parameters (velocity of lid, prandtl, and grashof numbers) */
ierr = PetscOptionsGetReal(PETSC_NULL,"-lidvelocity",&user.lidvelocity,PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetReal(PETSC_NULL,"-prandtl",&user.prandtl,PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetReal(PETSC_NULL,"-grashof",&user.grashof,PETSC_NULL);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create user context, set problem data, create vector data structures.
Also, compute the initial guess.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Setup Physics 1:
- Lap(U) - Grad_y(Omega) = 0
- Lap(V) + Grad_x(Omega) = 0
- Lap(Omega) + Div([U*Omega,V*Omega]) - GR*Grad_x(T) = 0
where T is given by the given x.temp
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
dof = 3;
ierr = DACreate2d(comm,DA_NONPERIODIC,DA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,dof,1,0,0,&da1);CHKERRQ(ierr);
ierr = DASetFieldName(da1,0,"x-velocity");CHKERRQ(ierr);
ierr = DASetFieldName(da1,1,"y-velocity");CHKERRQ(ierr);
ierr = DASetFieldName(da1,2,"Omega");CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Setup Physics 2:
- Lap(T) + PR*Div([U*T,V*T]) = 0
where U and V are given by the given x.u and x.v
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
dof = 1;
ierr = DACreate2d(comm,DA_NONPERIODIC,DA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,dof,1,0,0,&da2);CHKERRQ(ierr);
ierr = DASetFieldName(da2,0,"temperature");CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create the DMComposite object to manage the two grids/physics.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = DMCompositeCreate(comm,&pack);CHKERRQ(ierr);
ierr = DMCompositeAddDM(pack,(DM)da1);CHKERRQ(ierr);
ierr = DMCompositeAddDM(pack,(DM)da2);CHKERRQ(ierr);
/* Create the solver object and attach the grid/physics info */
ierr = DMMGCreate(comm,nlevels,&user,&dmmg_comp);CHKERRQ(ierr);
ierr = DMMGSetDM(dmmg_comp,(DM)pack);CHKERRQ(ierr);
ierr = DMMGSetISColoringType(dmmg_comp,IS_COLORING_GLOBAL);CHKERRQ(ierr);
ierr = DMMGSetInitialGuess(dmmg_comp,FormInitialGuessComp);CHKERRQ(ierr);
ierr = DMMGSetSNES(dmmg_comp,FormFunctionComp,0);CHKERRQ(ierr);
ierr = DAGetInfo(da1,PETSC_NULL,&mx,&my,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
user.lidvelocity = 1.0/(mx*my);
user.prandtl = 1.0;
user.grashof = 1.0;
/* Solve the nonlinear system */
ierr = DMMGSolve(dmmg_comp);CHKERRQ(ierr);
snes = DMMGGetSNES(dmmg_comp);
ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr);
ierr = PetscPrintf(comm,"Composite Physics: Number of Newton iterations = %D\n\n", its);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Free spaces
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = DMCompositeDestroy(pack);CHKERRQ(ierr);
ierr = DADestroy(da1);CHKERRQ(ierr);
ierr = DADestroy(da2);CHKERRQ(ierr);
ierr = DMMGDestroy(dmmg_comp);CHKERRQ(ierr);
ierr = PetscFinalize();CHKERRQ(ierr);
return 0;
}
#undef __FUNCT__
#define __FUNCT__ "FormInitialGuessComp"
/*
FormInitialGuessComp -
Forms the initial guess for the composite model
Unwraps the global solution vector and passes its local pieces into the user functions
*/
PetscErrorCode FormInitialGuessComp(DMMG dmmg,Vec X)
{
PetscErrorCode ierr;
AppCtx *user = (AppCtx*)dmmg->user;
DMComposite dm = (DMComposite)dmmg->dm;
Vec X1,X2;
Field1 **x1;
Field2 **x2;
DALocalInfo info1,info2;
DA da1,da2;
PetscFunctionBegin;
ierr = DMCompositeGetEntries(dm,&da1,&da2);CHKERRQ(ierr);
/* Access the subvectors in X */
ierr = DMCompositeGetAccess(dm,X,&X1,&X2);CHKERRQ(ierr);
/* Access the arrays inside the subvectors of X */
ierr = DAVecGetArray(da1,X1,(void**)&x1);CHKERRQ(ierr);
ierr = DAVecGetArray(da2,X2,(void**)&x2);CHKERRQ(ierr);
ierr = DAGetLocalInfo(da1,&info1);CHKERRQ(ierr);
ierr = DAGetLocalInfo(da2,&info2);CHKERRQ(ierr);
/* Evaluate local user provided function */
ierr = FormInitialGuessLocal1(&info1,x1);CHKERRQ(ierr);
ierr = FormInitialGuessLocal2(&info2,x2,user);CHKERRQ(ierr);
ierr = DAVecRestoreArray(da1,X1,(void**)&x1);CHKERRQ(ierr);
ierr = DAVecRestoreArray(da2,X2,(void**)&x2);CHKERRQ(ierr);
ierr = DMCompositeRestoreAccess(dm,X,&X1,&X2);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "FormFunctionComp"
/*
FormFunctionComp - Unwraps the input vector and passes its local ghosted pieces into the user function
*/
PetscErrorCode FormFunctionComp(SNES snes,Vec X,Vec F,void *ctx)
{
PetscErrorCode ierr;
DMMG dmmg = (DMMG)ctx;
AppCtx *user = (AppCtx*)dmmg->user;
DMComposite dm = (DMComposite)dmmg->dm;
DALocalInfo info1,info2;
DA da1,da2;
Field1 **x1,**f1;
Field2 **x2,**f2;
Vec X1,X2,F1,F2;
PetscFunctionBegin;
ierr = DMCompositeGetEntries(dm,&da1,&da2);CHKERRQ(ierr);
ierr = DAGetLocalInfo(da1,&info1);CHKERRQ(ierr);
ierr = DAGetLocalInfo(da2,&info2);CHKERRQ(ierr);
/* Get local vectors to hold ghosted parts of X */
ierr = DMCompositeGetLocalVectors(dm,&X1,&X2);CHKERRQ(ierr);
ierr = DMCompositeScatter(dm,X,X1,X2);CHKERRQ(ierr);
/* Access the arrays inside the subvectors of X */
ierr = DAVecGetArray(da1,X1,(void**)&x1);CHKERRQ(ierr);
ierr = DAVecGetArray(da2,X2,(void**)&x2);CHKERRQ(ierr);
/* Access the subvectors in F.
These are not ghosted so directly access the memory locations in F */
ierr = DMCompositeGetAccess(dm,F,&F1,&F2);CHKERRQ(ierr);
/* Access the arrays inside the subvectors of F */
ierr = DAVecGetArray(da1,F1,(void**)&f1);CHKERRQ(ierr);
ierr = DAVecGetArray(da2,F2,(void**)&f2);CHKERRQ(ierr);
/* Evaluate local user provided function */
ierr = FormFunctionLocal1(&info1,x1,x2,f1,(void**)user);CHKERRQ(ierr);
ierr = FormFunctionLocal2(&info2,x1,x2,f2,(void**)user);CHKERRQ(ierr);
ierr = DAVecRestoreArray(da1,F1,(void**)&f1);CHKERRQ(ierr);
ierr = DAVecRestoreArray(da2,F2,(void**)&f2);CHKERRQ(ierr);
ierr = DMCompositeRestoreAccess(dm,F,&F1,&F2);CHKERRQ(ierr);
ierr = DAVecRestoreArray(da1,X1,(void**)&x1);CHKERRQ(ierr);
ierr = DAVecRestoreArray(da2,X2,(void**)&x2);CHKERRQ(ierr);
ierr = DMCompositeRestoreLocalVectors(dm,&X1,&X2);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
|
