source: CIVL/examples/fortran/nek5000/core/cmt/bc.f

main
Last change on this file was ea777aa, checked in by Alex Wilton <awilton@…>, 3 years ago

Moved examples, include, build_default.properties, common.xml, and README out from dev.civl.com into the root of the repo.

git-svn-id: ​svn://vsl.cis.udel.edu/civl/trunk@5704 fb995dde-84ed-4084-dfe6-e5aef3e2452c
  • Property mode set to 100755
File size: 12.1 KB
Line 
1C> @file bc.f Boundary condition routines
2
3C> \ingroup bcond
4C> @{
5C> Determining rind state for Dirichlet boundary conditions
6 subroutine InviscidBC(wminus,wbc,nstate)
7!-------------------------------------------------------------------------------
8! JH091514 A fading copy of RFLU_ModAUSM.F90 from RocFlu
9!-------------------------------------------------------------------------------
10
11!#ifdef SPEC
12! USE ModSpecies, ONLY: t_spec_type
13!#endif
14 include 'SIZE'
15 include 'INPUT' ! do we need this?
16 include 'GEOM' ! for unx
17 include 'CMTDATA' ! do we need this without outflsub?
18! include 'TSTEP' ! for ifield?
19 include 'DG'
20
21! ==============================================================================
22! Arguments
23! ==============================================================================
24 integer nstate,nflux
25 real wminus(lx1*lz1,2*ldim,nelt,nstate),
26 > wbc(lx1*lz1,2*ldim,nelt,nstate)
27
28! ==============================================================================
29! Locals
30! ==============================================================================
31
32 integer e,f,fdim,i,k,nxz,nface,ifield
33 parameter (lfd=lxd*lzd)
34! JH111815 legacy rocflu names.
35!
36! nx,ny,nz : outward facing unit normal components
37! fs : face speed. zero until we have moving grid
38! jaco_c : fdim-D GLL grid Jacobian
39! nm : jaco_c, fine grid
40!
41! State on the interior (-, "left") side of the face
42! rl : density
43! ul,vl,wl : velocity
44! tl : temperature
45! al : sound speed
46! pl : pressure, then phi
47! cpl : rho*cp
48! State on the exterior (+, "right") side of the face
49! rr : density
50! ur,vr,wr : velocity
51! tr : temperature
52! ar : sound speed
53! pr : pressure
54! cpr : rho*cp
55
56 COMMON /SCRNS/ nx(lfd), ny(lfd), nz(lfd), rl(lfd), ul(lfd),
57 > vl(lfd), wl(lfd), pl(lfd), tl(lfd), al(lfd),
58 > cpl(lfd),rr(lfd), ur(lfd), vr(lfd), wr(lfd),
59 > pr(lfd),tr(lfd), ar(lfd),cpr(lfd),phl(lfd),fs(lfd),
60 > jaco_f(lfd),flx(lfd,toteq),jaco_c(lx1*lz1)
61 real nx, ny, nz, rl, ul, vl, wl, pl, tl, al, cpl, rr, ur, vr, wr,
62 > pr,tr, ar,cpr,phl,fs,jaco_f,flx,jaco_c
63
64! REAL vf(3)
65 real nTol
66 character*132 deathmessage
67 common /nekcb/ cb
68 character*3 cb
69
70 nTol = 1.0E-14
71
72 fdim=ldim-1
73 nface = 2*ldim
74 nxz = lx1*lz1
75 nxzd = lxd*lzd
76 ifield= 1 ! You need to figure out the best way of dealing with
77 ! this variable
78
79! if (outflsub)then
80! call maxMachnumber
81! endif
82 do e=1,nelt
83 do f=1,nface
84
85 cb=cbc(f,e,ifield)
86 if (cb.ne.'E '.and.cb.ne.'P ') then ! cbc bndy
87
88!-----------------------------------------------------------------------
89! compute flux for weakly-enforced boundary condition
90!-----------------------------------------------------------------------
91
92 do j=1,nstate
93 do i=1,nxz
94 if (abs(wbc(i,f,e,j)) .gt. ntol) then
95 write(6,*) nid,j,i,wbc(i,f,e,j),wminus(i,f,e,j),cb,
96 > nstate
97 write(deathmessage,*) 'GS hit a bndy,f,e=',f,e,'$'
98! Make sure you are not abusing this error handler
99 call exitti(deathmessage,f)
100 endif
101 enddo
102 enddo
103
104! JH060215 added SYM bc. Just use it as a slip wall hopefully.
105! JH021717 OK I just realized that this way doubles my userbc calls.
106! not sure if face loop and if-block for cb is a better way
107! to do it or not.
108 if (cb.eq.'v ' .or. cb .eq. 'V ') then
109 call inflow(nstate,f,e,wminus,wbc)
110 elseif (cb.eq.'O ') then
111 call outflow(nstate,f,e,wminus,wbc)
112 elseif (cb .eq. 'W ' .or. cb .eq.'I '.or.cb .eq.'SYM')then
113 call wallbc_inviscid(nstate,f,e,wminus,wbc)
114 endif
115
116 endif
117 enddo
118 enddo
119
120C> @}
121 return
122 end
123
124!-----------------------------------------------------------------------
125
126C> \ingroup bcond
127C> @{
128C> Determining IGU contribution to boundary flux. 0 for artificial
129C> viscosity, and strictly interior for physical viscosity.
130 subroutine bcflux(flux,agradu,qminus)
131! Need proper indexing and nekasgn & cmtasgn calls
132 include 'SIZE'
133 include 'INPUT'
134 include 'DG'
135! include 'NEKUSE'
136 include 'TSTEP' ! wait how do we know what ifield is?
137 integer e,eq,f
138 real flux (lx1*lz1,2*ldim,nelt,toteq)
139 real agradu(lx1*lz1,2*ldim,nelt,toteq)
140! real qminus(lx1*lz1,2*ldim,nelt,nqq) ! include CMTDATA?
141 real qminus(*) ! 'scuse me. comin' through
142 common /nekcb/ cb
143 character*3 cb
144
145 nfaces=2*ldim
146 nxz=lx1*lz1
147 ifield=1
148
149 do e=1,nelt
150 do f=1,nfaces
151 if (cbc(f, e, ifield).ne.'E '.and.
152 > cbc(f, e, ifield).ne.'P ') then ! cbc bndy
153 cb=cbc(f,e,ifield)
154 if (cb .eq. 'I ') then ! NEUMANN CONDITIONS GO HERE
155!-------------------------------------------------------------
156! JH112216 HARDCODING ADIABATIC WALL. DO SMARTER SOON
157 call rzero(flux(1,f,e,1),nxz)
158 do eq=2,ldim+1
159 call copy(flux(1,f,e,eq),agradu(1,f,e,eq),nxz)
160 enddo
161! METHOD "B", ADIABATIC NO-SLIP
162 call rzero(flux(1,f,e,toteq),nxz)
163! METHOD "A", ADIABATIC NO-SLIP augments with viscous work. triage below
164! because, predictably, NOW I need to computate AgradU on surfaces and I don't
165! have general code for that.
166 call a5adiabatic_wall(flux(1,1,1,toteq),f,e,agradu,
167 > qminus)
168! JH112216 HARDCODING ADIABATIC WALL. DO SMARTER SOON
169!-------------------------------------------------------------
170! cbu=cb
171! do eq=1,toteq
172! call userflux(flux(1,f,e,eq)) ! replace this with userbc
173! enddo
174 else ! if (cb .eq. 'SYM') then ! NEED PHYSICAL VISC TEST
175! JH031617 But this code block basically guarantees that artificial viscosity
176! does not contribute to viscous fluxes at boundaries.
177 do eq=1,toteq
178 call rzero(flux(1,f,e,eq),nxz)
179 enddo
180 endif
181 endif
182 enddo
183 enddo
184
185C> @}
186 return
187 end
188
189 subroutine a5adiabatic_wall(eflx,f,e,dU,wstate)
190 include 'SIZE'
191 include 'INPUT'
192 include 'GEOM' ! for UNX under ADIABATIC WALL METHOD "A"
193 include 'CMTDATA'
194 real eflx (lx1*lz1,2*ldim,nelt) ! better be zero on entry
195 real dU (lx1*lz1,2*ldim,nelt,toteq)
196 real wstate(lx1*lz1,2*ldim,nelt,nqq)
197 common /scrns/ flxscr(lx1*lz1)
198 real flxscr
199 integer e,f
200
201 nxz=lx1*lz1
202
203 call rzero(eflx(1,f,e),nxz)
204 call rzero(hface,nxz)
205
206 call a51dUadia(flxscr,f,e,dU,wstate)
207 call add2col2(eflx(1,f,e),flxsxcr,unx(1,1,f,e),nxz)
208 call a52dUadia(flxscr,f,e,dU,wstate)
209 call add2col2(eflx(1,f,e),flxsxcr,uny(1,1,f,e),nxz)
210 if (if3d) then
211 call a53dUadia(flxscr,f,e,dU,wstate)
212 call add2col2(eflx(1,f,e),flxsxcr,unz(1,1,f,e),nxz)
213 endif
214 return
215 end
216
217 subroutine a51dUadia(flux,f,ie,dU,wstate)
218! same as A51 for volume flux, but
219! 1. uses surface storage of quantities in wstate <-qminus (intent(in))
220! 2. SETS K=0. ADIABATIC WALLS HAVE VISCOUS HEATING, BUT DON'T CONDUCT
221 include 'SIZE'
222 include 'CMTDATA'
223 real wstate(lx1*lz1,2*ldim,nelt,nqq)
224 real dU (lx1*lz1,2*ldim,nelt,toteq,3)
225 real flux (lx1*ly1*lz1)
226 real K,E,kmcvmu,lambdamu
227 integer f
228 npt=lx1*lz1
229
230 do i=1,npt
231 dU1x=dU(i,f,ie,1,1)
232 dU2x=dU(i,f,ie,2,1)
233 dU3x=dU(i,f,ie,3,1)
234 dU4x=dU(i,f,ie,4,1)
235 dU5x=dU(i,f,ie,5,1)
236 dU1y=dU(i,f,ie,1,2)
237 dU2y=dU(i,f,ie,2,2)
238 dU3y=dU(i,f,ie,3,2)
239 dU4y=dU(i,f,ie,4,2)
240 dU5y=dU(i,f,ie,5,2)
241 dU1z=dU(i,f,ie,1,3)
242 dU2z=dU(i,f,ie,2,3)
243 dU3z=dU(i,f,ie,3,3)
244 dU4z=dU(i,f,ie,4,3)
245 dU5z=dU(i,f,ie,5,3)
246 rho =wstate(i,f,ie,irho)
247 cv =wstate(i,f,ie,icvf)/rho
248 lambda=wstate(i,f,ie,ilamf)
249 mu =wstate(i,f,ie,imuf)
250 K =0.0 ! ADIABATIC HARDCODING
251 u1 =wstate(i,f,ie,iux)
252 u2 =wstate(i,f,ie,iuy)
253 u3 =wstate(i,f,ie,iuz)
254 E =wstate(i,f,ie,iu5)/rho
255 lambdamu=lambda+mu
256 kmcvmu=K-cv*mu
257 flux(i)=
258 >(K*dU5x+cv*lambda*u1*dU4z-kmcvmu*u3*dU4x+cv*lambda*u1*dU3y
259 1 -kmcvmu*u2*dU3x+cv*mu*u3*dU2z+cv*mu*u2*dU2y+(cv*lambda-
260 2 K+2*cv*mu)*u1*dU2x-cv*lambdamu*u1*u3*dU1z-cv*lambdamu
261 3 *u1*u2*dU1y+(K*u3**2-cv*mu*u3**2+K*u2**2-cv*mu*u2**2-cv*la
262 4 mbda*u1**2+K*u1**2-2*cv*mu*u1**2-E*K)*dU1x)/(cv*rho)
263 enddo
264 return
265 end
266
267 subroutine a52dUadia(flux,f,ie,dU,wstate)
268! same as A52 for volume flux, but
269! 1. uses surface storage of quantities in wstate <-qminus (intent(in))
270! 2. SETS K=0. ADIABATIC WALLS HAVE VISCOUS HEATING, BUT DON'T CONDUCT
271 include 'SIZE'
272 include 'CMTDATA'
273 real wstate(lx1*lz1,2*ldim,nelt,nqq)
274 real dU (lx1*lz1,2*ldim,nelt,toteq,3)
275 real flux (lx1*ly1*lz1)
276 real K,E,kmcvmu,lambdamu
277 integer f
278 npt=lx1*lz1
279
280 do i=1,npt
281 dU1x=dU(i,f,ie,1,1)
282 dU2x=dU(i,f,ie,2,1)
283 dU3x=dU(i,f,ie,3,1)
284 dU4x=dU(i,f,ie,4,1)
285 dU5x=dU(i,f,ie,5,1)
286 dU1y=dU(i,f,ie,1,2)
287 dU2y=dU(i,f,ie,2,2)
288 dU3y=dU(i,f,ie,3,2)
289 dU4y=dU(i,f,ie,4,2)
290 dU5y=dU(i,f,ie,5,2)
291 dU1z=dU(i,f,ie,1,3)
292 dU2z=dU(i,f,ie,2,3)
293 dU3z=dU(i,f,ie,3,3)
294 dU4z=dU(i,f,ie,4,3)
295 dU5z=dU(i,f,ie,5,3)
296 rho =wstate(i,f,ie,irho)
297 cv =wstate(i,f,ie,icvf)/rho
298 lambda=wstate(i,f,ie,ilamf)
299 mu =wstate(i,f,ie,imuf)
300 K =0.0 ! ADIABATIC HARDCODING
301 u1 =wstate(i,f,ie,iux)
302 u2 =wstate(i,f,ie,iuy)
303 u3 =wstate(i,f,ie,iuz)
304 E =wstate(i,f,ie,iu5)/rho
305 lambdamu=lambda+mu
306 kmcvmu=K-cv*mu
307 flux(i)=
308 >(K*dU5y+cv*lambda*u2*dU4z-kmcvmu*u3*dU4y+cv*mu*u3*dU3z+(cv
309 1 *lambda-K+2*cv*mu)*u2*dU3y+cv*mu*u1*dU3x-kmcvmu*u1*dU2y+
310 2 cv*lambda*u2*dU2x-cv*lambdamu*u2*u3*dU1z+(K*u3**2-cv*mu
311 3 *u3**2-cv*lambda*u2**2+K*u2**2-2*cv*mu*u2**2+K*u1**2-cv*mu*
312 4 u1**2-E*K)*dU1y-cv*lambdamu*u1*u2*dU1x)/(cv*rho)
313 enddo
314 return
315 end
316
317 subroutine a53dUadia(flux,f,ie,dU,wstate)
318! same as A53 for volume flux, but
319! 1. uses surface storage of quantities in wstate <-qminus (intent(in))
320! 2. SETS K=0. ADIABATIC WALLS HAVE VISCOUS HEATING, BUT DON'T CONDUCT
321 include 'SIZE'
322 include 'CMTDATA'
323 real wstate(lx1*lz1,2*ldim,nelt,nqq)
324 real dU (lx1*lz1,2*ldim,nelt,toteq,3)
325 real flux (lx1*ly1*lz1)
326 real K,E,kmcvmu,lambdamu
327 integer f
328 npt=lx1*lz1
329
330 do i=1,npt
331 dU1x=dU(i,f,ie,1,1)
332 dU2x=dU(i,f,ie,2,1)
333 dU3x=dU(i,f,ie,3,1)
334 dU4x=dU(i,f,ie,4,1)
335 dU5x=dU(i,f,ie,5,1)
336 dU1y=dU(i,f,ie,1,2)
337 dU2y=dU(i,f,ie,2,2)
338 dU3y=dU(i,f,ie,3,2)
339 dU4y=dU(i,f,ie,4,2)
340 dU5y=dU(i,f,ie,5,2)
341 dU1z=dU(i,f,ie,1,3)
342 dU2z=dU(i,f,ie,2,3)
343 dU3z=dU(i,f,ie,3,3)
344 dU4z=dU(i,f,ie,4,3)
345 dU5z=dU(i,f,ie,5,3)
346 rho =wstate(i,f,ie,irho)
347 cv =wstate(i,f,ie,icvf)/rho
348 lambda=wstate(i,f,ie,ilamf)
349 mu =wstate(i,f,ie,imuf)
350 K =0.0 ! ADIABATIC HARDCODING
351 u1 =wstate(i,f,ie,iux)
352 u2 =wstate(i,f,ie,iuy)
353 u3 =wstate(i,f,ie,iuz)
354 E =wstate(i,f,ie,iu5)/rho
355 lambdamu=lambda+mu
356 kmcvmu=K-cv*mu
357 flux(i)=
358 >(K*(dU5z-E*dU1z)+cv*u3*(lambda*dU4z+2*mu*dU4z+lambda*dU3y+lambda
359 1 *dU2x)-K*u3*dU4z+cv*mu*u2*(dU4y+dU3z)+cv*mu*u1*(dU4x+dU2z)-
360 2 K*u2*dU3z-K*u1*dU2z-cv*(lambda+2*mu)*u3**2*dU1z+K*u3**2*dU1z+
361 3 K*u2**2*dU1z-cv*mu*u2**2*dU1z+K*u1**2*dU1z-cv*mu*u1**2*dU1z-c
362 4 _v*(lambda+mu)*u2*u3*dU1y-cv*(lambda+mu)*u1*u3*dU1x)/(cv*rho)
363 enddo
364 return
365 end
Note: See TracBrowser for help on using the repository browser.