Page MenuHomec4science
Files F76754334

compute_temperature.cc
No OneTemporary
Actions

File Metadata

Created
Sat, Aug 10, 05:44

compute_temperature.cc
View Options

#include "compute_temperature.hh"
#include "fft.hh"
#include "material_point.hh"
#include <cmath>
/* -------------------------------------------------------------------------- */
ComputeTemperature::ComputeTemperature(Real dt) : dt(dt) {}
void ComputeTemperature::setDeltaT(Real dt) {
this->dt = dt;
}
/*Compute takes a list of particles creates matrices of T and h solve for T using spectral discr. in space and finite
differences in t*/
void ComputeTemperature::compute(System& system) {
UInt N = system.getNbParticles();
UInt Msize = sqrt(N);
Real K;
Matrix<complex> Temp(Msize);
Matrix<complex> hrate(Msize);
Real L = 2.;
Real dx = L/Msize;
//Fill the matrix
for (auto & par : system) {
auto &mpoint = static_cast<MaterialPoint &>(par);
auto &T = mpoint.getTemperature();
auto &h = mpoint.getHeatRate();
auto &x = mpoint.getPosition();//vector x,y
int i = int((x[0]+L/2) / dx);
int j = int((x[1]+L/2) / dx);
Temp(i, j) = T;
hrate(i, j) = h;
}
//FFT forward
FFT matFFT;
auto Temp_sp = matFFT.transform(Temp);
auto hrate_sp = matFFT.transform(hrate);
auto freq = matFFT.computeFrequencies(Msize);
Matrix<complex> rhs_sp(Msize);
//compute right-hand-side in Spectral Space
auto q2 = abs2(freq);
Real fac = (K*4*M_PI*M_PI)/(L*L);
Temp_sp *= fac;
Temp_sp = Temp_sp*q2;
rhs_sp = hrate_sp -Temp_sp;
//FFT backward
auto rhs = matFFT.itransform(rhs_sp);
//update T
rhs *= dt;
Temp = Temp + rhs;
//update list of particles
for (auto & par : system) {
auto &mpoint = static_cast<MaterialPoint &>(par);
auto &x = mpoint.getPosition();//vector x,y
int i = int((x[0]+L/2) / dx);
int j = int((x[1]+L/2) / dx);
mpoint.getTemperature()= Temp(i,j).real();
}
}
/* -------------------------------------------------------------------------- */

Event Timeline

c4science · Help