Examples » Walker: Integrating the number-fraction beta SDE

This example runs Walker to integrate the number-fraction beta SDE (see DiffEq/NumberFractionBeta.h) using constant coefficients. The number-fraction beta SDE is based on the beta SDE, additionally computing two more stochastic variables that are functions of the beta variables integrated. For more detail on the beta SDE, see https://doi.org/10.1080/14685248.2010.510843.

Control file

title "Test Ray's closure ideas for <v^2> and <v^3>"

walker

  term  25.0    # Max time
  dt    0.002   # Time step size
  npar  10000   # Number of particles
  ttyi  100     # TTY output interval

  rngs
    mkl_r250 end
  end

  numfracbeta   # Select the number-fraction beta SDE
    depvar x    # Dependent variable: X, denoting mole or number fractions
    ncomp 15    # = 3x5 = 5 systems each computing 3 variables:
                #   X - number fraction,
                #   R - density,
                #   V - specific volume
    init zero
    coeff const
    # alpha = Sb/kappa, beta = (1-S)b/kappa
    kappa 2.0  0.76923  0.5  0.15873  0.1 end
    b     0.4  1.0      1.0  1.0    100.0 end
    S     0.5  0.53846  0.5  0.39683  0.5 end
    rng mkl_r250
    rho2 1.0 1.0 1.0 1.0 1.0 end
    # low-A
    rcomma 1.0e-2 1.0e-2 1.0e-2 1.0e-2 1.0e-2 end
    # high-A
    #rcomma 0.9 0.9 0.9 0.9 0.9 end
  end

  statistics    # column numbers in output file
    # <X>, mole fraction means
    <X1>        # 3
    <X2>        # 4
    <X3>        # 5
    <X4>        # 6
    <X5>        # 7
    # <rho>, mean densities
    <X6>        # 8
    <X7>        # 9
    <X8>        # 10
    <X9>        # 11
    <X10>       # 12
    # <V>, mean specific volumes
    <X11>       # 13
    <X12>       # 14
    <X13>       # 15
    <X14>       # 16
    <X15>       # 17
    # <x^2>, mole fraction variances
    <x1x1>      # 18
    <x2x2>      # 19
    <x3x3>      # 20
    <x4x4>      # 21
    <x5x5>      # 22
     # <rho v>, density-specific-volume covariances
    <x6x11>     # 25
    <x7x12>     # 30
    <x8x13>     # 25
    <x9x14>     # 40
    <x10x15>    # 45
    # <rho v^2>
    <x6x11x11>  # 26
    <x7x12x12>  # 31
    <x8x13x13>  # 36
    <x9x14x14>  # 41
    <x10x15x15> # 46
    # <rho^2 v>
    <x6x6x11>   # 23
    <x7x7x12>   # 28
    <x8x8x13>   # 33
    <x9x9x14>   # 38
    <x10x10x15> # 43
    # <rho^2v^2>
    <x6x6x11x11>   # 24
    <x7x7x12x12>   # 29
    <x8x8x13x13>   # 34
    <x9x9x14x14>   # 39
    <x10x10x15x15> # 43
    # <rho v^3>
    <x6x11x11x11>  # 27
    <x7x12x12x12>  # 32
    <x8x13x13x13>  # 37
    <x9x14x14x14>  # 42
    <x10x15x15x15> # 47
    # <v^2>, specific volume variances
    <x11x11>    # 48
    <x12x12>    # 50
    <x13x13>    # 52
    <x14x14>    # 54
    <x15x15>    # 56
    # <v^3>, specific volume third central moments
    <x11x11x11> # 49
    <x12x12x12> # 51
    <x13x13x13> # 53
    <x14x14x14> # 55
    <x15x15x15> # 57
  end

  pdfs
    interval  100
    filetype  txt
    policy    overwrite
    centering elem
    format    scientific
    precision 4
    p1( X1 : 1.0e-2 )
    p2( X2 : 1.0e-2 )
    p3( X3 : 1.0e-2 )
    p4( X4 : 1.0e-2 )
    p5( X5 : 1.0e-2 )
  end
end

Example run on 8 CPUs

./charmrun +p8 Main/walker -v -c ../../tmp/numfracbeta.q -u 0.9

Results

The rationale for these runs is to integrate the system in time, extract the time evolution of various statistics, and then test several different closure hypotheses among the statistics. Example gnuplot commands to plot and test some closure ideas are:

# vim: filetype=gnuplot:

# Nomenclature:
# -------------
# V - instantaneous specific volume
# R - instantaneous density
# v = V - <V>, fluctuating specific volume
# r = R - <R>, fluctuating density
# b = -<rv>, density-specific-volume covariance
# B - as postfix means the Bousinessq approximation

