Examples » Walker: Integrating the mass-fraction beta SDE

This example runs Walker to integrate the mass-fraction beta SDE (see DiffEq/MassFractionBeta.h) using constant coefficients. The mass-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 <y^2> and <rho v>"

walker

  #nstep 1      # Max number of time steps
  term  25.0    # Max time
  dt    0.002    # Time step size
  npar  10000   # Number of particles
  ttyi  500      # TTY output interval

  rngs
    mkl_r250 end
  end

  massfracbeta  # Select the mass-fraction beta SDE
    depvar Y    # Dependent variable: Y, denoting mass fractions
    ncomp 15    # = 3x5 = 5 systems each computing 3 variables:
                #   Y - mass 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
    #r 0.0101 0.0101 0.0101 0.0101 0.0101 end # low-A
    r 9.0 9.0 9.0 9.0 9.0 end   # high-A
  end

  statistics
    #precision 5
    #format    default
    # <Y>, mass fraction means
    <Y1>        # 3
    <Y2>        # 4
    <Y3>        # 5
    <Y4>        # 6
    <Y5>        # 7
    # <rho>, mean densities
    <Y6>        # 8
    <Y7>        # 9
    <Y8>        # 10
    <Y9>        # 11
    <Y10>       # 12
    # <V>, mean specific volumes
    <Y11>       # 13
    <Y12>       # 14
    <Y13>       # 15
    <Y14>       # 16
    <Y15>       # 17
    # <y^2>, mass fraction variances
    <y1y1>      # 23
    <y2y2>      # 24
    <y3y3>      # 25
    <y4y4>      # 26
    <y5y5>      # 27
     # <rho^2>, density variances
    <y6y6>      # 28
    <y7y7>      # 30
    <y8y8>      # 32
    <y9y9>      # 34
    <y10y10>    # 36
    # <v^2>, specific volume variances
    <y11y11>    # 38
    <y12y12>    # 39
    <y13y13>    # 40
    <y14y14>    # 41
    <y15y15>    # 42
     # <rho v>, density-specific-volume covariances
    <y6y11>     # 29
    <y7y12>     # 31
    <y8y13>     # 33
    <y9y14>     # 35
    <y10y15>    # 37
     # <rho v^2>
    <Y6y11y11>  # 18
    <Y7y12y12>  # 19
    <Y8y13y13>  # 20
    <Y9y14y14>  # 21
    <Y10y15y15> # 22
  end

  pdfs
    interval  100
    filetype  txt
    policy    overwrite
    centering elem
    format    scientific
    precision 4
    p1( Y1 : 1.0e-2 )
    p2( Y2 : 1.0e-2 )
    p3( Y3 : 1.0e-2 )
    p4( Y4 : 1.0e-2 )
    p5( Y5 : 1.0e-2 )
  end
end

Example run on 8 CPUs

./charmrun +p8 Main/walker -v -c ../../tmp/massfracbeta.q

Output

Running on 8 processors:  Main/walker -v -c massfracbeta.q
charmrun>  /usr/bin/setarch x86_64 -R  mpirun -np 8  Main/walker -v -c massfracbeta.q
Charm++> Running on MPI version: 3.0
Charm++> level of thread support used: MPI_THREAD_SINGLE (desired: MPI_THREAD_SINGLE)
Charm++> Running in non-SMP mode: numPes 8
Converse/Charm++ Commit ID: 63927de
CharmLB> Load balancer assumes all CPUs are same.
Charm++> Running on 1 unique compute nodes (8-way SMP).
Charm++> cpu topology info is gathered in 0.002 seconds.

      ,::,`                                                            `.
   .;;;'';;;:                                                          ;;#
  ;;;@+   +;;;  ;;;;;,   ;;;;. ;;;;;, ;;;;      ;;;;   `;;;;;;:        ;;;
 :;;@`     :;;' .;;;@,    ,;@, ,;;;@: .;;;'     .;+;. ;;;@#:';;;      ;;;;'
 ;;;#       ;;;: ;;;'      ;:   ;;;'   ;;;;;     ;#  ;;;@     ;;;     ;+;;'
