Peano 4
Loading...
Searching...
No Matches
TwoBodyProblem.cpph
Go to the documentation of this file.
1#include "swift2/kernels/kernel_hydro.h"
2#include "swift2/kernels/equation_of_state.h"
3
4
5template <typename ParticleContainer>
7 const ParticleContainer& localParticles,
8 const ParticleContainer& activeParticles
9) {
10 for (auto* localParticle: localParticles ) {
11 if ( not localParticle->getCellHasUpdatedParticle() ) {
12 for (auto* activeParticle: activeParticles ) {
13
14 #if Dimensions==2
15 tarch::la::Vector<Dimensions,double> pos_sun = {benchmarks::swift2::planetorbit::SUN_X_COORD,benchmarks::swift2::planetorbit::SUN_Y_COORD};
16 #else
17 tarch::la::Vector<Dimensions,double> pos_sun = {benchmarks::swift2::planetorbit::SUN_X_COORD,benchmarks::swift2::planetorbit::SUN_Y_COORD,0.0};
18 #endif
19
20 /* G_NEWTON in internal units */
21 double const L3 = benchmarks::swift2::planetorbit::LENGTH_UNIT * benchmarks::swift2::planetorbit::LENGTH_UNIT * benchmarks::swift2::planetorbit::LENGTH_UNIT;
22 double const T2 = benchmarks::swift2::planetorbit::PERIOD * benchmarks::swift2::planetorbit::PERIOD;
23 double const GN = benchmarks::swift2::planetorbit::G_NEWTON * benchmarks::swift2::planetorbit::M_SUN * T2 / L3;
24
25 tarch::la::Vector<Dimensions,double> earthSunDistance = localParticle->getX() - pos_sun;
26 const double r = tarch::la::norm2(earthSunDistance);
27 if (tarch::la::greater(r,0.0) ) {
28 const double r_inv = 1.0 / r;
29 localParticle->setA( - GN * earthSunDistance * (r_inv * r_inv * r_inv) );
30
31 /* Update Energy */
32 const double vel = tarch::la::norm2( localParticle->getV() );
33 const double E_kin = 0.5 * vel * vel;
34 const double E_pot = - GN * r_inv;
35
36 /* Update total energy (overall minus sign) */
37 localParticle->setEnergyKin( E_kin );
38 localParticle->setEnergyPot(-E_pot );
39 localParticle->setEnergyTot( - ( E_kin + E_pot ) );
40 }
41
42
43 }
44 }
45 }
46}
47
void computeGravitationalForce(const ParticleContainer &localParticles, const ParticleContainer &activeParticles)
Hard-coded two-body gravitational interaction for Earth-Sun system.
bool greater(double lhs, double rhs, double tolerance=NUMERICAL_ZERO_DIFFERENCE)
Scalar norm2(const Vector< Size, Scalar > &vector)
Computes the 2-norm of the vector, i.e.
Simple vector class.
Definition Vector.h:134