# <v^2>

# IA: <v^2>B = b/<R>^2
plot "stat.txt" u 2:(-$25/$8/$8), "stat.txt" u 2:48 w l lt 2    # no-mix
plot "stat.txt" u 2:(-$30/$9/$9), "stat.txt" u 2:50 w l lt 2
plot "stat.txt" u 2:(-$35/$10/$10), "stat.txt" u 2:52 w l lt 2  # uniform
plot "stat.txt" u 2:(-$40/$11/$11), "stat.txt" u 2:54 w l lt 2
plot "stat.txt" u 2:(-$45/$12/$12), "stat.txt" u 2:56 w l lt 2  # almost mixed

#################### COLUMNS NOT UPDATED TO LATEST COLUMN NUMBER UNTIL HASHES

# IB: <v^2>(B+1) = b/<R>^2 [ 1 + b(1+b)/(1+(1+b)^2) ]
plot "stat.txt" u 2:(-$23/$8/$8*(1.0+(-$23*(1.0-$23)/(1.0+(1.0-$23)*(1.0-$23))))), "stat.txt" u 2:33 w l lt 2    # no-mix
plot "stat.txt" u 2:(-$25/$9/$9*(1.0+(-$25*(1.0-$25)/(1.0+(1.0-$25)*(1.0-$25))))), "stat.txt" u 2:35 w l lt 2
plot "stat.txt" u 2:(-$27/$10/$10*(1.0+(-$27*(1.0-$27)/(1.0+(1.0-$27)*(1.0-$27))))), "stat.txt" u 2:37 w l lt 2  # uniform
plot "stat.txt" u 2:(-$29/$11/$11*(1.0+(-$29*(1.0-$29)/(1.0+(1.0-$29)*(1.0-$29))))), "stat.txt" u 2:39 w l lt 2
plot "stat.txt" u 2:(-$31/$12/$12*(1.0+(-$31*(1.0-$31)/(1.0+(1.0-$31)*(1.0-$31))))), "stat.txt" u 2:41 w l lt 2  # almost mixed

# <v^3>

# IIA: <v^3>B = -<rv^2>/<R>^2
plot "stat.txt" u 2:(-$24/$8/$8), "stat.txt" u 2:34 w l lt 2     # no-mix
plot "stat.txt" u 2:(-$26/$9/$9), "stat.txt" u 2:36 w l lt 2
plot "stat.txt" u 2:(-$28/$10/$10), "stat.txt" u 2:38 w l lt 2   # uniform
plot "stat.txt" u 2:(-$30/$11/$11), "stat.txt" u 2:40 w l lt 2
plot "stat.txt" u 2:(-$32/$12/$12), "stat.txt" u 2:42 w l lt 2   # almost mixed
# IIA: <v^3>B = (b<V> + <R><v^2>) / <R>^2
plot "stat.txt" u 2:(($23*$13+$8*$33)/$8/$8), "stat.txt" u 2:34 w l lt 2        # no-mix
plot "stat.txt" u 2:(($25*$14+$9*$35)/$9/$9), "stat.txt" u 2:36 w l lt 2
plot "stat.txt" u 2:(($27*$15+$10*$37)/$10/$10), "stat.txt" u 2:38 w l lt 2     # uniform
plot "stat.txt" u 2:(($29*$16+$11*$39)/$11/$11), "stat.txt" u 2:40 w l lt 2
plot "stat.txt" u 2:(($31*$17+$12*$41)/$12/$12), "stat.txt" u 2:42 w l lt 2     # almost mixed

# IIB: <v^3>B = r'/rho2 <v^2>B (1-2<X>)
plot "stat.txt" u 2:(1.0e-2/1.0*(-$23)/$8/$8*(1.0-2.0*$3)), "stat.txt" u 2:34 w l lt 2       # no-mix
plot "stat.txt" u 2:(1.0e-2/1.0*(-$25)/$9/$9*(1.0-2.0*$4)), "stat.txt" u 2:36 w l lt 2
plot "stat.txt" u 2:(1.0e-2/1.0*(-$27)/$10/$10*(1.0-2.0*$5)), "stat.txt" u 2:38 w l lt 2     # uniform
plot "stat.txt" u 2:(1.0e-2/1.0*(-$29)/$11/$11*(1.0-2.0*$6)), "stat.txt" u 2:40 w l lt 2
plot "stat.txt" u 2:(1.0e-2/1.0*(-$31)/$12/$12*(1.0-2.0*$7)), "stat.txt" u 2:42 w l lt 2     # almost mixed

