Homec4science
Diffusion pNbody (master) Browse

pNbody/examples/movies/movie08/Gadget/src/TODO
b9f64bd8b5f6master

Last Change
Raw File

TODO
View Options

This document is not UTF8. It was detected as ISO-8859-1 (Latin 1) and converted to UTF8 for display.
multiple IMF :
----------
- add parameters
- add code
- add maxlivetime
---> !!! si toutes les étoiles explosent,
la masse de certaines étoiles devient très petite,
voir négative.
--> il faut voir les tables, regarder la perte de masse pour une masse donnée.
new cooling :
----------
--- tester le cooling avec M=4, on devrait avoir Tcool << Tdyn à faire
- ajouter implicit cooling integration ok
+ tester ok, semble bien
- verifier ComovingIntegrationOn à tester
- tester sur cluster à faire
- tester sur boite cosmo à faire
on peut enlever :
DtEgySpecRadSph
DtEntropyRadSph
new sfr :
----------
- crée des étoiles de plus petite masse, et donc, garder la part. de gaz. ok
- utiliser tous les critères, jeans aussi... à faire
- tester un model statique --> loi exp décroissante ok
- transformer toute la particule lorsque elle est trops petite ok
- ecrire les output de StP ok
- prendre en compte l'énergie des étoiles... ok
- conservation de l'énergie ok
1 proc, StarFormationNStarsFromGas = 1 snap00
4 proc, StarFormationNStarsFromGas = 1 snap01 mini diff.
1 proc, StarFormationNStarsFromGas = 4 snap02 sfr correcte, mais
!!! il y a des oscillations...
1 proc, StarFormationNStarsFromGas = 4 snap03 ok, on corrige avec TreeDomainUpdateFrequency = 0
4 proc, StarFormationNStarsFromGas = 4
- permettre : TreeDomainUpdateFrequency > 0 à faire
domain :
----------
- forcer les particules stellaire a etre arrangées dans le block 1 ok
domain.c à faire (difficile)
maxloadsph -> maxloadstars
pour PEANOHILBERT : trier les particules star à faire
peano.c à faire
domain_sumCost -> add local_DomainCountSt à faire (difficile ?)
- modifier StP ok
All.TotN_stars ok
All.MaxPartStars = All.StarFormationNStarsFromGas*All.MaxPartSph ok
+ All.PartAllocFactor * (All.TotN_stars / NTask)
allocate.c ok
All.MaxPartStars
restart.c ok
init.c ok
domain.c ok
starformation.c ok
rearrange_particle_sequence ok
io.c
read_ic.c ok
ok
- vérifier les output ok
int flag_sfr
int flag_feedback
int flag_cooling
int flag_stellarage
int flag_metals
typelist[6] ok
--> 0 si aucune particule de ce type est présente dans le block
get_particles_in_block ok
- vérifier les input ok
- partir avec un fichier qui contiend déjà des étoiles
--> conditions initiales (lit qu'une partie...) ok
--> redémarrage
- partir avec un fichier qui contiend déjà la métalicité du gaz
--> conditions initiales
--> redémarrage
- tester All.TotN_stars ok
- N_stars ok
chimie :
----------
- pour une part. d'étoile, trouver les plus proches particules de gaz ok
- il faut calculer hsml/density (gaz) pour les étoiles... ok
- dans do_chimie, commencer la boucle par la première étoile ok
et non particule
- verifier le bon nombre de voisins pour les étoiles ok
stars_density.c ok
chimie.c ok
tests :
- vérifier la densité (écrire) ok
- vérifier hsml (écrire) ok
-> on écrit dans le fichier de sortie Density,Hsml ok
- vérifier les plus proches voisins ok
--> dans chimie,
-> vérifier : Sum m/rho wij = 1, par example... ok
- initialisation
utiliser get_nelts pour avoir NELEMENTS, ou alors, check.... à faire
- pour une particule stellaire, calculer à faire
- la masse totale ejectée ok
- les métaux éjectés ok
- l'énergie éjectée ok
- injecter la masse au plus proches voisins ok
- enlever mass et elt ejecté par une étoile ok
!!! conservation de l'énergie kin,pot, int
!!! conservation de l'impulsion...
!!! conservation de la masse ok
ok, a 1%pres... a vérifier par la suite...
dans des conditions moins difficiles...
<-------- Wed Jul 22 15:09:25 CEST 2009
- mieux conserver l'énergie lors de feedback thermqie (rendre part. active ?) à faire !!! (ou feedback cinétique...)
- faire un fichier statistique qui comptes ce qui est crée... à faire
- injecter l'énergie thermique ok
- injecter l'énergie cinetique à faire
- utiliser la metallicité pour le cooling ok
- utiliser TreeDomainUpdateFrequency != 0.0 à faire
- vérifier dans chimie.c -> utilisation de vel ou velpred...
- unifier à faire
SolarAbun_Fe = 0.001771 pNbody
SolarAbun_Mg = 0.00091245 pNbody
#define FEH_SOLAR 0.001771 Gadget /* 0.00181 */
#define MGH_SOLAR 0.00091245 Gadget
FeHSolar = 0.00181 Gadget cooling !!!
- vérifier restart -> particulièrement les parametres à faire
??? All.ChimieSupernovaEnergy = all.ChimieSupernovaEnergy
- revoir sticky (faire attention au leap frog, par ex.) à faire
- cooling : on peut completement sortie docooling de timestep, non ?
!! initialisation correcte de : StarEnergyInt,StarEnergyRadSph, à faire
StarEnergyRadSticky...
--> aussi lors d'un restart !!!
--> should be in All.
--------------------------------------------------------------------------------
CPU
--------------------------------------------------------------------------------
timediff(t0, t1) = t1-t0
All.CPU_Total,
run.c : on somme les diff de la boucle principale
All.CPU_Gravity,
accel.c :
gravity_tree()
All.CPU_Hydro,
accel.c :
density()
hydro_force();
All.CPU_Domain,
domain.c :
domain_Decomposition() sans peano_hilbert_order
All.CPU_Potential, potential.c :
compute_potential()
All.CPU_Predict, accel.c :
force_update_hmax(); ???
predic.c :
move_particles()
run.c :
find_next_sync_point_and_drift()
All.CPU_TimeLine, timestep.c :
advance_and_find_timesteps()
All.CPU_Snapshot, io.c
savepositions()
All.CPU_TreeConstruction potential.c
force_treebuild()
gravtree.c
force_treebuild()
All.CPU_TreeWalk, gravtree.c !!!
All.CPU_CommSum, gravtree.c
All.CPU_Imbalance, gravtree.c
All.CPU_TreeWalk += sumt / NTask;
All.CPU_Imbalance += sumimbalance / NTask;
All.CPU_CommSum += sumcomm / NTask;
All.CPU_HydCompWalk, chimie.c density.c hydra.c stars_density.c !!!
All.CPU_HydCommSumm, chimie.c density.c hydra.c stars_density.c
All.CPU_HydImbalance, chimie.c density.c hydra.c stars_density.c
All.CPU_HydImbalance += sumimbalance / NTask;
All.CPU_HydCommSumm += sumcomm / NTask;
All.CPU_HydCompWalk += sumt / NTask;
All.CPU_EnsureNgb, density.c !!!
All.CPU_EnsureNgb += sumtimengb / NTask;
stars_density.c
All.CPU_EnsureNgb += sumtimengb / NTask;
All.CPU_PM, accel.c
long_range_force()
All.CPU_Peano
domain.c : peano_hilbert_order()
on pourait rajouter :
-------------------
All.CPU_Accel
if faut rajouter :
All.CPU_Chimie
All.CPU_StarFormation
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
- sticky bournaud
- faire des tests
- choix du sticky --> pour cluster, mieux dans rayon sph...
- éviter que la densité soit calculée ??? vraiment utiles ?
- conservation de l'énergie pour sticky
----------------------------------------
- improve multiphase computation : only with MULTIPHASE switched on
- test sticky only
* for each particle, find one particle to collide with
* try to do that in the same function than the main hydra (???)
- test sph only
- test dark gas
- test combination of all
check
* indep. of number of proc
##########################################3
MULTIPHASE
##########################################
- part of particles behave differently (sticky)
GAS_SPH -> normal sph
GAS_STICKY -> sticky
GAS_DARK -> no (weak) interaction
functions
multi_hydro_evaluate
sticky_evaluate
ngb_treefind_phase_pairs
io.c ok
init.c ok
read_ic.c ok
run.c ok
density.c ok (rien de particulier)
ngb.c
ngb_treefind_pairs hydra.c
every body sees every body
ngb_treefind_variable density.c
find density and hsml !!!! hsml is then needed to find colliding particules
!!!! but may be problematic if only sph part. are present in hsml !!!
ngb_treefind_phase_pairs bondi_accretion.c
hydra.c !!!!!!!!!!!!
multi_hydro_evaluate(i, 0);
hydro_evaluate(i, 0);
sticky_evaluate(i, 0);
timestep.c
global.c
cooling.c
starformation.c
bubbles.c
##########################################
PHASE_MIXING
##########################################
proto.h
phase_mixing
run.c
phase.c
##########################################
COLDGAS_CYCLE (need MULTIPHASE)
##########################################
- compute cycling or not
---------------------------------------
COLDGAS_CYCLE :parameters
-------------------------
ColdGasCycleTransitionTime;
ColdGasCycleTransitionParameter;
allvars.h ok
begrun.c ok
phase.c ok
##########################################
EXTERNAL_FLUX : really only external flux
##########################################
begrun.c
cooling.c
phase.c
allvars.h
HeatingPeSolarEnergyDensity
HeatingExternalFLuxEnergyDensity
##########################################
STELLAR_FLUX : only stellar flux
##########################################
allvars.h
init.c
gravtree.c
forcetree.c
begrun.c
HeatingPeSolarEnergyDensity;
HeatingPeLMRatioGas;
HeatingPeLMRatioHalo;
HeatingPeLMRatioDisk;
HeatingPeLMRatioBulge;
HeatingPeLMRatioStars;
HeatingPeLMRatioBndry;
HeatingPeLMRatio[6];
####################################################
# sticky_evaluate(i, 0)
####################################################
1) check who can inteact with who ?
- loop over active particles
- ngb_treefind_pairs ! find all particles in h_i ... check !!!
here we could use different fct, depending on the type
- if(P[j].Ti_endstep == All.Ti_Current) only active particles ! ensure symetry,
really necessary ???
- if(SphP[j].Phase == GAS_STICKY) ok, but may be done with "ngb_treefind_pairs"
- if(SphP[j].StickyFlag) SphP[i].StickyFlag = 1; in init.c
SphP[i].StickyFlag is determined in phase.c
2) what is modified by sticky_evaluate
P[target].Vel[k] <------ here, we change the velocity !!!!!!!
P[j].Vel[k]
SphP[target].EgySpecRadSticky <------ here, we count the energy !!!!!!!
SphP[j].EgySpecRadSticky
SphP[target].StickyCollisionNumber++;
SphP[j].StickyCollisionNumber++;
SphP[target].HydroAccel[k] = 0;
SphP[target].StickyFlag = 0;
SphP[target].DtEntropy = 0; /* force dt entropy to zero */
SphP[j].DtEgySpecFeedback = 0; /* should not be there */
!!!! we change the velocity and count energy further
tests:
- only sticky : be sure that particles interact symetrically...
####################################################
# problemes
####################################################
(!!) si Entropy < 0, le gaz est considéré GAS_DARK (output)
-> bien vérifier...
(!!) si une particule oscille frequemment entre dark et visible, elle risque
de colisonner trops souvent... StickyFlag=1 automatique lors du retour au sticky
####################################################
#
# comoving integration
#
####################################################
velocity
-------
gadget-1 : w=sqrt(a)*xp
gadget-2 : u=a*a*xp
################################################
# conversion (comobile -> physical')
################################################
r_p = x_c*a = Pos*a
v_p = x_c*H(a)*a + v_c * a = Pos*H(a)*a + Vel/a (in the code)
v_p = x_c*H(a)*a + v_c * a = Pos*H(a)*a + Vel*sqrt(a) (in the snapshot)
m_p = m_c = Mass
u_p = u_c (need to compute from entropy)
rho_p = rho_c / a^3 = Density/a^3
A = Ap
Pp = Pc * a^(-3 * gamma)
c_c = c_p a^(3(gamma-1))/2
################################################
# conversion (physical' -> physical)
################################################
pos_p = pos_p'/h
mass_p= mass_p'/h
t_p = t_p'/h
vel_p = vel_p'
u_p = u_p'
rho_p = rho_p'*h^2
* where :
starformation.c
---------------
1) .Density
if (SphP[i].Density*a3inv > All.ThresholdDensity)
-> All.ThresholdDensity en terme physical
2) .Vel
!!! FEEDBACK_WIND use Vel ---> bad !
chimie.c
---------------
1) .Density
LocalSysState.EnergyInt1 += P[i].Mass * SphP[i].EntropyPred / (GAMMA_MINUS1) * pow(SphP[i].Density*a3inv, GAMMA_MINUS1);
LocalSysState.EnergyInt2 += P[i].Mass * SphP[i].EntropyPred / (GAMMA_MINUS1) * pow(SphP[i].Density*a3inv, GAMMA_MINUS1);
EgySpec = SphP[i].EntropyPred / GAMMA_MINUS1 * pow(SphP[i].Density*a3inv, GAMMA_MINUS1);
DeltaEntropy = GAMMA_MINUS1*NewEgySpec/pow(SphP[i].Density*a3inv, GAMMA_MINUS1) - SphP[i].EntropyPred;
EgySpec = SphP[j].EntropyPred / GAMMA_MINUS1 * pow(SphP[j].Density*a3inv, GAMMA_MINUS1);
DeltaEntropy = GAMMA_MINUS1*NewEgySpec/pow(SphP[j].Density*a3inv, GAMMA_MINUS1) - SphP[j].EntropyPred;
2) .Vel
!!! chimie_compute_energy_kin
!!! chimie_apply_wind
!!! vj2 += SphP[j].VelPred[k]*SphP[j].VelPred[k];
!!! vi2 += vel[k]*vel[k];
timestep.c
---------------
a bien refaire
cooling.c
---------------
1) Density
CoolingForOne
DoCooling
lambda
2) Vel
no
#################################################
# steps for the entropy
predict.c
dt_drift = get_drift_factor(time0, time1);
dt_gravkick = get_gravkick_factor(time0, time1); -> pour v = a_grav * dt
dt_hydrokick = get_hydrokick_factor(time0, time1); -> pour v = a_hydro * dt
dt_entr = (time1 - (P[i].Ti_begstep + P[i].Ti_endstep) / 2) * All.Timebase_interval;
SphP[i].EntropyPred = (SphP[i].Entropy + SphP[i].DtEntropy * dt_entr);
timestep.c
dt_gravkick = get_gravkick_factor(tstart, tend);
dt_hydrokick = get_hydrokick_factor(tstart, tend);
dt_entr = (tend - tstart) * All.Timebase_interval;
SphP[i].Entropy += SphP[i].DtEntropy * dt_entr;
- que vaut DtEntropy
hydra.c
dtEntropy += 0.5 * hfc_visc * vdotr2;
SphP[i].DtEntropy *= GAMMA_MINUS1 / (hubble_a2 * pow(SphP[i].Density, GAMMA_MINUS1));
et c'est tout !
-------------------------------------------------------
Tue Dec 21 22:37:07 CET 2010
-------------------------------------------------------
dA/dt = SphP[i].DtEntropy * hubble_a
==>
dA/dt * ds/hubble_a
dans le code est ecrit par
SphP[i].DtEntropy * (tend - tstart) * All.Timebase_interval
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
j'ai suivit l'evolution de
SphP[i].DtEntropy * hubble_a
dans le cas d'un système isolé
ou du meme système en comobile,
le resultat est identique !!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
/* Conversion into program time unit */
Cp->coeff_z[2][2] = Cp->coeff_z[2][2] + u_lt;
double u_lt;
double UnitLength_in_kpc; ok
double UnitMass_in_Msol; ok
u_lt = -log10( 4.7287e11*sqrt(pow(UnitLength_in_kpc,3)/UnitMass_in_Msol));
treeasph
unit_mass 1.e10 Msol
unit_length 1. kpc
u_lt = -log10( 4.7287d11*SQRT(unit_length**3/unit_mass))
u_lt : unit_local_time ? calculé a partir de G=1 ???
u_lt = 4728700.0 unit time in yr
d'ou vient :
4.7287d11 =
time unit = 4.72e6 ans
All.G = GRAVITY / pow(All.UnitLength_in_cm, 3) * All.UnitMass_in_g * pow(All.UnitTime_in_s, 2);
GRAVITY 6.672e-8
UnitTime_in_s = sqrt ( pow(UnitLength_in_cm, 3)/GRAVITY /UnitMass_in_g )
SEC_PER_MEGAYEAR=3.155e13
UnitLength_in_cm =3.085678e+21
UnitVelocity_in_cm_per_s =100000.0
UnitTime_in_s = UnitLength_in_cm / UnitVelocity_in_cm_per_s;
UnitTime_in_Megayears = UnitTime_in_s / SEC_PER_MEGAYEAR;
u_lt = -log10(UnitTime_in_Megayears/1e6)
5.3262507580106977
c Conversion into program time unit
coeff_z(3,3) = coeff_z(3,3) + u_lt
minlivetime maxlivetime
unite treeasph : 0.828579 5.64403e+08
4.72
==> Myr : 3.9108 2.6639e+09
première étoile : 16.8457 = 79.511 Myr
première sn : 17.6926 = 83.50 Myr
_________________________________________________________________
minlivetime maxlivetime
unite gadget : 0.00400612 2.72885e+06
978.027
==> Myr 3.918 2.668e+08
6.96579
première étoile : 7.04278 = 75.28
première sn : 7.06
isolated
--------
StP[m].FormationTime : myr formation de la ssp
minlivetime : myr temps de vie minimum
maxlivetime : myr temps de vie maximum
tmin : myr temps minimum d'exection
tmax : myr temps maximum d'exection
t1 : myr debut du pas de temps
t2 : myr find du pas de temps
il y a chimie, lorsque
1) tmin < t1 && tmax >=t1 debut du pas de temps dans l'interval [tmin,tmax]
tmin <=t2 && tmax > t2 find du pas de temps dans l'interval [tmin,tmax]
tmin >=t1 && tmax <=t2 t1,t2 hord de [tmin,tmax]
starf of sfr
0.708
################################################################################
# passage unite en h -> unités physiques
################################################################################
pos_p = pos_p'/h
mass_p= mass_p'/h
t_p = t_p'/h
rho_p = rho_p'*h^2
e_p = e_p'/h
vel_p = vel_p'
u_p = u_p'
# parties a modifier :
* starformation.c
#################
All.ThresholdDensity -> unite de h
All.StarFormationDensity -> unite de h
--> rien a faire
* cooling.c
#################
--> lambda
nHcgs -> nHcgs*h^2 ok
All.mumh/All.Boltzmann ok indep de h
T= All.mumh/All.Boltzmann * GAMMA_MINUS1 * Entropy ok indep de h
l = cooling_function_with_metals ok ne dép que de T et Fe
!! l doit etre converti en unité h
l = l * h**3
l = l * nH2 si nH2 est en unité h oui
l = l * h**3
l = l * nH2h**2 si nH2 a été exprimé en unité normales ok ICI
* chimie.c
#################
* !!! toutes les grandeurs liées au temps
minlivetime = star_lifetime(StP[m].Metal[NELEMENTS-1],Cp->Mmax*SOLAR_MASS/All.UnitMass_in_g);
maxlivetime = star_lifetime(StP[m].Metal[NELEMENTS-1],Cp->Mmin*SOLAR_MASS/All.UnitMass_in_g);
-> doit etre converti en h
minlivetime = minlivetime*h
maxlivetime = maxlivetime*h
suite :
comparaison t01 t02 -> ok
# star_mass_from_age prend des valeur non h
m2 = star_mass_from_age(StP[m].Metal[NELEMENTS-1],t01/h); !!!
m1 = star_mass_from_age(StP[m].Metal[NELEMENTS-1],t02/h); !!!
m1, n2 seront en non h
Total_mass_ejection -> M0/h
EjectedMass = EjectedMass*h
le reste devrait etre bon...
__________> a vérifier... et à réfléchir...
# tests
- prendre une simulation avec unités sans h,
-> transformer en unité avec h
pos = pos*h
vel = vel
mass = mass*h
u = u
rho = rho/h**2
-> transformer les parametes (qui doivent etre en unité de h dans le code)
pos = pos*h
rho = rho/h**2
t = t/h
e = e*h
HubbleParam ok
BoxSize ok
TimeBetSnapshot ok
TimeOfFirstSnapshot ok
TimeBetStatistics ok
MaxSizeTimestep ok
MinSizeTimestep ok
All.HubbleParam
StickyDensity non -> convertir dans le code ok
StarFormationTime non -> convertir dans le code ok
StarFormationDensity non -> convertir dans le code ok
ChimieMaxSizeTimestep non -> convertir dans le code ok
ChimieWindTime non -> convertir dans le code ok
ChimieThermalTime non -> convertir dans le code ok
ChimieSupernovaEnergy non -> convertir dans le code ok
SofteningGas ok
SofteningHalo
SofteningDisk
SofteningBulge
SofteningStars
SofteningBndry
SofteningGasMaxPhys ok
SofteningHaloMaxPhys
SofteningDiskMaxPhys
SofteningBulgeMaxPhys
SofteningStarsMaxPhys
SofteningBndryMaxPhys
################################################################################
# only the master read chimie
################################################################################
pour chaque fichier chimie à lire:
---------------------------------
read_chimie(filename,nf)
-> Cps[it]
-> Elts[it]
allocation de
MassFracSNIIs[it]
EjectedMasss[it]
SingleMassFracSNIIs[it]
SingleEjectedMasss[it]
Cps
alloué par la suite
Elts
MassFracSNIIs[it] ->
EjectedMasss[it]
SingleMassFracSNIIs[it]
SingleEjectedMasss[it]
variables initializée durant read_chime
Cps
!!!
* allocate memory for elts */
Elts
!!!
/* allocate memory */
MassFracSNIIs[it] = malloc((Cps[it].nelts+2) * sizeof(double));
EjectedMasss[it] = malloc((Cps[it].nelts+2) * sizeof(double));
SingleMassFracSNIIs[it] = malloc((Cps[it].nelts+2) * sizeof(double));
SingleEjectedMasss[it] = malloc((Cps[it].nelts+2) * sizeof(double));
/* injected metals */
-- elts
c4science · Help