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
|
<center><a href="ff1.c">Actual source code: ff1.c</a></center><br>
<html>
<head>
<title></title>
<meta name="generator" content="c2html 0.9.5">
<meta name="date" content="2010-04-08T19:58:01+00:00">
</head>
<body bgcolor="#FFFFFF">
<pre width="80">
<a name="line2"> 2: </a> #include <A href="mp.h.html">mp.h</A>
<a name="line3"> 3: </a><font color="#B22222">/*</font>
<a name="line4"> 4: </a><font color="#B22222"> Defines the u, v, omega physics for a given tempature</font>
<a name="line5"> 5: </a><font color="#B22222">*/</font>
<a name="line8"> 8: </a><strong><font color="#4169E1"><a name="FormInitialGuessLocal1"></a><A href="tutorials/../../docs/manualpages/Sys/PetscErrorCode.html#PetscErrorCode">PetscErrorCode</A> FormInitialGuessLocal1(<A href="tutorials/../../docs/manualpages/DA/DALocalInfo.html#DALocalInfo">DALocalInfo</A> *info,Field1 **x)</font></strong>
<a name="line9"> 9: </a>{
<a name="line10"> 10: </a> <A href="tutorials/../../docs/manualpages/Sys/PetscInt.html#PetscInt">PetscInt</A> i,j;
<a name="line11"> 11: </a>
<a name="line12"> 12: </a> <font color="#4169E1">for</font> (j=info->ys; j<info->ys+info->ym; j++) {
<a name="line13"> 13: </a> <font color="#4169E1">for</font> (i=info->xs; i<info->xs+info->xm; i++) {
<a name="line14"> 14: </a> x[j][i].u = 0.0;
<a name="line15"> 15: </a> x[j][i].v = 0.0;
<a name="line16"> 16: </a> x[j][i].omega = 0.0;
<a name="line17"> 17: </a> }
<a name="line18"> 18: </a> }
<a name="line19"> 19: </a> <font color="#4169E1">return</font> 0;
<a name="line20"> 20: </a>}
<a name="line22"> 22: </a><A href="tutorials/../../docs/manualpages/Sys/PetscLogEvent.html#PetscLogEvent">PetscLogEvent</A> EVENT_FORMFUNCTIONLOCAL1;
<a name="line26"> 26: </a><font color="#B22222">/* </font>
<a name="line27"> 27: </a><font color="#B22222"> x2 contains given tempature field</font>
<a name="line28"> 28: </a><font color="#B22222">*/</font>
<a name="line29"> 29: </a><strong><font color="#4169E1"><a name="FormFunctionLocal1"></a><A href="tutorials/../../docs/manualpages/Sys/PetscErrorCode.html#PetscErrorCode">PetscErrorCode</A> FormFunctionLocal1(<A href="tutorials/../../docs/manualpages/DA/DALocalInfo.html#DALocalInfo">DALocalInfo</A> *info,Field1 **x,Field2 **x2,Field1 **f,void *ptr)</font></strong>
<a name="line30"> 30: </a> {
<a name="line31"> 31: </a> AppCtx *user = (AppCtx*)ptr;
<a name="line32"> 32: </a> <A href="tutorials/../../docs/manualpages/Sys/PetscInt.html#PetscInt">PetscInt</A> xints,xinte,yints,yinte,i,j;
<a name="line33"> 33: </a> <A href="tutorials/../../docs/manualpages/Sys/PetscReal.html#PetscReal">PetscReal</A> hx,hy,dhx,dhy,hxdhy,hydhx;
<a name="line34"> 34: </a> <A href="tutorials/../../docs/manualpages/Sys/PetscReal.html#PetscReal">PetscReal</A> grashof,prandtl,lid;
<a name="line35"> 35: </a> <A href="tutorials/../../docs/manualpages/Sys/PetscScalar.html#PetscScalar">PetscScalar</A> u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;
<a name="line39"> 39: </a> <A href="tutorials/../../docs/manualpages/Profiling/PetscLogEventBegin.html#PetscLogEventBegin">PetscLogEventBegin</A>(EVENT_FORMFUNCTIONLOCAL1,0,0,0,0);
<a name="line40"> 40: </a> grashof = user->grashof;
<a name="line41"> 41: </a> prandtl = user->prandtl;
<a name="line42"> 42: </a> lid = user->lidvelocity;
<a name="line44"> 44: </a> dhx = (<A href="tutorials/../../docs/manualpages/Sys/PetscReal.html#PetscReal">PetscReal</A>)(info->mx-1); dhy = (<A href="tutorials/../../docs/manualpages/Sys/PetscReal.html#PetscReal">PetscReal</A>)(info->my-1);
<a name="line45"> 45: </a> hx = 1.0/dhx; hy = 1.0/dhy;
<a name="line46"> 46: </a> hxdhy = hx*dhy; hydhx = hy*dhx;
<a name="line48"> 48: </a> xints = info->xs; xinte = info->xs+info->xm; yints = info->ys; yinte = info->ys+info->ym;
<a name="line50"> 50: </a> <font color="#B22222">/* Test whether we are on the bottom edge of the global array */</font>
<a name="line51"> 51: </a> <font color="#4169E1">if</font> (yints == 0) {
<a name="line52"> 52: </a> j = 0;
<a name="line53"> 53: </a> yints = yints + 1;
<a name="line54"> 54: </a> <font color="#B22222">/* bottom edge */</font>
<a name="line55"> 55: </a> <font color="#4169E1">for</font> (i=info->xs; i<info->xs+info->xm; i++) {
<a name="line56"> 56: </a> f[j][i].u = x[j][i].u;
<a name="line57"> 57: </a> f[j][i].v = x[j][i].v;
<a name="line58"> 58: </a> f[j][i].omega = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy;
<a name="line59"> 59: </a> }
<a name="line60"> 60: </a> }
<a name="line62"> 62: </a> <font color="#B22222">/* Test whether we are on the top edge of the global array */</font>
<a name="line63"> 63: </a> <font color="#4169E1">if</font> (yinte == info->my) {
<a name="line64"> 64: </a> j = info->my - 1;
<a name="line65"> 65: </a> yinte = yinte - 1;
<a name="line66"> 66: </a> <font color="#B22222">/* top edge */</font>
<a name="line67"> 67: </a> <font color="#4169E1">for</font> (i=info->xs; i<info->xs+info->xm; i++) {
<a name="line68"> 68: </a> f[j][i].u = x[j][i].u - lid;
<a name="line69"> 69: </a> f[j][i].v = x[j][i].v;
<a name="line70"> 70: </a> f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy;
<a name="line71"> 71: </a> }
<a name="line72"> 72: </a> }
<a name="line74"> 74: </a> <font color="#B22222">/* Test whether we are on the left edge of the global array */</font>
<a name="line75"> 75: </a> <font color="#4169E1">if</font> (xints == 0) {
<a name="line76"> 76: </a> i = 0;
<a name="line77"> 77: </a> xints = xints + 1;
<a name="line78"> 78: </a> <font color="#B22222">/* left edge */</font>
<a name="line79"> 79: </a> <font color="#4169E1">for</font> (j=info->ys; j<info->ys+info->ym; j++) {
<a name="line80"> 80: </a> f[j][i].u = x[j][i].u;
<a name="line81"> 81: </a> f[j][i].v = x[j][i].v;
<a name="line82"> 82: </a> f[j][i].omega = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx;
<a name="line83"> 83: </a> }
<a name="line84"> 84: </a> }
<a name="line86"> 86: </a> <font color="#B22222">/* Test whether we are on the right edge of the global array */</font>
<a name="line87"> 87: </a> <font color="#4169E1">if</font> (xinte == info->mx) {
<a name="line88"> 88: </a> i = info->mx - 1;
<a name="line89"> 89: </a> xinte = xinte - 1;
<a name="line90"> 90: </a> <font color="#B22222">/* right edge */</font>
<a name="line91"> 91: </a> <font color="#4169E1">for</font> (j=info->ys; j<info->ys+info->ym; j++) {
<a name="line92"> 92: </a> f[j][i].u = x[j][i].u;
<a name="line93"> 93: </a> f[j][i].v = x[j][i].v;
<a name="line94"> 94: </a> f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx;
<a name="line95"> 95: </a> }
<a name="line96"> 96: </a> }
<a name="line98"> 98: </a> <font color="#B22222">/* Compute over the interior points */</font>
<a name="line99"> 99: </a> <font color="#4169E1">for</font> (j=yints; j<yinte; j++) {
<a name="line100">100: </a> <font color="#4169E1">for</font> (i=xints; i<xinte; i++) {
<a name="line101">101: </a> <font color="#B22222">/* convective coefficients for upwinding */</font>
<a name="line102">102: </a> vx = x[j][i].u; avx = PetscAbsScalar(vx);
<a name="line103">103: </a> vxp = .5*(vx+avx); vxm = .5*(vx-avx);
<a name="line104">104: </a> vy = x[j][i].v; avy = PetscAbsScalar(vy);
<a name="line105">105: </a> vyp = .5*(vy+avy); vym = .5*(vy-avy);
<a name="line107">107: </a> <font color="#B22222">/* U velocity */</font>
<a name="line108">108: </a> u = x[j][i].u;
<a name="line109">109: </a> uxx = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
<a name="line110">110: </a> uyy = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
<a name="line111">111: </a> f[j][i].u = uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx;
<a name="line113">113: </a> <font color="#B22222">/* V velocity */</font>
<a name="line114">114: </a> u = x[j][i].v;
<a name="line115">115: </a> uxx = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
<a name="line116">116: </a> uyy = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
<a name="line117">117: </a> f[j][i].v = uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy;
<a name="line119">119: </a> <font color="#B22222">/* Omega */</font>
<a name="line120">120: </a> u = x[j][i].omega;
<a name="line121">121: </a> uxx = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
<a name="line122">122: </a> uyy = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
<a name="line123">123: </a> 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;
<a name="line124">124: </a> <font color="#4169E1">if</font> (x2) {
<a name="line125">125: </a> f[j][i].omega += - .5 * grashof * (x2[j][i+1].temp - x2[j][i-1].temp) * hy;
<a name="line126">126: </a> }
<a name="line127">127: </a> }
<a name="line128">128: </a> }
<a name="line129">129: </a> <A href="tutorials/../../docs/manualpages/Profiling/PetscLogEventEnd.html#PetscLogEventEnd">PetscLogEventEnd</A>(EVENT_FORMFUNCTIONLOCAL1,0,0,0,0);
<a name="line130">130: </a> <font color="#4169E1">return</font>(0);
<a name="line131">131: </a>}
</pre>
</body>
</html>
|
