From 3ae844b3312ba8871775c0082f9315ed82a400e3 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Fri, 29 Dec 2023 16:34:19 -0500 Subject: [PATCH 1/6] fix continues --- mfc.sh | 4 +- src/common/m_eigen_solver.f90 | 854 +++++++++++++++++----------------- 2 files changed, 429 insertions(+), 429 deletions(-) diff --git a/mfc.sh b/mfc.sh index 542acd863..ff52cd176 100755 --- a/mfc.sh +++ b/mfc.sh @@ -161,8 +161,8 @@ elif [ "$1" == "format" ]; then fprettify ${@:-src} --exclude "src/*/autogen" --recursive --silent \ --indent 4 --c-relations --enable-replacements --enable-decl \ - --whitespace-comma 1 --whitespace-multdiv 1 --whitespace-plusminus 1 \ - --case 1 1 1 1 --strict-indent + --whitespace-comma 0 --whitespace-multdiv 0 --whitespace-plusminus 1 \ + --case 1 1 1 1 --strict-indent --line-length 1000 ret="$?" if [ "$ret" != '0' ]; then diff --git a/src/common/m_eigen_solver.f90 b/src/common/m_eigen_solver.f90 index 03a6d88f7..b255893fa 100644 --- a/src/common/m_eigen_solver.f90 +++ b/src/common/m_eigen_solver.f90 @@ -11,7 +11,7 @@ module m_eigen_solver implicit none private; public :: cg,cbal,corth,comqr2,csroot,cdiv,pythag - + contains subroutine cg(nm,nl,ar,ai,wr,wi,zr,zi,fv1,fv2,fv3,ierr) @@ -50,22 +50,22 @@ subroutine cg(nm,nl,ar,ai,wr,wi,zr,zi,fv1,fv2,fv3,ierr) ! ! this version dated august 1983. ! -! ------------------------------------------------------------------ +! ------------------------------------------------------------------ integer nm,nl,is1,is2,ierr - real(kind(0d0)), dimension(nm,nl) :: ar,ai,zr,zi - real(kind(0d0)), dimension(nl) :: wr,wi,fv1,fv2,fv3 + real(kind(0d0)),dimension(nm,nl) :: ar,ai,zr,zi + real(kind(0d0)),dimension(nl) :: wr,wi,fv1,fv2,fv3 - if (nl .le. nm) go to 10 + if (nl <= nm) go to 10 ierr = 10*nl go to 50 10 call cbal(nm,nl,ar,ai,is1,is2,fv1) call corth(nm,nl,is1,is2,ar,ai,fv2,fv3) call comqr2(nm,nl,is1,is2,fv2,fv3,ar,ai,wr,wi,zr,zi,ierr) - if (ierr .ne. 0) go to 50 + if (ierr /= 0) go to 50 call cbabk2(nm,nl,is1,is2,fv1,nl,zr,zi) -50 return +50 return end subroutine cg subroutine cbal(nm,nl,ar,ai,low,igh,scale) @@ -126,130 +126,130 @@ subroutine cbal(nm,nl,ar,ai,low,igh,scale) ! ! ------------------------------------------------------------------ integer i,j,k,l,ml,nl,jj,nm,igh,low,iexc - real(kind(0d0)), dimension(nm,nl) :: ar,ai - real(kind(0d0)), dimension(nl) :: scale + real(kind(0d0)),dimension(nm,nl) :: ar,ai + real(kind(0d0)),dimension(nl) :: scale real(kind(0d0)) :: c,f,g,r,s,b2,radix logical noconv radix = 16.0d0 - b2 = radix * radix + b2 = radix*radix k = 1 l = nl go to 100 ! .......... in-line procedure for row and ! column exchange .......... 20 scale(ml) = j - if (j .eq. ml) go to 50 + if (j == ml) go to 50 - do 30 i = 1, l + do 30 i = 1,l f = ar(i,j) ar(i,j) = ar(i,ml) ar(i,ml) = f f = ai(i,j) ai(i,j) = ai(i,ml) ai(i,ml) = f -30 continue +30 end do - do 40 i = k, nl + do 40 i = k,nl f = ar(j,i) ar(j,i) = ar(ml,i) ar(ml,i) = f f = ai(j,i) ai(j,i) = ai(ml,i) ai(ml,i) = f -40 continue +40 end do -50 go to (80,130), iexc +50 go to(80,130),iexc ! .......... search for rows isolating an eigenvalue ! and push them down .......... -80 if (l .eq. 1) go to 280 - l = l - 1 +80 if (l == 1) go to 280 + l = l - 1 ! .......... for j=l step -1 until 1 do -- .......... -100 do 120 jj = 1, l +100 do 120 jj = 1,l j = l + 1 - jj - do 110 i = 1, l - if (i .eq. j) go to 110 - if (ar(j,i) .ne. 0.0d0 .or. ai(j,i) .ne. 0.0d0) go to 120 -110 continue + do 110 i = 1,l + if (i == j) go to 110 + if (ar(j,i) /= 0.0d0 .or. ai(j,i) /= 0.0d0) go to 120 +110 end do - ml = l - iexc = 1 - go to 20 -120 continue + ml = l + iexc = 1 + go to 20 +120 end do go to 140 ! .......... search for columns isolating an eigenvalue ! and push them left .......... 130 k = k + 1 -140 do 170 j = k, l +140 do 170 j = k,l - do 150 i = k, l - if (i .eq. j) go to 150 - if (ar(i,j) .ne. 0.0d0 .or. ai(i,j) .ne. 0.0d0) go to 170 -150 continue + do 150 i = k,l + if (i == j) go to 150 + if (ar(i,j) /= 0.0d0 .or. ai(i,j) /= 0.0d0) go to 170 +150 end do - ml = k - iexc = 2 - go to 20 -170 continue + ml = k + iexc = 2 + go to 20 +170 end do ! .......... now balance the submatrix in rows k to l .......... - do 180 i = k, l - scale(i) = 1.0d0 -180 continue + do 180 i = k,l + scale(i) = 1.0d0 +180 end do ! .......... iterative loop for norm reduction .......... 190 noconv = .false. - do 270 i = k, l + do 270 i = k,l c = 0.0d0 r = 0.0d0 - do 200 j = k, l - if (j .eq. i) go to 200 - c = c + dabs(ar(j,i)) + dabs(ai(j,i)) - r = r + dabs(ar(i,j)) + dabs(ai(i,j)) -200 continue + do 200 j = k,l + if (j == i) go to 200 + c = c + dabs(ar(j,i)) + dabs(ai(j,i)) + r = r + dabs(ar(i,j)) + dabs(ai(i,j)) +200 end do ! .......... guard against zero c or r due to underflow .......... - if (c .eq. 0.0d0 .or. r .eq. 0.0d0) go to 270 - g = r / radix - f = 1.0d0 - s = c + r -210 if (c .ge. g) go to 220 - f = f * radix - c = c * b2 - go to 210 -220 g = r * radix -230 if (c .lt. g) go to 240 - f = f / radix - c = c / b2 - go to 230 + if (c == 0.0d0 .or. r == 0.0d0) go to 270 + g = r/radix + f = 1.0d0 + s = c + r +210 if (c >= g) go to 220 + f = f*radix + c = c*b2 + go to 210 +220 g = r*radix +230 if (c < g) go to 240 + f = f/radix + c = c/b2 + go to 230 ! .......... now balance .......... -240 if ((c + r) / f .ge. 0.95d0 * s) go to 270 - g = 1.0d0 / f - scale(i) = scale(i) * f - noconv = .true. +240 if ((c + r)/f >= 0.95d0*s) go to 270 + g = 1.0d0/f + scale(i) = scale(i)*f + noconv = .true. - do 250 j = k, nl - ar(i,j) = ar(i,j) * g - ai(i,j) = ai(i,j) * g -250 continue + do 250 j = k,nl + ar(i,j) = ar(i,j)*g + ai(i,j) = ai(i,j)*g +250 end do - do 260 j = 1, l - ar(j,i) = ar(j,i) * f - ai(j,i) = ai(j,i) * f -260 continue + do 260 j = 1,l + ar(j,i) = ar(j,i)*f + ai(j,i) = ai(j,i)*f +260 end do -270 continue +270 end do if (noconv) go to 190 280 low = k igh = l - return + return end subroutine cbal - + subroutine corth(nm,nl,low,igh,ar,ai,ortr,orti) ! this subroutine is a translation of a complex analogue of ! the algol procedure orthes, num. math. 12, 349-368(1968) @@ -302,91 +302,91 @@ subroutine corth(nm,nl,low,igh,ar,ai,ortr,orti) integer mll mll = 6 - + la = igh - 1 kp1 = low + 1 - if (la .lt. kp1) go to 200 + if (la < kp1) go to 200 - do 180 ml = kp1, la + do 180 ml = kp1,la h = 0.0d0 ortr(ml) = 0.0d0 orti(ml) = 0.0d0 scale = 0.0d0 ! .......... scale column (algol tol then not needed) .......... - do 90 i = ml, igh - scale = scale + dabs(ar(i,ml-1)) + dabs(ai(i,ml-1)) -90 continue - if (scale .eq. 0d0) go to 180 - mp = ml + igh + do 90 i = ml,igh + scale = scale + dabs(ar(i,ml - 1)) + dabs(ai(i,ml - 1)) +90 end do + if (scale == 0d0) go to 180 + mp = ml + igh ! .......... for i=igh step -1 until ml do -- .......... - do 100 ii = ml, igh - i = mp - ii - ortr(i) = ar(i,ml-1) / scale - orti(i) = ai(i,ml-1) / scale - h = h + ortr(i) * ortr(i) + orti(i) * orti(i) -100 continue -! - g = dsqrt(h) - call pythag(ortr(ml),orti(ml),f) - if (f .eq. 0d0) go to 103 - h = h + f * g - g = g / f - ortr(ml) = (1.0d0 + g) * ortr(ml) - orti(ml) = (1.0d0 + g) * orti(ml) - go to 105 - -103 ortr(ml) = g - ar(ml,ml-1) = scale + do 100 ii = ml,igh + i = mp - ii + ortr(i) = ar(i,ml - 1)/scale + orti(i) = ai(i,ml - 1)/scale + h = h + ortr(i)*ortr(i) + orti(i)*orti(i) +100 end do +! + g = dsqrt(h) + call pythag(ortr(ml),orti(ml),f) + if (f == 0d0) go to 103 + h = h + f*g + g = g/f + ortr(ml) = (1.0d0 + g)*ortr(ml) + orti(ml) = (1.0d0 + g)*orti(ml) + go to 105 + +103 ortr(ml) = g + ar(ml,ml - 1) = scale ! .......... form (i-(u*ut)/h) * a .......... -105 do 130 j = ml, nl - fr = 0.0d0 - fi = 0.0d0 +105 do 130 j = ml,nl + fr = 0.0d0 + fi = 0.0d0 ! .......... for i=igh step -1 until ml do -- .......... - do 110 ii = ml, igh - i = mp - ii - fr = fr + ortr(i) * ar(i,j) + orti(i) * ai(i,j) - fi = fi + ortr(i) * ai(i,j) - orti(i) * ar(i,j) -110 continue + do 110 ii = ml,igh + i = mp - ii + fr = fr + ortr(i)*ar(i,j) + orti(i)*ai(i,j) + fi = fi + ortr(i)*ai(i,j) - orti(i)*ar(i,j) +110 end do ! - fr = fr / h - fi = fi / h + fr = fr/h + fi = fi/h ! - do 120 i = ml, igh - ar(i,j) = ar(i,j) - fr * ortr(i) + fi * orti(i) - ai(i,j) = ai(i,j) - fr * orti(i) - fi * ortr(i) -120 continue + do 120 i = ml,igh + ar(i,j) = ar(i,j) - fr*ortr(i) + fi*orti(i) + ai(i,j) = ai(i,j) - fr*orti(i) - fi*ortr(i) +120 end do ! -130 continue +130 end do ! .......... form (i-(u*ut)/h)*a*(i-(u*ut)/h) .......... - do 160 i = 1, igh - fr = 0.0d0 - fi = 0.0d0 + do 160 i = 1,igh + fr = 0.0d0 + fi = 0.0d0 ! .......... for j=igh step -1 until ml do -- .......... - do 140 jj = ml, igh - j = mp - jj - fr = fr + ortr(j) * ar(i,j) - orti(j) * ai(i,j) - fi = fi + ortr(j) * ai(i,j) + orti(j) * ar(i,j) -140 continue + do 140 jj = ml,igh + j = mp - jj + fr = fr + ortr(j)*ar(i,j) - orti(j)*ai(i,j) + fi = fi + ortr(j)*ai(i,j) + orti(j)*ar(i,j) +140 end do ! - fr = fr / h - fi = fi / h + fr = fr/h + fi = fi/h ! - do 150 j = ml, igh - ar(i,j) = ar(i,j) - fr * ortr(j) - fi * orti(j) - ai(i,j) = ai(i,j) + fr * orti(j) - fi * ortr(j) -150 continue + do 150 j = ml,igh + ar(i,j) = ar(i,j) - fr*ortr(j) - fi*orti(j) + ai(i,j) = ai(i,j) + fr*orti(j) - fi*ortr(j) +150 end do ! -160 continue +160 end do ! - ortr(ml) = scale * ortr(ml) - orti(ml) = scale * orti(ml) - ar(ml,ml-1) = -g * ar(ml,ml-1) - ai(ml,ml-1) = -g * ai(ml,ml-1) -180 continue + ortr(ml) = scale*ortr(ml) + orti(ml) = scale*orti(ml) + ar(ml,ml - 1) = -g*ar(ml,ml - 1) + ai(ml,ml - 1) = -g*ai(ml,ml - 1) +180 end do ! -200 return +200 return end subroutine corth - + subroutine comqr2(nm,nl,low,igh,ortr,orti,hr,hi,wr,wi,zr,zi,ierr) ! MESHED overflow control WITH vectors of isolated roots (10/19/89 BSG) ! MESHED overflow control WITH triangular multiply (10/30/89 BSG) @@ -460,147 +460,147 @@ subroutine comqr2(nm,nl,low,igh,ortr,orti,hr,hi,wr,wi,zr,zi,ierr) ! this version dated october 1989. ! ! ------------------------------------------------------------------ - integer i,j,k,l,ml,nl,en,ii,jj,ll,nm,nn,igh,ip1,& - itn,its,low,lp1,enm1,iend,ierr + integer i,j,k,l,ml,nl,en,ii,jj,ll,nm,nn,igh,ip1, & + itn,its,low,lp1,enm1,iend,ierr real(kind(0d0)),dimension(nm,nl) :: hr,hi,zr,zi real(kind(0d0)),dimension(nl) :: wr,wi real(kind(0d0)),dimension(igh) :: ortr,orti - real(kind(0d0)) :: si,sr,ti,tr,xi,xr,yi,yr,zzi,zzr,& - norm,tst1,tst2,c,d + real(kind(0d0)) :: si,sr,ti,tr,xi,xr,yi,yr,zzi,zzr, & + norm,tst1,tst2,c,d ! ierr = 0 ! .......... initialize eigenvector matrix .......... - do 101 j = 1, nl -! - do 100 i = 1, nl - zr(i,j) = 0.0d0 - zi(i,j) = 0.0d0 -100 continue - zr(j,j) = 1.0d0 -101 continue + do 101 j = 1,nl +! + do 100 i = 1,nl + zr(i,j) = 0.0d0 + zi(i,j) = 0.0d0 +100 end do + zr(j,j) = 1.0d0 +101 end do ! .......... form the matrix of accumulated transformations ! from the information left by corth .......... iend = igh - low - 1 - if (iend .lt. 0) go to 180 - if (iend .eq. 0) go to 150 - if (iend .gt. 0) go to 105 + if (iend < 0) go to 180 + if (iend == 0) go to 150 + if (iend > 0) go to 105 ! .......... for i=igh-1 step -1 until low+1 do -- .......... -105 do 140 ii = 1, iend - i = igh - ii - if (dabs(ortr(i)) .eq. 0d0 .and. dabs(orti(i)) .eq. 0d0) go to 140 - if (dabs(hr(i,i-1)) .eq. 0d0 .and. dabs(hi(i,i-1)) .eq. 0d0) go to 140 +105 do 140 ii = 1,iend + i = igh - ii + if (dabs(ortr(i)) == 0d0 .and. dabs(orti(i)) == 0d0) go to 140 + if (dabs(hr(i,i - 1)) == 0d0 .and. dabs(hi(i,i - 1)) == 0d0) go to 140 ! .......... norm below is negative of h formed in corth .......... - norm = hr(i,i-1) * ortr(i) + hi(i,i-1) * orti(i) - ip1 = i + 1 + norm = hr(i,i - 1)*ortr(i) + hi(i,i - 1)*orti(i) + ip1 = i + 1 - do 110 k = ip1, igh - ortr(k) = hr(k,i-1) - orti(k) = hi(k,i-1) -110 continue + do 110 k = ip1,igh + ortr(k) = hr(k,i - 1) + orti(k) = hi(k,i - 1) +110 end do ! - do 130 j = i, igh - sr = 0.0d0 - si = 0.0d0 + do 130 j = i,igh + sr = 0.0d0 + si = 0.0d0 ! - do 115 k = i, igh - sr = sr + ortr(k) * zr(k,j) + orti(k) * zi(k,j) - si = si + ortr(k) * zi(k,j) - orti(k) * zr(k,j) -115 continue + do 115 k = i,igh + sr = sr + ortr(k)*zr(k,j) + orti(k)*zi(k,j) + si = si + ortr(k)*zi(k,j) - orti(k)*zr(k,j) +115 end do ! - sr = sr / norm - si = si / norm + sr = sr/norm + si = si/norm ! - do 120 k = i, igh - zr(k,j) = zr(k,j) + sr * ortr(k) - si * orti(k) - zi(k,j) = zi(k,j) + sr * orti(k) + si * ortr(k) -120 continue + do 120 k = i,igh + zr(k,j) = zr(k,j) + sr*ortr(k) - si*orti(k) + zi(k,j) = zi(k,j) + sr*orti(k) + si*ortr(k) +120 end do ! -130 continue +130 end do ! -140 continue +140 end do ! .......... create real subdiagonal elements .......... 150 l = low + 1 ! - do 170 i = l, igh - ll = min0(i+1,igh) - if (dabs(hi(i,i-1)) .eq. 0d0) go to 170 - call pythag(hr(i,i-1),hi(i,i-1),norm) - yr = hr(i,i-1) / norm - yi = hi(i,i-1) / norm - hr(i,i-1) = norm - hi(i,i-1) = 0.0d0 -! - do 155 j = i, nl - si = yr * hi(i,j) - yi * hr(i,j) - hr(i,j) = yr * hr(i,j) + yi * hi(i,j) - hi(i,j) = si -155 continue -! - do 160 j = 1, ll - si = yr * hi(j,i) + yi * hr(j,i) - hr(j,i) = yr * hr(j,i) - yi * hi(j,i) - hi(j,i) = si -160 continue -! - do 165 j = low, igh - si = yr * zi(j,i) + yi * zr(j,i) - zr(j,i) = yr * zr(j,i) - yi * zi(j,i) - zi(j,i) = si -165 continue - -170 continue + do 170 i = l,igh + ll = min0(i + 1,igh) + if (dabs(hi(i,i - 1)) == 0d0) go to 170 + call pythag(hr(i,i - 1),hi(i,i - 1),norm) + yr = hr(i,i - 1)/norm + yi = hi(i,i - 1)/norm + hr(i,i - 1) = norm + hi(i,i - 1) = 0.0d0 +! + do 155 j = i,nl + si = yr*hi(i,j) - yi*hr(i,j) + hr(i,j) = yr*hr(i,j) + yi*hi(i,j) + hi(i,j) = si +155 end do +! + do 160 j = 1,ll + si = yr*hi(j,i) + yi*hr(j,i) + hr(j,i) = yr*hr(j,i) - yi*hi(j,i) + hi(j,i) = si +160 end do +! + do 165 j = low,igh + si = yr*zi(j,i) + yi*zr(j,i) + zr(j,i) = yr*zr(j,i) - yi*zi(j,i) + zi(j,i) = si +165 end do + +170 end do ! .......... store roots isolated by cbal .......... -180 do 200 i = 1, nl - if (i .ge. low .and. i .le. igh) go to 200 - wr(i) = hr(i,i) - wi(i) = hi(i,i) -200 continue +180 do 200 i = 1,nl + if (i >= low .and. i <= igh) go to 200 + wr(i) = hr(i,i) + wi(i) = hi(i,i) +200 end do ! en = igh tr = 0.0d0 ti = 0.0d0 itn = 30*nl ! .......... search for next eigenvalue .......... -220 if (en .lt. low) go to 680 +220 if (en < low) go to 680 its = 0 enm1 = en - 1 ! .......... look for single small sub-diagonal element ! for l=en step -1 until low do -- .......... -240 do 260 ll = low, en - l = en + low - ll - if (l .eq. low) go to 300 - tst1 = dabs(hr(l-1,l-1)) + dabs(hi(l-1,l-1)) & - + dabs(hr(l,l)) + dabs(hi(l,l)) - tst2 = tst1 + dabs(hr(l,l-1)) - if (tst2 .eq. tst1) go to 300 -260 continue +240 do 260 ll = low,en + l = en + low - ll + if (l == low) go to 300 + tst1 = dabs(hr(l - 1,l - 1)) + dabs(hi(l - 1,l - 1)) & + + dabs(hr(l,l)) + dabs(hi(l,l)) + tst2 = tst1 + dabs(hr(l,l - 1)) + if (tst2 == tst1) go to 300 +260 end do ! .......... form shift .......... -300 if (l .eq. en) go to 660 - if (itn .eq. 0) go to 1000 - if (its .eq. 10 .or. its .eq. 20) go to 320 +300 if (l == en) go to 660 + if (itn == 0) go to 1000 + if (its == 10 .or. its == 20) go to 320 sr = hr(en,en) si = hi(en,en) - xr = hr(enm1,en) * hr(en,enm1) - xi = hi(enm1,en) * hr(en,enm1) - if (xr .eq. 0.0d0 .and. xi .eq. 0.0d0) go to 340 - yr = (hr(enm1,enm1) - sr) / 2.0d0 - yi = (hi(enm1,enm1) - si) / 2.0d0 - call csroot(yr**2-yi**2+xr,2.0d0*yr*yi+xi,zzr,zzi) - if (yr * zzr + yi * zzi .ge. 0.0d0) go to 310 + xr = hr(enm1,en)*hr(en,enm1) + xi = hi(enm1,en)*hr(en,enm1) + if (xr == 0.0d0 .and. xi == 0.0d0) go to 340 + yr = (hr(enm1,enm1) - sr)/2.0d0 + yi = (hi(enm1,enm1) - si)/2.0d0 + call csroot(yr**2 - yi**2 + xr,2.0d0*yr*yi + xi,zzr,zzi) + if (yr*zzr + yi*zzi >= 0.0d0) go to 310 zzr = -zzr zzi = -zzi -310 call cdiv(xr,xi,yr+zzr,yi+zzi,xr,xi) +310 call cdiv(xr,xi,yr + zzr,yi + zzi,xr,xi) sr = sr - xr si = si - xi go to 340 ! .......... form exceptional shift .......... -320 sr = dabs(hr(en,enm1)) + dabs(hr(enm1,en-2)) +320 sr = dabs(hr(en,enm1)) + dabs(hr(enm1,en - 2)) si = 0.0d0 ! -340 do 360 i = low, en - hr(i,i) = hr(i,i) - sr - hi(i,i) = hi(i,i) - si -360 continue +340 do 360 i = low,en + hr(i,i) = hr(i,i) - sr + hi(i,i) = hi(i,i) - si +360 end do ! tr = tr + sr ti = ti + si @@ -609,94 +609,94 @@ subroutine comqr2(nm,nl,low,igh,ortr,orti,hr,hi,wr,wi,zr,zi,ierr) ! .......... reduce to triangle (rows) .......... lp1 = l + 1 ! - do 500 i = lp1, en - sr = hr(i,i-1) - hr(i,i-1) = 0.0d0 - call pythag(hr(i-1,i-1),hi(i-1,i-1),c) - call pythag(c,sr,norm) - xr = hr(i-1,i-1) / norm - wr(i-1) = xr - xi = hi(i-1,i-1) / norm - wi(i-1) = xi - hr(i-1,i-1) = norm - hi(i-1,i-1) = 0.0d0 - hi(i,i-1) = sr / norm -! - do 490 j = i, nl - yr = hr(i-1,j) - yi = hi(i-1,j) - zzr = hr(i,j) - zzi = hi(i,j) - hr(i-1,j) = xr * yr + xi * yi + hi(i,i-1) * zzr - hi(i-1,j) = xr * yi - xi * yr + hi(i,i-1) * zzi - hr(i,j) = xr * zzr - xi * zzi - hi(i,i-1) * yr - hi(i,j) = xr * zzi + xi * zzr - hi(i,i-1) * yi -490 continue -! -500 continue + do 500 i = lp1,en + sr = hr(i,i - 1) + hr(i,i - 1) = 0.0d0 + call pythag(hr(i - 1,i - 1),hi(i - 1,i - 1),c) + call pythag(c,sr,norm) + xr = hr(i - 1,i - 1)/norm + wr(i - 1) = xr + xi = hi(i - 1,i - 1)/norm + wi(i - 1) = xi + hr(i - 1,i - 1) = norm + hi(i - 1,i - 1) = 0.0d0 + hi(i,i - 1) = sr/norm +! + do 490 j = i,nl + yr = hr(i - 1,j) + yi = hi(i - 1,j) + zzr = hr(i,j) + zzi = hi(i,j) + hr(i - 1,j) = xr*yr + xi*yi + hi(i,i - 1)*zzr + hi(i - 1,j) = xr*yi - xi*yr + hi(i,i - 1)*zzi + hr(i,j) = xr*zzr - xi*zzi - hi(i,i - 1)*yr + hi(i,j) = xr*zzi + xi*zzr - hi(i,i - 1)*yi +490 end do +! +500 end do ! si = hi(en,en) - if (dabs(si) .eq. 0d0) go to 540 + if (dabs(si) == 0d0) go to 540 call pythag(hr(en,en),si,norm) - sr = hr(en,en) / norm - si = si / norm + sr = hr(en,en)/norm + si = si/norm hr(en,en) = norm hi(en,en) = 0.0d0 - if (en .eq. nl) go to 540 + if (en == nl) go to 540 ip1 = en + 1 ! - do 520 j = ip1, nl - yr = hr(en,j) - yi = hi(en,j) - hr(en,j) = sr * yr + si * yi - hi(en,j) = sr * yi - si * yr -520 continue + do 520 j = ip1,nl + yr = hr(en,j) + yi = hi(en,j) + hr(en,j) = sr*yr + si*yi + hi(en,j) = sr*yi - si*yr +520 end do ! .......... inverse operation (columns) .......... -540 do 600 j = lp1, en - xr = wr(j-1) - xi = wi(j-1) -! - do 580 i = 1, j - yr = hr(i,j-1) - yi = 0.0d0 - zzr = hr(i,j) - zzi = hi(i,j) - if (i .eq. j) go to 560 - yi = hi(i,j-1) - hi(i,j-1) = xr * yi + xi * yr + hi(j,j-1) * zzi -560 hr(i,j-1) = xr * yr - xi * yi + hi(j,j-1) * zzr - hr(i,j) = xr * zzr + xi * zzi - hi(j,j-1) * yr - hi(i,j) = xr * zzi - xi * zzr - hi(j,j-1) * yi -580 continue -! - do 590 i = low, igh - yr = zr(i,j-1) - yi = zi(i,j-1) - zzr = zr(i,j) - zzi = zi(i,j) - zr(i,j-1) = xr * yr - xi * yi + hi(j,j-1) * zzr - zi(i,j-1) = xr * yi + xi * yr + hi(j,j-1) * zzi - zr(i,j) = xr * zzr + xi * zzi - hi(j,j-1) * yr - zi(i,j) = xr * zzi - xi * zzr - hi(j,j-1) * yi -590 continue - -600 continue -! - if (dabs(si) .eq. 0d0) go to 240 -! - do 630 i = 1, en - yr = hr(i,en) - yi = hi(i,en) - hr(i,en) = sr * yr - si * yi - hi(i,en) = sr * yi + si * yr -630 continue -! - do 640 i = low, igh - yr = zr(i,en) - yi = zi(i,en) - zr(i,en) = sr * yr - si * yi - zi(i,en) = sr * yi + si * yr -640 continue +540 do 600 j = lp1,en + xr = wr(j - 1) + xi = wi(j - 1) +! + do 580 i = 1,j + yr = hr(i,j - 1) + yi = 0.0d0 + zzr = hr(i,j) + zzi = hi(i,j) + if (i == j) go to 560 + yi = hi(i,j - 1) + hi(i,j - 1) = xr*yi + xi*yr + hi(j,j - 1)*zzi +560 hr(i,j - 1) = xr*yr - xi*yi + hi(j,j - 1)*zzr + hr(i,j) = xr*zzr + xi*zzi - hi(j,j - 1)*yr + hi(i,j) = xr*zzi - xi*zzr - hi(j,j - 1)*yi +580 end do +! + do 590 i = low,igh + yr = zr(i,j - 1) + yi = zi(i,j - 1) + zzr = zr(i,j) + zzi = zi(i,j) + zr(i,j - 1) = xr*yr - xi*yi + hi(j,j - 1)*zzr + zi(i,j - 1) = xr*yi + xi*yr + hi(j,j - 1)*zzi + zr(i,j) = xr*zzr + xi*zzi - hi(j,j - 1)*yr + zi(i,j) = xr*zzi - xi*zzr - hi(j,j - 1)*yi +590 end do + +600 end do +! + if (dabs(si) == 0d0) go to 240 +! + do 630 i = 1,en + yr = hr(i,en) + yi = hi(i,en) + hr(i,en) = sr*yr - si*yi + hi(i,en) = sr*yi + si*yr +630 end do +! + do 640 i = low,igh + yr = zr(i,en) + yi = zi(i,en) + zr(i,en) = sr*yr - si*yi + zi(i,en) = sr*yi + si*yr +640 end do ! go to 240 ! .......... a root found .......... @@ -710,88 +710,88 @@ subroutine comqr2(nm,nl,low,igh,ortr,orti,hr,hi,wr,wi,zr,zi,ierr) ! vectors of upper triangular form .......... 680 norm = 0.0d0 ! - do i = 1, nl - do j = i, nl - tr = dabs(hr(i,j)) + dabs(hi(i,j)) - if (tr .gt. norm) norm = tr - end do + do i = 1,nl + do j = i,nl + tr = dabs(hr(i,j)) + dabs(hi(i,j)) + if (tr > norm) norm = tr + end do end do ! - if (nl .eq. 1 .or. norm .eq. 0d0) go to 1001 + if (nl == 1 .or. norm == 0d0) go to 1001 ! .......... for en=nl step -1 until 2 do -- .......... - do 800 nn = 2, nl - en = nl + 2 - nn - xr = wr(en) - xi = wi(en) - hr(en,en) = 1.0d0 - hi(en,en) = 0.0d0 - enm1 = en - 1 + do 800 nn = 2,nl + en = nl + 2 - nn + xr = wr(en) + xi = wi(en) + hr(en,en) = 1.0d0 + hi(en,en) = 0.0d0 + enm1 = en - 1 ! .......... for i=en-1 step -1 until 1 do -- .......... - do 780 ii = 1, enm1 - i = en - ii - zzr = 0.0d0 - zzi = 0.0d0 - ip1 = i + 1 - - do 740 j = ip1, en - zzr = zzr + hr(i,j) * hr(j,en) - hi(i,j) * hi(j,en) - zzi = zzi + hr(i,j) * hi(j,en) + hi(i,j) * hr(j,en) -740 continue -! - yr = xr - wr(i) - yi = xi - wi(i) - if (yr .ne. 0.0d0 .or. yi .ne. 0.0d0) go to 765 - tst1 = norm - yr = tst1 -760 yr = 0.01d0 * yr - tst2 = norm + yr - if (tst2 .gt. tst1) go to 760 -765 continue - call cdiv(zzr,zzi,yr,yi,hr(i,en),hi(i,en)) + do 780 ii = 1,enm1 + i = en - ii + zzr = 0.0d0 + zzi = 0.0d0 + ip1 = i + 1 + + do 740 j = ip1,en + zzr = zzr + hr(i,j)*hr(j,en) - hi(i,j)*hi(j,en) + zzi = zzi + hr(i,j)*hi(j,en) + hi(i,j)*hr(j,en) +740 end do +! + yr = xr - wr(i) + yi = xi - wi(i) + if (yr /= 0.0d0 .or. yi /= 0.0d0) go to 765 + tst1 = norm + yr = tst1 +760 yr = 0.01d0*yr + tst2 = norm + yr + if (tst2 > tst1) go to 760 +765 continue + call cdiv(zzr,zzi,yr,yi,hr(i,en),hi(i,en)) ! .......... overflow control .......... - tr = dabs(hr(i,en)) + dabs(hi(i,en)) - if (tr .eq. 0.0d0) go to 780 - tst1 = tr - tst2 = tst1 + 1.0d0/tst1 - if (tst2 .gt. tst1) go to 780 - do 770 j = i, en - hr(j,en) = hr(j,en)/tr - hi(j,en) = hi(j,en)/tr -770 continue -! -780 continue -! -800 continue + tr = dabs(hr(i,en)) + dabs(hi(i,en)) + if (tr == 0.0d0) go to 780 + tst1 = tr + tst2 = tst1 + 1.0d0/tst1 + if (tst2 > tst1) go to 780 + do 770 j = i,en + hr(j,en) = hr(j,en)/tr + hi(j,en) = hi(j,en)/tr +770 end do +! +780 end do +! +800 end do ! .......... end backsubstitution .......... ! .......... vectors of isolated roots .......... - do 840 i = 1, nl - if (i .ge. low .and. i .le. igh) go to 840 + do 840 i = 1,nl + if (i >= low .and. i <= igh) go to 840 ! - do 820 j = I, nl - zr(i,j) = hr(i,j) - zi(i,j) = hi(i,j) -820 continue + do 820 j = I,nl + zr(i,j) = hr(i,j) + zi(i,j) = hi(i,j) +820 end do ! -840 continue +840 end do ! .......... multiply by transformation matrix to give ! vectors of original full matrix. ! for j=nl step -1 until low do -- .......... - do jj = low, nl - j = nl + low - jj - ml = min0(j,igh) -! - do i = low, igh - zzr = 0.0d0 - zzi = 0.0d0 -! - do 860 k = low, ml - zzr = zzr + zr(i,k) * hr(k,j) - zi(i,k) * hi(k,j) - zzi = zzi + zr(i,k) * hi(k,j) + zi(i,k) * hr(k,j) -860 continue -! - zr(i,j) = zzr - zi(i,j) = zzi - end do + do jj = low,nl + j = nl + low - jj + ml = min0(j,igh) +! + do i = low,igh + zzr = 0.0d0 + zzi = 0.0d0 +! + do 860 k = low,ml + zzr = zzr + zr(i,k)*hr(k,j) - zi(i,k)*hi(k,j) + zzi = zzi + zr(i,k)*hi(k,j) + zi(i,k)*hr(k,j) +860 end do +! + zr(i,j) = zzr + zi(i,j) = zzi + end do end do ! go to 1001 @@ -843,51 +843,51 @@ subroutine cbabk2(nm,nl,low,igh,scale,ml,zr,zi) ! ! ------------------------------------------------------------------ ! - integer i,j,k,ml,nl,ii,nm,igh,low - double precision scale(nl),zr(nm,ml),zi(nm,ml) - double precision s + integer i,j,k,ml,nl,ii,nm,igh,low + double precision scale(nl),zr(nm,ml),zi(nm,ml) + double precision s - if (ml .eq. 0) go to 200 - if (igh .eq. low) go to 120 + if (ml == 0) go to 200 + if (igh == low) go to 120 ! - do 110 i = low, igh - s = scale(i) + do 110 i = low,igh + s = scale(i) ! .......... left hand eigenvectors are back transformed ! if the foregoing statement is replaced by ! s=1.0d0/scale(i). .......... - do 100 j = 1, ml - zr(i,j) = zr(i,j) * s - zi(i,j) = zi(i,j) * s -100 continue + do 100 j = 1,ml + zr(i,j) = zr(i,j)*s + zi(i,j) = zi(i,j)*s +100 end do ! -110 continue +110 end do ! .......... for i=low-1 step -1 until 1, ! igh+1 step 1 until nl do -- .......... -120 do 140 ii = 1, nl - i = ii - if (i .ge. low .and. i .le. igh) go to 140 - if (i .lt. low) i = low - ii - k = scale(i) - if (k .eq. i) go to 140 -! - do 130 j = 1, ml - s = zr(i,j) - zr(i,j) = zr(k,j) - zr(k,j) = s - s = zi(i,j) - zi(i,j) = zi(k,j) - zi(k,j) = s -130 continue -! -140 continue -! -200 return +120 do 140 ii = 1,nl + i = ii + if (i >= low .and. i <= igh) go to 140 + if (i < low) i = low - ii + k = scale(i) + if (k == i) go to 140 +! + do 130 j = 1,ml + s = zr(i,j) + zr(i,j) = zr(k,j) + zr(k,j) = s + s = zi(i,j) + zi(i,j) = zi(k,j) + zi(k,j) = s +130 end do +! +140 end do +! +200 return end subroutine cbabk2 - + subroutine csroot(xr,xi,yr,yi) real(kind(0d0)) :: xr,xi,yr,yi ! -! (yr,yi) = complex dsqrt(xr,xi) +! (yr,yi) = complex dsqrt(xr,xi) ! branch chosen so that yr .ge. 0.0 and sign(yi) .eq. sign(xi) ! real(kind(0d0)) :: s,tr,ti,c @@ -895,14 +895,14 @@ subroutine csroot(xr,xi,yr,yi) ti = xi call pythag(tr,ti,c) s = dsqrt(0.5d0*(c + dabs(tr))) - if (tr .ge. 0.0d0) yr = s - if (ti .lt. 0.0d0) s = -s - if (tr .le. 0.0d0) yi = s - if (tr .lt. 0.0d0) yr = 0.5d0*(ti/yi) - if (tr .gt. 0.0d0) yi = 0.5d0*(ti/yr) + if (tr >= 0.0d0) yr = s + if (ti < 0.0d0) s = -s + if (tr <= 0.0d0) yi = s + if (tr < 0.0d0) yr = 0.5d0*(ti/yi) + if (tr > 0.0d0) yi = 0.5d0*(ti/yr) return end subroutine csroot - + subroutine cdiv(ar,ai,br,bi,cr,ci) real(kind(0d0)) :: ar,ai,br,bi,cr,ci ! @@ -927,18 +927,18 @@ subroutine pythag(a,b,c) ! real(kind(0d0)) :: p,r,s,t,u p = dmax1(dabs(a),dabs(b)) - if (p .eq. 0.0d0) go to 20 + if (p == 0.0d0) go to 20 r = (dmin1(dabs(a),dabs(b))/p)**2 10 continue t = 4.0d0 + r - if (t .eq. 4.0d0) go to 20 + if (t == 4.0d0) go to 20 s = r/t u = 1.0d0 + 2.0d0*s p = u*p - r = (s/u)**2 * r + r = (s/u)**2*r go to 10 20 c = p return end subroutine pythag -end module m_eigen_solver \ No newline at end of file +end module m_eigen_solver From 9741608682431f6d614bcd436226a9e1f6360266 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Fri, 29 Dec 2023 16:59:25 -0500 Subject: [PATCH 2/6] add diff and err message --- mfc.sh | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/mfc.sh b/mfc.sh index ff52cd176..b7ef089a4 100755 --- a/mfc.sh +++ b/mfc.sh @@ -159,10 +159,26 @@ elif [ "$1" == "format" ]; then pip3 install --upgrade fprettify fi + if [ "$1" == "diff" ]; then + shift + out=$(fprettify ${@:-src} --exclude "src/*/autogen" --recursive --silent \ + --indent 4 --c-relations --enable-replacements --enable-decl \ + --whitespace-comma 0 --whitespace-multdiv 0 --whitespace-plusminus 1 \ + --case 1 1 1 1 --strict-indent --line-length 1000 -s -d) + if [ -z "${out}" ]; then + echo "Already pretty!" + exit 0 + else + error "Not pretty, run ./mfc.sh format" + + exit 1 + fi + fi + fprettify ${@:-src} --exclude "src/*/autogen" --recursive --silent \ --indent 4 --c-relations --enable-replacements --enable-decl \ --whitespace-comma 0 --whitespace-multdiv 0 --whitespace-plusminus 1 \ - --case 1 1 1 1 --strict-indent --line-length 1000 + --case 1 1 1 1 --strict-indent --line-length 1000 ret="$?" if [ "$ret" != '0' ]; then From a0c48e4854ae7648e964fcd66ccb32affefec9a5 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Fri, 29 Dec 2023 19:08:18 -0500 Subject: [PATCH 3/6] DRY --- mfc.sh | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/mfc.sh b/mfc.sh index b7ef089a4..9a7fd771d 100755 --- a/mfc.sh +++ b/mfc.sh @@ -159,12 +159,10 @@ elif [ "$1" == "format" ]; then pip3 install --upgrade fprettify fi + opts="--recursive --silent --indent 4 --c-relations --enable-replacements --enable-decl --whitespace-comma 0 --whitespace-multdiv 0 --whitespace-plusminus 1 --case 1 1 1 1 --strict-indent --line-length 1000" if [ "$1" == "diff" ]; then shift - out=$(fprettify ${@:-src} --exclude "src/*/autogen" --recursive --silent \ - --indent 4 --c-relations --enable-replacements --enable-decl \ - --whitespace-comma 0 --whitespace-multdiv 0 --whitespace-plusminus 1 \ - --case 1 1 1 1 --strict-indent --line-length 1000 -s -d) + out=$(fprettify ${@:-src} --exclude "src/*/autogen" $opts -s -d) if [ -z "${out}" ]; then echo "Already pretty!" exit 0 @@ -175,10 +173,7 @@ elif [ "$1" == "format" ]; then fi fi - fprettify ${@:-src} --exclude "src/*/autogen" --recursive --silent \ - --indent 4 --c-relations --enable-replacements --enable-decl \ - --whitespace-comma 0 --whitespace-multdiv 0 --whitespace-plusminus 1 \ - --case 1 1 1 1 --strict-indent --line-length 1000 + fprettify ${@:-src} --exclude "src/*/autogen" $opts ret="$?" if [ "$ret" != '0' ]; then From 154657fac9b617fd5b4da3fb42ec3eb2fafdc806 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Fri, 29 Dec 2023 19:09:56 -0500 Subject: [PATCH 4/6] add CI --- .github/workflows/pretty.yml | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) create mode 100644 .github/workflows/pretty.yml diff --git a/.github/workflows/pretty.yml b/.github/workflows/pretty.yml new file mode 100644 index 000000000..f2c1fd761 --- /dev/null +++ b/.github/workflows/pretty.yml @@ -0,0 +1,19 @@ +name: Pretty + +on: + push: + + pull_request: + + workflow_dispatch: + +jobs: + docs: + name: Code formatting + runs-on: ubuntu-latest + + steps: + - uses: actions/checkout@v3 + + - name: Count + run: ./mfc.sh format diff From 650ce5f2ca9c6da0665de558fbcdcd0421bd2f38 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Sat, 30 Dec 2023 23:53:54 +0000 Subject: [PATCH 5/6] Update run-phoenix-release-gpu.sh --- misc/run-phoenix-release-gpu.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/misc/run-phoenix-release-gpu.sh b/misc/run-phoenix-release-gpu.sh index bb27af6b5..cdd994a50 100644 --- a/misc/run-phoenix-release-gpu.sh +++ b/misc/run-phoenix-release-gpu.sh @@ -19,6 +19,6 @@ set -x gpu_count=$(nvidia-smi -L | wc -l) # number of GPUs on node gpu_ids=$(seq -s ' ' 0 $(($gpu_count-1))) # 0,1,2,...,gpu_count-1 -./mfc.sh test -a -b mpirun -j $(nproc) \ +./mfc.sh test -a -b mpirun -j 1 \ --gpu -g $gpu_ids From c57893e4cfea97af1779281419385ed77a7f31ba Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Fri, 5 Jan 2024 14:28:32 -0500 Subject: [PATCH 6/6] fast forward --- src/common/m_derived_types.fpp | 6 + src/common/m_eigen_solver.f90 | 854 +++++++++++------------ src/common/m_phase_change.fpp | 795 +++++++++++++++++++++ src/common/m_variables_conversion.fpp | 142 ++-- src/post_process/m_global_parameters.fpp | 7 + src/post_process/m_mpi_proxy.fpp | 7 +- src/post_process/m_start_up.f90 | 3 +- src/pre_process/m_assign_variables.f90 | 12 +- src/pre_process/m_checker.f90 | 25 +- src/pre_process/m_data_output.fpp | 7 +- src/pre_process/m_global_parameters.fpp | 14 + src/pre_process/m_mpi_proxy.fpp | 11 +- src/pre_process/m_start_up.fpp | 17 +- src/simulation/m_cbc.fpp | 10 +- src/simulation/m_checker.fpp | 23 + src/simulation/m_data_output.fpp | 43 +- src/simulation/m_global_parameters.fpp | 17 +- src/simulation/m_mpi_proxy.fpp | 8 +- src/simulation/m_riemann_solvers.fpp | 66 +- src/simulation/m_start_up.fpp | 21 +- src/simulation/m_time_steppers.fpp | 20 +- 21 files changed, 1563 insertions(+), 545 deletions(-) create mode 100644 src/common/m_phase_change.fpp diff --git a/src/common/m_derived_types.fpp b/src/common/m_derived_types.fpp index 163bf0392..5d44d12c9 100644 --- a/src/common/m_derived_types.fpp +++ b/src/common/m_derived_types.fpp @@ -152,6 +152,9 @@ module m_derived_types real(kind(0d0)), dimension(num_fluids_max) :: alpha real(kind(0d0)) :: gamma real(kind(0d0)) :: pi_inf !< + real(kind(0d0)) :: cv !< + real(kind(0d0)) :: qv !< + real(kind(0d0)) :: qvp !< !! Primitive variables associated with the patch. In order, these include @@ -178,6 +181,9 @@ module m_derived_types real(kind(0d0)) :: gamma !< Sp. heat ratio real(kind(0d0)) :: pi_inf !< Liquid stiffness real(kind(0d0)), dimension(2) :: Re !< Reynolds number + REAL(KIND(0d0)) :: cv !< heat capacity + REAL(KIND(0d0)) :: qv !< reference energy per unit mass for SGEOS, q (see Le Metayer (2004)) + REAL(KIND(0d0)) :: qvp !< reference entropy per unit mass for SGEOS, q' (see Le Metayer (2004)) real(kind(0d0)) :: mul0 !< Bubble viscosity real(kind(0d0)) :: ss !< Bubble surface tension real(kind(0d0)) :: pv !< Bubble vapour pressure diff --git a/src/common/m_eigen_solver.f90 b/src/common/m_eigen_solver.f90 index b255893fa..03a6d88f7 100644 --- a/src/common/m_eigen_solver.f90 +++ b/src/common/m_eigen_solver.f90 @@ -11,7 +11,7 @@ module m_eigen_solver implicit none private; public :: cg,cbal,corth,comqr2,csroot,cdiv,pythag - + contains subroutine cg(nm,nl,ar,ai,wr,wi,zr,zi,fv1,fv2,fv3,ierr) @@ -50,22 +50,22 @@ subroutine cg(nm,nl,ar,ai,wr,wi,zr,zi,fv1,fv2,fv3,ierr) ! ! this version dated august 1983. ! -! ------------------------------------------------------------------ +! ------------------------------------------------------------------ integer nm,nl,is1,is2,ierr - real(kind(0d0)),dimension(nm,nl) :: ar,ai,zr,zi - real(kind(0d0)),dimension(nl) :: wr,wi,fv1,fv2,fv3 + real(kind(0d0)), dimension(nm,nl) :: ar,ai,zr,zi + real(kind(0d0)), dimension(nl) :: wr,wi,fv1,fv2,fv3 - if (nl <= nm) go to 10 + if (nl .le. nm) go to 10 ierr = 10*nl go to 50 10 call cbal(nm,nl,ar,ai,is1,is2,fv1) call corth(nm,nl,is1,is2,ar,ai,fv2,fv3) call comqr2(nm,nl,is1,is2,fv2,fv3,ar,ai,wr,wi,zr,zi,ierr) - if (ierr /= 0) go to 50 + if (ierr .ne. 0) go to 50 call cbabk2(nm,nl,is1,is2,fv1,nl,zr,zi) -50 return +50 return end subroutine cg subroutine cbal(nm,nl,ar,ai,low,igh,scale) @@ -126,130 +126,130 @@ subroutine cbal(nm,nl,ar,ai,low,igh,scale) ! ! ------------------------------------------------------------------ integer i,j,k,l,ml,nl,jj,nm,igh,low,iexc - real(kind(0d0)),dimension(nm,nl) :: ar,ai - real(kind(0d0)),dimension(nl) :: scale + real(kind(0d0)), dimension(nm,nl) :: ar,ai + real(kind(0d0)), dimension(nl) :: scale real(kind(0d0)) :: c,f,g,r,s,b2,radix logical noconv radix = 16.0d0 - b2 = radix*radix + b2 = radix * radix k = 1 l = nl go to 100 ! .......... in-line procedure for row and ! column exchange .......... 20 scale(ml) = j - if (j == ml) go to 50 + if (j .eq. ml) go to 50 - do 30 i = 1,l + do 30 i = 1, l f = ar(i,j) ar(i,j) = ar(i,ml) ar(i,ml) = f f = ai(i,j) ai(i,j) = ai(i,ml) ai(i,ml) = f -30 end do +30 continue - do 40 i = k,nl + do 40 i = k, nl f = ar(j,i) ar(j,i) = ar(ml,i) ar(ml,i) = f f = ai(j,i) ai(j,i) = ai(ml,i) ai(ml,i) = f -40 end do +40 continue -50 go to(80,130),iexc +50 go to (80,130), iexc ! .......... search for rows isolating an eigenvalue ! and push them down .......... -80 if (l == 1) go to 280 - l = l - 1 +80 if (l .eq. 1) go to 280 + l = l - 1 ! .......... for j=l step -1 until 1 do -- .......... -100 do 120 jj = 1,l +100 do 120 jj = 1, l j = l + 1 - jj - do 110 i = 1,l - if (i == j) go to 110 - if (ar(j,i) /= 0.0d0 .or. ai(j,i) /= 0.0d0) go to 120 -110 end do + do 110 i = 1, l + if (i .eq. j) go to 110 + if (ar(j,i) .ne. 0.0d0 .or. ai(j,i) .ne. 0.0d0) go to 120 +110 continue - ml = l - iexc = 1 - go to 20 -120 end do + ml = l + iexc = 1 + go to 20 +120 continue go to 140 ! .......... search for columns isolating an eigenvalue ! and push them left .......... 130 k = k + 1 -140 do 170 j = k,l +140 do 170 j = k, l - do 150 i = k,l - if (i == j) go to 150 - if (ar(i,j) /= 0.0d0 .or. ai(i,j) /= 0.0d0) go to 170 -150 end do + do 150 i = k, l + if (i .eq. j) go to 150 + if (ar(i,j) .ne. 0.0d0 .or. ai(i,j) .ne. 0.0d0) go to 170 +150 continue - ml = k - iexc = 2 - go to 20 -170 end do + ml = k + iexc = 2 + go to 20 +170 continue ! .......... now balance the submatrix in rows k to l .......... - do 180 i = k,l - scale(i) = 1.0d0 -180 end do + do 180 i = k, l + scale(i) = 1.0d0 +180 continue ! .......... iterative loop for norm reduction .......... 190 noconv = .false. - do 270 i = k,l + do 270 i = k, l c = 0.0d0 r = 0.0d0 - do 200 j = k,l - if (j == i) go to 200 - c = c + dabs(ar(j,i)) + dabs(ai(j,i)) - r = r + dabs(ar(i,j)) + dabs(ai(i,j)) -200 end do + do 200 j = k, l + if (j .eq. i) go to 200 + c = c + dabs(ar(j,i)) + dabs(ai(j,i)) + r = r + dabs(ar(i,j)) + dabs(ai(i,j)) +200 continue ! .......... guard against zero c or r due to underflow .......... - if (c == 0.0d0 .or. r == 0.0d0) go to 270 - g = r/radix - f = 1.0d0 - s = c + r -210 if (c >= g) go to 220 - f = f*radix - c = c*b2 - go to 210 -220 g = r*radix -230 if (c < g) go to 240 - f = f/radix - c = c/b2 - go to 230 + if (c .eq. 0.0d0 .or. r .eq. 0.0d0) go to 270 + g = r / radix + f = 1.0d0 + s = c + r +210 if (c .ge. g) go to 220 + f = f * radix + c = c * b2 + go to 210 +220 g = r * radix +230 if (c .lt. g) go to 240 + f = f / radix + c = c / b2 + go to 230 ! .......... now balance .......... -240 if ((c + r)/f >= 0.95d0*s) go to 270 - g = 1.0d0/f - scale(i) = scale(i)*f - noconv = .true. +240 if ((c + r) / f .ge. 0.95d0 * s) go to 270 + g = 1.0d0 / f + scale(i) = scale(i) * f + noconv = .true. - do 250 j = k,nl - ar(i,j) = ar(i,j)*g - ai(i,j) = ai(i,j)*g -250 end do + do 250 j = k, nl + ar(i,j) = ar(i,j) * g + ai(i,j) = ai(i,j) * g +250 continue - do 260 j = 1,l - ar(j,i) = ar(j,i)*f - ai(j,i) = ai(j,i)*f -260 end do + do 260 j = 1, l + ar(j,i) = ar(j,i) * f + ai(j,i) = ai(j,i) * f +260 continue -270 end do +270 continue if (noconv) go to 190 280 low = k igh = l - return + return end subroutine cbal - + subroutine corth(nm,nl,low,igh,ar,ai,ortr,orti) ! this subroutine is a translation of a complex analogue of ! the algol procedure orthes, num. math. 12, 349-368(1968) @@ -302,91 +302,91 @@ subroutine corth(nm,nl,low,igh,ar,ai,ortr,orti) integer mll mll = 6 - + la = igh - 1 kp1 = low + 1 - if (la < kp1) go to 200 + if (la .lt. kp1) go to 200 - do 180 ml = kp1,la + do 180 ml = kp1, la h = 0.0d0 ortr(ml) = 0.0d0 orti(ml) = 0.0d0 scale = 0.0d0 ! .......... scale column (algol tol then not needed) .......... - do 90 i = ml,igh - scale = scale + dabs(ar(i,ml - 1)) + dabs(ai(i,ml - 1)) -90 end do - if (scale == 0d0) go to 180 - mp = ml + igh + do 90 i = ml, igh + scale = scale + dabs(ar(i,ml-1)) + dabs(ai(i,ml-1)) +90 continue + if (scale .eq. 0d0) go to 180 + mp = ml + igh ! .......... for i=igh step -1 until ml do -- .......... - do 100 ii = ml,igh - i = mp - ii - ortr(i) = ar(i,ml - 1)/scale - orti(i) = ai(i,ml - 1)/scale - h = h + ortr(i)*ortr(i) + orti(i)*orti(i) -100 end do -! - g = dsqrt(h) - call pythag(ortr(ml),orti(ml),f) - if (f == 0d0) go to 103 - h = h + f*g - g = g/f - ortr(ml) = (1.0d0 + g)*ortr(ml) - orti(ml) = (1.0d0 + g)*orti(ml) - go to 105 - -103 ortr(ml) = g - ar(ml,ml - 1) = scale + do 100 ii = ml, igh + i = mp - ii + ortr(i) = ar(i,ml-1) / scale + orti(i) = ai(i,ml-1) / scale + h = h + ortr(i) * ortr(i) + orti(i) * orti(i) +100 continue +! + g = dsqrt(h) + call pythag(ortr(ml),orti(ml),f) + if (f .eq. 0d0) go to 103 + h = h + f * g + g = g / f + ortr(ml) = (1.0d0 + g) * ortr(ml) + orti(ml) = (1.0d0 + g) * orti(ml) + go to 105 + +103 ortr(ml) = g + ar(ml,ml-1) = scale ! .......... form (i-(u*ut)/h) * a .......... -105 do 130 j = ml,nl - fr = 0.0d0 - fi = 0.0d0 +105 do 130 j = ml, nl + fr = 0.0d0 + fi = 0.0d0 ! .......... for i=igh step -1 until ml do -- .......... - do 110 ii = ml,igh - i = mp - ii - fr = fr + ortr(i)*ar(i,j) + orti(i)*ai(i,j) - fi = fi + ortr(i)*ai(i,j) - orti(i)*ar(i,j) -110 end do + do 110 ii = ml, igh + i = mp - ii + fr = fr + ortr(i) * ar(i,j) + orti(i) * ai(i,j) + fi = fi + ortr(i) * ai(i,j) - orti(i) * ar(i,j) +110 continue ! - fr = fr/h - fi = fi/h + fr = fr / h + fi = fi / h ! - do 120 i = ml,igh - ar(i,j) = ar(i,j) - fr*ortr(i) + fi*orti(i) - ai(i,j) = ai(i,j) - fr*orti(i) - fi*ortr(i) -120 end do + do 120 i = ml, igh + ar(i,j) = ar(i,j) - fr * ortr(i) + fi * orti(i) + ai(i,j) = ai(i,j) - fr * orti(i) - fi * ortr(i) +120 continue ! -130 end do +130 continue ! .......... form (i-(u*ut)/h)*a*(i-(u*ut)/h) .......... - do 160 i = 1,igh - fr = 0.0d0 - fi = 0.0d0 + do 160 i = 1, igh + fr = 0.0d0 + fi = 0.0d0 ! .......... for j=igh step -1 until ml do -- .......... - do 140 jj = ml,igh - j = mp - jj - fr = fr + ortr(j)*ar(i,j) - orti(j)*ai(i,j) - fi = fi + ortr(j)*ai(i,j) + orti(j)*ar(i,j) -140 end do + do 140 jj = ml, igh + j = mp - jj + fr = fr + ortr(j) * ar(i,j) - orti(j) * ai(i,j) + fi = fi + ortr(j) * ai(i,j) + orti(j) * ar(i,j) +140 continue ! - fr = fr/h - fi = fi/h + fr = fr / h + fi = fi / h ! - do 150 j = ml,igh - ar(i,j) = ar(i,j) - fr*ortr(j) - fi*orti(j) - ai(i,j) = ai(i,j) + fr*orti(j) - fi*ortr(j) -150 end do + do 150 j = ml, igh + ar(i,j) = ar(i,j) - fr * ortr(j) - fi * orti(j) + ai(i,j) = ai(i,j) + fr * orti(j) - fi * ortr(j) +150 continue ! -160 end do +160 continue ! - ortr(ml) = scale*ortr(ml) - orti(ml) = scale*orti(ml) - ar(ml,ml - 1) = -g*ar(ml,ml - 1) - ai(ml,ml - 1) = -g*ai(ml,ml - 1) -180 end do + ortr(ml) = scale * ortr(ml) + orti(ml) = scale * orti(ml) + ar(ml,ml-1) = -g * ar(ml,ml-1) + ai(ml,ml-1) = -g * ai(ml,ml-1) +180 continue ! -200 return +200 return end subroutine corth - + subroutine comqr2(nm,nl,low,igh,ortr,orti,hr,hi,wr,wi,zr,zi,ierr) ! MESHED overflow control WITH vectors of isolated roots (10/19/89 BSG) ! MESHED overflow control WITH triangular multiply (10/30/89 BSG) @@ -460,147 +460,147 @@ subroutine comqr2(nm,nl,low,igh,ortr,orti,hr,hi,wr,wi,zr,zi,ierr) ! this version dated october 1989. ! ! ------------------------------------------------------------------ - integer i,j,k,l,ml,nl,en,ii,jj,ll,nm,nn,igh,ip1, & - itn,its,low,lp1,enm1,iend,ierr + integer i,j,k,l,ml,nl,en,ii,jj,ll,nm,nn,igh,ip1,& + itn,its,low,lp1,enm1,iend,ierr real(kind(0d0)),dimension(nm,nl) :: hr,hi,zr,zi real(kind(0d0)),dimension(nl) :: wr,wi real(kind(0d0)),dimension(igh) :: ortr,orti - real(kind(0d0)) :: si,sr,ti,tr,xi,xr,yi,yr,zzi,zzr, & - norm,tst1,tst2,c,d + real(kind(0d0)) :: si,sr,ti,tr,xi,xr,yi,yr,zzi,zzr,& + norm,tst1,tst2,c,d ! ierr = 0 ! .......... initialize eigenvector matrix .......... - do 101 j = 1,nl -! - do 100 i = 1,nl - zr(i,j) = 0.0d0 - zi(i,j) = 0.0d0 -100 end do - zr(j,j) = 1.0d0 -101 end do + do 101 j = 1, nl +! + do 100 i = 1, nl + zr(i,j) = 0.0d0 + zi(i,j) = 0.0d0 +100 continue + zr(j,j) = 1.0d0 +101 continue ! .......... form the matrix of accumulated transformations ! from the information left by corth .......... iend = igh - low - 1 - if (iend < 0) go to 180 - if (iend == 0) go to 150 - if (iend > 0) go to 105 + if (iend .lt. 0) go to 180 + if (iend .eq. 0) go to 150 + if (iend .gt. 0) go to 105 ! .......... for i=igh-1 step -1 until low+1 do -- .......... -105 do 140 ii = 1,iend - i = igh - ii - if (dabs(ortr(i)) == 0d0 .and. dabs(orti(i)) == 0d0) go to 140 - if (dabs(hr(i,i - 1)) == 0d0 .and. dabs(hi(i,i - 1)) == 0d0) go to 140 +105 do 140 ii = 1, iend + i = igh - ii + if (dabs(ortr(i)) .eq. 0d0 .and. dabs(orti(i)) .eq. 0d0) go to 140 + if (dabs(hr(i,i-1)) .eq. 0d0 .and. dabs(hi(i,i-1)) .eq. 0d0) go to 140 ! .......... norm below is negative of h formed in corth .......... - norm = hr(i,i - 1)*ortr(i) + hi(i,i - 1)*orti(i) - ip1 = i + 1 + norm = hr(i,i-1) * ortr(i) + hi(i,i-1) * orti(i) + ip1 = i + 1 - do 110 k = ip1,igh - ortr(k) = hr(k,i - 1) - orti(k) = hi(k,i - 1) -110 end do + do 110 k = ip1, igh + ortr(k) = hr(k,i-1) + orti(k) = hi(k,i-1) +110 continue ! - do 130 j = i,igh - sr = 0.0d0 - si = 0.0d0 + do 130 j = i, igh + sr = 0.0d0 + si = 0.0d0 ! - do 115 k = i,igh - sr = sr + ortr(k)*zr(k,j) + orti(k)*zi(k,j) - si = si + ortr(k)*zi(k,j) - orti(k)*zr(k,j) -115 end do + do 115 k = i, igh + sr = sr + ortr(k) * zr(k,j) + orti(k) * zi(k,j) + si = si + ortr(k) * zi(k,j) - orti(k) * zr(k,j) +115 continue ! - sr = sr/norm - si = si/norm + sr = sr / norm + si = si / norm ! - do 120 k = i,igh - zr(k,j) = zr(k,j) + sr*ortr(k) - si*orti(k) - zi(k,j) = zi(k,j) + sr*orti(k) + si*ortr(k) -120 end do + do 120 k = i, igh + zr(k,j) = zr(k,j) + sr * ortr(k) - si * orti(k) + zi(k,j) = zi(k,j) + sr * orti(k) + si * ortr(k) +120 continue ! -130 end do +130 continue ! -140 end do +140 continue ! .......... create real subdiagonal elements .......... 150 l = low + 1 ! - do 170 i = l,igh - ll = min0(i + 1,igh) - if (dabs(hi(i,i - 1)) == 0d0) go to 170 - call pythag(hr(i,i - 1),hi(i,i - 1),norm) - yr = hr(i,i - 1)/norm - yi = hi(i,i - 1)/norm - hr(i,i - 1) = norm - hi(i,i - 1) = 0.0d0 -! - do 155 j = i,nl - si = yr*hi(i,j) - yi*hr(i,j) - hr(i,j) = yr*hr(i,j) + yi*hi(i,j) - hi(i,j) = si -155 end do -! - do 160 j = 1,ll - si = yr*hi(j,i) + yi*hr(j,i) - hr(j,i) = yr*hr(j,i) - yi*hi(j,i) - hi(j,i) = si -160 end do -! - do 165 j = low,igh - si = yr*zi(j,i) + yi*zr(j,i) - zr(j,i) = yr*zr(j,i) - yi*zi(j,i) - zi(j,i) = si -165 end do - -170 end do + do 170 i = l, igh + ll = min0(i+1,igh) + if (dabs(hi(i,i-1)) .eq. 0d0) go to 170 + call pythag(hr(i,i-1),hi(i,i-1),norm) + yr = hr(i,i-1) / norm + yi = hi(i,i-1) / norm + hr(i,i-1) = norm + hi(i,i-1) = 0.0d0 +! + do 155 j = i, nl + si = yr * hi(i,j) - yi * hr(i,j) + hr(i,j) = yr * hr(i,j) + yi * hi(i,j) + hi(i,j) = si +155 continue +! + do 160 j = 1, ll + si = yr * hi(j,i) + yi * hr(j,i) + hr(j,i) = yr * hr(j,i) - yi * hi(j,i) + hi(j,i) = si +160 continue +! + do 165 j = low, igh + si = yr * zi(j,i) + yi * zr(j,i) + zr(j,i) = yr * zr(j,i) - yi * zi(j,i) + zi(j,i) = si +165 continue + +170 continue ! .......... store roots isolated by cbal .......... -180 do 200 i = 1,nl - if (i >= low .and. i <= igh) go to 200 - wr(i) = hr(i,i) - wi(i) = hi(i,i) -200 end do +180 do 200 i = 1, nl + if (i .ge. low .and. i .le. igh) go to 200 + wr(i) = hr(i,i) + wi(i) = hi(i,i) +200 continue ! en = igh tr = 0.0d0 ti = 0.0d0 itn = 30*nl ! .......... search for next eigenvalue .......... -220 if (en < low) go to 680 +220 if (en .lt. low) go to 680 its = 0 enm1 = en - 1 ! .......... look for single small sub-diagonal element ! for l=en step -1 until low do -- .......... -240 do 260 ll = low,en - l = en + low - ll - if (l == low) go to 300 - tst1 = dabs(hr(l - 1,l - 1)) + dabs(hi(l - 1,l - 1)) & - + dabs(hr(l,l)) + dabs(hi(l,l)) - tst2 = tst1 + dabs(hr(l,l - 1)) - if (tst2 == tst1) go to 300 -260 end do +240 do 260 ll = low, en + l = en + low - ll + if (l .eq. low) go to 300 + tst1 = dabs(hr(l-1,l-1)) + dabs(hi(l-1,l-1)) & + + dabs(hr(l,l)) + dabs(hi(l,l)) + tst2 = tst1 + dabs(hr(l,l-1)) + if (tst2 .eq. tst1) go to 300 +260 continue ! .......... form shift .......... -300 if (l == en) go to 660 - if (itn == 0) go to 1000 - if (its == 10 .or. its == 20) go to 320 +300 if (l .eq. en) go to 660 + if (itn .eq. 0) go to 1000 + if (its .eq. 10 .or. its .eq. 20) go to 320 sr = hr(en,en) si = hi(en,en) - xr = hr(enm1,en)*hr(en,enm1) - xi = hi(enm1,en)*hr(en,enm1) - if (xr == 0.0d0 .and. xi == 0.0d0) go to 340 - yr = (hr(enm1,enm1) - sr)/2.0d0 - yi = (hi(enm1,enm1) - si)/2.0d0 - call csroot(yr**2 - yi**2 + xr,2.0d0*yr*yi + xi,zzr,zzi) - if (yr*zzr + yi*zzi >= 0.0d0) go to 310 + xr = hr(enm1,en) * hr(en,enm1) + xi = hi(enm1,en) * hr(en,enm1) + if (xr .eq. 0.0d0 .and. xi .eq. 0.0d0) go to 340 + yr = (hr(enm1,enm1) - sr) / 2.0d0 + yi = (hi(enm1,enm1) - si) / 2.0d0 + call csroot(yr**2-yi**2+xr,2.0d0*yr*yi+xi,zzr,zzi) + if (yr * zzr + yi * zzi .ge. 0.0d0) go to 310 zzr = -zzr zzi = -zzi -310 call cdiv(xr,xi,yr + zzr,yi + zzi,xr,xi) +310 call cdiv(xr,xi,yr+zzr,yi+zzi,xr,xi) sr = sr - xr si = si - xi go to 340 ! .......... form exceptional shift .......... -320 sr = dabs(hr(en,enm1)) + dabs(hr(enm1,en - 2)) +320 sr = dabs(hr(en,enm1)) + dabs(hr(enm1,en-2)) si = 0.0d0 ! -340 do 360 i = low,en - hr(i,i) = hr(i,i) - sr - hi(i,i) = hi(i,i) - si -360 end do +340 do 360 i = low, en + hr(i,i) = hr(i,i) - sr + hi(i,i) = hi(i,i) - si +360 continue ! tr = tr + sr ti = ti + si @@ -609,94 +609,94 @@ subroutine comqr2(nm,nl,low,igh,ortr,orti,hr,hi,wr,wi,zr,zi,ierr) ! .......... reduce to triangle (rows) .......... lp1 = l + 1 ! - do 500 i = lp1,en - sr = hr(i,i - 1) - hr(i,i - 1) = 0.0d0 - call pythag(hr(i - 1,i - 1),hi(i - 1,i - 1),c) - call pythag(c,sr,norm) - xr = hr(i - 1,i - 1)/norm - wr(i - 1) = xr - xi = hi(i - 1,i - 1)/norm - wi(i - 1) = xi - hr(i - 1,i - 1) = norm - hi(i - 1,i - 1) = 0.0d0 - hi(i,i - 1) = sr/norm -! - do 490 j = i,nl - yr = hr(i - 1,j) - yi = hi(i - 1,j) - zzr = hr(i,j) - zzi = hi(i,j) - hr(i - 1,j) = xr*yr + xi*yi + hi(i,i - 1)*zzr - hi(i - 1,j) = xr*yi - xi*yr + hi(i,i - 1)*zzi - hr(i,j) = xr*zzr - xi*zzi - hi(i,i - 1)*yr - hi(i,j) = xr*zzi + xi*zzr - hi(i,i - 1)*yi -490 end do -! -500 end do + do 500 i = lp1, en + sr = hr(i,i-1) + hr(i,i-1) = 0.0d0 + call pythag(hr(i-1,i-1),hi(i-1,i-1),c) + call pythag(c,sr,norm) + xr = hr(i-1,i-1) / norm + wr(i-1) = xr + xi = hi(i-1,i-1) / norm + wi(i-1) = xi + hr(i-1,i-1) = norm + hi(i-1,i-1) = 0.0d0 + hi(i,i-1) = sr / norm +! + do 490 j = i, nl + yr = hr(i-1,j) + yi = hi(i-1,j) + zzr = hr(i,j) + zzi = hi(i,j) + hr(i-1,j) = xr * yr + xi * yi + hi(i,i-1) * zzr + hi(i-1,j) = xr * yi - xi * yr + hi(i,i-1) * zzi + hr(i,j) = xr * zzr - xi * zzi - hi(i,i-1) * yr + hi(i,j) = xr * zzi + xi * zzr - hi(i,i-1) * yi +490 continue +! +500 continue ! si = hi(en,en) - if (dabs(si) == 0d0) go to 540 + if (dabs(si) .eq. 0d0) go to 540 call pythag(hr(en,en),si,norm) - sr = hr(en,en)/norm - si = si/norm + sr = hr(en,en) / norm + si = si / norm hr(en,en) = norm hi(en,en) = 0.0d0 - if (en == nl) go to 540 + if (en .eq. nl) go to 540 ip1 = en + 1 ! - do 520 j = ip1,nl - yr = hr(en,j) - yi = hi(en,j) - hr(en,j) = sr*yr + si*yi - hi(en,j) = sr*yi - si*yr -520 end do + do 520 j = ip1, nl + yr = hr(en,j) + yi = hi(en,j) + hr(en,j) = sr * yr + si * yi + hi(en,j) = sr * yi - si * yr +520 continue ! .......... inverse operation (columns) .......... -540 do 600 j = lp1,en - xr = wr(j - 1) - xi = wi(j - 1) -! - do 580 i = 1,j - yr = hr(i,j - 1) - yi = 0.0d0 - zzr = hr(i,j) - zzi = hi(i,j) - if (i == j) go to 560 - yi = hi(i,j - 1) - hi(i,j - 1) = xr*yi + xi*yr + hi(j,j - 1)*zzi -560 hr(i,j - 1) = xr*yr - xi*yi + hi(j,j - 1)*zzr - hr(i,j) = xr*zzr + xi*zzi - hi(j,j - 1)*yr - hi(i,j) = xr*zzi - xi*zzr - hi(j,j - 1)*yi -580 end do -! - do 590 i = low,igh - yr = zr(i,j - 1) - yi = zi(i,j - 1) - zzr = zr(i,j) - zzi = zi(i,j) - zr(i,j - 1) = xr*yr - xi*yi + hi(j,j - 1)*zzr - zi(i,j - 1) = xr*yi + xi*yr + hi(j,j - 1)*zzi - zr(i,j) = xr*zzr + xi*zzi - hi(j,j - 1)*yr - zi(i,j) = xr*zzi - xi*zzr - hi(j,j - 1)*yi -590 end do - -600 end do -! - if (dabs(si) == 0d0) go to 240 -! - do 630 i = 1,en - yr = hr(i,en) - yi = hi(i,en) - hr(i,en) = sr*yr - si*yi - hi(i,en) = sr*yi + si*yr -630 end do -! - do 640 i = low,igh - yr = zr(i,en) - yi = zi(i,en) - zr(i,en) = sr*yr - si*yi - zi(i,en) = sr*yi + si*yr -640 end do +540 do 600 j = lp1, en + xr = wr(j-1) + xi = wi(j-1) +! + do 580 i = 1, j + yr = hr(i,j-1) + yi = 0.0d0 + zzr = hr(i,j) + zzi = hi(i,j) + if (i .eq. j) go to 560 + yi = hi(i,j-1) + hi(i,j-1) = xr * yi + xi * yr + hi(j,j-1) * zzi +560 hr(i,j-1) = xr * yr - xi * yi + hi(j,j-1) * zzr + hr(i,j) = xr * zzr + xi * zzi - hi(j,j-1) * yr + hi(i,j) = xr * zzi - xi * zzr - hi(j,j-1) * yi +580 continue +! + do 590 i = low, igh + yr = zr(i,j-1) + yi = zi(i,j-1) + zzr = zr(i,j) + zzi = zi(i,j) + zr(i,j-1) = xr * yr - xi * yi + hi(j,j-1) * zzr + zi(i,j-1) = xr * yi + xi * yr + hi(j,j-1) * zzi + zr(i,j) = xr * zzr + xi * zzi - hi(j,j-1) * yr + zi(i,j) = xr * zzi - xi * zzr - hi(j,j-1) * yi +590 continue + +600 continue +! + if (dabs(si) .eq. 0d0) go to 240 +! + do 630 i = 1, en + yr = hr(i,en) + yi = hi(i,en) + hr(i,en) = sr * yr - si * yi + hi(i,en) = sr * yi + si * yr +630 continue +! + do 640 i = low, igh + yr = zr(i,en) + yi = zi(i,en) + zr(i,en) = sr * yr - si * yi + zi(i,en) = sr * yi + si * yr +640 continue ! go to 240 ! .......... a root found .......... @@ -710,88 +710,88 @@ subroutine comqr2(nm,nl,low,igh,ortr,orti,hr,hi,wr,wi,zr,zi,ierr) ! vectors of upper triangular form .......... 680 norm = 0.0d0 ! - do i = 1,nl - do j = i,nl - tr = dabs(hr(i,j)) + dabs(hi(i,j)) - if (tr > norm) norm = tr - end do + do i = 1, nl + do j = i, nl + tr = dabs(hr(i,j)) + dabs(hi(i,j)) + if (tr .gt. norm) norm = tr + end do end do ! - if (nl == 1 .or. norm == 0d0) go to 1001 + if (nl .eq. 1 .or. norm .eq. 0d0) go to 1001 ! .......... for en=nl step -1 until 2 do -- .......... - do 800 nn = 2,nl - en = nl + 2 - nn - xr = wr(en) - xi = wi(en) - hr(en,en) = 1.0d0 - hi(en,en) = 0.0d0 - enm1 = en - 1 + do 800 nn = 2, nl + en = nl + 2 - nn + xr = wr(en) + xi = wi(en) + hr(en,en) = 1.0d0 + hi(en,en) = 0.0d0 + enm1 = en - 1 ! .......... for i=en-1 step -1 until 1 do -- .......... - do 780 ii = 1,enm1 - i = en - ii - zzr = 0.0d0 - zzi = 0.0d0 - ip1 = i + 1 - - do 740 j = ip1,en - zzr = zzr + hr(i,j)*hr(j,en) - hi(i,j)*hi(j,en) - zzi = zzi + hr(i,j)*hi(j,en) + hi(i,j)*hr(j,en) -740 end do -! - yr = xr - wr(i) - yi = xi - wi(i) - if (yr /= 0.0d0 .or. yi /= 0.0d0) go to 765 - tst1 = norm - yr = tst1 -760 yr = 0.01d0*yr - tst2 = norm + yr - if (tst2 > tst1) go to 760 -765 continue - call cdiv(zzr,zzi,yr,yi,hr(i,en),hi(i,en)) + do 780 ii = 1, enm1 + i = en - ii + zzr = 0.0d0 + zzi = 0.0d0 + ip1 = i + 1 + + do 740 j = ip1, en + zzr = zzr + hr(i,j) * hr(j,en) - hi(i,j) * hi(j,en) + zzi = zzi + hr(i,j) * hi(j,en) + hi(i,j) * hr(j,en) +740 continue +! + yr = xr - wr(i) + yi = xi - wi(i) + if (yr .ne. 0.0d0 .or. yi .ne. 0.0d0) go to 765 + tst1 = norm + yr = tst1 +760 yr = 0.01d0 * yr + tst2 = norm + yr + if (tst2 .gt. tst1) go to 760 +765 continue + call cdiv(zzr,zzi,yr,yi,hr(i,en),hi(i,en)) ! .......... overflow control .......... - tr = dabs(hr(i,en)) + dabs(hi(i,en)) - if (tr == 0.0d0) go to 780 - tst1 = tr - tst2 = tst1 + 1.0d0/tst1 - if (tst2 > tst1) go to 780 - do 770 j = i,en - hr(j,en) = hr(j,en)/tr - hi(j,en) = hi(j,en)/tr -770 end do -! -780 end do -! -800 end do + tr = dabs(hr(i,en)) + dabs(hi(i,en)) + if (tr .eq. 0.0d0) go to 780 + tst1 = tr + tst2 = tst1 + 1.0d0/tst1 + if (tst2 .gt. tst1) go to 780 + do 770 j = i, en + hr(j,en) = hr(j,en)/tr + hi(j,en) = hi(j,en)/tr +770 continue +! +780 continue +! +800 continue ! .......... end backsubstitution .......... ! .......... vectors of isolated roots .......... - do 840 i = 1,nl - if (i >= low .and. i <= igh) go to 840 + do 840 i = 1, nl + if (i .ge. low .and. i .le. igh) go to 840 ! - do 820 j = I,nl - zr(i,j) = hr(i,j) - zi(i,j) = hi(i,j) -820 end do + do 820 j = I, nl + zr(i,j) = hr(i,j) + zi(i,j) = hi(i,j) +820 continue ! -840 end do +840 continue ! .......... multiply by transformation matrix to give ! vectors of original full matrix. ! for j=nl step -1 until low do -- .......... - do jj = low,nl - j = nl + low - jj - ml = min0(j,igh) -! - do i = low,igh - zzr = 0.0d0 - zzi = 0.0d0 -! - do 860 k = low,ml - zzr = zzr + zr(i,k)*hr(k,j) - zi(i,k)*hi(k,j) - zzi = zzi + zr(i,k)*hi(k,j) + zi(i,k)*hr(k,j) -860 end do -! - zr(i,j) = zzr - zi(i,j) = zzi - end do + do jj = low, nl + j = nl + low - jj + ml = min0(j,igh) +! + do i = low, igh + zzr = 0.0d0 + zzi = 0.0d0 +! + do 860 k = low, ml + zzr = zzr + zr(i,k) * hr(k,j) - zi(i,k) * hi(k,j) + zzi = zzi + zr(i,k) * hi(k,j) + zi(i,k) * hr(k,j) +860 continue +! + zr(i,j) = zzr + zi(i,j) = zzi + end do end do ! go to 1001 @@ -843,51 +843,51 @@ subroutine cbabk2(nm,nl,low,igh,scale,ml,zr,zi) ! ! ------------------------------------------------------------------ ! - integer i,j,k,ml,nl,ii,nm,igh,low - double precision scale(nl),zr(nm,ml),zi(nm,ml) - double precision s + integer i,j,k,ml,nl,ii,nm,igh,low + double precision scale(nl),zr(nm,ml),zi(nm,ml) + double precision s - if (ml == 0) go to 200 - if (igh == low) go to 120 + if (ml .eq. 0) go to 200 + if (igh .eq. low) go to 120 ! - do 110 i = low,igh - s = scale(i) + do 110 i = low, igh + s = scale(i) ! .......... left hand eigenvectors are back transformed ! if the foregoing statement is replaced by ! s=1.0d0/scale(i). .......... - do 100 j = 1,ml - zr(i,j) = zr(i,j)*s - zi(i,j) = zi(i,j)*s -100 end do + do 100 j = 1, ml + zr(i,j) = zr(i,j) * s + zi(i,j) = zi(i,j) * s +100 continue ! -110 end do +110 continue ! .......... for i=low-1 step -1 until 1, ! igh+1 step 1 until nl do -- .......... -120 do 140 ii = 1,nl - i = ii - if (i >= low .and. i <= igh) go to 140 - if (i < low) i = low - ii - k = scale(i) - if (k == i) go to 140 -! - do 130 j = 1,ml - s = zr(i,j) - zr(i,j) = zr(k,j) - zr(k,j) = s - s = zi(i,j) - zi(i,j) = zi(k,j) - zi(k,j) = s -130 end do -! -140 end do -! -200 return +120 do 140 ii = 1, nl + i = ii + if (i .ge. low .and. i .le. igh) go to 140 + if (i .lt. low) i = low - ii + k = scale(i) + if (k .eq. i) go to 140 +! + do 130 j = 1, ml + s = zr(i,j) + zr(i,j) = zr(k,j) + zr(k,j) = s + s = zi(i,j) + zi(i,j) = zi(k,j) + zi(k,j) = s +130 continue +! +140 continue +! +200 return end subroutine cbabk2 - + subroutine csroot(xr,xi,yr,yi) real(kind(0d0)) :: xr,xi,yr,yi ! -! (yr,yi) = complex dsqrt(xr,xi) +! (yr,yi) = complex dsqrt(xr,xi) ! branch chosen so that yr .ge. 0.0 and sign(yi) .eq. sign(xi) ! real(kind(0d0)) :: s,tr,ti,c @@ -895,14 +895,14 @@ subroutine csroot(xr,xi,yr,yi) ti = xi call pythag(tr,ti,c) s = dsqrt(0.5d0*(c + dabs(tr))) - if (tr >= 0.0d0) yr = s - if (ti < 0.0d0) s = -s - if (tr <= 0.0d0) yi = s - if (tr < 0.0d0) yr = 0.5d0*(ti/yi) - if (tr > 0.0d0) yi = 0.5d0*(ti/yr) + if (tr .ge. 0.0d0) yr = s + if (ti .lt. 0.0d0) s = -s + if (tr .le. 0.0d0) yi = s + if (tr .lt. 0.0d0) yr = 0.5d0*(ti/yi) + if (tr .gt. 0.0d0) yi = 0.5d0*(ti/yr) return end subroutine csroot - + subroutine cdiv(ar,ai,br,bi,cr,ci) real(kind(0d0)) :: ar,ai,br,bi,cr,ci ! @@ -927,18 +927,18 @@ subroutine pythag(a,b,c) ! real(kind(0d0)) :: p,r,s,t,u p = dmax1(dabs(a),dabs(b)) - if (p == 0.0d0) go to 20 + if (p .eq. 0.0d0) go to 20 r = (dmin1(dabs(a),dabs(b))/p)**2 10 continue t = 4.0d0 + r - if (t == 4.0d0) go to 20 + if (t .eq. 4.0d0) go to 20 s = r/t u = 1.0d0 + 2.0d0*s p = u*p - r = (s/u)**2*r + r = (s/u)**2 * r go to 10 20 c = p return end subroutine pythag -end module m_eigen_solver +end module m_eigen_solver \ No newline at end of file diff --git a/src/common/m_phase_change.fpp b/src/common/m_phase_change.fpp new file mode 100644 index 000000000..c805ce82d --- /dev/null +++ b/src/common/m_phase_change.fpp @@ -0,0 +1,795 @@ +!> +!! @file m_phase_change.fpp +!! @brief Contains module m_phasechange + +#:include 'macros.fpp' + +!> @brief This module is used to relax the model equations (6-eqn model) +!> towards pressure and temperature (6-eqn to 4-eqn), and (if wanted) Gibbs free +!> energies (6-eqn to 4-eqn) equilibrium through an infinitely fast (algebraic) +!> procedure. +module m_phase_change + +#ifndef MFC_POST_PROCESS + + ! Dependencies ============================================================= + + use m_derived_types !< Definitions of the derived types + + use m_global_parameters !< Definitions of the global parameters + + use m_mpi_proxy !< Message passing interface (MPI) module proxy + + use m_variables_conversion !< State variables type conversion procedures + + use ieee_arithmetic + + ! ========================================================================== + + implicit none + + private; public :: s_initialize_phasechange_module, & + s_relaxation_solver, & + s_infinite_relaxation_k, & + s_finalize_relaxation_solver_module + + !> @name Abstract interface for creating function pointers + !> @{ + abstract interface + + !> @name Abstract subroutine for the infinite relaxation solver + !> @{ + subroutine s_abstract_relaxation_solver(q_cons_vf) ! ------- + import :: scalar_field, sys_size + type(scalar_field), dimension(sys_size), intent(INOUT) :: q_cons_vf + end subroutine + !> @} + + end interface + !> @} + + !> @name Parameters for the first order transition phase change + !> @{ + integer, parameter :: max_iter = 1e8 !< max # of iterations + real(kind(0d0)), parameter :: pCr = 4.94d7 !< Critical water pressure + real(kind(0d0)), parameter :: TCr = 385.05+273.15 !< Critical water temperature + real(kind(0d0)), parameter :: mixM = 1.0d-8 !< threshold for 'mixture cell'. If Y < mixM, phase change does not happen + integer, parameter :: lp = 1 !< index for the liquid phase of the reacting fluid + integer, parameter :: vp = 2 !< index for the vapor phase of the reacting fluid + !> @} + + !> @name Gibbs free energy phase change parameters + !> @{ + real(kind(0d0)) :: A, B, C, D + !> @} + + !$acc declare create(max_iter,pCr,TCr,mixM,lp,vp,A,B,C,D) + + procedure(s_abstract_relaxation_solver), pointer :: s_relaxation_solver => null() + +contains + + !> The purpose of this subroutine is to initialize the phase change module + !! by setting the parameters needed for phase change and + !! selecting the phase change module that will be used + !! (pT- or pTg-equilibrium) + subroutine s_initialize_phasechange_module() + + ! variables used in the calculation of the saturation curves for fluids 1 and 2 + A = (gs_min(lp)*cvs(lp) - gs_min(vp)*cvs(vp) & + + qvps(vp) - qvps(lp))/((gs_min(vp) - 1.0d0)*cvs(vp)) + + B = (qvs(lp) - qvs(vp))/((gs_min(vp) - 1.0d0)*cvs(vp)) + + C = (gs_min(vp)*cvs(vp) - gs_min(lp)*cvs(lp)) & + /((gs_min(vp) - 1.0d0)*cvs(vp)) + + D = ((gs_min(lp) - 1.0d0)*cvs(lp)) & + /((gs_min(vp) - 1.0d0)*cvs(vp)) + + ! Associating procedural pointer to the subroutine that will be + ! utilized to calculate the solution to the selected relaxation system + if ((relax_model == 5) .or. (relax_model == 6)) then + s_relaxation_solver => s_infinite_relaxation_k + else + call s_mpi_abort('relaxation solver was not set!') + end if + + end subroutine s_initialize_phasechange_module !------------------------------- + + !> This subroutine is created to activate either the pT- (N fluids) or the + !! pTg-equilibrium (2 fluids for g-equilibrium) + !! model, also considering mass depletion, depending on the incoming + !! state conditions. + !! @param q_cons_vf Cell-average conservative variables + subroutine s_infinite_relaxation_k(q_cons_vf) ! ---------------- + type(scalar_field), dimension(sys_size), intent(INOUT) :: q_cons_vf + real(kind(0.0d0)) :: pS, pSOV, pSSL !< equilibrium pressure for mixture, overheated vapor, and subcooled liquid + real(kind(0.0d0)) :: TS, TSOV, TSSL, TSatOV, TSatSL !< equilibrium temperature for mixture, overheated vapor, and subcooled liquid. Saturation Temperatures at overheated vapor and subcooled liquid + real(kind(0.0d0)) :: rhoe, dynE, rhos !< total internal energy, kinetic energy, and total entropy + real(kind(0.0d0)) :: rho, rM, m1, m2, MCT !< total density, total reacting mass, individual reacting masses + real(kind(0.0d0)) :: TvF !< total volume fraction + + !$acc declare create(pS, pSOV, pSSL, TS, TSOV, TSatOV, TSatSL, TSSL, rhoe, dynE, rhos, rho, rM, m1, m2, MCT, TvF) + + real(kind(0d0)), dimension(num_fluids) :: p_infOV, p_infpT, p_infSL, sk, hk, gk, ek, rhok + + !< Generic loop iterators + integer :: i, j, k, l + + !$acc declare create(p_infOV, p_infpT, p_infSL, sk, hk, gk, ek, rhok) + + ! starting equilibrium solver + !$acc parallel loop collapse(3) gang vector default(present) private(p_infOV, p_infpT, p_infSL, sk, hk, gk, ek, rhok,pS, pSOV, pSSL, TS, TSOV, TSatOV, TSatSL, TSSL, rhoe, dynE, rhos, rho, rM, m1, m2, MCT, TvF) + do j = 0, m + do k = 0, n + do l = 0, p + + rho = 0.0d0; TvF = 0.0d0 +!$acc loop seq + do i = 1, num_fluids + + ! Mixture density + rho = rho + q_cons_vf(i + contxb - 1)%sf(j, k, l) + + ! Total Volume Fraction + TvF = TvF + q_cons_vf(i + advxb - 1)%sf(j, k, l) + + end do + + ! calculating the total reacting mass for the phase change process. By hypothesis, this should not change + ! throughout the phase-change process. + rM = q_cons_vf(lp + contxb - 1)%sf(j, k, l) + q_cons_vf(vp + contxb - 1)%sf(j, k, l) + + ! correcting negative (recating) mass fraction values in case they happen + call s_correct_partial_densities(MCT, q_cons_vf, rM, j, k, l) + + ! fixing m1 and m2 AFTER correcting the partial densities. Note that these values must be stored for the phase + ! change process that will happen a posteriori + m1 = q_cons_vf(lp + contxb - 1)%sf(j, k, l) + + m2 = q_cons_vf(vp + contxb - 1)%sf(j, k, l) + + ! kinetic energy as an auxiliary variable to the calculation of the total internal energy + dynE = 0.0d0 +!$acc loop seq + do i = momxb, momxe + + dynE = dynE + 5.0d-1*q_cons_vf(i)%sf(j, k, l)**2/rho + + end do + + ! calculating the total energy that MUST be preserved throughout the pT- and pTg-relaxation procedures + ! at each of the cells. The internal energy is calculated as the total energy minus the kinetic + ! energy to preserved its value at sharp interfaces + rhoe = q_cons_vf(E_idx)%sf(j, k, l) - dynE + + ! Calling pT-equilibrium for either finishing phase-change module, or as an IC for the pTg-equilibrium + ! for this case, MFL cannot be either 0 or 1, so I chose it to be 2 + call s_infinite_pt_relaxation_k(j, k, l, 2, pS, p_infpT, rM, q_cons_vf, rhoe, TS) + + ! check if pTg-equilibrium is required + ! NOTE that NOTHING else needs to be updated OTHER than the individual partial densities + ! given the outputs from the pT- and pTg-equilibrium solvers are just p and one of the partial masses + ! (pTg- case) + if ((relax_model == 6) .and. ((q_cons_vf(lp + contxb - 1)%sf(j, k, l) > mixM*rM) & + .and. (q_cons_vf(vp + contxb - 1)%sf(j, k, l) > mixM*rM)) & + .and. (pS < pCr) .and. (TS < TCr)) then + + ! Checking if phase change is needed, by checking whether the final solution is either subcoooled + ! liquid or overheated vapor. + + ! overheated vapor case + ! depleting the mass of liquid + q_cons_vf(lp + contxb - 1)%sf(j, k, l) = mixM*rM + + ! tranferring the total mass to vapor + q_cons_vf(vp + contxb - 1)%sf(j, k, l) = (1.0d0 - mixM)*rM + + ! calling pT-equilibrium for overheated vapor, which is MFL = 0 + call s_infinite_pt_relaxation_k(j, k, l, 0, pSOV, p_infOV, rM, q_cons_vf, rhoe, TSOV) + + ! calculating Saturation temperature + call s_TSat(pSOV, TSatOV, TSOV) + + ! subcooled liquid case + ! tranferring the total mass to liquid + q_cons_vf(lp + contxb - 1)%sf(j, k, l) = (1.0d0 - mixM)*rM + + ! depleting the mass of vapor + q_cons_vf(vp + contxb - 1)%sf(j, k, l) = mixM*rM + + ! calling pT-equilibrium for subcooled liquid, which is MFL = 1 + call s_infinite_pt_relaxation_k(j, k, l, 1, pSSL, p_infSL, rM, q_cons_vf, rhoe, TSSL) + + ! calculating Saturation temperature + call s_TSat(pSSL, TSatSL, TSSL) + + ! checking the conditions for overheated vapor and subcooled liquide + if (TSOV > TSatOV) then + + ! Assigning pressure + pS = pSOV + + ! Assigning Temperature + TS = TSOV + + ! correcting the liquid partial density + q_cons_vf(lp + contxb - 1)%sf(j, k, l) = mixM*rM + + ! correcting the vapor partial density + q_cons_vf(vp + contxb - 1)%sf(j, k, l) = (1.0d0 - mixM)*rM + + elseif (TSSL < TSatSL) then + + ! Assigning pressure + pS = pSSL + + ! Assigning Temperature + TS = TSSL + + ! correcting the liquid partial density + q_cons_vf(lp + contxb - 1)%sf(j, k, l) = (1.0d0 - mixM)*rM + + ! correcting the vapor partial density + q_cons_vf(vp + contxb - 1)%sf(j, k, l) = mixM*rM + + else + + ! returning partial pressures to what they were from the homogeneous solver + ! liquid + q_cons_vf(lp + contxb - 1)%sf(j, k, l) = m1 + + ! vapor + q_cons_vf(vp + contxb - 1)%sf(j, k, l) = m2 + + ! calling the pTg-equilibrium solver + call s_infinite_ptg_relaxation_k(j, k, l, pS, p_infpT, rhoe, q_cons_vf, TS) + + end if + + end if + + ! Calculations AFTER equilibrium + + ! entropy + sk(1:num_fluids) = cvs(1:num_fluids)*DLOG((TS**gs_min(1:num_fluids)) & + /((pS + ps_inf(1:num_fluids))**(gs_min(1:num_fluids) - 1.0d0))) + qvps(1:num_fluids) + + ! enthalpy + hk(1:num_fluids) = gs_min(1:num_fluids)*cvs(1:num_fluids)*TS & + + qvs(1:num_fluids) + + ! Gibbs-free energy + gk(1:num_fluids) = hk(1:num_fluids) - TS*sk(1:num_fluids) + + ! densities + rhok(1:num_fluids) = (pS + ps_inf(1:num_fluids)) & + /((gs_min(1:num_fluids) - 1)*cvs(1:num_fluids)*TS) + + ! internal energy + ek(1:num_fluids) = (pS + gs_min(1:num_fluids) & + *ps_inf(1:num_fluids))/(pS + ps_inf(1:num_fluids)) & + *cvs(1:num_fluids)*TS + qvs(1:num_fluids) + + ! calculating volume fractions, internal energies, and total entropy + rhos = 0.0d0 +!$acc loop seq + do i = 1, num_fluids + + ! volume fractions + q_cons_vf(i + advxb - 1)%sf(j, k, l) = q_cons_vf(i + contxb - 1)%sf(j, k, l)/rhok(i) + + ! alpha*rho*e + q_cons_vf(i + intxb - 1)%sf(j, k, l) = q_cons_vf(i + contxb - 1)%sf(j, k, l)*ek(i) + + ! Total entropy + rhos = rhos + q_cons_vf(i + contxb - 1)%sf(j, k, l)*sk(i) + + end do + end do + end do + end do + + end subroutine s_infinite_relaxation_k ! ---------------- + + !> This auxiliary subroutine is created to activate the pT-equilibrium for N fluids + !! @param j generic loop iterator for x direction + !! @param k generic loop iterator for y direction + !! @param l generic loop iterator for z direction + !! @param MFL flag that tells whether the fluid is pure gas (0), pure liquid (1), or a mixture (2) + !! @param pS equilibrium pressure at the interface + !! @param p_infpT stiffness for the participating fluids under pT-equilibrium + !! @param rM sum of the reacting masses + !! @param q_cons_vf Cell-average conservative variables + !! @param rhoe mixture energy + !! @param TS equilibrium temperature at the interface + subroutine s_infinite_pt_relaxation_k(j, k, l, MFL, pS, p_infpT, rM, q_cons_vf, rhoe, TS) +!$acc routine seq + + ! initializing variables + type(scalar_field), dimension(sys_size), intent(IN) :: q_cons_vf + real(kind(0.0d0)), intent(OUT) :: pS, TS + real(kind(0.0d0)), dimension(num_fluids), intent(OUT) :: p_infpT + real(kind(0.0d0)), intent(IN) :: rM, rhoe + integer, intent(IN) :: j, k, l, MFL + real(kind(0.0d0)), dimension(num_fluids) :: pk !< individual initial pressures + integer, dimension(num_fluids) :: ig !< flags to toggle the inclusion of fluids for the pT-equilibrium + + real(kind(0.0d0)) :: gp, gpp, hp, pO, mCP, mQ !< variables for the Newton Solver + + integer :: i, ns !< generic loop iterators + + ! auxiliary variables for the pT-equilibrium solver + mCP = 0.0d0; mQ = 0.0d0; p_infpT = ps_inf; pk(1:num_fluids) = 0.0d0 + + ig(1:num_fluids) = 0 + + ! Performing tests before initializing the pT-equilibrium +!$acc loop seq + do i = 1, num_fluids + + ! sum of the total alpha*rho*cp of the system + mCP = mCP + q_cons_vf(i + contxb - 1)%sf(j, k, l)*cvs(i)*gs_min(i) + + ! sum of the total alpha*rho*q of the system + mQ = mQ + q_cons_vf(i + contxb - 1)%sf(j, k, l)*qvs(i) + + end do + + ! Checking energy constraint + if ((rhoe - mQ - minval(p_infpT)) < 0.0d0) then + + if ((MFL == 0) .or. (MFL == 1)) then + + ! Assigning zero values for mass depletion cases + ! pressure + pS = 0.0d0 + + ! temperature + TS = 0.0d0 + + return + end if + + end if + + ! calculating initial estimate for pressure in the pT-relaxation procedure. I will also use this variable to + ! iterate over the Newton's solver + pO = 0.0d0 + + ! Maybe improve this condition afterwards. As long as the initial guess is in between -min(ps_inf) + ! and infinity, a solution should be able to be found. + pS = 1.0d4 + + ! Newton Solver for the pT-equilibrium + ns = 0 + ! change this relative error metric. 1E4 is just arbitrary + do while ((DABS(pS - pO) > palpha_eps) .and. (DABS((pS - pO)/pO) > palpha_eps/1e4) .or. (ns == 0)) + + ! increasing counter + ns = ns + 1 + + ! updating old pressure + pO = pS + + ! updating functions used in the Newton's solver + gpp = 0.0d0; gp = 0.0d0; hp = 0.0d0 +!$acc loop seq + do i = 1, num_fluids + + gp = gp + (gs_min(i) - 1.0d0)*q_cons_vf(i + contxb - 1)%sf(j, k, l)*cvs(i) & + *(rhoe + pS - mQ)/(mCP*(pS + p_infpT(i))) + + gpp = gpp + (gs_min(i) - 1.0d0)*q_cons_vf(i + contxb - 1)%sf(j, k, l)*cvs(i) & + *(p_infpT(i) - rhoe + mQ)/(mCP*(pS + p_infpT(i))**2) + + end do + + hp = 1.0d0/(rhoe + pS - mQ) + 1.0d0/(pS + minval(p_infpT)) + + ! updating common pressure for the newton solver + pS = pO + ((1.0d0 - gp)/gpp)/(1.0d0 - (1.0d0 - gp + DABS(1.0d0 - gp)) & + /(2.0d0*gpp)*hp) + end do + + ! common temperature + TS = (rhoe + pS - mQ)/mCP + + end subroutine s_infinite_pt_relaxation_k ! ----------------------- + + !> This auxiliary subroutine is created to activate the pTg-equilibrium for N fluids under pT + !! and 2 fluids under pTg-equilibrium. There is a final common p and T during relaxation + !! @param j generic loop iterator for x direction + !! @param k generic loop iterator for y direction + !! @param l generic loop iterator for z direction + !! @param pS equilibrium pressure at the interface + !! @param p_infpT stiffness for the participating fluids under pT-equilibrium + !! @param rhoe mixture energy + !! @param q_cons_vf Cell-average conservative variables + !! @param TS equilibrium temperature at the interface + subroutine s_infinite_ptg_relaxation_k(j, k, l, pS, p_infpT, rhoe, q_cons_vf, TS) +!$acc routine seq + + type(scalar_field), dimension(sys_size), intent(INOUT) :: q_cons_vf + real(kind(0.0d0)), dimension(num_fluids), intent(IN) :: p_infpT + real(kind(0.0d0)), intent(INOUT) :: pS, TS + real(kind(0.0d0)), intent(IN) :: rhoe + integer, intent(IN) :: j, k, l + real(kind(0.0d0)), dimension(num_fluids) :: p_infpTg !< stiffness for the participating fluids for pTg-equilibrium + real(kind(0.0d0)), dimension(2, 2) :: Jac, InvJac, TJac !< matrices for the Newton Solver + real(kind(0.0d0)), dimension(2) :: R2D, DeltamP !< residual and correction array + real(kind(0.0d0)) :: Om ! underrelaxation factor + real(kind(0.0d0)) :: mCP, mCPD, mCVGP, mCVGP2, mQ, mQD ! auxiliary variables for the pTg-solver + + !< Generic loop iterators + integer :: i, ns + + ! pTg-equilibrium solution procedure + ! Newton Solver parameters + ! counter + ns = 0 + + ! Relaxation factor + Om = 1.0d-3 + + p_infpTg = p_infpT + + if (((pS < 0.0d0) .and. ((q_cons_vf(lp + contxb - 1)%sf(j, k, l) & + + q_cons_vf(vp + contxb - 1)%sf(j, k, l)) > ((rhoe & + - gs_min(lp)*ps_inf(lp)/(gs_min(lp) - 1))/qvs(lp)))) .or. & + ((pS >= 0.0d0) .and. (pS < 1.0d-1))) then + + ! improve this initial condition + pS = 1.0d4 + + end if + + ! Loop until the solution for F(X) is satisfied + ! Check whether I need to use both absolute and relative values + ! for the residual, and how to do it adequately. + ! Dummy guess to start the pTg-equilibrium problem. + ! improve this initial condition + R2D(1) = 0.0d0; R2D(2) = 0.0d0 + DeltamP(1) = 0.0d0; DeltamP(2) = 0.0d0 + do while (((DSQRT(R2D(1)**2 + R2D(2)**2) > ptgalpha_eps) & + .and. ((DSQRT(R2D(1)**2 + R2D(2)**2)/rhoe) > (ptgalpha_eps/1d6))) & + .or. (ns == 0)) + + ! Updating counter for the iterative procedure + ns = ns + 1 + + ! Auxiliary variables to help in the calculation of the residue + mCP = 0.0d0; mCPD = 0.0d0; mCVGP = 0.0d0; mCVGP2 = 0.0d0; mQ = 0.0d0; mQD = 0.0d0 + ! Those must be updated through the iterations, as they either depend on + ! the partial masses for all fluids, or on the equilibrium pressure +!$acc loop seq + do i = 1, num_fluids + + ! sum of the total alpha*rho*cp of the system + mCP = mCP + q_cons_vf(i + contxb - 1)%sf(j, k, l) & + *cvs(i)*gs_min(i) + + ! sum of the total alpha*rho*q of the system + mQ = mQ + q_cons_vf(i + contxb - 1)%sf(j, k, l)*qvs(i) + + ! These auxiliary variables now need to be updated, as the partial densities now + ! vary at every iteration + if ((i /= lp) .and. (i /= vp)) then + + mCVGP = mCVGP + q_cons_vf(i + contxb - 1)%sf(j, k, l) & + *cvs(i)*(gs_min(i) - 1)/(pS + ps_inf(i)) + + mCVGP2 = mCVGP2 + q_cons_vf(i + contxb - 1)%sf(j, k, l) & + *cvs(i)*(gs_min(i) - 1)/((pS + ps_inf(i))**2) + + mQD = mQD + q_cons_vf(i + contxb - 1)%sf(j, k, l)*qvs(i) + + ! sum of the total alpha*rho*cp of the system + mCPD = mCPD + q_cons_vf(i + contxb - 1)%sf(j, k, l)*cvs(i) & + *gs_min(i) + + end if + + end do + + ! calculating the (2D) Jacobian Matrix used in the solution of the pTg-quilibrium model + call s_compute_jacobian_matrix(InvJac, j, Jac, k, l, mCPD, mCVGP, mCVGP2, pS, q_cons_vf, TJac) + + ! calculating correction array for Newton's method + DeltamP = -1.0d0*matmul(InvJac, R2D) + + ! updating two reacting 'masses'. Recall that inert 'masses' do not change during the phase change + ! liquid + q_cons_vf(lp + contxb - 1)%sf(j, k, l) = q_cons_vf(lp + contxb - 1)%sf(j, k, l) + Om*DeltamP(1) + + ! gas + q_cons_vf(vp + contxb - 1)%sf(j, k, l) = q_cons_vf(vp + contxb - 1)%sf(j, k, l) - Om*DeltamP(1) + + ! updating pressure + pS = pS + Om*DeltamP(2) + + ! calculating residuals, which are (i) the difference between the Gibbs Free energy of the gas and the liquid + ! and (ii) the energy before and after the phase-change process. + call s_compute_pTg_residue(j, k, l, mCPD, mCVGP, mQD, q_cons_vf, pS, rhoe, R2D) + + end do + + ! common temperature + TS = (rhoe + pS - mQ)/mCP + + end subroutine s_infinite_ptg_relaxation_k ! ----------------------- + + !> This auxiliary subroutine corrects the partial densities of the REACTING fluids in case one of them is negative + !! but their sum is positive. Inert phases are not corrected at this moment + !! @param MCT partial density correction parameter + !! @param q_cons_vf Cell-average conservative variables + !! @param rM sum of the reacting masses + !! @param j generic loop iterator for x direction + !! @param k generic loop iterator for y direction + !! @param l generic loop iterator for z direction + subroutine s_correct_partial_densities(MCT, q_cons_vf, rM, j, k, l) +!$acc routine seq + + !> @name variables for the correction of the reacting partial densities + !> @{ + type(scalar_field), dimension(sys_size), intent(INOUT) :: q_cons_vf + real(kind(0.0d0)), intent(INOUT) :: rM + real(kind(0.0d0)), intent(OUT) :: MCT + integer, intent(IN) :: j, k, l + !> @} + + if (rM < 0.0d0) then + + if ( (q_cons_vf(lp + contxb - 1)%sf(j, k, l) .ge. -1.0d0*mixM) .and. & + (q_cons_vf(vp + contxb - 1)%sf(j, k, l) .ge. -1.0d0*mixM) ) then + + q_cons_vf(lp + contxb - 1)%sf(j, k, l) = 0.0d0 + + q_cons_vf(vp + contxb - 1)%sf(j, k, l) = 0.0d0 + + rM = q_cons_vf(lp + contxb - 1)%sf(j, k, l) + q_cons_vf(vp + contxb - 1)%sf(j, k, l) + + end if + + end if + + ! Defining the correction in terms of an absolute value might not be the best practice. + ! Maybe a good way to do this is to partition the partial densities, giving a small percentage of the total reacting density + MCT = 2*mixM + + ! correcting the partial densities of the reacting fluids. What to do for the nonreacting ones? + if (q_cons_vf(lp + contxb - 1)%sf(j, k, l) < 0.0d0) then + + q_cons_vf(lp + contxb - 1)%sf(j, k, l) = MCT*rM + + q_cons_vf(vp + contxb - 1)%sf(j, k, l) = (1.0d0 - MCT)*rM + + elseif (q_cons_vf(vp + contxb - 1)%sf(j, k, l) < 0.0d0) then + + q_cons_vf(lp + contxb - 1)%sf(j, k, l) = (1.0d0 - MCT)*rM + + q_cons_vf(vp + contxb - 1)%sf(j, k, l) = MCT*rM + + end if + + end subroutine s_correct_partial_densities + + !> This auxiliary subroutine calculates the 2 x 2 Jacobian and, its inverse and transpose + !! to be used in the pTg-equilibirium procedure + !! @param InvJac Inverse of the Jacobian Matrix + !! @param j generic loop iterator for x direction + !! @param Jac Jacobian Matrix + !! @param k generic loop iterator for y direction + !! @param l generic loop iterator for z direction + !! @param mCPD sum of the total alpha*rho*cp + !! @param mCVGP auxiliary variable for the calculation of the matrices: alpha*rho*cv*(g-1)/press + !! @param mCVGP2 auxiliary variable for the calculation of the matrices: alpha*rho*cv*(g-1)/press^2 + !! @param pS equilibrium pressure at the interface + !! @param q_cons_vf Cell-average conservative variables + !! @param TJac Transpose of the Jacobian Matrix + subroutine s_compute_jacobian_matrix(InvJac, j, Jac, k, l, mCPD, mCVGP, mCVGP2, pS, q_cons_vf, TJac) +!$acc routine seq + + type(scalar_field), dimension(sys_size), intent(IN) :: q_cons_vf + real(kind(0.0d0)), intent(IN) :: pS, mCPD, mCVGP, mCVGP2 + integer, intent(IN) :: j, k, l + real(kind(0.0d0)), dimension(2, 2), intent(OUT) :: Jac, InvJac, TJac + real(kind(0.0d0)) :: ml, mT, TS, dFdT, dTdm, dTdp ! mass of the reacting fluid, total reacting mass, and auxiliary variables + + ! mass of the reacting liquid + ml = q_cons_vf(lp + contxb - 1)%sf(j, k, l) + + ! mass of the two participating fluids + mT = q_cons_vf(lp + contxb - 1)%sf(j, k, l) & + + q_cons_vf(vp + contxb - 1)%sf(j, k, l) + + TS = 1/(mT*cvs(vp)*(gs_min(vp) - 1)/(pS + ps_inf(vp)) & + + ml*(cvs(lp)*(gs_min(lp) - 1)/(pS + ps_inf(lp)) & + - cvs(vp)*(gs_min(vp) - 1)/(pS + ps_inf(vp))) & + + mCVGP) + + dFdT = & + -(cvs(lp)*gs_min(lp) - cvs(vp)*gs_min(vp))*DLOG(TS) & + - (qvps(lp) - qvps(vp)) & + + cvs(lp)*(gs_min(lp) - 1)*DLOG(pS + ps_inf(lp)) & + - cvs(vp)*(gs_min(vp) - 1)*DLOG(pS + ps_inf(vp)) + + dTdm = -(cvs(lp)*(gs_min(lp) - 1)/(pS + ps_inf(lp)) & + - cvs(vp)*(gs_min(vp) - 1)/(pS + ps_inf(vp)))*TS**2 + + dTdp = (mT*cvs(vp)*(gs_min(vp) - 1)/(pS + ps_inf(vp))**2 & + + ml*(cvs(lp)*(gs_min(lp) - 1)/(pS + ps_inf(lp))**2 & + - cvs(vp)*(gs_min(vp) - 1)/(pS + ps_inf(vp))**2) & + + mCVGP2)*TS**2 + + ! F = (F1,F2) is the function whose roots we are looking for + ! x = (m1, p) are the independent variables. m1 = mass of the first participant fluid, p = pressure + ! F1 = 0 is the Gibbs free energy quality + ! F2 = 0 is the enforcement of the thermodynamic (total - kinectic) energy + ! dF1dm + Jac(1, 1) = dFdT*dTdm + + ! dF1dp + Jac(1, 2) = dFdT*dTdp + TS & + *(cvs(lp)*(gs_min(lp) - 1)/(pS + ps_inf(lp)) & + - cvs(vp)*(gs_min(vp) - 1)/(pS + ps_inf(vp))) + + ! dF2dm + Jac(2, 1) = (qvs(vp) - qvs(lp) & + + (cvs(vp)*gs_min(vp) - cvs(lp)*gs_min(lp)) & + /(ml*(cvs(lp)*(gs_min(lp) - 1)/(pS + ps_inf(lp)) & + - cvs(vp)*(gs_min(vp) - 1)/(pS + ps_inf(vp))) & + + mT*cvs(vp)*(gs_min(vp) - 1)/(pS + ps_inf(vp)) + mCVGP) & + - (ml*(cvs(vp)*gs_min(vp) - cvs(lp)*gs_min(lp)) & + - mT*cvs(vp)*gs_min(vp) - mCPD) & + *(cvs(lp)*(gs_min(lp) - 1)/(pS + ps_inf(lp)) & + - cvs(vp)*(gs_min(vp) - 1)/(pS + ps_inf(vp))) & + /((ml*(cvs(lp)*(gs_min(lp) - 1)/(pS + ps_inf(lp)) & + - cvs(vp)*(gs_min(vp) - 1)/(pS + ps_inf(vp))) & + + mT*cvs(vp)*(gs_min(vp) - 1)/(pS + ps_inf(vp)) + mCVGP)**2))/1 + ! dF2dp + Jac(2, 2) = (1 + (ml*(cvs(vp)*gs_min(vp) - cvs(lp)*gs_min(lp)) & + - mT*cvs(vp)*gs_min(vp) - mCPD) & + *(ml*(cvs(lp)*(gs_min(lp) - 1)/(pS + ps_inf(lp))**2 & + - cvs(vp)*(gs_min(vp) - 1)/(pS + ps_inf(vp))**2) & + + mT*cvs(vp)*(gs_min(vp) - 1)/(pS + ps_inf(vp))**2 + mCVGP2) & + /(ml*(cvs(lp)*(gs_min(lp) - 1)/(pS + ps_inf(lp)) & + - cvs(vp)*(gs_min(vp) - 1)/(pS + ps_inf(vp))) & + + mT*cvs(vp)*(gs_min(vp) - 1)/(pS + ps_inf(vp)) + mCVGP)**2)/1 + + ! intermediate elements of J^{-1} + InvJac(1, 1) = Jac(2, 2) + InvJac(1, 2) = -1.0d0*Jac(1, 2) + InvJac(2, 1) = -1.0d0*Jac(2, 1) + InvJac(2, 2) = Jac(1, 1) + + ! elements of J^{T} + TJac(1, 1) = Jac(1, 1) + TJac(1, 2) = Jac(2, 1) + TJac(2, 1) = Jac(1, 2) + TJac(2, 2) = Jac(2, 2) + + ! dividing by det(J) + InvJac = InvJac/(Jac(1, 1)*Jac(2, 2) - Jac(1, 2)*Jac(2, 1)) + + end subroutine s_compute_jacobian_matrix + + !> This auxiliary subroutine computes the residue of the pTg-equilibrium procedure + !! @param j generic loop iterator for x direction + !! @param k generic loop iterator for y direction + !! @param l generic loop iterator for z direction + !! @param mCPD sum of the total alpha*rho*cp + !! @param mCVGP auxiliary variable for the calculation of the matrices: alpha*rho*cv*(g-1)/press + !! @param mQD sum of the total alpha*rho*qv + !! @param q_cons_vf Cell-average conservative variables + !! @param pS equilibrium pressure at the interface + !! @param rhoe mixture energy + !! @param R2D (2D) residue array + subroutine s_compute_pTg_residue(j, k, l, mCPD, mCVGP, mQD, q_cons_vf, pS, rhoe, R2D) +!$acc routine seq + + type(scalar_field), dimension(sys_size), intent(IN) :: q_cons_vf + real(kind(0.0d0)), intent(IN) :: pS, rhoe, mCPD, mCVGP, mQD + integer, intent(IN) :: j, k, l + real(kind(0.0d0)), dimension(2), intent(OUT) :: R2D + real(kind(0.0d0)) :: ml, mT, TS !< mass of the reacting liquid, total reacting mass, equilibrium temperature + + ! mass of the reacting liquid + ml = q_cons_vf(lp + contxb - 1)%sf(j, k, l) + + ! mass of the two participating fluids + mT = q_cons_vf(lp + contxb - 1)%sf(j, k, l) & + + q_cons_vf(vp + contxb - 1)%sf(j, k, l) + + TS = 1/(mT*cvs(vp)*(gs_min(vp) - 1)/(pS + ps_inf(vp)) & + + ml*(cvs(lp)*(gs_min(lp) - 1)/(pS + ps_inf(lp)) & + - cvs(vp)*(gs_min(vp) - 1)/(pS + ps_inf(vp))) & + + mCVGP) + + ! Gibbs Free Energy Equality condition (DG) + R2D(1) = TS*((cvs(lp)*gs_min(lp) - cvs(vp)*gs_min(vp)) & + *(1 - DLOG(TS)) - (qvps(lp) - qvps(vp)) & + + cvs(lp)*(gs_min(lp) - 1)*DLOG(pS + ps_inf(lp)) & + - cvs(vp)*(gs_min(vp) - 1)*DLOG(pS + ps_inf(vp))) & + + qvs(lp) - qvs(vp) + + ! Constant Energy Process condition (DE) + R2D(2) = (rhoe + pS & + + ml*(qvs(vp) - qvs(lp)) - mT*qvs(vp) - mQD & + + (ml*(gs_min(vp)*cvs(vp) - gs_min(lp)*cvs(lp)) & + - mT*gs_min(vp)*cvs(vp) - mCPD) & + /(ml*(cvs(lp)*(gs_min(lp) - 1)/(pS + ps_inf(lp)) & + - cvs(vp)*(gs_min(vp) - 1)/(pS + ps_inf(vp))) & + + mT*cvs(vp)*(gs_min(vp) - 1)/(pS + ps_inf(vp)) + mCVGP))/1 + + end subroutine s_compute_pTg_residue + + !> This auxiliary subroutine finds the Saturation temperature for a given + !! saturation pressure through a newton solver + !! @param pSat Saturation Pressure + !! @param TSat Saturation Temperature + !! @param TSIn equilibrium Temperature + subroutine s_TSat(pSat, TSat, TSIn) +!$acc routine seq + + real(kind(0.0d0)), intent(OUT) :: TSat + real(kind(0.0d0)), intent(IN) :: pSat, TSIn + real(kind(0.0d0)) :: dFdT, FT, Om !< auxiliary variables + + ! Generic loop iterators + integer :: ns + + if ((pSat == 0.0d0) .and. (TSIn == 0.0d0)) then + + ! assigning Saturation temperature + TSat = 0.0d0 + + else + + ! calculating initial estimate for temperature in the TSat procedure. I will also use this variable to + ! iterate over the Newton's solver + TSat = TSIn + + ! iteration counter + ns = 0 + + ! underrelaxation factor + Om = 1.0d-3 + do while ((DABS(FT) > ptgalpha_eps) .or. (ns == 0)) + ! increasing counter + ns = ns + 1 + + ! calculating residual + FT = TSat*((cvs(lp)*gs_min(lp) - cvs(vp)*gs_min(vp)) & + *(1 - DLOG(TSat)) - (qvps(lp) - qvps(vp)) & + + cvs(lp)*(gs_min(lp) - 1)*DLOG(pSat + ps_inf(lp)) & + - cvs(vp)*(gs_min(vp) - 1)*DLOG(pSat + ps_inf(vp))) & + + qvs(lp) - qvs(vp) + + ! calculating the jacobian + dFdT = & + -(cvs(lp)*gs_min(lp) - cvs(vp)*gs_min(vp))*DLOG(TSat) & + - (qvps(lp) - qvps(vp)) & + + cvs(lp)*(gs_min(lp) - 1)*DLOG(pSat + ps_inf(lp)) & + - cvs(vp)*(gs_min(vp) - 1)*DLOG(pSat + ps_inf(vp)) + + ! updating saturation temperature + TSat = TSat - Om * FT/dFdT + + end do + + end if + + end subroutine s_TSat + + !> This subroutine finalizes the phase change module + subroutine s_finalize_relaxation_solver_module() + + s_relaxation_solver => null() + + end subroutine + +#endif + +end module m_phase_change \ No newline at end of file diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index ba8c79726..5b1056449 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -53,8 +53,9 @@ module m_variables_conversion !! @param rho Density !! @param gamma Specific heat ratio function !! @param pi_inf Liquid stiffness function + !! @param qv Fluid reference energy subroutine s_convert_xxxxx_to_mixture_variables(q_vf, i, j, k, & - rho, gamma, pi_inf, Re_K, G_K, G) + rho, gamma, pi_inf, qv, Re_K, G_K, G) ! Importing the derived type scalar_field from m_derived_types.f90 ! and global variable sys_size, from m_global_variables.f90, as @@ -68,6 +69,7 @@ module m_variables_conversion real(kind(0d0)), intent(OUT), target :: rho real(kind(0d0)), intent(OUT), target :: gamma real(kind(0d0)), intent(OUT), target :: pi_inf + real(kind(0d0)), intent(OUT), target :: qv real(kind(0d0)), optional, dimension(2), intent(OUT) :: Re_K @@ -81,10 +83,10 @@ module m_variables_conversion integer, public :: ixb, ixe, iyb, iye, izb, ize !$acc declare create(ixb, ixe, iyb, iye, izb, ize) - !! In simulation, gammas and pi_infs is already declared in m_global_variables + !! In simulation, gammas, pi_infs, and qvs are already declared in m_global_variables #ifndef MFC_SIMULATION - real(kind(0d0)), allocatable, public, dimension(:) :: gammas, pi_infs - !$acc declare create(gammas, pi_infs) + real(kind(0d0)), allocatable, public, dimension(:) :: gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps + !$acc declare create(gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps) #endif real(kind(0d0)), allocatable, dimension(:) :: Gs @@ -97,7 +99,8 @@ module m_variables_conversion real(kind(0d0)), allocatable, dimension(:, :, :), public :: rho_sf !< Scalar density function real(kind(0d0)), allocatable, dimension(:, :, :), public :: gamma_sf !< Scalar sp. heat ratio function - real(kind(0d0)), allocatable, dimension(:, :, :), public :: pi_inf_sf !< Scalar liquid stiffness function + real(kind(0d0)), allocatable, dimension(:, :, :), public :: pi_inf_sf !< Scalar liquid stiffness function + real(kind(0d0)), allocatable, dimension(:, :, :), public :: qv_sf !< Scalar liquid energy reference function procedure(s_convert_xxxxx_to_mixture_variables), & pointer :: s_convert_to_mixture_variables => null() !< @@ -113,9 +116,10 @@ contains !! @param mom Momentum !! @param dyn_p Dynamic Pressure !! @param pi_inf Liquid Stiffness + !! @param qv fluid reference energy !! @param gamma Specific Heat Ratio !! @param pres Pressure to calculate - subroutine s_compute_pressure(energy, alf, dyn_p, pi_inf, gamma, rho, pres, stress, mom, G) + subroutine s_compute_pressure(energy, alf, dyn_p, pi_inf, gamma, rho, qv, pres, stress, mom, G) !$acc routine seq real(kind(0d0)), intent(IN) :: energy, alf @@ -124,7 +128,7 @@ contains real(kind(0d0)), intent(IN) :: dyn_p real(kind(0d0)), intent(OUT) :: pres - real(kind(0d0)), intent(IN) :: pi_inf, gamma, rho + real(kind(0d0)), intent(IN) :: pi_inf, gamma, rho, qv real(kind(0d0)) :: E_e @@ -134,9 +138,9 @@ contains ! for computing pressure is targeted by the procedure pointer if ((model_eqns /= 4) .and. (bubbles .neqv. .true.)) then - pres = (energy - dyn_p - pi_inf)/gamma + pres = (energy - dyn_p - pi_inf - qv)/gamma else if ((model_eqns /= 4) .and. bubbles) then - pres = ((energy - dyn_p)/(1.d0 - alf) - pi_inf)/gamma + pres = ((energy - dyn_p)/(1.d0 - alf) - pi_inf - qv)/gamma else pres = (pref + pi_inf)* & (energy/ & @@ -162,7 +166,7 @@ contains pres = ( & energy - & 0.5d0*(mom**2.d0)/rho - & - pi_inf - E_e & + pi_inf - qv - E_e & )/gamma @@ -185,7 +189,7 @@ contains !! @param gamma specific heat ratio function !! @param pi_inf liquid stiffness subroutine s_convert_mixture_to_mixture_variables(q_vf, i, j, k, & - rho, gamma, pi_inf, Re_K, G_K, G) + rho, gamma, pi_inf, qv, Re_K, G_K, G) type(scalar_field), dimension(sys_size), intent(IN) :: q_vf @@ -194,6 +198,7 @@ contains real(kind(0d0)), intent(OUT), target :: rho real(kind(0d0)), intent(OUT), target :: gamma real(kind(0d0)), intent(OUT), target :: pi_inf + real(kind(0d0)), intent(OUT), target :: qv real(kind(0d0)), optional, dimension(2), intent(OUT) :: Re_K @@ -205,13 +210,15 @@ contains rho = q_vf(1)%sf(i, j, k) gamma = q_vf(gamma_idx)%sf(i, j, k) pi_inf = q_vf(pi_inf_idx)%sf(i, j, k) + qv = 0d0 ! keep this value nill for now. For future adjustment - ! Post process requires rho_sf/gamma_sf/pi_inf_sf to also be updated + ! Post process requires rho_sf/gamma_sf/pi_inf_sf/qv_sf to also be updated #ifdef MFC_POST_PROCESS rho_sf (i, j, k) = rho gamma_sf (i, j, k) = gamma pi_inf_sf(i, j, k) = pi_inf + qv_sf(i, j, k) = qv #endif end subroutine s_convert_mixture_to_mixture_variables ! ---------------- @@ -226,11 +233,12 @@ contains !! @param rho_K density !! @param gamma_K specific heat ratio !! @param pi_inf_K liquid stiffness + !! @param qv_K fluid reference energy !! @param j Cell index !! @param k Cell index !! @param l Cell index subroutine s_convert_species_to_mixture_variables_bubbles(q_vf, j, k, l, & - rho, gamma, pi_inf, Re_K, G_K, G) + rho, gamma, pi_inf, qv, Re_K, G_K, G) type(scalar_field), dimension(sys_size), intent(IN) :: q_vf @@ -239,6 +247,7 @@ contains real(kind(0d0)), intent(OUT), target :: rho real(kind(0d0)), intent(OUT), target :: gamma real(kind(0d0)), intent(OUT), target :: pi_inf + real(kind(0d0)), intent(OUT), target :: qv real(kind(0d0)), dimension(num_fluids) :: alpha_rho_K, alpha_K real(kind(0d0)), optional, dimension(2), intent(OUT) :: Re_K @@ -275,25 +284,29 @@ contains rho = q_vf(1)%sf(j, k, l) gamma = fluid_pp(1)%gamma !qK_vf(gamma_idx)%sf(i,j,k) pi_inf = fluid_pp(1)%pi_inf !qK_vf(pi_inf_idx)%sf(i,j,k) + qv = fluid_pp(1)%qv else if ((model_eqns == 2) .and. bubbles) then - rho = 0d0; gamma = 0d0; pi_inf = 0d0 + rho = 0d0; gamma = 0d0; pi_inf = 0d0; qv = 0d0 if (mpp_lim .and. (num_fluids > 2)) then do i = 1, num_fluids rho = rho + q_vf(i)%sf(j, k, l) gamma = gamma + q_vf(i + E_idx)%sf(j, k, l)*fluid_pp(i)%gamma pi_inf = pi_inf + q_vf(i + E_idx)%sf(j, k, l)*fluid_pp(i)%pi_inf + qv = qv + q_vf(i)%sf(j, k, l)*fluid_pp(i)%qv end do else if (num_fluids == 2) then rho = q_vf(1)%sf(j, k, l) gamma = fluid_pp(1)%gamma pi_inf = fluid_pp(1)%pi_inf + qv = fluid_pp(1)%qv else if (num_fluids > 2) then !TODO: This may need fixing for hypo + bubbles do i = 1, num_fluids - 1 !leave out bubble part of mixture rho = rho + q_vf(i)%sf(j, k, l) gamma = gamma + q_vf(i + E_idx)%sf(j, k, l)*fluid_pp(i)%gamma pi_inf = pi_inf + q_vf(i + E_idx)%sf(j, k, l)*fluid_pp(i)%pi_inf + qv = qv + q_vf(i)%sf(j, k, l)*fluid_pp(i)%qv end do !rho = qK_vf(1)%sf(j,k,l) !gamma_K = fluid_pp(1)%gamma @@ -302,6 +315,7 @@ contains rho = q_vf(1)%sf(j, k, l) gamma = fluid_pp(1)%gamma pi_inf = fluid_pp(1)%pi_inf + qv = fluid_pp(1)%qv end if end if @@ -325,11 +339,12 @@ contains end if #endif - ! Post process requires rho_sf/gamma_sf/pi_inf_sf to also be updated + ! Post process requires rho_sf/gamma_sf/pi_inf_sf/qv_sf to also be updated #ifdef MFC_POST_PROCESS rho_sf (j, k, l) = rho gamma_sf (j, k, l) = gamma pi_inf_sf(j, k, l) = pi_inf + qv_sf(j, k, l) = qv #endif end subroutine s_convert_species_to_mixture_variables_bubbles ! ---------------- @@ -343,11 +358,12 @@ contains !! @param rho density !! @param gamma specific heat ratio !! @param pi_inf liquid stiffness + !! @param fluid reference energy !! @param j Cell index !! @param k Cell index !! @param l Cell index - subroutine s_convert_species_to_mixture_variables(q_vf, k, l, r, & - rho, gamma, pi_inf, Re_K, G_K, G) + subroutine s_convert_species_to_mixture_variables(q_vf, k, l, r, rho, & + gamma, pi_inf, qv, Re_K, G_K, G) type(scalar_field), dimension(sys_size), intent(IN) :: q_vf @@ -356,6 +372,7 @@ contains real(kind(0d0)), intent(OUT), target :: rho real(kind(0d0)), intent(OUT), target :: gamma real(kind(0d0)), intent(OUT), target :: pi_inf + real(kind(0d0)), intent(OUT), target :: qv real(kind(0d0)), optional, dimension(2), intent(OUT) :: Re_K @@ -386,14 +403,16 @@ contains end if - ! Calculating the density, the specific heat ratio function and the - ! liquid stiffness function, respectively, from the species analogs - rho = 0d0; gamma = 0d0; pi_inf = 0d0 + ! Calculating the density, the specific heat ratio function, the + ! liquid stiffness function, and the energy reference function, + ! respectively, from the species analogs + rho = 0d0; gamma = 0d0; pi_inf = 0d0; qv = 0d0 do i = 1, num_fluids rho = rho + alpha_rho_K(i) gamma = gamma + alpha_K(i)*gammas(i) pi_inf = pi_inf + alpha_K(i)*pi_infs(i) + qv = qv + alpha_rho_K(i)*qvs(i) end do #ifdef MFC_SIMULATION @@ -420,22 +439,23 @@ contains G_K = max(0d0, G_K) end if - ! Post process requires rho_sf/gamma_sf/pi_inf_sf to also be updated + ! Post process requires rho_sf/gamma_sf/pi_inf_sf/qv_sf to also be updated #ifdef MFC_POST_PROCESS rho_sf (k, l, r) = rho gamma_sf (k, l, r) = gamma pi_inf_sf(k, l, r) = pi_inf + qv_sf(k, l, r) = qv #endif end subroutine s_convert_species_to_mixture_variables ! ---------------- subroutine s_convert_species_to_mixture_variables_acc(rho_K, & - gamma_K, pi_inf_K, & + gamma_K, pi_inf_K, qv_K, & alpha_K, alpha_rho_K, Re_K, k, l, r, & G_K, G) !$acc routine seq - real(kind(0d0)), intent(OUT) :: rho_K, gamma_K, pi_inf_K + real(kind(0d0)), intent(OUT) :: rho_K, gamma_K, pi_inf_K, qv_K real(kind(0d0)), dimension(num_fluids), intent(INOUT) :: alpha_rho_K, alpha_K !< real(kind(0d0)), dimension(2), intent(OUT) :: Re_K @@ -457,6 +477,7 @@ contains rho_K = 0d0 gamma_K = 0d0 pi_inf_K = 0d0 + qv_K = 0d0 alpha_K_sum = 0d0 @@ -475,6 +496,7 @@ contains rho_K = rho_K + alpha_rho_K(i) gamma_K = gamma_K + alpha_K(i)*gammas(i) pi_inf_K = pi_inf_K + alpha_K(i)*pi_infs(i) + qv_K = qv_K + alpha_rho_K(i)*qvs(i) end do if (present(G_K)) then @@ -507,11 +529,11 @@ contains end subroutine s_convert_species_to_mixture_variables_acc ! ---------------- subroutine s_convert_species_to_mixture_variables_bubbles_acc(rho_K, & - gamma_K, pi_inf_K, & + gamma_K, pi_inf_K, qv_K, & alpha_K, alpha_rho_K, Re_K, k, l, r) !$acc routine seq - real(kind(0d0)), intent(INOUT) :: rho_K, gamma_K, pi_inf_K + real(kind(0d0)), intent(INOUT) :: rho_K, gamma_K, pi_inf_K, qv_K real(kind(0d0)), dimension(num_fluids), intent(IN) :: alpha_rho_K, alpha_K !< !! Partial densities and volume fractions @@ -525,23 +547,27 @@ contains rho_K = 0d0 gamma_K = 0d0 pi_inf_K = 0d0 + qv_K = 0d0 if (mpp_lim .and. (model_eqns == 2) .and. (num_fluids > 2)) then do i = 1, num_fluids rho_K = rho_K + alpha_rho_K(i) gamma_K = gamma_K + alpha_K(i)*gammas(i) pi_inf_K = pi_inf_K + alpha_K(i)*pi_infs(i) + qv_K = qv_K + alpha_rho_K(i)*qvs(i) end do else if ((model_eqns == 2) .and. (num_fluids > 2)) then do i = 1, num_fluids - 1 rho_K = rho_K + alpha_rho_K(i) gamma_K = gamma_K + alpha_K(i)*gammas(i) pi_inf_K = pi_inf_K + alpha_K(i)*pi_infs(i) + qv_K = qv_K + alpha_rho_K(i)*qvs(i) end do else rho_K = alpha_rho_K(1) gamma_K = gammas(1) pi_inf_K = pi_infs(1) + qv_K = qvs(1) end if if (any(Re_size > 0)) then @@ -593,15 +619,25 @@ contains !$acc update device(ixb, ixe, iyb, iye, izb, ize) @:ALLOCATE(gammas (1:num_fluids)) + @:ALLOCATE(gs_min (1:num_fluids)) @:ALLOCATE(pi_infs(1:num_fluids)) + @:ALLOCATE(ps_inf(1:num_fluids)) + @:ALLOCATE(cvs (1:num_fluids)) + @:ALLOCATE(qvs (1:num_fluids)) + @:ALLOCATE(qvps (1:num_fluids)) @:ALLOCATE(Gs (1:num_fluids)) do i = 1, num_fluids gammas(i) = fluid_pp(i)%gamma + gs_min(i) = 1.0d0/gammas(i) + 1.0d0 pi_infs(i) = fluid_pp(i)%pi_inf Gs(i) = fluid_pp(i)%G + ps_inf(i) = pi_infs(i)/(1.0d0 + gammas(i)) + cvs(i) = fluid_pp(i)%cv + qvs(i) = fluid_pp(i)%qv + qvps(i) = fluid_pp(i)%qvp end do - !$acc update device(gammas, pi_infs, Gs) + !$acc update device(gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps, Gs) #ifdef MFC_SIMULATION @@ -647,6 +683,9 @@ contains allocate (pi_inf_sf(-buff_size:m + buff_size, & -buff_size:n + buff_size, & -buff_size:p + buff_size)) + allocate (qv_sf(-buff_size:m + buff_size, & + -buff_size:n + buff_size, & + -buff_size:p + buff_size)) ! Simulation is 2D else @@ -660,7 +699,9 @@ contains allocate (pi_inf_sf(-buff_size:m + buff_size, & -buff_size:n + buff_size, & 0:0)) - + allocate (qv_sf(-buff_size:m + buff_size, & + -buff_size:n + buff_size, & + 0:0)) end if ! Simulation is 1D @@ -675,6 +716,9 @@ contains allocate (pi_inf_sf(-buff_size:m + buff_size, & 0:0, & 0:0)) + allocate (qv_sf(-buff_size:m + buff_size, & + 0:0, & + 0:0)) end if #endif @@ -785,7 +829,7 @@ contains real(kind(0d0)), dimension(num_fluids) :: alpha_K, alpha_rho_K real(kind(0d0)), dimension(2) :: Re_K - real(kind(0d0)) :: rho_K, gamma_K, pi_inf_K, dyn_pres_K + real(kind(0d0)) :: rho_K, gamma_K, pi_inf_K, qv_K, dyn_pres_K #:if MFC_CASE_OPTIMIZATION #ifndef MFC_SIMULATION @@ -823,7 +867,7 @@ contains endif #:endif - !$acc parallel loop collapse(3) gang vector default(present) private(alpha_K, alpha_rho_K, Re_K, nRtmp, rho_K, gamma_K, pi_inf_K, dyn_pres_K, R3tmp) + !$acc parallel loop collapse(3) gang vector default(present) private(alpha_K, alpha_rho_K, Re_K, nRtmp, rho_K, gamma_K, pi_inf_K, qv_K, dyn_pres_K, R3tmp) do l = izb, ize do k = iyb, iye do j = ixb, ixe @@ -846,23 +890,23 @@ contains #ifdef MFC_SIMULATION ! If in simulation, use acc mixture subroutines if (hypoelasticity) then - call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, alpha_K, & + call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, qv_K, alpha_K, & alpha_rho_K, Re_K, j, k, l, G_K, Gs) else if (bubbles) then - call s_convert_species_to_mixture_variables_bubbles_acc(rho_K, gamma_K, pi_inf_K, & + call s_convert_species_to_mixture_variables_bubbles_acc(rho_K, gamma_K, pi_inf_K, qv_K, & alpha_K, alpha_rho_K, Re_K, j, k, l) else - call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, & + call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, qv_K, & alpha_K, alpha_rho_K, Re_K, j, k, l) end if #else ! If pre-processing, use non acc mixture subroutines if (hypoelasticity) then call s_convert_to_mixture_variables(qK_cons_vf, j, k, l, & - rho_K, gamma_K, pi_inf_K, Re_K, G_K, fluid_pp(:)%G) + rho_K, gamma_K, pi_inf_K, qv_K, Re_K, G_K, fluid_pp(:)%G) else call s_convert_to_mixture_variables(qK_cons_vf, j, k, l, & - rho_K, gamma_K, pi_inf_K) + rho_K, gamma_K, pi_inf_K, qv_K) end if #endif end if @@ -886,7 +930,7 @@ contains call s_compute_pressure(qK_cons_vf(E_idx)%sf(j, k, l), & qK_cons_vf(alf_idx)%sf(j, k, l), & - dyn_pres_K, pi_inf_K, gamma_K, rho_K, pres) + dyn_pres_K, pi_inf_K, gamma_K, rho_K, qv_K, pres) qK_prim_vf(E_idx)%sf(j, k, l) = pres @@ -978,11 +1022,12 @@ contains real(kind(0d0)) :: rho real(kind(0d0)) :: gamma real(kind(0d0)) :: pi_inf + real(kind(0d0)) :: qv real(kind(0d0)) :: dyn_pres real(kind(0d0)) :: nbub, R3, vftmp, R3tmp real(kind(0d0)), dimension(nb) :: Rtmp - real(kind(0d0)) :: G + real(kind(0d0)), dimension(2) :: Re_K integer :: i, j, k, l, q !< Generic loop iterators @@ -995,7 +1040,7 @@ contains ! Obtaining the density, specific heat ratio function ! and the liquid stiffness function, respectively call s_convert_to_mixture_variables(q_prim_vf, j, k, l, & - rho, gamma, pi_inf) + rho, gamma, pi_inf, qv, Re_K, G, fluid_pp(:)%G) ! Transferring the continuity equation(s) variable(s) do i = 1, contxe @@ -1020,14 +1065,15 @@ contains ! Computing the energy from the pressure if ((model_eqns /= 4) .and. (bubbles .neqv. .true.)) then - ! E = Gamma*P + \rho u u /2 + \pi_inf + ! E = Gamma*P + \rho u u /2 + \pi_inf + (\alpha\rho qv) q_cons_vf(E_idx)%sf(j, k, l) = & - gamma*q_prim_vf(E_idx)%sf(j, k, l) + dyn_pres + pi_inf + gamma*q_prim_vf(E_idx)%sf(j, k, l) + dyn_pres + pi_inf & + + qv else if ((model_eqns /= 4) .and. (bubbles)) then ! \tilde{E} = dyn_pres + (1-\alf)(\Gamma p_l + \Pi_inf) q_cons_vf(E_idx)%sf(j, k, l) = dyn_pres + & (1.d0 - q_prim_vf(alf_idx)%sf(j, k, l))* & - (gamma*q_prim_vf(E_idx)%sf(j, k, l) + pi_inf) + (gamma*q_prim_vf(E_idx)%sf(j, k, l) + pi_inf) else !Tait EOS, no conserved energy variable q_cons_vf(E_idx)%sf(j, k, l) = 0. @@ -1040,7 +1086,8 @@ contains q_cons_vf(i + internalEnergies_idx%beg - 1)%sf(j, k, l) = & q_cons_vf(i + adv_idx%beg - 1)%sf(j, k, l)* & (fluid_pp(i)%gamma*q_prim_vf(E_idx)%sf(j, k, l) + & - fluid_pp(i)%pi_inf) + fluid_pp(i)%pi_inf) + & + q_cons_vf(i + cont_idx%beg - 1)%sf(j, k, l)*fluid_pp(i)%qv end do end if @@ -1135,6 +1182,7 @@ contains real(kind(0d0)) :: E_K real(kind(0d0)) :: gamma_K real(kind(0d0)) :: pi_inf_K + real(kind(0d0)) :: qv_K real(kind(0d0)), dimension(2) :: Re_K real(kind(0d0)) :: G_K @@ -1176,20 +1224,20 @@ contains pres_K = qK_prim_vf(j, k, l, E_idx) if (hypoelasticity) then - call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, & + call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, qv_K, & alpha_K, alpha_rho_K, Re_K, & j, k, l, G_K, Gs) else if (bubbles) then call s_convert_species_to_mixture_variables_bubbles_acc(rho_K, gamma_K, & - pi_inf_K, alpha_K, alpha_rho_K, Re_K, j, k, l) + pi_inf_K, qv_K, alpha_K, alpha_rho_K, Re_K, j, k, l) else - call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, & + call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, qv_K, & alpha_K, alpha_rho_K, Re_K, j, k, l) end if ! Computing the energy from the pressure E_K = gamma_K*pres_K + pi_inf_K & - + 5d-1*rho_K*vel_K_sum + + 5d-1*rho_K*vel_K_sum + qv_K ! mass flux, this should be \alpha_i \rho_i u_i !$acc loop seq @@ -1241,10 +1289,10 @@ contains ! Deallocating the density, the specific heat ratio function and the ! liquid stiffness function #ifdef MFC_POST_PROCESS - deallocate(rho_sf, gamma_sf, pi_inf_sf) + deallocate(rho_sf, gamma_sf, pi_inf_sf, qv_sf) #endif - @:DEALLOCATE(gammas, pi_infs, Gs) + @:DEALLOCATE(gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps, Gs) if (bubbles) then @:DEALLOCATE(bubrs) diff --git a/src/post_process/m_global_parameters.fpp b/src/post_process/m_global_parameters.fpp index 87584fc21..ad1bfa1b1 100644 --- a/src/post_process/m_global_parameters.fpp +++ b/src/post_process/m_global_parameters.fpp @@ -82,6 +82,8 @@ module m_global_parameters !> @{ integer :: model_eqns !< Multicomponent flow model integer :: num_fluids !< Number of different fluids present in the flow + logical :: relax !< phase change + integer :: relax_model !< Phase change relaxation model logical :: adv_alphan !< Advection of the last volume fraction logical :: mpp_lim !< Maximum volume fraction limiter integer :: sys_size !< Number of unknowns in the system of equations @@ -270,6 +272,8 @@ contains weno_order = dflt_int mixture_err = .false. alt_soundspeed = .false. + relax = .false. + relax_model = dflt_int hypoelasticity = .false. bc_x%beg = dflt_int; bc_x%end = dflt_int @@ -280,6 +284,9 @@ contains do i = 1, num_fluids_max fluid_pp(i)%gamma = dflt_real fluid_pp(i)%pi_inf = dflt_real + fluid_pp(i)%cv = 0d0 + fluid_pp(i)%qv = 0d0 + fluid_pp(i)%qvp = 0d0 fluid_pp(i)%G = dflt_real end do diff --git a/src/post_process/m_mpi_proxy.fpp b/src/post_process/m_mpi_proxy.fpp index 4d2b8ebce..8250fb47d 100644 --- a/src/post_process/m_mpi_proxy.fpp +++ b/src/post_process/m_mpi_proxy.fpp @@ -159,7 +159,7 @@ contains & 't_step_start', 't_step_stop', 't_step_save', 'weno_order', & & 'model_eqns', 'num_fluids', 'bc_x%beg', 'bc_x%end', 'bc_y%beg', & & 'bc_y%end', 'bc_z%beg', 'bc_z%end', 'flux_lim', 'format', & - & 'precision', 'fd_order', 'thermal', 'nb' ] + & 'precision', 'fd_order', 'thermal', 'nb', 'relax_model' ] call MPI_BCAST(${VAR}$, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) #:endfor @@ -168,7 +168,7 @@ contains & 'E_wrt', 'pres_wrt', 'gamma_wrt', & & 'heat_ratio_wrt', 'pi_inf_wrt', 'pres_inf_wrt', 'cons_vars_wrt', & & 'prim_vars_wrt', 'c_wrt', 'qm_wrt','schlieren_wrt', 'bubbles', & - & 'polytropic', 'polydisperse', 'file_per_process' ] + & 'polytropic', 'polydisperse', 'file_per_process', 'relax' ] call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) #:endfor @@ -182,6 +182,9 @@ contains do i = 1, num_fluids_max call MPI_BCAST(fluid_pp(i)%gamma, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) call MPI_BCAST(fluid_pp(i)%pi_inf, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(fluid_pp(i)%cv, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(fluid_pp(i)%qv, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(fluid_pp(i)%qvp, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) call MPI_BCAST(fluid_pp(i)%G, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) end do diff --git a/src/post_process/m_start_up.f90 b/src/post_process/m_start_up.f90 index 120156a61..4d5e07e18 100644 --- a/src/post_process/m_start_up.f90 +++ b/src/post_process/m_start_up.f90 @@ -71,7 +71,8 @@ subroutine s_read_input_file() ! --------------------------------------- flux_lim, flux_wrt, cyl_coord, & parallel_io, rhoref, pref, bubbles, qbmm, sigR, & R0ref, nb, polytropic, thermal, Ca, Web, Re_inv, & - polydisperse, poly_sigma, file_per_process + polydisperse, poly_sigma, file_per_process, relax, & + relax_model ! Inquiring the status of the post_process.inp file file_loc = 'post_process.inp' diff --git a/src/pre_process/m_assign_variables.f90 b/src/pre_process/m_assign_variables.f90 index 8390731ad..c11f9fd4e 100644 --- a/src/pre_process/m_assign_variables.f90 +++ b/src/pre_process/m_assign_variables.f90 @@ -255,9 +255,11 @@ subroutine s_assign_patch_species_primitive_variables(patch_id, j, k, l, & real(kind(0d0)) :: gamma real(kind(0d0)) :: lit_gamma !< specific heat ratio real(kind(0d0)) :: pi_inf !< stiffness from SEOS + real(kind(0d0)) :: qv !< reference energy from SEOS real(kind(0d0)) :: orig_rho real(kind(0d0)) :: orig_gamma real(kind(0d0)) :: orig_pi_inf + real(kind(0d0)) :: orig_qv real(kind(0d0)) :: muR, muV real(kind(0d0)), dimension(int(E_idx - mom_idx%beg)) :: vel !< velocity @@ -298,7 +300,7 @@ subroutine s_assign_patch_species_primitive_variables(patch_id, j, k, l, & q_prim_vf, j, k, l, & orig_rho, & orig_gamma, & - orig_pi_inf) + orig_pi_inf, orig_qv) ! Computing Mixture Variables of Current Patch ===================== @@ -333,7 +335,8 @@ subroutine s_assign_patch_species_primitive_variables(patch_id, j, k, l, & q_prim_vf, j, k, l, & patch_icpp(patch_id)%rho, & patch_icpp(patch_id)%gamma, & - patch_icpp(patch_id)%pi_inf) + patch_icpp(patch_id)%pi_inf, & + patch_icpp(patch_id)%qv ) ! ================================================================== @@ -403,7 +406,8 @@ subroutine s_assign_patch_species_primitive_variables(patch_id, j, k, l, & q_prim_vf, j, k, l, & patch_icpp(smooth_patch_id)%rho, & patch_icpp(smooth_patch_id)%gamma, & - patch_icpp(smooth_patch_id)%pi_inf) + patch_icpp(smooth_patch_id)%pi_inf, & + patch_icpp(smooth_patch_id)%qv ) ! ================================================================== @@ -464,7 +468,7 @@ subroutine s_assign_patch_species_primitive_variables(patch_id, j, k, l, & ! Density and the specific heat ratio and liquid stiffness functions ! call s_convert_species_to_mixture_variables(q_prim_vf, j, k, l, & call s_convert_to_mixture_variables(q_prim_vf, j, k, l, & - rho, gamma, pi_inf) + rho, gamma, pi_inf, qv) ! Velocity do i = 1, E_idx - mom_idx%beg diff --git a/src/pre_process/m_checker.f90 b/src/pre_process/m_checker.f90 index 9579e6c0b..7f8d9b34c 100644 --- a/src/pre_process/m_checker.f90 +++ b/src/pre_process/m_checker.f90 @@ -91,8 +91,25 @@ subroutine s_check_inputs() call s_mpi_abort('hypoelasticity requires model_eqns = 2'// & 'exiting ...') end if - - ! Constraints on the use of a preexisting grid and initial condition + ! phase change checkers. + if ( relax ) then + if ( model_eqns /= 3 ) then + call s_mpi_abort( 'phase change requires model_eqns = 3. ' // & + 'Exiting ...') + elseif ( ( relax_model .lt. 0 ) .or. ( relax_model .gt. 6 ) ) then + call s_mpi_abort( 'relax_model should be in between 0 and 6. ' // & + 'Exiting ...' ) + elseif ( ( palpha_eps .le. 0d0 ) .or. ( palpha_eps .ge. 1d0 ) .or. & + ( ptgalpha_eps .le. 0d0 ) .or. ( ptgalpha_eps .ge. 1d0 ) ) then + call s_mpi_abort( 'both palpha_eps and ptgalpha_eps must & + be in (0,1). ' // 'Exiting ...') + end if + + elseif ( ( relax_model /= dflt_int ) .or. ( palpha_eps /= dflt_real ) & + .or. ( ptgalpha_eps /= dflt_real ) ) then + call s_mpi_abort( 'relax is not set as true, but other phase change parameters have & + been modified. Either activate phase change or set the values to default. ' // 'Exiting ...') + end if if ((old_grid .neqv. .true.) .and. old_ic) then call s_mpi_abort('Unsupported choice of the combination of '// & 'values for old_grid and old_ic. Exiting ...') @@ -641,6 +658,10 @@ subroutine s_check_inputs() 'of values of num_fluids '// & 'and fluid_pp('//trim(iStr)//')%'// & 'pi_inf. Exiting ...') + elseif ( fluid_pp(i)%cv < 0d0 ) then + call s_mpi_abort('Unsupported value of '// & + 'fluid_pp('//trim(iStr)//')%'// & + 'cv. Make sure cv is positive. Exiting ...') end if end do diff --git a/src/pre_process/m_data_output.fpp b/src/pre_process/m_data_output.fpp index c10fc76e5..3782494f7 100644 --- a/src/pre_process/m_data_output.fpp +++ b/src/pre_process/m_data_output.fpp @@ -93,7 +93,7 @@ contains real(kind(0d0)), dimension(nb) :: nRtmp !< Temporary bubble concentration real(kind(0d0)) :: nbub !< Temporary bubble number density - real(kind(0d0)) :: gamma, lit_gamma, pi_inf !< Temporary EOS params + real(kind(0d0)) :: gamma, lit_gamma, pi_inf, qv !< Temporary EOS params real(kind(0d0)) :: rho !< Temporary density real(kind(0d0)) :: pres !< Temporary pressure @@ -178,6 +178,7 @@ contains gamma = fluid_pp(1)%gamma lit_gamma = 1d0/fluid_pp(1)%gamma + 1d0 pi_inf = fluid_pp(1)%pi_inf + qv = fluid_pp(1)%qv if (precision == 1) then FMT = "(2F30.3)" @@ -200,7 +201,7 @@ contains open (2, FILE=trim(file_loc)) do j = 0, m - call s_convert_to_mixture_variables(q_cons_vf, j, 0, 0, rho, gamma, pi_inf) + call s_convert_to_mixture_variables(q_cons_vf, j, 0, 0, rho, gamma, pi_inf, qv ) lit_gamma = 1d0/gamma + 1d0 @@ -218,7 +219,7 @@ contains q_cons_vf(E_idx)%sf(j, 0, 0), & q_cons_vf(alf_idx)%sf(j, 0, 0), & 0.5d0*(q_cons_vf(mom_idx%beg)%sf(j, 0, 0)**2.d0)/rho, & - pi_inf, gamma, rho, pres) + pi_inf, gamma, rho, qv, pres) write (2, FMT) x_cb(j), pres else if ((i >= bub_idx%beg) .and. (i <= bub_idx%end) .and. bubbles) then diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index 8d23c22ac..84df94b63 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/src/pre_process/m_global_parameters.fpp @@ -72,6 +72,10 @@ module m_global_parameters ! Simulation Algorithm Parameters ========================================== integer :: model_eqns !< Multicomponent flow model + logical :: relax !< activate phase change + integer :: relax_model !< Relax Model + real(kind(0d0)) :: palpha_eps !< trigger parameter for the p relaxation procedure, phase change model + real(kind(0d0)) :: ptgalpha_eps !< trigger parameter for the pTg relaxation procedure, phase change model integer :: num_fluids !< Number of different fluids present in the flow logical :: adv_alphan !< Advection of the last volume fraction logical :: mpp_lim !< Alpha limiter @@ -240,6 +244,10 @@ contains ! Simulation algorithm parameters model_eqns = dflt_int + relax = .false. + relax_model= dflt_int + palpha_eps = dflt_real + ptgalpha_eps = dflt_real num_fluids = dflt_int adv_alphan = .false. weno_order = dflt_int @@ -294,6 +302,9 @@ contains patch_icpp(i)%alpha = dflt_real patch_icpp(i)%gamma = dflt_real patch_icpp(i)%pi_inf = dflt_real + patch_icpp(i)%cv = 0d0 + patch_icpp(i)%qv = 0d0 + patch_icpp(i)%qvp = 0d0 patch_icpp(i)%tau_e = 0d0 !should get all of r0's and v0's patch_icpp(i)%r0 = dflt_real @@ -349,6 +360,9 @@ contains fluid_pp(i)%M_v = dflt_real fluid_pp(i)%mu_v = dflt_real fluid_pp(i)%k_v = dflt_real + fluid_pp(i)%cv = 0d0 + fluid_pp(i)%qv = 0d0 + fluid_pp(i)%qvp = 0d0 fluid_pp(i)%G = 0d0 end do diff --git a/src/pre_process/m_mpi_proxy.fpp b/src/pre_process/m_mpi_proxy.fpp index bda163ca5..b8a091029 100644 --- a/src/pre_process/m_mpi_proxy.fpp +++ b/src/pre_process/m_mpi_proxy.fpp @@ -48,12 +48,12 @@ contains & 'loops_x', 'loops_y', 'loops_z', 'model_eqns', 'num_fluids', & & 'weno_order', 'precision', 'perturb_flow_fluid', & & 'perturb_sph_fluid', 'num_patches', 'thermal', 'nb', 'dist_type',& - & 'R0_type' ] + & 'R0_type', 'relax_model' ] call MPI_BCAST(${VAR}$, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) #:endfor #:for VAR in [ 'old_grid','old_ic','stretch_x','stretch_y','stretch_z',& - & 'cyl_coord','adv_alphan','mpp_lim','hypoelasticity', & + & 'cyl_coord','adv_alphan','mpp_lim','hypoelasticity', 'relax', & & 'parallel_io', 'perturb_flow', 'vel_profile', 'instability_wave', 'perturb_sph', & 'bubbles', 'polytropic', 'polydisperse', 'qbmm', 'file_per_process' ] call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) @@ -66,7 +66,8 @@ contains & 'a_z', 'x_a', 'x_b', 'y_a', 'y_b', 'z_a', 'z_b', 'bc_x%beg', & & 'bc_x%end', 'bc_y%beg', 'bc_y%end', 'bc_z%beg', 'bc_z%end', & & 'perturb_flow_mag', 'pref', 'rhoref', 'poly_sigma', 'R0ref', & - & 'Web', 'Ca', 'Re_inv', 'sigR', 'sigV', 'rhoRV' ] + & 'Web', 'Ca', 'Re_inv', 'sigR', 'sigV', 'rhoRV', 'palpha_eps', & + & 'ptgalpha_eps' ] call MPI_BCAST(${VAR}$, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) #:endfor @@ -81,7 +82,7 @@ contains #:for VAR in [ 'x_centroid', 'y_centroid', 'z_centroid', & & 'length_x', 'length_y', 'length_z', 'radius', 'epsilon', & & 'beta', 'smooth_coeff', 'rho', 'p0', 'm0', 'r0', 'v0', & - & 'pres', 'gamma', 'pi_inf', 'hcid' ] + & 'pres', 'gamma', 'pi_inf', 'hcid', 'cv', 'qv', 'qvp' ] call MPI_BCAST(patch_icpp(i)%${VAR}$, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) #:endfor @@ -98,7 +99,7 @@ contains ! Fluids physical parameters do i = 1, num_fluids_max #:for VAR in [ 'gamma','pi_inf','mul0','ss','pv','gamma_v','M_v', & - & 'mu_v','k_v', 'G' ] + & 'mu_v','k_v', 'G', 'qv' ] call MPI_BCAST(fluid_pp(i)%${VAR}$, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) #:endfor end do diff --git a/src/pre_process/m_start_up.fpp b/src/pre_process/m_start_up.fpp index 4df15c5e5..452dfb994 100644 --- a/src/pre_process/m_start_up.fpp +++ b/src/pre_process/m_start_up.fpp @@ -31,6 +31,8 @@ module m_start_up use m_assign_variables + use m_phase_change !< Phase-change module + use m_helper #ifdef MFC_MPI @@ -124,7 +126,8 @@ contains polytropic, thermal, Ca, Web, Re_inv, & polydisperse, poly_sigma, qbmm, & sigR, sigV, dist_type, rhoRV, R0_type, & - file_per_process + file_per_process, relax, relax_model, & + palpha_eps, ptgalpha_eps ! Inquiring the status of the pre_process.inp file file_loc = 'pre_process.inp' @@ -703,6 +706,7 @@ contains call s_initialize_grid_module() call s_initialize_initial_condition_module() call s_initialize_assign_variables_module() + if (relax) call s_initialize_phasechange_module() ! Associate pointers for serial or parallel I/O if (parallel_io .neqv. .true.) then @@ -763,6 +767,15 @@ contains call s_generate_initial_condition() + if (relax) then + if (proc_rank == 0) then + print *, 'initial condition might have been altered due to enforcement of & + pTg-equilirium (relax = "T" activated)' + end if + + call s_relaxation_solver(q_cons_vf) + end if + call s_write_data_files(q_cons_vf) call cpu_time(finish) @@ -840,6 +853,8 @@ contains call s_finalize_data_output_module() call s_finalize_global_parameters_module() call s_finalize_assign_variables_module() + if (relax) call s_finalize_relaxation_solver_module() + ! Finalization of the MPI environment call s_mpi_finalize() diff --git a/src/simulation/m_cbc.fpp b/src/simulation/m_cbc.fpp index bbe2191b4..de2ba373b 100644 --- a/src/simulation/m_cbc.fpp +++ b/src/simulation/m_cbc.fpp @@ -597,6 +597,7 @@ contains real(kind(0d0)), dimension(num_fluids) :: dadv_dt real(kind(0d0)) :: dgamma_dt real(kind(0d0)) :: dpi_inf_dt + real(kind(0d0)) :: dqv_dt real(kind(0d0)), dimension(contxe) :: alpha_rho, dalpha_rho_ds, mf real(kind(0d0)), dimension(2) :: Re_cbc real(kind(0d0)), dimension(num_dims) :: vel, dvel_ds @@ -610,6 +611,7 @@ contains real(kind(0d0)) :: H !< Cell averaged enthalpy real(kind(0d0)) :: gamma !< Cell averaged specific heat ratio real(kind(0d0)) :: pi_inf !< Cell averaged liquid stiffness + real(kind(0d0)) :: qv !< Cell averaged fluid reference energy real(kind(0d0)) :: c real(kind(0d0)) :: vel_K_sum, vel_dv_dt_sum @@ -749,10 +751,10 @@ contains end do if (bubbles) then - call s_convert_species_to_mixture_variables_bubbles_acc(rho, gamma, pi_inf, adv, alpha_rho, Re_cbc, 0, k, r) + call s_convert_species_to_mixture_variables_bubbles_acc(rho, gamma, pi_inf, qv, adv, alpha_rho, Re_cbc, 0, k, r) else - call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, adv, alpha_rho, Re_cbc, 0, k, r) + call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, adv, alpha_rho, Re_cbc, 0, k, r) end if !$acc loop seq @@ -878,7 +880,7 @@ contains end do end if - drho_dt = 0d0; dgamma_dt = 0d0; dpi_inf_dt = 0d0 + drho_dt = 0d0; dgamma_dt = 0d0; dpi_inf_dt = 0d0; dqv_dt = 0d0 if (model_eqns == 1) then drho_dt = dalpha_rho_dt(1) @@ -890,6 +892,7 @@ contains drho_dt = drho_dt + dalpha_rho_dt(i) dgamma_dt = dgamma_dt + dadv_dt(i)*gammas(i) dpi_inf_dt = dpi_inf_dt + dadv_dt(i)*pi_infs(i) + dqv_dt = dqv_dt + dalpha_rho_dt(i)*qvs(i) end do end if ! ============================================================ @@ -912,6 +915,7 @@ contains + ds(0)*(pres*dgamma_dt & + gamma*dpres_dt & + dpi_inf_dt & + + dqv_dt & + rho*vel_dv_dt_sum & + 5d-1*drho_dt*vel_K_sum) diff --git a/src/simulation/m_checker.fpp b/src/simulation/m_checker.fpp index 4ea6a2954..de98dd60c 100644 --- a/src/simulation/m_checker.fpp +++ b/src/simulation/m_checker.fpp @@ -135,6 +135,25 @@ contains end if end if + ! phase change checkers. + if ( relax ) then + if ( model_eqns /= 3 ) then + call s_mpi_abort( 'phase change requires model_eqns = 3. ' // & + 'Exiting ...') + elseif ( ( relax_model .lt. 0 ) .or. ( relax_model .gt. 6 ) ) then + call s_mpi_abort( 'relax_model should be in between 0 and 6. ' // & + 'Exiting ...' ) + elseif ( ( palpha_eps .le. 0d0 ) .or. ( palpha_eps .ge. 1d0 ) .or. & + ( ptgalpha_eps .le. 0d0 ) .or. ( ptgalpha_eps .ge. 1d0 ) ) then + call s_mpi_abort( 'both palpha_eps and ptgalpha_eps must & + be in (0,1). ' // 'Exiting ...') + end if + elseif ( ( relax_model /= dflt_int ) .or. ( palpha_eps /= dflt_real ) & + .or. ( ptgalpha_eps /= dflt_real ) ) then + call s_mpi_abort( 'relax is not set as true, but other phase change parameters have & + been modified. Either activate phase change or set the values to default. ' // 'Exiting ...') + end if + if (num_fluids /= dflt_int & .and. & (num_fluids < 1 .or. num_fluids > num_fluids)) then @@ -339,6 +358,10 @@ contains 'of values of num_fluids '// & 'and fluid_pp('//trim(iStr)//')%'// & 'pi_inf. Exiting ...') + elseif ( fluid_pp(i)%cv < 0d0 ) then + call s_mpi_abort('Unsupported value of '// & + 'fluid_pp('//trim(iStr)//')%'// & + 'cv. Make sure cv is positive. Exiting ...') end if do j = 1, 2 diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index 1d594f847..71774e5aa 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -227,6 +227,7 @@ contains real(kind(0d0)), dimension(num_fluids) :: alpha !< Cell-avg. volume fraction real(kind(0d0)) :: gamma !< Cell-avg. sp. heat ratio real(kind(0d0)) :: pi_inf !< Cell-avg. liquid stiffness function + real(kind(0d0)) :: qv !< Cell-avg. fluid reference energy real(kind(0d0)) :: c !< Cell-avg. sound speed real(kind(0d0)) :: E !< Cell-avg. energy real(kind(0d0)) :: H !< Cell-avg. enthalpy @@ -257,9 +258,9 @@ contains end do if (bubbles) then - call s_convert_species_to_mixture_variables_bubbles_acc(rho, gamma, pi_inf, alpha, alpha_rho, Re, j, k, l) + call s_convert_species_to_mixture_variables_bubbles_acc(rho, gamma, pi_inf, qv, alpha, alpha_rho, Re, j, k, l) else - call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, alpha, alpha_rho, Re, j, k, l) + call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, alpha, alpha_rho, Re, j, k, l) end if do i = 1, num_dims @@ -273,7 +274,7 @@ contains pres = q_prim_vf(E_idx)%sf(j, k, l) - E = gamma*pres + pi_inf + 5d-1*rho*vel_sum + E = gamma*pres + pi_inf + 5d-1*rho*vel_sum + qv H = (E + pres)/rho @@ -455,7 +456,7 @@ contains real(kind(0d0)), dimension(nb) :: nRtmp !< Temporary bubble concentration real(kind(0d0)) :: nbub, nR3, vftmp !< Temporary bubble number density - real(kind(0d0)) :: gamma, lit_gamma, pi_inf !< Temporary EOS params + real(kind(0d0)) :: gamma, lit_gamma, pi_inf, qv !< Temporary EOS params real(kind(0d0)) :: rho !< Temporary density real(kind(0d0)), dimension(2) :: Re !< Temporary Reynolds number real(kind(0d0)) :: E_e !< Temp. elastic energy contribution @@ -546,6 +547,7 @@ contains gamma = fluid_pp(1)%gamma lit_gamma = 1d0/fluid_pp(1)%gamma + 1d0 pi_inf = fluid_pp(1)%pi_inf + qv = fluid_pp(1)%qv if (precision == 1) then FMT = "(2F30.3)" @@ -992,6 +994,7 @@ contains real(kind(0d0)), dimension(num_fluids) :: alpha real(kind(0d0)) :: gamma real(kind(0d0)) :: pi_inf + real(kind(0d0)) :: qv real(kind(0d0)) :: c real(kind(0d0)) :: M00, M10, M01, M20, M11, M02 real(kind(0d0)) :: varR, varV @@ -1039,6 +1042,7 @@ contains pres = 0d0 gamma = 0d0 pi_inf = 0d0 + qv = 0d0 c = 0d0 accel = 0d0 nR = 0d0; R = 0d0 @@ -1072,11 +1076,11 @@ contains ! Computing/Sharing necessary state variables if(hypoelasticity) then call s_convert_to_mixture_variables(q_cons_vf, j - 2, k, l, & - rho, gamma, pi_inf, & + rho, gamma, pi_inf, qv, & Re, G, fluid_pp(:)%G) else call s_convert_to_mixture_variables(q_cons_vf, j - 2, k, l, & - rho, gamma, pi_inf) + rho, gamma, pi_inf, qv) end if do s = 1, num_dims vel(s) = q_cons_vf(cont_idx%end + s)%sf(j - 2, k, l)/rho @@ -1088,7 +1092,7 @@ contains q_cons_vf(alf_idx)%sf(j - 2, k, l), & 0.5d0*(q_cons_vf(2)%sf(j - 2, k, l)**2.d0)/ & q_cons_vf(1)%sf(j - 2, k, l), & - pi_inf, gamma, rho, pres, & + pi_inf, gamma, rho, qv, pres, & q_cons_vf(stress_idx%beg)%sf(j - 2, k, l), & q_cons_vf(mom_idx%beg)%sf(j - 2, k, l), G) else @@ -1097,7 +1101,7 @@ contains q_cons_vf(alf_idx)%sf(j - 2, k, l), & 0.5d0*(q_cons_vf(2)%sf(j - 2, k, l)**2.d0)/ & q_cons_vf(1)%sf(j - 2, k, l), & - pi_inf, gamma, rho, pres) + pi_inf, gamma, rho, qv, pres) end if if (model_eqns == 4) then @@ -1176,7 +1180,7 @@ contains ! Computing/Sharing necessary state variables call s_convert_to_mixture_variables(q_cons_vf, j - 2, k - 2, l, & - rho, gamma, pi_inf, & + rho, gamma, pi_inf, qv, & Re, G, fluid_pp(:)%G) do s = 1, num_dims vel(s) = q_cons_vf(cont_idx%end + s)%sf(j - 2, k - 2, l)/rho @@ -1187,7 +1191,7 @@ contains q_cons_vf(alf_idx)%sf(j - 2, k - 2, l), & 0.5d0*(q_cons_vf(2)%sf(j - 2, k - 2, l)**2.d0)/ & q_cons_vf(1)%sf(j - 2, k - 2, l), & - pi_inf, gamma, pres, rho, & + pi_inf, gamma, rho, qv, pres, & q_cons_vf(stress_idx%beg)%sf(j - 2, k - 2, l), & q_cons_vf(mom_idx%beg)%sf(j - 2, k - 2, l), G) @@ -1250,14 +1254,14 @@ contains ! Computing/Sharing necessary state variables call s_convert_to_mixture_variables(q_cons_vf, j - 2, k - 2, l - 2, & - rho, gamma, pi_inf, & + rho, gamma, pi_inf, qv, & Re, G, fluid_pp(:)%G) do s = 1, num_dims vel(s) = q_cons_vf(cont_idx%end + s)%sf(j - 2, k - 2, l - 2)/rho end do call s_compute_pressure(q_cons_vf(E_idx)%sf(j - 2, k - 2, l - 2), & - 0d0, 0.5d0*rho*dot_product(vel, vel), pi_inf, gamma, rho, pres) + 0d0, 0.5d0*rho*dot_product(vel, vel), pi_inf, gamma, rho, qv, pres) ! Compute mixture sound speed call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, & @@ -1270,7 +1274,7 @@ contains end if if (num_procs > 1) then - #:for VAR in ['rho','pres','gamma','pi_inf','c','accel'] + #:for VAR in ['rho','pres','gamma','pi_inf','qv','c','accel'] tmp = ${VAR}$ call s_mpi_allreduce_sum(tmp, ${VAR}$) #:endfor @@ -1408,7 +1412,7 @@ contains end if else write (i + 30, '(6X,F12.6,F24.8,F24.8,F24.8,F24.8,'// & - 'F24.8,F24.8,F24.8,F24.8,'// & + 'F24.8,F24.8,F24.8,F24.8,F24.8,'// & 'F24.8)') & nondim_time, & rho, & @@ -1418,6 +1422,7 @@ contains pres, & gamma, & pi_inf, & + qv, & c, & accel end if @@ -1440,11 +1445,12 @@ contains pres = 0d0 gamma = 0d0 pi_inf = 0d0 + qv = 0d0 if ((integral(i)%xmin <= x_cb(j)) .and. (integral(i)%xmax >= x_cb(j))) then npts = npts + 1 call s_convert_to_mixture_variables(q_cons_vf, j, k, l, & - rho, gamma, pi_inf, Re) + rho, gamma, pi_inf, qv, Re) do s = 1, num_dims vel(s) = q_cons_vf(cont_idx%end + s)%sf(j, k, l)/rho end do @@ -1453,7 +1459,7 @@ contains (q_cons_vf(E_idx)%sf(j, k, l) - & 0.5d0*(q_cons_vf(mom_idx%beg)%sf(j, k, l)**2.d0)/rho)/ & (1.d0 - q_cons_vf(alf_idx)%sf(j, k, l)) - & - pi_inf & + pi_inf - qv & )/gamma int_pres = int_pres + (pres - 1.d0)**2.d0 end if @@ -1511,11 +1517,12 @@ contains pres = 0d0 gamma = 0d0 pi_inf = 0d0 + qv = 0d0 if (trigger) then npts = npts + 1 call s_convert_to_mixture_variables(q_cons_vf, j, k, l, & - rho, gamma, pi_inf, Re) + rho, gamma, pi_inf, qv, Re) do s = 1, num_dims vel(s) = q_cons_vf(cont_idx%end + s)%sf(j, k, l)/rho end do @@ -1524,7 +1531,7 @@ contains (q_cons_vf(E_idx)%sf(j, k, l) - & 0.5d0*(q_cons_vf(mom_idx%beg)%sf(j, k, l)**2.d0)/rho)/ & (1.d0 - q_cons_vf(alf_idx)%sf(j, k, l)) - & - pi_inf & + pi_inf - qv & )/gamma int_pres = int_pres + abs(pres - 1.d0) max_pres = max(max_pres, abs(pres - 1.d0)) diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index 439533ca3..494c209cf 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -114,6 +114,10 @@ module m_global_parameters integer :: riemann_solver !< Riemann solver algorithm integer :: wave_speeds !< Wave speeds estimation method integer :: avg_state !< Average state evaluation method + logical :: relax !< activate phase change + integer :: relax_model !< Relaxation model + real(kind(0d0)) :: palpha_eps !< trigger parameter for the p relaxation procedure, phase change model + real(kind(0d0)) :: ptgalpha_eps !< trigger parameter for the pTg relaxation procedure, phase change model logical :: alt_soundspeed !< Alternate mixture sound speed logical :: null_weights !< Null undesired WENO weights logical :: mixture_err !< Mixture properties correction @@ -126,7 +130,7 @@ module m_global_parameters !$acc declare create(num_dims, weno_polyn, weno_order) #:endif -!$acc declare create(mpp_lim, num_fluids, model_eqns, mixture_err, alt_soundspeed, avg_state, mapped_weno, mp_weno, weno_eps, hypoelasticity) +!$acc declare create(mpp_lim, num_fluids, model_eqns, mixture_err, alt_soundspeed, avg_state, mapped_weno, mp_weno, weno_eps, hypoelasticity, relax, palpha_eps,ptgalpha_eps) !> @name Boundary conditions (BC) in the x-, y- and z-directions, respectively !> @{ @@ -317,8 +321,8 @@ module m_global_parameters integer :: strxb, strxe !$acc declare create(momxb, momxe, advxb, advxe, contxb, contxe, intxb, intxe, bubxb, bubxe, strxb, strxe) - real(kind(0d0)), allocatable, dimension(:) :: gammas, pi_infs - !$acc declare create(gammas, pi_infs) + real(kind(0d0)), allocatable, dimension(:) :: gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps + !$acc declare create(gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps) real(kind(0d0)) :: mytime !< Current simulation time @@ -376,6 +380,10 @@ contains parallel_io = .false. file_per_process = .false. precision = 2 + relax = .false. + relax_model = dflt_int + palpha_eps = dflt_real + ptgalpha_eps = dflt_real hypoelasticity = .false. weno_flat = .true. riemann_flat = .true. @@ -393,6 +401,9 @@ contains do i = 1, num_fluids_max fluid_pp(i)%gamma = dflt_real fluid_pp(i)%pi_inf = dflt_real + fluid_pp(i)%cv = 0d0 + fluid_pp(i)%qv = 0d0 + fluid_pp(i)%qvp = 0d0 fluid_pp(i)%Re(:) = dflt_real fluid_pp(i)%mul0 = dflt_real fluid_pp(i)%ss = dflt_real diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp index 73ece6f0a..5b7b89caf 100644 --- a/src/simulation/m_mpi_proxy.fpp +++ b/src/simulation/m_mpi_proxy.fpp @@ -127,7 +127,7 @@ contains & 'wave_speeds', 'avg_state', 'precision', 'bc_x%beg', 'bc_x%end', & & 'bc_y%beg', 'bc_y%end', 'bc_z%beg', 'bc_z%end', 'fd_order', & & 'num_probes', 'num_integrals', 'bubble_model', 'thermal', & - & 'R0_type', 'num_mono' ] + & 'R0_type', 'num_mono', 'relax_model'] call MPI_BCAST(${VAR}$, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) #:endfor @@ -136,12 +136,12 @@ contains & 'weno_Re_flux', 'alt_soundspeed', 'null_weights', 'mixture_err', & & 'parallel_io', 'hypoelasticity', 'bubbles', 'polytropic', & & 'polydisperse', 'qbmm', 'monopole', 'probe_wrt', 'integral_wrt', & - & 'prim_vars_wrt', 'weno_avg', 'file_per_process'] + & 'prim_vars_wrt', 'weno_avg', 'file_per_process', 'relax' ] call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) #:endfor #:for VAR in [ 'dt','weno_eps','pref','rhoref','R0ref','Web','Ca', & - & 'Re_inv','poly_sigma' ] + & 'Re_inv','poly_sigma', 'palpha_eps', 'ptgalpha_eps' ] call MPI_BCAST(${VAR}$, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) #:endfor @@ -152,7 +152,7 @@ contains do i = 1, num_fluids_max #:for VAR in [ 'gamma','pi_inf','mul0','ss','pv','gamma_v','M_v', & - & 'mu_v','k_v','G' ] + & 'mu_v','k_v','G', 'cv', 'qv', 'qvp' ] call MPI_BCAST(fluid_pp(i)%${VAR}$, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) #:endfor call MPI_BCAST(fluid_pp(i)%Re(1), 2, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index 1e36dce6a..4319f0c7a 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -257,6 +257,7 @@ contains real(kind(0d0)) :: Y_L, Y_R real(kind(0d0)) :: gamma_L, gamma_R real(kind(0d0)) :: pi_inf_L, pi_inf_R + real(kind(0d0)) :: qv_L, qv_R real(kind(0d0)) :: c_L, c_R real(kind(0d0)), dimension(6) :: tau_e_L, tau_e_R real(kind(0d0)) :: G_L, G_R @@ -341,10 +342,12 @@ contains rho_L = 0d0 gamma_L = 0d0 pi_inf_L = 0d0 + qv_L = 0d0 rho_R = 0d0 gamma_R = 0d0 pi_inf_R = 0d0 + qv_R = 0d0 alpha_L_sum = 0d0 alpha_R_sum = 0d0 @@ -374,10 +377,12 @@ contains rho_L = rho_L + alpha_rho_L(i) gamma_L = gamma_L + alpha_L(i)*gammas(i) pi_inf_L = pi_inf_L + alpha_L(i)*pi_infs(i) + qv_L = qv_L + alpha_rho_L(i)*qvs(i) rho_R = rho_R + alpha_rho_R(i) gamma_R = gamma_R + alpha_R(i)*gammas(i) pi_inf_R = pi_inf_R + alpha_R(i)*pi_infs(i) + qv_R = qv_R + alpha_rho_R(i)*qvs(i) end do if (any(Re_size > 0)) then @@ -413,8 +418,8 @@ contains end do end if - E_L = gamma_L*pres_L + pi_inf_L + 5d-1*rho_L*vel_L_rms - E_R = gamma_R*pres_R + pi_inf_R + 5d-1*rho_R*vel_R_rms + E_L = gamma_L*pres_L + pi_inf_L + 5d-1*rho_L*vel_L_rms + qv_L + E_R = gamma_R*pres_R + pi_inf_R + 5d-1*rho_R*vel_R_rms + qv_R H_L = (E_L + pres_L)/rho_L H_R = (E_R + pres_R)/rho_R @@ -793,6 +798,7 @@ contains real(kind(0d0)) :: Y_L, Y_R real(kind(0d0)) :: gamma_L, gamma_R real(kind(0d0)) :: pi_inf_L, pi_inf_R + real(kind(0d0)) :: qv_L, qv_R real(kind(0d0)) :: c_L, c_R real(kind(0d0)), dimension(2) :: Re_L, Re_R @@ -879,10 +885,12 @@ contains rho_L = 0d0 gamma_L = 0d0 pi_inf_L = 0d0 + qv_L = 0d0 rho_R = 0d0 gamma_R = 0d0 pi_inf_R = 0d0 + qv_R = 0d0 alpha_L_sum = 0d0 alpha_R_sum = 0d0 @@ -918,10 +926,12 @@ contains rho_L = rho_L + qL_prim_rs${XYZ}$_vf(j, k, l, i) gamma_L = gamma_L + qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i)*gammas(i) pi_inf_L = pi_inf_L + qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i)*pi_infs(i) + qv_L = qv_L + qL_prim_rs${XYZ}$_vf(j, k, l, i)*qvs(i) rho_R = rho_R + qR_prim_rs${XYZ}$_vf(j + 1, k, l, i) gamma_R = gamma_R + qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + i)*gammas(i) pi_inf_R = pi_inf_R + qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + i)*pi_infs(i) + qv_R = qv_R + qR_prim_rs${XYZ}$_vf(j + 1, k, l, i)*qvs(i) alpha_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, advxb + i - 1) alpha_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, advxb + i - 1) @@ -960,9 +970,9 @@ contains end do end if - E_L = gamma_L*pres_L + pi_inf_L + 5d-1*rho_L*vel_L_rms + E_L = gamma_L*pres_L + pi_inf_L + 5d-1*rho_L*vel_L_rms + qv_L - E_R = gamma_R*pres_R + pi_inf_R + 5d-1*rho_R*vel_R_rms + E_R = gamma_R*pres_R + pi_inf_R + 5d-1*rho_R*vel_R_rms + qv_R H_L = (E_L + pres_L)/rho_L H_R = (E_R + pres_R)/rho_R @@ -1031,8 +1041,10 @@ contains qL_prim_rs${XYZ}$_vf(j, k, l, i + contxb - 1)*vel_L(dir_idx(1)) flux_rs${XYZ}$_vf(j, k, l, i + intxb - 1) = & - qL_prim_rs${XYZ}$_vf(j, k, l, i + advxb - 1)* & - (gammas(i)*pres_L + pi_infs(i))*vel_L(dir_idx(1)) + ( qL_prim_rs${XYZ}$_vf(j, k, l, i + advxb - 1) * & + ( gammas(i)*pres_L + pi_infs(i) ) + & + qL_prim_rs${XYZ}$_vf(j, k, l, i + contxb - 1) * & + qvs(i) ) * vel_L(dir_idx(1) ) end do !$acc loop seq do i = 1, num_dims @@ -1058,8 +1070,10 @@ contains qR_prim_rs${XYZ}$_vf(j + 1, k, l, i + contxb - 1)*vel_R(dir_idx(1)) flux_rs${XYZ}$_vf(j, k, l, i + intxb - 1) = & - qR_prim_rs${XYZ}$_vf(j + 1, k, l, i + advxb - 1)* & - (gammas(i)*pres_R + pi_infs(i))*vel_R(dir_idx(1)) + ( qR_prim_rs${XYZ}$_vf(j + 1, k, l, i + advxb - 1) * & + ( gammas(i)*pres_R + pi_infs(i) ) + & + qR_prim_rs${XYZ}$_vf(j + 1, k, l, i + contxb - 1) * & + qvs(i) ) * vel_R(dir_idx(1)) end do !$acc loop seq do i = 1, num_dims @@ -1091,8 +1105,10 @@ contains qL_prim_rs${XYZ}$_vf(j, k, l, i + contxb - 1)*xi_L*s_S flux_rs${XYZ}$_vf(j, k, l, i + intxb - 1) = & - qL_prim_rs${XYZ}$_vf(j, k, l, i + advxb - 1)* & - (gammas(i)*p_K_Star + pi_infs(i))*s_S + ( qL_prim_rs${XYZ}$_vf(j, k, l, i + advxb - 1) * & + ( gammas(i)*p_K_Star + pi_infs(i) ) + & + qL_prim_rs${XYZ}$_vf(j, k, l, i + contxb - 1) * & + qvs(i) ) * s_S end do !$acc loop seq do i = 1, num_dims @@ -1128,8 +1144,10 @@ contains qR_prim_rs${XYZ}$_vf(j + 1, k, l, i + contxb - 1)*xi_R*s_S flux_rs${XYZ}$_vf(j, k, l, i + intxb - 1) = & - qR_prim_rs${XYZ}$_vf(j + 1, k, l, i + advxb - 1)* & - (gammas(i)*p_K_Star + pi_infs(i))*s_S + ( qR_prim_rs${XYZ}$_vf(j + 1, k, l, i + advxb - 1) * & + ( gammas(i)*p_K_Star + pi_infs(i) ) + & + qR_prim_rs${XYZ}$_vf(j + 1, k, l, i + contxb - 1) * & + qvs(i) ) * s_S end do !$acc loop seq do i = 1, num_dims @@ -1210,26 +1228,30 @@ contains rho_L = 0d0 gamma_L = 0d0 pi_inf_L = 0d0 + qv_L = 0d0 !$acc loop seq do i = 1, num_fluids rho_L = rho_L + alpha_rho_L(i) gamma_L = gamma_L + alpha_L(i)*gammas(i) pi_inf_L = pi_inf_L + alpha_L(i)*pi_infs(i) + qv_L = qv_L + alpha_rho_L(i)*qvs(i) end do rho_R = 0d0 gamma_R = 0d0 pi_inf_R = 0d0 + qv_R = 0d0 !$acc loop seq do i = 1, num_fluids rho_R = rho_R + alpha_rho_R(i) gamma_R = gamma_R + alpha_R(i)*gammas(i) pi_inf_R = pi_inf_R + alpha_R(i)*pi_infs(i) + qv_R = qv_R + alpha_rho_R(i)*qvs(i) end do - E_L = gamma_L*pres_L + pi_inf_L + 5d-1*rho_L*vel_L_rms + E_L = gamma_L*pres_L + pi_inf_L + 5d-1*rho_L*vel_L_rms + qv_L - E_R = gamma_R*pres_R + pi_inf_R + 5d-1*rho_R*vel_R_rms + E_R = gamma_R*pres_R + pi_inf_R + 5d-1*rho_R*vel_R_rms + qv_R H_L = (E_L + pres_L)/rho_L H_R = (E_R + pres_R)/rho_R @@ -1447,6 +1469,7 @@ contains rho_L = 0d0 gamma_L = 0d0 pi_inf_L = 0d0 + qv_L = 0d0 if (mpp_lim .and. (num_fluids > 2)) then !$acc loop seq @@ -1454,6 +1477,7 @@ contains rho_L = rho_L + qL_prim_rs${XYZ}$_vf(j, k, l, i) gamma_L = gamma_L + qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i)*gammas(i) pi_inf_L = pi_inf_L + qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i)*pi_infs(i) + qv_L = qv_L + qL_prim_rs${XYZ}$_vf(j, k, l, i)*qvs(i) end do else if (num_fluids > 2) then !$acc loop seq @@ -1461,16 +1485,19 @@ contains rho_L = rho_L + qL_prim_rs${XYZ}$_vf(j, k, l, i) gamma_L = gamma_L + qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i)*gammas(i) pi_inf_L = pi_inf_L + qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i)*pi_infs(i) + qv_L = qv_L + qL_prim_rs${XYZ}$_vf(j, k, l, i)*qvs(i) end do else rho_L = qL_prim_rs${XYZ}$_vf(j, k, l, 1) gamma_L = gammas(1) pi_inf_L = pi_infs(1) + qv_L = qvs(1) end if rho_R = 0d0 gamma_R = 0d0 pi_inf_R = 0d0 + qv_R = 0d0 if (mpp_lim .and. (num_fluids > 2)) then !$acc loop seq @@ -1478,6 +1505,7 @@ contains rho_R = rho_R + qR_prim_rs${XYZ}$_vf(j + 1, k, l, i) gamma_R = gamma_R + qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + i)*gammas(i) pi_inf_R = pi_inf_R + qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + i)*pi_infs(i) + qv_R = qv_R + qR_prim_rs${XYZ}$_vf(j + 1, k, l, i)*qvs(i) end do else if (num_fluids > 2) then !$acc loop seq @@ -1485,11 +1513,13 @@ contains rho_R = rho_R + qR_prim_rs${XYZ}$_vf(j + 1, k, l, i) gamma_R = gamma_R + qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + i)*gammas(i) pi_inf_R = pi_inf_R + qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + i)*pi_infs(i) + qv_R = qv_R + qR_prim_rs${XYZ}$_vf(j + 1, k, l, i)*qvs(i) end do else rho_R = qR_prim_rs${XYZ}$_vf(j + 1, k, l, 1) gamma_R = gammas(1) pi_inf_R = pi_infs(1) + qv_R = qvs(1) end if if (any(Re_size > 0)) then @@ -1884,10 +1914,12 @@ contains rho_L = 0d0 gamma_L = 0d0 pi_inf_L = 0d0 + qv_L = 0d0 rho_R = 0d0 gamma_R = 0d0 pi_inf_R = 0d0 + qv_R = 0d0 alpha_L_sum = 0d0 alpha_R_sum = 0d0 @@ -1923,10 +1955,12 @@ contains rho_L = rho_L + qL_prim_rs${XYZ}$_vf(j, k, l, i) gamma_L = gamma_L + qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i)*gammas(i) pi_inf_L = pi_inf_L + qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + i)*pi_infs(i) + qv_L = qv_L + qL_prim_rs${XYZ}$_vf(j, k, l, i)*qvs(i) rho_R = rho_R + qR_prim_rs${XYZ}$_vf(j + 1, k, l, i) gamma_R = gamma_R + qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + i)*gammas(i) pi_inf_R = pi_inf_R + qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + i)*pi_infs(i) + qv_R = qv_R + qR_prim_rs${XYZ}$_vf(j + 1, k, l, i)*qvs(i) end do if (any(Re_size > 0)) then @@ -1962,9 +1996,9 @@ contains end do end if - E_L = gamma_L*pres_L + pi_inf_L + 5d-1*rho_L*vel_L_rms + E_L = gamma_L*pres_L + pi_inf_L + 5d-1*rho_L*vel_L_rms + qv_L - E_R = gamma_R*pres_R + pi_inf_R + 5d-1*rho_R*vel_R_rms + E_R = gamma_R*pres_R + pi_inf_R + 5d-1*rho_R*vel_R_rms + qv_R H_L = (E_L + pres_L)/rho_L H_R = (E_R + pres_R)/rho_R diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 82b35173b..eac47bcd3 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -49,6 +49,8 @@ module m_start_up use m_hypoelastic + use m_phase_change !< Phase-change module + use m_viscous use m_bubbles @@ -143,7 +145,8 @@ contains polytropic, thermal, & integral, integral_wrt, num_integrals, & polydisperse, poly_sigma, qbmm, & - R0_type, file_per_process + R0_type, file_per_process, relax, relax_model, & + palpha_eps, ptgalpha_eps ! Checking that an input file has been provided by the user. If it ! has, then the input file is read in, otherwise, simulation exits. @@ -832,6 +835,7 @@ contains real(kind(0d0)) :: dyn_pres real(kind(0d0)) :: gamma real(kind(0d0)) :: pi_inf + real(kind(0d0)) :: qv real(kind(0d0)), dimension(2) :: Re real(kind(0d0)) :: pres @@ -841,7 +845,7 @@ contains do k = 0, n do l = 0, p - call s_convert_to_mixture_variables(v_vf, j, k, l, rho, gamma, pi_inf, Re) + call s_convert_to_mixture_variables(v_vf, j, k, l, rho, gamma, pi_inf, qv, Re) dyn_pres = 0d0 do i = mom_idx%beg, mom_idx%end @@ -850,11 +854,12 @@ contains end do call s_compute_pressure(v_vf(E_idx)%sf(j, k, l), 0d0, & - dyn_pres, pi_inf, gamma, rho, pres) + dyn_pres, pi_inf, gamma, rho, qv, pres) do i = 1, num_fluids v_vf(i + internalEnergies_idx%beg - 1)%sf(j, k, l) = v_vf(i + adv_idx%beg - 1)%sf(j, k, l)* & - (fluid_pp(i)%gamma*pres + fluid_pp(i)%pi_inf) + (fluid_pp(i)%gamma*pres + fluid_pp(i)%pi_inf) & + + v_vf(i + cont_idx%beg - 1)%sf(j, k, l)*fluid_pp(i)%qv end do end do @@ -905,6 +910,8 @@ contains call s_3rd_order_tvd_rk(t_step, time_avg) end if + if ( relax ) call s_relaxation_solver(q_cons_ts(1)%vf) + ! Time-stepping loop controls if ((mytime + dt) >= finaltime) dt = finaltime - mytime t_step = t_step + 1 @@ -1056,6 +1063,7 @@ contains #endif if (hypoelasticity) call s_initialize_hypoelastic_module() + if (relax) call s_initialize_phasechange_module() call s_initialize_data_output_module() call s_initialize_derived_variables_module() call s_initialize_time_steppers_module() @@ -1169,6 +1177,10 @@ contains !$acc update device(nb, R0ref, Ca, Web, Re_inv, weight, R0, V0, bubbles, polytropic, polydisperse, qbmm, R0_type, ptil, bubble_model, thermal, poly_sigma) !$acc update device(R_n, R_v, phi_vn, phi_nv, Pe_c, Tw, pv, M_n, M_v, k_n, k_v, pb0, mass_n0, mass_v0, Pe_T, Re_trans_T, Re_trans_c, Im_trans_T, Im_trans_c, omegaN , mul0, ss, gamma_v, mu_v, gamma_m, gamma_n, mu_n, gam) !$acc update device(monopole, num_mono) + !$acc update device(relax) + if(relax) then + !$acc update device(palpha_eps, ptgalpha_eps) + end if end subroutine s_initialize_gpu_vars @@ -1188,6 +1200,7 @@ contains if (grid_geometry == 3) call s_finalize_fftw_module call s_finalize_mpi_proxy_module() call s_finalize_global_parameters_module() + if (relax) call s_finalize_relaxation_solver_module() if (any(Re_size > 0)) then call s_finalize_viscous_module() diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index e782d3aef..679ed254c 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -395,7 +395,9 @@ contains if (grid_geometry == 3) call s_apply_fourier_filter(q_cons_ts(2)%vf) - if (model_eqns == 3) call s_pressure_relaxation_procedure(q_cons_ts(2)%vf) + if ( model_eqns == 3 .and. (.not.relax ) ) then + call s_pressure_relaxation_procedure(q_cons_ts(2)%vf) + end if ! ================================================================== ! Stage 2 of 2 ===================================================== @@ -454,7 +456,9 @@ contains if (grid_geometry == 3) call s_apply_fourier_filter(q_cons_ts(1)%vf) - if (model_eqns == 3) call s_pressure_relaxation_procedure(q_cons_ts(1)%vf) + if ( model_eqns == 3 .and. (.not.relax ) ) then + call s_pressure_relaxation_procedure(q_cons_ts(1)%vf) + end if call nvtxEndRange @@ -548,7 +552,9 @@ contains if (grid_geometry == 3) call s_apply_fourier_filter(q_cons_ts(2)%vf) - if (model_eqns == 3) call s_pressure_relaxation_procedure(q_cons_ts(2)%vf) + if ( model_eqns == 3 .and. (.not.relax ) ) then + call s_pressure_relaxation_procedure(q_cons_ts(2)%vf) + end if ! ================================================================== @@ -607,7 +613,9 @@ contains end if if (grid_geometry == 3) call s_apply_fourier_filter(q_cons_ts(2)%vf) - if (model_eqns == 3) call s_pressure_relaxation_procedure(q_cons_ts(2)%vf) + if ( model_eqns == 3 .and. (.not.relax ) ) then + call s_pressure_relaxation_procedure(q_cons_ts(2)%vf) + end if ! ================================================================== @@ -666,7 +674,9 @@ contains if (grid_geometry == 3) call s_apply_fourier_filter(q_cons_ts(1)%vf) - if (model_eqns == 3) call s_pressure_relaxation_procedure(q_cons_ts(1)%vf) + if ( model_eqns == 3 .and. (.not.relax ) ) then + call s_pressure_relaxation_procedure(q_cons_ts(1)%vf) + end if call nvtxEndRange