Viscosity Coefficients
Contents
Viscosity Coefficients#
Prelude#
In this notebook we will calculate the viscosity of the Yukawa OCP.
The YAML input file can be found at input_file and this notebook at notebook.
[1]:
# Import the usual libraries
%pylab
%matplotlib inline
import os
plt.style.use('MSUstyle')
# Import sarkas
from sarkas.processes import Simulation, PostProcess, PreProcess
# Create the file path to the YAML input file
input_file_name = os.path.join('input_files', 'yocp_transport.yaml')
Using matplotlib backend: <object object at 0x7f56d85b3ff0>
%pylab is deprecated, use %matplotlib inline and import the required libraries.
Populating the interactive namespace from numpy and matplotlib
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
File ~/checkouts/readthedocs.org/user_builds/sarkas/conda/dev/lib/python3.9/site-packages/matplotlib/style/core.py:137, in use(style)
136 try:
--> 137 style = _rc_params_in_file(style)
138 except OSError as err:
File ~/checkouts/readthedocs.org/user_builds/sarkas/conda/dev/lib/python3.9/site-packages/matplotlib/__init__.py:866, in _rc_params_in_file(fname, transform, fail_on_error)
865 rc_temp = {}
--> 866 with _open_file_or_url(fname) as fd:
867 try:
File ~/checkouts/readthedocs.org/user_builds/sarkas/conda/dev/lib/python3.9/contextlib.py:119, in _GeneratorContextManager.__enter__(self)
118 try:
--> 119 return next(self.gen)
120 except StopIteration:
File ~/checkouts/readthedocs.org/user_builds/sarkas/conda/dev/lib/python3.9/site-packages/matplotlib/__init__.py:843, in _open_file_or_url(fname)
842 fname = os.path.expanduser(fname)
--> 843 with open(fname, encoding='utf-8') as f:
844 yield f
FileNotFoundError: [Errno 2] No such file or directory: 'MSUstyle'
The above exception was the direct cause of the following exception:
OSError Traceback (most recent call last)
Cell In[1], line 6
3 get_ipython().run_line_magic('matplotlib', 'inline')
5 import os
----> 6 plt.style.use('MSUstyle')
8 # Import sarkas
9 from sarkas.processes import Simulation, PostProcess, PreProcess
File ~/checkouts/readthedocs.org/user_builds/sarkas/conda/dev/lib/python3.9/site-packages/matplotlib/style/core.py:139, in use(style)
137 style = _rc_params_in_file(style)
138 except OSError as err:
--> 139 raise OSError(
140 f"{style!r} is not a valid package style, path of style "
141 f"file, URL of style file, or library style name (library "
142 f"styles are listed in `style.available`)") from err
143 filtered = {}
144 for k in style: # don't trigger RcParams.__getitem__('backend')
OSError: 'MSUstyle' is not a valid package style, path of style file, URL of style file, or library style name (library styles are listed in `style.available`)
[2]:
# pre = PreProcess(input_file_name)
# pre.setup(read_yaml=True)
# pre.run(loops=150)
[3]:
# sim = Simulation(input_file_name)
# sim.setup(read_yaml=True)
# sim.run()
[4]:
postproc = PostProcess(input_file_name)
postproc.setup(read_yaml=True)
# postproc.parameters.verbose = True
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[4], line 1
----> 1 postproc = PostProcess(input_file_name)
2 postproc.setup(read_yaml=True)
3 # postproc.parameters.verbose = True
NameError: name 'PostProcess' is not defined
[5]:
from sarkas.tools.observables import Thermodynamics, PressureTensor, HeatFlux, RadialDistributionFunction, VelocityAutoCorrelationFunction
[6]:
therm = Thermodynamics()
therm.setup(postproc.parameters, no_slices=3 )
therm.compute()
therm.grab_sim_data()
therm.temp_energy_plot(postproc)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[6], line 2
1 therm = Thermodynamics()
----> 2 therm.setup(postproc.parameters, no_slices=3 )
3 therm.compute()
4 therm.grab_sim_data()
NameError: name 'postproc' is not defined
Pair Distribution Function#
The first observable to calculate is always the RDF.
[7]:
rdf = RadialDistributionFunction()
rdf.setup(postproc.parameters, no_slices=3 )
rdf.parse()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[7], line 2
1 rdf = RadialDistributionFunction()
----> 2 rdf.setup(postproc.parameters, no_slices=3 )
3 rdf.parse()
NameError: name 'postproc' is not defined
[8]:
rdf.plot(scaling = rdf.a_ws,
y = ('C-C RDF', 'Mean'),
xlabel = r'$r /a$')
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In[8], line 1
----> 1 rdf.plot(scaling = rdf.a_ws,
2 y = ('C-C RDF', 'Mean'),
3 xlabel = r'$r /a$')
AttributeError: 'RadialDistributionFunction' object has no attribute 'a_ws'
Pressure Tensor#
The viscosity is obtained from the autocorrelation function of the Pressure Tensor \(\overleftrightarrow{\mathcal P}\) whose elements are
\begin{equation} \mathcal P_{\alpha\gamma}(t) = \frac{1}{V} \sum_{i}^{N} \left [ m_i v^{\alpha}_{i} v^{\gamma}_{i} - \sum_{j > i} \frac{r_{ij}^{\alpha} r_{ij}^{\gamma} }{r_{ij}} \frac{d}{dr}\phi(r) \right ], \end{equation}
where \(r_{ij}^{\alpha}\) is the \(\alpha\) component of the distance between particles \(i\) and \(j\). The first term is the kinetic term and the second term is the virial term, but it is often referred to as the potential contribution. The virial is calculated during the simulation phase and saved together with particles corrdinates.
In order to check that our code are correct, let’s verify some laws.
The pressure of the system is calculated from \(\mathcal P(t)= \frac1{3} {\rm Tr} \overleftrightarrow{\mathcal P}(t)\) and also from
\begin{equation} P = \frac{n}{\beta} - \frac{2\pi}{3} n^2 \int_0^{\infty} dr \, r^3 \frac{d\phi(r)}{dr} g(r) \end{equation}
where \(g(r)\) is the pair distribution function that we have already calculated.
Let’s calculate the Pressure tensor and the pressure \(\mathcal P\).
[9]:
pt = PressureTensor()
pt.setup(postproc.parameters, no_slices = 3)
# pt.compute()
pt.parse(acf_data=True)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[9], line 2
1 pt = PressureTensor()
----> 2 pt.setup(postproc.parameters, no_slices = 3)
3 # pt.compute()
4 pt.parse(acf_data=True)
NameError: name 'postproc' is not defined
As usual the data is saved in several dataframes. In this case we have 4 dataframes
A dataframe for the values of each of the elements of the pressure tensor for each of the slices,
pt.dataframe_slicesA dataframe for the mean and std values of each of the elements of the pressure tensor,
pt.dataframeA dataframe for the ACF of each pair \(\langle \mathcal P_{\alpha\beta}(t)\mathcal P_{\mu\nu}(0) \rangle\) for each slice,
pt.dataframe_acf_slicesA dataframe for the mean and std of the ACF of each pair \(\langle \mathcal P_{\alpha\beta}(t)\mathcal P_{\mu\nu}(0) \rangle\),
pt.dataframe_acf
Let’s look at pt.dataframe and at its columns
[10]:
pt.dataframe
Note that the Pressure \(\mathcal P(t)\) is readily calculated and provided as a column of the dataframe.
Note also that there is a multitude of columns. This is because in dense plasmas it is useful to know the contribution of both the kinetic term and potential term separately, as such the columns of each dataframe contain the kinetic, the potential, and the total value of each \(\mathcal P_{\alpha\beta}\) and their ACFs.
Let’s plot the Pressure as a function of time
[11]:
# Let's plot it
t_wp = 2.0*np.pi/therm.total_plasma_frequency
p_id = pt.total_num_density / therm.beta_slice.mean()
ax = pt.plot(
scaling = (t_wp, p_id),
y = ("Total","Pressure", "Mean"),
xlabel = "Plasma cycles",
ylabel = r"$ \beta P(t)/n$"
)
ax.plot(pt.dataframe[("Species", "Quantity", 'Time')]/t_wp, pt.dataframe[("Total",'Pressure','Mean')].expanding().mean()/p_id )
ax.legend(['Pressure', 'Moving avg'])
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In[11], line 2
1 # Let's plot it
----> 2 t_wp = 2.0*np.pi/therm.total_plasma_frequency
4 p_id = pt.total_num_density / therm.beta_slice.mean()
5 ax = pt.plot(
6 scaling = (t_wp, p_id),
7 y = ("Total","Pressure", "Mean"),
8 xlabel = "Plasma cycles",
9 ylabel = r"$ \beta P(t)/n$"
10 )
AttributeError: 'Thermodynamics' object has no attribute 'total_plasma_frequency'
[12]:
from sarkas.tools.observables import make_gaussian_plot
[13]:
pt.dataframe.columns
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In[13], line 1
----> 1 pt.dataframe.columns
AttributeError: 'NoneType' object has no attribute 'columns'
[14]:
make_gaussian_plot(pt.dataframe[("Species", "Quantity", 'Time')]/t_wp, pt.dataframe[("C",'Pressure Tensor XY','Mean')], "Pressure Tensor XY", "mks")
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[14], line 1
----> 1 make_gaussian_plot(pt.dataframe[("Species", "Quantity", 'Time')]/t_wp, pt.dataframe[("C",'Pressure Tensor XY','Mean')], "Pressure Tensor XY", "mks")
TypeError: 'NoneType' object is not subscriptable
Pressure from RDF#
Let’s now calculate the pressure from the integral of the RDF. This is obtained from the method compute_from_rdf of the Thermodynamics object.
Looking at the documentation of this method we notice that it returns five values: the Hartree and correlational terms between species :math:A and :math:B and the ideal pressure \(n k_B T\).
The total pressure is given from the sum of the three terms and should be equal to the
[15]:
nkT, _, _, p_h, p_c = therm.compute_from_rdf(rdf, postproc.potential)
P_rdf = nkT + p_h + p_c
P_trace = pt.dataframe[("Total","Pressure", "Mean")].mean()
print("The relative difference between the two methods is = {:.2f} %".format((P_rdf[0] - P_trace)*100/P_rdf[0] ) )
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[15], line 1
----> 1 nkT, _, _, p_h, p_c = therm.compute_from_rdf(rdf, postproc.potential)
3 P_rdf = nkT + p_h + p_c
4 P_trace = pt.dataframe[("Total","Pressure", "Mean")].mean()
NameError: name 'postproc' is not defined
It seems that we have done a good job!
Sum rule#
Let’s now check that we have calculated the ACF correctly. The equal time ACFs of the elements of \(\overleftrightarrow{\mathcal P}(t)\) obey the following sum rules
where
Notice that all three equal time ACF satisfy
Let’s look at the dataframe of the ACF first
[16]:
pt.dataframe_acf
Notice that in this case we have many more columns since now we have the ACF of the kinetic-kinetic, kinetic-potential, potential-kinetic, potential-potential, and the total ACF of each pair of elements.
Let’s verify the sum rules.
[17]:
# Diagonal terms
column_zzzz = [
('Pressure Tensor ACF XXXX', 'Mean'),
('Pressure Tensor ACF YYYY', 'Mean'),
('Pressure Tensor ACF ZZZZ', 'Mean'),
]
J_zzzz_0 = pt.dataframe_acf[column_zzzz].iloc[0].mean()
# Cross-Diagonal Terms
column_zzxx = [
('Pressure Tensor ACF XXYY', 'Mean'),
('Pressure Tensor ACF XXZZ', 'Mean'),
('Pressure Tensor ACF YYZZ', 'Mean'),
]
J_zzxx_0 = pt.dataframe_acf[column_zzxx].iloc[0].mean()
# Cross Off Diagonal terms
column_xyxy = [
('Pressure Tensor ACF XYXY', 'Mean'),
('Pressure Tensor ACF XZXZ', 'Mean'),
('Pressure Tensor ACF YZYZ', 'Mean'),
]
J_xyxy_0 = pt.dataframe_acf[column_xyxy].iloc[0].mean()
# The units of J's are [Density * Energy]^2
print('The isotropy condition : (J_zzzz_0 - J_zzxx_0 )/( 2*J_xyxy_0 ) = {:.4f}'.format( (J_zzzz_0 - J_zzxx_0)/(2.0 * J_xyxy_0) ))
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[17], line 7
1 # Diagonal terms
2 column_zzzz = [
3 ('Pressure Tensor ACF XXXX', 'Mean'),
4 ('Pressure Tensor ACF YYYY', 'Mean'),
5 ('Pressure Tensor ACF ZZZZ', 'Mean'),
6 ]
----> 7 J_zzzz_0 = pt.dataframe_acf[column_zzzz].iloc[0].mean()
9 # Cross-Diagonal Terms
10 column_zzxx = [
11 ('Pressure Tensor ACF XXYY', 'Mean'),
12 ('Pressure Tensor ACF XXZZ', 'Mean'),
13 ('Pressure Tensor ACF YYZZ', 'Mean'),
14 ]
TypeError: 'NoneType' object is not subscriptable
Not exactly 1 but pretty close.
Let’s now verify the sum rules. These are calculated from the pt.sum_rule method
[18]:
h_r, c_r = rdf.compute_sum_rule_integrals(postproc.potential)
sigma_zzzz, sigma_zzxx, sigma_xyxy = pt.sum_rule(therm.beta_slice.mean(), rdf, postproc.potential)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[18], line 1
----> 1 h_r, c_r = rdf.compute_sum_rule_integrals(postproc.potential)
2 sigma_zzzz, sigma_zzxx, sigma_xyxy = pt.sum_rule(therm.beta_slice.mean(), rdf, postproc.potential)
NameError: name 'postproc' is not defined
[19]:
G_inf = J_xyxy_0*therm.beta_slice.mean()*rdf.box_volume
K_inf = 1.0/3.0*(J_zzzz_0 + 2.0* J_zzxx_0)*therm.beta_slice.mean()*rdf.box_volume
K_sr = (sigma_zzzz + 2.0*sigma_zzxx)/3.0
const = 1.0 #postproc.species[0].sigma**3/postproc.species[0].epsilon
print("G_inf = {:2.1f}, sum_rule = {:2.1f}, {:2.2f} %".format( G_inf* const, sigma_xyxy * const, 100*abs(G_inf - sigma_xyxy) /G_inf) )
print("K_inf = {:2.1f}, sum_rule = {:2.1f}, {:2.2f} %".format( K_inf * const, K_sr * const, 100*abs(K_inf - K_sr) /K_inf) )
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[19], line 1
----> 1 G_inf = J_xyxy_0*therm.beta_slice.mean()*rdf.box_volume
2 K_inf = 1.0/3.0*(J_zzzz_0 + 2.0* J_zzxx_0)*therm.beta_slice.mean()*rdf.box_volume
4 K_sr = (sigma_zzzz + 2.0*sigma_zzxx)/3.0
NameError: name 'J_xyxy_0' is not defined
Viscosity#
The shear viscosity is calculated from the Green-Kubo relation
\begin{equation} \eta = \frac{\beta V}{3} \sum_{\alpha} \sum_{\gamma \neq \alpha} \int_0^{\infty} dt \, \left \langle \mathcal P_{\alpha\gamma}(t) \mathcal P_{\alpha\gamma}(0) \right \rangle, \end{equation}
where \(\beta = 1/k_B T\), \(\alpha,\gamma = {x, y, z}\).
The bulk viscosity is given by a similar relation
\begin{equation} \eta_V = \beta V \int_0^{\infty}dt \, \left \langle \delta \mathcal P(t) \delta \mathcal P(0) \right \rangle, \end{equation}
where
\begin{equation} \delta \mathcal P(t) = \mathcal P(t) - \left \langle \mathcal P \right \rangle \end{equation}
is the deviation of the scalar pressure.
[20]:
from sarkas.tools.transport import Viscosity
[21]:
tc = Viscosity()
tc.setup(postproc.parameters, observable=pt)
tc.compute(observable = pt)
# tc.parse(observable = pt, tc_name = "Viscosities")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[21], line 2
1 tc = Viscosity()
----> 2 tc.setup(postproc.parameters, observable=pt)
3 tc.compute(observable = pt)
4 # tc.parse(observable = pt, tc_name = "Viscosities")
NameError: name 'postproc' is not defined
[22]:
acf_str = "Delta Pressure ACF"
acf_avg = pt.dataframe_acf[("Pressure Bulk ACF", "Mean")]
acf_std = pt.dataframe_acf[("Pressure Bulk ACF", "Std")]
pq = "Bulk Viscosity"
tc_avg = tc.dataframe[(pq, "Mean")]
tc_std = tc.dataframe[(pq, "Std")]
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[22], line 2
1 acf_str = "Delta Pressure ACF"
----> 2 acf_avg = pt.dataframe_acf[("Pressure Bulk ACF", "Mean")]
3 acf_std = pt.dataframe_acf[("Pressure Bulk ACF", "Std")]
5 pq = "Bulk Viscosity"
TypeError: 'NoneType' object is not subscriptable
[23]:
tc.dataframe
[24]:
fig, axes = tc.plot_tc(
time = tc.dataframe[("Integration","Interval")].to_numpy(),
acf_data=np.column_stack((acf_avg, acf_std)),
tc_data=np.column_stack((tc_avg, tc_std)),
acf_name=acf_str,
tc_name="Bulk Viscosity",
figname="{}_Plot.png".format("Bulk Viscosity"),
show=False
)
axes[0].set(ylim = (-1, 1.05))
# axes[1].set(ylim = (-0.5, 1000 ) )
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[24], line 2
1 fig, axes = tc.plot_tc(
----> 2 time = tc.dataframe[("Integration","Interval")].to_numpy(),
3 acf_data=np.column_stack((acf_avg, acf_std)),
4 tc_data=np.column_stack((tc_avg, tc_std)),
5 acf_name=acf_str,
6 tc_name="Bulk Viscosity",
7 figname="{}_Plot.png".format("Bulk Viscosity"),
8 show=False
9 )
10 axes[0].set(ylim = (-1, 1.05))
11 # axes[1].set(ylim = (-0.5, 1000 ) )
TypeError: 'NoneType' object is not subscriptable
[25]:
def murillo_yvm(kappa, gamma):
Ak = 0.46 *kappa**4/(1 + 0.44 * kappa**4)
Bk = 1.01*np.exp(-0.92 * kappa)
Ck = -3.7e-5 + 9.0e-4 * kappa - 2.9e-4*kappa**2
gamma_ocp = Ak + Bk * gamma + Ck*gamma**2
lambda_yvm = 4.0 * np.pi/3.0 * (3.0 * gamma_ocp)**(3/2)
I1 = 1.0/ (180 * gamma_ocp * np.pi **(3/2) )
I2 = (0.49 - 2.23 * gamma_ocp**(-1/3) )/ (60 *np.pi**2)
I3 = 0.241 * gamma_ocp**(1/9)/np.pi**(3/2)
eta = lambda_yvm * I1 + (1 + lambda_yvm * I2)**2/(lambda_yvm * I3)
return eta
[26]:
pq = "Shear Viscosity"
tc_avg = tc.dataframe[(pq, "Mean")]
tc_std = tc.dataframe[(pq, "Std")]
rescale = pt.total_plasma_frequency * pt.a_ws**2 * pt.species_masses[0] * pt.total_num_density
fig, ax = plt.subplots(1,1)
ax.plot(tc.dataframe[("Integration","Interval")].to_numpy()*1e12,
tc_avg / rescale,
label = r'$\mu$')
ax.fill_between(
tc.dataframe[("Integration","Interval")].to_numpy()*1e12,
(tc_avg - tc_std) / rescale,
(tc_avg + tc_std) / rescale,
alpha = 0.2)
ax.plot(tc.dataframe[("Integration","Interval")].to_numpy()*1e12,
tc_avg.expanding().mean()/rescale,
label = r'Moving avg')
ax.set(xlabel = r'Time lag $\tau$ [ps]',
ylabel = r"Shear viscosity $\eta$",
xscale= 'log'
)
eta_yvm = murillo_yvm(postproc.potential.kappa, postproc.potential.coupling_constant)
ax.axhline(0.0654, ls = '--', c = 'r', label = "Daligault MD")
ax.axhline(eta_yvm, ls = ':', c = 'r', label = "Murillo YVM")
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[26], line 2
1 pq = "Shear Viscosity"
----> 2 tc_avg = tc.dataframe[(pq, "Mean")]
3 tc_std = tc.dataframe[(pq, "Std")]
6 rescale = pt.total_plasma_frequency * pt.a_ws**2 * pt.species_masses[0] * pt.total_num_density
TypeError: 'NoneType' object is not subscriptable
Thermal Conductivity#
[27]:
from sarkas.tools.observables import HeatFlux
from sarkas.tools.transport import ThermalConductivity
[28]:
ht = HeatFlux()
ht.setup(postproc.parameters, no_slices=3)
ht.compute(calculate_acf=True)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[28], line 2
1 ht = HeatFlux()
----> 2 ht.setup(postproc.parameters, no_slices=3)
3 ht.compute(calculate_acf=True)
NameError: name 'postproc' is not defined
[29]:
thc = ThermalConductivity()
thc.setup(postproc.parameters, ht)
thc.compute(ht)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[29], line 2
1 thc = ThermalConductivity()
----> 2 thc.setup(postproc.parameters, ht)
3 thc.compute(ht)
NameError: name 'postproc' is not defined
[ ]: