SPH 2014/08/06 2014
1. ( ) 2. 3. 4. SPH 1. 1: 2: 2.
1. ( ) 2. 3. 4.
MPI SIMD
:
1. MPI : 2. ( ) MPI :
( ) : DRY (Don t Repeat Yourself)
( ) : DRY (Don t Repeat Yourself) ( 2014 7 )
(SPH MPS MLS ) ( : Tree, FMM, PME ) ( ) ( ) ( )
API C++ ( ) (I/O )
SPH MPI
Recursive Multisection (JM 2004) Tofu ( 2009 2012)
1. #include<particle_simulator.h> // class ResultForce{ // public: void clear(){ // acc = 0.0; pot = 0.0; } // PS::F64vec acc; PS::F64 pot; };
(1) class RealPtcl{ // public: PS::F64vec getpos() const { // return this->pos; } PS::F64 getcharge() const { // return this->mass; } void copyfromforce(const ResultForce & force){ // this->acc = force.acc; this->pot = force.pot; }
(2) // PS::S64 id; PS::F64 mass; PS::F64vec pos; PS::F64vec vel; PS::F64vec acc; PS::F64 pot; PS::S32 loadoneparticle(file * fp) { PS::S32 ret = 0; ret = fscanf(fp, "%lf%lf%lf%lf%lf%lf%lf", &this->mass, &this->pos[0], &this->pos[1], &this->pos[2 &this->vel[0], &this->vel[1], &this->vel[2 std::cout<<"this->mass"<<this->mass<<std::endl; return ret; }
(3) void dumponeparticle(file * fp){ fprintf(fp, "%lf %lf %lf %lf %lf %lf %lf this->mass, this->pos[0], this->pos[1], this->pos[2], this->vel[0], this->vel[1], this->vel[2]); } }; class EssentialPtclI{ public: void copyfromrp(const RealPtcl & rp){ // pos = rp.pos; id = rp.id; } }; // PS::F64vec pos; PS::S64 id;
(4) class EssentialPtclJ{ public: void copyfromrp(const RealPtcl & rp){ // mass = rp.mass; pos = rp.pos; id = rp.id; } }; // PS::S64 id; PS::F64 mass; PS::F64vec pos;
void CalcForceEpEp(const EssentialPtclI * ep_i, const PS::S32 n_ip, const EssentialPtclJ * ep_j, const PS::S32 n_jp, ResultForce * force){ for(ps::s32 i=0; i<n_ip; i++){ PS::F64vec xi = ep_i[i].pos; PS::F64vec ai = 0.0; PS::F64 poti = 0.0; PS::S64 idi = ep_i[i].id; for(ps::s32 j=0; j<n_jp; j++){ if( idi == ep_j[j].id ) continue; PS::F64vec rij = xi - ep_j[j].pos; PS::F64 r3_inv = rij * rij; PS::F64 r_inv = 1.0/sqrt(r3_inv); r3_inv = r_inv * r_inv; r_inv *= ep_j[j].mass; r3_inv *= r_inv; ai -= r3_inv * rij; poti -= r_inv;
} } ( ) } force[i].acc = ai; force[i].pot = poti;
( ) // (Leapfrog) void Kick(PS::ParticleSystem<RealPtcl> system, const PS::F64 dt){ RealPtcl * rp = system.getparticlepointer(); PS::S32 ni = system.getnumberofparticlelocal(); for(int i=0; i<ni; i++){ rp[i].vel = rp[i].acc * dt; } } void Drift(PS::ParticleSystem<RealPtcl> system, const PS::F64 dt){ RealPtcl * rp = system.getparticlepointer(); PS::S32 ni = system.getnumberofparticlelocal(); for(int i=0; i<ni; i++){ rp[i].pos = rp[i].vel * dt; } }
int main(int argc, char *argv[]){ PS::Initialize(argc, argv); // PS::ParticleSystem<RealPtcl> nbody_system; // nbody_system.loadparticlesingle(argv[1], "r", &RealPtcl::lo PS::DomainInfo dinfo; dinfo.initialize("domain_info.para"); dinfo.decomposedomainall(nbody_system); // PS::TreeType<ResultForce, EssentialPtclI, EssentialPtclJ>:: PS::F32 theta = 0.5; PS::S32 n_leaf_max = 8; PS::S32 n_grp_max = 64; // grav_tree.initialize(nbody_system.getnumberofparticletotal(
// LT LET GT nbody_system grav_tree.calcforceallandwriteback(calcforceepep, dinfo, nb PS::F64 dt = 0.125; PS::F64 tend = 100.0; PS::F64 tsys = 0.0; Kick(nbody_system, dt*0.5); while(tsys < tend){ tsys += dt; Drift(nbody_system, dt); if( fmod(tsys, 1.0) == 0){ // dinfo.decomposedomainall(nbody_system); } grav_tree.calcforceallandwriteback(calcforceepep, dinfo Kick(nbody_system, dt); }
} PS::Finalize(); return 0;
150 ( )
1: 2:
: : v v = 0 t + (v )v = p + 1 Re 2 v CFL
1 2 3 CG ( )
: : CFL 1 : ( )
( ) Explicit-MPS ( 2011, MPS S semiimplicit ) (Fan et al. 2003, Hotta et al. 2012)
0 1 ξ ξ
: : ( ) ( )
: : 1
ρ dv dt = p ρge z + Π ρc p dt dt = λ dp dt = (k T ) + (others) ξρ dv dt = p ρge z + 1 c Π ρc p dt dt = λ dp dt = (ck T ) + (others)
ξ: c: Q: 1 2 A: 1 Q: A: R R 1/3
( )
SPH SPH ( ) ( )3
1000
1000 1/1000 1
2
( ) ( )
( ) ( )
SPH
SPH f f ( x) = f( x )W ( x x )d x. (1) ρ( x) = j m j W ( x x j ), (2) SPH f = j m j f j ρ( x) W ( x x j). (3)
SPH (1) f : f = f 1 = j f ( x) j m j 1 ρ( x) W ( x x j). (4) m j f( x j ) ρ( x j ) W ( x x j). (5)
SPH (2) 1 P ρ 1 ρ P = P ρ 2 ρ + P ρ2. (6) v i = j m j ( Pi ρ 2 i + P j ρ 2 j ) W (x i x j ), (7) x i
SPH KH Agertz et al (MN 2007, 380,963) SPH Grid Brob test Kelvin-Helmholtz ( 3 ) SPH
(1) ( 1/10 10 ) 3 Grid 2 Gasoline ( 10M ) SPH
(2) SPH KH
(3) 2
SPH 2 ρ 2 1 = 1 m j ρ( x) W ( x x j). (8) j 1 ρ P = P ρ 2 ρ + P ρ2. (9) SPH
: ρ u ( ) u u u ( ) ( )
( ) : f ( x) = j m j f( x j ) ρ( x j ) W ( x x j). (10) d x m j /ρ( x j )
(Saito and Makino 2013, Hosono et al. 2013) (SM13, Hosono et al. 2014) (Yamamoto et al. submitted) ( )
(1) : : =1 : =
(2) KHI : : =1 : =
( ) SPH 0