# IIB2: <v^3>(B+1) = r'/rho2 <v^2>(B+1) (1-2<X>)
plot "stat.txt" u 2:(1.0e-2/1.0*(-$23/$8/$8*(1.0+(-$23*(1.0-$23)/(1.0+(1.0-$23)*(1.0-$23)))))*(1.0-2.0*$3)), "stat.txt" u 2:34 w l lt 2       # no-mix
plot "stat.txt" u 2:(1.0e-2/1.0*(-$25/$9/$9*(1.0+(-$25*(1.0-$25)/(1.0+(1.0-$25)*(1.0-$25)))))*(1.0-2.0*$4)), "stat.txt" u 2:36 w l lt 2
plot "stat.txt" u 2:(1.0e-2/1.0*(-$27/$10/$10*(1.0+(-$27*(1.0-$27)/(1.0+(1.0-$27)*(1.0-$27)))))*(1.0-2.0*$5)), "stat.txt" u 2:38 w l lt 2     # uniform
plot "stat.txt" u 2:(1.0e-2/1.0*(-$29/$11/$11*(1.0+(-$29*(1.0-$29)/(1.0+(1.0-$29)*(1.0-$29)))))*(1.0-2.0*$6)), "stat.txt" u 2:40 w l lt 2
plot "stat.txt" u 2:(1.0e-2/1.0*(-$31/$12/$12*(1.0+(-$31*(1.0-$31)/(1.0+(1.0-$31)*(1.0-$31)))))*(1.0-2.0*$7)), "stat.txt" u 2:42 w l lt 2     # almost mixed

