From: "Lee, Bo-sung" Newsgroups: comp.parallel.mpi Subject: Re: [Q] Out-of-order calls in MPI vs. PVM Date: Sun, 27 Jun 1999 16:43:18 +0900 Organization: ETRI Supercomputer Center / Seoul National Univ. Sender: bslee6@210.216.149.243 Message-Id: <3775D615.1B9018BC@pvmcube4.snu.ac.kr> References: <7l0u34$q9g$1@mozo.cc.purdue.edu> Reply-To: bslee@pvmcube4.snu.ac.kr Mime-Version: 1.0 Content-Type: multipart/mixed; boundary="------------BB2C3A587E7F9B68300F0821" Xref: ukc comp.parallel.mpi:5251 This is a multi-part message in MIME format. --------------BB2C3A587E7F9B68300F0821 Content-Type: text/plain; charset=EUC-KR Content-Transfer-Encoding: 7bit "Brian R. Jones" wrote: > Hello! > > [Short version at the end] > > Background: > > I am in the process of converting a Fortran 77 3D CFD code from > PVM to MPI. The ability to compile with PVM or MPI is a requirement, > thus I can not (extensively) modify the structure of the code. I > realize that this limits the extent to which I can exploit the > added capability of MPI. Eventually, PVM will be dropped and a > more effective MPI rewrite will occur. > > That said, I currently have a direct translation completed (meaning > pvmfsend mapped to MPI_Send, etc.) and even managed to sidestep the > PVM spawning process in MPI. I am having problems, however, with > the data exchange between "blocks" of the 3D computational grid. > In the PVM version of the code, the process for each node looks > like this: > > 1) perform computation on current block > 2) pass data (via pvmfsend/MPI_Send) at interface between > current block and neighboring blocks to corresponding > nodes. > 3) wait for interface data from neighboring nodes/blocks > (via pvmfrecv/MPI_Recv) > > Note that there is no synchronization here.... if all the nodes had > an equal computational burden, they would all send their data at the > same time and then all receive data at the same time. > > While I have had a bit of exposure to MPI, I'm new to PVM. From what I > can tell, however, pvmfsend and pvmfrecv equate to MPI_Send and MPI_Recv > (i.e. both are blocking). Dismissing that possibility, I'm now working > on why the program hangs. I'm guessing that the PVM daemon is handling > the out-of-order communication that I will need to explicitly address > in MPI. > > The short version: > > A node has a message pending for receipt, but can't post the receive > until after sending off another message. The buffering supplied by > the PVM daemon will allow this behavior, while the MPI code must be > structured to handle the MPI commands in the order that they occur. > Am I interpreting this correctly? Is the simple (but ugly) solution > to add MPI_Probe commands before each send to test for a pending > receive? Two years ago, I translated my pvm cfd code to mpi cfd code like you. First I changed pvm routine to mpi routine line-by-line. After succeeding in running in MPI code, I rewrited my MPI code using derived data type to simplify the stride data packing in PVM. I had no problem like you. I think your problem results from MPI_SEND-RECV mode. In PVM, pvmfsend always returns whether the message is arrived to receiver or not. i.e. in pvm_send, message always send to send buffer not to receiver and never waits for receive acknoledgement. But in mpi, MPI_BSEND sends message to message buffer and returns after messages are sent to buffer. and MPI_SSEND always waits for receiver's acknowledgements. And MPI_SEND is either MPI_BSEND or MPI_SSEND in according to the mpi implementations in your system. If your code always hang up in send-recv routines, i think your MPI_SEND routine operates as synchoronous mode(MPI_SSEND). So I think you need to use MPI_BSEND rather MPI_SEND to run like pvm code. I attached my two version of codes below this message. Two versions are nearly the same. you can compare the send-recv routines easily. And I didn't use group function in pvm code because the group fuctions in pvm are not well operated. Sorry for my poor English :) If i can explain by Korean, it would be more helpful to you. if you want full version of the code, send mail to me. -- Lee, Bo-sung Ph.D Candidate in Aerospace Eng. Seoul National Univ. Invited Reseacher in ETRI/Supercomputing Center mailto:bslee@pvmcube4.snu.ac.kr mailto:bslee@bslee.hpcnet.ne.kr http://www.shinbiro.com/~bslee6 --------------BB2C3A587E7F9B68300F0821 Content-Type: text/plain; charset=EUC-KR; name="ns2d-pvm.f" Content-Transfer-Encoding: 7bit Content-Disposition: inline; filename="ns2d-pvm.f" program NSSPMD2D include 'nsspmd2d.inc' include 'fpvm3.h' integer task_ids(maxpe,maxpe) INTEGERX = INTEGER4 IREALX = REAL4 call pvmfmytid(mytid) call pvmfjoingroup('PVMALL',mype) call pvmfsetopt(PvmRoute,PvmRouteDirect,ioldval) call pvmfgsize('PVMALL',npes) if (mype.eq.0) then open(12,file='init.inp') open(13,file='grid.grd') read(12,500) read(12,*) alpa, re, amach, tinf read(12,500) read(12,*) nx, ny, tol, itmax, isubmax read(12,500) read(12,*) ifread, delt, intwrt, cfl read(12,500) read(12,*) cappa, eee, isweep, init read(12,500) read(12,*) iminbc, imaxbc, jminbc, jmaxbc read(12,500) read(12,*) i_proc, j_proc nproc = i_proc*j_proc do i = 1, i_proc do j = 1, j_proc call pvmfgettid('PVMALL',(i-1)*j_proc+j-1,task_ids(i,j)) enddo enddo open(21,file='result.dat') open(22,file='wnstep.dat') open(23,file='wn1step.dat') open(30,file='error.dat') call GRID if (ifread.lt.0) then ptime = 0.0 itbegin = 1 else read(22,*) ptime, itbegin, errbegin do i = 1, nx+1 do j = 1, ny+1 read(22,*) qq(1,i,j),qq(2,i,j),qq(3,i,j), . qq(4,i,j),qq(5,i,j),qq(6,i,j) read(23,*) qm(1,i,j),qm(2,i,j),qm(3,i,j), . qm(4,i,j),qm(5,i,j),qm(6,i,j) enddo enddo endif do i = 1, i_proc do j = 1, j_proc if (i.eq.1) then if (iminbc.eq.PERIODIC) then imtask = task_ids(i_proc,j) else imtask = iminbc endif else imtask = task_ids(i-1,j) endif if (i.eq.i_proc) then if (imaxbc.eq.PERIODIC) then iptask = task_ids(1,j) else iptask = imaxbc endif else iptask = task_ids(i+1,j) endif if (j.eq.1) then jmtask = jminbc else jmtask = task_ids(i,j-1) endif if (j.eq.j_proc) then jptask = jmaxbc else jptask = task_ids(i,j+1) endif i_cell = (nx-1)/i_proc j_cell = (ny-1)/j_proc i_remain = (nx-1)-(nx-1)/i_proc*i_proc j_remain = (ny-1)-(ny-1)/j_proc*j_proc if (i.le.i_remain) i_cell = i_cell+1 if (j.le.j_remain) j_cell = j_cell+1 i_end = i*i_cell+1 j_end = j*j_cell+1 if (i.gt.i_remain) i_end = i_end+i_remain if (j.gt.j_remain) j_end = j_end+j_remain i_begin = i_end-i_cell j_begin = j_end-j_cell nfx_slave = i_cell+1 nfy_slave = j_cell+1 call pvmfinitsend(PvmDataRaw,ibufid) call pvmfpack(INTEGERX,imtask,1,1,info) call pvmfpack(INTEGERX,iptask,1,1,info) call pvmfpack(INTEGERX,jmtask,1,1,info) call pvmfpack(INTEGERX,jptask,1,1,info) call pvmfpack(INTEGERX,nfx_slave,1,1,info) call pvmfpack(INTEGERX,nfy_slave,1,1,info) call pvmfpack(INTEGERX,itbegin,1,1,info) call pvmfpack(INTEGERX,itmax,1,1,info) call pvmfpack(INTEGERX,intwrt,1,1,info) call pvmfpack(INTEGERX,isweep,1,1,info) call pvmfpack(INTEGERX,ifread,1,1,info) call pvmfpack(INTEGERX,init,1,1,info) call pvmfpack(IREALX,alpa,1,1,info) call pvmfpack(IREALX,re,1,1,info) call pvmfpack(IREALX,amach,1,1,info) call pvmfpack(IREALX,tinf,1,1,info) call pvmfpack(IREALX,cfl,1,1,info) call pvmfpack(IREALX,delt,1,1,info) call pvmfpack(IREALX,cappa,1,1,info) call pvmfpack(IREALX,eee,1,1,info) do jj = j_begin+1,j_end+1 call pvmfpack(IREALX,gx(i_begin+1,jj), . i_end-i_begin+1,1,info) call pvmfpack(IREALX,gy(i_begin+1,jj), . i_end-i_begin+1,1,info) enddo if (ifread.ge.0) then do jj = j_begin,j_end+1 call pvmfpack(IREALX,qq(1,i_begin,jj), . (i_end-i_begin+2)*6,1,info) call pvmfpack(IREALX,qm(1,i_begin,jj), . (i_end-i_begin+2)*6,1,info) enddo endif do jj = j_begin+1,j_end call pvmfpack(INTEGERX,nb(i_begin+1,jj), . i_end-i_begin,1,info) call pvmfpack(IREALX,yyn(i_begin+1,jj), . i_end-i_begin,1,info) enddo call pvmfsend(task_ids(i,j),666,info) enddo enddo endif call pvmfgettid('PVMALL',0,master_id) call pvmfrecv(master_id,666,ibufid) call pvmfunpack(INTEGERX,imtask,1,1,info) call pvmfunpack(INTEGERX,iptask,1,1,info) call pvmfunpack(INTEGERX,jmtask,1,1,info) call pvmfunpack(INTEGERX,jptask,1,1,info) call pvmfunpack(INTEGERX,nfx,1,1,info) call pvmfunpack(INTEGERX,nfy,1,1,info) call pvmfunpack(INTEGERX,itbegin,1,1,info) call pvmfunpack(INTEGERX,itmax,1,1,info) call pvmfunpack(INTEGERX,intwrt,1,1,info) call pvmfunpack(INTEGERX,isweep,1,1,info) call pvmfunpack(INTEGERX,ifread,1,1,info) call pvmfunpack(INTEGERX,init,1,1,info) call pvmfunpack(IREALX,alpa,1,1,info) call pvmfunpack(IREALX,re,1,1,info) call pvmfunpack(IREALX,amach,1,1,info) call pvmfunpack(IREALX,tinf,1,1,info) call pvmfunpack(IREALX,cfl,1,1,info) call pvmfunpack(IREALX,delt,1,1,info) call pvmfunpack(IREALX,cappa,1,1,info) call pvmfunpack(IREALX,eee,1,1,info) do j = 2, nfy+1 call pvmfunpack(IREALX,x(2,j),nfx,1,info) call pvmfunpack(IREALX,y(2,j),nfx,1,info) enddo pi = 4.0*atan(1.0) alpa = alpa*pi/180.0 rho0 = 1.0 u0 = amach*cos(alpa) v0 = amach*sin(alpa) e0 = 1./gamma/gam1+0.5*amach**2 p0 = 1.0/gamma c0 = sqrt(gamma*p0/rho0) ra = re/amach tcoe = 110.4/tinf rk0 = 1.e-7 rw0 = 1.0 tgamma1 = beta1/betaA-sigmao1*0.41**2/sqrt(betaA) tgamma2 = beta2/betaA-sigmao2*0.41**2/sqrt(betaA) nfx_start = 1 nfx_end = nfx nfy_start = 1 nfy_end = nfy if (imtask.lt.0) nfx_start = 2 if (iptask.lt.0) nfx_end = nfx-1 if (jmtask.lt.0) nfy_start = 2 if (jptask.lt.0) nfy_end = nfy-1 do j = 1, nfy+1 do i = 1, nfx+1 ib(i,j) = 1 zac(i,j) = 1.0 ff(i,j) = 1.0 enddo enddo if (ifread.ge.0) then do j = 1, nfy+1 call pvmfunpack(IREALX,w_n(1,1,j),(nfx+1)*6,1,info) call pvmfunpack(IREALX,w_nm(1,1,j),(nfx+1)*6,1,info) do i = 1, nfx+1 w(1,i,j) = w_n(1,i,j) w(2,i,j) = w_n(2,i,j) w(3,i,j) = w_n(3,i,j) w(4,i,j) = w_n(4,i,j) w(5,i,j) = w_n(5,i,j)/w(1,i,j) w(6,i,j) = w_n(6,i,j)/w(1,i,j) p(i,j) = gam1*(w(4,i,j) . -0.5*(w(2,i,j)**2+w(3,i,j)**2)/w(1,i,j)) c(i,j) = sqrt(gamma*p(i,j)/w(1,i,j)) enddo enddo else if (init.le.0) then do j = 1, nfy+1 do i = 1, nfx+1 w(2,i,j) = 0 w(3,i,j) = 0 enddo enddo else if (init.eq.1) then do j = 1, nfy+1 do i = 1, nfx+1 w(2,i,j) = rho0*u0 w(3,i,j) = rho0*v0 enddo enddo else if (init.eq.2) then do j = 1, nfy+1 do i = 1, nfx+1 if (y(i,j).ge.0) then esp = -0.01 else esp = +0.01 endif w(2,i,j) = rho0*u0*(1+esp) w(3,i,j) = rho0*v0 enddo enddo endif do j = 1, nfy+1 do i = 1, nfx+1 w(1,i,j) = rho0 w(4,i,j) = rho0*e0 w(5,i,j) = rk0 w(6,i,j) = rw0 p(i,j) = p0 c(i,j) = c0 w_n(1,i,j) = w(1,i,j) w_n(2,i,j) = w(2,i,j) w_n(3,i,j) = w(3,i,j) w_n(4,i,j) = w(4,i,j) w_n(5,i,j) = w(5,i,j)*w(1,i,j) w_n(6,i,j) = w(6,i,j)*w(1,i,j) w_nm(1,i,j) = w(1,i,j) w_nm(2,i,j) = w(2,i,j) w_nm(3,i,j) = w(3,i,j) w_nm(4,i,j) = w(4,i,j) w_nm(5,i,j) = w(5,i,j)*w(1,i,j) w_nm(6,i,j) = w(6,i,j)*w(1,i,j) enddo enddo endif do j = 2, nfy call pvmfunpack(INTEGERX,ib(2,j),nfx-1,1,info) call pvmfunpack(IREALX,yn(2,j),nfx-1,1,info) enddo call HOLE_CHECK call METRIC call METDUM call BC call RENEW call F1CALC call VISBC i_total = 0 iconv = 0 do it = itbegin, itmax itsub = 0 ptime = ptime + delt call STEP call NODE call PHYSICS call SOLVE call UPGRADE call BC call RENEW call F1CALC call VISBC do while (iconv.eq.0) itsub = itsub + 1 call STEP call NODE call PHYSICS call SOLVE call UPGRADE call BC call RENEW call F1CALC call VISBC err = 0.0 do j = 2, nfy do i = 2, nfx err = err + dw(1,i,j)**2 enddo enddo call pvmfinitsend(PvmDataRaw,ibufid) call pvmfpack(IREALX,err,1,1,info) call pvmfpack(IREALX,dtmin,1,1,info) call pvmfpack(IREALX,dtmax,1,1,info) call pvmfsend(master_id,777,info) if (mype.eq.0) then i_total = i_total + 1 errtotal = 0.0 dtmin = 100 dtmax = 0 do ii = 1, i_proc do jj = 1, j_proc call pvmfrecv(task_ids(ii,jj),777,ibufid) call pvmfunpack(IREALX,err,1,1,info) call pvmfunpack(IREALX,dtmi,1,1,info) call pvmfunpack(IREALX,dtma,1,1,info) errtotal = errtotal + err dtmin = min(dtmin, dtmi) dtmax = max(dtmax, dtma) enddo enddo errtotal = log10(sqrt(errtotal/((nx-1)*(ny-1)))) if (it.eq.1.and.i_total.eq.1) errbegin = errtotal errtotal = errtotal - errbegin write(0,400) i_total, errtotal,dtmin,dtmax write(30,400) i_total, errtotal,dtmin,dtmax if (errtotal.lt.tol.or.itsub.eq.isubmax) then iconv = 1 else iconv = 0 endif call pvmfinitsend(PvmDataRaw,ibufid) call pvmfpack(INTEGERX,iconv,1,1,info) do ii = 1, i_proc do jj = 1, j_proc call pvmfsend(task_ids(ii,jj),888,info) enddo enddo endif call pvmfrecv(master_id,888,ibufid) call pvmfunpack(INTEGERX,iconv,1,1,info) enddo do j = 1, nfy+1 do i = 1, nfx+1 do n = 1, 4 w_nm(n,i,j) = w_n(n,i,j) w_n(n,i,j) = w(n,i,j) enddo w_nm(5,i,j) = w_n(5,i,j) w_nm(6,i,j) = w_n(6,i,j) w_n(5,i,j) = w(5,i,j)*w(1,i,j) w_n(6,i,j) = w(6,i,j)*w(1,i,j) enddo enddo call FORCE(rLp,rDp,rLf,rDf) call pvmfinitsend(PvmDataRaw,ibufid) call pvmfpack(IREALX,rLp,1,1,info) call pvmfpack(IREALX,rDp,1,1,info) call pvmfpack(IREALX,rLf,1,1,info) call pvmfpack(IREALX,rDf,1,1,info) call pvmfsend(master_id,7777,info) if (myPE.eq.0) then rLp = 0.0 rDp = 0.0 rLf = 0.0 rDf = 0.0 do ii = 1, i_proc do jj = 1, j_proc call pvmfrecv(task_ids(ii,jj),7777,ibufid) call pvmfunpack(IREALX,rLpp,1,1,info) call pvmfunpack(IREALX,rDpp,1,1,info) call pvmfunpack(IREALX,rLff,1,1,info) call pvmfunpack(IREALX,rDff,1,1,info) rLp = rLp + rLpp rDp = rDp + rDpp rLf = rLf + rLff rDf = rDf + rDff enddo enddo Clp = rLp/(0.5*rho0*amach**2) Cdp = rDp/(0.5*rho0*amach**2) Cl = (rLp+rLf)/(0.5*rho0*amach**2) Cd = (rDp+rDf)/(0.5*rho0*amach**2) write(21,340) ptime, Clp, Cdp, Cl, Cd, itsub call flush(21) endif if (myPE.ne.0.and.mod(it-itbegin+1,intwrt).eq.0) then call pvmfinitsend(PvmDataRaw,ibufid) do j = 1, nfy+1 call pvmfpack(IREALX,w_n(1,1,j),(nfx+1)*6,1,info) call pvmfpack(IREALX,w_nm(1,1,j),(nfx+1)*6,1,info) call pvmfpack(IREALX,vis(1,j),nfx+1,1,info) enddo call pvmfsend(master_id,999,info) endif if (myPE.eq.0.and.mod(it-itbegin+1,intwrt).eq.0) then do ii = 1, i_proc do jj = 1, j_proc i_cell = (nx-1)/i_proc j_cell = (ny-1)/j_proc i_remain = (nx-1)-(nx-1)/i_proc*i_proc j_remain = (ny-1)-(ny-1)/j_proc*j_proc if (ii.le.i_remain) i_cell = i_cell+1 if (jj.le.j_remain) j_cell = j_cell+1 i_end = ii*i_cell+1 j_end = jj*j_cell+1 if (ii.gt.i_remain) i_end = i_end+i_remain if (jj.gt.j_remain) j_end = j_end+j_remain i_begin = i_end-i_cell j_begin = j_end-j_cell if ((ii.eq.1).and.(jj.eq.1)) then do j = j_begin,j_end+1 do i = i_begin,i_end+1 do n = 1, 6 qq(n,i,j) = w_n(n,i,j) qm(n,i,j) = w_nm(n,i,j) enddo yyn(i,j) = vis(i,j) enddo enddo else call pvmfrecv(task_ids(ii,jj),999,ibufid) do j = j_begin,j_end+1 call pvmfunpack(IREALX,qq(1,i_begin,j), . (i_end-i_begin+2)*6,1,info) call pvmfunpack(IREALX,qm(1,i_begin,j), . (i_end-i_begin+2)*6,1,info) call pvmfunpack(IREALX,yyn(i_begin,j), . i_end-i_begin+2,1,info) enddo endif enddo enddo call AFTER endif iconv = 0 enddo 300 format('it =',i5,' itsub =',i5,' log(RES) =',e12.6, ' Dt =',e12.6) 400 format(1x,i6,3(1x,e11.5)) 500 format(1x) 340 format(1x,e12.6,4(' ',e12.6),' ',i5) stop end subroutine metric include 'nsspmd2d.inc' include 'fpvm3.h' do j = 2, nfy jp = j + 1 do i = 2, nfx ip = i + 1 chx(i,j) = 0.5*( y(i,jp)+y(ip,jp)-y(i,j)-y(ip,j) ) chy(i,j) = 0.5*( x(i,jp)+x(ip,jp)-x(i,j)-x(ip,j) ) etx(i,j) = 0.5*( y(ip,j)+y(ip,jp)-y(i,j)-y(i,jp) ) ety(i,j) = 0.5*( x(ip,j)+x(ip,jp)-x(i,j)-x(i,jp) ) zac(i,j) = 1.0/ ( chx(i,j)*ety(i,j)-etx(i,j)*chy(i,j) ) chx(i,j) = zac(i,j)*chx(i,j) chy(i,j) = -zac(i,j)*chy(i,j) etx(i,j) = -zac(i,j)*etx(i,j) ety(i,j) = zac(i,j)*ety(i,j) enddo enddo if (jmtask.ge.0) then call pvmfinitsend(PvmDataRaw,ibufid) call pvmfpack(IREALX,chx(1,2),nfx+1,1,info) call pvmfpack(IREALX,chy(1,2),nfx+1,1,info) call pvmfpack(IREALX,etx(1,2),nfx+1,1,info) call pvmfpack(IREALX,ety(1,2),nfx+1,1,info) call pvmfpack(IREALX,zac(1,2),nfx+1,1,info) call pvmfsend(jmtask,3221,info) else do i = 1, nfx+1 chx(i,1) = chx(i,2) chy(i,1) = chy(i,2) etx(i,1) = etx(i,2) ety(i,1) = ety(i,2) zac(i,1) = zac(i,2) enddo endif if (jptask.ge.0) then call pvmfrecv(jptask,3221,ibufid) call pvmfunpack(IREALX,chx(1,nfy+1),nfx+1,1,info) call pvmfunpack(IREALX,chy(1,nfy+1),nfx+1,1,info) call pvmfunpack(IREALX,etx(1,nfy+1),nfx+1,1,info) call pvmfunpack(IREALX,ety(1,nfy+1),nfx+1,1,info) call pvmfunpack(IREALX,zac(1,nfy+1),nfx+1,1,info) call pvmfinitsend(PvmDataRaw,ibufid) call pvmfpack(IREALX,chx(1,nfy),nfx+1,1,info) call pvmfpack(IREALX,chy(1,nfy),nfx+1,1,info) call pvmfpack(IREALX,etx(1,nfy),nfx+1,1,info) call pvmfpack(IREALX,ety(1,nfy),nfx+1,1,info) call pvmfpack(IREALX,zac(1,nfy),nfx+1,1,info) call pvmfsend(jptask,3222,info) else do i = 1, nfx+1 chx(i,nfy+1) = chx(i,nfy) chy(i,nfy+1) = chy(i,nfy) etx(i,nfy+1) = etx(i,nfy) ety(i,nfy+1) = ety(i,nfy) zac(i,nfy+1) = zac(i,nfy) enddo endif if (jmtask.ge.0) then call pvmfrecv(jmtask,3222,ibufid) call pvmfunpack(IREALX,chx(1,1),nfx+1,1,info) call pvmfunpack(IREALX,chy(1,1),nfx+1,1,info) call pvmfunpack(IREALX,etx(1,1),nfx+1,1,info) call pvmfunpack(IREALX,ety(1,1),nfx+1,1,info) call pvmfunpack(IREALX,zac(1,1),nfx+1,1,info) endif if (imtask.ge.0) then call pvmfinitsend(PvmDataRaw,ibufid) call pvmfpack(IREALX,chx(2,1),nfy+1,imx,info) call pvmfpack(IREALX,chy(2,1),nfy+1,imx,info) call pvmfpack(IREALX,etx(2,1),nfy+1,imx,info) call pvmfpack(IREALX,ety(2,1),nfy+1,imx,info) call pvmfpack(IREALX,zac(2,1),nfy+1,imx,info) call pvmfsend(imtask,3231,info) else do j = 1, nfy+1 chx(1,j) = chx(2,j) chy(1,j) = chy(2,j) etx(1,j) = etx(2,j) ety(1,j) = ety(2,j) zac(1,j) = zac(2,j) enddo endif if (iptask.ge.0) then call pvmfrecv(iptask,3231,ibufid) call pvmfunpack(IREALX,chx(nfx+1,1),nfy+1,imx,info) call pvmfunpack(IREALX,chy(nfx+1,1),nfy+1,imx,info) call pvmfunpack(IREALX,etx(nfx+1,1),nfy+1,imx,info) call pvmfunpack(IREALX,ety(nfx+1,1),nfy+1,imx,info) call pvmfunpack(IREALX,zac(nfx+1,1),nfy+1,imx,info) call pvmfinitsend(PvmDataRaw,ibufid) call pvmfpack(IREALX,chx(nfx,1),nfy+1,imx,info) call pvmfpack(IREALX,chy(nfx,1),nfy+1,imx,info) call pvmfpack(IREALX,etx(nfx,1),nfy+1,imx,info) call pvmfpack(IREALX,ety(nfx,1),nfy+1,imx,info) call pvmfpack(IREALX,zac(nfx,1),nfy+1,imx,info) call pvmfsend(iptask,3232,info) else do j = 1, nfy+1 chx(nfx+1,j) = chx(nfx,j) chy(nfx+1,j) = chy(nfx,j) etx(nfx+1,j) = etx(nfx,j) ety(nfx+1,j) = ety(nfx,j) zac(nfx+1,j) = zac(nfx,j) enddo endif if (imtask.ge.0) then call pvmfrecv(imtask,3232,ibufid) call pvmfunpack(IREALX,chx(1,1),nfy+1,imx,info) call pvmfunpack(IREALX,chy(1,1),nfy+1,imx,info) call pvmfunpack(IREALX,etx(1,1),nfy+1,imx,info) call pvmfunpack(IREALX,ety(1,1),nfy+1,imx,info) call pvmfunpack(IREALX,zac(1,1),nfy+1,imx,info) endif return end subroutine metdum include 'nsspmd2d.inc' do j = 2, nfy jp = j + 1 do i = 1, nfx ip = i + 1 if ((ib(i,j).eq.0).and.(ib(ip,j).eq.1)) then zac(i,j) = zac(ip,j) else if ((ib(i,j).eq.1).and.(ib(ip,j).eq.0)) then zac(ip,j) = zac(i,j) endif xzac(i,j) = 0.5*( zac(i,j)+zac(ip,j) ) cx_x(i,j) = xzac(i,j)*(y(ip,jp)-y(ip,j) ) cy_x(i,j) = xzac(i,j)*(x(ip,j)-x(ip,jp) ) ex_x(i,j) = 0.5*( etx(i,j)+etx(ip,j) ) ey_x(i,j) = 0.5*( ety(i,j)+ety(ip,j) ) enddo enddo do j = 1, nfy jp = j + 1 do i = 2, nfx ip = i + 1 if ((ib(i,j).eq.0).and.(ib(i,jp).eq.1)) then zac(i,j) = zac(i,jp) else if ((ib(i,j).eq.1).and.(ib(i,jp).eq.0)) then zac(i,jp) = zac(i,j) endif yzac(i,j) = 0.5*( zac(i,j)+zac(i,jp) ) cx_y(i,j) = 0.5*( chx(i,j)+chx(i,jp) ) cy_y(i,j) = 0.5*( chy(i,j)+chy(i,jp) ) ex_y(i,j) = yzac(i,j)*(y(i,jp)-y(ip,jp) ) ey_y(i,j) = yzac(i,j)*(x(ip,jp)-x(i,jp) ) enddo enddo return end subroutine step include 'nsspmd2d.inc' dtmin = 100. dtmax = 0.0 do j = 1, nfy+1 do i = 1, nfx+1 uuu = (w(2,i,j)*chx(i,j)+w(3,i,j)*chy(i,j))/w(1,i,j) vvv = (w(2,i,j)*etx(i,j)+w(3,i,j)*ety(i,j))/w(1,i,j) cuu = c(i,j)*sqrt(chx(i,j)**2+chy(i,j)**2) cvv = c(i,j)*sqrt(etx(i,j)**2+ety(i,j)**2) chuv = abs(uuu)+cuu etuv = abs(vvv)+cvv dtc = 1./(chuv+etuv) t = c(i,j)**2 rmul = t**1.5*(1.+tcoe)/(t+tcoe) rmu = rmul/ra+vis(i,j) dtd = 1./(2.5*rmu*(chx(i,j)**2+chy(i,j)**2 . +etx(i,j)**2+ety(i,j)**2)) dt(i,j) = cfl*dtc*dtd/(dtc+dtd) dtmin = min(dtmin, dt(i,j)) dtmax = max(dtmax, dt(i,j)) enddo enddo return end subroutine physics include 'nsspmd2d.inc' common/flux/rh(2),uu(2),vv(2),pp(2),hh(2),cc(2),oo(2) dimension eflux(6,imx,jmx),fflux(6,imx,jmx) dimension ev(6,imx,jmx),fv(6,imx,jmx) do j = 2, nfy jp = j + 1 do i = 1, nfx_start - 1 ip = i + 1 if ((ib(i,j).eq.0).and.(ib(ip,j).eq.1)) then w(1,i,j) = w(1,ip,j) w(2,i,j) = -w(2,ip,j) w(3,i,j) = -w(3,ip,j) w(4,i,j) = w(4,ip,j) w(5,i,j) = -w(5,ip,j) p(i,j) = p(ip,j) vis(i,j) = -vis(ip,j) ff(i,j) = ff(ip,j) xx = 0.25*(x(ip,j)+x(ip,jp)+x(ip+1,j)+x(ip+1,jp)) yy = 0.25*(y(ip,j)+y(ip,jp)+y(ip+1,j)+y(ip+1,jp)) xc = 0.5*(x(ip,j)+x(ip,jp)) yc = 0.5*(y(ip,j)+y(ip,jp)) dyn = (xx-xc)**2+(yy-yc)**2 c ut = (w(2,ip,j)**2+w(3,ip,j)**2)/w(1,ip,j)**2 c owall = 2500*sqrt(ut/dyn) ccc = sqrt(gamma*p(ip,j)/w(1,ip,j)) t = ccc**2 rmul = t**1.5*(1.+tcoe)/(t+tcoe) owall = 60*rmul/(beta1*dyn*ra) w(6,i,j) = -w(6,ip,j)+2*owall else if ((ib(i,j).eq.1).and.(ib(ip,j).eq.0)) then w(1,ip,j) = w(1,i,j) w(2,ip,j) = -w(2,i,j) w(3,ip,j) = -w(3,i,j) w(4,ip,j) = w(4,i,j) w(5,ip,j) = -w(5,i,j) p(ip,j) = p(i,j) vis(ip,j) = -vis(i,j) ff(ip,j) = ff(i,j) xx = 0.25*(x(ip,j)+x(ip,jp)+x(i,j)+x(i,jp)) yy = 0.25*(y(ip,j)+y(ip,jp)+y(i,j)+y(i,jp)) xc = 0.5*(x(ip,j)+x(ip,jp)) yc = 0.5*(y(ip,j)+y(ip,jp)) dyn = (xx-xc)**2+(yy-yc)**2 c ut = (w(2,i,j)**2+w(3,i,j)**2)/w(1,i,j)**2 c owall = 2500*sqrt(ut/dyn) ccc = sqrt(gamma*p(i,j)/w(1,i,j)) t = ccc**2 rmul = t**1.5*(1.+tcoe)/(t+tcoe) owall = 60*rmul/(beta1*dyn*ra) w(6,ip,j) = -w(6,i,j)+2*owall endif rh(1) = w(1,i,j) uu(1) = w(2,i,j)/w(1,i,j) vv(1) = w(3,i,j)/w(1,i,j) pp(1) = p(i,j) cc(1) = w(5,i,j) oo(1) = w(6,i,j) rh(2) = w(1,ip,j) uu(2) = w(2,ip,j)/w(1,ip,j) vv(2) = w(3,ip,j)/w(1,ip,j) pp(2) = p(ip,j) cc(2) = w(5,ip,j) oo(2) = w(6,ip,j) hh(1) = gamma/gam1*pp(1)/rh(1) + 0.5*(uu(1)**2+vv(1)**2) hh(2) = gamma/gam1*pp(2)/rh(2) + 0.5*(uu(2)**2+vv(2)**2) cx = y(ip,jp) - y(ip,j) cy = x(ip,j) - x(ip,jp) call roeflux(eflux,cx,cy,i,j) enddo do i = nfx_end + 1, nfx ip = i + 1 if ((ib(i,j).eq.0).and.(ib(ip,j).eq.1)) then w(1,i,j) = w(1,ip,j) w(2,i,j) = -w(2,ip,j) w(3,i,j) = -w(3,ip,j) w(4,i,j) = w(4,ip,j) w(5,i,j) = -w(5,ip,j) p(i,j) = p(ip,j) vis(i,j) = -vis(ip,j) ff(i,j) = ff(ip,j) xx = 0.25*(x(ip,j)+x(ip,jp)+x(ip+1,j)+x(ip+1,jp)) yy = 0.25*(y(ip,j)+y(ip,jp)+y(ip+1,j)+y(ip+1,jp)) xc = 0.5*(x(ip,j)+x(ip,jp)) yc = 0.5*(y(ip,j)+y(ip,jp)) dyn = (xx-xc)**2+(yy-yc)**2 c ut = (w(2,ip,j)**2+w(3,ip,j)**2)/w(1,ip,j)**2 c owall = 2500*sqrt(ut/dyn) ccc = sqrt(gamma*p(ip,j)/w(1,ip,j)) t = ccc**2 rmul = t**1.5*(1.+tcoe)/(t+tcoe) owall = 60*rmul/(beta1*dyn*ra) w(6,i,j) = -w(6,ip,j)+2*owall else if ((ib(i,j).eq.1).and.(ib(ip,j).eq.0)) then w(1,ip,j) = w(1,i,j) w(2,ip,j) = -w(2,i,j) w(3,ip,j) = -w(3,i,j) w(4,ip,j) = w(4,i,j) w(5,ip,j) = -w(5,i,j) p(ip,j) = p(i,j) vis(ip,j) = -vis(i,j) ff(ip,j) = ff(i,j) xx = 0.25*(x(ip,j)+x(ip,jp)+x(i,j)+x(i,jp)) yy = 0.25*(y(ip,j)+y(ip,jp)+y(i,j)+y(i,jp)) xc = 0.5*(x(ip,j)+x(ip,jp)) yc = 0.5*(y(ip,j)+y(ip,jp)) dyn = (xx-xc)**2+(yy-yc)**2 c ut = (w(2,i,j)**2+w(3,i,j)**2)/w(1,i,j)**2 c owall = 2500*sqrt(ut/dyn) ccc = sqrt(gamma*p(i,j)/w(1,i,j)) t = ccc**2 rmul = t**1.5*(1.+tcoe)/(t+tcoe) owall = 60*rmul/(beta1*dyn*ra) w(6,ip,j) = -w(6,i,j)+2*owall endif rh(1) = w(1,i,j) uu(1) = w(2,i,j)/w(1,i,j) vv(1) = w(3,i,j)/w(1,i,j) pp(1) = p(i,j) cc(1) = w(5,i,j) oo(1) = w(6,i,j) rh(2) = w(1,ip,j) uu(2) = w(2,ip,j)/w(1,ip,j) vv(2) = w(3,ip,j)/w(1,ip,j) pp(2) = p(ip,j) cc(2) = w(5,ip,j) oo(2) = w(6,ip,j) hh(1) = gamma/gam1*pp(1)/rh(1) + 0.5*(uu(1)**2+vv(1)**2) hh(2) = gamma/gam1*pp(2)/rh(2) + 0.5*(uu(2)**2+vv(2)**2) cx = y(ip,jp) - y(ip,j) cy = x(ip,j) - x(ip,jp) call roeflux(eflux,cx,cy,i,j) enddo do i = nfx_start, nfx_end ip = i + 1 im = i - 1 if ((ib(i,j).eq.0).and.(ib(ip,j).eq.1)) then w(1,i,j) = w(1,ip,j) w(2,i,j) = -w(2,ip,j) w(3,i,j) = -w(3,ip,j) w(4,i,j) = w(4,ip,j) w(5,i,j) = -w(5,ip,j) p(i,j) = p(ip,j) vis(i,j) = -vis(ip,j) ff(i,j) = ff(ip,j) xx = 0.25*(x(ip,j)+x(ip,jp)+x(ip+1,j)+x(ip+1,jp)) yy = 0.25*(y(ip,j)+y(ip,jp)+y(ip+1,j)+y(ip+1,jp)) xc = 0.5*(x(ip,j)+x(ip,jp)) yc = 0.5*(y(ip,j)+y(ip,jp)) dyn = (xx-xc)**2+(yy-yc)**2 c ut = (w(2,ip,j)**2+w(3,ip,j)**2)/w(1,ip,j)**2 c owall = 2500*sqrt(ut/dyn) ccc = sqrt(gamma*p(ip,j)/w(1,ip,j)) t = ccc**2 rmul = t**1.5*(1.+tcoe)/(t+tcoe) owall = 60*rmul/(beta1*dyn*ra) w(6,i,j) = -w(6,ip,j)+2*owall else if ((ib(i,j).eq.1).and.(ib(ip,j).eq.0)) then w(1,ip,j) = w(1,i,j) w(2,ip,j) = -w(2,i,j) w(3,ip,j) = -w(3,i,j) w(4,ip,j) = w(4,i,j) w(5,ip,j) = -w(5,i,j) p(ip,j) = p(i,j) vis(ip,j) = -vis(i,j) ff(ip,j) = ff(i,j) xx = 0.25*(x(ip,j)+x(ip,jp)+x(i,j)+x(i,jp)) yy = 0.25*(y(ip,j)+y(ip,jp)+y(i,j)+y(i,jp)) xc = 0.5*(x(ip,j)+x(ip,jp)) yc = 0.5*(y(ip,j)+y(ip,jp)) dyn = (xx-xc)**2+(yy-yc)**2 c ut = (w(2,i,j)**2+w(3,i,j)**2)/w(1,i,j)**2 c owall = 2500*sqrt(ut/dyn) ccc = sqrt(gamma*p(i,j)/w(1,i,j)) t = ccc**2 rmul = t**1.5*(1.+tcoe)/(t+tcoe) owall = 60*rmul/(beta1*dyn*ra) w(6,ip,j) = -w(6,i,j)+2*owall endif rh(1) = w(1,i,j) uu(1) = w(2,i,j)/w(1,i,j) vv(1) = w(3,i,j)/w(1,i,j) pp(1) = p(i,j) cc(1) = w(5,i,j) oo(1) = w(6,i,j) rh(2) = w(1,ip,j) uu(2) = w(2,ip,j)/w(1,ip,j) vv(2) = w(3,ip,j)/w(1,ip,j) pp(2) = p(ip,j) cc(2) = w(5,ip,j) oo(2) = w(6,ip,j) call muscl(w(1,im,j),rh(1),rh(2),w(1,ip+1,j),ib(i,j),ib(ip,j)) call muscl(w(2,im,j)/w(1,im,j),uu(1),uu(2), . w(2,ip+1,j)/w(1,ip+1,j),ib(i,j),ib(ip,j)) call muscl(w(3,im,j)/w(1,im,j),vv(1),vv(2), . w(3,ip+1,j)/w(1,ip+1,j),ib(i,j),ib(ip,j)) call muscl(w(5,im,j),cc(1),cc(2),w(5,ip+1,j),ib(i,j),ib(ip,j)) call muscl(w(6,im,j),oo(1),oo(2),w(6,ip+1,j),ib(i,j),ib(ip,j)) call muscl(p(im,j),pp(1),pp(2),p(ip+1,j),ib(i,j),ib(ip,j)) hh(1) = gamma/gam1*pp(1)/rh(1) + 0.5*(uu(1)**2+vv(1)**2) hh(2) = gamma/gam1*pp(2)/rh(2) + 0.5*(uu(2)**2+vv(2)**2) cx = y(ip,jp) - y(ip,j) cy = x(ip,j) - x(ip,jp) call roeflux(eflux,cx,cy,i,j) enddo enddo do j = 1, nfy_start-1 jp = j + 1 do i = 2, nfx ip = i + 1 if ((ib(i,j).eq.0).and.(ib(i,jp).eq.1)) then w(1,i,j) = w(1,i,jp) w(2,i,j) = -w(2,i,jp) w(3,i,j) = -w(3,i,jp) w(4,i,j) = w(4,i,jp) w(5,i,j) = -w(5,i,jp) p(i,j) = p(i,jp) vis(i,j) = -vis(i,jp) ff(i,j) = ff(i,jp) xx = 0.25*(x(i,jp)+x(i+1,jp)+x(i,jp+1)+x(i+1,jp+1)) yy = 0.25*(y(i,jp)+y(i+1,jp)+y(i,jp+1)+y(i+1,jp+1)) xc = 0.5*(x(i,jp)+x(i+1,jp)) yc = 0.5*(y(i,jp)+y(i+1,jp)) dyn = (xx-xc)**2+(yy-yc)**2 c ut = (w(2,i,jp)**2+w(3,i,jp)**2)/w(1,i,jp)**2 c owall = 2500*sqrt(ut/dyn) ccc = sqrt(gamma*p(i,jp)/w(1,i,jp)) t = ccc**2 rmul = t**1.5*(1.+tcoe)/(t+tcoe) owall = 60*rmul/(beta1*dyn*ra) w(6,i,j) = -w(6,i,jp)+2*owall else if ((ib(i,j).eq.1).and.(ib(i,jp).eq.0)) then w(1,i,jp) = w(1,i,j) w(2,i,jp) = -w(2,i,j) w(3,i,jp) = -w(3,i,j) w(4,i,jp) = w(4,i,j) w(5,i,jp) = -w(5,i,j) p(i,jp) = p(i,j) vis(i,jp) = -vis(i,j) ff(i,jp) = ff(i,j) xx = 0.25*(x(i,jp)+x(i+1,jp)+x(i,j)+x(i+1,j)) yy = 0.25*(y(i,jp)+y(i+1,jp)+y(i,j)+y(i+1,j)) xc = 0.5*(x(i,jp)+x(i+1,jp)) yc = 0.5*(y(i,jp)+y(i+1,jp)) dyn = (xx-xc)**2+(yy-yc)**2 c ut = (w(2,i,j)**2+w(3,i,j)**2)/w(1,i,j)**2 c owall = 2500*sqrt(ut/dyn) ccc = sqrt(gamma*p(i,j)/w(1,i,j)) t = ccc**2 rmul = t**1.5*(1.+tcoe)/(t+tcoe) owall = 60*rmul/(beta1*dyn*ra) w(6,i,jp) = -w(6,i,j)+2*owall endif rh(1) = w(1,i,j) uu(1) = w(2,i,j)/w(1,i,j) vv(1) = w(3,i,j)/w(1,i,j) pp(1) = p(i,j) cc(1) = w(5,i,j) oo(1) = w(6,i,j) rh(2) = w(1,i,jp) uu(2) = w(2,i,jp)/w(1,i,jp) vv(2) = w(3,i,jp)/w(1,i,jp) pp(2) = p(i,jp) cc(2) = w(5,i,jp) oo(2) = w(6,i,jp) hh(1) = gamma/gam1*pp(1)/rh(1) + 0.5*(uu(1)**2+vv(1)**2) hh(2) = gamma/gam1*pp(2)/rh(2) + 0.5*(uu(2)**2+vv(2)**2) ex = y(i,jp) - y(ip,jp) ey = x(ip,jp) - x(i,jp) call roeflux(fflux,ex,ey,i,j) enddo enddo do j = nfy_end+1, nfy jp = j + 1 do i = 2, nfx ip = i + 1 if ((ib(i,j).eq.0).and.(ib(i,jp).eq.1)) then w(1,i,j) = w(1,i,jp) w(2,i,j) = -w(2,i,jp) w(3,i,j) = -w(3,i,jp) w(4,i,j) = w(4,i,jp) w(5,i,j) = -w(5,i,jp) p(i,j) = p(i,jp) vis(i,j) = -vis(i,jp) ff(i,j) = ff(i,jp) xx = 0.25*(x(i,jp)+x(i+1,jp)+x(i,jp+1)+x(i+1,jp+1)) yy = 0.25*(y(i,jp)+y(i+1,jp)+y(i,jp+1)+y(i+1,jp+1)) xc = 0.5*(x(i,jp)+x(i+1,jp)) yc = 0.5*(y(i,jp)+y(i+1,jp)) dyn = (xx-xc)**2+(yy-yc)**2 c ut = (w(2,i,jp)**2+w(3,i,jp)**2)/w(1,i,jp)**2 c owall = 2500*sqrt(ut/dyn) ccc = sqrt(gamma*p(i,jp)/w(1,i,jp)) t = ccc**2 rmul = t**1.5*(1.+tcoe)/(t+tcoe) owall = 60*rmul/(beta1*dyn*ra) w(6,i,j) = -w(6,i,jp)+2*owall else if ((ib(i,j).eq.1).and.(ib(i,jp).eq.0)) then w(1,i,jp) = w(1,i,j) w(2,i,jp) = -w(2,i,j) w(3,i,jp) = -w(3,i,j) w(4,i,jp) = w(4,i,j) w(5,i,jp) = -w(5,i,j) p(i,jp) = p(i,j) vis(i,jp) = -vis(i,j) ff(i,jp) = ff(i,j) xx = 0.25*(x(i,jp)+x(i+1,jp)+x(i,j)+x(i+1,j)) yy = 0.25*(y(i,jp)+y(i+1,jp)+y(i,j)+y(i+1,j)) xc = 0.5*(x(i,jp)+x(i+1,jp)) yc = 0.5*(y(i,jp)+y(i+1,jp)) dyn = (xx-xc)**2+(yy-yc)**2 c ut = (w(2,i,j)**2+w(3,i,j)**2)/w(1,i,j)**2 c owall = 2500*sqrt(ut/dyn) ccc = sqrt(gamma*p(i,j)/w(1,i,j)) t = ccc**2 rmul = t**1.5*(1.+tcoe)/(t+tcoe) owall = 60*rmul/(beta1*dyn*ra) w(6,i,jp) = -w(6,i,j)+2*owall endif rh(1) = w(1,i,j) uu(1) = w(2,i,j)/w(1,i,j) vv(1) = w(3,i,j)/w(1,i,j) pp(1) = p(i,j) cc(1) = w(5,i,j) oo(1) = w(6,i,j) rh(2) = w(1,i,jp) uu(2) = w(2,i,jp)/w(1,i,jp) vv(2) = w(3,i,jp)/w(1,i,jp) pp(2) = p(i,jp) cc(2) = w(5,i,jp) oo(2) = w(6,i,jp) hh(1) = gamma/gam1*pp(1)/rh(1) + 0.5*(uu(1)**2+vv(1)**2) hh(2) = gamma/gam1*pp(2)/rh(2) + 0.5*(uu(2)**2+vv(2)**2) ex = y(i,jp) - y(ip,jp) ey = x(ip,jp) - x(i,jp) call roeflux(fflux,ex,ey,i,j) enddo enddo do j = nfy_start, nfy_end jp = j + 1 jm = j - 1 do i = 2, nfx ip = i + 1 if ((ib(i,j).eq.0).and.(ib(i,jp).eq.1)) then w(1,i,j) = w(1,i,jp) w(2,i,j) = -w(2,i,jp) w(3,i,j) = -w(3,i,jp) w(4,i,j) = w(4,i,jp) w(5,i,j) = -w(5,i,jp) p(i,j) = p(i,jp) vis(i,j) = -vis(i,jp) ff(i,j) = ff(i,jp) xx = 0.25*(x(i,jp)+x(i+1,jp)+x(i,jp+1)+x(i+1,jp+1)) yy = 0.25*(y(i,jp)+y(i+1,jp)+y(i,jp+1)+y(i+1,jp+1)) xc = 0.5*(x(i,jp)+x(i+1,jp)) yc = 0.5*(y(i,jp)+y(i+1,jp)) dyn = (xx-xc)**2+(yy-yc)**2 c ut = (w(2,i,jp)**2+w(3,i,jp)**2)/w(1,i,jp)**2 c owall = 2500*sqrt(ut/dyn) ccc = sqrt(gamma*p(i,jp)/w(1,i,jp)) t = ccc**2 rmul = t**1.5*(1.+tcoe)/(t+tcoe) owall = 60*rmul/(beta1*dyn*ra) w(6,i,j) = -w(6,i,jp)+2*owall else if ((ib(i,j).eq.1).and.(ib(i,jp).eq.0)) then w(1,i,jp) = w(1,i,j) w(2,i,jp) = -w(2,i,j) w(3,i,jp) = -w(3,i,j) w(4,i,jp) = w(4,i,j) w(5,i,jp) = -w(5,i,j) p(i,jp) = p(i,j) vis(i,jp) = -vis(i,j) ff(i,jp) = ff(i,j) xx = 0.25*(x(i,jp)+x(i+1,jp)+x(i,j)+x(i+1,j)) yy = 0.25*(y(i,jp)+y(i+1,jp)+y(i,j)+y(i+1,j)) xc = 0.5*(x(i,jp)+x(i+1,jp)) yc = 0.5*(y(i,jp)+y(i+1,jp)) dyn = (xx-xc)**2+(yy-yc)**2 c ut = (w(2,i,j)**2+w(3,i,j)**2)/w(1,i,j)**2 c owall = 2500*sqrt(ut/dyn) ccc = sqrt(gamma*p(i,j)/w(1,i,j)) t = ccc**2 rmul = t**1.5*(1.+tcoe)/(t+tcoe) owall = 60*rmul/(beta1*dyn*ra) w(6,i,jp) = -w(6,i,j)+2*owall endif rh(1) = w(1,i,j) uu(1) = w(2,i,j)/w(1,i,j) vv(1) = w(3,i,j)/w(1,i,j) pp(1) = p(i,j) cc(1) = w(5,i,j) oo(1) = w(6,i,j) rh(2) = w(1,i,jp) uu(2) = w(2,i,jp)/w(1,i,jp) vv(2) = w(3,i,jp)/w(1,i,jp) pp(2) = p(i,jp) cc(2) = w(5,i,jp) oo(2) = w(6,i,jp) call muscl(w(1,i,jm),rh(1),rh(2),w(1,i,jp+1),ib(i,j),ib(i,jp)) call muscl(w(2,i,jm)/w(1,i,jm),uu(1),uu(2), . w(2,i,jp+1)/w(1,i,jp+1),ib(i,j),ib(i,jp)) call muscl(w(3,i,jm)/w(1,i,jm),vv(1),vv(2), . w(3,i,jp+1)/w(1,i,jp+1),ib(i,j),ib(i,jp)) call muscl(w(5,i,jm),cc(1),cc(2),w(5,i,jp+1),ib(i,j),ib(i,jp)) call muscl(w(6,i,jm),oo(1),oo(2),w(6,i,jp+1),ib(i,j),ib(i,jp)) call muscl(p(i,jm),pp(1),pp(2),p(i,jp+1),ib(i,j),ib(i,jp)) hh(1) = gamma/gam1*pp(1)/rh(1) + 0.5*(uu(1)**2+vv(1)**2) hh(2) = gamma/gam1*pp(2)/rh(2) + 0.5*(uu(2)**2+vv(2)**2) ex = y(i,jp) - y(ip,jp) ey = x(ip,jp) - x(i,jp) call roeflux(fflux,ex,ey,i,j) enddo enddo do j = 2, nfy jp = j + 1 do i = 1, nfx ip = i + 1 t = 0.5*(c(i,j)**2+c(ip,j)**2) rmul = t**1.5*(1.+tcoe)/(t+tcoe) rmut = 0.5*ra*(vis(i,j)+vis(ip,j)) rmu = rmul+rmut prcoef = rmul/prl + rmut/prt pr = rmu/prcoef f2 = 0.5*(ff(i,j)+ff(ip,j)) uip = w(2,ip,j)/w(1,ip,j) ui = w(2,i,j)/w(1,i,j) vip = w(3,ip,j)/w(1,ip,j) vi = w(3,i,j)/w(1,i,j) c2ip = c(ip,j)**2 c2i = c(i,j)**2 uface = 0.5*(uip+ui) vface = 0.5*(vip+vi) uup = unod(ip,jp) udown = unod(ip,j) vup = vnod(ip,jp) vdown = vnod(ip,j) c2up = cnod(ip,jp)**2 c2down = cnod(ip,j)**2 uch = uip-ui vch = vip-vi uet = uup-udown vet = vup-vdown c2ch = c2ip-c2i c2et = c2up-c2down rkup = rknod(ip,jp) rkdown = rknod(ip,j) rwup = rwnod(ip,jp) rwdown = rwnod(ip,j) rkch = w(5,ip,j) - w(5,i,j) rwch = w(6,ip,j) - w(6,i,j) rket = rkup - rkdown rwet = rwup - rwdown a1 = 4./3.*cx_x(i,j)**2+cy_x(i,j)**2 a2 = 4./3.*cx_x(i,j)*ex_x(i,j)+cy_x(i,j)*ey_x(i,j) a3 = 1./3.*cx_x(i,j)*cy_x(i,j) a4 = cy_x(i,j)*ex_x(i,j)-2./3.*cx_x(i,j)*ey_x(i,j) b1 = cx_x(i,j)**2+4./3.*cy_x(i,j)**2 b2 = cx_x(i,j)*ex_x(i,j)+4./3.*cy_x(i,j)*ey_x(i,j) b3 = 1./3.*cx_x(i,j)*cy_x(i,j) b4 = cx_x(i,j)*ey_x(i,j)-2./3.*cy_x(i,j)*ex_x(i,j) c1 = cx_x(i,j)**2+cy_x(i,j)**2 c2 = cx_x(i,j)*ex_x(i,j)+cy_x(i,j)*ey_x(i,j) omet = cx_x(i,j)**2+cy_x(i,j)**2 pmet = cx_x(i,j)*ex_x(i,j)+cy_x(i,j)*ey_x(i,j) sigmak = f2*sigmak1+(1.-f2)*sigmak2 sigmao = f2*sigmao1+(1.-f2)*sigmao2 aa = rmul+sigmak*rmut bb = rmul+sigmao*rmut ev(1,i,j) = 0.0 ev(2,i,j) = a1*uch+a2*uet+a3*vch+a4*vet ev(3,i,j) = b1*vch+b2*vet+b3*uch+b4*uet ev(4,i,j) = uface*ev(2,i,j)+vface*ev(3,i,j) . +(c1*c2ch+c2*c2et)/(gam1*pr) ev(5,i,j) = aa*(omet*rkch+pmet*rket) ev(6,i,j) = bb*(omet*rwch+pmet*rwet) do n = 2, 4 ev(n,i,j) = rmu*ev(n,i,j)/xzac(i,j) enddo do n = 5, 6 ev(n,i,j) = ev(n,i,j)/xzac(i,j) enddo enddo enddo do j = 1, nfy jp = j + 1 do i = 2, nfx ip = i + 1 t = 0.5*(c(i,j)**2+c(i,jp)**2) rmul = t**1.5*(1.+tcoe)/(t+tcoe) rmut = 0.5*ra*(vis(i,j)+vis(i,jp)) rmu = rmul + rmut prcoef = rmul/prl + rmut/prt pr = rmu/prcoef f2 = 0.5*(ff(i,j)+ff(i,jp)) ujp = w(2,i,jp)/w(1,i,jp) uj = w(2,i,j)/w(1,i,j) vjp = w(3,i,jp)/w(1,i,jp) vj = w(3,i,j)/w(1,i,j) c2jp = c(i,jp)**2 c2j = c(i,j)**2 uface = 0.5*(ujp+uj) vface = 0.5*(vjp+vj) uup = unod(ip,jp) udown = unod(i,jp) vup = vnod(ip,jp) vdown = vnod(i,jp) c2up = cnod(ip,jp)**2 c2down = cnod(i,jp)**2 uch = uup-udown vch = vup-vdown uet = ujp-uj vet = vjp-vj c2ch = c2up-c2down c2et = c2jp-c2j rkup = rknod(ip,jp) rkdown = rknod(i,jp) rwup = rwnod(ip,jp) rwdown = rwnod(i,jp) rkch = rkup - rkdown rwch = rwup - rwdown rket = w(5,i,jp) - w(5,i,j) rwet = w(6,i,jp) - w(6,i,j) a1 = 4./3.*ex_y(i,j)**2+ey_y(i,j)**2 a2 = 4./3.*cx_y(i,j)*ex_y(i,j)+cy_y(i,j)*ey_y(i,j) a3 = 1./3.*ex_y(i,j)*ey_y(i,j) a4 = cx_y(i,j)*ey_y(i,j)-2./3.*cy_y(i,j)*ex_y(i,j) b1 = ex_y(i,j)**2+4./3.*ey_y(i,j)**2 b2 = cx_y(i,j)*ex_y(i,j)+4./3.*cy_y(i,j)*ey_y(i,j) b3 = 1./3.*ex_y(i,j)*ey_y(i,j) b4 = cy_y(i,j)*ex_y(i,j)-2./3.*cx_y(i,j)*ey_y(i,j) c1 = ex_y(i,j)**2+ey_y(i,j)**2 c2 = cx_y(i,j)*ex_y(i,j)+cy_y(i,j)*ey_y(i,j) omet = cx_y(i,j)*ex_y(i,j)+cy_y(i,j)*ey_y(i,j) pmet = ex_y(i,j)**2+ey_y(i,j)**2 sigmak = f2*sigmak1+(1.-f2)*sigmak2 sigmao = f2*sigmao1+(1.-f2)*sigmao2 aa = rmul+sigmak*rmut bb = rmul+sigmao*rmut fv(1,i,j) = 0.0 fv(2,i,j) = a1*uet+a2*uch+a3*vet+a4*vch fv(3,i,j) = b1*vet+b2*vch+b3*uet+b4*uch fv(4,i,j) = uface*fv(2,i,j)+vface*fv(3,i,j) . +(c1*c2et+c2*c2ch)/(gam1*pr) fv(5,i,j) = aa*(omet*rkch+pmet*rket) fv(6,i,j) = bb*(omet*rwch+pmet*rwet) do n = 2, 4 fv(n,i,j) = rmu*fv(n,i,j)/yzac(i,j) enddo do n = 5, 6 fv(n,i,j) = fv(n,i,j)/yzac(i,j) enddo enddo enddo do j = 2, nfy jm = j - 1 do i = 2, nfx im = i - 1 do n = 1, 6 rhs(n,i,j) = eflux(n,i,j)-eflux(n,im,j) . +fflux(n,i,j)-fflux(n,i,jm) . -(ev(n,i,j)-ev(n,im,j) . + fv(n,i,j)-fv(n,i,jm))/ra enddo enddo enddo return end subroutine roeflux(eflux,cex,cey,i,j) include 'nsspmd2d.inc' common/flux/rh(2),uu(2),vv(2),pp(2),hh(2),cc(2),oo(2) dimension eflux(6,imx,jmx),fl(4),fr(4), . df1(4),df2(4),df3(4),df4(4) dimension ram(4) dss = sqrt(cey**2+cex**2) chxn = cex/dss chyn = cey/dss rrr = sqrt(rh(2)/rh(1)) rrrp1 = rrr + 1.0 rhoav = rrr*rh(1) uav = ( rrr*uu(2)+uu(1) )/rrrp1 vav = ( rrr*vv(2)+vv(1) )/rrrp1 hav = ( rrr*hh(2)+hh(1) )/rrrp1 cavsq = gam1*( hav-0.5*(uav**2+vav**2) ) cav = sqrt(cavsq) qavsq = uav**2+vav**2 dqsq = uu(2)**2+vv(2)**2-uu(1)**2-vv(1)**2 contra= chxn*uav+chyn*vav dp = pp(2)-pp(1) drho = rh(2)-rh(1) du = uu(2)-uu(1) dv = vv(2)-vv(1) ram(1) = abs( contra ) ram(2) = ram(1) ram(3) = abs( contra+cav ) ram(4) = abs( contra-cav ) do n= 1, 4 if(ram(n).lt.eee)then ram(n) = 0.5*(ram(n)**2/eee + eee) endif enddo ram(1) = ram(1)*dss ram(2) = ram(2)*dss ram(3) = ram(3)*dss ram(4) = ram(4)*dss chduv = chxn*du +chyn*dv coeff1 = ram(1)*(drho-dp/cavsq) coeff2 = ram(2)*rhoav coeff3 = ram(3)*0.5*(dp+rhoav*cav*chduv)/cavsq coeff4 = ram(4)*0.5*(dp-rhoav*cav*chduv)/cavsq df1(1) = coeff1 df1(2) = coeff1*uav df1(3) = coeff1*vav df1(4) = 0.5*coeff1*qavsq df2(1) = 0.0 df2(2) = coeff2*(du - chxn*chduv) df2(3) = coeff2*(dv - chyn*chduv) df2(4) = coeff2*(0.5*dqsq - contra*chduv) df3(1) = coeff3 df3(2) = coeff3*(uav + chxn*cav) df3(3) = coeff3*(vav + chyn*cav) df3(4) = coeff3*(hav + cav*contra) df4(1) = coeff4 df4(2) = coeff4*(uav - chxn*cav) df4(3) = coeff4*(vav - chyn*cav) df4(4) = coeff4*(hav - cav*contra) commo = rh(1)*( cex*uu(1) + cey*vv(1) ) fl(1) = commo fl(2) = uu(1)*commo+cex*pp(1) fl(3) = vv(1)*commo+cey*pp(1) fl(4) = hh(1)*commo commo = rh(2)*( cex*uu(2) + cey*vv(2) ) fr(1) = commo fr(2) = uu(2)*commo+cex*pp(2) fr(3) = vv(2)*commo+cey*pp(2) fr(4) = hh(2)*commo do n= 1, 4 eflux(n,i,j)=0.5*(fl(n)+fr(n)) . +0.5*(-df1(n)-df2(n)-df3(n)-df4(n)) enddo eflux(5,i,j)=0.5*(cc(1)*fl(1)+cc(2)*fr(1)) eflux(6,i,j)=0.5*(oo(1)*fl(1)+oo(2)*fr(1)) delk = (cc(2) - cc(1))*rhoav delo = (oo(2) - oo(1))*rhoav eflux(5,i,j) = eflux(5,i,j)-0.5*ram(1)*delk eflux(6,i,j) = eflux(6,i,j)-0.5*ram(1)*delo return end subroutine node include 'nsspmd2d.inc' do i = 2, nfx+1 im = i - 1 do j = 2, nfy+1 jm = j - 1 uij = w(2,i,j)/w(1,i,j) uimj = w(2,im,j)/w(1,im,j) uimjm = w(2,im,jm)/w(1,im,jm) uijm = w(2,i,jm)/w(1,i,jm) unod(i,j) = 0.25*(uij+uimj+uimjm+uijm) vij = w(3,i,j)/w(1,i,j) vimj = w(3,im,j)/w(1,im,j) vimjm = w(3,im,jm)/w(1,im,jm) vijm = w(3,i,jm)/w(1,i,jm) vnod(i,j) = 0.25*(vij+vimj+vimjm+vijm) cnod(i,j) = 0.25*(c(i,j)+c(im,j)+c(im,jm)+c(i,jm)) rknod(i,j) = 0.25*(w(5,i,j)+w(5,im,j)+w(5,im,jm)+w(5,i,jm)) rwnod(i,j) = 0.25*(w(6,i,j)+w(6,im,j)+w(6,im,jm)+w(6,i,jm)) enddo enddo return end subroutine solve include 'nsspmd2d.inc' dimension roa(imx,jmx),rob(imx,jmx),rov(imx,jmx) dimension apdw1(6),amdw1(6),bpdw1(6),bmdw1(6) dimension amdw(6,imx,jmx),bmdw(6,imx,jmx) dimension sr(6,imx,jmx), scr(6,imx,jmx) if (it.eq.1) then factor = 1./1.5 factor1 = 1.0 else factor = 1.0 factor1 = 1.5 endif if (delt.lt.0) then factor = 0.0 factor1 = 0.0 endif do j = 1, nfy+1 do i = 1, nfx+1 cx = chx(i,j)/zac(i,j) cy = chy(i,j)/zac(i,j) uuu = (cx*w(2,i,j)+cy*w(3,i,j))/w(1,i,j) cuu = c(i,j)*sqrt(cx**2+cy**2) roa(i,j) = cappa*(abs(uuu)+cuu) ex = etx(i,j)/zac(i,j) ey = ety(i,j)/zac(i,j) vvv = (ex*w(2,i,j)+ey*w(3,i,j))/w(1,i,j) cvv = c(i,j)*sqrt(ex**2+ey**2) rob(i,j) = cappa*(abs(vvv)+cvv) t = c(i,j)**2 rmul = t**1.5*(1.+tcoe)/(t+tcoe) rmut = ra*vis(i,j) rmu = rmul+rmut prcoef = rmul/prl + rmut/prt pr = rmu/prcoef rov(i,j) = gamma*rmu/(w(1,i,j)*pr*ra*zac(i,j)) do n = 1, 4 scr(n,i,j) = 1./((1./dt(i,j)+factor1/delt)/zac(i,j) . +roa(i,j)+rob(i,j)+4*rov(i,j)) enddo scr(5,i,j) = 1./((1./dt(i,j)+factor1/delt+2*betaA*w(6,i,j)) . /zac(i,j)+roa(i,j)+rob(i,j)+4*rov(i,j)) betao = ff(i,j)*beta1+(1.-ff(i,j))*beta2 scr(6,i,j) = 1./((1./dt(i,j)+factor1/delt+2*betao*w(6,i,j)) . /zac(i,j)+roa(i,j)+rob(i,j)+4*rov(i,j)) enddo enddo do j = 2, nfy jp = j + 1 jm = j - 1 do i = 2, nfx ip = i + 1 im = i - 1 do n = 1, 4 sr(n,i,j) = factor*(1.5*w(n,i,j)-2.0*w_n(n,i,j) . +0.5*w_nm(n,i,j))/(zac(i,j)*delt) enddo do n = 5, 6 sr(n,i,j) = factor*(1.5*w(n,i,j)*w(1,i,j)-2.0*w_n(n,i,j) . +0.5*w_nm(n,i,j))/(zac(i,j)*delt) enddo uch = 0.5*(w(2,ip,j)/w(1,ip,j)-w(2,im,j)/w(1,im,j)) uet = 0.5*(w(2,i,jp)/w(1,i,jp)-w(2,i,jm)/w(1,i,jm)) vch = 0.5*(w(3,ip,j)/w(1,ip,j)-w(3,im,j)/w(1,im,j)) vet = 0.5*(w(3,i,jp)/w(1,i,jp)-w(3,i,jm)/w(1,i,jm)) rkch = 0.5*(w(5,ip,j)-w(5,im,j)) rket = 0.5*(w(5,i,jp)-w(5,i,jm)) rwch = 0.5*(w(6,ip,j)-w(6,im,j)) rwet = 0.5*(w(6,i,jp)-w(6,i,jm)) ux = uch*chx(i,j)+uet*etx(i,j) uy = uch*chy(i,j)+uet*ety(i,j) vx = vch*chx(i,j)+vet*etx(i,j) vy = vch*chy(i,j)+vet*ety(i,j) rkx = rkch*chx(i,j)+rket*etx(i,j) rky = rkch*chy(i,j)+rket*ety(i,j) rwx = rwch*chx(i,j)+rwet*etx(i,j) rwy = rwch*chy(i,j)+rwet*ety(i,j) Pr = 4./3.*(ux**2+vy**2-ux*vy)+(uy+vx)**2 . -2./3.*w(1,i,j)*w(5,i,j)*(ux+vy)/vis(i,j) Pk = vis(i,j)*Pr Dk = betaA*w(1,i,j)*w(5,i,j)*w(6,i,j) Pk = min(Pk,20.*Dk) sr(5,i,j) = sr(5,i,j)-(Pk-Dk)/zac(i,j) betao = ff(i,j)*beta1+(1.-ff(i,j))*beta2 tgamma = ff(i,j)*tgamma1+(1.-ff(i,j))*tgamma2 sr(6,i,j) = sr(6,i,j) . -(tgamma*Pr-betao*w(6,i,j)**2)*w(1,i,j)/zac(i,j) . -2.*w(1,i,j)*(1.-ff(i,j))*sigmao2*(rkx*rwx+rky*rwy) . /(w(6,i,j)*zac(i,j)) do n = 1, 6 dw(n,i,j) = ib(i,j)*(-rhs(n,i,j)-sr(n,i,j))*scr(n,i,j) enddo enddo enddo do l = 1, isweep call DW_CHANGE do j = 2, nfy jm = j - 1 jp = j + 1 do i = 2, nfx im = i - 1 ip = i + 1 ox = chx(im,j)/zac(im,j) oy = chy(im,j)/zac(im,j) uu = w(2,im,j)/w(1,im,j) vv = w(3,im,j)/w(1,im,j) Q = 0.5*gam1*(uu**2+vv**2) H = (w(4,im,j)+p(im,j))/w(1,im,j) U = ox*uu+oy*vv V = oy*uu+ox*vv dU = uu*dw1(2,im,j)+vv*dw1(3,im,j) dUU = U*dw1(2,im,j)+V*dw1(3,im,j) dVV = V*dw1(2,im,j)+U*dw1(3,im,j) dww = ox*dw1(2,im,j)+oy*dw1(3,im,j) apdw1(1) = (roa(im,j)-2*rov(im,j))*dw1(1,im,j)+dww apdw1(2) = (-U*uu+ox*Q)*dw1(1,im,j)+dUU-gamma*ox*dU . +(2*ox*uu+roa(im,j)-2*rov(im,j))*dw1(2,im,j) . +gam1*ox*dw1(4,im,j) apdw1(3) = (-U*vv+oy*Q)*dw1(1,im,j)+dVV-gamma*oy*dU . +(2*oy*vv+roa(im,j)-2*rov(im,j))*dw1(3,im,j) . +gam1*oy*dw1(4,im,j) apdw1(4) = U*(Q-H)*dw1(1,im,j)+H*dww-gam1*U*dU . +(gamma*U+roa(im,j)-2*rov(im,j))*dw1(4,im,j) apdw1(5) = (U+roa(im,j)-2*rov(im,j))*dw1(5,im,j) apdw1(6) = (U+roa(im,j)-2*rov(im,j))*dw1(6,im,j) ox = etx(i,jm)/zac(i,jm) oy = ety(i,jm)/zac(i,jm) uu = w(2,i,jm)/w(1,i,jm) vv = w(3,i,jm)/w(1,i,jm) Q = 0.5*gam1*(uu**2+vv**2) H = (w(4,i,jm)+p(i,jm))/w(1,i,jm) U = ox*uu+oy*vv V = oy*uu+ox*vv dU = uu*dw1(2,i,jm)+vv*dw1(3,i,jm) dUU = U*dw1(2,i,jm)+V*dw1(3,i,jm) dVV = V*dw1(2,i,jm)+U*dw1(3,i,jm) dww = ox*dw1(2,i,jm)+oy*dw1(3,i,jm) bpdw1(1) = (rob(i,jm)-2*rov(i,jm))*dw1(1,i,jm)+dww bpdw1(2) = (-U*uu+ox*Q)*dw1(1,i,jm)+dUU-gamma*ox*dU . +(2*ox*uu+rob(i,jm)-2*rov(i,jm))*dw1(2,i,jm) . +gam1*ox*dw1(4,i,jm) bpdw1(3) = (-U*vv+oy*Q)*dw1(1,i,jm)+dVV-gamma*oy*dU . +(2*oy*vv+rob(i,jm)-2*rov(i,jm))*dw1(3,i,jm) . +gam1*oy*dw1(4,i,jm) bpdw1(4) = U*(Q-H)*dw1(1,i,jm)+H*dww-gam1*U*dU . +(gamma*U+rob(i,jm)-2*rov(i,jm))*dw1(4,i,jm) bpdw1(5) = (U+rob(i,jm)-2*rov(i,jm))*dw1(5,i,jm) bpdw1(6) = (U+rob(i,jm)-2*rov(i,jm))*dw1(6,i,jm) ox = chx(ip,j)/zac(ip,j) oy = chy(ip,j)/zac(ip,j) uu = w(2,ip,j)/w(1,ip,j) vv = w(3,ip,j)/w(1,ip,j) Q = 0.5*gam1*(uu**2+vv**2) H = (w(4,ip,j)+p(ip,j))/w(1,ip,j) U = ox*uu+oy*vv V = oy*uu+ox*vv dU = uu*dw(2,ip,j)+vv*dw(3,ip,j) dUU = U*dw(2,ip,j)+V*dw(3,ip,j) dVV = V*dw(2,ip,j)+U*dw(3,ip,j) dww = ox*dw(2,ip,j)+oy*dw(3,ip,j) amdw(1,i,j) = (-roa(ip,j)+2*rov(ip,j))*dw(1,ip,j)+dww amdw(2,i,j) = (-U*uu+ox*Q)*dw(1,ip,j)+dUU-gamma*ox*dU . +(2*ox*uu-roa(ip,j)+2*rov(ip,j))*dw(2,ip,j) . +gam1*ox*dw(4,ip,j) amdw(3,i,j) = (-U*vv+oy*Q)*dw(1,ip,j)+dVV-gamma*oy*dU . +(2*oy*vv-roa(ip,j)+2*rov(ip,j))*dw(3,ip,j) . +gam1*oy*dw(4,ip,j) amdw(4,i,j) = U*(Q-H)*dw(1,ip,j)+H*dww-gam1*U*dU . +(gamma*U-roa(ip,j)+2*rov(ip,j))*dw(4,ip,j) amdw(5,i,j) = (U-roa(ip,j)+2*rov(ip,j))*dw(5,ip,j) amdw(6,i,j) = (U-roa(ip,j)+2*rov(ip,j))*dw(6,ip,j) ox = etx(i,jp)/zac(i,jp) oy = ety(i,jp)/zac(i,jp) uu = w(2,i,jp)/w(1,i,jp) vv = w(3,i,jp)/w(1,i,jp) Q = 0.5*gam1*(uu**2+vv**2) H = (w(4,i,jp)+p(i,jp))/w(1,i,jp) U = ox*uu+oy*vv V = oy*uu+ox*vv dU = uu*dw(2,i,jp)+vv*dw(3,i,jp) dUU = U*dw(2,i,jp)+V*dw(3,i,jp) dVV = V*dw(2,i,jp)+U*dw(3,i,jp) dww = ox*dw(2,i,jp)+oy*dw(3,i,jp) bmdw(1,i,j) = (-rob(i,jp)+2*rov(i,jp))*dw(1,i,jp)+dww bmdw(2,i,j) = (-U*uu+ox*Q)*dw(1,i,jp)+dUU-gamma*ox*dU . +(2*ox*uu-rob(i,jp)+2*rov(i,jp))*dw(2,i,jp) . +gam1*ox*dw(4,i,jp) bmdw(3,i,j) = (-U*vv+oy*Q)*dw(1,i,jp)+dVV-gamma*oy*dU . +(2*oy*vv-rob(i,jp)+2*rov(i,jp))*dw(3,i,jp) . +gam1*oy*dw(4,i,jp) bmdw(4,i,j) = U*(Q-H)*dw(1,i,jp)+H*dww-gam1*U*dU . +(gamma*U-rob(i,jp)+2*rov(i,jp))*dw(4,i,jp) bmdw(5,i,j) = (U-rob(i,jp)+2*rov(i,jp))*dw(5,i,jp) bmdw(6,i,j) = (U-rob(i,jp)+2*rov(i,jp))*dw(6,i,jp) do n = 1, 6 dw1(n,i,j) = (-rhs(n,i,j) - sr(n,i,j) . + 0.5*(apdw1(n)+bpdw1(n)-amdw(n,i,j)-bmdw(n,i,j))) . * ib(i,j)*scr(n,i,j) enddo enddo enddo do j = nfy, 2, -1 jp = j + 1 do i = nfx, 2, -1 ip = i + 1 ox = chx(ip,j)/zac(ip,j) oy = chy(ip,j)/zac(ip,j) uu = w(2,ip,j)/w(1,ip,j) vv = w(3,ip,j)/w(1,ip,j) Q = 0.5*gam1*(uu**2+vv**2) H = (w(4,ip,j)+p(ip,j))/w(1,ip,j) U = ox*uu+oy*vv V = oy*uu+ox*vv dU = uu*dw1(2,ip,j)+vv*dw1(3,ip,j) dUU = U*dw1(2,ip,j)+V*dw1(3,ip,j) dVV = V*dw1(2,ip,j)+U*dw1(3,ip,j) dw1w = ox*dw1(2,ip,j)+oy*dw1(3,ip,j) amdw1(1) = (-roa(ip,j)+2*rov(ip,j))*dw1(1,ip,j)+dw1w amdw1(2) = (-U*uu+ox*Q)*dw1(1,ip,j)+dUU-gamma*ox*dU . +(2*ox*uu-roa(ip,j)+2*rov(ip,j))*dw1(2,ip,j) . +gam1*ox*dw1(4,ip,j) amdw1(3) = (-U*vv+oy*Q)*dw1(1,ip,j)+dVV-gamma*oy*dU . +(2*oy*vv-roa(ip,j)+2*rov(ip,j))*dw1(3,ip,j) . +gam1*oy*dw1(4,ip,j) amdw1(4) = U*(Q-H)*dw1(1,ip,j)+H*dw1w-gam1*U*dU . +(gamma*U-roa(ip,j)+2*rov(ip,j))*dw1(4,ip,j) amdw1(5) = (U-roa(ip,j)+2*rov(ip,j))*dw1(5,ip,j) amdw1(6) = (U-roa(ip,j)+2*rov(ip,j))*dw1(6,ip,j) ox = etx(i,jp)/zac(i,jp) oy = ety(i,jp)/zac(i,jp) uu = w(2,i,jp)/w(1,i,jp) vv = w(3,i,jp)/w(1,i,jp) Q = 0.5*gam1*(uu**2+vv**2) H = (w(4,i,jp)+p(i,jp))/w(1,i,jp) U = ox*uu+oy*vv V = oy*uu+ox*vv dU = uu*dw1(2,i,jp)+vv*dw1(3,i,jp) dUU = U*dw1(2,i,jp)+V*dw1(3,i,jp) dVV = V*dw1(2,i,jp)+U*dw1(3,i,jp) dw1w = ox*dw1(2,i,jp)+oy*dw1(3,i,jp) bmdw1(1) = (-rob(i,jp)+2*rov(i,jp))*dw1(1,i,jp)+dw1w bmdw1(2) = (-U*uu+ox*Q)*dw1(1,i,jp)+dUU-gamma*ox*dU . +(2*ox*uu-rob(i,jp)+2*rov(i,jp))*dw1(2,i,jp) . +gam1*ox*dw1(4,i,jp) bmdw1(3) = (-U*vv+oy*Q)*dw1(1,i,jp)+dVV-gamma*oy*dU . +(2*oy*vv-rob(i,jp)+2*rov(i,jp))*dw1(3,i,jp) . +gam1*oy*dw1(4,i,jp) bmdw1(4) = U*(Q-H)*dw1(1,i,jp)+H*dw1w-gam1*U*dU . +(gamma*U-rob(i,jp)+2*rov(i,jp))*dw1(4,i,jp) bmdw1(5) = (U-rob(i,jp)+2*rov(i,jp))*dw1(5,i,jp) bmdw1(6) = (U-rob(i,jp)+2*rov(i,jp))*dw1(6,i,jp) do n = 1, 6 dw1(n,i,j) = (dw1(n,i,j) . - (0.5*(amdw1(n)+bmdw1(n)-amdw(n,i,j)-bmdw(n,i,j))) . *scr(n,i,j))*ib(i,j) enddo enddo enddo do i = 2, nfx do j = 2, nfy do n = 1, 6 dw(n,i,j) = ib(i,j)*dw1(n,i,j) enddo enddo enddo enddo return end subroutine dw_change include 'nsspmd2d.inc' include 'fpvm3.h' if (jmtask.ge.0) then call pvmfinitsend(PvmDataRaw,ibufid) call pvmfpack(IREALX,dw(1,2,2),(nfx-1)*6,1,info) call pvmfsend(jmtask,1021,info) endif if (jptask.ge.0) then call pvmfrecv(jptask,1021,ibufid) call pvmfunpack(IREALX,dw(1,2,nfy+1),(nfx-1)*6,1,info) do i = 2, nfx dw1(1,i,nfy+1) = dw(1,i,nfy+1) dw1(2,i,nfy+1) = dw(2,i,nfy+1) dw1(3,i,nfy+1) = dw(3,i,nfy+1) dw1(4,i,nfy+1) = dw(4,i,nfy+1) dw1(5,i,nfy+1) = dw(5,i,nfy+1) dw1(6,i,nfy+1) = dw(6,i,nfy+1) enddo call pvmfinitsend(PvmDataRaw,ibufid) call pvmfpack(IREALX,dw(1,2,nfy),(nfx-1)*6,1,info) call pvmfsend(jptask,1022,info) endif if (jmtask.ge.0) then call pvmfrecv(jmtask,1022,ibufid) call pvmfunpack(IREALX,dw1(1,2,1),(nfx-1)*6,1,info) endif if (imtask.ge.0) then call pvmfinitsend(PvmDataRaw,ibufid) call pvmfpack(IREALX,dw(1,2,2),nfy-1,imx*6,info) call pvmfpack(IREALX,dw(2,2,2),nfy-1,imx*6,info) call pvmfpack(IREALX,dw(3,2,2),nfy-1,imx*6,info) call pvmfpack(IREALX,dw(4,2,2),nfy-1,imx*6,info) call pvmfpack(IREALX,dw(5,2,2),nfy-1,imx*6,info) call pvmfpack(IREALX,dw(6,2,2),nfy-1,imx*6,info) call pvmfsend(imtask,1031,info) endif if (iptask.ge.0) then call pvmfrecv(iptask,1031,ibufid) call pvmfunpack(IREALX,dw(1,nfx+1,2),nfy-1,imx*6,info) call pvmfunpack(IREALX,dw(2,nfx+1,2),nfy-1,imx*6,info) call pvmfunpack(IREALX,dw(3,nfx+1,2),nfy-1,imx*6,info) call pvmfunpack(IREALX,dw(4,nfx+1,2),nfy-1,imx*6,info) call pvmfunpack(IREALX,dw(5,nfx+1,2),nfy-1,imx*6,info) call pvmfunpack(IREALX,dw(6,nfx+1,2),nfy-1,imx*6,info) do j = 2, nfy dw1(1,nfx+1,j) = dw(1,nfx+1,j) dw1(2,nfx+1,j) = dw(2,nfx+1,j) dw1(3,nfx+1,j) = dw(3,nfx+1,j) dw1(4,nfx+1,j) = dw(4,nfx+1,j) dw1(5,nfx+1,j) = dw(5,nfx+1,j) dw1(6,nfx+1,j) = dw(6,nfx+1,j) enddo call pvmfinitsend(PvmDataRaw,ibufid) call pvmfpack(IREALX,dw(1,nfx,2),nfy-1,imx*6,info) call pvmfpack(IREALX,dw(2,nfx,2),nfy-1,imx*6,info) call pvmfpack(IREALX,dw(3,nfx,2),nfy-1,imx*6,info) call pvmfpack(IREALX,dw(4,nfx,2),nfy-1,imx*6,info) call pvmfpack(IREALX,dw(5,nfx,2),nfy-1,imx*6,info) call pvmfpack(IREALX,dw(6,nfx,2),nfy-1,imx*6,info) call pvmfsend(iptask,1032,info) endif if (imtask.ge.0) then call pvmfrecv(imtask,1032,ibufid) call pvmfunpack(IREALX,dw1(1,1,2),nfy-1,imx*6,info) call pvmfunpack(IREALX,dw1(2,1,2),nfy-1,imx*6,info) call pvmfunpack(IREALX,dw1(3,1,2),nfy-1,imx*6,info) call pvmfunpack(IREALX,dw1(4,1,2),nfy-1,imx*6,info) call pvmfunpack(IREALX,dw1(5,1,2),nfy-1,imx*6,info) call pvmfunpack(IREALX,dw1(6,1,2),nfy-1,imx*6,info) endif return end subroutine upgrade include 'nsspmd2d.inc' do i = 2, nfx do j = 2, nfy w(1,i,j) = w(1,i,j) + dw(1,i,j) w(2,i,j) = w(2,i,j) + dw(2,i,j) w(3,i,j) = w(3,i,j) + dw(3,i,j) w(4,i,j) = w(4,i,j) + dw(4,i,j) w(5,i,j) = w(5,i,j) + dw(5,i,j)/w(1,i,j) w(6,i,j) = w(6,i,j) + dw(6,i,j)/w(1,i,j) if (w(5,i,j).lt.0) w(5,i,j) = . 0.25*(rknod(i,j)+rknod(i+1,j)+rknod(i,j+1)+rknod(i+1,j+1)) if (w(6,i,j).le.0) w(6,i,j) = . 0.25*(rwnod(i,j)+rwnod(i+1,j)+rwnod(i,j+1)+rwnod(i+1,j+1)) enddo enddo return end subroutine renew include 'nsspmd2d.inc' do j = nfy_start-1, 0 do i = 2,nfx p(i,j) = gam1*(w(4,i,j) . -0.5*(w(2,i,j)**2+w(3,i,j)**2)/w(1,i,j)) enddo enddo do j = nfy+2, nfy_end+2 do i = 2,nfx p(i,j) = gam1*(w(4,i,j) . -0.5*(w(2,i,j)**2+w(3,i,j)**2)/w(1,i,j)) enddo enddo do j = 2, nfy do i = nfx_start-1, 0 p(i,j) = gam1*(w(4,i,j) . -0.5*(w(2,i,j)**2+w(3,i,j)**2)/w(1,i,j)) enddo do i = nfx+2,nfx_end+2 p(i,j) = gam1*(w(4,i,j) . -0.5*(w(2,i,j)**2+w(3,i,j)**2)/w(1,i,j)) enddo enddo do j = 1, nfy+1 do i = 1, nfx+1 p(i,j) = gam1*(w(4,i,j) . -0.5*(w(2,i,j)**2+w(3,i,j)**2)/w(1,i,j)) c(i,j) = sqrt(gamma*p(i,j)/w(1,i,j)) enddo enddo return end subroutine bc include 'nsspmd2d.inc' include 'fpvm3.h' if (jmtask.ge.0) then call pvmfinitsend(PvmDataRaw,ibufid) call pvmfpack(IREALX,w(1,1,2),(nfx+1)*6,1,info) call pvmfpack(IREALX,w(1,1,3),(nfx+1)*6,1,info) call pvmfsend(jmtask,21,info) else call BCJmin endif if (jptask.ge.0) then call pvmfrecv(jptask,21,ibufid) call pvmfunpack(IREALX,w(1,1,nfy+1),(nfx+1)*6,1,info) call pvmfunpack(IREALX,w(1,1,nfy+2),(nfx+1)*6,1,info) call pvmfinitsend(PvmDataRaw,ibufid) call pvmfpack(IREALX,w(1,1,nfy-1),(nfx+1)*6,1,info) call pvmfpack(IREALX,w(1,1,nfy),(nfx+1)*6,1,info) call pvmfsend(jptask,22,info) else call BCJmax endif if (jmtask.ge.0) then call pvmfrecv(jmtask,22,ibufid) call pvmfunpack(IREALX,w(1,1,0),(nfx+1)*6,1,info) call pvmfunpack(IREALX,w(1,1,1),(nfx+1)*6,1,info) endif if (imtask.ge.0) then call pvmfinitsend(PvmDataRaw,ibufid) call pvmfpack(IREALX,w(1,2,1),nfy+1,(imx+1)*6,info) call pvmfpack(IREALX,w(2,2,1),nfy+1,(imx+1)*6,info) call pvmfpack(IREALX,w(3,2,1),nfy+1,(imx+1)*6,info) call pvmfpack(IREALX,w(4,2,1),nfy+1,(imx+1)*6,info) call pvmfpack(IREALX,w(5,2,1),nfy+1,(imx+1)*6,info) call pvmfpack(IREALX,w(6,2,1),nfy+1,(imx+1)*6,info) call pvmfpack(IREALX,w(1,3,1),nfy+1,(imx+1)*6,info) call pvmfpack(IREALX,w(2,3,1),nfy+1,(imx+1)*6,info) call pvmfpack(IREALX,w(3,3,1),nfy+1,(imx+1)*6,info) call pvmfpack(IREALX,w(4,3,1),nfy+1,(imx+1)*6,info) call pvmfpack(IREALX,w(5,3,1),nfy+1,(imx+1)*6,info) call pvmfpack(IREALX,w(6,3,1),nfy+1,(imx+1)*6,info) call pvmfsend(imtask,31,info) else call BCImin endif if (iptask.ge.0) then call pvmfrecv(iptask,31,ibufid) call pvmfunpack(IREALX,w(1,nfx+1,1),nfy+1,(imx+1)*6,info) call pvmfunpack(IREALX,w(2,nfx+1,1),nfy+1,(imx+1)*6,info) call pvmfunpack(IREALX,w(3,nfx+1,1),nfy+1,(imx+1)*6,info) call pvmfunpack(IREALX,w(4,nfx+1,1),nfy+1,(imx+1)*6,info) call pvmfunpack(IREALX,w(5,nfx+1,1),nfy+1,(imx+1)*6,info) call pvmfunpack(IREALX,w(6,nfx+1,1),nfy+1,(imx+1)*6,info) call pvmfunpack(IREALX,w(1,nfx+2,1),nfy+1,(imx+1)*6,info) call pvmfunpack(IREALX,w(2,nfx+2,1),nfy+1,(imx+1)*6,info) call pvmfunpack(IREALX,w(3,nfx+2,1),nfy+1,(imx+1)*6,info) call pvmfunpack(IREALX,w(4,nfx+2,1),nfy+1,(imx+1)*6,info) call pvmfunpack(IREALX,w(5,nfx+2,1),nfy+1,(imx+1)*6,info) call pvmfunpack(IREALX,w(6,nfx+2,1),nfy+1,(imx+1)*6,info) call pvmfinitsend(PvmDataRaw,ibufid) call pvmfpack(IREALX,w(1,nfx-1,1),nfy+1,(imx+1)*6,info) call pvmfpack(IREALX,w(2,nfx-1,1),nfy+1,(imx+1)*6,info) call pvmfpack(IREALX,w(3,nfx-1,1),nfy+1,(imx+1)*6,info) call pvmfpack(IREALX,w(4,nfx-1,1),nfy+1,(imx+1)*6,info) call pvmfpack(IREALX,w(5,nfx-1,1),nfy+1,(imx+1)*6,info) call pvmfpack(IREALX,w(6,nfx-1,1),nfy+1,(imx+1)*6,info) call pvmfpack(IREALX,w(1,nfx,1),nfy+1,(imx+1)*6,info) call pvmfpack(IREALX,w(2,nfx,1),nfy+1,(imx+1)*6,info) call pvmfpack(IREALX,w(3,nfx,1),nfy+1,(imx+1)*6,info) call pvmfpack(IREALX,w(4,nfx,1),nfy+1,(imx+1)*6,info) call pvmfpack(IREALX,w(5,nfx,1),nfy+1,(imx+1)*6,info) call pvmfpack(IREALX,w(6,nfx,1),nfy+1,(imx+1)*6,info) call pvmfsend(iptask,32,info) else call BCImax endif if (imtask.ge.0) then call pvmfrecv(imtask,32,ibufid) call pvmfunpack(IREALX,w(1,0,1),nfy+1,(imx+1)*6,info) call pvmfunpack(IREALX,w(2,0,1),nfy+1,(imx+1)*6,info) call pvmfunpack(IREALX,w(3,0,1),nfy+1,(imx+1)*6,info) call pvmfunpack(IREALX,w(4,0,1),nfy+1,(imx+1)*6,info) call pvmfunpack(IREALX,w(5,0,1),nfy+1,(imx+1)*6,info) call pvmfunpack(IREALX,w(6,0,1),nfy+1,(imx+1)*6,info) call pvmfunpack(IREALX,w(1,1,1),nfy+1,(imx+1)*6,info) call pvmfunpack(IREALX,w(2,1,1),nfy+1,(imx+1)*6,info) call pvmfunpack(IREALX,w(3,1,1),nfy+1,(imx+1)*6,info) call pvmfunpack(IREALX,w(4,1,1),nfy+1,(imx+1)*6,info) call pvmfunpack(IREALX,w(5,1,1),nfy+1,(imx+1)*6,info) call pvmfunpack(IREALX,w(6,1,1),nfy+1,(imx+1)*6,info) endif return end subroutine BCImin include 'nsspmd2d.inc' if (imtask.eq.RIEMANN) then call IminRiemann else if (imtask.eq.SYMMETRIC) then call IminSymmetric else if (imtask.eq.FAR) then call IminFar endif return end subroutine BCImax include 'nsspmd2d.inc' if (iptask.eq.RIEMANN) then call ImaxRiemann else if (iptask.eq.OUTFLOW) then call ImaxOutflow else if (iptask.eq.FAR) then call ImaxFar endif return end subroutine IminSymmetric include 'nsspmd2d.inc' do j = 2, nfy w(1,1,j) = w(1,2,j) w(2,1,j) = w(2,2,j) w(3,1,j) = -w(3,2,j) w(4,1,j) = w(4,2,j) w(5,1,j) = w(5,2,j) w(6,1,j) = w(6,2,j) vis(1,j) = vis(2,j) ff(1,j) = ff(2,j) enddo return end subroutine IminFar include 'nsspmd2d.inc' do j = 2, nfy w(1,1,j) = rho0 w(2,1,j) = rho0*u0 w(3,1,j) = rho0*v0 w(4,1,j) = rho0*e0 w(5,1,j) = rk0 w(6,1,j) = rw0 vis(1,j) = vis(2,j) ff(1,j) = ff(2,j) enddo return end subroutine ImaxOutflow include 'nsspmd2d.inc' do j = 2, nfy w(1,nfx+1,j) = w(1,nfx,j) w(2,nfx+1,j) = w(2,nfx,j) w(3,nfx+1,j) = w(3,nfx,j) w(4,nfx+1,j) = w(4,nfx,j) w(5,nfx+1,j) = w(5,nfx,j) w(6,nfx+1,j) = w(6,nfx,j) vis(nfx+1,j) = vis(nfx,j) ff(nfx+1,j) = ff(nfx,j) enddo return end subroutine ImaxFar include 'nsspmd2d.inc' do j = 2, nfy w(1,nfx+1,j) = rho0 w(2,nfx+1,j) = rho0*u0 w(3,nfx+1,j) = rho0*v0 w(4,nfx+1,j) = rho0*e0 w(5,nfx+1,j) = rk0 w(6,nfx+1,j) = rw0 vis(nfx+1,j) = vis(nfx,j) ff(nfx+1,j) = ff(nfx,j) enddo return end subroutine IminRiemann include 'nsspmd2d.inc' do j = 2, nfy jp = j + 1 cx = y(2,jp)-y(2,j) cy = x(2,j)-x(2,jp) dd = sqrt(cx**2+cy**2) cx = cx/dd cy = cy/dd qne = (w(2,2,j)*cx+w(3,2,j)*cy)/w(1,2,j) qte = (-w(2,2,j)*cy+w(3,2,j)*cx)/w(1,2,j) qnf = u0*cx+v0*cy qtf = -u0*cy+v0*cx R2 = qnf+2*c0/gam1 R1 = qne-2*c(2,j)/gam1 qn = 0.5*(R1+R2) c(1,j) = 0.25*gam1*(R2-R1) if (qn.ge.0) then qt = qtf w(1,1,j) = . (rho0**gamma*c(1,j)**2/p0/gamma)**(1./gam1) w(5,1,j) = rk0 w(6,1,j) = rw0 else qt = qte w(1,1,j) = . (w(1,2,j)**gamma*c(1,j)**2/p(2,j)/gamma)**(1./gam1) w(5,1,j) = w(5,2,j) w(6,1,j) = w(6,2,j) endif p(1,j) = w(1,1,j)*c(1,j)**2/gamma w(2,1,j) = w(1,1,j)*(qn*cx - qt*cy) w(3,1,j) = w(1,1,j)*(qn*cy + qt*cx) w(4,1,j) = p(1,j)/gam1 . +0.5*(w(2,1,j)**2+w(3,1,j)**2)/w(1,1,j) vis(1,j) = vis(2,j) ff(1,j) = ff(2,j) enddo return end subroutine ImaxRiemann include 'nsspmd2d.inc' do j = 2, nfy jp = j + 1 cx = y(nfx+1,jp)-y(nfx+1,j) cy = x(nfx+1,j)-x(nfx+1,jp) dd = sqrt(cx**2+cy**2) cx = cx/dd cy = cy/dd qne = (w(2,nfx,j)*cx+w(3,nfx,j)*cy)/w(1,nfx,j) qte = (-w(2,nfx,j)*cy+w(3,nfx,j)*cx)/w(1,nfx,j) qnf = u0*cx+v0*cy qtf = -u0*cy+v0*cx R2 = qne+2*c(nfx,j)/gam1 R1 = qnf-2*c0/gam1 qn = 0.5*(R1+R2) c(nfx+1,j) = 0.25*gam1*(R2-R1) if (qn.le.0) then qt = qtf w(1,nfx+1,j) = . (rho0**gamma*c(nfx+1,j)**2/p0/gamma)**(1./gam1) w(5,nfx+1,j) = rk0 w(6,nfx+1,j) = rw0 else qt = qte w(1,nfx+1,j) = . (w(1,nfx,j)**gamma*c(nfx+1,j)**2/p(nfx,j)/gamma)**(1./gam1) w(5,nfx+1,j) = w(5,nfx,j) w(6,nfx+1,j) = w(6,nfx,j) endif p(nfx+1,j) = w(1,nfx+1,j)*c(nfx+1,j)**2/gamma w(2,nfx+1,j) = w(1,nfx+1,j)*(qn*cx - qt*cy) w(3,nfx+1,j) = w(1,nfx+1,j)*(qn*cy + qt*cx) w(4,nfx+1,j) = p(nfx+1,j)/gam1 . +0.5*(w(2,nfx+1,j)**2+w(3,nfx+1,j)**2)/w(1,nfx+1,j) vis(nfx+1,j) = vis(nfx,j) ff(nfx+1,j) = ff(nfx,j) enddo return end subroutine BCJmin include 'nsspmd2d.inc' if (jmtask.eq.RIEMANN) then call JminRiemann else if (jmtask.eq.WALL) then call JminWall else if (jmtask.eq.FAR) then call JminFar else if (jmtask.eq.OUTFLOW) then call JminOutflow else if (jmtask.eq.SYMMETRIC) then call JminSymmetric endif return end subroutine BCJmax include 'nsspmd2d.inc' if (jptask.eq.RIEMANN) then call JmaxRiemann else if (jptask.eq.WALL) then call JmaxWall else if (jptask.eq.FAR) then call JmaxFar else if (jptask.eq.OUTFLOW) then call JmaxOutflow else if (jptask.eq.SYMMETRIC) then call JmaxSymmetric endif return end subroutine JminSymmetric include 'nsspmd2d.inc' do i = 2, nfx w(1,i,1) = w(1,i,2) w(2,i,1) = w(2,i,2) w(3,i,1) = -w(3,i,2) w(4,i,1) = w(4,i,2) w(5,i,1) = w(5,i,2) w(6,i,1) = w(6,i,2) vis(i,1) = vis(i,2) ff(i,1) = ff(i,2) enddo return end subroutine JminWall include 'nsspmd2d.inc' do i = 2, nfx w(1,i,1) = w(1,i,2) w(2,i,1) = -w(2,i,2) w(3,i,1) = -w(3,i,2) w(4,i,1) = w(4,i,2) xx = 0.25*(x(i,2)+x(i+1,2)+x(i,3)+x(i+1,3)) yy = 0.25*(y(i,2)+y(i+1,2)+y(i,3)+y(i+1,3)) xc = 0.5*(x(i,2)+x(i+1,2)) yc = 0.5*(y(i,2)+y(i+1,2)) dyn = (xx-xc)**2+(yy-yc)**2 c ut = (w(2,i,2)**2+w(3,i,2)**2)/w(1,i,2)**2 c owall = 2500*sqrt(ut/dyn) ccc = sqrt(gamma*p(i,2)/w(1,i,2)) t = ccc**2 rmul = t**1.5*(1.+tcoe)/(t+tcoe) owall = 60*rmul/(beta1*dyn*ra) w(5,i,1) = -w(5,i,2) w(6,i,1) = -w(6,i,2)+2*owall vis(i,1) = -vis(i,2) ff(i,1) = ff(i,2) enddo return end subroutine JminFar include 'nsspmd2d.inc' do i = 2, nfx w(1,i,1) = rho0 w(2,i,1) = rho0*u0 w(3,i,1) = rho0*v0 w(4,i,1) = rho0*e0 w(5,i,1) = rk0 w(6,i,1) = rw0 vis(i,1) = vis(i,2) ff(i,1) = ff(i,2) enddo return end subroutine JminOutflow include 'nsspmd2d.inc' do i = 2, nfx w(1,i,1) = w(1,i,2) w(2,i,1) = w(2,i,2) w(3,i,1) = w(3,i,2) w(4,i,1) = w(4,i,2) w(5,i,1) = w(5,i,2) w(6,i,1) = w(6,i,2) vis(i,1) = vis(i,2) ff(i,1) = ff(i,2) enddo return end subroutine JmaxSymmetric include 'nsspmd2d.inc' do i = 2, nfx w(1,i,nfy+1) = w(1,i,nfy) w(2,i,nfy+1) = w(2,i,nfy) w(3,i,nfy+1) = -w(3,i,nfy) w(4,i,nfy+1) = w(4,i,nfy) w(5,i,nfy+1) = w(5,i,nfy) w(6,i,nfy+1) = w(6,i,nfy) vis(i,nfy+1) = vis(i,nfy) ff(i,nfy+1) = ff(i,nfy) enddo return end subroutine JminRiemann include 'nsspmd2d.inc' do i = 2, nfx ip = i + 1 ex = y(i,2)-y(ip,2) ey = x(ip,2)-x(i,2) dd = sqrt(ex**2+ey**2) ex = ex/dd ey = ey/dd qne = (w(2,i,2)*ex+w(3,i,2)*ey)/w(1,i,2) qte = (-w(2,i,2)*ey+w(3,i,2)*ex)/w(1,i,2) qnf = u0*ex+v0*ey qtf = -u0*ey+v0*ex R2 = qnf+2*c0/gam1 R1 = qne-2*c(i,2)/gam1 qn = 0.5*(R1+R2) c(i,1) = 0.25*gam1*(R2-R1) if (qn.ge.0) then qt = qtf w(1,i,1) = . (rho0**gamma*c(i,1)**2/p0/gamma)**(1./gam1) w(5,i,1) = rk0 w(6,i,1) = rw0 else qt = qte w(1,i,1) = . (w(1,i,2)**gamma*c(i,1)**2/p(i,2)/gamma)**(1./gam1) w(5,i,1) = w(5,i,2) w(6,i,1) = w(6,i,2) endif p(i,1) = w(1,i,1)*c(i,1)**2/gamma w(2,i,1) = w(1,i,1)*(qn*ex - qt*ey) w(3,i,1) = w(1,i,1)*(qn*ey + qt*ex) w(4,i,1) = p(i,1)/gam1 . +0.5*(w(2,i,1)**2+w(3,i,1)**2)/w(1,i,1) vis(i,1) = vis(i,2) ff(i,1) = ff(i,2) enddo return end subroutine JmaxRiemann include 'nsspmd2d.inc' do i = 2, nfx ip = i + 1 ex = y(i,nfy+1)-y(ip,nfy+1) ey = x(ip,nfy+1)-x(i,nfy+1) dd = sqrt(ex**2+ey**2) ex = ex/dd ey = ey/dd qne = (w(2,i,nfy)*ex+w(3,i,nfy)*ey)/w(1,i,nfy) qte = (-w(2,i,nfy)*ey+w(3,i,nfy)*ex)/w(1,i,nfy) qnf = u0*ex+v0*ey qtf = -u0*ey+v0*ex R2 = qne+2*c(i,nfy)/gam1 R1 = qnf-2*c0/gam1 qn = 0.5*(R1+R2) c(i,nfy+1) = 0.25*gam1*(R2-R1) if (qn.le.0) then qt = qtf w(1,i,nfy+1) = . (rho0**gamma*c(i,nfy+1)**2/p0/gamma)**(1./gam1) w(5,i,nfy+1) = rk0 w(6,i,nfy+1) = rw0 else qt = qte w(1,i,nfy+1) = . (w(1,i,nfy)**gamma*c(i,nfy+1)**2/p(i,nfy)/gamma)**(1./gam1) w(5,i,nfy+1) = w(5,i,nfy) w(6,i,nfy+1) = w(6,i,nfy) endif p(i,nfy+1) = w(1,i,nfy+1)*c(i,nfy+1)**2/gamma w(2,i,nfy+1) = w(1,i,nfy+1)*(qn*ex - qt*ey) w(3,i,nfy+1) = w(1,i,nfy+1)*(qn*ey + qt*ex) w(4,i,nfy+1) = p(i,nfy+1)/gam1 . +0.5*(w(2,i,nfy+1)**2+w(3,i,nfy+1)**2)/w(1,i,nfy+1) vis(i,nfy+1) = vis(i,nfy) ff(i,nfy+1) = ff(i,nfy) enddo return end subroutine JmaxWall include 'nsspmd2d.inc' do i = 2, nfx w(1,i,nfy+1) = w(1,i,nfy) w(2,i,nfy+1) = -w(2,i,nfy) w(3,i,nfy+1) = -w(3,i,nfy) w(4,i,nfy+1) = w(4,i,nfy) xx = 0.25*(x(i,nfy+1)+x(i+1,nfy+1)+x(i,nfy)+x(i+1,nfy)) yy = 0.25*(y(i,nfy+1)+y(i+1,nfy+1)+y(i,nfy)+y(i+1,nfy)) xc = 0.5*(x(i,nfy+1)+x(i+1,nfy+1)) yc = 0.5*(y(i,nfy+1)+y(i+1,nfy+1)) dyn = (xx-xc)**2+(yy-yc)**2 c ut = (w(2,i,nfy)**2+w(3,i,nfy)**2)/w(1,i,nfy)**2 c owall = 2500*sqrt(ut/dyn) ccc = sqrt(gamma*p(i,nfy)/w(1,i,nfy)) t = ccc**2 rmul = t**1.5*(1.+tcoe)/(t+tcoe) owall = 60*rmul/(beta1*dyn*ra) w(5,i,nfy+1) = -w(5,i,nfy) w(6,i,nfy+1) = -w(6,i,nfy)+2*owall vis(i,nfy+1) = -vis(i,nfy) ff(i,nfy+1) = ff(i,nfy) enddo return end subroutine JmaxFar include 'nsspmd2d.inc' do i = 2, nfx w(1,i,nfy+1) = rho0 w(2,i,nfy+1) = rho0*u0 w(3,i,nfy+1) = rho0*v0 w(4,i,nfy+1) = rho0*e0 w(5,i,nfy+1) = rk0 w(6,i,nfy+1) = rw0 vis(i,nfy+1) = vis(i,nfy) ff(i,nfy+1) = ff(i,nfy) enddo return end subroutine JmaxOutflow include 'nsspmd2d.inc' do i = 2, nfx w(1,i,nfy+1) = w(1,i,nfy) w(2,i,nfy+1) = w(2,i,nfy) w(3,i,nfy+1) = w(3,i,nfy) w(4,i,nfy+1) = w(4,i,nfy) w(5,i,nfy+1) = w(5,i,nfy) w(6,i,nfy+1) = w(6,i,nfy) vis(i,nfy+1) = vis(i,nfy) ff(i,nfy+1) = ff(i,nfy) enddo return end subroutine hole_check include 'nsspmd2d.inc' include 'fpvm3.h' if (jmtask.ge.0) then call pvmfinitsend(PvmDataRaw,ibufid) call pvmfpack(INTEGERX,ib(2,2),nfx-1,1,info) call pvmfpack(IREALX,x(2,3),nfx,1,info) call pvmfpack(IREALX,y(2,3),nfx,1,info) call pvmfsend(jmtask,221,info) else do i = 2, nfx if (ib(i,2).eq.0) ib(i,1) = 0 enddo endif if (jptask.ge.0) then call pvmfrecv(jptask,221,ibufid) call pvmfunpack(INTEGERX,ib(2,nfy+1),nfx-1,1,info) call pvmfunpack(IREALX,x(2,nfy+2),nfx,1,info) call pvmfunpack(IREALX,y(2,nfy+2),nfx,1,info) call pvmfinitsend(PvmDataRaw,ibufid) call pvmfpack(INTEGERX,ib(2,nfy),nfx-1,1,info) call pvmfpack(IREALX,x(2,nfy),nfx,1,info) call pvmfpack(IREALX,y(2,nfy),nfx,1,info) call pvmfsend(jptask,222,info) else do i = 2, nfx if (ib(i,nfy).eq.0) ib(i,nfy+1) = 0 enddo endif if (jmtask.ge.0) then call pvmfrecv(jmtask,222,ibufid) call pvmfunpack(INTEGERX,ib(2,1),nfx-1,1,info) call pvmfunpack(IREALX,x(2,1),nfx,1,info) call pvmfunpack(IREALX,y(2,1),nfx,1,info) endif if (imtask.ge.0) then call pvmfinitsend(PvmDataRaw,ibufid) call pvmfpack(INTEGERX,ib(2,2),nfy-1,imx,info) call pvmfpack(IREALX,x(3,2),nfy,imx,info) call pvmfpack(IREALX,y(3,2),nfy,imx,info) call pvmfsend(imtask,231,info) else do j = 2, nfy if (ib(2,j).eq.0) ib(1,j) = 0 enddo endif if (iptask.ge.0) then call pvmfrecv(iptask,231,ibufid) call pvmfunpack(INTEGERX,ib(nfx+1,2),nfy-1,imx,info) call pvmfunpack(IREALX,x(nfx+2,2),nfy,imx,info) call pvmfunpack(IREALX,y(nfx+2,2),nfy,imx,info) call pvmfinitsend(PvmDataRaw,ibufid) call pvmfpack(INTEGERX,ib(nfx,2),nfy-1,imx,info) call pvmfpack(IREALX,x(nfx,2),nfy,imx,info) call pvmfpack(IREALX,y(nfx,2),nfy,imx,info) call pvmfsend(iptask,232,info) else do j = 2, nfy if (ib(nfx,j).eq.0) ib(nfx+1,j) = 0 enddo endif if (imtask.ge.0) then call pvmfrecv(imtask,232,ibufid) call pvmfunpack(INTEGERX,ib(1,2),nfy-1,imx,info) call pvmfunpack(IREALX,x(1,2),nfy,imx,info) call pvmfunpack(IREALX,y(1,2),nfy,imx,info) endif return end subroutine muscl(qll,ql,qr,qrr,ib1,ib2) capa = 1./3. dql = ql - qll dqc = qr - ql dqr = qrr - qr ss = (2.*dql*dqc+1.e-6)/(dqc**2+dql**2+1.e-6) omks = 1.0 - capa*ss opks = 1.0 + capa*ss ql = ql + 0.25*ss*(omks*dql+opks*dqc)*ib1*ib2 ss = (2.*dqc*dqr+1.e-6)/(dqr**2+dqc**2+1.e-6) omks = 1.0 - capa*ss opks = 1.0 + capa*ss qr = qr - 0.25*ss*(omks*dqr+opks*dqc)*ib1*ib2 return end subroutine grid include 'nsspmd2d.inc' dimension yt(imx,jmx) read(13,*) nfxc, nfyc, i_hole if(nx.ne.nfxc.or.ny.ne.nfyc) then call pvmfexit(ierror) endif do n = 1, i_hole read(13,*) i_hole_start(n), j_hole_start(n), . i_hole_end(n), j_hole_end(n) enddo do i = 2, nx+1 do j = 2, ny+1 read(13,*) gx(i,j), gy(i,j) enddo enddo call sstgrid(yt) do j = 2, ny do i = 2, nx yyn(i,j) = 0.25*(yt(i,j)+yt(i+1,j)+yt(i,j+1)+yt(i+1,j+1)) enddo enddo do j = 1, ny+1 do i = 1, nx+1 nb(i,j) = 1 enddo enddo do n = 1, i_hole do i = i_hole_start(n)+1, i_hole_end(n) do j = j_hole_start(n)+1, j_hole_end(n) nb(i,j) = 0 enddo enddo enddo return end subroutine F1calc include 'nsspmd2d.inc' do j = 2, nfy jp = j + 1 jm = j - 1 do i = 2, nfx ip = i + 1 im = i - 1 uch = 0.5*(w(2,ip,j)/w(1,ip,j)-w(2,im,j)/w(1,im,j)) uet = 0.5*(w(2,i,jp)/w(1,i,jp)-w(2,i,jm)/w(1,i,jm)) vch = 0.5*(w(3,ip,j)/w(1,ip,j)-w(3,im,j)/w(1,im,j)) vet = 0.5*(w(3,i,jp)/w(1,i,jp)-w(3,i,jm)/w(1,i,jm)) rkch = 0.5*(w(5,ip,j)-w(5,im,j)) rket = 0.5*(w(5,i,jp)-w(5,i,jm)) rwch = 0.5*(w(6,ip,j)-w(6,im,j)) rwet = 0.5*(w(6,i,jp)-w(6,i,jm)) uy = uch*chy(i,j)+uet*ety(i,j) vx = vch*chx(i,j)+vet*etx(i,j) rkx = rkch*chx(i,j)+rket*etx(i,j) rky = rkch*chy(i,j)+rket*ety(i,j) rwx = rwch*chx(i,j)+rwet*etx(i,j) rwy = rwch*chy(i,j)+rwet*ety(i,j) CDkw = max(2*w(1,i,j)*sigmao2*(rkx*rwx+rky*rwy) . /w(6,i,j),1.e-20) term1 = sqrt(w(5,i,j))/(0.09*w(6,i,j)*yn(i,j)) t = c(i,j)**2 rmu = t**1.5*(1.+tcoe)/(t+tcoe) term2 = 500.*rmu/(yn(i,j)**2*w(6,i,j)*ra*w(1,i,j)) term3 = 4*w(1,i,j)*sigmao2*w(5,i,j)/(CDkw*yn(i,j)**2) arg1 = min(max(term1,term2),term3) ff(i,j) = tanh(arg1**4) arg2 = max(2*term1,term2) f2 = tanh(arg2**2) vor = abs(uy-vx) vis(i,j) = aa1*w(5,i,j)*w(1,i,j)/max(aa1*w(6,i,j),vor*f2) enddo enddo return end subroutine visbc include 'nsspmd2d.inc' include 'fpvm3.h' if (jmtask.ge.0) then call pvmfinitsend(PvmDataRaw,ibufid) call pvmfpack(IREALX,vis(1,2),nfx+1,1,info) call pvmfpack(IREALX,ff(1,2),nfx+1,1,info) call pvmfsend(jmtask,921,info) endif if (jptask.ge.0) then call pvmfrecv(jptask,921,ibufid) call pvmfunpack(IREALX,vis(1,nfy+1),nfx+1,1,info) call pvmfunpack(IREALX,ff(1,nfy+1),nfx+1,1,info) call pvmfinitsend(PvmDataRaw,ibufid) call pvmfpack(IREALX,vis(1,nfy),nfx+1,1,info) call pvmfpack(IREALX,ff(1,nfy),nfx+1,1,info) call pvmfsend(jptask,922,info) endif if (jmtask.ge.0) then call pvmfrecv(jmtask,922,ibufid) call pvmfunpack(IREALX,vis(1,1),nfx+1,1,info) call pvmfunpack(IREALX,ff(1,1),nfx+1,1,info) endif if (imtask.ge.0) then call pvmfinitsend(PvmDataRaw,ibufid) call pvmfpack(IREALX,vis(2,1),nfy+1,imx,info) call pvmfpack(IREALX,ff(2,1),nfy+1,imx,info) call pvmfsend(imtask,931,info) endif if (iptask.ge.0) then call pvmfrecv(iptask,931,ibufid) call pvmfunpack(IREALX,vis(nfx+1,1),nfy+1,imx,info) call pvmfunpack(IREALX,ff(nfx+1,1),nfy+1,imx,info) call pvmfinitsend(PvmDataRaw,ibufid) call pvmfpack(IREALX,vis(nfx,1),nfy+1,imx,info) call pvmfpack(IREALX,ff(nfx,1),nfy+1,imx,info) call pvmfsend(iptask,932,info) endif if (imtask.ge.0) then call pvmfrecv(imtask,932,ibufid) call pvmfunpack(IREALX,vis(1,1),nfy+1,imx,info) call pvmfunpack(IREALX,ff(1,1),nfy+1,imx,info) endif return end subroutine sstgrid(yt) include 'nsspmd2d.inc' dimension yt(imx,jmx),mb(imx,jmx) do j = 2, ny+1 do i = 2, nx+1 mb(i,j) = 1 enddo enddo if (iminbc.eq.WALL) then do j = 2, ny+1 mb(2,j) = 0 enddo endif if (imaxbc.eq.WALL) then do j = 2, ny+1 mb(nx+1,j) = 0 enddo endif if (jminbc.eq.WALL) then do i = 2, nx+1 mb(i,2) = 0 enddo endif if (jmaxbc.eq.WALL) then do i = 2, nx+1 mb(i,ny+1) = 0 enddo endif do n = 1, i_hole do i = i_hole_start(n)+2,i_hole_end(n) do j = j_hole_start(n)+2,j_hole_end(n) mb(i,j) = -1 enddo enddo do j = j_hole_start(n)+1,j_hole_end(n)+1 mb(i_hole_start(n)+1,j) = 0 mb(i_hole_end(n)+1,j) = 0 enddo do i = i_hole_start(n)+1,i_hole_end(n)+1 mb(i,j_hole_start(n)+1) = 0 mb(i,j_hole_end(n)+1) = 0 enddo enddo do j = 2, ny+1 if (mb(3,j).eq.-1) mb(2,j) = -1 if (mb(nx,j).eq.-1) mb(nx+1,j) = -1 enddo do i = 2, nx+1 if (mb(i,3).eq.-1) mb(i,2) = -1 if (mb(i,ny).eq.-1) mb(i,ny+1) = -1 enddo do i = 3, nx im = i - 1 ip = i + 1 do j = 2, ny+1 if ((mb(i,j).eq.0).and.(mb(im,j).eq.-1).and.(mb(ip,j).eq.-1)) . mb(i,j) = -1 enddo enddo do j = 3, ny jm = j - 1 jp = j + 1 do i = 2, nx+1 if ((mb(i,j).eq.0).and.(mb(i,jm).eq.-1).and.(mb(i,jp).eq.-1)) . mb(i,j) = -1 enddo enddo do j = 2, ny+1 do i = 2, nx+1 yt(i,j) = 1.e10 do n = 2, ny+1 do m = 2, nx+1 if ((mb(i,j).eq.1).and.(mb(m,n).eq.0)) then ytmp = sqrt((gx(i,j)-gx(m,n))**2+(gy(i,j)-gy(m,n))**2) if (ytmp.le.yt(i,j)) yt(i,j) = ytmp endif if (mb(i,j).eq.0) yt(i,j) = 0.0 if (mb(i,j).eq.-1) yt(i,j) = -1 enddo enddo enddo enddo return end --------------BB2C3A587E7F9B68300F0821 Content-Type: text/plain; charset=EUC-KR; name="ns2d-mpi.f" Content-Transfer-Encoding: 7bit Content-Disposition: inline; filename="ns2d-mpi.f" program mpins2d include 'mpins2d.inc' include 'mpif.h' integer pe(maxpe,maxpe),ibuf(12) real rbuf(8) call MPI_INIT(ierr) call MPI_COMM_SIZE(MPI_COMM_WORLD,npes,ierr) call MPI_COMM_RANK(MPI_COMM_WORLD,myPE,ierr) call MPI_TYPE_CONTIGUOUS(6,MPI_REAL,MPI_VEC,ierr) call MPI_TYPE_COMMIT(MPI_VEC,ierr) if (myPE.eq.0) then open(12,file='init.inp') open(13,file='grid.grd') read(12,500) read(12,*) alpa, re, amach, tinf read(12,500) read(12,*) nx, ny, tol, itmax, isubmax read(12,500) read(12,*) ifread, delt, intwrt, cfl read(12,500) read(12,*) cappa, eee, isweep, init read(12,500) read(12,*) iminbc, imaxbc, jminbc, jmaxbc read(12,500) read(12,*) i_proc, j_proc nproc = i_proc*j_proc do i = 1, i_proc do j = 1, j_proc pe(i,j) = (i-1)*j_proc+j-1 enddo enddo open(21,file='result.dat') open(22,file='wnstep.dat') open(23,file='wn1step.dat') open(30,file='error.dat') call GRID if (ifread.lt.0) then ptime = 0.0 itbegin = 1 else read(22,*) ptime, itbegin, errbegin do i = 1, nx+1 do j = 1, ny+1 read(22,*) qq(1,i,j),qq(2,i,j),qq(3,i,j), . qq(4,i,j),qq(5,i,j),qq(6,i,j) read(23,*) qm(1,i,j),qm(2,i,j),qm(3,i,j), . qm(4,i,j),qm(5,i,j),qm(6,i,j) enddo enddo endif do i = 1, i_proc do j = 1, j_proc if (i.eq.1) then if (iminbc.eq.PERIODIC) then imPE = pe(i_proc,j) else imPE = iminbc endif else imPE = pe(i-1,j) endif if (i.eq.i_proc) then if (imaxbc.eq.PERIODIC) then ipPE = pe(1,j) else ipPE = imaxbc endif else ipPE = pe(i+1,j) endif if (j.eq.1) then jmPE = jminbc else jmPE = pe(i,j-1) endif if (j.eq.j_proc) then jpPE = jmaxbc else jpPE = pe(i,j+1) endif i_cell = (nx-1)/i_proc j_cell = (ny-1)/j_proc i_remain = (nx-1)-(nx-1)/i_proc*i_proc j_remain = (ny-1)-(ny-1)/j_proc*j_proc if (i.le.i_remain) i_cell = i_cell+1 if (j.le.j_remain) j_cell = j_cell+1 i_end = i*i_cell+1 j_end = j*j_cell+1 if (i.gt.i_remain) i_end = i_end+i_remain if (j.gt.j_remain) j_end = j_end+j_remain i_begin = i_end-i_cell j_begin = j_end-j_cell nfx_slave = i_cell+1 nfy_slave = j_cell+1 ibuf(1) = imPE ibuf(2) = ipPE ibuf(3) = jmPE ibuf(4) = jpPE ibuf(5) = nfx_slave ibuf(6) = nfy_slave ibuf(7) = itbegin ibuf(8) = itmax ibuf(9) = intwrt ibuf(10) = isweep ibuf(11) = ifread ibuf(12) = init rbuf(1) = alpa rbuf(2) = re rbuf(3) = amach rbuf(4) = tinf rbuf(5) = cfl rbuf(6) = delt rbuf(7) = cappa rbuf(8) = eee call MPI_SEND(ibuf,12,MPI_INTEGER,pe(i,j),100, . MPI_COMM_WORLD,ierr) call MPI_SEND(rbuf,8,MPI_REAL,pe(i,j),200, . MPI_COMM_WORLD,ierr) do jj = j_begin+1,j_end+1 call MPI_SEND(gx(i_begin+1,jj),i_end-i_begin+1, . MPI_REAL,pe(i,j),11,MPI_COMM_WORLD,ierr) call MPI_SEND(gy(i_begin+1,jj),i_end-i_begin+1, . MPI_REAL,pe(i,j),22,MPI_COMM_WORLD,ierr) enddo if (ifread.ge.0) then do jj = j_begin,j_end+1 call MPI_SEND(qq(1,i_begin,jj),i_end-i_begin+2, . MPI_VEC,pe(i,j),33,MPI_COMM_WORLD,ierr) call MPI_SEND(qm(1,i_begin,jj),i_end-i_begin+2, . MPI_VEC,pe(i,j),44,MPI_COMM_WORLD,ierr) enddo endif do jj = j_begin+1,j_end call MPI_SEND(nb(i_begin+1,jj),i_end-i_begin, . MPI_INTEGER,pe(i,j),55,MPI_COMM_WORLD,ierr) call MPI_SEND(yyn(i_begin+1,jj),i_end-i_begin, . MPI_REAL,pe(i,j),66,MPI_COMM_WORLD,ierr) enddo enddo enddo endif call MPI_RECV(ibuf,12,MPI_INTEGER,0,100, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(rbuf,8,MPI_REAL,0,200, . MPI_COMM_WORLD,istatus,ierr) imPE = ibuf(1) ipPE = ibuf(2) jmPE = ibuf(3) jpPE = ibuf(4) nfx = ibuf(5) nfy = ibuf(6) itbegin = ibuf(7) itmax = ibuf(8) intwrt = ibuf(9) isweep = ibuf(10) ifread = ibuf(11) init = ibuf(12) alpa = rbuf(1) re = rbuf(2) amach = rbuf(3) tinf = rbuf(4) cfl = rbuf(5) delt = rbuf(6) cappa = rbuf(7) eee = rbuf(8) do j = 2, nfy+1 call MPI_RECV(x(2,j),nfx,MPI_REAL,0,11, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(y(2,j),nfx,MPI_REAL,0,22, . MPI_COMM_WORLD,istatus,ierr) enddo pi = 4.0*atan(1.0) alpa = alpa*pi/180.0 rho0 = 1.0 u0 = amach*cos(alpa) v0 = amach*sin(alpa) e0 = 1./gamma/gam1+0.5*amach**2 p0 = 1.0/gamma c0 = sqrt(gamma*p0/rho0) ra = re/amach tcoe = 110.4/tinf rk0 = 1.e-7 rw0 = 1.0 tgamma1 = beta1/betaA-sigmao1*0.41**2/sqrt(betaA) tgamma2 = beta2/betaA-sigmao2*0.41**2/sqrt(betaA) nfx_start = 1 nfx_end = nfx nfy_start = 1 nfy_end = nfy if (imPE.lt.0) nfx_start = 2 if (ipPE.lt.0) nfx_end = nfx-1 if (jmPE.lt.0) nfy_start = 2 if (jpPE.lt.0) nfy_end = nfy-1 do j = 1, nfy+1 do i = 1, nfx+1 ib(i,j) = 1 zac(i,j) = 1.0 ff(i,j) = 1.0 enddo enddo if (ifread.ge.0) then do j = 1, nfy+1 call MPI_RECV(w_n(1,1,j),nfx+1,MPI_VEC,0,33, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(w_nm(1,1,j),nfx+1,MPI_VEC,0,44, . MPI_COMM_WORLD,istatus,ierr) do i = 1, nfx+1 w(1,i,j) = w_n(1,i,j) w(2,i,j) = w_n(2,i,j) w(3,i,j) = w_n(3,i,j) w(4,i,j) = w_n(4,i,j) w(5,i,j) = w_n(5,i,j)/w(1,i,j) w(6,i,j) = w_n(6,i,j)/w(1,i,j) p(i,j) = gam1*(w(4,i,j) . -0.5*(w(2,i,j)**2+w(3,i,j)**2)/w(1,i,j)) c(i,j) = sqrt(gamma*p(i,j)/w(1,i,j)) enddo enddo else if (init.le.0) then do j = 1, nfy+1 do i = 1, nfx+1 w(2,i,j) = 0 w(3,i,j) = 0 enddo enddo else if (init.eq.1) then do j = 1, nfy+1 do i = 1, nfx+1 w(2,i,j) = rho0*u0 w(3,i,j) = rho0*v0 enddo enddo else if (init.eq.2) then do j = 1, nfy+1 do i = 1, nfx+1 if (y(i,j).ge.0) then esp = -0.01 else esp = +0.01 endif w(2,i,j) = rho0*u0*(1+esp) w(3,i,j) = rho0*v0 enddo enddo endif do j = 1, nfy+1 do i = 1, nfx+1 w(1,i,j) = rho0 w(4,i,j) = rho0*e0 w(5,i,j) = rk0 w(6,i,j) = rw0 p(i,j) = p0 c(i,j) = c0 w_n(1,i,j) = w(1,i,j) w_n(2,i,j) = w(2,i,j) w_n(3,i,j) = w(3,i,j) w_n(4,i,j) = w(4,i,j) w_n(5,i,j) = w(5,i,j)*w(1,i,j) w_n(6,i,j) = w(6,i,j)*w(1,i,j) w_nm(1,i,j) = w(1,i,j) w_nm(2,i,j) = w(2,i,j) w_nm(3,i,j) = w(3,i,j) w_nm(4,i,j) = w(4,i,j) w_nm(5,i,j) = w(5,i,j)*w(1,i,j) w_nm(6,i,j) = w(6,i,j)*w(1,i,j) enddo enddo endif do j = 2, nfy call MPI_RECV(ib(2,j),nfx-1,MPI_INTEGER,0,55, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(yn(2,j),nfx-1,MPI_REAL,0,66, . MPI_COMM_WORLD,istatus,ierr) enddo call HOLE_CHECK call METRIC call METDUM call BC call RENEW call F1CALC call VISBC i_total = 0 iconv = 0 do it = itbegin, itmax itsub = 0 ptime = ptime + delt call STEP call NODE call PHYSICS call SOLVE call UPGRADE call BC call RENEW call F1CALC call VISBC do while (iconv.eq.0) itsub = itsub + 1 call STEP call NODE call PHYSICS call SOLVE call UPGRADE call BC call RENEW call F1CALC call VISBC err = 0.0 do j = 2, nfy do i = 2, nfx err = err + dw(1,i,j)**2 enddo enddo call MPI_REDUCE(err,err,1,MPI_REAL,MPI_SUM,0, . MPI_COMM_WORLD,ierr) call MPI_REDUCE(dtmin,dtmin,1,MPI_REAL,MPI_MIN,0, . MPI_COMM_WORLD,ierr) call MPI_REDUCE(dtmax,dtmax,1,MPI_REAL,MPI_MAX,0, . MPI_COMM_WORLD,ierr) call FORCE(rLp,rDp,rLf,rDf) call MPI_REDUCE(rLp,rLp,1,MPI_REAL,MPI_SUM,0, . MPI_COMM_WORLD,ierr) call MPI_REDUCE(rDp,rDp,1,MPI_REAL,MPI_SUM,0, . MPI_COMM_WORLD,ierr) call MPI_REDUCE(rLf,rLf,1,MPI_REAL,MPI_SUM,0, . MPI_COMM_WORLD,ierr) call MPI_REDUCE(rDf,rDf,1,MPI_REAL,MPI_SUM,0, . MPI_COMM_WORLD,ierr) if (myPE.eq.0) then i_total = i_total + 1 errtotal = log10(sqrt(err/((nx-1)*(ny-1)))) if (it.eq.1.and.i_total.eq.1) errbegin = errtotal c errtotal = errtotal - errbegin if (i_total-i_total/10*10.eq.1) then write(30,400) i_total, errtotal, dtmin, dtmax endif if (errtotal.lt.tol.or.itsub.eq.isubmax) then iconv = 1 else iconv = 0 endif Clp = rLp/(0.5*rho0*amach**2) Cdp = rDp/(0.5*rho0*amach**2) Cl = (rLp+rLf)/(0.5*rho0*amach**2) Cd = (rDp+rDf)/(0.5*rho0*amach**2) endif call MPI_BCAST(iconv,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(Clp,1,MPI_REAL,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(Cdp,1,MPI_REAL,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(Cl,1,MPI_REAL,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(Cd,1,MPI_REAL,0,MPI_COMM_WORLD,ierr) enddo do j = 1, nfy+1 do i = 1, nfx+1 do n = 1, 4 w_nm(n,i,j) = w_n(n,i,j) w_n(n,i,j) = w(n,i,j) enddo w_nm(5,i,j) = w_n(5,i,j) w_nm(6,i,j) = w_n(6,i,j) w_n(5,i,j) = w(5,i,j)*w(1,i,j) w_n(6,i,j) = w(6,i,j)*w(1,i,j) enddo enddo if (mype.eq.0) write(21,340) ptime, Clp, Cdp, Cl, Cd, itsub if (myPE.ne.0.and.mod(it-itbegin+1,intwrt).eq.0) then do j = 1, nfy+1 call MPI_SSEND(w_n(1,1,j),nfx+1,MPI_VEC,0,777, . MPI_COMM_WORLD,ierr) call MPI_SSEND(w_nm(1,1,j),nfx+1,MPI_VEC,0,888, . MPI_COMM_WORLD,ierr) call MPI_SSEND(vis(1,j),nfx+1,MPI_REAL,0,999, . MPI_COMM_WORLD,ierr) enddo endif if (myPE.eq.0.and.mod(it-itbegin+1,intwrt).eq.0) then do ii = 1, i_proc do jj = 1, j_proc i_cell = (nx-1)/i_proc j_cell = (ny-1)/j_proc i_remain = (nx-1)-(nx-1)/i_proc*i_proc j_remain = (ny-1)-(ny-1)/j_proc*j_proc if (ii.le.i_remain) i_cell = i_cell+1 if (jj.le.j_remain) j_cell = j_cell+1 i_end = ii*i_cell+1 j_end = jj*j_cell+1 if (ii.gt.i_remain) i_end = i_end+i_remain if (jj.gt.j_remain) j_end = j_end+j_remain i_begin = i_end-i_cell j_begin = j_end-j_cell if ((ii.eq.1).and.(jj.eq.1)) then do j = j_begin,j_end+1 do i = i_begin,i_end+1 do n = 1, 6 qq(n,i,j) = w_n(n,i,j) qm(n,i,j) = w_nm(n,i,j) enddo yyn(i,j) = vis(i,j) enddo enddo else do j = j_begin,j_end+1 call MPI_RECV(qq(1,i_begin,j),i_end-i_begin+2, . MPI_VEC,pe(ii,jj),777, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(qm(1,i_begin,j),i_end-i_begin+2, . MPI_VEC,pe(ii,jj),888, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(yyn(i_begin,j),i_end-i_begin+2, . MPI_REAL,pe(ii,jj),999, . MPI_COMM_WORLD,istatus,ierr) enddo endif enddo enddo call AFTER endif iconv = 0 enddo 400 format(1x,i6,3(1x,e11.5)) 500 format(1x) 340 format(1x,e12.6,4(' ',e12.6),' ',i5) call MPI_TYPE_FREE(MPI_VEC,ierr) call MPI_FINALIZE(ierr) end subroutine metric include 'mpins2d.inc' include 'mpif.h' call MPI_TYPE_VECTOR(nfy+1,1,imx,MPI_REAL,MPI_REAL_ROW,ierr) call MPI_TYPE_COMMIT(MPI_REAL_ROW,ierr) do j = 2, nfy jp = j + 1 do i = 2, nfx ip = i + 1 chx(i,j) = 0.5*( y(i,jp)+y(ip,jp)-y(i,j)-y(ip,j) ) chy(i,j) = 0.5*( x(i,jp)+x(ip,jp)-x(i,j)-x(ip,j) ) etx(i,j) = 0.5*( y(ip,j)+y(ip,jp)-y(i,j)-y(i,jp) ) ety(i,j) = 0.5*( x(ip,j)+x(ip,jp)-x(i,j)-x(i,jp) ) zac(i,j) = 1.0/ ( chx(i,j)*ety(i,j)-etx(i,j)*chy(i,j) ) chx(i,j) = zac(i,j)*chx(i,j) chy(i,j) = -zac(i,j)*chy(i,j) etx(i,j) = -zac(i,j)*etx(i,j) ety(i,j) = zac(i,j)*ety(i,j) enddo enddo if (jmPE.ge.0) then call MPI_SEND(chx(1,2),nfx+1,MPI_REAL,jmPE,32211, . MPI_COMM_WORLD,ierr) call MPI_SEND(chy(1,2),nfx+1,MPI_REAL,jmPE,32212, . MPI_COMM_WORLD,ierr) call MPI_SEND(etx(1,2),nfx+1,MPI_REAL,jmPE,32213, . MPI_COMM_WORLD,ierr) call MPI_SEND(ety(1,2),nfx+1,MPI_REAL,jmPE,32214, . MPI_COMM_WORLD,ierr) call MPI_SEND(zac(1,2),nfx+1,MPI_REAL,jmPE,32215, . MPI_COMM_WORLD,ierr) else do i = 1, nfx+1 chx(i,1) = chx(i,2) chy(i,1) = chy(i,2) etx(i,1) = etx(i,2) ety(i,1) = ety(i,2) zac(i,1) = zac(i,2) enddo endif if (jpPE.ge.0) then call MPI_RECV(chx(1,nfy+1),nfx+1,MPI_REAL,jpPE,32211, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(chy(1,nfy+1),nfx+1,MPI_REAL,jpPE,32212, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(etx(1,nfy+1),nfx+1,MPI_REAL,jpPE,32213, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(ety(1,nfy+1),nfx+1,MPI_REAL,jpPE,32214, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(zac(1,nfy+1),nfx+1,MPI_REAL,jpPE,32215, . MPI_COMM_WORLD,istatus,ierr) call MPI_SEND(chx(1,nfy),nfx+1,MPI_REAL,jpPE,32221, . MPI_COMM_WORLD,ierr) call MPI_SEND(chy(1,nfy),nfx+1,MPI_REAL,jpPE,32222, . MPI_COMM_WORLD,ierr) call MPI_SEND(etx(1,nfy),nfx+1,MPI_REAL,jpPE,32223, . MPI_COMM_WORLD,ierr) call MPI_SEND(ety(1,nfy),nfx+1,MPI_REAL,jpPE,32224, . MPI_COMM_WORLD,ierr) call MPI_SEND(zac(1,nfy),nfx+1,MPI_REAL,jpPE,32225, . MPI_COMM_WORLD,ierr) else do i = 1, nfx+1 chx(i,nfy+1) = chx(i,nfy) chy(i,nfy+1) = chy(i,nfy) etx(i,nfy+1) = etx(i,nfy) ety(i,nfy+1) = ety(i,nfy) zac(i,nfy+1) = zac(i,nfy) enddo endif if (jmPE.ge.0) then call MPI_RECV(chx(1,1),nfx+1,MPI_REAL,jmPE,32221, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(chy(1,1),nfx+1,MPI_REAL,jmPE,32222, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(etx(1,1),nfx+1,MPI_REAL,jmPE,32223, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(ety(1,1),nfx+1,MPI_REAL,jmPE,32224, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(zac(1,1),nfx+1,MPI_REAL,jmPE,32225, . MPI_COMM_WORLD,istatus,ierr) endif if (imPE.ge.0) then call MPI_SEND(chx(2,1),1,MPI_REAL_ROW,imPE,32311, . MPI_COMM_WORLD,ierr) call MPI_SEND(chy(2,1),1,MPI_REAL_ROW,imPE,32312, . MPI_COMM_WORLD,ierr) call MPI_SEND(etx(2,1),1,MPI_REAL_ROW,imPE,32313, . MPI_COMM_WORLD,ierr) call MPI_SEND(ety(2,1),1,MPI_REAL_ROW,imPE,32314, . MPI_COMM_WORLD,ierr) call MPI_SEND(zac(2,1),1,MPI_REAL_ROW,imPE,32315, . MPI_COMM_WORLD,ierr) else do j = 1, nfy+1 chx(1,j) = chx(2,j) chy(1,j) = chy(2,j) etx(1,j) = etx(2,j) ety(1,j) = ety(2,j) zac(1,j) = zac(2,j) enddo endif if (ipPE.ge.0) then call MPI_RECV(chx(nfx+1,1),1,MPI_REAL_ROW,ipPE,32311, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(chy(nfx+1,1),1,MPI_REAL_ROW,ipPE,32312, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(etx(nfx+1,1),1,MPI_REAL_ROW,ipPE,32313, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(ety(nfx+1,1),1,MPI_REAL_ROW,ipPE,32314, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(zac(nfx+1,1),1,MPI_REAL_ROW,ipPE,32315, . MPI_COMM_WORLD,istatus,ierr) call MPI_SEND(chx(nfx,1),1,MPI_REAL_ROW,ipPE,32321, . MPI_COMM_WORLD,ierr) call MPI_SEND(chy(nfx,1),1,MPI_REAL_ROW,ipPE,32322, . MPI_COMM_WORLD,ierr) call MPI_SEND(etx(nfx,1),1,MPI_REAL_ROW,ipPE,32323, . MPI_COMM_WORLD,ierr) call MPI_SEND(ety(nfx,1),1,MPI_REAL_ROW,ipPE,32324, . MPI_COMM_WORLD,ierr) call MPI_SEND(zac(nfx,1),1,MPI_REAL_ROW,ipPE,32325, . MPI_COMM_WORLD,ierr) else do j = 1, nfy+1 chx(nfx+1,j) = chx(nfx,j) chy(nfx+1,j) = chy(nfx,j) etx(nfx+1,j) = etx(nfx,j) ety(nfx+1,j) = ety(nfx,j) zac(nfx+1,j) = zac(nfx,j) enddo endif if (imPE.ge.0) then call MPI_RECV(chx(1,1),1,MPI_REAL_ROW,imPE,32321, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(chy(1,1),1,MPI_REAL_ROW,imPE,32322, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(etx(1,1),1,MPI_REAL_ROW,imPE,32323, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(ety(1,1),1,MPI_REAL_ROW,imPE,32324, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(zac(1,1),1,MPI_REAL_ROW,imPE,32325, . MPI_COMM_WORLD,istatus,ierr) endif call MPI_TYPE_FREE(MPI_REAL_ROW,ierr) return end subroutine metdum include 'mpins2d.inc' do j = 2, nfy jp = j + 1 do i = 1, nfx ip = i + 1 if ((ib(i,j).eq.0).and.(ib(ip,j).eq.1)) then zac(i,j) = zac(ip,j) else if ((ib(i,j).eq.1).and.(ib(ip,j).eq.0)) then zac(ip,j) = zac(i,j) endif xzac(i,j) = 0.5*( zac(i,j)+zac(ip,j) ) cx_x(i,j) = xzac(i,j)*(y(ip,jp)-y(ip,j) ) cy_x(i,j) = xzac(i,j)*(x(ip,j)-x(ip,jp) ) ex_x(i,j) = 0.5*( etx(i,j)+etx(ip,j) ) ey_x(i,j) = 0.5*( ety(i,j)+ety(ip,j) ) enddo enddo do j = 1, nfy jp = j + 1 do i = 2, nfx ip = i + 1 if ((ib(i,j).eq.0).and.(ib(i,jp).eq.1)) then zac(i,j) = zac(i,jp) else if ((ib(i,j).eq.1).and.(ib(i,jp).eq.0)) then zac(i,jp) = zac(i,j) endif yzac(i,j) = 0.5*( zac(i,j)+zac(i,jp) ) cx_y(i,j) = 0.5*( chx(i,j)+chx(i,jp) ) cy_y(i,j) = 0.5*( chy(i,j)+chy(i,jp) ) ex_y(i,j) = yzac(i,j)*(y(i,jp)-y(ip,jp) ) ey_y(i,j) = yzac(i,j)*(x(ip,jp)-x(i,jp) ) enddo enddo return end subroutine step include 'mpins2d.inc' dtmin = 100. dtmax = 0.0 do j = 1, nfy+1 do i = 1, nfx+1 uuu = (w(2,i,j)*chx(i,j)+w(3,i,j)*chy(i,j))/w(1,i,j) vvv = (w(2,i,j)*etx(i,j)+w(3,i,j)*ety(i,j))/w(1,i,j) cuu = c(i,j)*sqrt(chx(i,j)**2+chy(i,j)**2) cvv = c(i,j)*sqrt(etx(i,j)**2+ety(i,j)**2) chuv = abs(uuu)+cuu etuv = abs(vvv)+cvv dtc = 1./(chuv+etuv) t = c(i,j)**2 rmul = t**1.5*(1.+tcoe)/(t+tcoe) rmu = rmul/ra+vis(i,j) dtd = 1./(2.5*rmu*(chx(i,j)**2+chy(i,j)**2 . +etx(i,j)**2+ety(i,j)**2)) dt(i,j) = cfl*dtc*dtd/(dtc+dtd) dtmin = min(dtmin, dt(i,j)) dtmax = max(dtmax, dt(i,j)) enddo enddo return end subroutine physics include 'mpins2d.inc' common/flux/rh(2),uu(2),vv(2),pp(2),hh(2),cc(2),oo(2) dimension eflux(6,imx,jmx),fflux(6,imx,jmx) dimension ev(6,imx,jmx),fv(6,imx,jmx) do j = 2, nfy jp = j + 1 do i = 1, nfx_start - 1 ip = i + 1 if ((ib(i,j).eq.0).and.(ib(ip,j).eq.1)) then w(1,i,j) = w(1,ip,j) w(2,i,j) = -w(2,ip,j) w(3,i,j) = -w(3,ip,j) w(4,i,j) = w(4,ip,j) w(5,i,j) = -w(5,ip,j) p(i,j) = p(ip,j) vis(i,j) = -vis(ip,j) ff(i,j) = ff(ip,j) xx = 0.25*(x(ip,j)+x(ip,jp)+x(ip+1,j)+x(ip+1,jp)) yy = 0.25*(y(ip,j)+y(ip,jp)+y(ip+1,j)+y(ip+1,jp)) xc = 0.5*(x(ip,j)+x(ip,jp)) yc = 0.5*(y(ip,j)+y(ip,jp)) dyn = (xx-xc)**2+(yy-yc)**2 c ut = (w(2,ip,j)**2+w(3,ip,j)**2)/w(1,ip,j)**2 c owall = 2500*sqrt(ut/dyn) ccc = sqrt(gamma*p(ip,j)/w(1,ip,j)) t = ccc**2 rmul = t**1.5*(1.+tcoe)/(t+tcoe) owall = 60*rmul/(beta1*dyn*ra) w(6,i,j) = -w(6,ip,j)+2*owall else if ((ib(i,j).eq.1).and.(ib(ip,j).eq.0)) then w(1,ip,j) = w(1,i,j) w(2,ip,j) = -w(2,i,j) w(3,ip,j) = -w(3,i,j) w(4,ip,j) = w(4,i,j) w(5,ip,j) = -w(5,i,j) p(ip,j) = p(i,j) vis(ip,j) = -vis(i,j) ff(ip,j) = ff(i,j) xx = 0.25*(x(ip,j)+x(ip,jp)+x(i,j)+x(i,jp)) yy = 0.25*(y(ip,j)+y(ip,jp)+y(i,j)+y(i,jp)) xc = 0.5*(x(ip,j)+x(ip,jp)) yc = 0.5*(y(ip,j)+y(ip,jp)) dyn = (xx-xc)**2+(yy-yc)**2 c ut = (w(2,i,j)**2+w(3,i,j)**2)/w(1,i,j)**2 c owall = 2500*sqrt(ut/dyn) ccc = sqrt(gamma*p(i,j)/w(1,i,j)) t = ccc**2 rmul = t**1.5*(1.+tcoe)/(t+tcoe) owall = 60*rmul/(beta1*dyn*ra) w(6,ip,j) = -w(6,i,j)+2*owall endif rh(1) = w(1,i,j) uu(1) = w(2,i,j)/w(1,i,j) vv(1) = w(3,i,j)/w(1,i,j) pp(1) = p(i,j) cc(1) = w(5,i,j) oo(1) = w(6,i,j) rh(2) = w(1,ip,j) uu(2) = w(2,ip,j)/w(1,ip,j) vv(2) = w(3,ip,j)/w(1,ip,j) pp(2) = p(ip,j) cc(2) = w(5,ip,j) oo(2) = w(6,ip,j) hh(1) = gamma/gam1*pp(1)/rh(1) + 0.5*(uu(1)**2+vv(1)**2) hh(2) = gamma/gam1*pp(2)/rh(2) + 0.5*(uu(2)**2+vv(2)**2) cx = y(ip,jp) - y(ip,j) cy = x(ip,j) - x(ip,jp) call roeflux(eflux,cx,cy,i,j) enddo do i = nfx_end + 1, nfx ip = i + 1 if ((ib(i,j).eq.0).and.(ib(ip,j).eq.1)) then w(1,i,j) = w(1,ip,j) w(2,i,j) = -w(2,ip,j) w(3,i,j) = -w(3,ip,j) w(4,i,j) = w(4,ip,j) w(5,i,j) = -w(5,ip,j) p(i,j) = p(ip,j) vis(i,j) = -vis(ip,j) ff(i,j) = ff(ip,j) xx = 0.25*(x(ip,j)+x(ip,jp)+x(ip+1,j)+x(ip+1,jp)) yy = 0.25*(y(ip,j)+y(ip,jp)+y(ip+1,j)+y(ip+1,jp)) xc = 0.5*(x(ip,j)+x(ip,jp)) yc = 0.5*(y(ip,j)+y(ip,jp)) dyn = (xx-xc)**2+(yy-yc)**2 c ut = (w(2,ip,j)**2+w(3,ip,j)**2)/w(1,ip,j)**2 c owall = 2500*sqrt(ut/dyn) ccc = sqrt(gamma*p(ip,j)/w(1,ip,j)) t = ccc**2 rmul = t**1.5*(1.+tcoe)/(t+tcoe) owall = 60*rmul/(beta1*dyn*ra) w(6,i,j) = -w(6,ip,j)+2*owall else if ((ib(i,j).eq.1).and.(ib(ip,j).eq.0)) then w(1,ip,j) = w(1,i,j) w(2,ip,j) = -w(2,i,j) w(3,ip,j) = -w(3,i,j) w(4,ip,j) = w(4,i,j) w(5,ip,j) = -w(5,i,j) p(ip,j) = p(i,j) vis(ip,j) = -vis(i,j) ff(ip,j) = ff(i,j) xx = 0.25*(x(ip,j)+x(ip,jp)+x(i,j)+x(i,jp)) yy = 0.25*(y(ip,j)+y(ip,jp)+y(i,j)+y(i,jp)) xc = 0.5*(x(ip,j)+x(ip,jp)) yc = 0.5*(y(ip,j)+y(ip,jp)) dyn = (xx-xc)**2+(yy-yc)**2 c ut = (w(2,i,j)**2+w(3,i,j)**2)/w(1,i,j)**2 c owall = 2500*sqrt(ut/dyn) ccc = sqrt(gamma*p(i,j)/w(1,i,j)) t = ccc**2 rmul = t**1.5*(1.+tcoe)/(t+tcoe) owall = 60*rmul/(beta1*dyn*ra) w(6,ip,j) = -w(6,i,j)+2*owall endif rh(1) = w(1,i,j) uu(1) = w(2,i,j)/w(1,i,j) vv(1) = w(3,i,j)/w(1,i,j) pp(1) = p(i,j) cc(1) = w(5,i,j) oo(1) = w(6,i,j) rh(2) = w(1,ip,j) uu(2) = w(2,ip,j)/w(1,ip,j) vv(2) = w(3,ip,j)/w(1,ip,j) pp(2) = p(ip,j) cc(2) = w(5,ip,j) oo(2) = w(6,ip,j) hh(1) = gamma/gam1*pp(1)/rh(1) + 0.5*(uu(1)**2+vv(1)**2) hh(2) = gamma/gam1*pp(2)/rh(2) + 0.5*(uu(2)**2+vv(2)**2) cx = y(ip,jp) - y(ip,j) cy = x(ip,j) - x(ip,jp) call roeflux(eflux,cx,cy,i,j) enddo do i = nfx_start, nfx_end ip = i + 1 im = i - 1 if ((ib(i,j).eq.0).and.(ib(ip,j).eq.1)) then w(1,i,j) = w(1,ip,j) w(2,i,j) = -w(2,ip,j) w(3,i,j) = -w(3,ip,j) w(4,i,j) = w(4,ip,j) w(5,i,j) = -w(5,ip,j) p(i,j) = p(ip,j) vis(i,j) = -vis(ip,j) ff(i,j) = ff(ip,j) xx = 0.25*(x(ip,j)+x(ip,jp)+x(ip+1,j)+x(ip+1,jp)) yy = 0.25*(y(ip,j)+y(ip,jp)+y(ip+1,j)+y(ip+1,jp)) xc = 0.5*(x(ip,j)+x(ip,jp)) yc = 0.5*(y(ip,j)+y(ip,jp)) dyn = (xx-xc)**2+(yy-yc)**2 c ut = (w(2,ip,j)**2+w(3,ip,j)**2)/w(1,ip,j)**2 c owall = 2500*sqrt(ut/dyn) ccc = sqrt(gamma*p(ip,j)/w(1,ip,j)) t = ccc**2 rmul = t**1.5*(1.+tcoe)/(t+tcoe) owall = 60*rmul/(beta1*dyn*ra) w(6,i,j) = -w(6,ip,j)+2*owall else if ((ib(i,j).eq.1).and.(ib(ip,j).eq.0)) then w(1,ip,j) = w(1,i,j) w(2,ip,j) = -w(2,i,j) w(3,ip,j) = -w(3,i,j) w(4,ip,j) = w(4,i,j) w(5,ip,j) = -w(5,i,j) p(ip,j) = p(i,j) vis(ip,j) = -vis(i,j) ff(ip,j) = ff(i,j) xx = 0.25*(x(ip,j)+x(ip,jp)+x(i,j)+x(i,jp)) yy = 0.25*(y(ip,j)+y(ip,jp)+y(i,j)+y(i,jp)) xc = 0.5*(x(ip,j)+x(ip,jp)) yc = 0.5*(y(ip,j)+y(ip,jp)) dyn = (xx-xc)**2+(yy-yc)**2 c ut = (w(2,i,j)**2+w(3,i,j)**2)/w(1,i,j)**2 c owall = 2500*sqrt(ut/dyn) ccc = sqrt(gamma*p(i,j)/w(1,i,j)) t = ccc**2 rmul = t**1.5*(1.+tcoe)/(t+tcoe) owall = 60*rmul/(beta1*dyn*ra) w(6,ip,j) = -w(6,i,j)+2*owall endif rh(1) = w(1,i,j) uu(1) = w(2,i,j)/w(1,i,j) vv(1) = w(3,i,j)/w(1,i,j) pp(1) = p(i,j) cc(1) = w(5,i,j) oo(1) = w(6,i,j) rh(2) = w(1,ip,j) uu(2) = w(2,ip,j)/w(1,ip,j) vv(2) = w(3,ip,j)/w(1,ip,j) pp(2) = p(ip,j) cc(2) = w(5,ip,j) oo(2) = w(6,ip,j) call muscl(w(1,im,j),rh(1),rh(2),w(1,ip+1,j),ib(i,j),ib(ip,j)) call muscl(w(2,im,j)/w(1,im,j),uu(1),uu(2), . w(2,ip+1,j)/w(1,ip+1,j),ib(i,j),ib(ip,j)) call muscl(w(3,im,j)/w(1,im,j),vv(1),vv(2), . w(3,ip+1,j)/w(1,ip+1,j),ib(i,j),ib(ip,j)) call muscl(w(5,im,j),cc(1),cc(2),w(5,ip+1,j),ib(i,j),ib(ip,j)) call muscl(w(6,im,j),oo(1),oo(2),w(6,ip+1,j),ib(i,j),ib(ip,j)) call muscl(p(im,j),pp(1),pp(2),p(ip+1,j),ib(i,j),ib(ip,j)) hh(1) = gamma/gam1*pp(1)/rh(1) + 0.5*(uu(1)**2+vv(1)**2) hh(2) = gamma/gam1*pp(2)/rh(2) + 0.5*(uu(2)**2+vv(2)**2) cx = y(ip,jp) - y(ip,j) cy = x(ip,j) - x(ip,jp) call roeflux(eflux,cx,cy,i,j) enddo enddo do j = 1, nfy_start-1 jp = j + 1 do i = 2, nfx ip = i + 1 if ((ib(i,j).eq.0).and.(ib(i,jp).eq.1)) then w(1,i,j) = w(1,i,jp) w(2,i,j) = -w(2,i,jp) w(3,i,j) = -w(3,i,jp) w(4,i,j) = w(4,i,jp) w(5,i,j) = -w(5,i,jp) p(i,j) = p(i,jp) vis(i,j) = -vis(i,jp) ff(i,j) = ff(i,jp) xx = 0.25*(x(i,jp)+x(i+1,jp)+x(i,jp+1)+x(i+1,jp+1)) yy = 0.25*(y(i,jp)+y(i+1,jp)+y(i,jp+1)+y(i+1,jp+1)) xc = 0.5*(x(i,jp)+x(i+1,jp)) yc = 0.5*(y(i,jp)+y(i+1,jp)) dyn = (xx-xc)**2+(yy-yc)**2 c ut = (w(2,i,jp)**2+w(3,i,jp)**2)/w(1,i,jp)**2 c owall = 2500*sqrt(ut/dyn) ccc = sqrt(gamma*p(i,jp)/w(1,i,jp)) t = ccc**2 rmul = t**1.5*(1.+tcoe)/(t+tcoe) owall = 60*rmul/(beta1*dyn*ra) w(6,i,j) = -w(6,i,jp)+2*owall else if ((ib(i,j).eq.1).and.(ib(i,jp).eq.0)) then w(1,i,jp) = w(1,i,j) w(2,i,jp) = -w(2,i,j) w(3,i,jp) = -w(3,i,j) w(4,i,jp) = w(4,i,j) w(5,i,jp) = -w(5,i,j) p(i,jp) = p(i,j) vis(i,jp) = -vis(i,j) ff(i,jp) = ff(i,j) xx = 0.25*(x(i,jp)+x(i+1,jp)+x(i,j)+x(i+1,j)) yy = 0.25*(y(i,jp)+y(i+1,jp)+y(i,j)+y(i+1,j)) xc = 0.5*(x(i,jp)+x(i+1,jp)) yc = 0.5*(y(i,jp)+y(i+1,jp)) dyn = (xx-xc)**2+(yy-yc)**2 c ut = (w(2,i,j)**2+w(3,i,j)**2)/w(1,i,j)**2 c owall = 2500*sqrt(ut/dyn) ccc = sqrt(gamma*p(i,j)/w(1,i,j)) t = ccc**2 rmul = t**1.5*(1.+tcoe)/(t+tcoe) owall = 60*rmul/(beta1*dyn*ra) w(6,i,jp) = -w(6,i,j)+2*owall endif rh(1) = w(1,i,j) uu(1) = w(2,i,j)/w(1,i,j) vv(1) = w(3,i,j)/w(1,i,j) pp(1) = p(i,j) cc(1) = w(5,i,j) oo(1) = w(6,i,j) rh(2) = w(1,i,jp) uu(2) = w(2,i,jp)/w(1,i,jp) vv(2) = w(3,i,jp)/w(1,i,jp) pp(2) = p(i,jp) cc(2) = w(5,i,jp) oo(2) = w(6,i,jp) hh(1) = gamma/gam1*pp(1)/rh(1) + 0.5*(uu(1)**2+vv(1)**2) hh(2) = gamma/gam1*pp(2)/rh(2) + 0.5*(uu(2)**2+vv(2)**2) ex = y(i,jp) - y(ip,jp) ey = x(ip,jp) - x(i,jp) call roeflux(fflux,ex,ey,i,j) enddo enddo do j = nfy_end+1, nfy jp = j + 1 do i = 2, nfx ip = i + 1 if ((ib(i,j).eq.0).and.(ib(i,jp).eq.1)) then w(1,i,j) = w(1,i,jp) w(2,i,j) = -w(2,i,jp) w(3,i,j) = -w(3,i,jp) w(4,i,j) = w(4,i,jp) w(5,i,j) = -w(5,i,jp) p(i,j) = p(i,jp) vis(i,j) = -vis(i,jp) ff(i,j) = ff(i,jp) xx = 0.25*(x(i,jp)+x(i+1,jp)+x(i,jp+1)+x(i+1,jp+1)) yy = 0.25*(y(i,jp)+y(i+1,jp)+y(i,jp+1)+y(i+1,jp+1)) xc = 0.5*(x(i,jp)+x(i+1,jp)) yc = 0.5*(y(i,jp)+y(i+1,jp)) dyn = (xx-xc)**2+(yy-yc)**2 c ut = (w(2,i,jp)**2+w(3,i,jp)**2)/w(1,i,jp)**2 c owall = 2500*sqrt(ut/dyn) ccc = sqrt(gamma*p(i,jp)/w(1,i,jp)) t = ccc**2 rmul = t**1.5*(1.+tcoe)/(t+tcoe) owall = 60*rmul/(beta1*dyn*ra) w(6,i,j) = -w(6,i,jp)+2*owall else if ((ib(i,j).eq.1).and.(ib(i,jp).eq.0)) then w(1,i,jp) = w(1,i,j) w(2,i,jp) = -w(2,i,j) w(3,i,jp) = -w(3,i,j) w(4,i,jp) = w(4,i,j) w(5,i,jp) = -w(5,i,j) p(i,jp) = p(i,j) vis(i,jp) = -vis(i,j) ff(i,jp) = ff(i,j) xx = 0.25*(x(i,jp)+x(i+1,jp)+x(i,j)+x(i+1,j)) yy = 0.25*(y(i,jp)+y(i+1,jp)+y(i,j)+y(i+1,j)) xc = 0.5*(x(i,jp)+x(i+1,jp)) yc = 0.5*(y(i,jp)+y(i+1,jp)) dyn = (xx-xc)**2+(yy-yc)**2 c ut = (w(2,i,j)**2+w(3,i,j)**2)/w(1,i,j)**2 c owall = 2500*sqrt(ut/dyn) ccc = sqrt(gamma*p(i,j)/w(1,i,j)) t = ccc**2 rmul = t**1.5*(1.+tcoe)/(t+tcoe) owall = 60*rmul/(beta1*dyn*ra) w(6,i,jp) = -w(6,i,j)+2*owall endif rh(1) = w(1,i,j) uu(1) = w(2,i,j)/w(1,i,j) vv(1) = w(3,i,j)/w(1,i,j) pp(1) = p(i,j) cc(1) = w(5,i,j) oo(1) = w(6,i,j) rh(2) = w(1,i,jp) uu(2) = w(2,i,jp)/w(1,i,jp) vv(2) = w(3,i,jp)/w(1,i,jp) pp(2) = p(i,jp) cc(2) = w(5,i,jp) oo(2) = w(6,i,jp) hh(1) = gamma/gam1*pp(1)/rh(1) + 0.5*(uu(1)**2+vv(1)**2) hh(2) = gamma/gam1*pp(2)/rh(2) + 0.5*(uu(2)**2+vv(2)**2) ex = y(i,jp) - y(ip,jp) ey = x(ip,jp) - x(i,jp) call roeflux(fflux,ex,ey,i,j) enddo enddo do j = nfy_start, nfy_end jp = j + 1 jm = j - 1 do i = 2, nfx ip = i + 1 if ((ib(i,j).eq.0).and.(ib(i,jp).eq.1)) then w(1,i,j) = w(1,i,jp) w(2,i,j) = -w(2,i,jp) w(3,i,j) = -w(3,i,jp) w(4,i,j) = w(4,i,jp) w(5,i,j) = -w(5,i,jp) p(i,j) = p(i,jp) vis(i,j) = -vis(i,jp) ff(i,j) = ff(i,jp) xx = 0.25*(x(i,jp)+x(i+1,jp)+x(i,jp+1)+x(i+1,jp+1)) yy = 0.25*(y(i,jp)+y(i+1,jp)+y(i,jp+1)+y(i+1,jp+1)) xc = 0.5*(x(i,jp)+x(i+1,jp)) yc = 0.5*(y(i,jp)+y(i+1,jp)) dyn = (xx-xc)**2+(yy-yc)**2 c ut = (w(2,i,jp)**2+w(3,i,jp)**2)/w(1,i,jp)**2 c owall = 2500*sqrt(ut/dyn) ccc = sqrt(gamma*p(i,jp)/w(1,i,jp)) t = ccc**2 rmul = t**1.5*(1.+tcoe)/(t+tcoe) owall = 60*rmul/(beta1*dyn*ra) w(6,i,j) = -w(6,i,jp)+2*owall else if ((ib(i,j).eq.1).and.(ib(i,jp).eq.0)) then w(1,i,jp) = w(1,i,j) w(2,i,jp) = -w(2,i,j) w(3,i,jp) = -w(3,i,j) w(4,i,jp) = w(4,i,j) w(5,i,jp) = -w(5,i,j) p(i,jp) = p(i,j) vis(i,jp) = -vis(i,j) ff(i,jp) = ff(i,j) xx = 0.25*(x(i,jp)+x(i+1,jp)+x(i,j)+x(i+1,j)) yy = 0.25*(y(i,jp)+y(i+1,jp)+y(i,j)+y(i+1,j)) xc = 0.5*(x(i,jp)+x(i+1,jp)) yc = 0.5*(y(i,jp)+y(i+1,jp)) dyn = (xx-xc)**2+(yy-yc)**2 c ut = (w(2,i,j)**2+w(3,i,j)**2)/w(1,i,j)**2 c owall = 2500*sqrt(ut/dyn) ccc = sqrt(gamma*p(i,j)/w(1,i,j)) t = ccc**2 rmul = t**1.5*(1.+tcoe)/(t+tcoe) owall = 60*rmul/(beta1*dyn*ra) w(6,i,jp) = -w(6,i,j)+2*owall endif rh(1) = w(1,i,j) uu(1) = w(2,i,j)/w(1,i,j) vv(1) = w(3,i,j)/w(1,i,j) pp(1) = p(i,j) cc(1) = w(5,i,j) oo(1) = w(6,i,j) rh(2) = w(1,i,jp) uu(2) = w(2,i,jp)/w(1,i,jp) vv(2) = w(3,i,jp)/w(1,i,jp) pp(2) = p(i,jp) cc(2) = w(5,i,jp) oo(2) = w(6,i,jp) call muscl(w(1,i,jm),rh(1),rh(2),w(1,i,jp+1),ib(i,j),ib(i,jp)) call muscl(w(2,i,jm)/w(1,i,jm),uu(1),uu(2), . w(2,i,jp+1)/w(1,i,jp+1),ib(i,j),ib(i,jp)) call muscl(w(3,i,jm)/w(1,i,jm),vv(1),vv(2), . w(3,i,jp+1)/w(1,i,jp+1),ib(i,j),ib(i,jp)) call muscl(w(5,i,jm),cc(1),cc(2),w(5,i,jp+1),ib(i,j),ib(i,jp)) call muscl(w(6,i,jm),oo(1),oo(2),w(6,i,jp+1),ib(i,j),ib(i,jp)) call muscl(p(i,jm),pp(1),pp(2),p(i,jp+1),ib(i,j),ib(i,jp)) hh(1) = gamma/gam1*pp(1)/rh(1) + 0.5*(uu(1)**2+vv(1)**2) hh(2) = gamma/gam1*pp(2)/rh(2) + 0.5*(uu(2)**2+vv(2)**2) ex = y(i,jp) - y(ip,jp) ey = x(ip,jp) - x(i,jp) call roeflux(fflux,ex,ey,i,j) enddo enddo do j = 2, nfy jp = j + 1 do i = 1, nfx ip = i + 1 t = 0.5*(c(i,j)**2+c(ip,j)**2) rmul = t**1.5*(1.+tcoe)/(t+tcoe) rmut = 0.5*ra*(vis(i,j)+vis(ip,j)) rmu = rmul+rmut prcoef = rmul/prl + rmut/prt pr = rmu/prcoef f2 = 0.5*(ff(i,j)+ff(ip,j)) uip = w(2,ip,j)/w(1,ip,j) ui = w(2,i,j)/w(1,i,j) vip = w(3,ip,j)/w(1,ip,j) vi = w(3,i,j)/w(1,i,j) c2ip = c(ip,j)**2 c2i = c(i,j)**2 uface = 0.5*(uip+ui) vface = 0.5*(vip+vi) uup = unod(ip,jp) udown = unod(ip,j) vup = vnod(ip,jp) vdown = vnod(ip,j) c2up = cnod(ip,jp)**2 c2down = cnod(ip,j)**2 uch = uip-ui vch = vip-vi uet = uup-udown vet = vup-vdown c2ch = c2ip-c2i c2et = c2up-c2down rkup = rknod(ip,jp) rkdown = rknod(ip,j) rwup = rwnod(ip,jp) rwdown = rwnod(ip,j) rkch = w(5,ip,j) - w(5,i,j) rwch = w(6,ip,j) - w(6,i,j) rket = rkup - rkdown rwet = rwup - rwdown a1 = 4./3.*cx_x(i,j)**2+cy_x(i,j)**2 a2 = 4./3.*cx_x(i,j)*ex_x(i,j)+cy_x(i,j)*ey_x(i,j) a3 = 1./3.*cx_x(i,j)*cy_x(i,j) a4 = cy_x(i,j)*ex_x(i,j)-2./3.*cx_x(i,j)*ey_x(i,j) b1 = cx_x(i,j)**2+4./3.*cy_x(i,j)**2 b2 = cx_x(i,j)*ex_x(i,j)+4./3.*cy_x(i,j)*ey_x(i,j) b3 = 1./3.*cx_x(i,j)*cy_x(i,j) b4 = cx_x(i,j)*ey_x(i,j)-2./3.*cy_x(i,j)*ex_x(i,j) c1 = cx_x(i,j)**2+cy_x(i,j)**2 c2 = cx_x(i,j)*ex_x(i,j)+cy_x(i,j)*ey_x(i,j) omet = cx_x(i,j)**2+cy_x(i,j)**2 pmet = cx_x(i,j)*ex_x(i,j)+cy_x(i,j)*ey_x(i,j) sigmak = f2*sigmak1+(1.-f2)*sigmak2 sigmao = f2*sigmao1+(1.-f2)*sigmao2 aa = rmul+sigmak*rmut bb = rmul+sigmao*rmut ev(1,i,j) = 0.0 ev(2,i,j) = a1*uch+a2*uet+a3*vch+a4*vet ev(3,i,j) = b1*vch+b2*vet+b3*uch+b4*uet ev(4,i,j) = uface*ev(2,i,j)+vface*ev(3,i,j) . +(c1*c2ch+c2*c2et)/(gam1*pr) ev(5,i,j) = aa*(omet*rkch+pmet*rket) ev(6,i,j) = bb*(omet*rwch+pmet*rwet) do n = 2, 4 ev(n,i,j) = rmu*ev(n,i,j)/xzac(i,j) enddo do n = 5, 6 ev(n,i,j) = ev(n,i,j)/xzac(i,j) enddo enddo enddo do j = 1, nfy jp = j + 1 do i = 2, nfx ip = i + 1 t = 0.5*(c(i,j)**2+c(i,jp)**2) rmul = t**1.5*(1.+tcoe)/(t+tcoe) rmut = 0.5*ra*(vis(i,j)+vis(i,jp)) rmu = rmul + rmut prcoef = rmul/prl + rmut/prt pr = rmu/prcoef f2 = 0.5*(ff(i,j)+ff(i,jp)) ujp = w(2,i,jp)/w(1,i,jp) uj = w(2,i,j)/w(1,i,j) vjp = w(3,i,jp)/w(1,i,jp) vj = w(3,i,j)/w(1,i,j) c2jp = c(i,jp)**2 c2j = c(i,j)**2 uface = 0.5*(ujp+uj) vface = 0.5*(vjp+vj) uup = unod(ip,jp) udown = unod(i,jp) vup = vnod(ip,jp) vdown = vnod(i,jp) c2up = cnod(ip,jp)**2 c2down = cnod(i,jp)**2 uch = uup-udown vch = vup-vdown uet = ujp-uj vet = vjp-vj c2ch = c2up-c2down c2et = c2jp-c2j rkup = rknod(ip,jp) rkdown = rknod(i,jp) rwup = rwnod(ip,jp) rwdown = rwnod(i,jp) rkch = rkup - rkdown rwch = rwup - rwdown rket = w(5,i,jp) - w(5,i,j) rwet = w(6,i,jp) - w(6,i,j) a1 = 4./3.*ex_y(i,j)**2+ey_y(i,j)**2 a2 = 4./3.*cx_y(i,j)*ex_y(i,j)+cy_y(i,j)*ey_y(i,j) a3 = 1./3.*ex_y(i,j)*ey_y(i,j) a4 = cx_y(i,j)*ey_y(i,j)-2./3.*cy_y(i,j)*ex_y(i,j) b1 = ex_y(i,j)**2+4./3.*ey_y(i,j)**2 b2 = cx_y(i,j)*ex_y(i,j)+4./3.*cy_y(i,j)*ey_y(i,j) b3 = 1./3.*ex_y(i,j)*ey_y(i,j) b4 = cy_y(i,j)*ex_y(i,j)-2./3.*cx_y(i,j)*ey_y(i,j) c1 = ex_y(i,j)**2+ey_y(i,j)**2 c2 = cx_y(i,j)*ex_y(i,j)+cy_y(i,j)*ey_y(i,j) omet = cx_y(i,j)*ex_y(i,j)+cy_y(i,j)*ey_y(i,j) pmet = ex_y(i,j)**2+ey_y(i,j)**2 sigmak = f2*sigmak1+(1.-f2)*sigmak2 sigmao = f2*sigmao1+(1.-f2)*sigmao2 aa = rmul+sigmak*rmut bb = rmul+sigmao*rmut fv(1,i,j) = 0.0 fv(2,i,j) = a1*uet+a2*uch+a3*vet+a4*vch fv(3,i,j) = b1*vet+b2*vch+b3*uet+b4*uch fv(4,i,j) = uface*fv(2,i,j)+vface*fv(3,i,j) . +(c1*c2et+c2*c2ch)/(gam1*pr) fv(5,i,j) = aa*(omet*rkch+pmet*rket) fv(6,i,j) = bb*(omet*rwch+pmet*rwet) do n = 2, 4 fv(n,i,j) = rmu*fv(n,i,j)/yzac(i,j) enddo do n = 5, 6 fv(n,i,j) = fv(n,i,j)/yzac(i,j) enddo enddo enddo do j = 2, nfy jm = j - 1 do i = 2, nfx im = i - 1 do n = 1, 6 rhs(n,i,j) = eflux(n,i,j)-eflux(n,im,j) . +fflux(n,i,j)-fflux(n,i,jm) . -(ev(n,i,j)-ev(n,im,j) . + fv(n,i,j)-fv(n,i,jm))/ra enddo enddo enddo return end subroutine roeflux(eflux,cex,cey,i,j) include 'mpins2d.inc' common/flux/rh(2),uu(2),vv(2),pp(2),hh(2),cc(2),oo(2) dimension eflux(6,imx,jmx),fl(4),fr(4), . df1(4),df2(4),df3(4),df4(4) dimension ram(4) dss = sqrt(cey**2+cex**2) chxn = cex/dss chyn = cey/dss rrr = sqrt(rh(2)/rh(1)) rrrp1 = rrr + 1.0 rhoav = rrr*rh(1) uav = ( rrr*uu(2)+uu(1) )/rrrp1 vav = ( rrr*vv(2)+vv(1) )/rrrp1 hav = ( rrr*hh(2)+hh(1) )/rrrp1 cavsq = gam1*( hav-0.5*(uav**2+vav**2) ) cav = sqrt(cavsq) qavsq = uav**2+vav**2 dqsq = uu(2)**2+vv(2)**2-uu(1)**2-vv(1)**2 contra= chxn*uav+chyn*vav dp = pp(2)-pp(1) drho = rh(2)-rh(1) du = uu(2)-uu(1) dv = vv(2)-vv(1) ram(1) = abs( contra ) ram(2) = ram(1) ram(3) = abs( contra+cav ) ram(4) = abs( contra-cav ) do n= 1, 4 if(ram(n).lt.eee)then ram(n) = 0.5*(ram(n)**2/eee + eee) endif enddo ram(1) = ram(1)*dss ram(2) = ram(2)*dss ram(3) = ram(3)*dss ram(4) = ram(4)*dss chduv = chxn*du +chyn*dv coeff1 = ram(1)*(drho-dp/cavsq) coeff2 = ram(2)*rhoav coeff3 = ram(3)*0.5*(dp+rhoav*cav*chduv)/cavsq coeff4 = ram(4)*0.5*(dp-rhoav*cav*chduv)/cavsq df1(1) = coeff1 df1(2) = coeff1*uav df1(3) = coeff1*vav df1(4) = 0.5*coeff1*qavsq df2(1) = 0.0 df2(2) = coeff2*(du - chxn*chduv) df2(3) = coeff2*(dv - chyn*chduv) df2(4) = coeff2*(0.5*dqsq - contra*chduv) df3(1) = coeff3 df3(2) = coeff3*(uav + chxn*cav) df3(3) = coeff3*(vav + chyn*cav) df3(4) = coeff3*(hav + cav*contra) df4(1) = coeff4 df4(2) = coeff4*(uav - chxn*cav) df4(3) = coeff4*(vav - chyn*cav) df4(4) = coeff4*(hav - cav*contra) commo = rh(1)*( cex*uu(1) + cey*vv(1) ) fl(1) = commo fl(2) = uu(1)*commo+cex*pp(1) fl(3) = vv(1)*commo+cey*pp(1) fl(4) = hh(1)*commo commo = rh(2)*( cex*uu(2) + cey*vv(2) ) fr(1) = commo fr(2) = uu(2)*commo+cex*pp(2) fr(3) = vv(2)*commo+cey*pp(2) fr(4) = hh(2)*commo do n= 1, 4 eflux(n,i,j)=0.5*(fl(n)+fr(n)) . +0.5*(-df1(n)-df2(n)-df3(n)-df4(n)) enddo eflux(5,i,j)=0.5*(cc(1)*fl(1)+cc(2)*fr(1)) eflux(6,i,j)=0.5*(oo(1)*fl(1)+oo(2)*fr(1)) delk = (cc(2) - cc(1))*rhoav delo = (oo(2) - oo(1))*rhoav eflux(5,i,j) = eflux(5,i,j)-0.5*ram(1)*delk eflux(6,i,j) = eflux(6,i,j)-0.5*ram(1)*delo return end subroutine node include 'mpins2d.inc' do i = 2, nfx+1 im = i - 1 do j = 2, nfy+1 jm = j - 1 uij = w(2,i,j)/w(1,i,j) uimj = w(2,im,j)/w(1,im,j) uimjm = w(2,im,jm)/w(1,im,jm) uijm = w(2,i,jm)/w(1,i,jm) unod(i,j) = 0.25*(uij+uimj+uimjm+uijm) vij = w(3,i,j)/w(1,i,j) vimj = w(3,im,j)/w(1,im,j) vimjm = w(3,im,jm)/w(1,im,jm) vijm = w(3,i,jm)/w(1,i,jm) vnod(i,j) = 0.25*(vij+vimj+vimjm+vijm) cnod(i,j) = 0.25*(c(i,j)+c(im,j)+c(im,jm)+c(i,jm)) rknod(i,j) = 0.25*(w(5,i,j)+w(5,im,j)+w(5,im,jm)+w(5,i,jm)) rwnod(i,j) = 0.25*(w(6,i,j)+w(6,im,j)+w(6,im,jm)+w(6,i,jm)) enddo enddo return end subroutine upgrade include 'mpins2d.inc' do i = 2, nfx do j = 2, nfy w(1,i,j) = w(1,i,j) + dw(1,i,j) w(2,i,j) = w(2,i,j) + dw(2,i,j) w(3,i,j) = w(3,i,j) + dw(3,i,j) w(4,i,j) = w(4,i,j) + dw(4,i,j) w(5,i,j) = w(5,i,j) + dw(5,i,j)/w(1,i,j) w(6,i,j) = w(6,i,j) + dw(6,i,j)/w(1,i,j) if (w(5,i,j).lt.0) w(5,i,j) = . 0.25*(rknod(i,j)+rknod(i+1,j)+rknod(i,j+1)+rknod(i+1,j+1)) if (w(6,i,j).le.0) w(6,i,j) = . 0.25*(rwnod(i,j)+rwnod(i+1,j)+rwnod(i,j+1)+rwnod(i+1,j+1)) enddo enddo return end subroutine renew include 'mpins2d.inc' do j = nfy_start-1, 0 do i = 2,nfx p(i,j) = gam1*(w(4,i,j) . -0.5*(w(2,i,j)**2+w(3,i,j)**2)/w(1,i,j)) enddo enddo do j = nfy+2, nfy_end+2 do i = 2,nfx p(i,j) = gam1*(w(4,i,j) . -0.5*(w(2,i,j)**2+w(3,i,j)**2)/w(1,i,j)) enddo enddo do j = 2, nfy do i = nfx_start-1, 0 p(i,j) = gam1*(w(4,i,j) . -0.5*(w(2,i,j)**2+w(3,i,j)**2)/w(1,i,j)) enddo do i = nfx+2,nfx_end+2 p(i,j) = gam1*(w(4,i,j) . -0.5*(w(2,i,j)**2+w(3,i,j)**2)/w(1,i,j)) enddo enddo do j = 1, nfy+1 do i = 1, nfx+1 p(i,j) = gam1*(w(4,i,j) . -0.5*(w(2,i,j)**2+w(3,i,j)**2)/w(1,i,j)) c(i,j) = sqrt(gamma*p(i,j)/w(1,i,j)) enddo enddo return end subroutine bc include 'mpins2d.inc' include 'mpif.h' call MPI_TYPE_VECTOR(nfy+1,1,imx+1,MPI_VEC,MPI_VEC_ROW,ierr) call MPI_TYPE_COMMIT(MPI_VEC_ROW,ierr) if (jmPE.ge.0) then call MPI_SEND(w(1,1,2),nfx+1,MPI_VEC,jmPE,211, . MPI_COMM_WORLD,ierr) call MPI_SEND(w(1,1,3),nfx+1,MPI_VEC,jmPE,212, . MPI_COMM_WORLD,ierr) else call BCJmin endif if (jpPE.ge.0) then call MPI_RECV(w(1,1,nfy+1),nfx+1,MPI_VEC,jpPE,211, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(w(1,1,nfy+2),nfx+1,MPI_VEC,jpPE,212, . MPI_COMM_WORLD,istatus,ierr) call MPI_SEND(w(1,1,nfy-1),nfx+1,MPI_VEC,jpPE,221, . MPI_COMM_WORLD,ierr) call MPI_SEND(w(1,1,nfy),nfx+1,MPI_VEC,jpPE,222, . MPI_COMM_WORLD,ierr) else call BCJmax endif if (jmPE.ge.0) then call MPI_RECV(w(1,1,0),nfx+1,MPI_VEC,jmPE,221, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(w(1,1,1),nfx+1,MPI_VEC,jmPE,222, . MPI_COMM_WORLD,istatus,ierr) endif if (imPE.ge.0) then call MPI_SEND(w(1,2,1),1,MPI_VEC_ROW,imPE,311, . MPI_COMM_WORLD,ierr) call MPI_SEND(w(1,3,1),1,MPI_VEC_ROW,imPE,312, . MPI_COMM_WORLD,ierr) else call BCImin endif if (ipPE.ge.0) then call MPI_RECV(w(1,nfx+1,1),1,MPI_VEC_ROW,ipPE,311, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(w(1,nfx+2,1),1,MPI_VEC_ROW,ipPE,312, . MPI_COMM_WORLD,istatus,ierr) call MPI_SEND(w(1,nfx-1,1),1,MPI_VEC_ROW,ipPE,321, . MPI_COMM_WORLD,ierr) call MPI_SEND(w(1,nfx,1),1,MPI_VEC_ROW,ipPE,322, . MPI_COMM_WORLD,ierr) else call BCImax endif if (imPE.ge.0) then call MPI_RECV(w(1,0,1),1,MPI_VEC_ROW,imPE,321, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(w(1,1,1),1,MPI_VEC_ROW,imPE,322, . MPI_COMM_WORLD,istatus,ierr) endif call MPI_TYPE_FREE(MPI_VEC_ROW,ierr) return end subroutine BCImin include 'mpins2d.inc' if (imPE.eq.RIEMANN) then call IminRiemann else if (imPE.eq.SYMMETRIC) then call IminSymmetric else if (imPE.eq.FAR) then call IminFar endif return end subroutine BCImax include 'mpins2d.inc' if (ipPE.eq.RIEMANN) then call ImaxRiemann else if (ipPE.eq.OUTFLOW) then call ImaxOutflow else if (ipPE.eq.FAR) then call ImaxFar endif return end subroutine IminSymmetric include 'mpins2d.inc' do j = 2, nfy w(1,1,j) = w(1,2,j) w(2,1,j) = w(2,2,j) w(3,1,j) = -w(3,2,j) w(4,1,j) = w(4,2,j) w(5,1,j) = w(5,2,j) w(6,1,j) = w(6,2,j) vis(1,j) = vis(2,j) ff(1,j) = ff(2,j) enddo return end subroutine IminFar include 'mpins2d.inc' do j = 2, nfy w(1,1,j) = rho0 w(2,1,j) = rho0*u0 w(3,1,j) = rho0*v0 w(4,1,j) = rho0*e0 w(5,1,j) = rk0 w(6,1,j) = rw0 vis(1,j) = vis(2,j) ff(1,j) = ff(2,j) enddo return end subroutine ImaxOutflow include 'mpins2d.inc' do j = 2, nfy w(1,nfx+1,j) = w(1,nfx,j) w(2,nfx+1,j) = w(2,nfx,j) w(3,nfx+1,j) = w(3,nfx,j) w(4,nfx+1,j) = w(4,nfx,j) w(5,nfx+1,j) = w(5,nfx,j) w(6,nfx+1,j) = w(6,nfx,j) vis(nfx+1,j) = vis(nfx,j) ff(nfx+1,j) = ff(nfx,j) enddo return end subroutine ImaxFar include 'mpins2d.inc' do j = 2, nfy w(1,nfx+1,j) = rho0 w(2,nfx+1,j) = rho0*u0 w(3,nfx+1,j) = rho0*v0 w(4,nfx+1,j) = rho0*e0 w(5,nfx+1,j) = rk0 w(6,nfx+1,j) = rw0 vis(nfx+1,j) = vis(nfx,j) ff(nfx+1,j) = ff(nfx,j) enddo return end subroutine IminRiemann include 'mpins2d.inc' do j = 2, nfy jp = j + 1 cx = y(2,jp)-y(2,j) cy = x(2,j)-x(2,jp) dd = sqrt(cx**2+cy**2) cx = cx/dd cy = cy/dd qne = (w(2,2,j)*cx+w(3,2,j)*cy)/w(1,2,j) qte = (-w(2,2,j)*cy+w(3,2,j)*cx)/w(1,2,j) qnf = u0*cx+v0*cy qtf = -u0*cy+v0*cx R2 = qnf+2*c0/gam1 R1 = qne-2*c(2,j)/gam1 qn = 0.5*(R1+R2) c(1,j) = 0.25*gam1*(R2-R1) if (qn.ge.0) then qt = qtf w(1,1,j) = . (rho0**gamma*c(1,j)**2/p0/gamma)**(1./gam1) w(5,1,j) = rk0 w(6,1,j) = rw0 else qt = qte w(1,1,j) = . (w(1,2,j)**gamma*c(1,j)**2/p(2,j)/gamma)**(1./gam1) w(5,1,j) = w(5,2,j) w(6,1,j) = w(6,2,j) endif p(1,j) = w(1,1,j)*c(1,j)**2/gamma w(2,1,j) = w(1,1,j)*(qn*cx - qt*cy) w(3,1,j) = w(1,1,j)*(qn*cy + qt*cx) w(4,1,j) = p(1,j)/gam1 . +0.5*(w(2,1,j)**2+w(3,1,j)**2)/w(1,1,j) vis(1,j) = vis(2,j) ff(1,j) = ff(2,j) enddo return end subroutine ImaxRiemann include 'mpins2d.inc' do j = 2, nfy jp = j + 1 cx = y(nfx+1,jp)-y(nfx+1,j) cy = x(nfx+1,j)-x(nfx+1,jp) dd = sqrt(cx**2+cy**2) cx = cx/dd cy = cy/dd qne = (w(2,nfx,j)*cx+w(3,nfx,j)*cy)/w(1,nfx,j) qte = (-w(2,nfx,j)*cy+w(3,nfx,j)*cx)/w(1,nfx,j) qnf = u0*cx+v0*cy qtf = -u0*cy+v0*cx R2 = qne+2*c(nfx,j)/gam1 R1 = qnf-2*c0/gam1 qn = 0.5*(R1+R2) c(nfx+1,j) = 0.25*gam1*(R2-R1) if (qn.le.0) then qt = qtf w(1,nfx+1,j) = . (rho0**gamma*c(nfx+1,j)**2/p0/gamma)**(1./gam1) w(5,nfx+1,j) = rk0 w(6,nfx+1,j) = rw0 else qt = qte w(1,nfx+1,j) = . (w(1,nfx,j)**gamma*c(nfx+1,j)**2/p(nfx,j)/gamma)**(1./gam1) w(5,nfx+1,j) = w(5,nfx,j) w(6,nfx+1,j) = w(6,nfx,j) endif p(nfx+1,j) = w(1,nfx+1,j)*c(nfx+1,j)**2/gamma w(2,nfx+1,j) = w(1,nfx+1,j)*(qn*cx - qt*cy) w(3,nfx+1,j) = w(1,nfx+1,j)*(qn*cy + qt*cx) w(4,nfx+1,j) = p(nfx+1,j)/gam1 . +0.5*(w(2,nfx+1,j)**2+w(3,nfx+1,j)**2)/w(1,nfx+1,j) vis(nfx+1,j) = vis(nfx,j) ff(nfx+1,j) = ff(nfx,j) enddo return end subroutine BCJmin include 'mpins2d.inc' if (jmPE.eq.RIEMANN) then call JminRiemann else if (jmPE.eq.WALL) then call JminWall else if (jmPE.eq.FAR) then call JminFar else if (jmPE.eq.OUTFLOW) then call JminOutflow else if (jmPE.eq.SYMMETRIC) then call JminSymmetric endif return end subroutine BCJmax include 'mpins2d.inc' if (jpPE.eq.RIEMANN) then call JmaxRiemann else if (jpPE.eq.WALL) then call JmaxWall else if (jpPE.eq.FAR) then call JmaxFar else if (jpPE.eq.OUTFLOW) then call JmaxOutflow else if (jpPE.eq.SYMMETRIC) then call JmaxSymmetric else if (jpPE.eq.VORTEXCOR) then call VortexCorrection endif return end subroutine JminSymmetric include 'mpins2d.inc' do i = 2, nfx w(1,i,1) = w(1,i,2) w(2,i,1) = w(2,i,2) w(3,i,1) = -w(3,i,2) w(4,i,1) = w(4,i,2) w(5,i,1) = w(5,i,2) w(6,i,1) = w(6,i,2) vis(i,1) = vis(i,2) ff(i,1) = ff(i,2) enddo return end subroutine JminWall include 'mpins2d.inc' do i = 2, nfx w(1,i,1) = w(1,i,2) w(2,i,1) = -w(2,i,2) w(3,i,1) = -w(3,i,2) w(4,i,1) = w(4,i,2) xx = 0.25*(x(i,2)+x(i+1,2)+x(i,3)+x(i+1,3)) yy = 0.25*(y(i,2)+y(i+1,2)+y(i,3)+y(i+1,3)) xc = 0.5*(x(i,2)+x(i+1,2)) yc = 0.5*(y(i,2)+y(i+1,2)) dyn = (xx-xc)**2+(yy-yc)**2 c ut = (w(2,i,2)**2+w(3,i,2)**2)/w(1,i,2)**2 c owall = 2500*sqrt(ut/dyn) ccc = sqrt(gamma*p(i,2)/w(1,i,2)) t = ccc**2 rmul = t**1.5*(1.+tcoe)/(t+tcoe) owall = 60*rmul/(beta1*dyn*ra) w(5,i,1) = -w(5,i,2) w(6,i,1) = -w(6,i,2)+2*owall vis(i,1) = -vis(i,2) ff(i,1) = ff(i,2) enddo return end subroutine JminFar include 'mpins2d.inc' do i = 2, nfx w(1,i,1) = rho0 w(2,i,1) = rho0*u0 w(3,i,1) = rho0*v0 w(4,i,1) = rho0*e0 w(5,i,1) = rk0 w(6,i,1) = rw0 vis(i,1) = vis(i,2) ff(i,1) = ff(i,2) enddo return end subroutine JminOutflow include 'mpins2d.inc' do i = 2, nfx w(1,i,1) = w(1,i,2) w(2,i,1) = w(2,i,2) w(3,i,1) = w(3,i,2) w(4,i,1) = w(4,i,2) w(5,i,1) = w(5,i,2) w(6,i,1) = w(6,i,2) vis(i,1) = vis(i,2) ff(i,1) = ff(i,2) enddo return end subroutine JmaxSymmetric include 'mpins2d.inc' do i = 2, nfx w(1,i,nfy+1) = w(1,i,nfy) w(2,i,nfy+1) = w(2,i,nfy) w(3,i,nfy+1) = -w(3,i,nfy) w(4,i,nfy+1) = w(4,i,nfy) w(5,i,nfy+1) = w(5,i,nfy) w(6,i,nfy+1) = w(6,i,nfy) vis(i,nfy+1) = vis(i,nfy) ff(i,nfy+1) = ff(i,nfy) enddo return end subroutine VortexCorrection include 'mpins2d.inc' x0 = 0.5 ei0 = p0/(gam1*rho0) h0 = gamma*ei0 + 0.5*(u0**2+v0**2) s0 = rho0**gamma/p0 gm = 1.0/gam1 gmg = gam1/gamma beta = sqrt(1.0 - amach**2) circ = 0.25*clp*amach*beta/pi do i=2,nfx xa = 0.5*( x(i+1,nfy+1) + x(i,nfy+1)) - x0 ya = 0.5*( y(i+1,nfy+1) + y(i,nfy+1)) r = sqrt( xa**2 + ya**2 ) angl = atan2(ya,xa) c1 = cos(angl) s = sin(angl) qv = circ/(r*(1.0 - (amach*(s*cos(alpa) . - c1*sin(alpa)))**2)) uf = u0 + qv*s vf = v0 - qv*c1 cf = sqrt(gam1*(h0 - 0.5*(uf**2 + vf**2))) xx = x(i+1,nfy+1) - x(i,nfy+1) yx = y(i+1,nfy+1) - y(i,nfy+1) dd = 1.0/sqrt(xx**2 + yx**2) xx = xx*dd yx = yx*dd u = w(2,i,nfy)/w(1,i,nfy) v = w(3,i,nfy)/w(1,i,nfy) qq1 = 0.5*(w(2,i,nfy)**2. + w(3,i,nfy)**2.)/w(1,i,nfy) p1 = gam1*(w(4,i,nfy) - qq1) c2 = sqrt(gamma*p1/w(1,i,nfy)) er = xx*v - yx*u + 2.0*gm*c2 fr = xx*vf - yx*uf - 2.0*gm*cf qn = 0.5*(er + fr) qt = xx*uf + yx*vf s = s0 if(qn.gt.0.0) then qt = xx*u + yx*v s = w(1,i,nfy)**gamma/p1 endif u = xx*qt - yx*qn v = xx*qn + yx*qt qq1 = u**2 + v**2 cc = gmg*(h0 - 0.5*qq1) w(1,i,nfy+1)= (s*cc)**gm w(2,i,nfy+1)= w(1,i,nfy+1)*u w(3,i,nfy+1)= w(1,i,nfy+1)*v p1 = w(1,i,nfy+1)*cc w(4,i,nfy+1)= gm*p1 + 0.5*w(1,i,nfy+1)*qq1 w(5,i,nfy+1) = w(5,i,nfy) w(6,i,nfy+1) = w(6,i,nfy) vis(i,nfy+1) = vis(i,nfy) ff(i,nfy+1) = ff(i,nfy) enddo return end subroutine JmaxWall include 'mpins2d.inc' do i = 2, nfx w(1,i,nfy+1) = w(1,i,nfy) w(2,i,nfy+1) = -w(2,i,nfy) w(3,i,nfy+1) = -w(3,i,nfy) w(4,i,nfy+1) = w(4,i,nfy) xx = 0.25*(x(i,nfy+1)+x(i+1,nfy+1)+x(i,nfy)+x(i+1,nfy)) yy = 0.25*(y(i,nfy+1)+y(i+1,nfy+1)+y(i,nfy)+y(i+1,nfy)) xc = 0.5*(x(i,nfy+1)+x(i+1,nfy+1)) yc = 0.5*(y(i,nfy+1)+y(i+1,nfy+1)) dyn = (xx-xc)**2+(yy-yc)**2 c ut = (w(2,i,nfy)**2+w(3,i,nfy)**2)/w(1,i,nfy)**2 c owall = 2500*sqrt(ut/dyn) ccc = sqrt(gamma*p(i,nfy)/w(1,i,nfy)) t = ccc**2 rmul = t**1.5*(1.+tcoe)/(t+tcoe) owall = 60*rmul/(beta1*dyn*ra) w(5,i,nfy+1) = -w(5,i,nfy) w(6,i,nfy+1) = -w(6,i,nfy)+2*owall vis(i,nfy+1) = -vis(i,nfy) ff(i,nfy+1) = ff(i,nfy) enddo return end subroutine JmaxFar include 'mpins2d.inc' do i = 2, nfx w(1,i,nfy+1) = rho0 w(2,i,nfy+1) = rho0*u0 w(3,i,nfy+1) = rho0*v0 w(4,i,nfy+1) = rho0*e0 w(5,i,nfy+1) = rk0 w(6,i,nfy+1) = rw0 vis(i,nfy+1) = vis(i,nfy) ff(i,nfy+1) = ff(i,nfy) enddo return end subroutine JmaxOutflow include 'mpins2d.inc' do i = 2, nfx w(1,i,nfy+1) = w(1,i,nfy) w(2,i,nfy+1) = w(2,i,nfy) w(3,i,nfy+1) = w(3,i,nfy) w(4,i,nfy+1) = w(4,i,nfy) w(5,i,nfy+1) = w(5,i,nfy) w(6,i,nfy+1) = w(6,i,nfy) vis(i,nfy+1) = vis(i,nfy) ff(i,nfy+1) = ff(i,nfy) enddo return end subroutine JminRiemann include 'mpins2d.inc' do i = 2, nfx ip = i + 1 ex = y(i,2)-y(ip,2) ey = x(ip,2)-x(i,2) dd = sqrt(ex**2+ey**2) ex = ex/dd ey = ey/dd qne = (w(2,i,2)*ex+w(3,i,2)*ey)/w(1,i,2) qte = (-w(2,i,2)*ey+w(3,i,2)*ex)/w(1,i,2) qnf = u0*ex+v0*ey qtf = -u0*ey+v0*ex R2 = qnf+2*c0/gam1 R1 = qne-2*c(i,2)/gam1 qn = 0.5*(R1+R2) c(i,1) = 0.25*gam1*(R2-R1) if (qn.ge.0) then qt = qtf w(1,i,1) = . (rho0**gamma*c(i,1)**2/p0/gamma)**(1./gam1) w(5,i,1) = rk0 w(6,i,1) = rw0 else qt = qte w(1,i,1) = . (w(1,i,2)**gamma*c(i,1)**2/p(i,2)/gamma)**(1./gam1) w(5,i,1) = w(5,i,2) w(6,i,1) = w(6,i,2) endif p(i,1) = w(1,i,1)*c(i,1)**2/gamma w(2,i,1) = w(1,i,1)*(qn*ex - qt*ey) w(3,i,1) = w(1,i,1)*(qn*ey + qt*ex) w(4,i,1) = p(i,1)/gam1 . +0.5*(w(2,i,1)**2+w(3,i,1)**2)/w(1,i,1) vis(i,1) = vis(i,2) ff(i,1) = ff(i,2) enddo return end subroutine JmaxRiemann include 'mpins2d.inc' do i = 2, nfx ip = i + 1 ex = y(i,nfy+1)-y(ip,nfy+1) ey = x(ip,nfy+1)-x(i,nfy+1) dd = sqrt(ex**2+ey**2) ex = ex/dd ey = ey/dd qne = (w(2,i,nfy)*ex+w(3,i,nfy)*ey)/w(1,i,nfy) qte = (-w(2,i,nfy)*ey+w(3,i,nfy)*ex)/w(1,i,nfy) qnf = u0*ex+v0*ey qtf = -u0*ey+v0*ex R2 = qne+2*c(i,nfy)/gam1 R1 = qnf-2*c0/gam1 qn = 0.5*(R1+R2) c(i,nfy+1) = 0.25*gam1*(R2-R1) if (qn.le.0) then qt = qtf w(1,i,nfy+1) = . (rho0**gamma*c(i,nfy+1)**2/p0/gamma)**(1./gam1) w(5,i,nfy+1) = rk0 w(6,i,nfy+1) = rw0 else qt = qte w(1,i,nfy+1) = . (w(1,i,nfy)**gamma*c(i,nfy+1)**2/p(i,nfy)/gamma)**(1./gam1) w(5,i,nfy+1) = w(5,i,nfy) w(6,i,nfy+1) = w(6,i,nfy) endif p(i,nfy+1) = w(1,i,nfy+1)*c(i,nfy+1)**2/gamma w(2,i,nfy+1) = w(1,i,nfy+1)*(qn*ex - qt*ey) w(3,i,nfy+1) = w(1,i,nfy+1)*(qn*ey + qt*ex) w(4,i,nfy+1) = p(i,nfy+1)/gam1 . +0.5*(w(2,i,nfy+1)**2+w(3,i,nfy+1)**2)/w(1,i,nfy+1) vis(i,nfy+1) = vis(i,nfy) ff(i,nfy+1) = ff(i,nfy) enddo return end subroutine hole_check include 'mpins2d.inc' include 'mpif.h' call MPI_TYPE_VECTOR(nfy-1,1,imx,MPI_INTEGER,MPI_INTEGER_ROW,ierr) call MPI_TYPE_VECTOR(nfy,1,imx,MPI_REAL,MPI_REAL_ROW,ierr) call MPI_TYPE_COMMIT(MPI_INTEGER_ROW,ierr) call MPI_TYPE_COMMIT(MPI_REAL_ROW,ierr) if (jmPE.ge.0) then call MPI_SEND(ib(2,2),nfx-1,MPI_INTEGER,jmPE,2211, . MPI_COMM_WORLD,ierr) call MPI_SEND(x(2,3),nfx,MPI_REAL,jmPE,12211, . MPI_COMM_WORLD,ierr) call MPI_SEND(y(2,3),nfx,MPI_REAL,jmPE,22211, . MPI_COMM_WORLD,ierr) else do i = 2, nfx if (ib(i,2).eq.0) ib(i,1) = 0 enddo endif if (jpPE.ge.0) then call MPI_RECV(ib(2,nfy+1),nfx-1,MPI_INTEGER,jpPE,2211, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(x(2,nfy+2),nfx,MPI_REAL,jpPE,12211, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(y(2,nfy+2),nfx,MPI_REAL,jpPE,22211, . MPI_COMM_WORLD,istatus,ierr) call MPI_SEND(ib(2,nfy),nfx-1,MPI_INTEGER,jpPE,2221, . MPI_COMM_WORLD,ierr) call MPI_SEND(x(2,nfy),nfx,MPI_REAL,jpPE,12221, . MPI_COMM_WORLD,ierr) call MPI_SEND(y(2,nfy),nfx,MPI_REAL,jpPE,22221, . MPI_COMM_WORLD,ierr) else do i = 2, nfx if (ib(i,nfy).eq.0) ib(i,nfy+1) = 0 enddo endif if (jmPE.ge.0) then call MPI_RECV(ib(2,1),nfx-1,MPI_INTEGER,jmPE,2221, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(x(2,1),nfx,MPI_REAL,jmPE,12221, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(y(2,1),nfx,MPI_REAL,jmPE,22221, . MPI_COMM_WORLD,istatus,ierr) endif if (imPE.ge.0) then call MPI_SEND(ib(2,2),1,MPI_INTEGER_ROW,imPE,2311, . MPI_COMM_WORLD,ierr) call MPI_SEND(x(3,2),1,MPI_REAL_ROW,imPE,12311, . MPI_COMM_WORLD,ierr) call MPI_SEND(y(3,2),1,MPI_REAL_ROW,imPE,22311, . MPI_COMM_WORLD,ierr) else do j = 2, nfy if (ib(2,j).eq.0) ib(1,j) = 0 enddo endif if (ipPE.ge.0) then call MPI_RECV(ib(nfx+1,2),1,MPI_INTEGER_ROW,ipPE,2311, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(x(nfx+2,2),1,MPI_REAL_ROW,ipPE,12311, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(y(nfx+2,2),1,MPI_REAL_ROW,ipPE,22311, . MPI_COMM_WORLD,istatus,ierr) call MPI_SEND(ib(nfx,2),1,MPI_INTEGER_ROW,ipPE,2321, . MPI_COMM_WORLD,ierr) call MPI_SEND(x(nfx,2),1,MPI_REAL_ROW,ipPE,12321, . MPI_COMM_WORLD,ierr) call MPI_SEND(y(nfx,2),1,MPI_REAL_ROW,ipPE,22321, . MPI_COMM_WORLD,ierr) else do j = 2, nfy if (ib(nfx,j).eq.0) ib(nfx+1,j) = 0 enddo endif if (imPE.ge.0) then call MPI_RECV(ib(1,2),1,MPI_INTEGER_ROW,imPE,2321, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(x(1,2),1,MPI_REAL_ROW,imPE,12321, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(y(1,2),1,MPI_REAL_ROW,imPE,22321, . MPI_COMM_WORLD,istatus,ierr) endif call MPI_TYPE_FREE(MPI_INTEGER_ROW,ierr) call MPI_TYPE_FREE(MPI_REAL_ROW,ierr) return end subroutine muscl(qll,ql,qr,qrr,ib1,ib2) capa = 1./3. dql = ql - qll dqc = qr - ql dqr = qrr - qr ss = (2.*dql*dqc+1.e-6)/(dqc**2+dql**2+1.e-6) omks = 1.0 - capa*ss opks = 1.0 + capa*ss ql = ql + 0.25*ss*(omks*dql+opks*dqc)*ib1*ib2 ss = (2.*dqc*dqr+1.e-6)/(dqr**2+dqc**2+1.e-6) omks = 1.0 - capa*ss opks = 1.0 + capa*ss qr = qr - 0.25*ss*(omks*dqr+opks*dqc)*ib1*ib2 return end subroutine grid include 'mpins2d.inc' dimension yt(imx,jmx) read(13,*) nfxc, nfyc, i_hole if(nx.ne.nfxc.or.ny.ne.nfyc) then call MPI_FINALIZE(ierr) endif do n = 1, i_hole read(13,*) i_hole_start(n), j_hole_start(n), . i_hole_end(n), j_hole_end(n) enddo do i = 2, nx+1 do j = 2, ny+1 read(13,*) gx(i,j), gy(i,j) enddo enddo call sstgrid(yt) do j = 2, ny do i = 2, nx yyn(i,j) = 0.25*(yt(i,j)+yt(i+1,j)+yt(i,j+1)+yt(i+1,j+1)) enddo enddo do j = 1, ny+1 do i = 1, nx+1 nb(i,j) = 1 enddo enddo do n = 1, i_hole do i = i_hole_start(n)+1, i_hole_end(n) do j = j_hole_start(n)+1, j_hole_end(n) nb(i,j) = 0 enddo enddo enddo return end subroutine F1calc include 'mpins2d.inc' do j = 2, nfy jp = j + 1 jm = j - 1 do i = 2, nfx ip = i + 1 im = i - 1 uch = 0.5*(w(2,ip,j)/w(1,ip,j)-w(2,im,j)/w(1,im,j)) uet = 0.5*(w(2,i,jp)/w(1,i,jp)-w(2,i,jm)/w(1,i,jm)) vch = 0.5*(w(3,ip,j)/w(1,ip,j)-w(3,im,j)/w(1,im,j)) vet = 0.5*(w(3,i,jp)/w(1,i,jp)-w(3,i,jm)/w(1,i,jm)) rkch = 0.5*(w(5,ip,j)-w(5,im,j)) rket = 0.5*(w(5,i,jp)-w(5,i,jm)) rwch = 0.5*(w(6,ip,j)-w(6,im,j)) rwet = 0.5*(w(6,i,jp)-w(6,i,jm)) uy = uch*chy(i,j)+uet*ety(i,j) vx = vch*chx(i,j)+vet*etx(i,j) rkx = rkch*chx(i,j)+rket*etx(i,j) rky = rkch*chy(i,j)+rket*ety(i,j) rwx = rwch*chx(i,j)+rwet*etx(i,j) rwy = rwch*chy(i,j)+rwet*ety(i,j) CDkw = max(2*w(1,i,j)*sigmao2*(rkx*rwx+rky*rwy) . /w(6,i,j),1.e-20) term1 = sqrt(w(5,i,j))/(0.09*w(6,i,j)*yn(i,j)) t = c(i,j)**2 rmu = t**1.5*(1.+tcoe)/(t+tcoe) term2 = 500.*rmu/(yn(i,j)**2*w(6,i,j)*ra*w(1,i,j)) term3 = 4*w(1,i,j)*sigmao2*w(5,i,j)/(CDkw*yn(i,j)**2) arg1 = min(max(term1,term2),term3) ff(i,j) = tanh(arg1**4) arg2 = max(2*term1,term2) f2 = tanh(arg2**2) vor = abs(uy-vx) vis(i,j) = aa1*w(5,i,j)*w(1,i,j)/max(aa1*w(6,i,j),vor*f2) enddo enddo return end subroutine visbc include 'mpins2d.inc' include 'mpif.h' call MPI_TYPE_VECTOR(nfy+1,1,imx,MPI_REAL,MPI_REAL_ROW,ierr) call MPI_TYPE_COMMIT(MPI_REAL_ROW,ierr) if (jmPE.ge.0) then call MPI_SEND(vis(1,2),nfx+1,MPI_REAL,jmPE,9211, . MPI_COMM_WORLD,ierr) call MPI_SEND(ff(1,2),nfx+1,MPI_REAL,jmPE,9212, . MPI_COMM_WORLD,ierr) endif if (jpPE.ge.0) then call MPI_RECV(vis(1,nfy+1),nfx+1,MPI_REAL,jpPE,9211, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(ff(1,nfy+1),nfx+1,MPI_REAL,jpPE,9212, . MPI_COMM_WORLD,istatus,ierr) call MPI_SEND(vis(1,nfy),nfx+1,MPI_REAL,jpPE,9221, . MPI_COMM_WORLD,ierr) call MPI_SEND(ff(1,nfy),nfx+1,MPI_REAL,jpPE,9222, . MPI_COMM_WORLD,ierr) endif if (jmPE.ge.0) then call MPI_RECV(vis(1,1),nfx+1,MPI_REAL,jmPE,9221, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(ff(1,1),nfx+1,MPI_REAL,jmPE,9222, . MPI_COMM_WORLD,istatus,ierr) endif if (imPE.ge.0) then call MPI_SEND(vis(2,1),1,MPI_REAL_ROW,imPE,9311, . MPI_COMM_WORLD,ierr) call MPI_SEND(ff(2,1),1,MPI_REAL_ROW,imPE,9312, . MPI_COMM_WORLD,ierr) endif if (ipPE.ge.0) then call MPI_RECV(vis(nfx+1,1),1,MPI_REAL_ROW,ipPE,9311, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(ff(nfx+1,1),1,MPI_REAL_ROW,ipPE,9312, . MPI_COMM_WORLD,istatus,ierr) call MPI_SEND(vis(nfx,1),1,MPI_REAL_ROW,ipPE,9321, . MPI_COMM_WORLD,ierr) call MPI_SEND(ff(nfx,1),1,MPI_REAL_ROW,ipPE,9322, . MPI_COMM_WORLD,ierr) endif if (imPE.ge.0) then call MPI_RECV(vis(1,1),1,MPI_REAL_ROW,imPE,9321, . MPI_COMM_WORLD,istatus,ierr) call MPI_RECV(ff(1,1),1,MPI_REAL_ROW,imPE,9322, . MPI_COMM_WORLD,istatus,ierr) endif call MPI_TYPE_FREE(MPI_REAL_ROW,ierr) return end subroutine sstgrid(yt) include 'mpins2d.inc' dimension yt(imx,jmx),mb(imx,jmx) do j = 2, ny+1 do i = 2, nx+1 mb(i,j) = 1 enddo enddo if (iminbc.eq.WALL) then do j = 2, ny+1 mb(2,j) = 0 enddo endif if (imaxbc.eq.WALL) then do j = 2, ny+1 mb(nx+1,j) = 0 enddo endif if (jminbc.eq.WALL) then do i = 2, nx+1 mb(i,2) = 0 enddo endif if (jmaxbc.eq.WALL) then do i = 2, nx+1 mb(i,ny+1) = 0 enddo endif do n = 1, i_hole do i = i_hole_start(n)+2,i_hole_end(n) do j = j_hole_start(n)+2,j_hole_end(n) mb(i,j) = -1 enddo enddo do j = j_hole_start(n)+1,j_hole_end(n)+1 mb(i_hole_start(n)+1,j) = 0 mb(i_hole_end(n)+1,j) = 0 enddo do i = i_hole_start(n)+1,i_hole_end(n)+1 mb(i,j_hole_start(n)+1) = 0 mb(i,j_hole_end(n)+1) = 0 enddo enddo do j = 2, ny+1 if (mb(3,j).eq.-1) mb(2,j) = -1 if (mb(nx,j).eq.-1) mb(nx+1,j) = -1 enddo do i = 2, nx+1 if (mb(i,3).eq.-1) mb(i,2) = -1 if (mb(i,ny).eq.-1) mb(i,ny+1) = -1 enddo do i = 3, nx im = i - 1 ip = i + 1 do j = 2, ny+1 if ((mb(i,j).eq.0).and.(mb(im,j).eq.-1).and.(mb(ip,j).eq.-1)) . mb(i,j) = -1 enddo enddo do j = 3, ny jm = j - 1 jp = j + 1 do i = 2, nx+1 if ((mb(i,j).eq.0).and.(mb(i,jm).eq.-1).and.(mb(i,jp).eq.-1)) . mb(i,j) = -1 enddo enddo do j = 2, ny+1 do i = 2, nx+1 yt(i,j) = 1.e10 do n = 2, ny+1 do m = 2, nx+1 if ((mb(i,j).eq.1).and.(mb(m,n).eq.0)) then ytmp = sqrt((gx(i,j)-gx(m,n))**2+(gy(i,j)-gy(m,n))**2) if (ytmp.le.yt(i,j)) yt(i,j) = ytmp endif if (mb(i,j).eq.0) yt(i,j) = 0.0 if (mb(i,j).eq.-1) yt(i,j) = -1 enddo enddo enddo enddo return end --------------BB2C3A587E7F9B68300F0821--