Page MenuHomec4science
Files F64981554

generate_input.py
No OneTemporary
Actions

File Metadata

Created
Thu, May 30, 21:06

generate_input.py
View Options

import numpy as np
import math
import argparse
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
def is_square(apositiveint):
"""Define if a value is square or not."""
x = apositiveint // 2
seen = set([x])
while x * x != apositiveint:
x = (x + (apositiveint // x)) // 2
if x in seen:
return False
seen.add(x)
return True
parser = argparse.ArgumentParser()
parser.add_argument("number", help="number of particles", type=int)
parser.add_argument("filename", help="name of generated input file")
parser.add_argument("temperature_distribution",
help="initial temperature distribution.\nType can be: 'homogeneous', 'random'")
parser.add_argument("heat_distribution_type",
help="type of heat distribution.\nType can be: 'null', 'linear', 'circular'")
parser.add_argument("plotFlag", help="Plot heat distribution flag (True or False)")
args = parser.parse_args()
# Get number of particles
N = int(args.number)
N_squared = int(math.sqrt(N))
if not(is_square(N)):
raise ValueError(
"Number of particles must be square")
positions = np.zeros((N, 3))
velocity = positions.copy()
force = positions.copy()
masses = 1e0*(np.random.random((N, 1)) + 1)
if args.temperature_distribution == "homogeneous":
temperature = 5*np.ones((N, 1))
elif args.temperature_distribution == "random":
temperature = 1e1*np.random.random((args.number, 1))
else:
raise ValueError(
"temperature_distribution can either be 'homogeneous' or 'random'")
# Set homogeneous heat_rate for all the particles
heat_rate = 3.14*np.ones((N, 1))
# Set heat distribution
if args.heat_distribution_type == "null":
heat_distribution_square = np.zeros((N_squared, N_squared))
index = 0
for i in range(0, N_squared):
for j in range(0, N_squared):
h = -1+2*float(i)/(N_squared-1)
k = -1+2*float(j)/(N_squared-1)
index += 1
positions[index-1, :] = [h, k, 0]
elif args.heat_distribution_type == "linear":
heat_distribution_square = np.zeros((N_squared, N_squared))
index = 0
for i in range(0, N_squared):
for j in range(0, N_squared):
h = -1+2*float(i)/(N_squared-1)
k = -1+2*float(j)/(N_squared-1)
index += 1
positions[index-1, :] = [h, k, 0]
if (j == int(N_squared/4)):
heat_distribution_square[i, j] = N
elif (j == int(3*N_squared/4)):
heat_distribution_square[i, j] = -N
elif args.heat_distribution_type == "circular":
heat_distribution_square = np.zeros((N_squared, N_squared))
# Setting the Radius length
R = 1./3
index = 0
for i in range(0, N_squared):
for j in range(0, N_squared):
h = -1+2*float(i)/(N_squared-1)
k = -1+2*float(j)/(N_squared-1)
index += 1
positions[index-1, :] = [h, k, 0]
if (h**2+k**2 < R):
heat_distribution_square[i, j] = N
else:
raise ValueError(
"heat_distribution_type can be 'null', 'linear', or 'circular'")
heat_distribution_line = np.zeros((N, 1))
heat_distribution = np.zeros((N_squared, N_squared))
L = 2
freqs = np.fft.fftfreq(N_squared)*2*np.pi/L*N_squared
freqs_2d = np.einsum('i,j->ij', np.ones(N_squared), freqs)
freqs = freqs_2d**2 + freqs_2d.T**2
freqs[0, 0] = 1.
heat_distribution = np.fft.fft2(heat_distribution_square)/freqs
heat_distribution[0, 0] = 0.
heat_distribution = np.fft.ifft2(heat_distribution)
# Ploting the Heat Field
if args.plotFlag == "true":
fig = plt.figure()
ax = fig.gca(projection='3d')
X = np.linspace(-1, 1, N_squared)
X_2d = np.einsum('i,j->ij', np.ones(N_squared), X)
surf = ax.plot_surface(X_2d, X_2d.T, heat_distribution.real,
cmap=cm.coolwarm, antialiased=False)
fig.suptitle('Heat distribution ')
plt.xlabel('X')
plt.ylabel('Y')
plt.show()
heat_distribution_line = np.reshape(
heat_distribution_square, (np.product(heat_distribution_square.shape), 1))
# Save data into csv file
file_data = np.hstack((positions, velocity, force, masses,
temperature, heat_rate, heat_distribution_line))
np.savetxt(args.filename, file_data, delimiter=" ")

Event Timeline

c4science · Help