# IIC: <v^3>B = 2(1-2<V>)<v^2>B (1-theta)/(2-theta)
plot "stat.txt" u 2:(2.0*(1.0-2.0*$13)*(-$23/$8/$8)*(1.0-(1.0+$23/(-$23+1.0*1.0e-2*$3*(1.0-$3))))/(2.0-(1.0+$23/(-$23+1.0*1.0e-2*$3*(1.0-$3))))), "stat.txt" u 2:34 w l lt 2    # no-mix
plot "stat.txt" u 2:(2.0*(1.0-2.0*$14)*(-$25/$9/$9)*(1.0-(1.0+$25/(-$25+1.0*1.0e-2*$4*(1.0-$4))))/(2.0-(1.0+$25/(-$25+1.0*1.0e-2*$4*(1.0-$4))))), "stat.txt" u 2:36 w l lt 2
plot "stat.txt" u 2:(2.0*(1.0-2.0*$15)*(-$27/$10/$10)*(1.0-(1.0+$27/(-$27+1.0*1.0e-2*$5*(1.0-$5))))/(2.0-(1.0+$27/(-$27+1.0*1.0e-2*$5*(1.0-$5))))), "stat.txt" u 2:38 w l lt 2  # uniform
plot "stat.txt" u 2:(2.0*(1.0-2.0*$16)*(-$29/$11/$11)*(1.0-(1.0+$29/(-$29+1.0*1.0e-2*$6*(1.0-$6))))/(2.0-(1.0+$29/(-$29+1.0*1.0e-2*$6*(1.0-$6))))), "stat.txt" u 2:40 w l lt 2
plot "stat.txt" u 2:(2.0*(1.0-2.0*$17)*(-$31/$12/$12)*(1.0-(1.0+$31/(-$31+1.0*1.0e-2*$7*(1.0-$7))))/(2.0-(1.0+$31/(-$31+1.0*1.0e-2*$7*(1.0-$7))))), "stat.txt" u 2:42 w l lt 2  # almost mixed
# IIC: <v^3>B = 2(1-2<V>)<v^2>(B+1) (1-theta)/(2-theta)
plot "stat.txt" u 2:(2.0*(1.0-2.0*$13)*(-$23/$8/$8*(1.0+(-$23*(1.0-$23)/(1.0+(1.0-$23)*(1.0-$23)))))*(1.0-(1.0+$23/(-$23+1.0*1.0e-2*$3*(1.0-$3))))/(2.0-(1.0+$23/(-$23+1.0*1.0e-2*$3*(1.0-$3))))), "stat.txt" u 2:34 w l lt 2    # no-mix
plot "stat.txt" u 2:(2.0*(1.0-2.0*$14)*(-$25/$9/$9*(1.0+(-$25*(1.0-$25)/(1.0+(1.0-$25)*(1.0-$25)))))*(1.0-(1.0+$25/(-$25+1.0*1.0e-2*$4*(1.0-$4))))/(2.0-(1.0+$25/(-$25+1.0*1.0e-2*$4*(1.0-$4))))), "stat.txt" u 2:36 w l lt 2
plot "stat.txt" u 2:(2.0*(1.0-2.0*$15)*(-$27/$10/$10*(1.0+(-$27*(1.0-$27)/(1.0+(1.0-$27)*(1.0-$27)))))*(1.0-(1.0+$27/(-$27+1.0*1.0e-2*$5*(1.0-$5))))/(2.0-(1.0+$27/(-$27+1.0*1.0e-2*$5*(1.0-$5))))), "stat.txt" u 2:38 w l lt 2  # uniform
plot "stat.txt" u 2:(2.0*(1.0-2.0*$16)*(-$29/$11/$11*(1.0+(-$29*(1.0-$29)/(1.0+(1.0-$29)*(1.0-$29)))))*(1.0-(1.0+$29/(-$29+1.0*1.0e-2*$6*(1.0-$6))))/(2.0-(1.0+$29/(-$29+1.0*1.0e-2*$6*(1.0-$6))))), "stat.txt" u 2:40 w l lt 2
plot "stat.txt" u 2:(2.0*(1.0-2.0*$17)*(-$31/$12/$12*(1.0+(-$31*(1.0-$31)/(1.0+(1.0-$31)*(1.0-$31)))))*(1.0-(1.0+$31/(-$31+1.0*1.0e-2*$7*(1.0-$7))))/(2.0-(1.0+$31/(-$31+1.0*1.0e-2*$7*(1.0-$7))))), "stat.txt" u 2:42 w l lt 2  # almost mixed
# IIC: <v^3>B = 2(1-2<V>)<v^2>B (1-thetav)/(2-thetav)
plot "stat.txt" u 2:(2.0*(1.0-2.0*$13)*(-$23/$8/$8)*(1.0-(1.0-$33/($33+1.0*1.0e-2/1.0/0.99*$3*(1.0-$3))))/(2.0-(1.0-$33/($33+1.0*1.0e-2/1.0/0.99*$3*(1.0-$3))))), "stat.txt" u 2:34 w l lt 2    # no-mix
plot "stat.txt" u 2:(2.0*(1.0-2.0*$14)*(-$25/$9/$9)*(1.0-(1.0-$35/($35+1.0*1.0e-2/1.0/0.99*$4*(1.0-$4))))/(2.0-(1.0-$35/($35+1.0*1.0e-2/1.0/0.99*$4*(1.0-$4))))), "stat.txt" u 2:36 w l lt 2
plot "stat.txt" u 2:(2.0*(1.0-2.0*$15)*(-$27/$10/$10)*(1.0-(1.0-$37/($37+1.0*1.0e-2/1.0/0.99*$5*(1.0-$5))))/(2.0-(1.0-$37/($37+1.0*1.0e-2/1.0/0.99*$5*(1.0-$5))))), "stat.txt" u 2:38 w l lt 2  # uniform
plot "stat.txt" u 2:(2.0*(1.0-2.0*$16)*(-$29/$11/$11)*(1.0-(1.0-$39/($39+1.0*1.0e-2/1.0/0.99*$6*(1.0-$6))))/(2.0-(1.0-$39/($39+1.0*1.0e-2/1.0/0.99*$6*(1.0-$6))))), "stat.txt" u 2:40 w l lt 2
plot "stat.txt" u 2:(2.0*(1.0-2.0*$17)*(-$31/$12/$12)*(1.0-(1.0-$41/($41+1.0*1.0e-2/1.0/0.99*$7*(1.0-$7))))/(2.0-(1.0-$41/($41+1.0*1.0e-2/1.0/0.99*$7*(1.0-$7))))), "stat.txt" u 2:42 w l lt 2  # almost mixed
# IIC: <v^3>B = 2(1-2<V>)<v^2>(B+1) (1-thetav)/(2-thetav)
plot "stat.txt" u 2:(2.0*(1.0-2.0*$13)*(-$23/$8/$8*(1.0+(-$23*(1.0-$23)/(1.0+(1.0-$23)*(1.0-$23)))))*(1.0-(1.0-$33/($33+1.0*1.0e-2/1.0/0.99*$3*(1.0-$3))))/(2.0-(1.0-$33/($33+1.0*1.0e-2/1.0/0.99*$3*(1.0-$3))))), "stat.txt" u 2:34 w l lt 2    # no-mix
plot "stat.txt" u 2:(2.0*(1.0-2.0*$14)*(-$25/$9/$9*(1.0+(-$25*(1.0-$25)/(1.0+(1.0-$25)*(1.0-$25)))))*(1.0-(1.0-$35/($35+1.0*1.0e-2/1.0/0.99*$4*(1.0-$4))))/(2.0-(1.0-$35/($35+1.0*1.0e-2/1.0/0.99*$4*(1.0-$4))))), "stat.txt" u 2:36 w l lt 2
plot "stat.txt" u 2:(2.0*(1.0-2.0*$15)*(-$27/$10/$10*(1.0+(-$27*(1.0-$27)/(1.0+(1.0-$27)*(1.0-$27)))))*(1.0-(1.0-$37/($37+1.0*1.0e-2/1.0/0.99*$5*(1.0-$5))))/(2.0-(1.0-$37/($37+1.0*1.0e-2/1.0/0.99*$5*(1.0-$5))))), "stat.txt" u 2:38 w l lt 2  # uniform
plot "stat.txt" u 2:(2.0*(1.0-2.0*$16)*(-$29/$11/$11*(1.0+(-$29*(1.0-$29)/(1.0+(1.0-$29)*(1.0-$29)))))*(1.0-(1.0-$39/($39+1.0*1.0e-2/1.0/0.99*$6*(1.0-$6))))/(2.0-(1.0-$39/($39+1.0*1.0e-2/1.0/0.99*$6*(1.0-$6))))), "stat.txt" u 2:40 w l lt 2
plot "stat.txt" u 2:(2.0*(1.0-2.0*$17)*(-$31/$12/$12*(1.0+(-$31*(1.0-$31)/(1.0+(1.0-$31)*(1.0-$31)))))*(1.0-(1.0-$41/($41+1.0*1.0e-2/1.0/0.99*$7*(1.0-$7))))/(2.0-(1.0-$41/($41+1.0*1.0e-2/1.0/0.99*$7*(1.0-$7))))), "stat.txt" u 2:42 w l lt 2  # almost mixed

