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
|
#include "mp.h"
/*
Defines the u, v, omega physics for a given tempature
*/
#undef __FUNCT__
#define __FUNCT__ "FormInitialGuessLocal1"
PetscErrorCode FormInitialGuessLocal1(DALocalInfo *info,Field1 **x)
{
PetscInt i,j;
for (j=info->ys; j<info->ys+info->ym; j++) {
for (i=info->xs; i<info->xs+info->xm; i++) {
x[j][i].u = 0.0;
x[j][i].v = 0.0;
x[j][i].omega = 0.0;
}
}
return 0;
}
PetscLogEvent EVENT_FORMFUNCTIONLOCAL1;
#undef __FUNCT__
#define __FUNCT__ "FormFunctionLocal1"
/*
x2 contains given tempature field
*/
PetscErrorCode FormFunctionLocal1(DALocalInfo *info,Field1 **x,Field2 **x2,Field1 **f,void *ptr)
{
AppCtx *user = (AppCtx*)ptr;
PetscInt xints,xinte,yints,yinte,i,j;
PetscReal hx,hy,dhx,dhy,hxdhy,hydhx;
PetscReal grashof,prandtl,lid;
PetscScalar u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;
PetscErrorCode ierr;
PetscFunctionBegin;
ierr = PetscLogEventBegin(EVENT_FORMFUNCTIONLOCAL1,0,0,0,0);CHKERRQ(ierr);
grashof = user->grashof;
prandtl = user->prandtl;
lid = user->lidvelocity;
dhx = (PetscReal)(info->mx-1); dhy = (PetscReal)(info->my-1);
hx = 1.0/dhx; hy = 1.0/dhy;
hxdhy = hx*dhy; hydhx = hy*dhx;
xints = info->xs; xinte = info->xs+info->xm; yints = info->ys; yinte = info->ys+info->ym;
/* Test whether we are on the bottom edge of the global array */
if (yints == 0) {
j = 0;
yints = yints + 1;
/* bottom edge */
for (i=info->xs; i<info->xs+info->xm; i++) {
f[j][i].u = x[j][i].u;
f[j][i].v = x[j][i].v;
f[j][i].omega = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
}
}
/* Test whether we are on the top edge of the global array */
if (yinte == info->my) {
j = info->my - 1;
yinte = yinte - 1;
/* top edge */
for (i=info->xs; i<info->xs+info->xm; i++) {
f[j][i].u = x[j][i].u - lid;
f[j][i].v = x[j][i].v;
f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
}
}
/* Test whether we are on the left edge of the global array */
if (xints == 0) {
i = 0;
xints = xints + 1;
/* left edge */
for (j=info->ys; j<info->ys+info->ym; j++) {
f[j][i].u = x[j][i].u;
f[j][i].v = x[j][i].v;
f[j][i].omega = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
}
}
/* Test whether we are on the right edge of the global array */
if (xinte == info->mx) {
i = info->mx - 1;
xinte = xinte - 1;
/* right edge */
for (j=info->ys; j<info->ys+info->ym; j++) {
f[j][i].u = x[j][i].u;
f[j][i].v = x[j][i].v;
f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
}
}
/* Compute over the interior points */
for (j=yints; j<yinte; j++) {
for (i=xints; i<xinte; i++) {
/* convective coefficients for upwinding */
vx = x[j][i].u; avx = PetscAbsScalar(vx);
vxp = .5*(vx+avx); vxm = .5*(vx-avx);
vy = x[j][i].v; avy = PetscAbsScalar(vy);
vyp = .5*(vy+avy); vym = .5*(vy-avy);
/* U velocity */
u = x[j][i].u;
uxx = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
uyy = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
f[j][i].u = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;
/* V velocity */
u = x[j][i].v;
uxx = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
uyy = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
f[j][i].v = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;
/* Omega */
u = x[j][i].omega;
uxx = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
uyy = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
f[j][i].omega = uxx + uyy + (vxp*(u - x[j][i-1].omega) + vxm*(x[j][i+1].omega - u))*hy + (vyp*(u - x[j-1][i].omega) + vym*(x[j+1][i].omega - u))*hx;
if (x2) {
f[j][i].omega += - .5 * grashof * (x2[j][i+1].temp - x2[j][i-1].temp) * hy;
}
}
}
ierr = PetscLogEventEnd(EVENT_FORMFUNCTIONLOCAL1,0,0,0,0);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
|
