1. ( ) SPH 1. 1: 2: 2.

Size: px

Start display at page:

Download "1. ( ) SPH 1. 1: 2: 2."

きみおあみおか
5 years ago
Views:

1 SPH 2014/08/

2 1. ( ) SPH 1. 1: 2: 2.

3 1. ( )

4 MPI SIMD

5 :

6 1. MPI : 2. ( ) MPI :

7 ( ) : DRY (Don t Repeat Yourself)

8 ( ) : DRY (Don t Repeat Yourself) ( )

9 (SPH MPS MLS ) ( : Tree, FMM, PME ) ( ) ( ) ( )

10 API C++ ( ) (I/O )

11 SPH MPI

12 Recursive Multisection (JM 2004) Tofu ( )

13 1. #include<particle_simulator.h> // class ResultForce{ // public: void clear(){ // acc = 0.0; pot = 0.0; } // PS::F64vec acc; PS::F64 pot; };

14 (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; }

15 (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; }

16 (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;

17 (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;

18 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;

19 } } ( ) } force[i].acc = ai; force[i].pot = poti;

20 ( ) // (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; } }

21 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(

22 // 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); }

23 } PS::Finalize(); return 0;

24 150 ( )

25 1: 2:

26 : : v v = 0 t + (v )v = p + 1 Re 2 v CFL

27 1 2 3 CG ( )

28 : : CFL 1 : ( )

29 ( ) Explicit-MPS ( 2011, MPS S semiimplicit ) (Fan et al. 2003, Hotta et al. 2012)

30 0 1 ξ ξ

31 : : ( ) ( )

32 : : 1

33 ρ 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)

34 ξ: c: Q: 1 2 A: 1 Q: A: R R 1/3

35 ( )

36 SPH SPH ( ) ( )3

37 1000

38 1000 1/1000 1

39 2

40 ( ) ( )

41 ( ) ( )

43 SPH

44 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)

45 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)

46 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

47 SPH KH Agertz et al (MN 2007, 380,963) SPH Grid Brob test Kelvin-Helmholtz ( 3 ) SPH

48 (1) ( 1/10 10 ) 3 Grid 2 Gasoline ( 10M ) SPH

49 (2) SPH KH

50 (3) 2

51 SPH 2 ρ 2 1 = 1 m j ρ( x) W ( x x j). (8) j 1 ρ P = P ρ 2 ρ + P ρ2. (9) SPH

52 : ρ u ( ) u u u ( ) ( )

53 ( ) : f ( x) = j m j f( x j ) ρ( x j ) W ( x x j). (10) d x m j /ρ( x j )

54 (Saito and Makino 2013, Hosono et al. 2013) (SM13, Hosono et al. 2014) (Yamamoto et al. submitted) ( )

55 (1) : : =1 : =

56 (2) KHI : : =1 : =

57 ( ) SPH 0

£Ã¥×¥í¥°¥é¥ß¥ó¥°(2018) - Âè11²ó – ½ÉÂê£²¤Î²òÀâ¡¤±é½¬£² –

£Ã¥×¥í¥°¥é¥ß¥ó¥°(2018) - Âè11²ó – ½ÉÂê£²¤Î²òÀâ¡¤±é½¬£² – (2018) 11 2018 12 13 2 g v dv x dt = bv x, dv y dt = g bv y (1) b v 0 θ x(t) = v 0 cos θ ( 1 e bt) (2) b y(t) = 1 ( v 0 sin θ + g ) ( 1 e bt) g b b b t (3) 11 ( ) p14 2 1 y 4 t m y > 0 y < 0 t m1 h = 0001