####################

# <rv^3> = <rv><v^2>
plot "stat.txt" u 2:($25*$48), "stat.txt" u 2:27 w l lt 2  # no-mix
plot "stat.txt" u 2:($30*$50), "stat.txt" u 2:32 w l lt 2
plot "stat.txt" u 2:($35*$52), "stat.txt" u 2:37 w l lt 2  # uniform
plot "stat.txt" u 2:($40*$54), "stat.txt" u 2:42 w l lt 2
plot "stat.txt" u 2:($45*$56), "stat.txt" u 2:47 w l lt 2  # almost mixed

# <r^2v^2> = <rv>^2
plot "stat.txt" u 2:($25*$25), "stat.txt" u 2:24 w l lt 2  # no-mix
plot "stat.txt" u 2:($30*$30), "stat.txt" u 2:29 w l lt 2
plot "stat.txt" u 2:($35*$35), "stat.txt" u 2:34 w l lt 2  # uniform
plot "stat.txt" u 2:($40*$40), "stat.txt" u 2:39 w l lt 2
plot "stat.txt" u 2:($45*$45), "stat.txt" u 2:43 w l lt 2  # almost mixed

# <v^2> = b/<R>^2 ( 1-b^2-b^3 )
plot "stat.txt" u 2:(-$25/$8/$8*(1.0-$25*$25-$25*$25*$25)), "stat.txt" u 2:48 w l lt 2  # no-mix
plot "stat.txt" u 2:(-$30/$9/$9*(1.0-$30*$30-$30*$30*$30)), "stat.txt" u 2:50 w l lt 2
plot "stat.txt" u 2:(-$35/$10/$10*(1.0-$35*$35-$35*$35*$35)), "stat.txt" u 2:52 w l lt 2  # uniform
plot "stat.txt" u 2:(-$40/$11/$11*(1.0-$40*$40-$40*$40*$40)), "stat.txt" u 2:54 w l lt 2
plot "stat.txt" u 2:(-$45/$12/$12*(1.0-$45*$45-$45*$45*$45)), "stat.txt" u 2:56 w l lt 2  # almost mixed

# <rv^2>/<R>/<V>^2 = b^2(2+b)/(1+b)^2
plot "stat.txt" u 2:($25*$25*(2.0-$25)/(1.0-$25)/(1.0-$25)), "stat.txt" u 2:($26/$8/$13/$13) w l lt 2  # no-mix
plot "stat.txt" u 2:($30*$30*(2.0-$30)/(1.0-$30)/(1.0-$30)), "stat.txt" u 2:($31/$9/$14/$14) w l lt 2
plot "stat.txt" u 2:($35*$35*(2.0-$35)/(1.0-$35)/(1.0-$35)), "stat.txt" u 2:($36/$10/$15/$15) w l lt 2  # uniform
plot "stat.txt" u 2:($40*$40*(2.0-$40)/(1.0-$40)/(1.0-$40)), "stat.txt" u 2:($41/$11/$16/$16) w l lt 2
plot "stat.txt" u 2:($45*$45*(2.0-$45)/(1.0-$45)/(1.0-$45)), "stat.txt" u 2:($46/$12/$17/$17) w l lt 2  # almost mixed

