Page MenuHomec4science
Files F62346076

compute_temperature.cc
No OneTemporary
Actions

File Metadata

Created
Sun, May 12, 13:53

compute_temperature.cc
View Options

#include "compute_temperature.hh"
#include "fft.hh"
#include "material_point.hh"
#include <cmath>
ComputeTemperature::ComputeTemperature(Real timestep,
Real Length,
Real rho_C,
Real kappa) : dt(timestep),
L(Length),
rhoC(rho_C),
k(kappa){}
/* -------------------------------------------------------------------------- */
void ComputeTemperature::setDeltaT(Real dt) {
this->dt = dt;
}
/* -------------------------------------------------------------------------- */
void ComputeTemperature::compute(System& system) {
UInt Nb = system.getNbParticles(); // Total num. of grid points (particles)
UInt dim = UInt(sqrt(Nb)); // Discretization in x and y directions
Matrix<complex> theta(dim); // Temperature matrix
Matrix<complex> h_v(dim); // Source term matrix
Matrix<complex> theta_fourier(dim); // Temperature in Fourier space
Matrix<complex> h_v_fourier(dim); // Source term in Fourier space
Matrix<std::complex<int>> wave_nb(dim); // Wave numbers corresponding to each node
Matrix<complex> dthetadt(dim); // Rate of change of temperature
// Initializing temperature and source term matrices from system of particles
for (auto&& par : index(theta)){
int i = std::get<0>(par);
int j = std::get<1>(par);
UInt id = i + j*dim;
auto& ic = dynamic_cast<MaterialPoint&>(system.getParticle(id));
theta(i,j) = complex(ic.getTemperature(),0);
h_v(i,j) = complex(ic.getHeatRate() ,0);
}
theta_fourier = FFT::transform(theta); // Transforming temperature
h_v_fourier = FFT::transform(h_v); // Transforming source term
wave_nb = FFT::computeFrequencies(dim); // Constructing wave numbers
// Calculating the rate of change of temperature in Fourier space
Real factor = pow(2.0*M_PI/L,2);
for(int j=0; j<dim; j++){
for(int i=0; i<dim; i++){
dthetadt(i, j) = (-k*theta_fourier(i,j)*factor*Real(norm(wave_nb(i,j)))
+h_v_fourier(i,j))/(rhoC);
}
}
// Transforming the rate of change of temperature to physical space
dthetadt = FFT::itransform(dthetadt);
// Calculating updated values of temperature
for (auto&& par : index(dthetadt)){
int i = std::get<0>(par);
int j = std::get<1>(par);
UInt id = i + j*dim;
auto& value = std::get<2>(par);
auto& BC = dynamic_cast<MaterialPoint&>(system.getParticle(id));
// Test if boundary --> zero temperature on the boundary of the domain
if (BC.getBoundary() == 1.0){
dynamic_cast<MaterialPoint&>(system.getParticle(id)).getTemperature() = 0.0;
} else {
dynamic_cast<MaterialPoint &>(system.getParticle(id)).getTemperature() += value.real() * this->dt;
}
}
}
/* -------------------------------------------------------------------------- */

Event Timeline

c4science · Help