Page MenuHomec4science
Files F62365475

TH_gen.py
No OneTemporary
Actions

File Metadata

Created
Sun, May 12, 17:03

TH_gen.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Dec 13 10:53:07 2018
@author: alessia
"""
import numpy as np
import math as m
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
#=======================================================================================
#
# This program generates the initial conditions for Temperature and Heat-rate to
# be used to solve the heat-equation with zero homogeneous boundary conditions.
# The heat-rate is uniformly distributed within a given radius R and zero elsewhere
# The Temperature can either have a zero homogeneous distribution or a Gaussian distrib.
#---------------------------------------------------------------------------------------
N2=10*10 #Number of particles
R=0.5 #Radius
L=1. #half-lenght of the domain
C=1.0 #Constant val for the uniform distrib. of h
TGauss=1 #Gaussian initial distr. for T
N=int(m.sqrt(N2))
dx=2*L/N
x, y = np.mgrid[-L:(L-dx):N*1j, -L:(L-dx):N*1j]
xy = np.column_stack([x.flat, y.flat])
if(TGauss == 1):
mu = np.array([0.0, 0.0])
sigma = np.array([.15, .15])
covariance = np.diag(sigma**2)
T = multivariate_normal.pdf(xy, mean=mu, cov=covariance)
#create the trivial initial condition for T
else:
T = (np.zeros(N2))
#homogeneous distr. for heat-rate within R
h = np.zeros(N2)
h = UniformR(R, N, L, C)
z=np.column_stack((xy,T.T,h.T))
#Write the output to file N2 lines 4 columns (xpos, ypos, T, hrate)
out_file = open("particles.txt","w")
for i in range (0,N2):
out_file.write('\n')
for j in range(0,4):
out_file.write('%.20f' % z[i,j] + '\t')
out_file.close()
#----------------------Visualization of T and h distributions---------------------------
h = h.reshape(x.shape)
T = T.reshape(x.shape)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
surf=ax.plot_surface(x,y,h.T)
#surf=ax.plot_surface(x,y,T.T,alpha=1)
plt.show()
#=======================================================================================
#---------------------------------------- END ------------------------------------------
#=======================================================================================
def UniformR(R, N, L, C):
H = np.zeros((N,N))
dx = 2*L/N
#Loop over 1/4 of the domain
for i in range(0, int(N/2)):
x2 = (i*dx-L)**2
for j in range (0, int(N/2)):
y2 = (L-j*dx)**2
if(x2+y2<=R):
H[i,j] = H[i,N-j-1] = C
H[N-i-1,j] = H[N-i-1,N-j-1] = C
hh=np.reshape(H,(N*N))
return hh

Event Timeline

c4science · Help