Page MenuHomec4science
Files F61356283

compute_temperature.cc
No OneTemporary
Actions

File Metadata

Created
Mon, May 6, 04:07

compute_temperature.cc
View Options

#include "compute_temperature.hh"
#include "fft.hh"
#include "material_point.hh"
#include <cmath>
#include "my_types.hh"
#include "matrix.hh"
/* -------------------------------------------------------------------------- */
void ComputeTemperature::compute(System& system){
// Matrices pre-allocation
UInt Sz = system.getNbParticles();
Matrix<complex> TemperatureMatrix(Sz);
Matrix<complex> FFTTemperatureMatrix(Sz);
Matrix<complex> FFTTemp_dt(Sz);
Matrix<complex> iFFTTemp_dt(Sz);
Matrix<std::complex<Real>> qx2_qy2;
Matrix<complex> HeatRate(Sz);
Matrix<complex> FFTHeatRate(Sz);
// filling up the matrices
// To access the particels mesh we ranked them and named each particel
// from left to right of each row the first row and and first element from left will hold
// an ID of 0, the second from left will be 1 ... etc.
// for accessing this we should use i*Number_of_particles_in_one_side + j.
// i is the x-axis (rows) and j is the y-axis.
for(UInt i =0; i < std::sqrt(Sz); i++){
for(UInt j =0; j < std::sqrt(Sz); j++){
MaterialPoint &par = static_cast<MaterialPoint&>(system.getParticle(i*std::sqrt(Sz)+j));
TemperatureMatrix(i, j) = par.getTemperature();
HeatRate(i, j) = par.getHeatRate();}}
FFTTemperatureMatrix = FFT::transform(TemperatureMatrix);
FFTHeatRate = FFT::transform(HeatRate);
qx2_qy2 = FFT::computeFrequencies(Sz);
// According to the comments on the moodle!
// Option 1:
//FFTTemp_dt = FFTHeatRate - 4*this->k * std::pow(M_PI,2) / (std::pow(this->Dim,2)) * qx2_qy2 * FFTTemperatureMatrix;
// Option 2 (as described in sujet):
auto c1 = FFTTemperatureMatrix * qx2_qy2;
FFTTemp_dt = FFTHeatRate - k * (c1) / (rho * C);
//iFFTTemp_dt = FFT::itransform(FFTTemp_dt);
//TemperatureMatrix += dt * iFFTTemp_dt;
//Updating particles
for(UInt i =0; i < std::sqrt(Sz); i++){
for(UInt j =0; j < std::sqrt(Sz); j++){
MaterialPoint &par = static_cast<MaterialPoint&>(system.getParticle(i*std::sqrt(Sz)+j));
par.getTemperature() = TemperatureMatrix(i, j).real();
par.getHeatRate() = HeatRate(i, j).real();}}
}

Event Timeline

c4science · Help