diff --git a/work/week12/starting_point/compute_temperature.cc b/work/week12/starting_point/compute_temperature.cc index b7b5f23..6311bea 100644 --- a/work/week12/starting_point/compute_temperature.cc +++ b/work/week12/starting_point/compute_temperature.cc @@ -1,94 +1,95 @@ #include "compute_temperature.hh" #include "fft.hh" #include "material_point.hh" #include /* -------------------------------------------------------------------------- */ ComputeTemperature::ComputeTemperature(Real dt) : dt(dt) {} void ComputeTemperature::compute(System& system) { // Get the number of material points auto N = system.getNbParticles(); // Create matrix of material points Matrix tempMtx(sqrt(N)); Matrix hv(sqrt(N)); for (auto&& entry : index(tempMtx)) { int i = std::get<0>(entry); int j = std::get<1>(entry); auto& val = std::get<2>(entry); Particle& par = system.getParticle(i*sqrt(N) + j); auto& mp = static_cast(par); val = mp.getTemperature(); // Set hv value hv(i,j) = mp.getHeatDistribution(); // std::cout << "particle(" << i*sqrt(N) + j << ")" << " "; // std::cout << "Init Temp val: " << val << std::endl; } // Apply FFT transformation Matrix tempMtx_fft = FFT::transform(tempMtx); Matrix hv_fft = FFT::transform(hv); // Get FFT coordinates Matrix> fft_coord = FFT::computeFrequencies(N); + fft_coord /= 2.; //parameters for gold Real pho= 19.32; //g/cm^3 Real C=0.129; //J/g*°C for (auto&& entry : index(tempMtx_fft)) { int i = std::get<0>(entry); int j = std::get<1>(entry); auto& val = std::get<2>(entry); auto hv_tmp = hv_fft(i,j); if (i == 0 && j == 0){ val = 0; hv_tmp = 0; } // std::cout << "Val is: " << val << ", and hv: " << hv_tmp << std::endl; // Get particle heat rate Particle& par = system.getParticle(i*sqrt(N) + j); auto& mp = static_cast(par); Real hr = mp.getHeatRate(); Real q_squared = pow(fft_coord(i,j).real(),2) + pow(fft_coord(j,i).real(),2); if (i == 0 && j == 0){ q_squared = 1; } // std::cout << "Before val is: " << val << ", " << hr << ", " << q_squared; val = (hv_tmp/q_squared - val * hr)/(pho*C); // std::cout << " After val is: " << val << std::endl; } // Inverse FFT transformation Matrix tempMtx_ifft = FFT::itransform(tempMtx_fft); // Multiply by timestep tempMtx_ifft *= this->dt; for (auto&& entry2 : index(tempMtx_ifft)) { int i = std::get<0>(entry2); int j = std::get<1>(entry2); auto& val = std::get<2>(entry2); // std::cout << " IFFT val is: " << val << std::endl; Particle& par = system.getParticle(i*sqrt(N) + j); auto& mp = static_cast(par); // Implement bourder condition // if ((i == 0 )|| (j == 0 ) || (j == sqrt(N)-1 ) || (i == sqrt(N)-1 ) ) { // mp.getTemperature() = 0; // } // else{ mp.getTemperature() += val.real(); // } // std::cout << "particle(" << i*sqrt(N) + j << ")" << " "; // std::cout << " New temp is: " << mp.getTemperature() << std::endl; } } /* -------------------------------------------------------------------------- */ diff --git a/work/week12/starting_point/generate_input.py b/work/week12/starting_point/generate_input.py index 822424a..78a64f6 100644 --- a/work/week12/starting_point/generate_input.py +++ b/work/week12/starting_point/generate_input.py @@ -1,123 +1,123 @@ 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] = 1 + heat_distribution_square[i, j] = N elif (j == int(3*N_squared/4)): - heat_distribution_square[i, j] = -1 + 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] = 1 + 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=" ") diff --git a/work/week12/starting_point/test_computeFrequency_1D.txt b/work/week12/starting_point/test_computeFrequency_1D.txt deleted file mode 100644 index 50203d9..0000000 --- a/work/week12/starting_point/test_computeFrequency_1D.txt +++ /dev/null @@ -1,10 +0,0 @@ -0.00 -6.28 -12.57 -18.85 -25.13 --31.42 --25.13 --18.85 --12.57 --6.28 diff --git a/work/week12/starting_point/test_computeFrequency_2D.txt b/work/week12/starting_point/test_computeFrequency_2D.txt deleted file mode 100644 index bf7b70f..0000000 --- a/work/week12/starting_point/test_computeFrequency_2D.txt +++ /dev/null @@ -1,10 +0,0 @@ -0.00 6.28 12.57 18.85 25.13 -31.42 -25.13 -18.85 -12.57 -6.28 -0.00 6.28 12.57 18.85 25.13 -31.42 -25.13 -18.85 -12.57 -6.28 -0.00 6.28 12.57 18.85 25.13 -31.42 -25.13 -18.85 -12.57 -6.28 -0.00 6.28 12.57 18.85 25.13 -31.42 -25.13 -18.85 -12.57 -6.28 -0.00 6.28 12.57 18.85 25.13 -31.42 -25.13 -18.85 -12.57 -6.28 -0.00 6.28 12.57 18.85 25.13 -31.42 -25.13 -18.85 -12.57 -6.28 -0.00 6.28 12.57 18.85 25.13 -31.42 -25.13 -18.85 -12.57 -6.28 -0.00 6.28 12.57 18.85 25.13 -31.42 -25.13 -18.85 -12.57 -6.28 -0.00 6.28 12.57 18.85 25.13 -31.42 -25.13 -18.85 -12.57 -6.28 -0.00 6.28 12.57 18.85 25.13 -31.42 -25.13 -18.85 -12.57 -6.28