.;;+        ;;;# ;;;'      ;:   ;;;'   ;#;;;`    ;#  ;;@      `;;+   .;#;;;.
;;;#        :;;' ;;;'      ;:   ;;;'   ;# ;;;    ;# ;;;@       ;;;   ;# ;;;+
;;;#        .;;; ;;;'      ;:   ;;;'   ;# ,;;;   ;# ;;;#       ;;;:  ;@  ;;;
;;;#        .;;' ;;;'      ;:   ;;;'   ;#  ;;;;  ;# ;;;'       ;;;+ ;',  ;;;@
;;;+        ,;;+ ;;;'      ;:   ;;;'   ;#   ;;;' ;# ;;;'       ;;;' ;':::;;;;
`;;;        ;;;@ ;;;'      ;:   ;;;'   ;#    ;;;';# ;;;@       ;;;:,;+++++;;;'
 ;;;;       ;;;@ ;;;#     .;.   ;;;'   ;#     ;;;;# `;;+       ;;# ;#     ;;;'
 .;;;      :;;@  ,;;+     ;+    ;;;'   ;#      ;;;#  ;;;      ;;;@ ;@      ;;;.
  ';;;    ;;;@,   ;;;;``.;;@    ;;;'   ;+      .;;#   ;;;    :;;@ ;;;      ;;;+
   :;;;;;;;+@`     ';;;;;'@    ;;;;;, ;;;;      ;;+    +;;;;;;#@ ;;;;.   .;;;;;;
     .;;#@'         `#@@@:     ;::::; ;::::      ;@      '@@@+   ;:::;    ;::::::
    :;;;;;;.      __      __        .__   __
  .;@+@';;;;;;'  /  \    /  \_____  |  | |  | __ ___________
    `     '#''@` \   \/\/   /\__  \ |  | |  |/ // __ \_  __ \
                  \        /  / __ \|  |_|    <\  ___/|  | \/
                   \__/\  /  (____  /____/__|_ \\___  >__|
                        \/        \/          \/    \/

                                 < ENVIRONMENT >
                                 ------ o ------

 * Build environment:
 --------------------
     Hostname                       : karman
     Executable                     : walker
     Version                        : 0.1
     Release                        : LA-CC-XX-XXX
     Revision                       : ab18ecb7bd1b27706d963643c4c1ae144a21bd22
     CMake build type               : DEBUG
     Asserts                        : on (turn off: CMAKE_BUILD_TYPE=RELEASE)
     Exception trace                : on (turn off: CMAKE_BUILD_TYPE=RELEASE)
     MPI C++ wrapper                : /opt/openmpi/1.8/clang/system/bin/mpicxx
     Underlying C++ compiler        : /usr/bin/clang++-3.5
     Build date                     : Mon Aug  3 15:17:35 MDT 2015

 * Run-time environment:
 -----------------------
     Date, time                     : Tue Aug  4 08:01:36 2015
     Work directory                 : /home/jbakosi/code/quinoa/build/clang
     Executable (rel. to work dir)  : Main/walker
     Command line arguments         : '-v -c massfracbeta.q'
     Output                         : verbose (quiet: omit -v)
     Control file                   : massfracbeta.q
     Parsed control file            : success


                                   < FACTORY >
                                   ---- o ----

 * Particle properties data layout policy (CMake: LAYOUT):
 ---------------------------------------------------------
     particle-major                

 * Registered differential equations:
 ------------------------------------
     Unique equation types          : 12
     With all policy combinations   : 56

     Legend: equation name : supported policies

     Policy codes:
      * i: initialization policy:
        R - raw
        Z - zero
        D - delta
        B - beta
      * c: coefficients policy:
        C - const
        D - decay
        H - homogeneous decay
        M - Monte Carlo homogeneous decay

     Beta                           : i:BDRZ, c:C
     Diagonal Ornstein-Uhlenbeck    : i:BDRZ, c:C
     Dirichlet                      : i:BDRZ, c:C
     Gamma                          : i:BDRZ, c:C
     Generalized Dirichlet          : i:BDRZ, c:C
     Mass-fraction beta             : i:BDRZ, c:C
     Mix mass-fraction beta         : i:BDRZ, c:DHM
     Mix number-fraction beta       : i:BDRZ, c:D
     Number-fraction beta           : i:BDRZ, c:C
     Ornstein-Uhlenbeck             : i:BDRZ, c:C
     Skew-Normal                    : i:BDRZ, c:C
     Wright-Fisher                  : i:BDRZ, c:C


                                   < PROBLEM >
                                   ---- o ----

 * Title: Test Ray's closure ideas for <y^2> and <rho v>
 -------------------------------------------------------

 * Differential equations integrated (1):
 ----------------------------------------
   < Mass-fraction beta >
     kind                           : stochastic
     dependent variable             : Y
     initialization policy          : Z
     coefficients policy            : C
     start offset in particle array : 0
     number of components           : 5
     random number generator        : MKL R250
     coeff b [5]                    : { 0.4 1 1 1 100 }
     coeff S [5]                    : { 0.5 0.53846 0.5 0.39683 0.5 }
     coeff kappa [5]                : { 2 0.76923 0.5 0.15873 0.1 }
     coeff rho2 [5]                 : { 1 1 1 1 1 }
     coeff r [5]                    : { 9 9 9 9 9 }

 * Output filenames:
 -------------------
     Statistics                     : stat.txt
     PDF                            : pdf

 * Discretization parameters:
 ----------------------------
     Number of time steps           : 18446744073709551615
     Terminate time                 : 25
     Initial time step size         : 0.002

 * Output intervals:
 -------------------
     TTY                            : 500
     Statistics                     : 1
     PDF                            : 100

 * Statistical moments and distributions:
 ----------------------------------------
     Estimated statistical moments  : <Y1> <Y2> <Y3> <Y4> <Y5> <Y6> <Y6y11y11> <Y7> <Y7y12y12> <Y8> <Y8y13y13> <Y9> <Y9y14y14> <Y10> <Y10y15y15> <Y11> <Y12> <Y13> <Y14> <Y15> <y1y1> <y2y2> <y3y3> <y4y4> <y5y5> <y6y6> <y6y11> <y7y7> <y7y12> <y8y8> <y8y13> <y9y9> <y9y14> <y10y10> <y10y15> <y11y11> <y12y12> <y13y13> <y14y14> <y15y15>
     Stats floating-point format    : default
     Stats text precision, digits   : 6
     Estimated PDFs                 : p1(Y1:0.01) p2(Y2:0.01) p3(Y3:0.01) p4(Y4:0.01) p5(Y5:0.01)
     PDF output file type           : txt
     PDF output file policy         : overwrite
     PDF output file centering      : elem
     PDF text floating-point format : scientific
     PDF text precision, digits     : 4

 * Load distribution:
 --------------------
     Virtualization [0.0...1.0]     : 0
     Number of processing elements  : 8
     Number of work units           : 8
     User load (# of particles)     : 10000
     Chunksize (load per work unit) : 1250
     Actual load (# of particles)   : 10000 (=8*1250)

 * Time integration: Differential equations testbed
 --------------------------------------------------
     Legend: it - iteration count
              t - time
             dt - time step size
            ETE - estimated time elapsed (h:m:s)
            ETA - estimated time for accomplishment (h:m:s)
            out - output-saved flags (S: statistics, P: PDFs)
     
      it             t            dt        ETE        ETA   out
 ---------------------------------------------------------------
     500  1.000000e+00  2.000000e-03  000:00:02  000:01:02  S
    1000  2.000000e+00  2.000000e-03  000:00:04  000:00:54  S
    1500  3.000000e+00  2.000000e-03  000:00:06  000:00:50  S
    2000  4.000000e+00  2.000000e-03  000:00:09  000:00:47  S
    2500  5.000000e+00  2.000000e-03  000:00:11  000:00:46  S
    3000  6.000000e+00  2.000000e-03  000:00:13  000:00:44  S
    3500  7.000000e+00  2.000000e-03  000:00:16  000:00:41  S
    4000  8.000000e+00  2.000000e-03  000:00:18  000:00:39  S
    4500  9.000000e+00  2.000000e-03  000:00:20  000:00:36  S
    5000  1.000000e+01  2.000000e-03  000:00:23  000:00:34  S
    5500  1.100000e+01  2.000000e-03  000:00:25  000:00:32  S
    6000  1.200000e+01  2.000000e-03  000:00:28  000:00:30  S
    6500  1.300000e+01  2.000000e-03  000:00:30  000:00:27  S
    7000  1.400000e+01  2.000000e-03  000:00:32  000:00:25  S
    7500  1.500000e+01  2.000000e-03  000:00:35  000:00:23  S
    8000  1.600000e+01  2.000000e-03  000:00:37  000:00:21  S
    8500  1.700000e+01  2.000000e-03  000:00:40  000:00:18  S
    9000  1.800000e+01  2.000000e-03  000:00:42  000:00:16  S
    9500  1.900000e+01  2.000000e-03  000:00:44  000:00:14  S
   10000  2.000000e+01  2.000000e-03  000:00:46  000:00:11  S
   10500  2.100000e+01  2.000000e-03  000:00:49  000:00:09  S
   11000  2.200000e+01  2.000000e-03  000:00:51  000:00:07  S
   11500  2.300000e+01  2.000000e-03  000:00:53  000:00:04  S
   12000  2.400000e+01  2.000000e-03  000:00:56  000:00:02  S
   12500  2.500000e+01  2.000000e-03  000:00:58  000:00:00  S

     Normal finish, maximum time reached: 25.000000

 * Timers (h:m:s):
 -----------------
     Migrate global-scope data                                                   : 0:0:0
     Initial conditions                                                          : 0:0:0
     Total runtime                                                               : 0:0:58

[Partition 0][Node 0] End of program

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
# nm - no-mix value

# <y^2> / <y^2>nm = <r^2> / <r^2>nm = b / bnm, low At
plot "stat.txt" u 2:($23/($3*(1.0-$3))) w l t "<y^2>/<y^2>nm", "stat.txt" u 2:($28/((0.0101**2)/((1.0+0.0101)**3)*$3*(1.0-$3))) w l t "<r^2>/<r^2>nm", "stat.txt" u 2:(-$29/($3*(1.0+$3)*0.0101*0.0101/(1.0+0.0101))) w l t "b/bnm" # no-mix
plot "stat.txt" u 2:($24/($4*(1.0-$4))) w l t "<y^2>/<y^2>nm", "stat.txt" u 2:($30/((0.0101**2)/((1.0+0.0101)**3)*$4*(1.0-$4))) w l t "<r^2>/<r^2>nm", "stat.txt" u 2:(-$31/($4*(1.0+$4)*0.0101*0.0101/(1.0+0.0101))) w l t "b/bnm"
plot "stat.txt" u 2:($25/($5*(1.0-$5))) w l t "<y^2>/<y^2>nm", "stat.txt" u 2:($32/((0.0101**2)/((1.0+0.0101)**3)*$5*(1.0-$5))) w l t "<r^2>/<r^2>nm", "stat.txt" u 2:(-$33/($5*(1.0+$5)*0.0101*0.0101/(1.0+0.0101))) w l t "b/bnm" # uniform
plot "stat.txt" u 2:($26/($6*(1.0-$6))) w l t "<y^2>/<y^2>nm", "stat.txt" u 2:($34/((0.0101**2)/((1.0+0.0101)**3)*$6*(1.0-$6))) w l t "<r^2>/<r^2>nm", "stat.txt" u 2:(-$35/($6*(1.0+$6)*0.0101*0.0101/(1.0+0.0101))) w l t "b/bnm"
plot "stat.txt" u 2:($27/($7*(1.0-$7))) w l t "<y^2>/<y^2>nm", "stat.txt" u 2:($36/((0.0101**2)/((1.0+0.0101)**3)*$7*(1.0-$7))) w l t "<r^2>/<r^2>nm", "stat.txt" u 2:(-$37/($7*(1.0+$7)*0.0101*0.0101/(1.0+0.0101))) w l t "b/bnm" # almost mixed

# <y^2> / <y^2>nm = <r^2> / <r^2>nm = b / bnm, high At
plot "stat.txt" u 2:($23/($3*(1.0-$3))) w l t "<y^2>/<y^2>nm", "stat.txt" u 2:($28/((9.0**2)/((1.0+9.0)**3)*$3*(1.0-$3))) w l t "<r^2>/<r^2>nm", "stat.txt" u 2:(-$29/($3*(1.0+$3)*9.0*9.0/(1.0+9.0))) w l t "b/bnm" # no-mix
plot "stat.txt" u 2:($24/($4*(1.0-$4))) w l t "<y^2>/<y^2>nm", "stat.txt" u 2:($30/((9.0**2)/((1.0+9.0)**3)*$4*(1.0-$4))) w l t "<r^2>/<r^2>nm", "stat.txt" u 2:(-$31/($4*(1.0+$4)*9.0*9.0/(1.0+9.0))) w l t "b/bnm"
plot "stat.txt" u 2:($25/($5*(1.0-$5))) w l t "<y^2>/<y^2>nm", "stat.txt" u 2:($32/((9.0**2)/((1.0+9.0)**3)*$5*(1.0-$5))) w l t "<r^2>/<r^2>nm", "stat.txt" u 2:(-$33/($5*(1.0+$5)*9.0*9.0/(1.0+9.0))) w l t "b/bnm" # uniform
plot "stat.txt" u 2:($26/($6*(1.0-$6))) w l t "<y^2>/<y^2>nm", "stat.txt" u 2:($34/((9.0**2)/((1.0+9.0)**3)*$6*(1.0-$6))) w l t "<r^2>/<r^2>nm", "stat.txt" u 2:(-$35/($6*(1.0+$6)*9.0*9.0/(1.0+9.0))) w l t "b/bnm"
plot "stat.txt" u 2:($27/($7*(1.0-$7))) w l t "<y^2>/<y^2>nm", "stat.txt" u 2:($36/((9.0**2)/((1.0+9.0)**3)*$7*(1.0-$7))) w l t "<r^2>/<r^2>nm", "stat.txt" u 2:(-$37/($7*(1.0+$7)*9.0*9.0/(1.0+9.0))) w l t "b/bnm" # almost mixed

# consistency test: <Rv^2> = b<V> = b(1+b)/<R>, both At
plot "stat.txt" u 2:18 w l t "<Rv^2>", "stat.txt" u 2:(-$29*$13) w p ps 1.0 pt 7 t "b<V>", "stat.txt" u 2:(-$29*(1.0-$29)/$8) w p ps 1.0 pt 7 t "b(1+b)/<R>" # no-mix
plot "stat.txt" u 2:19 w l t "<Rv^2>", "stat.txt" u 2:(-$31*$14) w p ps 1.0 pt 7 t "b<V>", "stat.txt" u 2:(-$31*(1.0-$31)/$9) w p ps 1.0 pt 7 t "b(1+b)/<R>"
plot "stat.txt" u 2:20 w l t "<Rv^2>", "stat.txt" u 2:(-$33*$15) w p ps 1.0 pt 7 t "b<V>", "stat.txt" u 2:(-$33*(1.0-$33)/$10) w p ps 1.0 pt 7 t "b(1+b)/<R>" # uniform
plot "stat.txt" u 2:21 w l t "<Rv^2>", "stat.txt" u 2:(-$35*$16) w p ps 1.0 pt 7 t "b<V>", "stat.txt" u 2:(-$35*(1.0-$35)/$11) w p ps 1.0 pt 7 t "b(1+b)/<R>"
plot "stat.txt" u 2:22 w l t "<Rv^2>", "stat.txt" u 2:(-$37*$17) w p ps 1.0 pt 7 t "b<V>", "stat.txt" u 2:(-$37*(1.0-$37)/$12) w p ps 1.0 pt 7 t "b(1+b)/<R>" # almost mixed

# b = <R>^2<v^2> [1 + <v^2>/<v^2>nm * (rho1rho2/<R>^2 - 1) ] where <v^2>nm = (r/rho_2)^2 <y^2>nm, <y^2>nm = <Y>(1-<Y>), low At
plot "stat.txt" u 2:(-$29) w l t "b", "stat.txt" u 2:($8*$8*$38*(1.0+$38/0.0101/0.0101/$3/(1.0+$3)*(0.99/$8/$8-1.0))) w p ps 1.0 pt 7 t "<R>^2<v^2> [1 + <v^2>/<v^2>nm * (rho_1rho_2/<R>^2 - 1) ]" # no-mix
plot "stat.txt" u 2:(-$31) w l t "b", "stat.txt" u 2:($9*$9*$39*(1.0+$39/0.0101/0.0101/$4/(1.0+$4)*(0.99/$9/$9-1.0))) w p ps 1.0 pt 7 t "<R>^2<v^2> [1 + <v^2>/<v^2>nm * (rho1_rho2/<R>^2 - 1) ]"
plot "stat.txt" u 2:(-$33) w l t "b", "stat.txt" u 2:($10*$10*$40*(1.0+$40/0.0101/0.0101/$5/(1.0+$5)*(0.99/$10/$10-1.0))) w p ps 1.0 pt 7 t "<R>^2<v^2> [1 + <v^2>/<v^2>nm * (rho_1rho2/<R>^2 - 1) ]" # uniform
plot "stat.txt" u 2:(-$35) w l t "b", "stat.txt" u 2:($11*$11*$41*(1.0+$41/0.0101/0.0101/$6/(1.0+$6)*(0.99/$11/$11-1.0))) w p ps 1.0 pt 7 t "<R>^2<v^2> [1 + <v^2>/<v^2>nm * (rho_1rho_2/<R>^2 - 1) ]"
plot "stat.txt" u 2:(-$37) w l t "b", "stat.txt" u 2:($12*$12*$42*(1.0+$42/0.0101/0.0101/$7/(1.0+$7)*(0.99/$12/$12-1.0))) w p ps 1.0 pt 7 t "<R>^2<v^2> [1 + <v^2>/<v^2>nm * (rho_1rho2/<R>^2 - 1) ]" # almost mixed

# b = <R>^2<v^2> [1 + <v^2>/<v^2>nm * (rho1rho2/<R>^2 - 1) ] where <v^2>nm = (r/rho_2)^2 <y^2>nm, <y^2>nm = <Y>(1-<Y>), high At
plot "stat.txt" u 2:(-$29) w l t "b", "stat.txt" u 2:($8*$8*$38*(1.0+$38/9.0/9.0/$3/(1.0+$3)*(0.1/$8/$8-1.0))) w p ps 1.0 pt 7 t "<R>^2<v^2> [1 + <v^2>/<v^2>nm * (rho_1rho_2/<R>^2 - 1) ]" # no-mix
plot "stat.txt" u 2:(-$31) w l t "b", "stat.txt" u 2:($9*$9*$39*(1.0+$39/9.0/9.0/$4/(1.0+$4)*(0.1/$9/$9-1.0))) w p ps 1.0 pt 7 t "<R>^2<v^2> [1 + <v^2>/<v^2>nm * (rho_1rho_2/<R>^2 - 1) ]"
plot "stat.txt" u 2:(-$33) w l t "b", "stat.txt" u 2:($10*$10*$40*(1.0+$40/9.0/9.0/$5/(1.0+$5)*(0.1/$10/$10-1.0))) w p ps 1.0 pt 7 t "<R>^2<v^2> [1 + <v^2>/<v^2>nm * (rho_1rho_2/<R>^2 - 1) ]" # uniform
plot "stat.txt" u 2:(-$35) w l t "b", "stat.txt" u 2:($11*$11*$41*(1.0+$41/9.0/9.0/$6/(1.0+$6)*(0.1/$11/$11-1.0))) w p ps 1.0 pt 7 t "<R>^2<v^2> [1 + <v^2>/<v^2>nm * (rho_1rho_2/<R>^2 - 1) ]"
plot "stat.txt" u 2:(-$37) w l t "b", "stat.txt" u 2:($12*$12*$42*(1.0+$42/9.0/9.0/$7/(1.0+$7)*(0.1/$12/$12-1.0))) w p ps 1.0 pt 7 t "<R>^2<v^2> [1 + <v^2>/<v^2>nm * (rho_1rho_2/<R>^2 - 1) ]" # almost mixed

# <r^2> = <R>^4<v^2> [1 + <v^2>/<v^2>nm * (rho1^3rho2/<R>^4 - 1) ] where <v^2>nm = (r/rho_2)^2 <y^2>nm, <y^2>nm = <Y>(1-<Y>), low At
plot "stat.txt" u 2:(-$29) w l t "b", "stat.txt" u 2:($8**4.0*$38*(1.0+$38/0.0101/0.0101/$3/(1.0+$3)*(0.99**2.0/$8**4.0-1.0))) w p ps 1.0 pt 7 t "<R>^4<v^2> [1 + <v^2>/<v^2>nm * (rho_1^3rho_2/<R>^4 - 1) ]" # no-mix
plot "stat.txt" u 2:(-$31) w l t "b", "stat.txt" u 2:($9**4.0*$39*(1.0+$39/0.0101/0.0101/$4/(1.0+$4)*(0.99**2.0/$9**4.0-1.0))) w p ps 1.0 pt 7 t "<R>^4<v^2> [1 + <v^2>/<v^2>nm * (rho_1^3rho_2/<R>^4 - 1) ]"
plot "stat.txt" u 2:(-$33) w l t "b", "stat.txt" u 2:($10**4.0*$40*(1.0+$40/0.0101/0.0101/$5/(1.0+$5)*(0.99**2.0/$10**4.0-1.0))) w p ps 1.0 pt 7 t "<R>^4<v^2> [1 + <v^2>/<v^2>nm * (rho_1^3rho_2/<R>^4 - 1) ]" # uniform
plot "stat.txt" u 2:(-$35) w l t "b", "stat.txt" u 2:($11**4.0*$41*(1.0+$41/0.0101/0.0101/$6/(1.0+$6)*(0.99**2.0/$11**4.0-1.0))) w p ps 1.0 pt 7 t "<R>^4<v^2> [1 + <v^2>/<v^2>nm * (rho_1^3rho_2/<R>^4 - 1) ]"
plot "stat.txt" u 2:(-$37) w l t "b", "stat.txt" u 2:($12**4.0*$42*(1.0+$42/0.0101/0.0101/$7/(1.0+$7)*(0.99**2.0/$12**4.0-1.0))) w p ps 1.0 pt 7 t "<R>^4<v^2> [1 + <v^2>/<v^2>nm * (rho_1^3rho_2/<R>^4 - 1) ]" # almost mixed

# <r^2> = <R>^4<v^2> [1 + <v^2>/<v^2>nm * (rho1^3rho2/<R>^4 - 1) ] where <v^2>nm = (r/rho_2)^2 <y^2>nm, <y^2>nm = <Y>(1-<Y>), high At
plot "stat.txt" u 2:(-$29) w l t "b", "stat.txt" u 2:($8**4.0*$38*(1.0+$38/9.0/9.0/$3/(1.0+$3)*(0.1**2.0/$8**4.0-1.0))) w p ps 1.0 pt 7 t "<R>^4<v^2> [1 + <v^2>/<v^2>nm * (rho_1^2rho_2^2/<R>^4 - 1) ]" # no-mix
plot "stat.txt" u 2:(-$31) w l t "b", "stat.txt" u 2:($9**4.0*$39*(1.0+$39/9.0/9.0/$4/(1.0+$4)*(0.1**2.0/$9**4.0-1.0))) w p ps 1.0 pt 7 t "<R>^4<v^2> [1 + <v^2>/<v^2>nm * (rho_1^2rho_2^2/<R>^4 - 1) ]"
plot "stat.txt" u 2:(-$33) w l t "b", "stat.txt" u 2:($10**4.0*$40*(1.0+$40/9.0/9.0/$5/(1.0+$5)*(0.1**2.0/$10**4.0-1.0))) w p ps 1.0 pt 7 t "<R>^4<v^2> [1 + <v^2>/<v^2>nm * (rho_1^2rho_2^2/<R>^4 - 1) ]" # uniform
plot "stat.txt" u 2:(-$35) w l t "b", "stat.txt" u 2:($11**4.0*$41*(1.0+$41/9.0/9.0/$6/(1.0+$6)*(0.1**2.0/$11**4.0-1.0))) w p ps 1.0 pt 7 t "<R>^4<v^2> [1 + <v^2>/<v^2>nm * (rho_1^2rho_2^2/<R>^4 - 1) ]"
plot "stat.txt" u 2:(-$37) w l t "b", "stat.txt" u 2:($12**4.0*$42*(1.0+$42/9.0/9.0/$7/(1.0+$7)*(0.1**2.0/$12**4.0-1.0))) w p ps 1.0 pt 7 t "<R>^4<v^2> [1 + <v^2>/<v^2>nm * (rho_1^2rho_2^2/<R>^4 - 1) ]" # almost mixed