Walker: Integrating the mix mass-fraction beta SDE
This example runs Walker to integrate the mix mass-fraction beta SDE (see DiffEq/MixMassFractionBeta.h) using constant coefficients. The mix mass-fraction beta SDE is based on the beta SDE, (1) additionally computing two more stochastic variables that are functions of the beta variables integrated, and (2) the coefficients b and kappa are specified via functions that constrain the beta SDE to be consistent with the turbulent mixing process. For more detail on the beta SDE, see https:/
Control file
# vim: filetype=sh: title "Test evolution of some stats in mass fraction mixing" walker #nstep 1 # Max number of time steps term 5.0 # Max time dt 0.01 # Time step size npar 10000 # Number of particles ttyi 10 # TTY output interval rngs mkl_r250 beta_method cja end end mixmassfracbeta depvar y ncomp 20 init jointbeta icbeta betapdf 0.01 0.01 0.0 1.0 end # light = heavy betapdf 0.2 0.8 0.0 1.0 end # light < heavy betapdf 0.8 0.2 0.0 1.0 end # light > heavy betapdf 0.116517612 0.2 0.0 1.0 end # <r^3> approximately 0 betapdf 2.0 5.0 0.0 1.0 end end # init jointdelta # icdelta # spike 0.01 0.5 0.99 0.5 end # spike 0.01 0.9 0.99 0.1 end # spike 0.01 0.1 0.99 0.9 end # spike 0.01 0.5 0.99 0.5 end # spike 0.01 0.5 0.99 0.5 end # end coeff homdecay # alpha = Sb/kappa, beta = (1-S)b/kappa # decay in <y^2> if bprime/kappaprime > 1/4 kappaprime 1.0 1.0 1.0 1.0 1.0 end bprime 1.9 1.9 1.9 1.9 1.9 end S 0.5 0.5 0.5 0.5 0.5 end rng mkl_r250 rho2 1.0 1.0 1.0 1.0 1.0 end #r 2.0 2.0 2.0 2.0 2.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 #r 4.0 4.0 4.0 4.0 4.0 end end statistics format scientific precision 12 # <Y>, mass fraction means <Y1> <Y2> <Y3> <Y4> <Y5> # <R>, mean densities <Y6> <Y7> <Y8> <Y9> <Y10> # <V>, mean specific volumes <Y11> <Y12> <Y13> <Y14> <Y15> # <y^2>, mass fraction variances <y1y1> <y2y2> <y3y3> <y4y4> <y5y5> # stats required for homdecay # <y^3>, mass fraction third moments <y1y1y1> <y2y2y2> <y3y3y3> <y4y4y4> <y5y5y5> # <r^2>, density variances <y6y6> <y7y7> <y8y8> <y9y9> <y10y10> # <r^3>, density third moments <y6y6y6> <y7y7y7> <y8y8y8> <y9y9y9> <y10y10y10> # <rv>, density-specific-volume covariances <y6y11> <y7y12> <y8y13> <y9y14> <y10y15> # <RY> <Y6Y1> <Y7Y2> <Y8Y3> <Y9Y4> <Y10Y5> # <Rv^2> <Y6y11y11> <Y7y12y12> <Y8y13y13> <Y9y14y14> <Y10y15y15> # <rv^2> <y6y11y11> <y7y12y12> <y8y13y13> <y9y14y14> <y10y15y15> # <v^2>, specific volume variances <y11y11> <y12y12> <y13y13> <y14y14> <y15y15> # # additional stats required for mchomdecay # # <R^2> # <Y6Y6> <Y7Y7> <Y8Y8> <Y9Y9> <Y10Y10> # # <YR^2> # <Y1Y6Y6> <Y2Y7Y7> <Y3Y8Y8> <Y4Y9Y9> <Y5Y10Y10> # # <Y(1-Y)R^3> # <Y1Y16Y6Y6Y6> # <Y2Y17Y7Y7Y7> # <Y3Y18Y8Y8Y8> # <Y4Y19Y9Y9Y9> # <Y5Y20Y10Y10Y10> end pdfs interval 1 filetype txt policy overwrite centering elem format scientific precision 4 p1( Y1 : 5.0e-3 ) p2( Y2 : 5.0e-3 ) p3( Y3 : 5.0e-3 ) p4( Y4 : 1.0e-2 ) # p5( Y5 : 1.0e-2 ) p6( Y6 : 5.0e-5 ) p7( Y7 : 5.0e-5 ) p8( Y8 : 5.0e-5 ) p9( Y9 : 1.0e-2 ) # p10( Y10 : 1.0e-2 ) end end
Example run on 8 CPUs
./charmrun +p8 Main/walker -v -c ../../tmp/mixmassfracbeta.q -v
Results
The rationale for these runs is to integrate the system in time, extract the time evolution of various statistics and PDFs, and then test if the evolution of the statistics allow us to qualitatively correctly model the physical process of mixing. Example questions these simulations help answer:
- Can we start with a symmetric initial mass fraction PDF and generate a nonzero skewness?
- Can we start with positive or negative initial mass fraction skewness?
- Does the time evolution ensure monotonic decay of the mass fraction variancess as well as the density variances at all times?
- Does the initially symmetric mass fraction PDF generate a smaller skewness for the low Atwood number than the large Atwood number case?
# 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 # For homdecay # <Y>, mass fraction means plot "stat.txt" u 2:3 w l t "<Y1>", "stat.txt" u 2:4 w l t "<Y2>", "stat.txt" u 2:5 w l t "<Y3>", "stat.txt" u 2:6 w l t "<Y4>", "stat.txt" u 2:7 w l t "<Y5>" # <R>, mean densities plot "stat.txt" u 2:8 w l t "<R1>", "stat.txt" u 2:10 w l t "<R2>", "stat.txt" u 2:12 w l t "<R3>", "stat.txt" u 2:14 w l t "<R4>", "stat.txt" u 2:16 w l t "<R5>" # <V>, mean specific volumes plot "stat.txt" u 2:18 w l t "<V1>", "stat.txt" u 2:19 w l t "<V2>", "stat.txt" u 2:20 w l t "<V3>", "stat.txt" u 2:21 w l t "<V4>", "stat.txt" u 2:22 w l t "<V5>" # <y^2>, mass fraction variances plot "stat.txt" u 2:28 w l t "<y1y1>", "stat.txt" u 2:30 w l t "<y2y2>", "stat.txt" u 2:32 w l t "<y3y3>", "stat.txt" u 2:34 w l t "<y4y4>", "stat.txt" u 2:36 w l t "<y5y5>" # <y^3>, mass fraction third moments plot "stat.txt" u 2:29 w l t "<y1y1y1>", "stat.txt" u 2:31 w l t "<y2y2y2>", "stat.txt" u 2:33 w l t "<y3y3y3>", "stat.txt" u 2:35 w l t "<y4y4y4>", "stat.txt" u 2:37 w l t "<y5y5y5>" # <r^2>, density variances plot "stat.txt" u 2:38 w l t "<r1r1>", "stat.txt" u 2:42 w l t "<r2r2>", "stat.txt" u 2:46 w l t "<r3r3>", "stat.txt" u 2:50 w l t "<r4r4>", "stat.txt" u 2:54 w l t "<r5r5>" # <r^3>, density third moments plot "stat.txt" u 2:39 w l t "<r1r1r1>", "stat.txt" u 2:43 w l t "<r2r2r2>", "stat.txt" u 2:47 w l t "<r3r3r3>", "stat.txt" u 2:51 w l t "<r4r4r4>", "stat.txt" u 2:55 w l t "<r5r5r5>" # <rv>, density-specific-volume covariances plot "stat.txt" u 2:40 w l t "<r1v1>", "stat.txt" u 2:44 w l t "<r2v2>", "stat.txt" u 2:48 w l t "<r3v3>", "stat.txt" u 2:52 w l t "<r4v4>", "stat.txt" u 2:56 w l t "<r5v5>" # <v^2>, specific-volume variances plot "stat.txt" u 2:58 w l t "<v1v1>", "stat.txt" u 2:59 w l t "<v2v2>", "stat.txt" u 2:60 w l t "<v3v3>", "stat.txt" u 2:61 w l t "<v4v4>", "stat.txt" u 2:62 w l t "<v5v5>" # <RY> plot "stat.txt" u 2:9 w l t "<R1Y1>", "stat.txt" u 2:11 w l t "<R2Y2>", "stat.txt" u 2:13 w l t "<R3Y3>", "stat.txt" u 2:15 w l t "<R4Y4>", "stat.txt" u 2:17 w l t "<R5Y5>" # <RY>/<R> plot "stat.txt" u 2:($9/$8) w l t "<R1Y1>/<R1>", "stat.txt" u 2:($11/$10) w l t "<R2Y2>/<R2>", "stat.txt" u 2:($13/$12) w l t "<R3Y3>/<R3>", "stat.txt" u 2:($15/$14) w l t "<R4Y4>/<R4>", "stat.txt" u 2:($17/$16) w l t "<R5Y5>/<R5>" # Normalized <y^2> = <y^2> / ( <Y> (1-<Y>) ) plot "stat.txt" u 2:($28/$3/(1.0-$3)) w l t "<y1y1>/<Y>/(1-<Y>), light = heavy", "" u 2:($30/$4/(1.0-$4)) w l t "<y1y1>/<Y>/(1-<Y>), light > heavy", "" u 2:($32/$5/(1.0-$5)) w l t "<y1y1>/<Y>/(1-<Y>), light < heavy", "" u 2:($34/$6/(1.0-$6)) w l t "<y1y1>/<Y>/(1-<Y>), <r^3>(t=0) = 0" # Normalized <y^2> = <y^2> / ( S (1-S) ) S = 0.5; plot "stat.txt" u 2:($28/S/(1.0-S)) w l t "<y1y1>/<Y>/(1-<Y>), light = heavy", "" u 2:($30/S/(1.0-S)) w l t "<y1y1>/<Y>/(1-<Y>), light > heavy", "" u 2:($32/S/(1.0-S)) w l t "<y1y1>/<Y>/(1-<Y>), light < heavy", "" u 2:($34/S/(1.0-S)) w l t "<y1y1>/<Y>/(1-<Y>), <r^3>(t=0) = 0" y1 = 4.961030229954e-01; y2 = 1.987758126229e-01; y3 = 7.939404938833e-01; y4 = 3.638634671405e-01; plot "stat.txt" u 2:($28/y1/(1.0-y1)) w l t "<y1y1>/<Y>/(1-<Y>), light = heavy", "" u 2:($30/y2/(1.0-y2)) w l t "<y1y1>/<Y>/(1-<Y>), light > heavy", "" u 2:($32/y3/(1.0-y3)) w l t "<y1y1>/<Y>/(1-<Y>), light < heavy", "" u 2:($34/y4/(1.0-y4)) w l t "<y1y1>/<Y>/(1-<Y>), <r^3>(t=0) = 0" # b(1+r)/r/r r = 9.0; plot "stat.txt" u 2:(-$40*(1.0+r)/r/r) w l t "b(1+r)/r^2", "" u 2:(-$44*(1.0+r)/r/r) w l t "b(1+r)/r^2", "" u 2:(-$48*(1.0+r)/r/r) w l t "b(1+r)/r^2", "" u 2:(-$52*(1.0+r)/r/r) w l t "b(1+r)/r^2" # <v^2> ( rho2 / r )^2 r = 9.0; plot "stat.txt" u 2:($58*(1.0/r)**2.0) w l t "<v^2>(rho2/r)^2", "" u 2:($59*(1.0/r)**2.0) w l t "<v^2>(rho2/r)^2", "" u 2:($60*(1.0/r)**2.0) w l t "<v^2>(rho2/r)^2", "" u 2:($61*(1.0/r)**2.0) w l t "<v^2>(rho2/r)^2" # <r^2> (1+r)^2 / (rho2 r)^2 r = 9.0; plot "stat.txt" u 2:($38*(1.0+r)**2.0/r/r) w l t "<rho^2>(1+r)^2/(rho2 r)^2", "" u 2:($42*(1.0+r)**2.0/r/r) w l t "<rho^2>(1+r)^2/(rho2 r)^2", "" u 2:($46*(1.0+r)**2.0/r/r) w l t "<rho^2>(1+r)^2/(rho2 r)^2", "" u 2:($50*(1.0+r)**2.0/r/r) w l t "<rho^2>(1+r)^2/(rho2 r)^2" # Test closure: 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:(-$40) w l t "b", "stat.txt" u 2:($8*$8*$58*(1.0+$58/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) ]" plot "stat.txt" u 2:(-$44) w l t "b", "stat.txt" u 2:($9*$9*$59*(1.0+$59/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 * (rho1_rho2/<R>^2 - 1) ]" plot "stat.txt" u 2:(-$48) w l t "b", "stat.txt" u 2:($10*$10*$60*(1.0+$60/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_1rho2/<R>^2 - 1) ]" plot "stat.txt" u 2:(-$52) w l t "b", "stat.txt" u 2:($11*$11*$61*(1.0+$61/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:(-$56) w l t "b", "stat.txt" u 2:($12*$12*$62*(1.0+$62/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) ]" # consistency test: <Rv^2> = b<V> = b(1+b)/<R> plot "stat.txt" u 2:23 w l t "<Rv^2>", "stat.txt" u 2:(-$40*$18) w p ps 1.0 pt 7 t "b<V>", "stat.txt" u 2:(-$40*(1.0-$40)/$8) w p ps 1.0 pt 6 t "b(1+b)/<R>" plot "stat.txt" u 2:24 w l t "<Rv^2>", "stat.txt" u 2:(-$44*$19) w p ps 1.0 pt 7 t "b<V>", "stat.txt" u 2:(-$44*(1.0-$44)/$10) w p ps 1.0 pt 6 t "b(1+b)/<R>" plot "stat.txt" u 2:25 w l t "<Rv^2>", "stat.txt" u 2:(-$48*$20) w p ps 1.0 pt 7 t "b<V>", "stat.txt" u 2:(-$48*(1.0-$48)/$12) w p ps 1.0 pt 6 t "b(1+b)/<R>" plot "stat.txt" u 2:26 w l t "<Rv^2>", "stat.txt" u 2:(-$52*$21) w p ps 1.0 pt 7 t "b<V>", "stat.txt" u 2:(-$52*(1.0-$52)/$14) w p ps 1.0 pt 6 t "b(1+b)/<R>" plot "stat.txt" u 2:27 w l t "<Rv^2>", "stat.txt" u 2:(-$56*$22) w p ps 1.0 pt 7 t "b<V>", "stat.txt" u 2:(-$56*(1.0-$56)/$16) w p ps 1.0 pt 6 t "b(1+b)/<R>" # test relative magnitude of <rv^2> compared to <R><v^2>, plot ratio plot "stat.txt" u 2:($41/$8/$58) w l t "<rv^2>/<R>/<v^2> light=heavy" plot "stat.txt" u 2:($45/$10/$59) w l t "<rv^2>/<R>/<v^2> light > heavy" plot "stat.txt" u 2:($49/$12/$60) w l t "<rv^2>/<R>/<v^2> light < heavy" plot "stat.txt" u 2:($53/$14/$61) w l t "<rv^2>/<R>/<v^2> <r^3>(t=0) = 0" plot "stat.txt" u 2:($57/$16/$62) w l t "<rv^2>/<R>/<v^2>" # first four in one plot plot "stat.txt" u 2:($41/$8/$58) w l t "<rv^2>/<R>/<v^2> light=heavy", "stat.txt" u 2:($45/$10/$59) w l t "<rv^2>/<R>/<v^2> light > heavy", "stat.txt" u 2:($49/$12/$60) w l t "<rv^2>/<R>/<v^2> light < heavy", "stat.txt" u 2:($53/$14/$61) w l t "<rv^2>/<R>/<v^2> <r^3>(t=0) = 0" # test closure model for <rv^2> = -r^3/rho_2 * 1/(1+r) * <y^2>/(1+r<RY>/<R>)^2 * [ 1 - 2<RY>/<R> - r(2+r)(<RY>/<R>)^2 ] =approx= <rv^2>nm, r = 9.0 -> high A, r = 0.0101 -> low A r = 9.0; plot "stat.txt" u 2:41 w l t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$28/(1.0+r*$9/$8)/(1.0+r*$9/$8)*(1.0-2.0*$9/$8-r*(2.0+r)*$9/$8*$9/$8)) w l t "<rv^2>nm light=heavy" r = 9.0; plot "stat.txt" u 2:45 w l t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$30/(1.0+r*$11/$10)/(1.0+r*$11/$10)*(1.0-2.0*$11/$10-r*(2.0+r)*$11/$10*$11/$10)) w l t "<rv^2>nm light > heavy" r = 9.0; plot "stat.txt" u 2:49 w l t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$32/(1.0+r*$13/$12)/(1.0+r*$13/$12)*(1.0-2.0*$13/$12-r*(2.0+r)*$13/$12*$13/$12)) w l t "<rv^2>nm light < heavy" r = 9.0; plot "stat.txt" u 2:53 w l t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$34/(1.0+r*$15/$14)/(1.0+r*$15/$14)*(1.0-2.0*$15/$14-r*(2.0+r)*$15/$14*$15/$14)) w l t "<rv^2>nm <r^3>(t=0) = 0" # all in on plot r = 9.0; plot "stat.txt" u 2:41 w l t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$28/(1.0+r*$9/$8)/(1.0+r*$9/$8)*(1.0-2.0*$9/$8-r*(2.0+r)*$9/$8*$9/$8)) w p pt 7 lt 1 t "<rv^2>nm light=heavy", "" u 2:45 w l lt 2 t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$30/(1.0+r*$11/$10)/(1.0+r*$11/$10)*(1.0-2.0*$11/$10-r*(2.0+r)*$11/$10*$11/$10)) w p pt 7 lt 2 t "<rv^2>nm light > heavy", "" u 2:49 w l lt 3 t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$32/(1.0+r*$13/$12)/(1.0+r*$13/$12)*(1.0-2.0*$13/$12-r*(2.0+r)*$13/$12*$13/$12)) w p pt 7 lt 3 t "<rv^2>nm light < heavy", "" u 2:53 w l lt 4 t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$34/(1.0+r*$15/$14)/(1.0+r*$15/$14)*(1.0-2.0*$15/$14-r*(2.0+r)*$15/$14*$15/$14)) w p pt 7 lt 4 t "<rv^2>nm <r^3>(t=0) = 0" # test closure model for <rv^2> = -r^3/rho_2 * 1/(1+r) * <y^2>/(1+r<RY>/<R>)^1 * [ 1 - 2<RY>/<R> - r^2(<RY>/<R>)^2 - b*(1+r<RY>/<R>)^2/r ], r = 9.0 -> high A, r = 0.0101 -> low A r = 9.0; plot "stat.txt" u 2:41 w l t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$28/(1.0+r*$9/$8)*(1.0-2.0*$9/$8-r**2.0*$9/$8*$9/$8+$40*(1.0+r*$9/$8)**2.0/r)) w l t "<rv^2>nm light=heavy" r = 9.0; plot "stat.txt" u 2:45 w l t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$30/(1.0+r*$11/$10)*(1.0-2.0*$11/$10-r**2.0*$11/$10*$11/$10+$44*(1.0+r*$11/$10)**2.0/r)) w l t "<rv^2>nm light > heavy" r = 9.0; plot "stat.txt" u 2:49 w l t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$32/(1.0+r*$13/$12)*(1.0-2.0*$13/$12-r**2.0*$13/$12*$13/$12+$48*(1.0+r*$13/$12)**2.0/r)) w l t "<rv^2>nm light < heavy" r = 9.0; plot "stat.txt" u 2:53 w l t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$34/(1.0+r*$15/$14)*(1.0-2.0*$15/$14-r**2.0*$15/$14*$15/$14+$52*(1.0+r*$15/$14)**2.0/r)) w l t "<rv^2>nm <r^3>(t=0) = 0" # all in on plot r = 9.0; plot "stat.txt" u 2:41 w l t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$28/(1.0+r*$9/$8)*(1.0-2.0*$9/$8-r**2.0*$9/$8*$9/$8+$40*(1.0+r*$9/$8)**2.0/r)) w p pt 7 lt 1 t "<rv^2>nm light=heavy", "" u 2:45 w l lt 2 t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$30/(1.0+r*$11/$10)*(1.0-2.0*$11/$10-r**2.0*$11/$10*$11/$10+$44*(1.0+r*$11/$10)**2.0/r)) w p pt 7 lt 2 t "<rv^2>nm light > heavy", "" u 2:49 w l lt 3 t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$32/(1.0+r*$13/$12)*(1.0-2.0*$13/$12-r**2.0*$13/$12*$13/$12+$48*(1.0+r*$13/$12)**2.0/r)) w p pt 7 lt 3 t "<rv^2>nm light < heavy", "" u 2:53 w l lt 4 t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$34/(1.0+r*$15/$14)*(1.0-2.0*$15/$14-r**2.0*$15/$14*$15/$14+$52*(1.0+r*$15/$14)**2.0/r)) w p pt 7 lt 4 t "<rv^2>nm <r^3>(t=0) = 0" # test closure model for <rv^2> = -r^3/rho_2 * 1/(1+r) * <y^2>/(1+r<RY>/<R>)^2 * [ 1 - 2<RY>/<R> - r(<RY>/<R>)^2 - b*(1+r<RY>/<R>)/r ], r = 9.0 -> high A, r = 0.0101 -> low A r = 9.0; plot "stat.txt" u 2:41 w l t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$28/(1.0+r*$9/$8)**2.0*(1.0-2.0*$9/$8-r*$9/$8*$9/$8+$40*(1.0+r*$9/$8)/r)) w l t "<rv^2>nm light=heavy" r = 9.0; plot "stat.txt" u 2:45 w l t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$30/(1.0+r*$11/$10)**2.0*(1.0-2.0*$11/$10-r*$11/$10*$11/$10+$44*(1.0+r*$11/$10)/r)) w l t "<rv^2>nm light > heavy" r = 9.0; plot "stat.txt" u 2:49 w l t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$32/(1.0+r*$13/$12)**2.0*(1.0-2.0*$13/$12-r*$13/$12*$13/$12+$48*(1.0+r*$13/$12)/r)) w l t "<rv^2>nm light < heavy" r = 9.0; plot "stat.txt" u 2:53 w l t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$34/(1.0+r*$15/$14)**2.0*(1.0-2.0*$15/$14-r*$15/$14*$15/$14+$52*(1.0+r*$15/$14)/r)) w l t "<rv^2>nm <r^3>(t=0) = 0" # all in on plot r = 9.0; plot "stat.txt" u 2:41 w l t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$28/(1.0+r*$9/$8)**2.0*(1.0-2.0*$9/$8-r*$9/$8*$9/$8+$40*(1.0+r*$9/$8)/r)) w p pt 7 lt 1 t "<rv^2>nm light=heavy", "" u 2:45 w l lt 2 t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$30/(1.0+r*$11/$10)**2.0*(1.0-2.0*$11/$10-r*$11/$10*$11/$10+$44*(1.0+r*$11/$10)/r)) w p pt 7 lt 2 t "<rv^2>nm light > heavy", "" u 2:49 w l lt 3 t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$32/(1.0+r*$13/$12)**2.0*(1.0-2.0*$13/$12-r*$13/$12*$13/$12+$48*(1.0+r*$13/$12)/r)) w p pt 7 lt 3 t "<rv^2>nm light < heavy", "" u 2:53 w l lt 4 t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$34/(1.0+r*$15/$14)**2.0*(1.0-2.0*$15/$14-r*$15/$14*$15/$14+$52*(1.0+r*$15/$14)/r)) w p pt 7 lt 4 t "<rv^2>nm <r^3>(t=0) = 0" ##### BASE MODEL v1 ##### # test closure model for <rv^2> = -r^3/rho_2 * 1/(1+r) * <y^2>/(1+r<RY>/<R>)^2 * [ 1 - 2<RY>/<R> - r(<RY>/<R>)^2 - b*(1+r<RY>/<R>)^2/r ], r = 9.0 -> high A, r = 0.0101 -> low A ##### BASE MODEL v1 ##### r = 9.0; plot "stat.txt" u 2:41 w l t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$28/(1.0+r*$9/$8)**2.0*(1.0-2.0*$9/$8-r*$9/$8*$9/$8+$40*(1.0+r*$9/$8)**2.0/r)) w l t "<rv^2>nm light=heavy" r = 9.0; plot "stat.txt" u 2:45 w l t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$30/(1.0+r*$11/$10)**2.0*(1.0-2.0*$11/$10-r*$11/$10*$11/$10+$44*(1.0+r*$11/$10)**2.0/r)) w l t "<rv^2>nm light > heavy" r = 9.0; plot "stat.txt" u 2:49 w l t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$32/(1.0+r*$13/$12)**2.0*(1.0-2.0*$13/$12-r*$13/$12*$13/$12+$48*(1.0+r*$13/$12)**2.0/r)) w l t "<rv^2>nm light < heavy" r = 9.0; plot "stat.txt" u 2:53 w l t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$34/(1.0+r*$15/$14)**2.0*(1.0-2.0*$15/$14-r*$15/$14*$15/$14+$52*(1.0+r*$15/$14)**2.0/r)) w l t "<rv^2>nm <r^3>(t=0) = 0" # all in on plot r = 9.0; plot "stat.txt" u 2:41 w l t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$28/(1.0+r*$9/$8)**2.0*(1.0-2.0*$9/$8-r*$9/$8*$9/$8+$40*(1.0+r*$9/$8)**2.0/r)) w p pt 7 lt 1 t "<rv^2>nm light=heavy", "" u 2:45 w l lt 2 t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$30/(1.0+r*$11/$10)**2.0*(1.0-2.0*$11/$10-r*$11/$10*$11/$10+$44*(1.0+r*$11/$10)**2.0/r)) w p pt 7 lt 2 t "<rv^2>nm light > heavy", "" u 2:49 w l lt 3 t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$32/(1.0+r*$13/$12)**2.0*(1.0-2.0*$13/$12-r*$13/$12*$13/$12+$48*(1.0+r*$13/$12)**2.0/r)) w p pt 7 lt 3 t "<rv^2>nm light < heavy", "" u 2:53 w l lt 4 t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$34/(1.0+r*$15/$14)**2.0*(1.0-2.0*$15/$14-r*$15/$14*$15/$14+$52*(1.0+r*$15/$14)**2.0/r)) w p pt 7 lt 4 t "<rv^2>nm <r^3>(t=0) = 0" ##### BASE MODEL v1 ##### # test closure model for <y^2> = b(1+b) / ( <R>^2 (r/rho2)^2 { 1 - r/rho2 <R>/(1+r) [ 1 - 2<RY>/<R> - r(<RY>/<R>)^2 - b/r(1+r<RY>/<R>)^2 ] } ) ##### BASE MODEL v1 ##### r = 9.0; plot "stat.txt" u 2:28 w l t "<y^2> MC", "" u 2:(-$40*(1.0-$40)/($8**2.0*r**2.0*(1.0-r*$8/(1.0+r)*(1.0-2.0*$9/$8-r*($9/$8)**2.0+$40/r*(1.0+r*$9/$8)**2.0)))) w l t "<y^2> model light=heavy" r = 9.0; plot "stat.txt" u 2:30 w l t "<y^2> MC", "" u 2:(-$44*(1.0-$44)/($10**2.0*r**2.0*(1.0-r*$10/(1.0+r)*(1.0-2.0*$11/$10-r*($11/$10)**2.0+$44/r*(1.0+r*$11/$10)**2.0)))) w l t "<y^2> model light > heavy" r = 9.0; plot "stat.txt" u 2:32 w l t "<y^2> MC", "" u 2:(-$48*(1.0-$48)/($12**2.0*r**2.0*(1.0-r*$12/(1.0+r)*(1.0-2.0*$13/$12-r*($13/$12)**2.0+$48/r*(1.0+r*$13/$12)**2.0)))) w l t "<y^2> model light < heavy" r = 9.0; plot "stat.txt" u 2:34 w l t "<y^2> MC", "" u 2:(-$52*(1.0-$52)/($14**2.0*r**2.0*(1.0-r*$14/(1.0+r)*(1.0-2.0*$15/$14-r*($15/$14)**2.0+$52/r*(1.0+r*$15/$14)**2.0)))) w l t "<y^2> model <r^3>(t=0) = 0" # all in on plot r = 9.0; plot "stat.txt" u 2:28 w l t "<y^2> MC", "" u 2:(-$40*(1.0-$40)/($8**2.0*r**2.0*(1.0-r*$8/(1.0+r)*(1.0-2.0*$9/$8-r*($9/$8)**2.0+$40/r*(1.0+r*$9/$8)**2.0)))) w p pt 7 lt 1 t "<y^2> model light=heavy", "" u 2:30 w l lt 2 t "<y^2> MC", "" u 2:(-$44*(1.0-$44)/($10**2.0*r**2.0*(1.0-r*$10/(1.0+r)*(1.0-2.0*$11/$10-r*($11/$10)**2.0+$44/r*(1.0+r*$11/$10)**2.0)))) w p pt 7 lt 2 t "<y^2> model light > heavy", "" u 2:32 w l lt 3 t "<y^2> MC", "" u 2:(-$48*(1.0-$48)/($12**2.0*r**2.0*(1.0-r*$12/(1.0+r)*(1.0-2.0*$13/$12-r*($13/$12)**2.0+$48/r*(1.0+r*$13/$12)**2.0)))) w p pt 7 lt 3 t "<y^2> model light < heavy", "" u 2:34 w l lt 4 t "<y^2> MC", "" u 2:(-$52*(1.0-$52)/($14**2.0*r**2.0*(1.0-r*$14/(1.0+r)*(1.0-2.0*$15/$14-r*($15/$14)**2.0+$52/r*(1.0+r*$15/$14)**2.0)))) w p pt 7 lt 4 t "<y^2> model <r^3>(t=0) = 0" # test linear closure model for <y^2> = (1-thetab)(1+r)<RY>/<R>(1-<RY>/<R>)/(1+r<RY>/<R>)^2 + thetab b(1+b)(1+r<RY>/<R>)^2/r^2, with thetab = 1.0-b/bnm, with bnm = r^2 <RY>/<R> (1-<RY>/<R>) / (1+<RY>/<R>)^2 # thetab = 1.0+$40/(-r*r*$9/$8*(1.0-$9/$8)/(1.0+r*$9/$8)**2.0) r = 9.0; plot "stat.txt" u 2:28 w l t "<y^2> MC", "" u 2:(((1.0-(1.0+$40/(-r*r*$9/$8*(1.0-$9/$8)/(1.0+r*$9/--with-gcc-toolchain$8)**2.0)))*(1.0+r)*$9/$8*(1.0-$9/$8)/(1.0+r*$9/$8)**2.0) + (1.0+$40/(-r*r*$9/$8*(1.0-$9/$8)/(1.0+r*$9/$8)**2.0))*(-$40*(1.0-$40)*(1.0+r*$9/$8)**2.0)/r/r) w l t "<y^2> linear model light=heavy" # thetab = 1.0+$44/(-r*r*$11/$10*(1.0-$11/$10)/(1.0+r*$11/$10)**2.0) r = 9.0; plot "stat.txt" u 2:30 w l t "<y^2> MC", "" u 2:(((1.0-(1.0+$44/(-r*r*$11/$10*(1.0-$11/$10)/(1.0+r*$11/$10)**2.0)))*(1.0+r)*$11/$10*(1.0-$11/$10)/(1.0+r*$11/$10)**2.0) + (1.0+$44/(-r*r*$11/$10*(1.0-$11/$10)/(1.0+r*$11/$10)**2.0))*(-$44*(1.0-$44)*(1.0+r*$11/$10)**2.0)/r/r) w l t "<y^2> linear model light>heavy" # thetab = 1.0+$48/(-r*r*$13/$12*(1.0-$13/$12)/(1.0+r*$13/$12)**2.0) r = 9.0; plot "stat.txt" u 2:32 w l t "<y^2> MC", "" u 2:(((1.0-(1.0+$48/(-r*r*$13/$12*(1.0-$13/$12)/(1.0+r*$13/$12)**2.0)))*(1.0+r)*$13/$12*(1.0-$13/$12)/(1.0+r*$13/$12)**2.0) + (1.0+$48/(-r*r*$13/$12*(1.0-$13/$12)/(1.0+r*$13/$12)**2.0))*(-$48*(1.0-$48)*(1.0+r*$13/$12)**2.0)/r/r) w l t "<y^2> linear model light<heavy" # thetab = 1.0+$52/(-r*r*$15/$14*(1.0-$15/$14)/(1.0+r*$15/$14)**2.0) r = 9.0; plot "stat.txt" u 2:34 w l t "<y^2> MC", "" u 2:(((1.0-(1.0+$52/(-r*r*$15/$14*(1.0-$15/$14)/(1.0+r*$15/$14)**2.0)))*(1.0+r)*$15/$14*(1.0-$15/$14)/(1.0+r*$15/$14)**2.0) + (1.0+$52/(-r*r*$15/$14*(1.0-$15/$14)/(1.0+r*$15/$14)**2.0))*(-$52*(1.0-$52)*(1.0+r*$15/$14)**2.0)/r/r) w l t "<y^2> linear model <r^3>(t=0) = 0" # Should we model b or <y^2> ? # b vs <y^2> r = 9.0; set xlabel "<y^2>"; set ylabel "b"; plot "stat.txt" u ($28):(-$40) w l t "light=heavy", "" u ($30):(-$44) w l t "light>heavy", "" u ($32):(-$48) w l t "light<heavy", "" u ($34):(-$52) w l t "<r^3>(t=0) = 0" # b/bnm vs <y^2> r = 9.0; set xlabel "<y^2>"; set ylabel "b/bnm"; plot "stat.txt" u ($28):(-$40/(r*r*$9/$8*(1.0-$9/$8)/(1.0+r*$9/$8)**2.0)) w l t "light=heavy", "" u ($30):(-$44/(r*r*$11/$10*(1.0-$11/$10)/(1.0+r*$11/$10)**2.0)) w l t "light>heavy", "" u ($32):(-$48/(r*r*$13/$12*(1.0-$13/$12)/(1.0+r*$13/$12)**2.0)) w l t "light<heavy", "" u ($34):(-$52/(r*r*$15/$14*(1.0-$15/$14)/(1.0+r*$15/$14)**2.0)) w l t "<r^3>(t=0) = 0" # b/bnm vs <y^2>/<y^2>nm r = 9.0; set xlabel "<y^2>/<y^2>nm"; set ylabel "b/bnm"; plot "stat.txt" u ($28/((1.0+r)*$9/$8*(1.0-$9/$8)/(1.0+r*$9/$8)**2.0)):(-$40/(r*r*$9/$8*(1.0-$9/$8)/(1.0+r*$9/$8)**2.0)) w l t "light=heavy", "" u ($30/((1.0+r)*$11/$10*(1.0-$11/$10)/(1.0+r*$11/$10)**2.0)):(-$44/(r*r*$11/$10*(1.0-$11/$10)/(1.0+r*$11/$10)**2.0)) w l t "light>heavy", "" u ($32/((1.0+r)*$13/$12*(1.0-$13/$12)/(1.0+r*$13/$12)**2.0)):(-$48/(r*r*$13/$12*(1.0-$13/$12)/(1.0+r*$13/$12)**2.0)) w l t "light<heavy", "" u ($34/((1.0+r)*$15/$14*(1.0-$15/$14)/(1.0+r*$15/$14)**2.0)):(-$52/(r*r*$15/$14*(1.0-$15/$14)/(1.0+r*$15/$14)**2.0)) w l t "<r^3>(t=0) = 0" # Evolution of b/bnm r = 9.0; plot "stat.txt" u 2:(-$40/(r*r*$9/$8*(1.0-$9/$8)/(1.0+r*$9/$8)**2.0)) w l t "b/bnm light=heavy", "" u 2:(-$44/(r*r*$11/$10*(1.0-$11/$10)/(1.0+r*$11/$10)**2.0)) w l t "b/bnm light>heavy", "" u 2:(-$48/(r*r*$13/$12*(1.0-$13/$12)/(1.0+r*$13/$12)**2.0)) w l t "b/bnm light<heavy", "" u 2:(-$52/(r*r*$15/$14*(1.0-$15/$14)/(1.0+r*$15/$14)**2.0)) w l t "b/bnm <r^3>(t=0) = 0" # Evolution of <y^2>/<y^2>nm r = 9.0; plot "stat.txt" u 2:($28/((1.0+r)*$9/$8*(1.0-$9/$8)/(1.0+r*$9/$8)/(1.0+r*$9/$8))) w l t "<y^2>/<y^2>nm light=heavy", "" u 2:($30/((1.0+r)*$11/$10*(1.0-$11/$10)/(1.0+r*$11/$10)/(1.0+r*$11/$10))) w l t "<y^2>/<y^2>nm light>heavy", "" u 2:($32/((1.0+r)*$13/$12*(1.0-$13/$12)/(1.0+r*$13/$12)/(1.0+r*$13/$12))) w l t "<y^2>/<y^2>nm light<heavy", "" u 2:($34/((1.0+r)*$15/$14*(1.0-$15/$14)/(1.0+r*$15/$14)/(1.0+r*$15/$14))) w l t "<y^2>/<y^2>nm <r^3>(t=0) = 0" r = 2.0; plot "stat.txt" u 2:($28/((1.0+r)*$9/$8*(1.0-$9/$8)/(1.0+r*$9/$8)/(1.0+r*$9/$8))) w l t "<y^2>/<y^2>nm light=heavy", "" u 2:($30/((1.0+r)*$11/$10*(1.0-$11/$10)/(1.0+r*$11/$10)/(1.0+r*$11/$10))) w l t "<y^2>/<y^2>nm light>heavy", "" u 2:($32/((1.0+r)*$13/$12*(1.0-$13/$12)/(1.0+r*$13/$12)/(1.0+r*$13/$12))) w l t "<y^2>/<y^2>nm light<heavy", "" u 2:($34/((1.0+r)*$15/$14*(1.0-$15/$14)/(1.0+r*$15/$14)/(1.0+r*$15/$14))) w l t "<y^2>/<y^2>nm <r^3>(t=0) = 0" # Evolution of <r^2>/<r^2>nm, <r^2>nm=(<R>-rho1)(rho2-<R>), high At: rho1=0.1, low At: rho1=0.99 r1=0.1; plot "stat.txt" u 2:($38/($8-r1)/(1.0-$8)) w l t "<r^2>/<r^2>nm light=heavy", "" u 2:($42/($10-r1)/(1.0-$10)) w l t "<r^2>/<r^2>nm light>heavy", "" u 2:($46/($12-r1)/(1.0-$12)) w l t "<r^2>/<r^2>nm light<heavy", "stat.txt" u 2:($50/($14-r1)/(1.0-$14)) w l t "<r^2>/<r^2>nm <r^3>(t=0) = 0" # <y^2>MC vs <y^2>linear_closure_model r = 9.0; set xlabel "<y^2>MC"; set ylabel "<y^2>linear_model"; plot "stat.txt" u 28:(((1.0-(1.0+$40/(-r*r*$9/$8*(1.0-$9/$8)/(1.0+r*$9/$8)**2.0)))*(1.0+r)*$9/$8*(1.0-$9/$8)/(1.0+r*$9/$8)**2.0) + (1.0+$40/(-r*r*$9/$8*(1.0-$9/$8)/(1.0+r*$9/$8)**2.0))*(-$40*(1.0-$40)*(1.0+r*$9/$8)**2.0)/r/r) w l t "light=heavy", "" u 30:(((1.0-(1.0+$44/(-r*r*$11/$10*(1.0-$11/$10)/(1.0+r*$11/$10)**2.0)))*(1.0+r)*$11/$10*(1.0-$11/$10)/(1.0+r*$11/$10)**2.0) + (1.0+$44/(-r*r*$11/$10*(1.0-$11/$10)/(1.0+r*$11/$10)**2.0))*(-$44*(1.0-$44)*(1.0+r*$11/$10)**2.0)/r/r) w l lt 2 t "light>heavy", "" u 32:(((1.0-(1.0+$48/(-r*r*$13/$12*(1.0-$13/$12)/(1.0+r*$13/$12)**2.0)))*(1.0+r)*$13/$12*(1.0-$13/$12)/(1.0+r*$13/$12)**2.0) + (1.0+$48/(-r*r*$13/$12*(1.0-$13/$12)/(1.0+r*$13/$12)**2.0))*(-$48*(1.0-$48)*(1.0+r*$13/$12)**2.0)/r/r) w l lt 3 t "light<heavy", "" u 34:(((1.0-(1.0+$52/(-r*r*$15/$14*(1.0-$15/$14)/(1.0+r*$15/$14)**2.0)))*(1.0+r)*$15/$14*(1.0-$15/$14)/(1.0+r*$15/$14)**2.0) + (1.0+$52/(-r*r*$15/$14*(1.0-$15/$14)/(1.0+r*$15/$14)**2.0))*(-$52*(1.0-$52)*(1.0+r*$15/$14)**2.0)/r/r) w l lt 4 t "<r^3>(t=0) = 0" # <Y>(1-<Y>) / <y^2>nm r = 9.0; plot "stat.txt" u 2:($3*(1.0-$3)/((1.0+r)*$9/$8*(1.0-$9/$8)/(1.0+r*$9/$8)/(1.0+r*$9/$8))) w l, "" u 2:($4*(1.0-$4)/((1.0+r)*$11/$10*(1.0-$11/$10)/(1.0+r*$11/$10)/(1.0+r*$11/$10))) w l, "" u 2:($6*(1.0-$6)/((1.0+r)*$13/$12*(1.0-$13/$12)/(1.0+r*$13/$12)/(1.0+r*$13/$12))) w l, "" u 2:($8*(1.0-$8)/((1.0+r)*$15/$14*(1.0-$15/$14)/(1.0+r*$15/$14)/(1.0+r*$15/$14))) w l ##### BASE MODEL v2 ##### # test closure model for <rv^2> = -r^3/rho_2 * 1/(1+r) * <y^2>/(1+r<RY>/<R>)^2 * [ 1 - (2-r)<RY>/<R> - 3r(<RY>/<R>)^2 - r^2(<RY>/<R>)^3 - b*(1+r<RY>/<R>)^3/r ], r = 9.0 -> high A, r = 0.0101 -> low A ##### BASE MODEL v2 ##### r = 9.0; plot "stat.txt" u 2:41 w l t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$28/(1.0+r*$9/$8)**2.0*(1.0-(2.0-r)*$9/$8-3.0*r*$9/$8*$9/$8-r*r*($9/$8)**3.0+$40*(1.0+r*$9/$8)**3.0/r)) w l t "<rv^2>nm light=heavy" r = 9.0; plot "stat.txt" u 2:45 w l t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$30/(1.0+r*$11/$10)**2.0*(1.0-(2.0-r)*$11/$10-3.0*r*$11/$10*$11/$10-r*r*($11/$10)**3.0+$44*(1.0+r*$11/$10)**3.0/r)) w l t "<rv^2>nm light=heavy" r = 9.0; plot "stat.txt" u 2:49 w l t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$32/(1.0+r*$13/$12)**2.0*(1.0-(2.0-r)*$13/$12-3.0*r*$13/$12*$13/$12-r*r*($13/$12)**3.0+$48*(1.0+r*$13/$12)**3.0/r)) w l t "<rv^2>nm light=heavy" r = 9.0; plot "stat.txt" u 2:53 w l t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$34/(1.0+r*$15/$14)**2.0*(1.0-(2.0-r)*$15/$14-3.0*r*$15/$14*$15/$14-r*r*($15/$14)**3.0+$52*(1.0+r*$15/$14)**3.0/r)) w l t "<rv^2>nm light=heavy" # all in on plot r = 9.0; plot "stat.txt" u 2:41 w l t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$28/(1.0+r*$9/$8)**2.0*(1.0-(2.0-r)*$9/$8-3.0*r*$9/$8*$9/$8-r*r*($9/$8)**3.0+$40*(1.0+r*$9/$8)**3.0/r)) w p pt 7 lt 1 t "<rv^2>nm light=heavy", "" u 2:45 w l lt 2 t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$30/(1.0+r*$11/$10)**2.0*(1.0-(2.0-r)*$11/$10-3.0*r*$11/$10*$11/$10-r*r*($11/$10)**3.0+$44*(1.0+r*$11/$10)**3.0/r)) w p pt 7 lt 2 t "<rv^2>nm light > heavy", "" u 2:49 w l lt 3 t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$32/(1.0+r*$13/$12)**2.0*(1.0-(2.0-r)*$13/$12-3.0*r*$13/$12*$13/$12-r*r*($13/$12)**3.0+$48*(1.0+r*$13/$12)**3.0/r)) w p pt 7 lt 3 t "<rv^2>nm light < heavy", "" u 2:53 w l lt 4 t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$34/(1.0+r*$15/$14)**2.0*(1.0-(2.0-r)*$15/$14-3.0*r*$15/$14*$15/$14-r*r*($15/$14)**3.0+$52*(1.0+r*$15/$14)**3.0/r)) w p pt 7 lt 4 t "<rv^2>nm <r^3>(t=0) = 0" ##### BASE MODEL v2 ##### # test closure model for <y^2> = b(1+b) / ( <R>^2 (r/rho2)^2 { 1 - r/rho2 <R>/(1+r) [ 1 - (2-r)<RY>/<R> - 3r(<RY>/<R>)^2 - r^2(<RY>/<R>)^3 - b*(1+r<RY>/<R>)^3/r ] } ) ##### BASE MODEL v2 ##### r = 9.0; plot "stat.txt" u 2:28 w l t "<y^2> MC", "" u 2:(-$40*(1.0-$40)/($8**2.0*r**2.0*(1.0-r*$8/(1.0+r)*(1.0-(2.0-r)*$9/$8-3.0*r*$9/$8*$9/$8-r*r*($9/$8)**3.0+$40*(1.0+r*$9/$8)**3.0/r)))) w l t "<y^2> model light=heavy" r = 9.0; plot "stat.txt" u 2:30 w l t "<y^2> MC", "" u 2:(-$44*(1.0-$44)/($10**2.0*r**2.0*(1.0-r*$10/(1.0+r)*(1.0-(2.0-r)*$11/$10-3.0*r*$11/$10*$11/$10-r*r*($11/$10)**3.0+$44*(1.0+r*$11/$10)**3.0/r)))) w l t "<y^2> model light > heavy" r = 9.0; plot "stat.txt" u 2:32 w l t "<y^2> MC", "" u 2:(-$48*(1.0-$48)/($12**2.0*r**2.0*(1.0-r*$12/(1.0+r)*(1.0-(2.0-r)*$13/$12-3.0*r*$13/$12*$13/$12-r*r*($13/$12)**3.0+$48*(1.0+r*$13/$12)**3.0/r)))) w l t "<y^2> model light < heavy" r = 9.0; plot "stat.txt" u 2:34 w l t "<y^2> MC", "" u 2:(-$52*(1.0-$52)/($14**2.0*r**2.0*(1.0-r*$14/(1.0+r)*(1.0-(2.0-r)*$15/$14-3.0*r*$15/$14*$15/$14-r*r*($15/$14)**3.0+$52*(1.0+r*$15/$14)**3.0/r)))) w l t "<y^2> model <r^3>(t=0) = 0" # all in on plot r = 9.0; plot "stat.txt" u 2:28 w l t "<y^2> MC", "" u 2:(-$40*(1.0-$40)/($8**2.0*r**2.0*(1.0-r*$8/(1.0+r)*(1.0-(2.0-r)*$9/$8-3.0*r*$9/$8*$9/$8-r*r*($9/$8)**3.0+$40*(1.0+r*$9/$8)**3.0/r)))) w p pt 7 lt 1 t "<y^2> model light=heavy", "" u 2:30 w l lt 2 t "<y^2> MC", "" u 2:(-$44*(1.0-$44)/($10**2.0*r**2.0*(1.0-r*$10/(1.0+r)*(1.0-(2.0-r)*$11/$10-3.0*r*$11/$10*$11/$10-r*r*($11/$10)**3.0+$44*(1.0+r*$11/$10)**3.0/r)))) w p pt 7 lt 2 t "<y^2> model light > heavy", "" u 2:32 w l lt 3 t "<y^2> MC", "" u 2:(-$48*(1.0-$48)/($12**2.0*r**2.0*(1.0-r*$12/(1.0+r)*(1.0-(2.0-r)*$13/$12-3.0*r*$13/$12*$13/$12-r*r*($13/$12)**3.0+$48*(1.0+r*$13/$12)**3.0/r)))) w p pt 7 lt 3 t "<y^2> model light < heavy", "" u 2:34 w l lt 4 t "<y^2> MC", "" u 2:(-$52*(1.0-$52)/($14**2.0*r**2.0*(1.0-r*$14/(1.0+r)*(1.0-(2.0-r)*$15/$14-3.0*r*$15/$14*$15/$14-r*r*($15/$14)**3.0+$52*(1.0+r*$15/$14)**3.0/r)))) w p pt 7 lt 4 t "<y^2> model <r^3>(t=0) = 0" ##### BASE MODEL v3 ##### # test closure model for <rv^2> = -r^3/rho_2 * 1/(1+r) * <y^2>/(1+r<RY>/<R>)^2 * [ 1 - (2-r)<RY>/<R> - 3r(<RY>/<R>)^2 - r^2(<RY>/<R>)^3 - r/(1+r)*<y^2>*(1+r<RY>/<R>)^3 ], r = 9.0 -> high A, r = 0.0101 -> low A ##### BASE MODEL v3 ##### r = 9.0; plot "stat.txt" u 2:41 w l t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$28/(1.0+r*$9/$8)**2.0*(1.0-(2.0-r)*$9/$8-3.0*r*$9/$8*$9/$8-r*r*($9/$8)**3.0-r/(1.0+r)*$28*(1.0+r*$9/$8)**3.0)) w l t "<rv^2>nm light=heavy" r = 9.0; plot "stat.txt" u 2:45 w l t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$30/(1.0+r*$11/$10)**2.0*(1.0-(2.0-r)*$11/$10-3.0*r*$11/$10*$11/$10-r*r*($11/$10)**3.0-r/(1.0+r)*$30*(1.0+r*$11/$10)**3.0)) w l t "<rv^2>nm light=heavy" r = 9.0; plot "stat.txt" u 2:49 w l t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$32/(1.0+r*$13/$12)**2.0*(1.0-(2.0-r)*$13/$12-3.0*r*$13/$12*$13/$12-r*r*($13/$12)**3.0-r/(1.0+r)*$32*(1.0+r*$13/$12)**3.0)) w l t "<rv^2>nm light=heavy" r = 9.0; plot "stat.txt" u 2:53 w l t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$34/(1.0+r*$15/$14)**2.0*(1.0-(2.0-r)*$15/$14-3.0*r*$15/$14*$15/$14-r*r*($15/$14)**3.0-r/(1.0+r)*$34*(1.0+r*$15/$14)**3.0)) w l t "<rv^2>nm light=heavy" # all in on plot r = 9.0; plot "stat.txt" u 2:41 w l t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$28/(1.0+r*$9/$8)**2.0*(1.0-(2.0-r)*$9/$8-3.0*r*$9/$8*$9/$8-r*r*($9/$8)**3.0-r/(1.0+r)*$28*(1.0+r*$9/$8)**3.0)) w p pt 7 lt 1 t "<rv^2>nm light=heavy", "" u 2:45 w l lt 2 t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$30/(1.0+r*$11/$10)**2.0*(1.0-(2.0-r)*$11/$10-3.0*r*$11/$10*$11/$10-r*r*($11/$10)**3.0-r/(1.0+r)*$30*(1.0+r*$11/$10)**3.0)) w p pt 7 lt 2 t "<rv^2>nm light > heavy", "" u 2:49 w l lt 3 t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$32/(1.0+r*$13/$12)**2.0*(1.0-(2.0-r)*$13/$12-3.0*r*$13/$12*$13/$12-r*r*($13/$12)**3.0-r/(1.0+r)*$32*(1.0+r*$13/$12)**3.0)) w p pt 7 lt 3 t "<rv^2>nm light < heavy", "" u 2:53 w l lt 4 t "<rv^2>", "" u 2:(-r**3.0/(1.0+r)*$34/(1.0+r*$15/$14)**2.0*(1.0-(2.0-r)*$15/$14-3.0*r*$15/$14*$15/$14-r*r*($15/$14)**3.0-r/(1.0+r)*$34*(1.0+r*$15/$14)**3.0)) w p pt 7 lt 4 t "<rv^2>nm <r^3>(t=0) = 0" ##### BASE MODEL v4 ##### # test closure model for <rv^2> = -(1+r)/rho2 * b/(1+r<RY>/<R>) * [ 1 - (1+r<RY>/<R>)^2*(1+b)/(1+r) ], r = 9.0 -> high A, r = 0.0101 -> low A ##### BASE MODEL v4 ##### r = 9.0; plot "stat.txt" u 2:41 w l t "<rv^2>", "" u 2:(-(1.0+r)*(-$40)/(1.0+r*$9/$8)*(1.0-(1.0+r*$9/$8)**2.0/(1.0+r)*(1.0-$40))) w l t "<rv^2>nm light=heavy" r = 9.0; plot "stat.txt" u 2:45 w l t "<rv^2>", "" u 2:(-(1.0+r)*(-$44)/(1.0+r*$11/$10)*(1.0-(1.0+r*$11/$10)**2.0/(1.0+r)*(1.0-$44))) w l t "<rv^2>nm light=heavy" r = 9.0; plot "stat.txt" u 2:49 w l t "<rv^2>", "" u 2:(-(1.0+r)*(-$48)/(1.0+r*$13/$12)*(1.0-(1.0+r*$13/$12)**2.0/(1.0+r)*(1.0-$48))) w l t "<rv^2>nm light=heavy" r = 9.0; plot "stat.txt" u 2:53 w l t "<rv^2>", "" u 2:(-(1.0+r)*(-$52)/(1.0+r*$15/$14)*(1.0-(1.0+r*$15/$14)**2.0/(1.0+r)*(1.0-$52))) w l t "<rv^2>nm light=heavy" # all in on plot r = 9.0; plot "stat.txt" u 2:41 w l t "<rv^2>", "" u 2:(-(1.0+r)*(-$40)/(1.0+r*$9/$8)*(1.0-(1.0+r*$9/$8)**2.0/(1.0+r)*(1.0-$40))) w p pt 7 lt 1 t "<rv^2>nm light=heavy", "" u 2:45 w l lt 2 t "<rv^2>", "" u 2:(-(1.0+r)*(-$44)/(1.0+r*$11/$10)*(1.0-(1.0+r*$11/$10)**2.0/(1.0+r)*(1.0-$44))) w p pt 7 lt 2 t "<rv^2>nm light > heavy", "" u 2:49 w l lt 3 t "<rv^2>", "" u 2:(-(1.0+r)*(-$48)/(1.0+r*$13/$12)*(1.0-(1.0+r*$13/$12)**2.0/(1.0+r)*(1.0-$48))) w p pt 7 lt 3 t "<rv^2>nm light < heavy", "" u 2:53 w l lt 4 t "<rv^2>", "" u 2:(-(1.0+r)*(-$52)/(1.0+r*$15/$14)*(1.0-(1.0+r*$15/$14)**2.0/(1.0+r)*(1.0-$52))) w p pt 7 lt 4 t "<rv^2>nm <r^3>(t=0) = 0" ##### BASE MODEL v5 ##### # test closure model for <rv^2> = -r^2/rho2 * <y^2>/(1+r<RY>/<R>) * [ 1 - (1+r<RY>/<R>)^2*(1+b)/(1+r) ], r = 9.0 -> high A, r = 0.0101 -> low A ##### BASE MODEL v5 ##### r = 9.0; plot "stat.txt" u 2:41 w l t "<rv^2>", "" u 2:(-r*r*$28/(1.0+r*$9/$8)*(1.0-(1.0+r*$9/$8)**2.0/(1.0+r)*(1.0-$40))) w l t "<rv^2>nm light=heavy" r = 9.0; plot "stat.txt" u 2:45 w l t "<rv^2>", "" u 2:(-r*r*$30/(1.0+r*$11/$10)*(1.0-(1.0+r*$11/$10)**2.0/(1.0+r)*(1.0-$44))) w l t "<rv^2>nm light=heavy" r = 9.0; plot "stat.txt" u 2:49 w l t "<rv^2>", "" u 2:(-r*r*$32/(1.0+r*$13/$12)*(1.0-(1.0+r*$13/$12)**2.0/(1.0+r)*(1.0-$48))) w l t "<rv^2>nm light=heavy" r = 9.0; plot "stat.txt" u 2:53 w l t "<rv^2>", "" u 2:(-r*r*$34/(1.0+r*$15/$14)*(1.0-(1.0+r*$15/$14)**2.0/(1.0+r)*(1.0-$52))) w l t "<rv^2>nm light=heavy" # all in on plot r = 9.0; plot "stat.txt" u 2:41 w l t "<rv^2>", "" u 2:(-r*r*$28/(1.0+r*$9/$8)*(1.0-(1.0+r*$9/$8)**2.0/(1.0+r)*(1.0-$40))) w p pt 7 lt 1 t "<rv^2>nm light=heavy", "" u 2:45 w l lt 2 t "<rv^2>", "" u 2:(-r*r*$30/(1.0+r*$11/$10)*(1.0-(1.0+r*$11/$10)**2.0/(1.0+r)*(1.0-$44))) w p pt 7 lt 2 t "<rv^2>nm light > heavy", "" u 2:49 w l lt 3 t "<rv^2>", "" u 2:(-r*r*$32/(1.0+r*$13/$12)*(1.0-(1.0+r*$13/$12)**2.0/(1.0+r)*(1.0-$48))) w p pt 7 lt 3 t "<rv^2>nm light < heavy", "" u 2:53 w l lt 4 t "<rv^2>", "" u 2:(-r*r*$34/(1.0+r*$15/$14)*(1.0-(1.0+r*$15/$14)**2.0/(1.0+r)*(1.0-$52))) w p pt 7 lt 4 t "<rv^2>nm <r^3>(t=0) = 0"