Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F73041113
sigvel.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
Thu, Jul 18, 05:42
Size
27 KB
Mime Type
text/x-c
Expires
Sat, Jul 20, 05:42 (2 d)
Engine
blob
Format
Raw Data
Handle
19133516
Attached To
rGEAR Gear
sigvel.c
View Options
/* ################################################################################## */
/* ### ### */
/* ### sigvel.c ### */
/* ### ### */
/* ### Original: hydra.c (public version of Gadget 2) ### */
/* ### Author: Volker Springel ### */
/* ### ### */
/* ### Modified: February 2011 ### */
/* ### Author: Fabrice Durier ### */
/* ### ### */
/* ################################################################################## */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
#include <gsl/gsl_math.h>
#include "allvars.h"
#include "proto.h"
/*! \file sigvel.c
* Re-Computation of SPH acceleration and maximum signal velocity
* for 'feedback' particles only.
*
* This file is a modified version of hydra.c, where the SPH forces are
* computed, and where the rate of change of entropy due to the shock heating
* (via artificial viscosity) is computed.
*/
#ifdef TIMESTEP_UPDATE_FOR_FEEDBACK
#ifdef CHIMIE_THERMAL_FEEDBACK
static
double
hubble_a
,
atime
,
hubble_a2
,
fac_mu
,
fac_vsic_fix
,
a3inv
,
fac_egy
;
#ifdef PERIODIC
static
double
boxSize
,
boxHalf
;
#ifdef LONG_X
static
double
boxSize_X
,
boxHalf_X
;
#else
#define boxSize_X boxSize
#define boxHalf_X boxHalf
#endif
#ifdef LONG_Y
static
double
boxSize_Y
,
boxHalf_Y
;
#else
#define boxSize_Y boxSize
#define boxHalf_Y boxHalf
#endif
#ifdef LONG_Z
static
double
boxSize_Z
,
boxHalf_Z
;
#else
#define boxSize_Z boxSize
#define boxHalf_Z boxHalf
#endif
#endif
/*! This function is the driver routine for the update of the calculation
* of hydrodynamical force due to the feedback impact
* from either SNe or BHs on active gas particles.
*/
void
get_sigvel
(
void
)
{
long
long
ntot
,
ntotleft
;
int
i
,
j
,
k
,
n
,
ngrp
,
maxfill
,
source
,
ndone
;
int
*
nbuffer
,
*
noffset
,
*
nsend_local
,
*
nsend
,
*
numlist
,
*
ndonelist
;
int
level
,
sendTask
,
recvTask
,
nexport
,
place
;
double
soundspeed_i
;
double
tstart
,
tend
,
sumt
,
sumcomm
;
double
timecomp
=
0
,
timecommsumm
=
0
,
timeimbalance
=
0
,
sumimbalance
;
MPI_Status
status
;
#ifdef PERIODIC
boxSize
=
All
.
BoxSize
;
boxHalf
=
0.5
*
All
.
BoxSize
;
#ifdef LONG_X
boxHalf_X
=
boxHalf
*
LONG_X
;
boxSize_X
=
boxSize
*
LONG_X
;
#endif
#ifdef LONG_Y
boxHalf_Y
=
boxHalf
*
LONG_Y
;
boxSize_Y
=
boxSize
*
LONG_Y
;
#endif
#ifdef LONG_Z
boxHalf_Z
=
boxHalf
*
LONG_Z
;
boxSize_Z
=
boxSize
*
LONG_Z
;
#endif
#endif
if
(
All
.
ComovingIntegrationOn
)
{
/* Factors for comoving integration of hydro */
hubble_a
=
All
.
Omega0
/
(
All
.
Time
*
All
.
Time
*
All
.
Time
)
+
(
1
-
All
.
Omega0
-
All
.
OmegaLambda
)
/
(
All
.
Time
*
All
.
Time
)
+
All
.
OmegaLambda
;
hubble_a
=
All
.
Hubble
*
sqrt
(
hubble_a
);
hubble_a2
=
All
.
Time
*
All
.
Time
*
hubble_a
;
fac_mu
=
pow
(
All
.
Time
,
3
*
(
GAMMA
-
1
)
/
2
)
/
All
.
Time
;
fac_egy
=
pow
(
All
.
Time
,
3
*
(
GAMMA
-
1
));
fac_vsic_fix
=
hubble_a
*
pow
(
All
.
Time
,
3
*
GAMMA_MINUS1
);
a3inv
=
1
/
(
All
.
Time
*
All
.
Time
*
All
.
Time
);
atime
=
All
.
Time
;
}
else
hubble_a
=
hubble_a2
=
atime
=
fac_mu
=
fac_vsic_fix
=
a3inv
=
fac_egy
=
1.0
;
/* `NumSphUpdate' gives the number of feedback particles on this processor that want a force update */
for
(
n
=
0
,
NumSphUpdate
=
0
;
n
<
N_gas
;
n
++
)
{
if
(
P
[
n
].
Ti_endstep
==
All
.
Ti_Current
&&
SphP
[
n
].
DeltaEgySpec
>
0
)
NumSphUpdate
++
;
for
(
j
=
0
;
j
<
3
;
j
++
)
SphP
[
n
].
FeedbackUpdatedAccel
[
j
]
=
0
;
}
numlist
=
malloc
(
NTask
*
sizeof
(
int
)
*
NTask
);
MPI_Allgather
(
&
NumSphUpdate
,
1
,
MPI_INT
,
numlist
,
1
,
MPI_INT
,
MPI_COMM_WORLD
);
for
(
i
=
0
,
ntot
=
0
;
i
<
NTask
;
i
++
)
ntot
+=
numlist
[
i
];
free
(
numlist
);
if
(
ThisTask
==
0
&&
ntot
)
printf
(
"
\t
---> Number of feedback particles = %ld
\n\n
"
,
ntot
);
noffset
=
malloc
(
sizeof
(
int
)
*
NTask
);
/* offsets of bunches in common list */
nbuffer
=
malloc
(
sizeof
(
int
)
*
NTask
);
nsend_local
=
malloc
(
sizeof
(
int
)
*
NTask
);
nsend
=
malloc
(
sizeof
(
int
)
*
NTask
*
NTask
);
ndonelist
=
malloc
(
sizeof
(
int
)
*
NTask
);
i
=
0
;
/* first particle for this task */
ntotleft
=
ntot
;
/* particles left for all tasks together */
while
(
ntotleft
>
0
)
{
for
(
j
=
0
;
j
<
NTask
;
j
++
)
nsend_local
[
j
]
=
0
;
/* do local particles and prepare export list */
tstart
=
second
();
for
(
nexport
=
0
,
ndone
=
0
;
i
<
N_gas
&&
nexport
<
All
.
BunchSizeHydro
-
NTask
;
i
++
)
if
(
P
[
i
].
Ti_endstep
==
All
.
Ti_Current
&&
SphP
[
i
].
DeltaEgySpec
>
0
)
{
ndone
++
;
for
(
j
=
0
;
j
<
NTask
;
j
++
)
Exportflag
[
j
]
=
0
;
get_sigvel_evaluate
(
i
,
0
);
for
(
j
=
0
;
j
<
NTask
;
j
++
)
{
if
(
Exportflag
[
j
])
{
for
(
k
=
0
;
k
<
3
;
k
++
)
{
HydroDataIn
[
nexport
].
Pos
[
k
]
=
P
[
i
].
Pos
[
k
];
HydroDataIn
[
nexport
].
Vel
[
k
]
=
SphP
[
i
].
VelPred
[
k
];
}
HydroDataIn
[
nexport
].
Hsml
=
SphP
[
i
].
Hsml
;
HydroDataIn
[
nexport
].
Mass
=
P
[
i
].
Mass
;
#ifdef DENSITY_INDEPENDENT_SPH
HydroDataIn
[
nexport
].
EgyRho
=
SphP
[
i
].
EgyWtDensity
;
HydroDataIn
[
nexport
].
EntVarPred
=
SphP
[
i
].
EntVarPred
;
HydroDataIn
[
nexport
].
DhsmlDensityFactor
=
SphP
[
i
].
DhsmlEgyDensityFactor
;
#else
HydroDataIn
[
nexport
].
DhsmlDensityFactor
=
SphP
[
i
].
DhsmlDensityFactor
;
#endif
HydroDataIn
[
nexport
].
Density
=
SphP
[
i
].
Density
;
#ifdef DENSITY_INDEPENDENT_SPH
HydroDataIn
[
nexport
].
Pressure
=
updated_pressure
(
SphP
[
i
].
EntropyPred
,
SphP
[
i
].
EgyWtDensity
,
SphP
[
i
].
DeltaEgySpec
);
#else
HydroDataIn
[
nexport
].
Pressure
=
updated_pressure
(
SphP
[
i
].
EntropyPred
,
SphP
[
i
].
Density
,
SphP
[
i
].
DeltaEgySpec
);
#endif
HydroDataIn
[
nexport
].
Timestep
=
P
[
i
].
Ti_endstep
-
P
[
i
].
Ti_begstep
;
/* calculation of F1 */
soundspeed_i
=
sqrt
(
GAMMA
*
updated_pressure
(
SphP
[
i
].
EntropyPred
,
SphP
[
i
].
Density
,
SphP
[
i
].
DeltaEgySpec
)
/
SphP
[
i
].
Density
);
HydroDataIn
[
nexport
].
F1
=
fabs
(
SphP
[
i
].
DivVel
)
/
(
fabs
(
SphP
[
i
].
DivVel
)
+
SphP
[
i
].
CurlVel
+
0.0001
*
soundspeed_i
/
SphP
[
i
].
Hsml
/
fac_mu
);
HydroDataIn
[
nexport
].
Index
=
i
;
HydroDataIn
[
nexport
].
Task
=
j
;
nexport
++
;
nsend_local
[
j
]
++
;
}
}
}
tend
=
second
();
timecomp
+=
timediff
(
tstart
,
tend
);
qsort
(
HydroDataIn
,
nexport
,
sizeof
(
struct
hydrodata_in
),
hydro_compare_key
);
for
(
j
=
1
,
noffset
[
0
]
=
0
;
j
<
NTask
;
j
++
)
noffset
[
j
]
=
noffset
[
j
-
1
]
+
nsend_local
[
j
-
1
];
tstart
=
second
();
MPI_Allgather
(
nsend_local
,
NTask
,
MPI_INT
,
nsend
,
NTask
,
MPI_INT
,
MPI_COMM_WORLD
);
tend
=
second
();
timeimbalance
+=
timediff
(
tstart
,
tend
);
/* now do the particles that need to be exported */
for
(
level
=
1
;
level
<
(
1
<<
PTask
);
level
++
)
{
tstart
=
second
();
for
(
j
=
0
;
j
<
NTask
;
j
++
)
nbuffer
[
j
]
=
0
;
for
(
ngrp
=
level
;
ngrp
<
(
1
<<
PTask
);
ngrp
++
)
{
maxfill
=
0
;
for
(
j
=
0
;
j
<
NTask
;
j
++
)
{
if
((
j
^
ngrp
)
<
NTask
)
if
(
maxfill
<
nbuffer
[
j
]
+
nsend
[(
j
^
ngrp
)
*
NTask
+
j
])
maxfill
=
nbuffer
[
j
]
+
nsend
[(
j
^
ngrp
)
*
NTask
+
j
];
}
if
(
maxfill
>=
All
.
BunchSizeHydro
)
break
;
sendTask
=
ThisTask
;
recvTask
=
ThisTask
^
ngrp
;
if
(
recvTask
<
NTask
)
{
if
(
nsend
[
ThisTask
*
NTask
+
recvTask
]
>
0
||
nsend
[
recvTask
*
NTask
+
ThisTask
]
>
0
)
{
/* get the particles */
MPI_Sendrecv
(
&
HydroDataIn
[
noffset
[
recvTask
]],
nsend_local
[
recvTask
]
*
sizeof
(
struct
hydrodata_in
),
MPI_BYTE
,
recvTask
,
TAG_HYDRO_A
,
&
HydroDataGet
[
nbuffer
[
ThisTask
]],
nsend
[
recvTask
*
NTask
+
ThisTask
]
*
sizeof
(
struct
hydrodata_in
),
MPI_BYTE
,
recvTask
,
TAG_HYDRO_A
,
MPI_COMM_WORLD
,
&
status
);
}
}
for
(
j
=
0
;
j
<
NTask
;
j
++
)
if
((
j
^
ngrp
)
<
NTask
)
nbuffer
[
j
]
+=
nsend
[(
j
^
ngrp
)
*
NTask
+
j
];
}
tend
=
second
();
timecommsumm
+=
timediff
(
tstart
,
tend
);
/* now do the imported particles */
tstart
=
second
();
for
(
j
=
0
;
j
<
nbuffer
[
ThisTask
];
j
++
)
get_sigvel_evaluate
(
j
,
1
);
tend
=
second
();
timecomp
+=
timediff
(
tstart
,
tend
);
/* do a block to measure imbalance */
tstart
=
second
();
MPI_Barrier
(
MPI_COMM_WORLD
);
tend
=
second
();
timeimbalance
+=
timediff
(
tstart
,
tend
);
/* get the result */
tstart
=
second
();
for
(
j
=
0
;
j
<
NTask
;
j
++
)
nbuffer
[
j
]
=
0
;
for
(
ngrp
=
level
;
ngrp
<
(
1
<<
PTask
);
ngrp
++
)
{
maxfill
=
0
;
for
(
j
=
0
;
j
<
NTask
;
j
++
)
{
if
((
j
^
ngrp
)
<
NTask
)
if
(
maxfill
<
nbuffer
[
j
]
+
nsend
[(
j
^
ngrp
)
*
NTask
+
j
])
maxfill
=
nbuffer
[
j
]
+
nsend
[(
j
^
ngrp
)
*
NTask
+
j
];
}
if
(
maxfill
>=
All
.
BunchSizeHydro
)
break
;
sendTask
=
ThisTask
;
recvTask
=
ThisTask
^
ngrp
;
if
(
recvTask
<
NTask
)
{
if
(
nsend
[
ThisTask
*
NTask
+
recvTask
]
>
0
||
nsend
[
recvTask
*
NTask
+
ThisTask
]
>
0
)
{
/* send the results */
MPI_Sendrecv
(
&
HydroDataResult
[
nbuffer
[
ThisTask
]],
nsend
[
recvTask
*
NTask
+
ThisTask
]
*
sizeof
(
struct
hydrodata_out
),
MPI_BYTE
,
recvTask
,
TAG_HYDRO_B
,
&
HydroDataPartialResult
[
noffset
[
recvTask
]],
nsend_local
[
recvTask
]
*
sizeof
(
struct
hydrodata_out
),
MPI_BYTE
,
recvTask
,
TAG_HYDRO_B
,
MPI_COMM_WORLD
,
&
status
);
/* add the result to the particles */
for
(
j
=
0
;
j
<
nsend_local
[
recvTask
];
j
++
)
{
source
=
j
+
noffset
[
recvTask
];
place
=
HydroDataIn
[
source
].
Index
;
for
(
k
=
0
;
k
<
3
;
k
++
)
SphP
[
place
].
FeedbackUpdatedAccel
[
k
]
+=
HydroDataPartialResult
[
source
].
Acc
[
k
];
if
(
SphP
[
place
].
MaxSignalVel
<
HydroDataPartialResult
[
source
].
MaxSignalVel
)
SphP
[
place
].
MaxSignalVel
=
HydroDataPartialResult
[
source
].
MaxSignalVel
;
}
}
}
for
(
j
=
0
;
j
<
NTask
;
j
++
)
if
((
j
^
ngrp
)
<
NTask
)
nbuffer
[
j
]
+=
nsend
[(
j
^
ngrp
)
*
NTask
+
j
];
}
tend
=
second
();
timecommsumm
+=
timediff
(
tstart
,
tend
);
level
=
ngrp
-
1
;
}
MPI_Allgather
(
&
ndone
,
1
,
MPI_INT
,
ndonelist
,
1
,
MPI_INT
,
MPI_COMM_WORLD
);
for
(
j
=
0
;
j
<
NTask
;
j
++
)
ntotleft
-=
ndonelist
[
j
];
}
free
(
ndonelist
);
free
(
nsend
);
free
(
nsend_local
);
free
(
nbuffer
);
free
(
noffset
);
/* collect some timing information */
MPI_Reduce
(
&
timecomp
,
&
sumt
,
1
,
MPI_DOUBLE
,
MPI_SUM
,
0
,
MPI_COMM_WORLD
);
MPI_Reduce
(
&
timecommsumm
,
&
sumcomm
,
1
,
MPI_DOUBLE
,
MPI_SUM
,
0
,
MPI_COMM_WORLD
);
MPI_Reduce
(
&
timeimbalance
,
&
sumimbalance
,
1
,
MPI_DOUBLE
,
MPI_SUM
,
0
,
MPI_COMM_WORLD
);
//if(ThisTask == 0)
// {
// All.CPU_HydCompWalk += sumt / NTask;
// All.CPU_HydCommSumm += sumcomm / NTask;
// All.CPU_HydImbalance += sumimbalance / NTask;
// }
}
/*! This function is the 'core' of the SPH force update. A target
* particle is specified which may either be local, or reside in the
* communication buffer.
* Note that only the acceleration due to feedback and the change
* of signal velocity are updated, not the actual variation of Entropy rate.
*/
void
get_sigvel_evaluate
(
int
target
,
int
mode
)
{
int
j
,
k
,
n
,
timestep
,
startnode
,
numngb
;
FLOAT
*
pos
,
*
vel
;
FLOAT
mass
,
h_i
,
dhsmlDensityFactor
,
rho
,
pressure
,
f1
,
f2
;
double
acc
[
3
],
dtEntropy
,
maxSignalVel
;
double
dx
,
dy
,
dz
,
dvx
,
dvy
,
dvz
;
double
h_i2
,
hinv
,
hinv4
;
double
p_over_rho2_i
,
p_over_rho2_j
,
soundspeed_i
,
soundspeed_j
;
double
hfc
,
dwk_i
,
vdotr
,
vdotr2
,
visc
,
mu_ij
,
rho_ij
,
vsig
;
double
h_j
,
dwk_j
,
r
,
r2
,
u
,
hfc_visc
;
int
phase
=
0
;
#ifndef NOVISCOSITYLIMITER
double
dt
;
#endif
#ifdef ART_CONDUCTIVITY
printf
(
"get_sigvel_evaluate is not implemented to run with the flag ART_CONDUCTIVITY"
);
endrun
(
674321
);
#endif
#if defined(ART_VISCO_MM) || defined(ART_VISCO_RO) || defined(ART_VISCO_CD)
printf
(
"get_sigvel_evaluate is not implemented to run with the flags ART_VISCO_MM, ART_VISCO_RO, ART_VISCO_CD"
);
endrun
(
674322
);
#endif
#ifdef DENSITY_INDEPENDENT_SPH
double
egyrho
,
entvarpred
;
#endif
if
(
mode
==
0
)
{
pos
=
P
[
target
].
Pos
;
vel
=
SphP
[
target
].
VelPred
;
h_i
=
SphP
[
target
].
Hsml
;
mass
=
P
[
target
].
Mass
;
dhsmlDensityFactor
=
SphP
[
target
].
DhsmlDensityFactor
;
rho
=
SphP
[
target
].
Density
;
#ifdef DENSITY_INDEPENDENT_SPH
pressure
=
updated_pressure
(
SphP
[
target
].
EntropyPred
,
SphP
[
target
].
EgyWtDensity
,
SphP
[
target
].
DeltaEgySpec
);
#else
pressure
=
updated_pressure
(
SphP
[
target
].
EntropyPred
,
SphP
[
target
].
Density
,
SphP
[
target
].
DeltaEgySpec
);
#endif
timestep
=
P
[
target
].
Ti_endstep
-
P
[
target
].
Ti_begstep
;
soundspeed_i
=
sqrt
(
GAMMA
*
pressure
/
rho
);
f1
=
fabs
(
SphP
[
target
].
DivVel
)
/
(
fabs
(
SphP
[
target
].
DivVel
)
+
SphP
[
target
].
CurlVel
+
0.0001
*
soundspeed_i
/
SphP
[
target
].
Hsml
/
fac_mu
);
#ifdef DENSITY_INDEPENDENT_SPH
egyrho
=
SphP
[
target
].
EgyWtDensity
;
entvarpred
=
SphP
[
target
].
EntVarPred
;
dhsmlDensityFactor
=
SphP
[
target
].
DhsmlEgyDensityFactor
;
#endif
}
else
{
pos
=
HydroDataGet
[
target
].
Pos
;
vel
=
HydroDataGet
[
target
].
Vel
;
h_i
=
HydroDataGet
[
target
].
Hsml
;
mass
=
HydroDataGet
[
target
].
Mass
;
dhsmlDensityFactor
=
HydroDataGet
[
target
].
DhsmlDensityFactor
;
rho
=
HydroDataGet
[
target
].
Density
;
pressure
=
HydroDataGet
[
target
].
Pressure
;
timestep
=
HydroDataGet
[
target
].
Timestep
;
soundspeed_i
=
sqrt
(
GAMMA
*
pressure
/
rho
);
f1
=
HydroDataGet
[
target
].
F1
;
#ifdef DENSITY_INDEPENDENT_SPH
egyrho
=
HydroDataGet
[
target
].
EgyRho
;
entvarpred
=
HydroDataGet
[
target
].
EntVarPred
;
#endif
}
/* initialize variables before SPH loop is started */
acc
[
0
]
=
acc
[
1
]
=
acc
[
2
]
=
dtEntropy
=
0
;
maxSignalVel
=
0
;
#ifdef DENSITY_INDEPENDENT_SPH
p_over_rho2_i
=
pressure
/
(
egyrho
*
egyrho
);
#else
p_over_rho2_i
=
pressure
/
(
rho
*
rho
)
*
dhsmlDensityFactor
;
#endif
h_i2
=
h_i
*
h_i
;
/* Now start the actual SPH computation for this particle */
startnode
=
All
.
MaxPart
;
do
{
numngb
=
ngb_treefind_pairs
(
&
pos
[
0
],
h_i
,
phase
,
&
startnode
);
for
(
n
=
0
;
n
<
numngb
;
n
++
)
{
j
=
Ngblist
[
n
];
dx
=
pos
[
0
]
-
P
[
j
].
Pos
[
0
];
dy
=
pos
[
1
]
-
P
[
j
].
Pos
[
1
];
dz
=
pos
[
2
]
-
P
[
j
].
Pos
[
2
];
#ifdef PERIODIC
/* find the closest image in the given box size */
if
(
dx
>
boxHalf_X
)
dx
-=
boxSize_X
;
if
(
dx
<
-
boxHalf_X
)
dx
+=
boxSize_X
;
if
(
dy
>
boxHalf_Y
)
dy
-=
boxSize_Y
;
if
(
dy
<
-
boxHalf_Y
)
dy
+=
boxSize_Y
;
if
(
dz
>
boxHalf_Z
)
dz
-=
boxSize_Z
;
if
(
dz
<
-
boxHalf_Z
)
dz
+=
boxSize_Z
;
#endif
r2
=
dx
*
dx
+
dy
*
dy
+
dz
*
dz
;
h_j
=
SphP
[
j
].
Hsml
;
if
(
r2
<
h_i2
||
r2
<
h_j
*
h_j
)
{
r
=
sqrt
(
r2
);
if
(
r
>
0
)
{
/*
here, we need to differenciate two cases, because DeltaEgySpec
is also used as a flag (see further) and may be equal to -1
*/
if
(
SphP
[
j
].
DeltaEgySpec
>
0
)
/* the particle is touch by feedback */
#ifdef DENSITY_INDEPENDENT_SPH
p_over_rho2_j
=
updated_pressure
(
SphP
[
j
].
EntropyPred
,
SphP
[
j
].
EgyWtDensity
,
SphP
[
j
].
DeltaEgySpec
)
/
(
SphP
[
j
].
EgyWtDensity
*
SphP
[
j
].
EgyWtDensity
);
#else
p_over_rho2_j
=
updated_pressure
(
SphP
[
j
].
EntropyPred
,
SphP
[
j
].
Density
,
SphP
[
j
].
DeltaEgySpec
)
/
(
SphP
[
j
].
Density
*
SphP
[
j
].
Density
);
#endif
else
#ifdef DENSITY_INDEPENDENT_SPH
p_over_rho2_j
=
SphP
[
j
].
Pressure
/
(
SphP
[
j
].
EgyWtDensity
*
SphP
[
j
].
EgyWtDensity
);
#else
p_over_rho2_j
=
SphP
[
j
].
Pressure
/
(
SphP
[
j
].
Density
*
SphP
[
j
].
Density
);
#endif
#ifdef DENSITY_INDEPENDENT_SPH
soundspeed_j
=
sqrt
(
GAMMA
*
SphP
[
j
].
Pressure
/
SphP
[
j
].
Density
);
#else
soundspeed_j
=
sqrt
(
GAMMA
*
p_over_rho2_j
*
SphP
[
j
].
Density
);
#endif
dvx
=
vel
[
0
]
-
SphP
[
j
].
VelPred
[
0
];
dvy
=
vel
[
1
]
-
SphP
[
j
].
VelPred
[
1
];
dvz
=
vel
[
2
]
-
SphP
[
j
].
VelPred
[
2
];
vdotr
=
dx
*
dvx
+
dy
*
dvy
+
dz
*
dvz
;
if
(
All
.
ComovingIntegrationOn
)
vdotr2
=
vdotr
+
hubble_a2
*
r2
;
else
vdotr2
=
vdotr
;
if
(
r2
<
h_i2
)
{
hinv
=
1.0
/
h_i
;
#ifndef TWODIMS
hinv4
=
hinv
*
hinv
*
hinv
*
hinv
;
#else
hinv4
=
hinv
*
hinv
*
hinv
/
boxSize_Z
;
#endif
u
=
r
*
hinv
;
if
(
u
<
0.5
)
dwk_i
=
hinv4
*
u
*
(
KERNEL_COEFF_3
*
u
-
KERNEL_COEFF_4
);
else
dwk_i
=
hinv4
*
KERNEL_COEFF_6
*
(
1.0
-
u
)
*
(
1.0
-
u
);
}
else
{
dwk_i
=
0
;
}
if
(
r2
<
h_j
*
h_j
)
{
hinv
=
1.0
/
h_j
;
#ifndef TWODIMS
hinv4
=
hinv
*
hinv
*
hinv
*
hinv
;
#else
hinv4
=
hinv
*
hinv
*
hinv
/
boxSize_Z
;
#endif
u
=
r
*
hinv
;
if
(
u
<
0.5
)
dwk_j
=
hinv4
*
u
*
(
KERNEL_COEFF_3
*
u
-
KERNEL_COEFF_4
);
else
dwk_j
=
hinv4
*
KERNEL_COEFF_6
*
(
1.0
-
u
)
*
(
1.0
-
u
);
}
else
{
dwk_j
=
0
;
}
if
(
soundspeed_i
+
soundspeed_j
>
maxSignalVel
)
maxSignalVel
=
soundspeed_i
+
soundspeed_j
;
if
(
vdotr2
<
0
)
/* ... artificial viscosity */
{
mu_ij
=
fac_mu
*
vdotr2
/
r
;
/* note: this is negative! */
vsig
=
soundspeed_i
+
soundspeed_j
-
3
*
mu_ij
;
if
(
vsig
>
maxSignalVel
)
maxSignalVel
=
vsig
;
if
(
P
[
j
].
Ti_endstep
==
All
.
Ti_Current
)
if
(
vsig
>
SphP
[
j
].
MaxSignalVel
)
SphP
[
j
].
MaxSignalVel
=
vsig
;
rho_ij
=
0.5
*
(
rho
+
SphP
[
j
].
Density
);
f2
=
fabs
(
SphP
[
j
].
DivVel
)
/
(
fabs
(
SphP
[
j
].
DivVel
)
+
SphP
[
j
].
CurlVel
+
0.0001
*
soundspeed_j
/
fac_mu
/
SphP
[
j
].
Hsml
);
visc
=
0.25
*
All
.
ArtBulkViscConst
*
vsig
*
(
-
mu_ij
)
/
rho_ij
*
(
f1
+
f2
);
/* .... end artificial viscosity evaluation */
#ifndef NOVISCOSITYLIMITER
/* make sure that viscous acceleration is not too large */
dt
=
imax
(
timestep
,
(
P
[
j
].
Ti_endstep
-
P
[
j
].
Ti_begstep
))
*
All
.
Timebase_interval
;
if
(
dt
>
0
&&
(
dwk_i
+
dwk_j
)
<
0
)
{
visc
=
dmin
(
visc
,
0.5
*
fac_vsic_fix
*
vdotr2
/
(
0.5
*
(
mass
+
P
[
j
].
Mass
)
*
(
dwk_i
+
dwk_j
)
*
r
*
dt
));
}
#endif
}
else
visc
=
0
;
#ifndef DENSITY_INDEPENDENT_SPH
p_over_rho2_j
*=
SphP
[
j
].
DhsmlDensityFactor
;
#endif
hfc_visc
=
0.5
*
P
[
j
].
Mass
*
visc
*
(
dwk_i
+
dwk_j
)
/
r
;
#ifdef DENSITY_INDEPENDENT_SPH
hfc
=
hfc_visc
;
/* leading-order term */
hfc
+=
P
[
j
].
Mass
*
(
dwk_i
*
p_over_rho2_i
*
SphP
[
j
].
EntVarPred
/
entvarpred
+
dwk_j
*
p_over_rho2_j
*
entvarpred
/
SphP
[
j
].
EntVarPred
)
/
r
;
/* grad-h corrections */
hfc
+=
P
[
j
].
Mass
*
(
dwk_i
*
p_over_rho2_i
*
egyrho
/
rho
*
dhsmlDensityFactor
+
dwk_j
*
p_over_rho2_j
*
SphP
[
j
].
EgyWtDensity
/
SphP
[
j
].
Density
*
SphP
[
j
].
DhsmlEgyDensityFactor
)
/
r
;
#else
hfc
=
hfc_visc
+
P
[
j
].
Mass
*
(
p_over_rho2_i
*
dwk_i
+
p_over_rho2_j
*
dwk_j
)
/
r
;
#endif
acc
[
0
]
-=
hfc
*
dx
;
acc
[
1
]
-=
hfc
*
dy
;
acc
[
2
]
-=
hfc
*
dz
;
if
(
P
[
j
].
Ti_endstep
==
All
.
Ti_Current
)
{
if
(
SphP
[
j
].
DeltaEgySpec
==
0
)
/* the particle is active but not affected by feedback */
SphP
[
j
].
DeltaEgySpec
=
-
1
;
if
(
maxSignalVel
>
SphP
[
j
].
MaxSignalVel
)
SphP
[
j
].
MaxSignalVel
=
maxSignalVel
;
SphP
[
j
].
FeedbackUpdatedAccel
[
0
]
-=
hfc
*
dx
;
SphP
[
j
].
FeedbackUpdatedAccel
[
1
]
-=
hfc
*
dy
;
SphP
[
j
].
FeedbackUpdatedAccel
[
2
]
-=
hfc
*
dz
;
}
}
}
}
fflush
(
stdout
);
}
while
(
startnode
>=
0
);
/* Now collect the result at the right place */
if
(
mode
==
0
)
{
SphP
[
target
].
MaxSignalVel
=
maxSignalVel
;
for
(
k
=
0
;
k
<
3
;
k
++
)
SphP
[
target
].
FeedbackUpdatedAccel
[
k
]
=
acc
[
k
];
}
else
{
HydroDataResult
[
target
].
MaxSignalVel
=
maxSignalVel
;
for
(
k
=
0
;
k
<
3
;
k
++
)
HydroDataResult
[
target
].
Acc
[
k
]
=
acc
[
k
];
}
}
/*! This function return the pressure as a function
* of entropy and density, but add the contribution
* of the feedback energy.
*/
FLOAT
updated_pressure
(
FLOAT
EntropyPred
,
FLOAT
Density
,
FLOAT
DeltaEgySpec
)
{
FLOAT
pressure
;
FLOAT
EgySpec
,
EgySpecUpdated
;
/* energy from entropy */
EgySpec
=
EntropyPred
/
GAMMA_MINUS1
*
pow
(
Density
*
a3inv
,
GAMMA_MINUS1
);
EgySpecUpdated
=
EgySpec
+
DeltaEgySpec
;
/* pressure */
pressure
=
GAMMA_MINUS1
*
(
Density
*
a3inv
)
*
EgySpecUpdated
;
return
pressure
;
}
/*! This function force the particle to be active
*/
void
make_particle_active
(
int
i
)
{
int
j
;
int
tstart
,
tend
;
double
dt_entr
;
double
dt_gravkick
;
double
dt_hydrokick
;
//printf("(%d) make particle %d active.\n",ThisTask,i);
#ifdef PMGRID
double
dt_gravkickA
,
dt_gravkickB
;
#endif
tstart
=
(
P
[
i
].
Ti_begstep
+
P
[
i
].
Ti_endstep
)
/
2
;
/* midpoint of old step */
tend
=
(
P
[
i
].
Ti_begstep
+
All
.
Ti_Current
)
/
2
;
/* midpoint of shrinked step */
/* kick the particle back */
kickback
(
i
,
tstart
,
tend
);
/* flag it as active */
P
[
i
].
Ti_endstep
=
All
.
Ti_Current
;
NumForceUpdate
++
;
#ifdef PMGRID
printf
(
"make_particle_active is not implemented to run with the flag PMGRID
\n
"
);
printf
(
"here we need to add the last part of advance_and_find_timesteps()
\n
"
);
endrun
(
674323
);
#endif
#ifdef MULTIPHASE
printf
(
"make_particle_active is not implemented to run with the flag MULTIPHASE
\n
"
);
endrun
(
674324
);
#endif
}
#endif
#endif
#ifdef SYNCHRONIZE_NGB_TIMESTEP
/*! In case we want to reduce the timestep of a particle,
* this function perform a negative kick on the particle.
* Positions, densities, predicted velocities and predicted
* entropies are not affected, as they are allready defined
* at the present step.
*/
void
kickback
(
int
i
,
int
tstart
,
int
tend
)
{
int
j
;
double
dt_entr
;
double
dt_gravkick
;
double
dt_hydrokick
;
double
DtEgySpec
;
double
a3inv
;
#ifdef PMGRID
double
dt_gravkickA
,
dt_gravkickB
;
!!!
this
is
probably
bad
,
please
check
!!!
if
(
All
.
ComovingIntegrationOn
)
dt_gravkickB
=
get_gravkick_factor
(
All
.
PM_Ti_begstep
,
All
.
Ti_Current
)
-
get_gravkick_factor
(
All
.
PM_Ti_begstep
,
(
All
.
PM_Ti_begstep
+
All
.
PM_Ti_endstep
)
/
2
);
else
dt_gravkickB
=
(
All
.
Ti_Current
-
(
All
.
PM_Ti_begstep
+
All
.
PM_Ti_endstep
)
/
2
)
*
All
.
Timebase_interval
;
#endif
if
(
All
.
ComovingIntegrationOn
)
{
dt_entr
=
(
tend
-
tstart
)
*
All
.
Timebase_interval
;
dt_gravkick
=
get_gravkick_factor
(
tstart
,
tend
);
dt_hydrokick
=
get_hydrokick_factor
(
tstart
,
tend
);
a3inv
=
1
/
(
All
.
Time
*
All
.
Time
*
All
.
Time
);
}
else
{
dt_entr
=
dt_gravkick
=
dt_hydrokick
=
(
tend
-
tstart
)
*
All
.
Timebase_interval
;
a3inv
=
1
;
}
for
(
j
=
0
;
j
<
3
;
j
++
)
{
P
[
i
].
Vel
[
j
]
+=
P
[
i
].
GravAccel
[
j
]
*
dt_gravkick
;
#ifdef PMGRID
P
[
i
].
Vel
[
j
]
+=
P
[
i
].
GravPM
[
j
]
*
dt_gravkickB
;
#endif
#ifdef AB_TURB
P
[
i
].
Vel
[
j
]
+=
SphP
[
i
].
TurbAccel
[
j
]
*
dt_hydrokick
;
#endif
#ifdef DISSIPATION_FORCES
P
[
i
].
Vel
[
j
]
+=
SphP
[
i
].
DissipationForcesAccel
[
j
]
*
dt_hydrokick
;
#endif
}
if
(
P
[
i
].
Type
==
0
)
/* SPH stuff */
{
for
(
j
=
0
;
j
<
3
;
j
++
)
P
[
i
].
Vel
[
j
]
+=
SphP
[
i
].
HydroAccel
[
j
]
*
dt_hydrokick
;
#ifdef COOLING
//DtEgySpec = - 1/GAMMA_MINUS1 * pow(SphP[i].Density * a3inv,GAMMA_MINUS1) * (SphP[i].DtEntropyRad);
DtEgySpec
=
SphP
[
i
].
DtEnergyRad
;
/* this is now independent of the density */
LocalSysState
.
RadiatedEnergy
+=
DtEgySpec
*
dt_entr
*
P
[
i
].
Mass
;
//printf("this is a check: dt_entr=%g RadiatedEnergy=%g\n",dt_entr,DtEgySpec * dt_entr * P[i].Mass);
#endif
SphP
[
i
].
Entropy
+=
SphP
[
i
].
DtEntropy
*
dt_entr
;
//#ifdef COOLING
// SphP[i].EntropyRad += -SphP[i].DtEntropyRad * dt_entr;
//#endif
}
#ifdef PMGRID
printf
(
"kickback is not implemented to run with the flag PMGRID
\n
"
);
printf
(
"here we need to add the last part of advance_and_find_timesteps()
\n
"
);
endrun
(
674324
);
#endif
#ifdef MULTIPHASE
printf
(
"kickback is not implemented to run with the flag MULTIPHASE
\n
"
);
endrun
(
674326
);
#endif
#ifdef LIMIT_DVEL
printf
(
"kickback is not implemented to run with the flag LIMIT_DVEL
\n
"
);
endrun
(
674327
)
#endif
}
#endif
/* SYNCHRONIZE_NGB_TIMESTEP */
Event Timeline
Log In to Comment