Homec4science
Diffusion SP4E-TB-TL-FR (master) Browse

R7571/work/week12/starting_point/compute_temperature.cc
2c37794e0152master

Last Change
Raw File

compute_temperature.cc
View Options

#include <cmath>
#include "compute_temperature.hh"
#include "fft.hh"
#include "material_point.hh"
/* -------------------------------------------------------------------------- */
ComputeTemperature::ComputeTemperature(Real dt) : dt(dt) {}
void ComputeTemperature::compute(System& system) {
// Get the number of material points
auto N = system.getNbParticles();
// Create matrix of temperature and heat source in the space domain
Matrix<complex> M_theta(sqrt(N));
Matrix<complex> M_hv(sqrt(N));
// Fill in temperature and heat source matrices,
// while computing min and max positions over the x-dimension
int minI = std::numeric_limits<int>::infinity();
int maxI = -std::numeric_limits<int>::infinity();
for (auto&& entry : index(M_theta)) {
int i = std::get<0>(entry);
int j = std::get<1>(entry);
auto& theta = std::get<2>(entry);
// Get corresponding system particle (material point cast)
Particle& par = system.getParticle(i * sqrt(N) + j);
auto& mp = static_cast<MaterialPoint&>(par);
// Assign temperautre and heat distribution to corresponding matrices
theta = mp.getTemperature();
M_hv(i, j) = mp.getHeatDistribution();
// std::cout << "particle(" << i*sqrt(N) + j << ")" << " ";
// std::cout << "Init Temp val: " << theta << std::endl;
// Update space boundaries if needed
auto pos = mp.getPosition();
if (pos[0] < minI) { minI = pos[0];}
if (pos[0] > maxI) { maxI = pos[0];}
}
// Get space dimension
float spaceDimension = maxI - minI;
// Apply FFT transformation to both matrices -> Fourier domain
Matrix<complex> M_theta_hat = FFT::transform(M_theta);
Matrix<complex> M_hv_hat = FFT::transform(M_hv);
// Get coordinates in the Fourier domain
Matrix<std::complex<int>> M_q = FFT::computeFrequencies(N);
M_q /= spaceDimension; // normalize by space dimension
// parameters for gold -> to be set outside of the method
Real pho = 19.32; // g/cm^3
Real C = 0.129; // J/g*°C
// Compute time derivative of temperature distribution in the Fourier domain
Matrix<complex> M_dtheta_hat_dt(sqrt(N));
for (auto&& entry : index(M_theta_hat)) {
int i = std::get<0>(entry);
int j = std::get<1>(entry);
auto& dtheta_hat_dt = std::get<2>(entry);
auto theta_hat = M_theta_hat(i, j);
auto hv_hat = M_hv_hat(i, j);
if (i == 0 && j == 0){
theta_hat = 0;
hv_hat = 0;
}
// std::cout << "theta_hat = " << theta_hat << ", hv_hat = " << hv_hat << std::endl;
// Get particle heat rate
Particle& par = system.getParticle(i * sqrt(N) + j);
auto& mp = static_cast<MaterialPoint&>(par);
Real kappa = mp.getHeatRate();
Real laplacian_q = pow(M_q(i, j).real(), 2) + pow(M_q(j, i).real(), 2);
if (i == 0 && j == 0){
laplacian_q = 1;
}
dtheta_hat_dt = (hv_hat - kappa * theta_hat * laplacian_q) / (pho * C);
// std::cout << "dtheta_hat_dt = " << dtheta_hat_dt << std::endl;
}
// Compute inverse FFT to get time derivative of temperature distribution in the space domain
Matrix<complex> M_dtheta_dt = FFT::itransform(M_dtheta_hat_dt);
// Update the temperature of all the particles in our system
for (auto&& entry2 : index(M_dtheta_dt)) {
int i = std::get<0>(entry2);
int j = std::get<1>(entry2);
auto& dtheta_dt = std::get<2>(entry2);
// std::cout << " dtheta_dt = " << dtheta_dt << std::endl;
Particle& par = system.getParticle(i*sqrt(N) + j);
auto& mp = static_cast<MaterialPoint&>(par);
// Set null temperature at the domain boundaries
if ((i == 0 ) || (j == 0 ) || (j == (int)sqrt(N) - 1) || (i == (int)sqrt(N) - 1) ) {
mp.getTemperature() = 0;
} else {
mp.getTemperature() += dtheta_dt.real() * this->dt;
}
// std::cout << "particle(" << i * sqrt(N) + j << "): ";
// std::cout << "theta = " << mp.getTemperature() << std::endl;
}
}
/* -------------------------------------------------------------------------- */
c4science · Help