...
 
Commits (2)
import numpy as np
import math
import random
import matplotlib.pyplot as plt
from init import *
from comp import *
from init import *
from post import *
from init.parameters import *
def conversion_natural_units(density_ph, h_ph, T_init_ph, T_final_ph, T_delta_ph, N_argon, N_steps, N_steps_set_T):
density = density_ph*sigma**3 # Particle density (dimensionless)
h = h_ph/tau # Timestep (dimensionless)
T_init = T_init_ph*kb/epsilon # Initial temperature (dimensionless)
T_final = T_final_ph*kb/epsilon # Final temperature (dimensionless)
T_delta = T_delta_ph*kb/epsilon # Temperature step for downscaling (dimensionless)
L_box = (N_argon/density)**(1/3) # Dimensions of box/Scaling factor (dimensionless)
N_steps = N_steps + 10*N_steps_set_T # Add time steps for temperature rescaling
return(density, h, T_init, T_final, T_delta, L_box, N_steps)
def execute_calculations(density_ph, h_ph, T_init_ph, T_final_ph, T_delta_ph, N_argon, N_steps, N_steps_set_T, N_steps_av_E, init_org, plot_org):
# CONVERSION TO NATURAL UNITS
density, h, T_init, T_final, T_delta, L_box, N_steps = conversion_natural_units(density_ph, h_ph, T_init_ph, T_final_ph, T_delta_ph, N_argon, N_steps, N_steps_set_T)
# INITIAL POSITIONS & VELOCITIES
x0, v0, density = initial_state.build(N_argon, L_box, T_init, topology, init_org)
print("Calculations started...")
# TIME EVOLUTION OF POSITIONS
X, V, E_kin, U, D, diff = movements.assemble(x0, v0, N_steps, N_argon,
L_box, h, N_steps_set_T,
T_init, T_final,
T_delta, N_steps_av_E)
# OBSERVABLES
P = observables.pressure_coll(D, V, L_box, T_init, h) # Pressure (time-dependent)
P_scalar = observables.pressure_scalar(D, L_box, T_init) # Pressure (time-averaged)
CV, n_specific_heat = observables.specific_heat(E_kin, N_argon) # Heat capacity
Cv = observables.specific_heat_total(E_kin, N_argon) # Total heat capacity
G,r,r_peak = observables.pair_correlation(L_box, N_argon, D) # Pair correlation
T = 2*E_kin/((N_argon-1)*3) # Temperature
# ERROR ESTIMATION
[std_energy, mean_energy] = errorest.error_combine(E_kin + U)
[std_kinetic, mean_kinetic] = errorest.error_combine(E_kin)
[std_potential,mean_potential]= errorest.error_combine(U)
[std_pressure, mean_pressure] = errorest.error_combine(P)
[std_CV, mean_CV] = errorest.error_combine(CV/N_argon)
print(f"Energy: std = {std_energy:.2e} | mean = {mean_energy:.2e}")
print(f"Kinetic: std = {std_kinetic:.2e} | mean = {mean_kinetic:.2e}")
print(f"Potential: std = {std_potential:.2e} | mean = {mean_potential:.2e}")
print(f"Pressure: std = {std_pressure:.2e} | mean = {mean_pressure:.2e}")
print(f"CV/N_argon: std = {std_CV:.2e} | mean = {mean_CV:.2e}")
print("\n... calculations successfully finished.")
""" OUTPUT """
#Energy
plotf.energy(E_kin, U, std_energy, std_kinetic, std_potential, error_alpha, h, N_argon, plot_org)
print(f"Mean total energy = {mean_energy:.2e} [-] = {epsilon*mean_energy:.2e} [J]")
print(f"Std total energy = {std_energy:.2e} [-] = {epsilon*std_energy:.2e} [J]")
print()
# PRESSURE
plotf.pressure(P, std_pressure, error_alpha, h, density, E_kin, N_argon, plot_org)
P_conv = density_ph*kb*np.mean(T_ref*2*E_kin/((N_argon-1)*3))
print(f"Mean pressure = {mean_pressure:.2e} [-] = {P_conv*mean_pressure:.2e} [Pa]")
print(f"Std pressure = {std_pressure:.2e} [-] = {P_conv*std_pressure:.2e} [Pa]")
print(f"Pressure according to time-average formula = {P_scalar:.2e} [-] = {P_conv*P_scalar:.2e} [Pa]")
print()
# SPECIFIC HEAT
plotf.specific_heat(CV, std_CV, error_alpha, n_specific_heat, h, N_argon, plot_org)
print(f"Mean specific heat (calculated differently) = {kb*Cv/N_argon:.2e} [-] = {kb*Cv/N_argon:.2e} [J/K]")
print(f"Mean specific heat = {kb*mean_CV:.2e} [-] = {kb*mean_CV:.2e} [J/K]")
print(f"Std specific heat = {kb*std_CV:.2e} [-] = {kb*std_CV:.2e} [J/K]")
print()
# PAIR CORRELATION
plotf.pair_correlation(r, G, N_argon, plot_org)
print()
# DIFFUSION
plotf.diffusion(diff, N_argon, N_steps, h, plot_org)
print()
# TEMPERATURE
plotf.temperature(T, h, N_argon, plot_org)
print()
# VALIDATE
validation_table = validate.validation_table(density, T_init, E_kin, P_scalar, N_argon, U, Cv)
return X, L_box
\ No newline at end of file