[FLASH-BUGS] "minor" glitch in flatten.F90
Tomasz Plewa
tomek at flash.uchicago.edu
Wed Jul 16 09:42:18 CDT 2003
Markus -
No such bias should be even formally present in the code. As a matter
of fact, some implementation of PPM do use modification similar to
yours. I simply do not think that C&W were considering existence of
possible left-right biases in their implementation. And, of course,
Your modification is reasonable and will produce unbiased result as
long as calculation of scrch3 is symmetric with respect to "i" index.
(This is actually the case.)
Thanks for pointing it out to us -
Tomek
--
On Wed, Jul 16, 2003 at 10:47:52AM +0100, Markus Gross wrote:
> Hi!
>
> I was scratching my head about this for the last couple of days and finally
> found a reason for it.
>
> Situation:
>
> assuming 1d, hydro, ppm, split. you have two blocks with identical data, the
> only difference is that the one block is an exact mirror copy of the other.
> ie:
>
> p_block1(:)= p_block2(16:1:-1)
> u_block1(:)=-u_block2(16:1:-1)
>
> for pressure and velocity.
>
> Expecting the same results for block1 and block2 would be reasonable but
> leads to suprises and a lot of scratching of the head.
>
> The reason for this was in the flattening which has a bias towards one
> direction if the pressure gradient is zero in a small section. There seems to
> be a similar bias in the monotonizing, just didn't get around to investigate
> it yet.
>
> The problem is:
>
> !this is the original
> !!$ if (dp(i) .LT. 0.e0) then
> !!$ scrch2(i) = scrch3(i+1)
> !!$ else
> !!$ scrch2(i) = scrch3(i-1)
> !!$ endif
>
> I replaced it with:
>
> !
> !
> ! this is the at.hoc modification
> if (dp(i) .LT. 0.e0) scrch2(i) = scrch3(i+1)
> if (dp(i) .GT. 0.e0) scrch2(i) = scrch3(i-1)
> if (dp(i) .eq. 0.e0) scrch2(i) = scrch3(i )
>
> I didn't check if it is reasonable yet, because my Colella &
> Woodward is buried here somewere at the moment and I cannot find it ...
>
> I attached an example with hardcoded values for pressure and velocity. This
> small programm runs twice, for the original and the mirror data.
>
> What do you think? Is a special case somewhat of course, but some bc's (as
> mine) or ic's may produce this sort of data and error subsequently.
>
>
> Markus.
>
>
>
>
>
>
> --
> _______________________________________________________________
>
> Markus Gross AMIMechE BEng (Hons.) Mechanical Engineering
>
> Heriot Watt University Edinburgh
> School of Engineering and Physical Sciences
> James Nasmyth Building
> Edinburgh
> EH14 4AS
> UK
>
> Member of IMechE, SPIE, CSME and VDI
> _______________________________________________________________
>
> further contact:
>
> Phone : +44 (0) 131 451 4737
>
> UNiX talk: talk markus at lasersim.mce.hw.ac.uk
>
> _______________________________________________________________
>
> "Plans are a place to begin," Grove said. "They rarely deliver
> you to where you expect. Make your plans knowing you are going
> to throw them away."
>
> _______________________________________________________________
>
>
>
>
> !rho: 0.10000000E-02 0.10000000E-02 0.14988374E-02 0.14988374E-02 0.14988374E-02 0.14988374E-02 0.10000000E-02
> !u: 0.00000000E+00 0.00000000E+00 0.00000000E+00 -0.32206959E+04 0.32206959E+04 0.00000000E+00 0.00000000E+00
> !p: 0.84087411E+06 0.84087411E+06 0.26033973E+07 0.26033973E+07 0.26033973E+07 0.26033973E+07 0.84087411E+06
>
>
>
> ! flaten zone structure in regions where shocks are too thin
>
> implicit none
> integer, parameter :: q=16
>
> integer:: nzn, ishkbn
> real, DIMENSION(q) :: u, p,flatn,flatn1,dp,du,shockd,shockd1
> real :: epsiln, omg1, omg2,smallu,smallp,igodu
>
> integer :: i, nzn5, nzn6, nzn7, nzn8
> real :: shkbrn, utest, dutest, dp2, dptest, ptest
> real, DIMENSION(q) :: scrch1(q), scrch2(q), scrch3(q)
> real :: sig1,sig2,ak1,ak2,wig1,wig2,wig3
> logical :: repeat = .FALSE.
>
> u=0.d0
> u(1) = 0.00000000E+00
> u(2) = 0.00000000E+00
> u(3) = 0.00000000E+00
> u(4) = -0.32206959E+04
> u(5) = 0.32206959E+04
> u(6) = 0.00000000E+00
> u(7) = 0.00000000E+00
>
> p = 0.84087411E+06
> p(1) = 0.84087411E+06
> p(2) = 0.84087411E+06
> p(3) = 0.26033973E+07
> p(4) = 0.26033973E+07
> p(5) = 0.26033973E+07
> p(6) = 0.26033973E+07
> p(7) = 0.84087411E+06
>
> 10 continue
>
> epsiln = 0.33
>
> ishkbn =0.
>
> igodu=0.d0
> smallu = 1.d-10
> smallp =1.d-20
>
> omg1 = 0.75
> omg2 = 10.0
> sig1 = 0.50
> sig2 = 1.00
> ak1 = 2.00
> ak2 = 0.01
>
> wig1 = 2.00
> wig2 = 0.00
> wig3 = 0.3333 - wig2
>
> nzn = 8
>
> nzn5 = nzn + 5
> nzn6 = nzn + 6
> nzn7 = nzn + 7
> nzn8 = nzn + 8
>
> do i = 1, nzn8
> flatn (i) = 0.0e00
> shockd(i) = 0.0e00
> end do
>
> do i = 2, nzn7
> shkbrn = ishkbn
>
> ! compute the w_j parameter in Eq. A.1 in Colella & Woodward. w_j is equal to
> ! 1 if the jth zone is inside a pressure and velocity jump in the sweep direction,
> ! in a manner consistent with a shock -- store it in shockd1
>
> dp(i) = p(i+1) - p(i-1)
> du(i) = u(i+1) - u(i-1)
> scrch1(i) = epsiln * min (p(i+1), p(i-1)) - abs( dp(i) )
> utest = smallu - abs (du(i))
>
> if (utest .LT. 0.e0) then
> dutest = du(i)
> else
> dutest = 0.e0
> endif
>
> if (scrch1(i) .LT. 0.e0) then
> scrch1(i) = 1.e0
> else
> scrch1(i) = 0.e0
> endif
>
> if (du(i) .GE. 0.e0) scrch1(i) = 0.e0
>
> if (dutest .EQ. 0.e0) scrch1(i) = 0.e0
>
> if (shkbrn .NE. 0.e0) then
> shockd1(i) = 0.e0
> else
> shockd1(i) = scrch1(i)
> endif
>
> end do
>
> do i = 3, nzn6
> shockd(i) = max (shockd1(i-1), shockd1(i), shockd1(i+1))
> end do
>
>
> ! compute the dissipative flux, using Eq. A.2 in Colella & Woodward
>
> do i = 3, nzn6
> dp2 = p(i+2) - p(i-2)
>
> if (dp2 .EQ. 0.e0) dp2 = smallp
>
> dptest = dp(i)
> ptest = dp2 - smallp
>
> if (ptest .EQ. 0.e0) dptest = 0.e0
>
> scrch2(i) = dptest / dp2 - omg1
> scrch3(i) = scrch1(i) * max (0.e00, scrch2(i) * omg2)
> end do
>
> do i = 2, nzn7
> !
> !
> ! this is the at.hoc modification
> if (dp(i) .LT. 0.e0) scrch2(i) = scrch3(i+1)
> if (dp(i) .GT. 0.e0) scrch2(i) = scrch3(i-1)
> if (dp(i) .eq. 0.e0) scrch2(i) = scrch3(i )
> !
> !this is the original
> !!$ if (dp(i) .LT. 0.e0) then
> !!$ scrch2(i) = scrch3(i+1)
> !!$ else
> !!$ scrch2(i) = scrch3(i-1)
> !!$ endif
> !
> end do
>
> do i = 1, nzn6
> flatn(i) = max (scrch3(i), scrch2(i))
> flatn(i) = max (0.0e00, min (1.0e00, flatn(i)))
> end do
>
> do i = 1, nzn8
> flatn (i) = flatn(i) * (1.0e00 - igodu) + igodu
> flatn1(i) = 1.0e00 - flatn(i)
> end do
>
> do i=1,16
> write(6,fmt='(I,7D10.3)') i,flatn(i),scrch1(i),scrch2(i),scrch3(i),dp(i),du(i),p(i)
> end do
>
> p(:)=p(16:1:-1)
> u(:)=-u(16:1:-1)
>
> if(.not.repeat) THEN
> repeat = .TRUE.
> write(6,*)
> goto 10
> end if
>
>
> end
--
More information about the flash-bugs
mailing list