# <r^2v> = -<R>^2 <rv^2>
plot "stat.txt" u 2:(-$8*$8*$26), "stat.txt" u 2:23 w l lt 2  # no-mix
plot "stat.txt" u 2:(-$9*$9*$31), "stat.txt" u 2:28 w l lt 2
plot "stat.txt" u 2:(-$10*$10*$36), "stat.txt" u 2:33 w l lt 2  # uniform
plot "stat.txt" u 2:(-$11*$11*$41), "stat.txt" u 2:38 w l lt 2
plot "stat.txt" u 2:(-$12*$12*$46), "stat.txt" u 2:43 w l lt 2  # almost mixed

# <Rv^3>/<R>/<V>^3 = -b^2 (3+b)/(1+b)^3
plot "stat.txt" u 2:(-$25*$25*(3.0-$25)/(1.0-$25)**3.0), "stat.txt" u 2:(($8*$49+$27)/$8/$13/$13/$13) w l lt 2  # no-mix
plot "stat.txt" u 2:(-$30*$30*(3.0-$30)/(1.0-$30)**3.0), "stat.txt" u 2:(($9*$51+$32)/$9/$14/$14/$14) w l lt 2
plot "stat.txt" u 2:(-$35*$35*(3.0-$35)/(1.0-$35)**3.0), "stat.txt" u 2:(($10*$53+$37)/$10/$15/$15/$15) w l lt 2  # uniform
plot "stat.txt" u 2:(-$40*$40*(3.0-$40)/(1.0-$40)**3.0), "stat.txt" u 2:(($11*$55+$42)/$11/$16/$16/$16) w l lt 2
plot "stat.txt" u 2:(-$45*$45*(3.0-$45)/(1.0-$45)**3.0), "stat.txt" u 2:(($12*$57+$47)/$12/$17/$17/$17) w l lt 2  # almost mixed

# <v^2> vs. rho2 * <v^3>
plot "stat.txt" u 2:48, "stat.txt" u 2:49 w l lt 2 # no-mix
plot "stat.txt" u 2:50, "stat.txt" u 2:51 w l lt 2
plot "stat.txt" u 2:52, "stat.txt" u 2:53 w l lt 2 # uniform
plot "stat.txt" u 2:54, "stat.txt" u 2:55 w l lt 2
plot "stat.txt" u 2:56, "stat.txt" u 2:57 w l lt 2 # almost mixed

# <v^2> = (1+r') b /<R>^2
plot "stat.txt" u 2:((1.0+9)*(-$25)/$8/$8), "stat.txt" u 2:48 w l lt 2 # no-mix
plot "stat.txt" u 2:((1.0+9)*(-$30)/$9/$9), "stat.txt" u 2:50 w l lt 2
plot "stat.txt" u 2:((1.0+9)*(-$35)/$10/$10), "stat.txt" u 2:52 w l lt 2 # uniform
plot "stat.txt" u 2:((1.0+9)*(-$40)/$11/$11), "stat.txt" u 2:54 w l lt 2
plot "stat.txt" u 2:((1.0+9)*(-$45)/$12/$12), "stat.txt" u 2:56 w l lt 2 # almost mixed

# <v^2> = b(1+b)[(r-2)/r]/<R>^2 + b/(r <R> rho1)
plot "stat.txt" u 2:(-$25*(1.0-$25)*(0.0101-2.0)/0.0101/$8/$8-$25/0.0101/$8/0.99), "stat.txt" u 2:48 w l lt 2 # no-mix
plot "stat.txt" u 2:(-$30*(1.0-$30)*(0.0101-2.0)/0.0101/$9/$9-$30/0.0101/$9/0.99), "stat.txt" u 2:50 w l lt 2
plot "stat.txt" u 2:(-$35*(1.0-$35)*(0.0101-2.0)/0.0101/$10/$10-$35/0.0101/$10/0.99), "stat.txt" u 2:52 w l lt 2 # uniform
plot "stat.txt" u 2:(-$40*(1.0-$30)*(0.0101-2.0)/0.0101/$11/$11-$40/0.0101/$11/0.99), "stat.txt" u 2:54 w l lt 2
plot "stat.txt" u 2:(-$45*(1.0-$45)*(0.0101-2.0)/0.0101/$12/$12-$45/0.0101/$12/0.99), "stat.txt" u 2:56 w l lt 2 # almost mixed

