Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F88722565
ab_turb.c
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
Sun, Oct 20, 08:56
Size
10 KB
Mime Type
text/x-c
Expires
Tue, Oct 22, 08:56 (2 d)
Engine
blob
Format
Raw Data
Handle
21812012
Attached To
rGEAR Gear
ab_turb.c
View Options
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "allvars.h"
#include "proto.h"
#ifdef AB_TURB
double
st_grn
(
void
);
void
st_init_ouseq
(
void
);
void
st_calc_phases
(
void
);
void
st_compute_injection
(
void
);
/*double StEnergyAcc;
double StEnergyDeacc;
double StLastStatTime;
FILE *FdTurb;*/
void
init_turb
()
{
if
(
ThisTask
==
0
)
{
printf
(
"Init turbulence routine
\n
"
);
}
int
ikx
,
iky
,
ikz
;
double
kx
,
ky
,
kz
,
k
;
double
ampl
;
int
ikxmax
=
All
.
BoxSize
*
All
.
StKmax
/
2.
/
M_PI
;
#ifndef ONEDIMS
int
ikymax
=
All
.
BoxSize
*
All
.
StKmax
/
2.
/
M_PI
;
#ifndef TWODIMS
int
ikzmax
=
All
.
BoxSize
*
All
.
StKmax
/
2.
/
M_PI
;
#else
int
ikzmax
=
0
;
#endif
#else
int
ikymax
=
0
;
int
ikzmax
=
0
;
#endif
StNModes
=
0
;
for
(
ikx
=
0
;
ikx
<=
ikxmax
;
ikx
++
)
{
kx
=
2.
*
M_PI
*
ikx
/
All
.
BoxSize
;
for
(
iky
=
0
;
iky
<=
ikymax
;
iky
++
)
{
ky
=
2.
*
M_PI
*
iky
/
All
.
BoxSize
;
for
(
ikz
=
0
;
ikz
<=
ikzmax
;
ikz
++
)
{
kz
=
2.
*
M_PI
*
ikz
/
All
.
BoxSize
;
k
=
sqrt
(
kx
*
kx
+
ky
*
ky
+
kz
*
kz
);
if
(
k
>=
All
.
StKmin
&&
k
<=
All
.
StKmax
)
{
#if NUMDIMS ==1
StNModes
+=
1
;
#endif
#if NUMDIMS == 2
StNModes
+=
2
;
#endif
#if NUMDIMS == 3
StNModes
+=
4
;
#endif
}
}
}
}
if
(
ThisTask
==
0
)
{
printf
(
"Using %d modes, %d %d %d
\n
"
,
StNModes
,
ikxmax
,
ikymax
,
ikzmax
);
}
// StMode = (double*) mymalloc_movable(&StMode,"StModes", StNModes * 3 * sizeof(double));
// StAka = (double*) mymalloc_movable(&StAka,"StAka", StNModes * 3 * sizeof(double));
// StAkb = (double*) mymalloc_movable(&StAkb,"StAkb", StNModes * 3 * sizeof(double));
// StAmpl = (double*) mymalloc_movable(&StAmpl,"StAmpl", StNModes * sizeof(double));
// StOUPhases = (double*) mymalloc_movable(StOUPhases,"StOUPhases", StNModes * 6 * sizeof(double));
StMode
=
(
double
*
)
malloc
(
StNModes
*
3
*
sizeof
(
double
));
StAka
=
(
double
*
)
malloc
(
StNModes
*
3
*
sizeof
(
double
));
StAkb
=
(
double
*
)
malloc
(
StNModes
*
3
*
sizeof
(
double
));
StAmpl
=
(
double
*
)
malloc
(
StNModes
*
sizeof
(
double
));
StOUPhases
=
(
double
*
)
malloc
(
StNModes
*
6
*
sizeof
(
double
));
StOUVar
=
sqrt
(
All
.
StEnergy
/
All
.
StDecay
);
double
kc
=
0.5
*
(
All
.
StKmin
+
All
.
StKmax
);
double
amin
=
0.
;
#if NUMDIMS == 3
StSolWeightNorm
=
sqrt
(
3.0
/
3.0
)
*
sqrt
(
3.0
)
*
1.0
/
sqrt
(
1.0
-
2.0
*
All
.
StSolWeight
+
3.0
*
All
.
StSolWeight
*
All
.
StSolWeight
);
#endif
#if NUMDIMS == 2
StSolWeightNorm
=
sqrt
(
3.0
/
2.0
)
*
sqrt
(
3.0
)
*
1.0
/
sqrt
(
1.0
-
2.0
*
All
.
StSolWeight
+
2.0
*
All
.
StSolWeight
*
All
.
StSolWeight
);
#endif
#if NUMDIMS == 1
StSolWeightNorm
=
sqrt
(
3.0
/
1.0
)
*
sqrt
(
3.0
)
*
1.0
/
sqrt
(
1.0
-
2.0
*
All
.
StSolWeight
+
1.0
*
All
.
StSolWeight
*
All
.
StSolWeight
);
#endif
StNModes
=
0
;
for
(
ikx
=
0
;
ikx
<=
ikxmax
;
ikx
++
)
{
kx
=
2.
*
M_PI
*
ikx
/
All
.
BoxSize
;
for
(
iky
=
0
;
iky
<=
ikymax
;
iky
++
)
{
ky
=
2.
*
M_PI
*
iky
/
All
.
BoxSize
;
for
(
ikz
=
0
;
ikz
<=
ikzmax
;
ikz
++
)
{
kz
=
2.
*
M_PI
*
ikz
/
All
.
BoxSize
;
k
=
sqrt
(
kx
*
kx
+
ky
*
ky
+
kz
*
kz
);
if
(
k
>=
All
.
StKmin
&&
k
<=
All
.
StKmax
)
{
if
(
All
.
StSpectForm
==
0
)
{
ampl
=
1.
;
}
else
if
(
All
.
StSpectForm
==
1
)
{
ampl
=
4.0
*
(
amin
-
1.0
)
/
((
All
.
StKmax
-
All
.
StKmin
)
*
(
All
.
StKmax
-
All
.
StKmin
))
*
((
k
-
kc
)
*
(
k
-
kc
))
+
1.0
;
}
else
if
(
All
.
StSpectForm
==
2
)
{
ampl
=
pow
(
All
.
StKmin
,
5.
/
3
)
/
pow
(
k
,
5.
/
3
);
}
else
if
(
All
.
StSpectForm
==
3
)
{
ampl
=
pow
(
All
.
StKmin
,
2.
)
/
pow
(
k
,
2.
);
}
else
{
//terminate("unkown spectral form");
endrun
(
564321
);
}
StAmpl
[
StNModes
+
0
]
=
ampl
;
StMode
[
3
*
StNModes
+
0
]
=
kx
;
StMode
[
3
*
StNModes
+
1
]
=
ky
;
StMode
[
3
*
StNModes
+
2
]
=
kz
;
if
(
ThisTask
==
0
)
{
printf
(
"Mode: %d, ikx=%d, iky=%d, ikz=%d, kx=%f, ky=%f, kz=%f, ampl=%f
\n
"
,
StNModes
,
ikx
,
iky
,
ikz
,
kx
,
ky
,
kz
,
ampl
);
}
StNModes
++
;
#if NUMDIMS > 1
StAmpl
[
StNModes
+
0
]
=
ampl
;
StMode
[
3
*
StNModes
+
0
]
=
kx
;
StMode
[
3
*
StNModes
+
1
]
=
-
ky
;
StMode
[
3
*
StNModes
+
2
]
=
kz
;
if
(
ThisTask
==
0
)
{
printf
(
"Mode: %d, ikx=%d, iky=%d, ikz=%d, kx=%f, ky=%f, kz=%f, ampl=%f
\n
"
,
StNModes
,
ikx
,
-
iky
,
ikz
,
kx
,
-
ky
,
kz
,
ampl
);
}
StNModes
++
;
#if NUMDIMS > 2
StAmpl
[
StNModes
+
0
]
=
ampl
;
StMode
[
3
*
StNModes
+
0
]
=
kx
;
StMode
[
3
*
StNModes
+
1
]
=
ky
;
StMode
[
3
*
StNModes
+
2
]
=
-
kz
;
if
(
ThisTask
==
0
)
{
printf
(
"Mode: %d, ikx=%d, iky=%d, ikz=%d, kx=%f, ky=%f, kz=%f, ampl=%f
\n
"
,
StNModes
,
ikx
,
iky
,
-
ikz
,
kx
,
ky
,
-
kz
,
ampl
);
}
StNModes
++
;
StAmpl
[
StNModes
+
0
]
=
ampl
;
StMode
[
3
*
StNModes
+
0
]
=
kx
;
StMode
[
3
*
StNModes
+
1
]
=
-
ky
;
StMode
[
3
*
StNModes
+
2
]
=
-
kz
;
if
(
ThisTask
==
0
)
{
printf
(
"Mode: %d, ikx=%d, iky=%d, ikz=%d, kx=%f, ky=%f, kz=%f, ampl=%f
\n
"
,
StNModes
,
ikx
,
-
iky
,
-
ikz
,
kx
,
-
ky
,
-
kz
,
ampl
);
}
StNModes
++
;
#endif
#endif
}
}
}
}
StTPrev
=
-
1
;
/* StEnergyDeacc = 0;
StEnergyAcc = 0;
StLastStatTime = All.TimeBegin - All.TimeBetStatistics;*/
StRng
=
gsl_rng_alloc
(
gsl_rng_ranlxd1
);
gsl_rng_set
(
StRng
,
All
.
StSeed
);
st_init_ouseq
();
st_calc_phases
();
if
(
ThisTask
==
0
)
printf
(
"calling set_turb_ampl in init_turb
\n
"
);
set_turb_ampl
();
StTPrev
=
All
.
Ti_Current
;
}
void
st_init_ouseq
(
void
)
{
int
i
;
for
(
i
=
0
;
i
<
6
*
StNModes
;
i
++
)
{
StOUPhases
[
i
]
=
st_grn
()
*
StOUVar
;
}
}
void
st_update_ouseq
(
void
)
{
int
i
;
double
damping
=
exp
(
-
All
.
StDtFreq
/
All
.
StDecay
);
for
(
i
=
0
;
i
<
6
*
StNModes
;
i
++
)
{
StOUPhases
[
i
]
=
StOUPhases
[
i
]
*
damping
+
StOUVar
*
sqrt
(
1.
-
damping
*
damping
)
*
st_grn
();
}
}
double
st_grn
(
void
)
{
double
r0
=
gsl_rng_uniform
(
StRng
);
double
r1
=
gsl_rng_uniform
(
StRng
);
return
sqrt
(
2.
*
log
(
1.
/
r0
)
)
*
cos
(
2.
*
M_PI
*
r1
);
}
void
st_calc_phases
(
void
)
{
int
i
,
j
;
for
(
i
=
0
;
i
<
StNModes
;
i
++
)
{
double
ka
=
0.
;
double
kb
=
0.
;
double
kk
=
0.
;
int
dim
=
NUMDIMS
;
for
(
j
=
0
;
j
<
dim
;
j
++
)
{
kk
+=
StMode
[
3
*
i
+
j
]
*
StMode
[
3
*
i
+
j
];
ka
+=
StMode
[
3
*
i
+
j
]
*
StOUPhases
[
6
*
i
+
2
*
j
+
1
];
kb
+=
StMode
[
3
*
i
+
j
]
*
StOUPhases
[
6
*
i
+
2
*
j
+
0
];
}
for
(
j
=
0
;
j
<
dim
;
j
++
)
{
double
diva
=
StMode
[
3
*
i
+
j
]
*
ka
/
kk
;
double
divb
=
StMode
[
3
*
i
+
j
]
*
kb
/
kk
;
double
curla
=
StOUPhases
[
6
*
i
+
2
*
j
+
0
]
-
divb
;
double
curlb
=
StOUPhases
[
6
*
i
+
2
*
j
+
1
]
-
diva
;
StAka
[
3
*
i
+
j
]
=
All
.
StSolWeight
*
curla
+
(
1.
-
All
.
StSolWeight
)
*
divb
;
StAkb
[
3
*
i
+
j
]
=
All
.
StSolWeight
*
curlb
+
(
1.
-
All
.
StSolWeight
)
*
diva
;
}
}
}
void
set_turb_ampl
(
void
)
{
int
i
;
double
delta
=
(
All
.
Ti_Current
-
StTPrev
)
*
All
.
Timebase_interval
;
if
(
delta
>=
All
.
StDtFreq
)
{
st_update_ouseq
();
st_calc_phases
();
StTPrev
=
StTPrev
+
All
.
StDtFreq
/
All
.
Timebase_interval
;
if
(
ThisTask
==
0
)
{
printf
(
"Updated turbulent stirring field at time %f.
\n
"
,
StTPrev
*
All
.
Timebase_interval
);
}
}
}
void
add_turb_accel
()
{
int
i
,
j
,
m
;
double
acc
[
3
];
// dt_kick;
set_turb_ampl
();
double
volume
=
0
;
//for(i = FirstActiveParticle; i >= 0; i = NextActiveParticle[i])
for
(
i
=
0
;
i
<
N_gas
;
i
++
)
if
(
P
[
i
].
Ti_endstep
==
All
.
Ti_Current
)
{
if
(
P
[
i
].
Type
==
0
)
{
double
fx
=
0
;
double
fy
=
0
;
double
fz
=
0
;
for
(
m
=
0
;
m
<
StNModes
;
m
++
)
//calc force
{
double
kxx
=
StMode
[
3
*
m
+
0
]
*
P
[
i
].
Pos
[
0
];
double
kyy
=
StMode
[
3
*
m
+
1
]
*
P
[
i
].
Pos
[
1
];
double
kzz
=
StMode
[
3
*
m
+
2
]
*
P
[
i
].
Pos
[
2
];
double
kdotx
=
kxx
+
kyy
+
kzz
;
double
ampl
=
StAmpl
[
m
];
double
realt
=
cos
(
kdotx
);
double
imagt
=
sin
(
kdotx
);
fx
+=
ampl
*
(
StAka
[
3
*
m
+
0
]
*
realt
-
StAkb
[
3
*
m
+
0
]
*
imagt
);
fy
+=
ampl
*
(
StAka
[
3
*
m
+
1
]
*
realt
-
StAkb
[
3
*
m
+
1
]
*
imagt
);
fz
+=
ampl
*
(
StAka
[
3
*
m
+
2
]
*
realt
-
StAkb
[
3
*
m
+
2
]
*
imagt
);
}
fx
*=
2.
*
All
.
StAmplFac
*
StSolWeightNorm
;
//*volume;
fy
*=
2.
*
All
.
StAmplFac
*
StSolWeightNorm
;
//*volume;
fz
*=
2.
*
All
.
StAmplFac
*
StSolWeightNorm
;
//*volume;
if
(
P
[
i
].
Mass
>
0.
)
{
//FIXME correct for mean force
acc
[
0
]
=
fx
;
///P[i].Mass;
acc
[
1
]
=
fy
;
///P[i].Mass;
#ifndef TWODIMS
acc
[
2
]
=
fz
;
///P[i].Mass;
#else
acc
[
2
]
=
0
;
#endif
for
(
j
=
0
;
j
<
3
;
j
++
)
{
SphP
[
i
].
TurbAccel
[
j
]
=
acc
[
j
];
}
}
}
}
if
(
ThisTask
==
0
)
{
printf
(
"Finished turbulent accel computation.
\n
"
);
}
}
// void do_turb_driving_step_first_half(void)
// {
// CPU_Step[CPU_MISC] += measure_time();
// int i, j;
// integertime ti_step, tstart, tend;
// double dvel[3], dt_gravkick, atime;
//
//
// add_turb_accel();
//
//
// if(All.ComovingIntegrationOn)
// atime = All.Time;
// else
// atime = 1.0;
//
//
// for(i = FirstActiveParticle; i >= 0; i = NextActiveParticle[i])
// {
// ti_step = P[i].TimeBin ? (((integertime) 1) << P[i].TimeBin) : 0;
//
// tstart = P[i].Ti_begstep; /* beginning of step */
// tend = P[i].Ti_begstep + ti_step / 2; /* midpoint of step */
//
// if(All.ComovingIntegrationOn)
// dt_gravkick = get_gravkick_factor(tstart, tend);
// else
// dt_gravkick = (tend - tstart) * All.Timebase_interval;
//
// if(P[i].Type == 0)
// {
//
//
// for(j = 0; j < 3; j++)
// {
// dvel[j] = SphP[i].TurbAccel[j] * dt_gravkick;
// P[i].Vel[j] += dvel[j];
// P[i].dp[j] += P[i].Mass * dvel[j];
// }
//
//
// }
// }
//
//
// CPU_Step[CPU_DRIFT] += measure_time();
// }
//
//
// void do_turb_driving_step_second_half(void)
// {
// CPU_Step[CPU_MISC] += measure_time();
// int i, j;
// integertime ti_step, tstart, tend;
// double dvel[3], dt_gravkick, atime;
//
// if(All.ComovingIntegrationOn)
// atime = All.Time;
// else
// atime = 1.0;
//
//
// for(i = FirstActiveParticle; i >= 0; i = NextActiveParticle[i])
// {
// ti_step = P[i].TimeBin ? (((integertime) 1) << P[i].TimeBin) : 0;
//
// tstart = P[i].Ti_begstep + ti_step / 2; /* midpoint of step */
// tend = P[i].Ti_begstep + ti_step; /* end of step */
//
// if(All.ComovingIntegrationOn)
// dt_gravkick = get_gravkick_factor(tstart, tend);
// else
// dt_gravkick = (tend - tstart) * All.Timebase_interval;
//
// if(P[i].Type == 0)
// {
//
// for(j = 0; j < 3; j++)
// {
// dvel[j] = SphP[i].TurbAccel[j] * dt_gravkick;
// P[i].Vel[j] += dvel[j];
// P[i].dp[j] += P[i].Mass * dvel[j];
//
// SphP[i].VelPred[j] = P[i].Vel[j];
// }
//
//
// }
// }
//
//
// CPU_Step[CPU_DRIFT] += measure_time();
// }
#endif
Event Timeline
Log In to Comment