Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F62245221
compute_temperature.cc
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sat, May 11, 21:59
Size
4 KB
Mime Type
text/x-c
Expires
Mon, May 13, 21:59 (2 d)
Engine
blob
Format
Raw Data
Handle
17624100
Attached To
R7571 SP4E-TB-TL-FR
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;
}
}
/* -------------------------------------------------------------------------- */
Event Timeline
Log In to Comment