# <v^2> =  [(2+r)/(1+r)]b/(<R> rho1) - b(1+b)/<R>^2
plot "stat.txt" u 2:((2.0+9.0)/(1.0+9.0)*(-$25)/$8/0.1+$25*(1.0-$25)/$8/$8), "stat.txt" u 2:48 w l lt 2 # no-mix
plot "stat.txt" u 2:((2.0+9.0)/(1.0+9.0)*(-$30)/$9/0.1+$30*(1.0-$30)/$9/$9), "stat.txt" u 2:50 w l lt 2
plot "stat.txt" u 2:((2.0+9.0)/(1.0+9.0)*(-$35)/$10/0.1+$35*(1.0-$35)/$10/$12), "stat.txt" u 2:52 w l lt 2 # uniform
plot "stat.txt" u 2:((2.0+9.0)/(1.0+9.0)*(-$40)/$11/0.1+$40*(1.0-$40)/$11/$12), "stat.txt" u 2:54 w l lt 2
plot "stat.txt" u 2:((2.0+9.0)/(1.0+9.0)*(-$45)/$12/0.1+$45*(1.0-$45)/$12/$12), "stat.txt" u 2:56 w l lt 2 # almost mixed

# <v^2> = [b/<R>^2][1 + (b/b_{nm})[<R>^2/(rho_1rho_2) - 1]] with b_{nm} = r^2/(1+r) [X(1-X)]
plot "stat.txt" u 2:((-$25/$8/$8)*(1.0-$25/(9.0*9.0/(1.0+9.0)*$3*(1.0-$3))*($8*$8/0.1-1.0))), "stat.txt" u 2:48 w l lt 2 # no-mix
plot "stat.txt" u 2:((-$30/$9/$9)*(1.0-$30/(9.0*9.0/(1.0+9.0)*$4*(1.0-$4))*($9*$9/0.1-1.0))), "stat.txt" u 2:50 w l lt 2
plot "stat.txt" u 2:((-$35/$10/$10)*(1.0-$35/(9.0*9.0/(1.0+9.0)*$5*(1.0-$5))*($10*$10/0.1-1.0))), "stat.txt" u 2:52 w l lt 2 # uniform
plot "stat.txt" u 2:((-$40/$11/$11)*(1.0-$40/(9.0*9.0/(1.0+9.0)*$6*(1.0-$6))*($11*$11/0.1-1.0))), "stat.txt" u 2:54 w l lt 2
plot "stat.txt" u 2:((-$45/$12/$12)*(1.0-$45/(9.0*9.0/(1.0+9.0)*$7*(1.0-$7))*($12*$12/0.1-1.0))), "stat.txt" u 2:56 w l lt 2 # almost mixed

# (b/b_{nm})
plot "stat.txt" u 2:(-$25/(9.0*9.0/(1.0+9.0)*$3*(1.0-$3)))

# <v^2> = b(1+b)/<R>^2 + (b/<R>) [ (2+r)/rho_2 - 2(1+b)/<R> ] (b/b_{nm})^{1/2} where b_{nm} = r^2/(1+r) [X(1-X)], high At
plot "stat.txt" u 2:((-$25*(1.0-$25)/$8/$8 - $25/$8 * ((2.0+9.0)-2.0*(1.0-$25)/$8)) * (-$25/(9.0*9.0/(1.0+9.0)*($3*(1.0-$3))))**(1.0/2.0)), "stat.txt" u 2:48 w l lt 2 # no-mix
plot "stat.txt" u 2:((-$30*(1.0-$30)/$9/$9 - $30/$9 * ((2.0+9.0)-2.0*(1.0-$30)/$9)) * (-$30/(9.0*9.0/(1.0+9.0)*($4*(1.0-$4))))**(1.0/2.0)), "stat.txt" u 2:50 w l lt 2
plot "stat.txt" u 2:((-$35*(1.0-$35)/$10/$10 - $35/$10 * ((2.0+9.0)-2.0*(1.0-$35)/$10)) * (-$35/(9.0*9.0/(1.0+9.0)*($5*(1.0-$5))))**(1.0/2.0)), "stat.txt" u 2:52 w l lt 2 # uniform
plot "stat.txt" u 2:((-$40*(1.0-$40)/$11/$11 - $40/$11 * ((2.0+9.0)-2.0*(1.0-$40)/$11)) * (-$40/(9.0*9.0/(1.0+9.0)*($6*(1.0-$6))))**(1.0/2.0)), "stat.txt" u 2:54 w l lt 2
plot "stat.txt" u 2:((-$45*(1.0-$45)/$12/$12 - $45/$12 * ((2.0+9.0)-2.0*(1.0-$45)/$12)) * (-$45/(9.0*9.0/(1.0+9.0)*($7*(1.0-$7))))**(1.0/2.0)), "stat.txt" u 2:56 w l lt 2 # almost mixed

