Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F84212811
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, Sep 21, 09:54
Size
3 KB
Mime Type
text/x-c
Expires
Mon, Sep 23, 09:54 (2 d)
Engine
blob
Format
Raw Data
Handle
20889161
Attached To
R9484 sp4e-homework-lars-bertil
compute_temperature.cc
View Options
#include "compute_temperature.hh"
#include "fft.hh"
#include "material_point.hh"
#include <cmath>
/* -------------------------------------------------------------------------- */
void
ComputeTemperature
::
setDeltaT
(
Real
dt
)
{
this
->
dt
=
dt
;
}
void
ComputeTemperature
::
compute
(
System
&
system
)
{
//std::cout << "Entered ComputeTemperature()" << std::endl;
//////// VALUES COULD ALSO BE SPECIFIED BY USER ////////////////////
Real
max_x
=
0
;
Real
min_x
=
0
;
for
(
auto
&
par
:
system
)
{
MaterialPoint
&
m_par
=
dynamic_cast
<
MaterialPoint
&>
(
par
);
if
(
m_par
.
getPosition
()[
0
]
>
max_x
)
{
max_x
=
m_par
.
getPosition
()[
0
];
}
if
(
m_par
.
getPosition
()[
0
]
<
min_x
)
{
min_x
=
m_par
.
getPosition
()[
0
];
}
}
// end for loop over particles
UInt
size
=
system
.
getNbParticles
();
UInt
N
=
std
::
sqrt
(
size
);
Real
L
=
max_x
-
min_x
;
Real
dx
=
L
/
(
N
-
1
);
// can also be extracted from distance between particle positions
Real
rho
=
1.0
;
Real
C
=
1.0
;
Real
kappa
=
1.0
;
////////////////////////////////////////////////////////////////////
Real
T
=
N
*
dx
;
// We have "size" particles aligned on a "sqrt(size) x sqrt(size)" grid
Matrix
<
complex
>
theta
(
N
);
Matrix
<
complex
>
hv
(
N
);
//UInt k = 0;
UInt
i
,
j
;
for
(
auto
&
par
:
system
)
{
MaterialPoint
&
m_par
=
dynamic_cast
<
MaterialPoint
&>
(
par
);
Vector
pos
=
m_par
.
getPosition
();
Real
x
=
pos
[
0
];
Real
y
=
pos
[
1
];
i
=
round
(((
x
+
L
/
2.0
)
/
dx
));
j
=
round
(((
y
+
L
/
2.0
)
/
dx
));
// Get 2d index from 1D index: k = i + j * N
//j = N / k;
//i = k % N;
// Fill the matrix entry (corresponding to the particle par) with its temperature and heat rate
theta
(
i
,
j
)
=
m_par
.
getTemperature
();
hv
(
i
,
j
)
=
m_par
.
getHeatRate
();
//theta(i,j) = 1.0;
//hv(i,j) = 1.0;
// update 1D index
//k++;
}
// end for loop over particles
Matrix
<
complex
>
theta_hat
=
FFT
::
transform
(
theta
);
Matrix
<
std
::
complex
<
int
>>
q_squared
=
FFT
::
computeFrequencies
(
N
);
/////// Temporary solution by converting from complex-int to complex-double
Matrix
<
std
::
complex
<
double
>>
q_squared_converted
(
N
);
for
(
auto
&&
entry
:
index
(
q_squared_converted
))
{
int
i
=
std
::
get
<
0
>
(
entry
);
int
j
=
std
::
get
<
1
>
(
entry
);
std
::
complex
<
double
>
in
=
(
std
::
real
(
q_squared
(
i
,
j
)),
std
::
imag
(
q_squared
(
i
,
j
))
);
auto
&
out
=
std
::
get
<
2
>
(
entry
);
out
=
in
;
}
//////////////////////////////////////////////////////////////////////////////
Matrix
<
complex
>
d_theta_hat_dt
(
N
);
for
(
auto
&&
entry
:
index
(
d_theta_hat_dt
))
{
int
i
=
std
::
get
<
0
>
(
entry
);
int
j
=
std
::
get
<
1
>
(
entry
);
auto
&
val
=
std
::
get
<
2
>
(
entry
);
val
=
(
1
/
(
rho
*
C
)
)
*
(
hv
(
i
,
j
)
-
kappa
*
theta_hat
(
i
,
j
)
*
(
2
*
M_PI
/
T
)
*
(
2
*
M_PI
/
T
)
*
q_squared_converted
(
i
,
j
)
);
}
Matrix
<
complex
>
d_theta_dt
=
FFT
::
itransform
(
d_theta_hat_dt
);
for
(
auto
&&
entry
:
index
(
theta
))
{
int
i
=
std
::
get
<
0
>
(
entry
);
int
j
=
std
::
get
<
1
>
(
entry
);
auto
&
val
=
std
::
get
<
2
>
(
entry
);
val
+=
dt
*
d_theta_dt
(
i
,
j
);
}
for
(
auto
&
par
:
system
)
{
MaterialPoint
&
m_par
=
dynamic_cast
<
MaterialPoint
&>
(
par
);
Vector
pos
=
m_par
.
getPosition
();
Real
x
=
pos
[
0
];
Real
y
=
pos
[
1
];
i
=
round
(((
x
+
L
/
2.0
)
/
dx
));
j
=
round
(((
y
+
L
/
2.0
)
/
dx
));
m_par
.
getTemperature
()
=
std
::
real
(
theta
(
i
,
j
));
//std::cout << "Material Point (" << i << ","<< j << ") : " << m_par << std::endl;
}
// end for loop over particles
}
/* -------------------------------------------------------------------------- */
Event Timeline
Log In to Comment