Lines Matching refs:c
2 c-----------------------------------------------------------------------
3 c
4 c this routine computes the body force term.
5 c
6 c currently this computes a swirl body with the axis alligned with
7 c the z-coordinate
8 c
9 c-----------------------------------------------------------------------
19 c
20 c This is the body force which will drive a swirl in a pipe flow
21 c
37 c$$$ nu=datmat(1,2,1)
38 c$$$ omag=datmat(3,5,1) ! frame rotation rate
39 c$$$ e1 = one/sqrt(3.0d0) ! axis of rotation
40 c$$$ e2 = e1
41 c$$$ e3 = e1
42 c$$$ om1=omag*e1
43 c$$$ om2=omag*e2
44 c$$$ om3=omag*e3
45 c$$$ w=0
46 c$$$
47 c$$$ do iel = 1, npro
48 c$$$
49 c$$$ x = xx(iel,1)
50 c$$$ y = xx(iel,2)
51 c$$$ z = xx(iel,3)
52 c$$$
53 c$$$c
54 c$$$c The following are MAPLE generated forcing functions for
55 c$$$c a contrived flow in a rotating reference frame.
56 c$$$c
57 c$$$
58 c$$$ t1 = x**2
59 c$$$ t2 = t1**2
60 c$$$ t5 = dexp(14.D0*x)
61 c$$$ t6 = t2*x*t5
62 c$$$ t7 = y**2
63 c$$$ t8 = t7**2
64 c$$$ t9 = t8*t7
65 c$$$ t12 = t2*t5
66 c$$$ t15 = nu*t1
67 c$$$ t17 = dexp(7.D0*x)
68 c$$$ t18 = t7*y
69 c$$$ t19 = t17*t18
70 c$$$ t22 = t1*x
71 c$$$ t23 = nu*t22
72 c$$$ t24 = t17*y
73 c$$$ t29 = t17*t7
74 c$$$ t34 = nu*x
75 c$$$ t43 = nu*t2
76 c$$$ t46 = -40.D0*t6*t9-6.D0*t12*t7+92.D0*t15*t19+132.D0*t23*t24+80.D0*
77 c$$$ #t6*t18-252.D0*t23*t29+168.D0*t23*t19+96.D0*t34*t29-64.D0*t34*t19+2
78 c$$$ #2.D0*t15*t24-138.D0*t15*t29+294.D0*t43*t29
79 c$$$ t50 = t2*t22*t5
80 c$$$ t57 = nu*t17
81 c$$$ t60 = t8*y
82 c$$$ t65 = t2**2
83 c$$$ t66 = t65*t5
84 c$$$ t73 = t22*t5
85 c$$$ t77 = t2*t1*t5
86 c$$$ t80 = -196.D0*t43*t19-96.D0*t50*t9-122.D0*t43*t24-32.D0*t34*t24-4.
87 c$$$ #D0*t57*y+36.D0*t12*t60+192.D0*t50*t18+98.D0*t66*t8+24.D0*t12*t18+2
88 c$$$ #8.D0*t66*t9-24.D0*t73*t60-336.D0*t77*t60
89 c$$$ t106 = -336.D0*t50*t8-12.D0*t12*t9+8.D0*t73*t9-140.D0*t6*t8+28.D0*
90 c$$$ #t73*t8+12.D0*t15*t17+12.D0*t43*t17+288.D0*t50*t60+112.D0*t77*t9-22
91 c$$$ #4.D0*t77*t18-56.D0*t66*t18+120.D0*t6*t60
92 c$$$ t131 = -8.D0*t57*t18-20.D0*t6*t7+56.D0*t77*t7-42.D0*t12*t8-24.D0*t
93 c$$$ #23*t17-48.D0*t50*t7+4.D0*t73*t7-16.D0*t73*t18+392.D0*t77*t8+14.D0*
94 c$$$ #t66*t7-84.D0*t66*t60+12.D0*t57*t7
95 c$$$ fx = t46+t80+t106+t131
96 c$$$
97 c$$$
98 c$$$ t1 = x**2
99 c$$$ t2 = t1**2
100 c$$$ t3 = nu*t2
101 c$$$ t5 = dexp(7.D0*x)
102 c$$$ t6 = t5*y
103 c$$$ t9 = y**2
104 c$$$ t10 = t9**2
105 c$$$ t11 = nu*t10
106 c$$$ t18 = t1*x
107 c$$$ t22 = nu*x
108 c$$$ t25 = nu*t18
109 c$$$ t32 = dexp(14.D0*x)
110 c$$$ t33 = t2*t32
111 c$$$ t39 = t2*t1*t32
112 c$$$ t40 = t10*t9
113 c$$$ t43 = t1*t32
114 c$$$ t46 = t9*y
115 c$$$ t47 = t10*t46
116 c$$$ t50 = -84.D0*t3*t6+343.D0*t11*t5*t2+66.D0*t11*t5*x-98.D0*t11*t5*t1
117 c$$$ #8-24.D0*t22*t6+120.D0*t25*t6-287.D0*t11*t5*t1-140.D0*t33*t10+14.D0
118 c$$$ #*t3*t5-56.D0*t39*t40-28.D0*t43*t40+8.D0*t43*t47
119 c$$$ t52 = t2*x*t32
120 c$$$ t55 = t18*t32
121 c$$$ t62 = t10*y
122 c$$$ t69 = nu*t5
123 c$$$ t80 = 120.D0*t52*t10-32.D0*t55*t47+56.D0*t33*t47+30.D0*t11*t5-216.
124 c$$$ #D0*t52*t62-144.D0*t55*t62+72.D0*t39*t62-60.D0*t69*t46+30.D0*t69*t9
125 c$$$ #-20.D0*t25*t5-20.D0*t43*t10+36.D0*t43*t62
126 c$$$ t86 = nu*t1
127 c$$$ t89 = t5*t9
128 c$$$ t92 = t5*t46
129 c$$$ t109 = 112.D0*t55*t40+80.D0*t55*t10-12.D0*t86*t6-275.D0*t86*t89+57
130 c$$$ #4.D0*t86*t92-218.D0*t25*t89+196.D0*t25*t92+90.D0*t22*t89-132.D0*t2
131 c$$$ #2*t92+427.D0*t3*t89+8.D0*t39*t46+4.D0*t43*t46
132 c$$$ t134 = 16.D0*t39*t47+28.D0*t33*t46-48.D0*t52*t47+252.D0*t33*t62+16
133 c$$$ #8.D0*t52*t40-24.D0*t52*t46-686.D0*t3*t92+2.D0*t86*t5-40.D0*t39*t10
134 c$$$ #-196.D0*t33*t40-16.D0*t55*t46+4.D0*t22*t5
135 c$$$ fy = t50+t80+t109+t134
136 c$$$
137 c$$$
138 c$$$ t3 = om3*x
139 c$$$ t5 = dexp(7.D0*x)
140 c$$$ t7 = x**2
141 c$$$ t12 = y**2
142 c$$$ t15 = (-1.D0+y)**2
143 c$$$ fxa = 2.D0*om2*w+2.D0*t3*t5*(2.D0+x-10.D0*t7+7.D0*t7*x)*t12*t15+om
144 c$$$ #2*(om1*y-om2*x)-om3*(t3-om1*z)
145 c$$$
146 c$$$ t1 = x**2
147 c$$$ t4 = (-1.D0+x)**2
148 c$$$ t7 = dexp(7.D0*x)
149 c$$$ t10 = y**2
150 c$$$ fya = 4.D0*om3*t1*t4*t7*y*(1.D0-3.D0*y+2.D0*t10)-2.D0*om1*w+om3*(o
151 c$$$ #m2*z-om3*y)-om1*(om1*y-om2*x)
152 c$$$
153 c$$$
154 c$$$ t3 = dexp(7.D0*x)
155 c$$$ t5 = x**2
156 c$$$ t10 = y**2
157 c$$$ t13 = (-1.D0+y)**2
158 c$$$ t19 = (-1.D0+x)**2
159 c$$$ fza = -2.D0*om1*x*t3*(2.D0+x-10.D0*t5+7.D0*t5*x)*t10*t13-4.D0*om2*
160 c$$$ #t5*t19*t3*y*(1.D0-3.D0*y+2.D0*t10)+om1*(om3*x-om1*z)-om2*(om2*z-om
161 c$$$ #3*y)
162 c$$$
163 c$$$ src(iel,1) = fx + fxa
164 c$$$ src(iel,2) = fy + fya
165 c$$$ src(iel,3) = fza
166 c$$$ enddo
174 c
175 c The following are MAPLE generated forcing functions for
176 c a Lid Driven cavity flow with an analytic solution
177 c
218 c-----------------------------------------------------------------------
219 c
220 c calculate the coefficients for the Spalart-Allmaras
221 c turbulence model and its jacobian
222 c
223 c input:
224 c Sclr : The scalar that is being transported in this pde.
225 c gradS : spatial gradient of Sclr
226 c rmu : kinematic viscosity
227 c
228 c output:
229 c rmu : diffusion for eddy viscosity equation
230 c srcR : residual terms for
231 c (srcR * turb)
232 c srcL : jacobian of srcR
233 c
234 c-----------------------------------------------------------------------
238 c coming in
246 c going out
249 c used locally
259 c Kay-Epsilon
264 c----------------------------------------------------------------------
265 c
266 c -- -- -- --
267 c | cb2 | 1 | |
268 c phi, + |u - -----phi, | phi, - -----|(nu+phi)phi, |
269 c t | i sigma i| i sigma| i|
270 c -- -- -- --,
271 c i
272 c
273 c ~ phi^2
274 c - cb1*S*phi + cw1*fw*-----
275 c d^2
276 c
277 c----------------------------------------------------------------------
278 c
280 c
281 c.... compute the global gradient of u
282 c
285 c
286 c du_i/dx_j
287 c
288 c i j indices match array where V is the velocity (u in our notes)
292 c
296 c
300 c a j u a i
301 c from our notes where we had N_{a,j} = dN_a/dx_j note that i is off by one because p was first in…
302 c
304 c
305 c.... magnitude of vorticity
306 c
316 c
317 c.... compute source and its jacobian
318 c
330 c
331 c Next 5 lines are DDES
332 c
340 c
341 c LHS
342 c
351 c DES97 dwall(e)=min(dmax,dwall(e))
411 c
412 c.... One source term has the form (beta_i)*(phi,_i). Since
413 c the convective term has (u_i)*(phi,_i), it is useful to treat
414 c beta_i as a "correction" to the velocity. In calculating the
415 c stabilization terms, the new "modified" velocity (u_i-beta_i) is
416 c then used in place of the pure velocity for stabilization terms,
417 c and the source term sneaks into the RHS and LHS. Here, the term
418 c is given by beta_i=(cb_2)*(phi,_i)/(sigma)
419 c
427 c.... compute the global gradient of u
428 c
431 c
432 c du_i/dx_j
433 c
434 c i j indices match array where V is the velocity (u in our notes)
438 c
442 c
446 c a j u a i
447 c from our notes where we had N_{a,j} = dN_a/dx_j note that i is off by one because p was first in…
448 c
465 c no source term in the LHS yet
479 c
480 c.... No source terms with the form (beta_i)*(phi,_i) for K or E
481 c