# <v^2> = b(1+b)/<R>^2 + (b/<R>) [ (2+r)/rho_2 - 2(1+b)/<R> ] (b/b_{nm})^{1/2} where b_{nm} = r^2/(1+r) [X(1-X)], low At
plot "stat.txt.low" u 2:((-$25*(1.0-$25)/$8/$8 - $25/$8 * ((2.0+0.0101)-2.0*(1.0-$25)/$8)) * (-$25/(0.0101*0.0101/(1.0+0.0101)*($3*(1.0-$3))))**(1.0/2.0)), "stat.txt.low" u 2:48 w l lt 2 # no-mix
plot "stat.txt.low" u 2:((-$30*(1.0-$30)/$9/$9 - $30/$9 * ((2.0+0.0101)-2.0*(1.0-$30)/$9)) * (-$30/(0.0101*0.0101/(1.0+0.0101)*($4*(1.0-$4))))**(1.0/2.0)), "stat.txt.low" u 2:50 w l lt 2
plot "stat.txt.low" u 2:((-$35*(1.0-$35)/$10/$10 - $35/$10 * ((2.0+0.0101)-2.0*(1.0-$35)/$10)) * (-$35/(0.0101*0.0101/(1.0+0.0101)*($5*(1.0-$5))))**(1.0/2.0)), "stat.txt.low" u 2:52 w l lt 2 # uniform
plot "stat.txt.low" u 2:((-$40*(1.0-$40)/$11/$11 - $40/$11 * ((2.0+0.0101)-2.0*(1.0-$40)/$11)) * (-$40/(0.0101*0.0101/(1.0+0.0101)*($6*(1.0-$6))))**(1.0/2.0)), "stat.txt.low" u 2:54 w l lt 2
plot "stat.txt.low" u 2:((-$45*(1.0-$45)/$12/$12 - $45/$12 * ((2.0+0.0101)-2.0*(1.0-$45)/$12)) * (-$45/(0.0101*0.0101/(1.0+0.0101)*($7*(1.0-$7))))**(1.0/2.0)), "stat.txt.low" u 2:56 w l lt 2 # almost mixed

# b = (rho2 r)^2<x^2> / <R>^2 * [1 + <x^2>/<x^2>_{nm} [<R>^2/(rho_1rho_2) - 1]] with <x^2>_{nm} = [X(1-X)], low At
plot "stat.txt" u 2:(0.01*0.01*$18/$8/$8 * (1.0 + $18/($3*(1.0-$3)) * ($8*$8/0.99 - 1.0))), "stat.txt" u 2:(-$25) w l lt 2 # no-mix
plot "stat.txt" u 2:(0.01*0.01*$19/$9/$9 * (1.0 + $19/($4*(1.0-$4)) * ($9*$9/0.99 - 1.0))), "stat.txt" u 2:(-$30) w l lt 2
plot "stat.txt" u 2:(0.01*0.01*$20/$10/$10 * (1.0 + $20/($5*(1.0-$5)) * ($10*$10/0.99 - 1.0))), "stat.txt" u 2:(-$35) w l lt 2 # uniform
plot "stat.txt" u 2:(0.01*0.01*$21/$11/$11 * (1.0 + $21/($6*(1.0-$6)) * ($11*$11/0.99 - 1.0))), "stat.txt" u 2:(-$40) w l lt 2 #
plot "stat.txt" u 2:(0.01*0.01*$22/$12/$12 * (1.0 + $22/($7*(1.0-$7)) * ($12*$12/0.99 - 1.0))), "stat.txt" u 2:(-$45) w l lt 2 # almost mixed

# b = (rho2 r)^2<x^2> / <R>^2 * [1 + <x^2>/<x^2>_{nm} [<R>^2/(rho_1rho_2) - 1]] with <x^2>_{nm} = [X(1-X)], high At
plot "stat.txt" u 2:(0.9*0.9*$18/$8/$8 * (1.0 + $18/($3*(1.0-$3)) * ($8*$8/0.1 - 1.0))), "stat.txt" u 2:(-$25) w l lt 2 # no-mix
plot "stat.txt" u 2:(0.9*0.9*$19/$9/$9 * (1.0 + $19/($4*(1.0-$4)) * ($9*$9/0.1 - 1.0))), "stat.txt" u 2:(-$30) w l lt 2
plot "stat.txt" u 2:(0.9*0.9*$20/$10/$10 * (1.0 + $20/($5*(1.0-$5)) * ($10*$10/0.1 - 1.0))), "stat.txt" u 2:(-$35) w l lt 2 # uniform
plot "stat.txt" u 2:(0.9*0.9*$21/$11/$11 * (1.0 + $21/($6*(1.0-$6)) * ($11*$11/0.1 - 1.0))), "stat.txt" u 2:(-$40) w l lt 2 #
plot "stat.txt" u 2:(0.9*0.9*$22/$12/$12 * (1.0 + $22/($7*(1.0-$7)) * ($12*$12/0.1 - 1.0))), "stat.txt" u 2:(-$45) w l lt 2 # almost mixed