From a01f8dcf98cd8bb234587ae091dc9f1280de8d59 Mon Sep 17 00:00:00 2001 From: christosneg Date: Thu, 2 Jan 2025 15:36:37 +0200 Subject: [PATCH 1/8] Add demo showcasing custom CelestialObject agent and Gravity behavior This commit introduces a new demo that includes: - A custom agent class `CelestialObject` to represent astrophysical objects such as stars, planets, and satellites. - A custom behavior class `Gravity` to simulate gravitational forces and their effects on celestial objects. - A simple demonstration `NewtonsLawTest` to validate the interaction between agents and gravitational physics. The purpose of this demo is to: - Demonstrate the creation of custom agents and behaviors in BioDynaMo. - Highlight the platform's ability to simulate environments beyond its original design. --- demo/newtons_law_test/CMakeLists.txt | 35 +++ demo/newtons_law_test/NewtonsLaw.png | Bin 0 -> 4581 bytes demo/newtons_law_test/README.md | 49 +++++ demo/newtons_law_test/bdm.toml | 6 + demo/newtons_law_test/src/CelestialObjects.h | 199 ++++++++++++++++++ demo/newtons_law_test/src/Gravity.h | 116 ++++++++++ demo/newtons_law_test/src/newtons_law_test.cc | 18 ++ demo/newtons_law_test/src/newtons_law_test.h | 68 ++++++ 8 files changed, 491 insertions(+) create mode 100644 demo/newtons_law_test/CMakeLists.txt create mode 100644 demo/newtons_law_test/NewtonsLaw.png create mode 100644 demo/newtons_law_test/README.md create mode 100644 demo/newtons_law_test/bdm.toml create mode 100644 demo/newtons_law_test/src/CelestialObjects.h create mode 100644 demo/newtons_law_test/src/Gravity.h create mode 100644 demo/newtons_law_test/src/newtons_law_test.cc create mode 100644 demo/newtons_law_test/src/newtons_law_test.h diff --git a/demo/newtons_law_test/CMakeLists.txt b/demo/newtons_law_test/CMakeLists.txt new file mode 100644 index 000000000..e09f08e2b --- /dev/null +++ b/demo/newtons_law_test/CMakeLists.txt @@ -0,0 +1,35 @@ +# ----------------------------------------------------------------------------- +# +# Copyright (C) 2021 CERN & University of Surrey for the benefit of the +# BioDynaMo collaboration. All Rights Reserved. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# +# See the LICENSE file distributed with this work for details. +# See the NOTICE file distributed with this work for additional information +# regarding copyright ownership. +# +# ----------------------------------------------------------------------------- +cmake_minimum_required(VERSION 3.19.3) +project(newtons_law_test) + +# BioDynaMo curretly uses the C++17 standard. +set(CMAKE_CXX_STANDARD 17) + +# Use BioDynaMo in this project. +find_package(BioDynaMo REQUIRED) + +# See UseBioDynaMo.cmake in your BioDynaMo build folder for details. +# Note that BioDynaMo provides gtest header/libraries in its include/lib dir. +include(${BDM_USE_FILE}) + +# Consider all files in src/ for BioDynaMo simulation. +include_directories("src") +file(GLOB_RECURSE PROJECT_HEADERS src/*.h) +file(GLOB_RECURSE PROJECT_SOURCES src/*.cc) + +bdm_add_executable(${CMAKE_PROJECT_NAME} + HEADERS ${PROJECT_HEADERS} + SOURCES ${PROJECT_SOURCES} + LIBRARIES ${BDM_REQUIRED_LIBRARIES}) diff --git a/demo/newtons_law_test/NewtonsLaw.png b/demo/newtons_law_test/NewtonsLaw.png new file mode 100644 index 0000000000000000000000000000000000000000..b6a02207511c8dffd20e28373c24da86994d43b4 GIT binary patch literal 4581 zcmdT|XHXN`w#H)thSiyz^lWnbK{&h8STkt7k4D#A+f_TKR~V#DIW zf%ikj!a=byq37d*?}vth;v!<>*?YUKMNa0C?%wgX_d^5YBV$1)+#_Q`MMg0(!dpT! z@$m@*vlE0;LN}rHgpq;Kg%bwmKw}f2k;qs=C^=zo@2(utCMPB)ZRBis$sJO{Bta3r z;SYP}`HaRm!;U#QnY#?FJ?2pdB`(=T{Ij*g$>ZcVPx%k;T=Lgl_RGn&>YRC(su3}E z@tSV_-N^Em>Wb^q@6KO1+e|8LJ*Dx-OFvLRTC4zY* zm|FB7#^GUra0edoO?O;_%h}UnPQf9Qx?(-T%)5I}Y0iobNB}Mx^&WnK_`5)3Q-=Jb zKY2r&yGQ|_q6}ptbu}~y9XnJr_SefnUtV#kEpY%`B>m3uwd1|b9@Tk9jQmG`u8I6N zzlSzu)}gn9Pt~WUo*CZN{4GJv(XQkHtHx}jwQ3oQ&>SWSaaJ6NfXI;(1<1Dv#xYVaj-Fl{iAEz&OEqrrB z-rOitdy`5et+q|JCpebKY2N$Qy}8jC!y;{|JFZutzp1|%QK}`{rG@*!?sc)6tmSG; z1v{tV`bSUBzQ6C5?S{qaL1V74NfZTUPK0=36(*pSt#7&v3P(kq(>3(cSg61bImrim zp^cIS0XCn`u@5pldh0iO+L)YXLJ}{{P-dkfz4@Ug)4vDGpJ^ZNpBQ$XgSOXmljJKu z$czS{ba|ke=AGD=CIhL(9`pfD+UE8OT)V$LCVRtBM%!Lq8|#Fyk9|T?y9Z`@W#)kC z;qH4YS9lYkymoR3-<{w_4-w)vN}==Gi64v1qsF!a2xGKMuCy9bpppW$!!+2lmZM; zOOylqF31moUHpYRp0X#)2yPn^xf7L4>^l6)BQ`zsf<{?~prPi9_kD-iD*ooCUDD1L z0GHy^j{8)fW7c6fqQ)Gu&o*_K#ed_Oy7C;T8L@8>EUm9$J;jx**JtigYv_IYtj(^o zvab+#bIIy{+i{$paf@IW0z2RVb?qb{HCJRM0;`S}0_*ta>W(RvPq^6ST?d9XX}B!) zC+Pv9h4G*d`0&1_7@>Bj#2=@wF zxc-}$gRK^rA7q&5-!SVOxAR2tFZJ|8W>qE_?rucB=aBsvpzj*4HGPDn-)V;yK8R#A zkSt2VR=+n6Z%$iI*6d2YjV$k0EYjrg9-B*pDUnz_=j3eF2D6i$=2?q|X)eeZdTofK z1a5=Cq0l+sQyGQ67184 zBeLTiHBfOB53oCl0D@sd1n=T2{I@E|ZIh{;4|iWKtSuZEX6=Brr*90WD+DaKn5(FZ z5^SqY8^;o}6N`eXjF+=qOI8fB9r)?#lYH>IWoD! z=^2e(Ns-nP<0|u7AJRxrDngN&gW7{x<3Bm~~`ApgU) z>*r&$kl~eRXb?2Zg6fn9A^+6(?R86OwL>~bQp-uC%q*BDGj-ALXt`aUj*G|E?|$N4 zR#!(PA3&q3f&K2zv1*-1k_Aa|SS#gfWBU<{;JJ`I2+Q0C6|XK4oqjiWWS0aDs`~!< zz2l4rBEHsA164`AZm;Wta39X(I^`)cf3(ZvN6~uL)T~A3FO$C3NUwo8@~KbPs7@`^ z_i;~)X|^IJM8dX(vOHAz7e}DxBH5|K3ey-&)N*@i0N zi85Y9rw#eEs|}btl63UFbMIA`nQlllxTvy>0%E0*R1Wg>-hRu7J6rnJ%Om(pNL+o$ zZscdk7$qS#zbjQmioL_wLYBnhW5l28LnB;QT|b$NtxUb3_sLd z=6xTCuya3+nkapbDn649lP~eSnN(l(@N^0JJoTwhn~G5CvB}^`c2)cs(N{s5dE2$X zQ`9JZiN&8)l6=|ohGxQ9YbUK~@vnJ=XMKc;seT%Zn|b!d=k1iwWoa_cFYEgAWEBo` zMuS=nsan+(!KY#x6b#orlI3P+2ispPQ+-eFPw(%Vtc8u!P)&BuB1N*g{TdfOqx!L< z3tnZ@d9kOqL57Z5#AWjxdyX#@#t8&%+TrA*lM(xLL00a zvWOJ3o5wOAT=2j4M!Q;xQFq4B;}*w#sBav{>yHe-Dp01k@)cLG z?day_m{UQ+k_Wsoo=oA_iep?&Dy$T_mBVB=sL`!^EqR4{9}1Nx{Z_|f zH=l<7RPnYLy7x)P(g>rUyy7u|rO{-W9gHZ>m+{MVo6{Qggp5RJ(YjuZ0mPx_v&|W< zTEXrF!6D7dmw$s(`$qVj#?>mfMcy<~Q=%wz3+kh(bDMSuPl4B`G(`#UdW}aBC+n$y zK%zJr?6wVDO%dPQ9mwrRF3GlIpL+Rmgz#X4i^!WEn2wDs$hOHuD_>D56-=4hIgmWz zl`_^C;IquG>ke$z{jpvFn_kG|a8v`*jRB8_rPofI|JZ%|EeVhTwxV%Wnc$#g?%ru z^Ij}@w4u8kp}w+nTG=!f0e*cv@U;Y>VWg6_O}*NLGCh*T_PKq(a!C29etZIdtD~>i zPu~W{?ccmydFNG5os4HS9R#!>a|V1+6G?YWKyik5JvJBp=GrkY33Jh3yV#L#nO#4Y z2knbP!?F1K2i#)vF$y`^%M!wdNm|<=iz&ZqvPU4k6q{h~6i2+qK~ZKt3=df)Z&RWC z9mgant{Ye8$H!r9rTTXJ^QoF!Zbl}GoGp{T*t#*0q%+0k(xCA!ffO@M&Wdz67`orPSPdFvw7c+$-=9Kz9<*+7~_eSvw0mBV28RgsnYMv%IcpuS=|6}WG z7Bo0&+3CtvYgOY>N9tluZP3-tfuw|?)#5X&QtwaQ^wd6d3?X~t1$&3bMW5! zs;nU}F=z-+z@lW`(x!5nMQIVKsdWR|v^B3Xz>n4^`Ty46TE#s>*MAH;2pD8FHo*bp&0==;^;@LS>K2W9 zokk3!k27ooryrd-;%Pk}b_SG4bfm6OcG#00(cqHqKf53U`{{keZq`@iW`%+jbnJGY zEa*`P!{e8~k2d>PyYJtiMQy)SW~-(v$D&O9A2MsHR`8^mk6z-(F8rk`^MCY%lYy8x zN(kbEkNQlvsWoHl@&1E=3oQ#r`pRp0SmzaDdgaPaN?!f|`$`1zCu{65f9M86X0GVv~ z#g_h0XU{GQgd2dh63TN*=yy`XkHAf=SIN=O(Jlt=z4&x+$bF+Mm1G^Y$X04w8cbRw z6k}>!yg8s`N>hJ!iM1@z^5;e0#Nnkc{;G<==UutIU&PB&ZGO(x()gex+meKZXMs&e s{*f + Simulation + + +## Contents + +- **Celestial Object Agents** +- **Gravity Behavior** + +## Assumptions + +To simplify the simulation, the following assumptions were made: + - All celestial objects (planets, stars, and moons) are modeled as **spheres**. + - Celestial objects are treated as **point masses**. + - Agents can't physically interact with each other. + +## Newton's Law of Universal Gravitation + +Newton’s law of universal gravitation governs the movement of all celestial objects in the simulation, describing the force of attraction between two masses. + +The force of attraction between two objects is proportional to both their masses and inversely proportional to the square of the distance between them. + +$$F = \frac{G * (m_1 * m_2)}{r^2}$$ + +where: +- $F$ is the gravitational force between the two masses, +- $G$ is the gravitational constant, +- $m_1$ and $m_2$ are the masses of the two objects, and +- $r$ is the distance between the centers of the two masses. + +By knowing the force of attraction and therefore the acceleration of a Celestial Object, the displacement can be predicted by solving the below differential equation: + +$$\frac{d^2x}{dt^2} = a$$ + +--- + +## Run the demos + +In order to run the demos, BioDynamo has to be correctly installed and sourced. + +```bash +cd [path_to_biodynamo]demos/NewtonsLawTest +bdm run +bdm view +``` diff --git a/demo/newtons_law_test/bdm.toml b/demo/newtons_law_test/bdm.toml new file mode 100644 index 000000000..2898a4993 --- /dev/null +++ b/demo/newtons_law_test/bdm.toml @@ -0,0 +1,6 @@ +[visualization] +export = true +interval = 1 + +[[visualize_agent]] +name = "CelestialObject" \ No newline at end of file diff --git a/demo/newtons_law_test/src/CelestialObjects.h b/demo/newtons_law_test/src/CelestialObjects.h new file mode 100644 index 000000000..10cddf869 --- /dev/null +++ b/demo/newtons_law_test/src/CelestialObjects.h @@ -0,0 +1,199 @@ +// ----------------------------------------------------------------------------- +// +// Copyright (C) 2021 CERN & University of Surrey for the benefit of the +// BioDynaMo collaboration. All Rights Reserved. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// +// See the LICENSE file distributed with this work for details. +// See the NOTICE file distributed with this work for additional information +// regarding copyright ownership. +// +// ----------------------------------------------------------------------------- +#ifndef CELESTIALOBJECTS_H_ +#define CELESTIALOBJECTS_H_ + +#include "biodynamo.h" + +namespace bdm { +namespace astrophysics { + +// creating the custom agent that +// defines celestial objects +class CelestialObject : public Agent { + BDM_AGENT_HEADER(CelestialObject, Agent, 1); + + public: + // constructors + CelestialObject() : diameter_(1.0), density_(1.0) { UpdateVolume(); } + + explicit CelestialObject(real_t diameter) + : diameter_(diameter), density_(1.0) { + UpdateVolume(); + } + + explicit CelestialObject(const Real3& position) + : position_(position), diameter_(1.0), density_(1.0) { + UpdateVolume(); + } + + // default destructor + ~CelestialObject() override = default; + + Real3 GetVelocity() const { return velocity_; } + + real_t GetDiameter() const override { return diameter_; } + + real_t GetMass() const { return density_ * volume_; } + + real_t GetDensity() const { return density_; } + + const Real3& GetPosition() const override { return position_; } + + real_t GetVolume() const { return volume_; } + + Shape GetShape() const override { return Shape::kSphere; } + + void SetVelocity(const Real3& vel) { velocity_ = vel; } + + void SetDiameter(real_t diameter) override { + if (diameter > diameter_) { + SetPropagateStaticness(); + } + diameter_ = diameter; + UpdateVolume(); + } + + void SetVolume(real_t volume) { + volume_ = volume; + UpdateDiameter(); + } + + void SetMass(real_t mass) { SetDensity(mass / volume_); } + + void SetDensity(real_t density) { + if (density > density_) { + SetPropagateStaticness(); + } + density_ = density; + } + + void SetPosition(const Real3& position) override { + position_ = position; + SetPropagateStaticness(); + } + + void ChangeVolume(real_t speed) { + // scaling for integration step + auto* param = Simulation::GetActive()->GetParam(); + real_t delta = speed * param->simulation_time_step; + volume_ += delta; + if (volume_ < real_t(5.2359877E-7)) { + volume_ = real_t(5.2359877E-7); + } + UpdateDiameter(); + } + + void UpdateDiameter() { + // V = (4/3)*pi*r^3 = (pi/6)*diameter^3 + real_t diameter = std::cbrt(volume_ * 6 / Math::kPi); + if (diameter > diameter_) { + Base::SetPropagateStaticness(); + } + diameter_ = diameter; + } + + void UpdateVolume() { + // V = (4/3)*pi*r^3 = (pi/6)*diameter^3 + volume_ = Math::kPi / real_t(6) * std::pow(diameter_, 3); + } + + void UpdatePosition(const Real3& delta) { + position_ += delta; + SetPropagateStaticness(); + } + + // empty unused functions + void ApplyDisplacement(const Real3& displacement) override{}; + + Real3 CalculateDisplacement(const InteractionForce* force, + real_t squared_radius, real_t dt) override { + return Real3{0, 0, 0}; + } + + private: + Real3 position_ = {0, 0, 0}; + real_t diameter_ = 0; + real_t volume_ = 0; + real_t density_ = 0; + Real3 velocity_ = {0, 0, 0}; +}; + +// Star subclass +class Star : public CelestialObject { + BDM_AGENT_HEADER(Star, CelestialObject, 1); + + public: + Star() : CelestialObject() {} + + explicit Star(real_t diameter) : CelestialObject(diameter) {} + + explicit Star(const Real3 position) : CelestialObject(position) {} + + ~Star() override = default; +}; + +// Planet subclass +class Planet : public CelestialObject { + BDM_AGENT_HEADER(Planet, CelestialObject, 1); + + public: + Planet() : CelestialObject() {} + + explicit Planet(real_t diameter) : CelestialObject(diameter) {} + + explicit Planet(const Real3 position) : CelestialObject(position) {} + + explicit Planet(const Star* main_star, real_t diameter, real_t distance, + real_t initial_speed) { + Real3 pos = main_star->GetPosition(); + pos[0] += distance; + Real3 initial_velocity = {0, initial_speed, 0}; + this->SetVelocity(main_star->GetVelocity() + initial_velocity); + this->SetPosition(pos); + this->SetDiameter(diameter); + } + + ~Planet() override = default; +}; + +// Satellite subclass +class Satellite : public CelestialObject { + BDM_AGENT_HEADER(Satellite, CelestialObject, 1); + + public: + Satellite() : CelestialObject() {} + + explicit Satellite(real_t diameter) : CelestialObject(diameter) {} + + explicit Satellite(const Real3 position) : CelestialObject(position) {} + + explicit Satellite(const Planet* main_planet, real_t diameter, + real_t distance, real_t initial_speed) { + Real3 pos = main_planet->GetPosition(); + pos[0] += distance; + Real3 initial_velocity = {0, initial_speed, 0}; + this->SetVelocity(main_planet->GetVelocity() + initial_velocity); + + this->SetPosition(pos); + this->SetDiameter(diameter); + } + + ~Satellite() override = default; +}; + +} // namespace astrophysics +} // namespace bdm + +#endif // CELESTIALOBJECT_H_ diff --git a/demo/newtons_law_test/src/Gravity.h b/demo/newtons_law_test/src/Gravity.h new file mode 100644 index 000000000..70117b706 --- /dev/null +++ b/demo/newtons_law_test/src/Gravity.h @@ -0,0 +1,116 @@ +// ----------------------------------------------------------------------------- +// +// Copyright (C) 2021 CERN & University of Surrey for the benefit of the +// BioDynaMo collaboration. All Rights Reserved. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// +// See the LICENSE file distributed with this work for details. +// See the NOTICE file distributed with this work for additional information +// regarding copyright ownership. +// +// ----------------------------------------------------------------------------- +#ifndef GRAVITY_H_ +#define GRAVITY_H_ + +#include "CelestialObjects.h" +#include "biodynamo.h" + +namespace bdm { +namespace astrophysics { + +// Gravitational constant in m^3/kg/s^2 +inline real_t gG = 6.674e-11; + +// set units of G constant +inline void SetUnits(char units = 'O', int times = 1) { + gG = 6.674e-11; + switch (units) { + case 'k': + gG *= pow(times * 1e-3, 3); + break; + case 'M': + gG *= pow(times * 1e-6, 3); + break; + case 'G': + gG *= pow(times * 1e-9, 3); + break; + default: + gG *= pow(times, 3); + break; + } +} + +// a behavior that applies Newtons Law for gravity and calculates the +// displacement for any celestial object +class Gravity : public Behavior { + Real3 acceleration_ = {0, 0, 0}; + real_t dt_ = Simulation::GetActive()->GetParam()->simulation_time_step; + + public: + Behavior* New() const override { return new Gravity(); } + + Behavior* NewCopy() const override { return new Gravity(*this); } + + void Run(Agent* agent) override { + CelestialObject* celestial = dynamic_cast(agent); + + Real3 celestial_pos = celestial->GetPosition(); + + auto* functor = new GravityFunctor(celestial, gG); + // for each agent calculate acceleration + Simulation::GetActive()->GetResourceManager()->ForEachAgentParallel( + *functor); + + dt_ = Simulation::GetActive()->GetParam()->simulation_time_step; + Real3 velocity = celestial->GetVelocity(); + Real3 acceleration = functor->acceleration; + Real3 position; + + // explicit Euler's method + velocity += acceleration * dt_; + position = celestial_pos + velocity * dt_; + + delete functor; + + celestial->SetVelocity(velocity); + celestial->SetPosition(position); + } + + class GravityFunctor : public Functor { + public: + CelestialObject* celestial; + real_t g; + Real3 acceleration; + + GravityFunctor(CelestialObject* celestial, real_t g) + : celestial(celestial), g(g), acceleration({0, 0, 0}) {} + + // override the function that overloads () operator to + // calculate acceleration based on given agent + void operator()(Agent* other_agent, AgentHandle) override { + CelestialObject* other_celestial = + dynamic_cast(other_agent); + + if (!other_celestial || other_celestial == celestial) { + return; + } + + Real3 dir = other_celestial->GetPosition() - celestial->GetPosition(); + real_t magnitude = dir.Norm(); + + dir.Normalize(); + real_t m1 = celestial->GetMass(); + real_t m2 = other_celestial->GetMass(); + real_t f = g * m1 * m2 / (magnitude * magnitude); + + acceleration += (f / m1) * dir; + }; + }; +}; +} // namespace astrophysics + +} // namespace bdm + +#endif // GRAVITY_H_ diff --git a/demo/newtons_law_test/src/newtons_law_test.cc b/demo/newtons_law_test/src/newtons_law_test.cc new file mode 100644 index 000000000..f64b690d5 --- /dev/null +++ b/demo/newtons_law_test/src/newtons_law_test.cc @@ -0,0 +1,18 @@ +// ----------------------------------------------------------------------------- +// +// Copyright (C) 2021 CERN & University of Surrey for the benefit of the +// BioDynaMo collaboration. All Rights Reserved. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// +// See the LICENSE file distributed with this work for details. +// See the NOTICE file distributed with this work for additional information +// regarding copyright ownership. +// +// ----------------------------------------------------------------------------- +#include "newtons_law_test.h" + +int main(int argc, const char** argv) { + return bdm::astrophysics::Simulate(argc, argv); +} diff --git a/demo/newtons_law_test/src/newtons_law_test.h b/demo/newtons_law_test/src/newtons_law_test.h new file mode 100644 index 000000000..8d74e2573 --- /dev/null +++ b/demo/newtons_law_test/src/newtons_law_test.h @@ -0,0 +1,68 @@ +// ----------------------------------------------------------------------------- +// +// Copyright (C) 2021 CERN & University of Surrey for the benefit of the +// BioDynaMo collaboration. All Rights Reserved. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// +// See the LICENSE file distributed with this work for details. +// See the NOTICE file distributed with this work for additional information +// regarding copyright ownership. +// +// ----------------------------------------------------------------------------- +#ifndef NEWTONSLAWTEST_H_ +#define NEWTONSLAWTEST_H_ + +#include "CelestialObjects.h" +#include "Gravity.h" +#include "biodynamo.h" + +namespace bdm { +namespace astrophysics { + +inline int Simulate(int argc, const char** argv) { + Simulation simulation(argc, argv); + auto* rm = simulation.GetResourceManager(); + + // create 4 celestial objects. Apply mass initial position and velocity + auto* object1 = new CelestialObject(10); + object1->SetMass(10e+16); + object1->SetPosition({50, 0, 0}); + object1->SetVelocity({100, 100, 100}); + + auto* object2 = new CelestialObject(5); + object2->SetMass(10e+16); + object2->SetPosition({-50, 0, 0}); + object2->SetVelocity({100, -100, -100}); + + auto* object3 = new CelestialObject(5); + object3->SetMass(5e+16); + object3->SetPosition({-50, -50, 0}); + object3->SetVelocity({200, -50, -20}); + + auto* object4 = new CelestialObject(50); + object4->SetMass(10e+17); + object4->SetPosition({100, 100, 0}); + object4->SetVelocity({-20, 50, 20}); + + // Add gravity behaviour to all agents + object1->AddBehavior(new Gravity); + object2->AddBehavior(new Gravity); + object3->AddBehavior(new Gravity); + object4->AddBehavior(new Gravity); + + // Add agents to simulation + rm->AddAgent(object1); + rm->AddAgent(object2); + rm->AddAgent(object3); + rm->AddAgent(object4); + + simulation.GetScheduler()->Simulate(200); + + return 0; +} + +} // namespace astrophysics +} // namespace bdm +#endif // NEWTONSLAWTEST_H_ From 7341a137a453d9477745340df6ffb69dbf27b236 Mon Sep 17 00:00:00 2001 From: christosneg Date: Thu, 2 Jan 2025 15:39:57 +0200 Subject: [PATCH 2/8] Add second demo showcasing Mars and its satellites Phobos and Deimos This commit introduces a second demo featuring Mars and its two satellites, Phobos and Deimos, modeled as agents influenced by the Gravity behavior. --- demo/mars/CMakeLists.txt | 35 ++++++ demo/mars/Mars.png | Bin 0 -> 4972 bytes demo/mars/README.md | 57 +++++++++ demo/mars/bdm.toml | 11 ++ demo/mars/src/CelestialObjects.h | 199 +++++++++++++++++++++++++++++++ demo/mars/src/Gravity.h | 116 ++++++++++++++++++ demo/mars/src/mars.cc | 48 ++++++++ demo/mars/src/mars.h | 66 ++++++++++ 8 files changed, 532 insertions(+) create mode 100644 demo/mars/CMakeLists.txt create mode 100644 demo/mars/Mars.png create mode 100644 demo/mars/README.md create mode 100644 demo/mars/bdm.toml create mode 100644 demo/mars/src/CelestialObjects.h create mode 100644 demo/mars/src/Gravity.h create mode 100644 demo/mars/src/mars.cc create mode 100644 demo/mars/src/mars.h diff --git a/demo/mars/CMakeLists.txt b/demo/mars/CMakeLists.txt new file mode 100644 index 000000000..cf002028a --- /dev/null +++ b/demo/mars/CMakeLists.txt @@ -0,0 +1,35 @@ +# ----------------------------------------------------------------------------- +# +# Copyright (C) 2021 CERN & University of Surrey for the benefit of the +# BioDynaMo collaboration. All Rights Reserved. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# +# See the LICENSE file distributed with this work for details. +# See the NOTICE file distributed with this work for additional information +# regarding copyright ownership. +# +# ----------------------------------------------------------------------------- +cmake_minimum_required(VERSION 3.19.3) +project(mars) + +# BioDynaMo curretly uses the C++17 standard. +set(CMAKE_CXX_STANDARD 17) + +# Use BioDynaMo in this project. +find_package(BioDynaMo REQUIRED) + +# See UseBioDynaMo.cmake in your BioDynaMo build folder for details. +# Note that BioDynaMo provides gtest header/libraries in its include/lib dir. +include(${BDM_USE_FILE}) + +# Consider all files in src/ for BioDynaMo simulation. +include_directories("src") +file(GLOB_RECURSE PROJECT_HEADERS src/*.h) +file(GLOB_RECURSE PROJECT_SOURCES src/*.cc) + +bdm_add_executable(${CMAKE_PROJECT_NAME} + HEADERS ${PROJECT_HEADERS} + SOURCES ${PROJECT_SOURCES} + LIBRARIES ${BDM_REQUIRED_LIBRARIES}) diff --git a/demo/mars/Mars.png b/demo/mars/Mars.png new file mode 100644 index 0000000000000000000000000000000000000000..c857ca9f730bdd336cd8f3f3975cfd554e142c81 GIT binary patch literal 4972 zcmeHL`8$;D-yTm{B9wTn8AbL`Sq35N*o7HdY+qLtr36F^0&6#Z& z9RKC>2BG88I%Pv+_s)+frgJyq+0N>V6iXZQUr6JTiQ+MgH>*=$;9-_A4X|6~jL%jzd41Q^? zC@&yir(Czn%HU-Nb)MaWgF&Dl9vrhE(2Fx%pr}*-HvaEPsPBuQK3pM^LvMkPmggtG z3Dnd4{G#*AT2^B~eoLA`-WK?+PntWO;~cF4Kigts%+mvv`rhKgUgBvxOMz=&yAo4* zq#pZ$&~Ra4w}`EC;IHLj@00~3gxsq~l+#d_aJIxNtKIw|J2``qCRPztEF7$mA~*|r zJvKRclK~kQIMeQ-BIMCOVi2QSp-c%N$UKqmu97vdsC?gXnC!9h9rI@d;1=$ccBC+Y zS{mE8-^4X_sM-`}o7y-Z0W0}>#tZqpZm=ZE6u{mZ`ZHQ#Cd6_O?}r+XOL5bKfaFz# zI0KaivhF7w)ZSzZk_MFjQLhvz(PK-r&NW9_cwiMss93+Y$7i`fW{BGWaJqlnh7cxW zm!3Z&@N6kmn8eybkwW(?WW&&T(X-1BWM;PM{-op>$&U`hkGx)?O;Og3F)6| zfD!<(NwR)q?Tx`y%qrcb+?y9Tj!T^UKH>7L|52sxbD%bwhYaT&4TZzu)$(gu@W$K4 z^q8Dosz;6X0}=5LOfJ=VYE)&RRNj(Nbf(wy>@9+ia{uD&KTl09N8}Zonsvr^U(w-U z)s(Zx90un4Xvn~L%eTVIBsRB5a+<$Q(+#e$SX5#c@+(Fin%_}_=6lTkC-hPwu6Uc@ z(kDNu38Q4orkie?i+)6wkb=Ws=s3~(QrcW?OrVq=XG|J(qZ5N#<-P7Y!WdWU#e;c(}Es=sN;`|sAK&j0UC{yOyp+yH-{szbTq4PBS|-$o+^Xz&8Y_V%4WyQo zm7Hb-;M%fLdq>@=#?TTwcH4#J=%grnJxJ8V^1JhH*PpLAIUCFP<23|G>WBpU zt6;B|MX&>Z@);Z^62Ye;XJ)Y#y0ckyw69P9u7kl2vu;41oIocZ&V^{d4ADW?e~qh- zy$y_oh;x!1P~Je8Wv+9%zV*{t&;8-T)?D6luKU=?-x(3~t{IH34oMxdhdhI3t_A*b zifMMnRrS^pOYc3GHse8+fH8#gmC1bT{ z?^9gR=GY{FN4TW;fsB?^dugSqZP1unj@XbOkx74S#v-Or$ANyfQc`4nqQ&V}Zm*by zvS`^IFrxKnS|(yK0}nW(e~+-^eO_i5(=y=TlfN71rh9H-_Hbav>In-W5?H*nnrFGN zF4}`J8oITWfD5f&x1qaQp%UDU4;zv977MBliCl-@>ha=<_+MGj5 zsfWdn z;?m)FS{ZDOa4*r2_4)W%n{t;tQN8IhuE@Qi$5>N7@P23Mwx>i#VOYE(HAZ_(vnyLm zk9-BA9SpAdFgJ4>+nr`2oEChzFd>TTSj}b)xljjZmiMO@UUtNL)J(Na;9&76YfIg9 z`^*K$2ik*X?MtN@B{j(UJ62&KwgMF6i>j5StqI#L?oW8!GDCu{Z)oEDOX^w0&T_d_ zuTdXXCt4d;!LSB!Hrbf^(%Iiza0KcZ@dohwzB4fN8|$DmUqif&35)(6zSw4?PEu6 z)(bhnB@UY1=VN%B>w$Dv)T`>MvQyAs&VZ?RaK734!k_RQU;3EuAe9={6n3NkOkSGt z1BI}^@F-%V8wYAo{Qqq`RU&rnAMC7nY)*HOS7qzp3V(r zI1(J!^cwKiW>?<{$BBg8V-2k>;;?>lLW)r>%*H;utj+AKwV%r>TjYh=;CvJ6@QJNCgX9v=y;pojKZ7UNBgl#Kfi%65@~N2AugclGuWq6@50o8VAr-OTN~9prpFlu z3?1AzEN?q4`}yVVYm?7feax}9>b(Iyl-xIgETei^jY|Uv-AQ#W#r9tT)4}e;m<|k6 z_ao`tpDl)bjFW!&dQo;()PKVO0?IO`jJnMb!+wk-djhOGIPk(HLXBw&&jbuE~S&m;nF#7(qM?w*vnY9Jj5UOSKl|J~NoTzgg+sABPzT*t;-#4H-Ygabj6>bft zlpNBK+LeM-OyoyzMd=o`K0deL#JkH zK(43#_Df1^$en3`+i1~V@ZIq-I)6I}$Th2P7Edf=AE=gXK9dqWb|DDYyI)y0rb!P% z3qBN#-YOZohg@ zbIsK)C`PJsHGA$s5$n%ru<^s%Atec5#K4OuE0$Ilp6>H}DRz&Ch|PkpYpLhngJRyV zhmUAYi>9;z$~k-kRZ?G|vN{akH_+bK?*yJMI0oC!v+{96ImXD3yF9A`qAN!AGFHrw z=mj435om~!75zs=U7^ma`BF z+X5teW`6g?n1jQ2bPGhevm~4nwgWnvIj$RVj<_GzS93z=rk{wzFu}sSiC%;C{_B- zc7(uyj0Abr1zkz<0)hkDbo3)VRp8{3SfavKcXCwNsKXKMRb>K;U92Ansx0ardhm*S z^$X4T_PNyE9oW&SQ;UBc1%?Dj_;oM~T2a~F~V-)KBJ zp(Mb4wSAw4VA#M4QswNJeA+5jpUBLyB%Y+t1+AExf1@O@^6EAFWYHpRPK^PW$|&<) z3>y+UlUMJBQV#!)qhaq0{|VpoCHw4Cr#)woCq{eR?x9>ziL0Njv`fJ=u6t^BsT_MJ z>btm`sXKf1^weOk3Su2OxamSQ-W0fds(f!Cr?uWo3NjWxzCK)~H^ia1&ciut3DZgWMr^mTA@7W?XT$7aNa&-GLhs1$X=mOtI6VgQBDg~~PCq&E zM<$syheR1VPU}v<#yg*z6vZnSFO!gMMBZ;R)+Ec!xHCaMSI6;Mds>Sse$Wddoofq) zqr#!&ygKj9@@-pdBE4(sXR9}{5Wv)i$+F2*hxRgEC%r#m+U$w|q z^z>C4LcahKt@kq=!WtcL897|{-&Z{9Tah~Wa!~{XI`@43FOWGKE2#5{AL!r4za;$s dC}H7 + Simulation + + +## Contents + +- **Celestial Object Agents** +- **Gravity Behavior** + +## Assumptions + +To simplify the simulation, the following assumptions were made: + - All celestial objects (Mars, Phobos, and Deimos) are modeled as **spheres**. + - Celestial objects are treated as **point masses** for computational efficiency. + - Agents can't physically interact with each other. + +## Details + +- **Mars**: Radius ~3,389.5 km, Mass ~6.417 × 10^23 kg. +- **Phobos**: Radius ~11.267 km, Mass ~1.0659 × 10^16 kg, Average orbital distance ~9,376 km. +- **Deimos**: Radius ~6.2 km, Mass ~1.4762 × 10^15 kg, Average orbital distance ~23,463 km. + +## Newton's Law of Universal Gravitation + +Newton’s law of universal gravitation governs the movement of all celestial objects in the simulation, describing the force of attraction between two masses. + +The force of attraction between two objects is proportional to both their masses and inversely proportional to the square of the distance between them. + +$$F = \frac{G * (m_1 * m_2)}{r^2}$$ + +where: +- $F$ is the gravitational force between the two masses, +- $G$ is the gravitational constant (6.674 × 10^-11 N·m²/kg²), +- $m_1$ and $m_2$ are the masses of the two objects, and +- $r$ is the distance between the centers of the two masses. + +By knowing the force of attraction and therefore the acceleration of a celestial object, the displacement can be predicted by solving the below differential equation: + +$$\frac{d^2x}{dt^2} = a$$ + +--- + +## Run the demos + +In order to run the demos, BioDynamo has to be correctly installed and sourced. + +```bash +cd [path_to_biodynamo]demos/NewtonsLawTest +bdm run +bdm view +``` +In order to clearly observe the satellites in relation to the planet, feel free to scale up the diameter or volume of the satellites in the simulation. + diff --git a/demo/mars/bdm.toml b/demo/mars/bdm.toml new file mode 100644 index 000000000..d6e7d9f5f --- /dev/null +++ b/demo/mars/bdm.toml @@ -0,0 +1,11 @@ +[visualization] +export = true +interval = 1 + +[simulation] +time_step = 360 + +[[visualize_agent]] +name = "Planet" +[[visualize_agent]] +name = "Satellite" \ No newline at end of file diff --git a/demo/mars/src/CelestialObjects.h b/demo/mars/src/CelestialObjects.h new file mode 100644 index 000000000..10cddf869 --- /dev/null +++ b/demo/mars/src/CelestialObjects.h @@ -0,0 +1,199 @@ +// ----------------------------------------------------------------------------- +// +// Copyright (C) 2021 CERN & University of Surrey for the benefit of the +// BioDynaMo collaboration. All Rights Reserved. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// +// See the LICENSE file distributed with this work for details. +// See the NOTICE file distributed with this work for additional information +// regarding copyright ownership. +// +// ----------------------------------------------------------------------------- +#ifndef CELESTIALOBJECTS_H_ +#define CELESTIALOBJECTS_H_ + +#include "biodynamo.h" + +namespace bdm { +namespace astrophysics { + +// creating the custom agent that +// defines celestial objects +class CelestialObject : public Agent { + BDM_AGENT_HEADER(CelestialObject, Agent, 1); + + public: + // constructors + CelestialObject() : diameter_(1.0), density_(1.0) { UpdateVolume(); } + + explicit CelestialObject(real_t diameter) + : diameter_(diameter), density_(1.0) { + UpdateVolume(); + } + + explicit CelestialObject(const Real3& position) + : position_(position), diameter_(1.0), density_(1.0) { + UpdateVolume(); + } + + // default destructor + ~CelestialObject() override = default; + + Real3 GetVelocity() const { return velocity_; } + + real_t GetDiameter() const override { return diameter_; } + + real_t GetMass() const { return density_ * volume_; } + + real_t GetDensity() const { return density_; } + + const Real3& GetPosition() const override { return position_; } + + real_t GetVolume() const { return volume_; } + + Shape GetShape() const override { return Shape::kSphere; } + + void SetVelocity(const Real3& vel) { velocity_ = vel; } + + void SetDiameter(real_t diameter) override { + if (diameter > diameter_) { + SetPropagateStaticness(); + } + diameter_ = diameter; + UpdateVolume(); + } + + void SetVolume(real_t volume) { + volume_ = volume; + UpdateDiameter(); + } + + void SetMass(real_t mass) { SetDensity(mass / volume_); } + + void SetDensity(real_t density) { + if (density > density_) { + SetPropagateStaticness(); + } + density_ = density; + } + + void SetPosition(const Real3& position) override { + position_ = position; + SetPropagateStaticness(); + } + + void ChangeVolume(real_t speed) { + // scaling for integration step + auto* param = Simulation::GetActive()->GetParam(); + real_t delta = speed * param->simulation_time_step; + volume_ += delta; + if (volume_ < real_t(5.2359877E-7)) { + volume_ = real_t(5.2359877E-7); + } + UpdateDiameter(); + } + + void UpdateDiameter() { + // V = (4/3)*pi*r^3 = (pi/6)*diameter^3 + real_t diameter = std::cbrt(volume_ * 6 / Math::kPi); + if (diameter > diameter_) { + Base::SetPropagateStaticness(); + } + diameter_ = diameter; + } + + void UpdateVolume() { + // V = (4/3)*pi*r^3 = (pi/6)*diameter^3 + volume_ = Math::kPi / real_t(6) * std::pow(diameter_, 3); + } + + void UpdatePosition(const Real3& delta) { + position_ += delta; + SetPropagateStaticness(); + } + + // empty unused functions + void ApplyDisplacement(const Real3& displacement) override{}; + + Real3 CalculateDisplacement(const InteractionForce* force, + real_t squared_radius, real_t dt) override { + return Real3{0, 0, 0}; + } + + private: + Real3 position_ = {0, 0, 0}; + real_t diameter_ = 0; + real_t volume_ = 0; + real_t density_ = 0; + Real3 velocity_ = {0, 0, 0}; +}; + +// Star subclass +class Star : public CelestialObject { + BDM_AGENT_HEADER(Star, CelestialObject, 1); + + public: + Star() : CelestialObject() {} + + explicit Star(real_t diameter) : CelestialObject(diameter) {} + + explicit Star(const Real3 position) : CelestialObject(position) {} + + ~Star() override = default; +}; + +// Planet subclass +class Planet : public CelestialObject { + BDM_AGENT_HEADER(Planet, CelestialObject, 1); + + public: + Planet() : CelestialObject() {} + + explicit Planet(real_t diameter) : CelestialObject(diameter) {} + + explicit Planet(const Real3 position) : CelestialObject(position) {} + + explicit Planet(const Star* main_star, real_t diameter, real_t distance, + real_t initial_speed) { + Real3 pos = main_star->GetPosition(); + pos[0] += distance; + Real3 initial_velocity = {0, initial_speed, 0}; + this->SetVelocity(main_star->GetVelocity() + initial_velocity); + this->SetPosition(pos); + this->SetDiameter(diameter); + } + + ~Planet() override = default; +}; + +// Satellite subclass +class Satellite : public CelestialObject { + BDM_AGENT_HEADER(Satellite, CelestialObject, 1); + + public: + Satellite() : CelestialObject() {} + + explicit Satellite(real_t diameter) : CelestialObject(diameter) {} + + explicit Satellite(const Real3 position) : CelestialObject(position) {} + + explicit Satellite(const Planet* main_planet, real_t diameter, + real_t distance, real_t initial_speed) { + Real3 pos = main_planet->GetPosition(); + pos[0] += distance; + Real3 initial_velocity = {0, initial_speed, 0}; + this->SetVelocity(main_planet->GetVelocity() + initial_velocity); + + this->SetPosition(pos); + this->SetDiameter(diameter); + } + + ~Satellite() override = default; +}; + +} // namespace astrophysics +} // namespace bdm + +#endif // CELESTIALOBJECT_H_ diff --git a/demo/mars/src/Gravity.h b/demo/mars/src/Gravity.h new file mode 100644 index 000000000..70117b706 --- /dev/null +++ b/demo/mars/src/Gravity.h @@ -0,0 +1,116 @@ +// ----------------------------------------------------------------------------- +// +// Copyright (C) 2021 CERN & University of Surrey for the benefit of the +// BioDynaMo collaboration. All Rights Reserved. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// +// See the LICENSE file distributed with this work for details. +// See the NOTICE file distributed with this work for additional information +// regarding copyright ownership. +// +// ----------------------------------------------------------------------------- +#ifndef GRAVITY_H_ +#define GRAVITY_H_ + +#include "CelestialObjects.h" +#include "biodynamo.h" + +namespace bdm { +namespace astrophysics { + +// Gravitational constant in m^3/kg/s^2 +inline real_t gG = 6.674e-11; + +// set units of G constant +inline void SetUnits(char units = 'O', int times = 1) { + gG = 6.674e-11; + switch (units) { + case 'k': + gG *= pow(times * 1e-3, 3); + break; + case 'M': + gG *= pow(times * 1e-6, 3); + break; + case 'G': + gG *= pow(times * 1e-9, 3); + break; + default: + gG *= pow(times, 3); + break; + } +} + +// a behavior that applies Newtons Law for gravity and calculates the +// displacement for any celestial object +class Gravity : public Behavior { + Real3 acceleration_ = {0, 0, 0}; + real_t dt_ = Simulation::GetActive()->GetParam()->simulation_time_step; + + public: + Behavior* New() const override { return new Gravity(); } + + Behavior* NewCopy() const override { return new Gravity(*this); } + + void Run(Agent* agent) override { + CelestialObject* celestial = dynamic_cast(agent); + + Real3 celestial_pos = celestial->GetPosition(); + + auto* functor = new GravityFunctor(celestial, gG); + // for each agent calculate acceleration + Simulation::GetActive()->GetResourceManager()->ForEachAgentParallel( + *functor); + + dt_ = Simulation::GetActive()->GetParam()->simulation_time_step; + Real3 velocity = celestial->GetVelocity(); + Real3 acceleration = functor->acceleration; + Real3 position; + + // explicit Euler's method + velocity += acceleration * dt_; + position = celestial_pos + velocity * dt_; + + delete functor; + + celestial->SetVelocity(velocity); + celestial->SetPosition(position); + } + + class GravityFunctor : public Functor { + public: + CelestialObject* celestial; + real_t g; + Real3 acceleration; + + GravityFunctor(CelestialObject* celestial, real_t g) + : celestial(celestial), g(g), acceleration({0, 0, 0}) {} + + // override the function that overloads () operator to + // calculate acceleration based on given agent + void operator()(Agent* other_agent, AgentHandle) override { + CelestialObject* other_celestial = + dynamic_cast(other_agent); + + if (!other_celestial || other_celestial == celestial) { + return; + } + + Real3 dir = other_celestial->GetPosition() - celestial->GetPosition(); + real_t magnitude = dir.Norm(); + + dir.Normalize(); + real_t m1 = celestial->GetMass(); + real_t m2 = other_celestial->GetMass(); + real_t f = g * m1 * m2 / (magnitude * magnitude); + + acceleration += (f / m1) * dir; + }; + }; +}; +} // namespace astrophysics + +} // namespace bdm + +#endif // GRAVITY_H_ diff --git a/demo/mars/src/mars.cc b/demo/mars/src/mars.cc new file mode 100644 index 000000000..bb0b79351 --- /dev/null +++ b/demo/mars/src/mars.cc @@ -0,0 +1,48 @@ +// ----------------------------------------------------------------------------- +// +// Copyright (C) 2021 CERN & University of Surrey for the benefit of the +// BioDynaMo collaboration. All Rights Reserved. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// +// See the LICENSE file distributed with this work for details. +// See the NOTICE file distributed with this work for additional information +// regarding copyright ownership. +// +// ----------------------------------------------------------------------------- +#include "mars.h" + +int main(int argc, const char** argv) { return bdm::astrophysics::Simulate(argc, argv); } + +namespace bdm { + namespace astrophysics { + + Real3 RotateVector(const Real3& Vec, real_t th, char Axis) { + th *= M_PI / 180; + real_t costh = std::cos(th), sinth = std::sin(th); + + Real3 RotVec; + switch (Axis) { + case 'x': + RotVec[0] = Vec[0]; + RotVec[1] = Vec[1] * costh - Vec[2] * sinth; + RotVec[2] = Vec[1] * sinth + Vec[2] * costh; + break; + case 'y': + RotVec[0] = Vec[0] * costh + Vec[2] * sinth; + RotVec[1] = Vec[1]; + RotVec[2] = -Vec[0] * sinth + Vec[2] * costh; + break; + case 'z': + RotVec[0] = Vec[0] * costh - Vec[1] * sinth; + RotVec[1] = Vec[0] * sinth + Vec[1] * costh; + RotVec[2] = Vec[2]; + break; + default: + throw std::invalid_argument("Invalid rotation axis. Must be 'x', 'y', or 'z'."); + } + return RotVec; + } + } // namespace astrophysics +} // namespace bdm diff --git a/demo/mars/src/mars.h b/demo/mars/src/mars.h new file mode 100644 index 000000000..17ae3c834 --- /dev/null +++ b/demo/mars/src/mars.h @@ -0,0 +1,66 @@ +// ----------------------------------------------------------------------------- +// +// Copyright (C) 2021 CERN & University of Surrey for the benefit of the +// BioDynaMo collaboration. All Rights Reserved. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// +// See the LICENSE file distributed with this work for details. +// See the NOTICE file distributed with this work for additional information +// regarding copyright ownership. +// +// ----------------------------------------------------------------------------- +#ifndef MARS_H_ +#define MARS_H_ + +#include "biodynamo.h" +#include "CelestialObjects.h" +#include "Gravity.h" + +namespace bdm { + namespace astrophysics { + + //Rotate a vector by a specified angle along a given axis + Real3 RotateVector(const Real3& vec, real_t th, char axis); + + inline int Simulate(int argc, const char** argv) { + Simulation simulation(argc, argv); + auto* rm = simulation.GetResourceManager(); + + SetUnits('k'); // Simulation units in kilometers. + + //Create Mars with distance from the Sun and initial orbital velocity. + auto* mars = new Planet(6779); + mars->SetPosition({0, 0, 0}); + mars->SetMass(0.64171e+24); + mars->AddBehavior(new Gravity); + rm->AddAgent(mars); + + //Create Phobos + auto* phobos = new Satellite(22.4); // Phobos diameter in kilometers + phobos->SetMass(1.0659e+16); // Phobos mass in kg + phobos->SetPosition({mars->GetPosition()[0] + 9376, mars->GetPosition()[1], mars->GetPosition()[2]}); // Position relative to Mars + Real3 v_phobos = mars->GetVelocity() + RotateVector({0, 2.138, 0}, 0, 'x'); // Phobos velocity relative to Mars + phobos->SetVelocity(v_phobos); + phobos->AddBehavior(new Gravity); + rm->AddAgent(phobos); + + //Create Deimos + auto* deimos = new Satellite(12.4); + deimos->SetMass(1.48e+15); + deimos->SetPosition({mars->GetPosition()[0] + 23459, mars->GetPosition()[1], mars->GetPosition()[2]}); // Position relative to Mars + Real3 v_deimos = mars->GetVelocity() + RotateVector({0, 1.351, 0}, 0, 'x'); // Deimos velocity relative to Mars + deimos->SetVelocity(v_deimos); + deimos->AddBehavior(new Gravity); + rm->AddAgent(deimos); + + simulation.GetScheduler()->Simulate(200); + + return 0; + } + + } // namespace astrophysics +} // namespace bdm + +#endif // MARS_H_ From d07e951f9070e9433ac7e12be0584e6434dd563e Mon Sep 17 00:00:00 2001 From: christosneg Date: Thu, 2 Jan 2025 15:50:01 +0200 Subject: [PATCH 3/8] Add third demo simulating the solar system This demo introduces the ability to read data from a .txt file to create a simulation of the solar system, including the Sun, its eight planets, and Pluto as agents influenced by gravity. --- demo/mars/src/mars.cc | 61 +++---- demo/mars/src/mars.h | 78 +++++---- demo/solar_system/CMakeLists.txt | 35 ++++ demo/solar_system/README.md | 65 ++++++++ demo/solar_system/SolarSystemSim.png | Bin 0 -> 10146 bytes demo/solar_system/bdm.toml | 12 ++ demo/solar_system/planets.txt | 10 ++ demo/solar_system/src/CelestialObjects.h | 199 +++++++++++++++++++++++ demo/solar_system/src/DataReader.h | 70 ++++++++ demo/solar_system/src/Gravity.h | 116 +++++++++++++ demo/solar_system/src/solar_system.cc | 51 ++++++ demo/solar_system/src/solar_system.h | 63 +++++++ 12 files changed, 695 insertions(+), 65 deletions(-) create mode 100644 demo/solar_system/CMakeLists.txt create mode 100644 demo/solar_system/README.md create mode 100644 demo/solar_system/SolarSystemSim.png create mode 100644 demo/solar_system/bdm.toml create mode 100644 demo/solar_system/planets.txt create mode 100644 demo/solar_system/src/CelestialObjects.h create mode 100644 demo/solar_system/src/DataReader.h create mode 100644 demo/solar_system/src/Gravity.h create mode 100644 demo/solar_system/src/solar_system.cc create mode 100644 demo/solar_system/src/solar_system.h diff --git a/demo/mars/src/mars.cc b/demo/mars/src/mars.cc index bb0b79351..d96e8d17b 100644 --- a/demo/mars/src/mars.cc +++ b/demo/mars/src/mars.cc @@ -13,36 +13,39 @@ // ----------------------------------------------------------------------------- #include "mars.h" -int main(int argc, const char** argv) { return bdm::astrophysics::Simulate(argc, argv); } +int main(int argc, const char** argv) { + return bdm::astrophysics::Simulate(argc, argv); +} namespace bdm { - namespace astrophysics { +namespace astrophysics { - Real3 RotateVector(const Real3& Vec, real_t th, char Axis) { - th *= M_PI / 180; - real_t costh = std::cos(th), sinth = std::sin(th); +Real3 RotateVector(const Real3& vec, real_t th, char axis) { + th *= M_PI / 180; + real_t costh = std::cos(th), sinth = std::sin(th); - Real3 RotVec; - switch (Axis) { - case 'x': - RotVec[0] = Vec[0]; - RotVec[1] = Vec[1] * costh - Vec[2] * sinth; - RotVec[2] = Vec[1] * sinth + Vec[2] * costh; - break; - case 'y': - RotVec[0] = Vec[0] * costh + Vec[2] * sinth; - RotVec[1] = Vec[1]; - RotVec[2] = -Vec[0] * sinth + Vec[2] * costh; - break; - case 'z': - RotVec[0] = Vec[0] * costh - Vec[1] * sinth; - RotVec[1] = Vec[0] * sinth + Vec[1] * costh; - RotVec[2] = Vec[2]; - break; - default: - throw std::invalid_argument("Invalid rotation axis. Must be 'x', 'y', or 'z'."); - } - return RotVec; - } - } // namespace astrophysics -} // namespace bdm + Real3 rot_vec; + switch (axis) { + case 'x': + rot_vec[0] = vec[0]; + rot_vec[1] = vec[1] * costh - vec[2] * sinth; + rot_vec[2] = vec[1] * sinth + vec[2] * costh; + break; + case 'y': + rot_vec[0] = vec[0] * costh + vec[2] * sinth; + rot_vec[1] = vec[1]; + rot_vec[2] = -vec[0] * sinth + vec[2] * costh; + break; + case 'z': + rot_vec[0] = vec[0] * costh - vec[1] * sinth; + rot_vec[1] = vec[0] * sinth + vec[1] * costh; + rot_vec[2] = vec[2]; + break; + default: + throw std::invalid_argument( + "Invalid rotation axis. Must be 'x', 'y', or 'z'."); + } + return rot_vec; +} +} // namespace astrophysics +} // namespace bdm diff --git a/demo/mars/src/mars.h b/demo/mars/src/mars.h index 17ae3c834..98887b435 100644 --- a/demo/mars/src/mars.h +++ b/demo/mars/src/mars.h @@ -14,53 +14,59 @@ #ifndef MARS_H_ #define MARS_H_ -#include "biodynamo.h" #include "CelestialObjects.h" #include "Gravity.h" +#include "biodynamo.h" namespace bdm { - namespace astrophysics { - - //Rotate a vector by a specified angle along a given axis - Real3 RotateVector(const Real3& vec, real_t th, char axis); +namespace astrophysics { + +// Rotate a vector by a specified angle along a given axis +Real3 RotateVector(const Real3& vec, real_t th, char axis); + +inline int Simulate(int argc, const char** argv) { + Simulation simulation(argc, argv); + auto* rm = simulation.GetResourceManager(); - inline int Simulate(int argc, const char** argv) { - Simulation simulation(argc, argv); - auto* rm = simulation.GetResourceManager(); - - SetUnits('k'); // Simulation units in kilometers. + SetUnits('k'); // Simulation units in kilometers. - //Create Mars with distance from the Sun and initial orbital velocity. - auto* mars = new Planet(6779); - mars->SetPosition({0, 0, 0}); - mars->SetMass(0.64171e+24); - mars->AddBehavior(new Gravity); - rm->AddAgent(mars); + // Create Mars with distance from the Sun and initial orbital velocity. + auto* mars = new Planet(6779); + mars->SetPosition({0, 0, 0}); + mars->SetMass(0.64171e+24); + mars->AddBehavior(new Gravity); + rm->AddAgent(mars); - //Create Phobos - auto* phobos = new Satellite(22.4); // Phobos diameter in kilometers - phobos->SetMass(1.0659e+16); // Phobos mass in kg - phobos->SetPosition({mars->GetPosition()[0] + 9376, mars->GetPosition()[1], mars->GetPosition()[2]}); // Position relative to Mars - Real3 v_phobos = mars->GetVelocity() + RotateVector({0, 2.138, 0}, 0, 'x'); // Phobos velocity relative to Mars - phobos->SetVelocity(v_phobos); - phobos->AddBehavior(new Gravity); - rm->AddAgent(phobos); + // Create Phobos + auto* phobos = new Satellite(22.4); // Phobos diameter in kilometers + phobos->SetMass(1.0659e+16); // Phobos mass in kg + phobos->SetPosition({mars->GetPosition()[0] + 9376, mars->GetPosition()[1], + mars->GetPosition()[2]}); // Position relative to Mars + Real3 v_phobos = + mars->GetVelocity() + + RotateVector({0, 2.138, 0}, 0, 'x'); // Phobos velocity relative to Mars + phobos->SetVelocity(v_phobos); + phobos->AddBehavior(new Gravity); + rm->AddAgent(phobos); - //Create Deimos - auto* deimos = new Satellite(12.4); - deimos->SetMass(1.48e+15); - deimos->SetPosition({mars->GetPosition()[0] + 23459, mars->GetPosition()[1], mars->GetPosition()[2]}); // Position relative to Mars - Real3 v_deimos = mars->GetVelocity() + RotateVector({0, 1.351, 0}, 0, 'x'); // Deimos velocity relative to Mars - deimos->SetVelocity(v_deimos); - deimos->AddBehavior(new Gravity); - rm->AddAgent(deimos); + // Create Deimos + auto* deimos = new Satellite(12.4); + deimos->SetMass(1.48e+15); + deimos->SetPosition({mars->GetPosition()[0] + 23459, mars->GetPosition()[1], + mars->GetPosition()[2]}); // Position relative to Mars + Real3 v_deimos = + mars->GetVelocity() + + RotateVector({0, 1.351, 0}, 0, 'x'); // Deimos velocity relative to Mars + deimos->SetVelocity(v_deimos); + deimos->AddBehavior(new Gravity); + rm->AddAgent(deimos); - simulation.GetScheduler()->Simulate(200); + simulation.GetScheduler()->Simulate(500); - return 0; - } + return 0; +} - } // namespace astrophysics +} // namespace astrophysics } // namespace bdm #endif // MARS_H_ diff --git a/demo/solar_system/CMakeLists.txt b/demo/solar_system/CMakeLists.txt new file mode 100644 index 000000000..8a692d976 --- /dev/null +++ b/demo/solar_system/CMakeLists.txt @@ -0,0 +1,35 @@ +# ----------------------------------------------------------------------------- +# +# Copyright (C) 2021 CERN & University of Surrey for the benefit of the +# BioDynaMo collaboration. All Rights Reserved. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# +# See the LICENSE file distributed with this work for details. +# See the NOTICE file distributed with this work for additional information +# regarding copyright ownership. +# +# ----------------------------------------------------------------------------- +cmake_minimum_required(VERSION 3.19.3) +project(solar_system) + +# BioDynaMo curretly uses the C++14 standard. +set(CMAKE_CXX_STANDARD 14) + +# Use BioDynaMo in this project. +find_package(BioDynaMo REQUIRED) + +# See UseBioDynaMo.cmake in your BioDynaMo build folder for details. +# Note that BioDynaMo provides gtest header/libraries in its include/lib dir. +include(${BDM_USE_FILE}) + +# Consider all files in src/ for BioDynaMo simulation. +include_directories("src") +file(GLOB_RECURSE PROJECT_HEADERS src/*.h) +file(GLOB_RECURSE PROJECT_SOURCES src/*.cc) + +bdm_add_executable(${CMAKE_PROJECT_NAME} + HEADERS ${PROJECT_HEADERS} + SOURCES ${PROJECT_SOURCES} + LIBRARIES ${BDM_REQUIRED_LIBRARIES}) diff --git a/demo/solar_system/README.md b/demo/solar_system/README.md new file mode 100644 index 000000000..a97578d50 --- /dev/null +++ b/demo/solar_system/README.md @@ -0,0 +1,65 @@ +# SolarSystem + +This simulation models the gravitational interactions within the Solar System, including the Sun, all eight planets, and Pluto, using realistic dimensions and masses. This demo consists of these celestial bodies as agents, with gravitational interactions modeled using the gravity behavior. + +
+ Simulation +
+ +## Contents + +- **Celestial Object Agents** +- **Gravity Behavior** + +## Assumptions + +To simplify the simulation, the following assumptions were made: + - All celestial objects (Sun, planets, and Pluto) are modeled as **spheres**. + - Celestial objects are treated as **point masses** for computational efficiency. + - Agents can't physically interact with each other. + +## Details + +The following celestial objects are included in the simulation with their approximate physical properties: + +- **Sun**: Radius ~696,340 km, Mass ~1.989 × 10^30 kg. +- **Mercury**: Radius ~2,439.7 km, Mass ~3.301 × 10^23 kg, Average orbital distance ~57.91 million km. +- **Venus**: Radius ~6,051.8 km, Mass ~4.867 × 10^24 kg, Average orbital distance ~108.2 million km. +- **Earth**: Radius ~6,371 km, Mass ~5.972 × 10^24 kg, Average orbital distance ~149.6 million km. +- **Mars**: Radius ~3,389.5 km, Mass ~6.417 × 10^23 kg, Average orbital distance ~227.9 million km. +- **Jupiter**: Radius ~69,911 km, Mass ~1.898 × 10^27 kg, Average orbital distance ~778.5 million km. +- **Saturn**: Radius ~58,232 km, Mass ~5.683 × 10^26 kg, Average orbital distance ~1.434 billion km. +- **Uranus**: Radius ~25,362 km, Mass ~8.681 × 10^25 kg, Average orbital distance ~2.871 billion km. +- **Neptune**: Radius ~24,622 km, Mass ~1.024 × 10^26 kg, Average orbital distance ~4.495 billion km. +- **Pluto**: Radius ~1,188.3 km, Mass ~1.309 × 10^22 kg, Average orbital distance ~5.906 billion km. + +## Newton's Law of Universal Gravitation + +Newton’s law of universal gravitation governs the movement of all celestial objects in the simulation, describing the force of attraction between two masses. + +The force of attraction between two objects is proportional to both their masses and inversely proportional to the square of the distance between them. + +$$F = \frac{G * (m_1 * m_2)}{r^2}$$ + +where: +- $F$ is the gravitational force between the two masses, +- $G$ is the gravitational constant (6.674 × 10^-11 N·m²/kg²), +- $m_1$ and $m_2$ are the masses of the two objects, and +- $r$ is the distance between the centers of the two masses. + +By knowing the force of attraction and therefore the acceleration of a celestial object, the displacement can be predicted by solving the below differential equation: + +$$\frac{d^2x}{dt^2} = a$$ + +--- + +## Run the demos + +In order to run the demos, BioDynamo has to be correctly installed and sourced. + +```bash +cd [path_to_biodynamo]demos/SolarSystem +bdm run +bdm view +``` +In order to clearly observe the planets and the sun, feel free to scale up the diameter or volume of the satellites in the simulation. diff --git a/demo/solar_system/SolarSystemSim.png b/demo/solar_system/SolarSystemSim.png new file mode 100644 index 0000000000000000000000000000000000000000..18d60e134e68f512088ecff08f7544a9d812273a GIT binary patch literal 10146 zcmeHMXIPVImkwB9ETE$!C@oS%R1!dn)G(qTAfPBMG%*AW5K2M|UHwEyN)Ql4gp5=j zNDxAiuEPfbX#;^ox=aEXLXQCg1on+P*Us+FH`iYCV}Irb+U)2X$~P9$g(&%Yk6*zJ0*< z$s(TSNd^V5xsinBUYdE#D)ij^k~_y21lqS31z6>1byW`$fQEQ@2HZwM!qI_%xgd}M zEIiNy;fo9g-$r_){0-$7>zn1lC{IIq2Q4c#t3VUv9n_Tw4Dv>V^-V;CFGAN-9%i)9 zAY2bnfJO#;fWy&#{#dD5uUmqB4dorJu7OPgFi5Z_L=&Q>dN~{w zsv&Q*4{U((^wP7vWcr&5;G3cRo#5a=Jt#CREDREM9uk1@hN|o8>O$2tpc)#gzz$U` z&Og{AT-6_|AW-q|IxZox2n;GP7!}|T7U=Z29S{<1C@&AJgMU*lh`@hn_s9MQ0sse8 zumV+ws6qep_F$CPpKKSb{JvfAvfedRIMUDY5(fJu62L&5FKzxiyZ(mjFCp-kHvgSn|F_`U_dB^m`U8wE44`f! z?;HRs3fhgmXns=^_(X|%J_3QH%grx+e=~e|erycy$@+d`c>-%~*WB}EE<23AN-7US zGr0T3PW|}tuxW#$@qup-UH;zU-rGm=M~^BlU^CfQY9Hlo?J+;9B5^&n1SCAO_poT_ z8wIJks>YZR*TtTNuDPo=VRT%s%M0vJKSk2rq=tMy^P+w%sw*jpG_|u`r1y@fn8QJ! zapRm!5a)Ia#-~tGE_-Eq(t9!pi_Af{8w~+qN(TlS@0&ZThlpPq#4sR{uvstnAs=s) z61#$rJ<~b$^L}BFFf+2u9YXCdWN$<<6&AVoK=0P#c7BL`Pl%q)ikfFbmCd(5q}@8Z&p7 z`96m^EAJHCQIWXjn)*^ntHa64pP0I{GRT~O<+#AO8$w{?buk&$v{-Hv^L0swiGD{6 z;haUqTxt%4)wymupoo5)GI;c4>{1d_2R5m-`xB;oFEeuPgFJ6OOetVo%FHmswg2gI z^=U+u-$Wgc0U2Va$wh;CPc|EqfsAk3B^QZZS@OqJemtOq7T2wsrwmW~Dj^=7vsi&n zoVOkUPT^!oV+ahl_S7@%1ozHyI_^M7gL~ETc4@L)+3*E!#q%b&kqfsmI~u$s<%V34 zi)`roC~?C#=S^EDZzujh)r-GUS1`D&t7T3oNMGi>L=$&q$rmhMZ2`{dx^{QMRcA}< z>Gfywh|g{9j+Nc=10ivq2k;&nO(P}C+a_yDJ2S|#vmoJr1b74tPht>S$+v@@UtTYA zu2a{xIBtT}N{Sg^G3r~qoEPdz`wjs`p$t(h$C#`Z#0 z%}XPMP|PBcfgmFhBJB?i0A)ppke( z*{sEDzdyReDf7S$zpj~B#q7@?-F9K%6wz8kMfWC#)*j=yCB22K31$JrwqED9^^a@{ z>nYhk%FlcFfMJ=Cla~ z?zXa>&zx-Rj=*Ckbvfm8z2jEf>vih_3|xU^UHSO_?GdjsxkY@jNvA(SsW&HaCg$sP zu?{H*PJ4KZt_V{W8yuVx`=vy#TWw&Ava{8UVbjVb_QdJDIGhsMMcE#o=tAHdRRT8e zKI9}cFCp9xftsL+aia(^ad8Yrub?U;(>=6WAILnNcZ4k#6*t{UzfqTp-+8ax2nI#D zf84VR)G%2pBqX?p0(nlmgtw|V>kRR!92X~|0%KpxI=6~c%R@pLMZE>H;wp5r3=q8D1fH)bey8+WBhVP<*>#9 zEoEhI9(8B)!B)-Ab~AQ!)fIX{^w6hVJf1(c^)n?#b8Rdvh2NMnx3KWlrtPSxV1!?( znXyG?KQb{jof}Y#Q&c1tTIZiEE-q;Z=i6`3H*d<-2DQA80_ykvz3XR@v9ZdqO)~67 zOUqbLG$pV1_E6HYI*}(yk4{|fd#bM2#>$-E>`w36ux~5lMrXF>^^T8Q%f@_$Z?^_) zUnI%MgfD*|2U^Q(Y|PpEViTh#vWbW~x3NqM<4cOU>cpGzj@2ww!E%d{)r64{rG(3b z%|-n#*iI6m_lc;GknkvCxxBv1EOv_z1k`EKkF^PLg(Z8YbARox*v{41xt@yp8QB0G zih9Ahn;P;Dpqfu0%N?1+(GPi>&08J;y8%4h%IoX99>TV6uC8V z6&d~++Q-_YrO`5S^FLGOqh3|+@3rX9cD?&9aqTfp6sL5L`=r;gJc32l#xk|Jsowow7=_r3gPGV_Vs-q!>7mS zIgsXXSip+RK?VFMVR$&g(3JjE{Zbu7#R(tYq zMh`7joCqov&CS__jSle$?5D`~s#|f3B_$54R3BURg){04ZdrSkhokPHa*idPDhwTR z8!LCB1#Oq}GOBSs8)_BvA&SoV$|_3woko$7oei80oyf0t^Z}N|!|6}TPx#DUyYZ07 zpgzp7EgFfKGPSbp*YGcs5yARSe?O*f1W(n7UgZ6#D+kU-4EGXMP^vZt%dgmN02zRc zeWncFSHDoi=<_AAH*pro(^|(NQHZlBq9r}me)@V+MH5t_XY1oeThP8DH))~ciw{|g zIg1l3(;|CBfpAs~jf`Y%BVb8)D!ZTKU((X;a9kZCBEPUYXg702JcX!|$@4~Hlb*pT zho24lG|0#l3-mDls$mA+o_IW>fC(5)#DoCdp-pAbd-DUUFI`bALVDM}{8Twf362 zk(`mP*V_6eygwDf^06&GoRXe_IjQYaFA2_W_H|71wT-vy-X2fZ*|%FXFl%lRzZ2w6 zjOVTnb8#{u&+olfLhb{6K~?5kepqiHDWb(vJY}b^oMq50j)ztI?KO6T8T&}l^A$U; zIg5#}VSuw9APIn>p)$2KjJJ2TH?g9JQS%0Z49(&(F%nnVZO2qPyx?99P^8NXJ(AE+ z*LxR(c0c`#)n+c!Po8;bS+|HDIs|^F6tuR0BDon#ic4D-X@MC;tNcR@;eO&K_ljvI zZ8qXPOIi*c5!KW`U(zc&GjGRx1ba42mjrS3r0pARlW~Zl&F;y z&Dm=;zaq{&n$J1MA2>pM+6{T8$PY~qRw6#-e(BURfh5K?oG~80u%S_^U0iJR<(I0c zWy@Y%iy%htS*$aKYz^6C_B2X49ktIZ5w8pWzIrYP6qH4)t$`nQ&{eJ9m zG7+e#+OijENe<=n7Upb9LRa^7+gX|unDaEah8{P-)xyw=b9?H8dV1oe2`{GC*EPeE zZC188&`7et#8$O|0b0dp?U#XT3r@!pTD`8%hQ^hesnmUO?GHg8NFf#Yy27*N*gZ0r zsttoR+z*FEf9s&wWxUKj5BIw7NJ)#+Rjx`Ka7yG`X8q)LgeXH+fhX>>9~ zmVd{*g!7FThrxinsOJoH=A@148g+2lH4Wxq-G$mh_ETsjZn_J)+k#a8F+k?XZdI7r7Pw{SD$$b88*#&fxV$7A%1dX+aZyF z_l32s->I5vl_2f_4WtX)eWRa@>%9R>nc8B%@-sD=R((sXeom5dME(k+7ZU4WPXH$8plsRZ|} zh;c570^Iz%Izbv{?Kk%`)LAC+q^l=^7R2%x|Ly6JWV~%#OZ$fnN;=)(p5gJD6_ULc zCCnOIeWGan&FU`TWV6QVaLQ#381F_BV{~J8Ycw2 zNE|q9oX~L|J9WgvHS>BB*ONP}>d&QQbE6L>c8;j9qt7~Co7YOZP z`qZyD*)*E=^>W+-Z4xmP_KCzm@C#YB#l-}Z@k8_-+ zHNz+T`(t1D7ryO}P7@sxOX+@|Jiq^!?p~?=o^I^xG^!p$j=g@tsJj@(ty~9S9wi@M z*;d8`zNruSi52SjVq!>RBB1QMT1PPFn^%k*)B8u2(0beab&_E7qm&$^(UOUcS6G;L zV>U9niBhtWM>NxRqei5AMt;dJ#k;`(=C!l=RkhmX_WE(IV=~E=mvOja++KgjG}7+J z(zBk83;9+e_{Oc3ZUa8A><4k6m7Q^|4^_e&t}Y#wXOt1=MYGI|Y0U#o9jDEQ9^)mN zf)w0iS6fE`GAP(QVgE?lUq-E3*Ns04t>&eav+$&|INsC^r&|9+#3?r85*)^m*- z|FMql>@sFsXOSaN1VC5yph7@50TK)~{-yZCq(iXt_lE)<|8njD3^hI-e6#c|2s8rJ z`IU8q`1CSGw#!6Zc4~I;{(-9LDXAJPr)^%)08eE>n^T-?l6U%J&CL5720(!Uk?WJ2 zD;YWygte_t&YX6IHDvAQG{uC=nItwf#C#&e+`LX>ym{rrm;q-hfQjE>a(G+RUZ8*b zG{TuSTkI8y{XCD`9Fgn6)V5>#yahx!v8HLUf6cXy)S0(vYmqhkve0TzU#bEBfj39L zd95RA0XyOL2=L{th>Y~ybUtjm$zYWx0_1B2d5!`CBGoD_7Y%S=D$&xx!omSK&Z;Oz zG+H%wTf)6}XKT!ci6v~-0rMzrAk9}eur-C+(BPD@F=S|1rc;lQ28EKRud{L8qW`Hz zTl-OPVR|~to2GKf$M?9s7J9n>_&3i^ct}A&qYA65 zp({sV%XjDI=Sz!kZg2-a9F^o~TG9*CgSsgj)!)e%N(Hz&*w%F)waSXIo(RKHP0V=n z6C174GBN@vB<2?&oEq55&CoCjfcvMbYkr{JX}gejoi=ek7G4nz!$hc>Ct?fB3#(Y#u}$V_*j`o64%*Ow z(t+ul8H25781Vj?4fh>l z&{3;lZc$NFpXSbgDAzq$vMG4+J2SdiKYZlG`~}#@zKS0n*M;~J3u)Ec!#ro(Vpk|< zZOY}1--iOWg>3jIi|rv8e+X!?l#M~3D(V-T&H`l%?i~LTZK{Qd(L~6|RTT|p*y}7` zFFsTo9!Z?(R#2$wA>enCbEf1sH}!e#6A@>A@O=+@SGf9amxD&BPy9P&@NF`D^Zn!Jxn_fF|1H%Z_OZ+mfIAx6p`eAIRLhoRO2n02EQmOy4J! z$FXlaJv9DvcC}bab#6=+Eny&oln!KwXFz~i9OxZ#aLv`n%D%)Z>OhIdMtJ+}LF8bF z%9z*3eDvZ!SOWgS9^P0gr=^l?HtcQL)1S)K-99`3p|Yz^g`7Ww+2K0~oUu>1+$?ra z*T~+a+egb&a-%fE3{ex<_UKCh3@HoU5|V(V2~^G;{UmY`T| z#7ZPKc_`}oGj8%TI9&nE&`VQ~A-Rnm0WfizpF2lsB}8Wq`m|$e9VryG0j~T2z_|62 z;|c~?cPdmg;YD28v7>U-z$-reb(Q`{PO%gEs|u=uDiQiu;ytz^@yTL)1YT9qmUj=Z ze#HPy*UwZO6QR}R>_3Ccmf>M{Uas8@;&&f)T;k=b>T6eX2g`nCRVBWr zhOLFU>Xh9~Tq*tTRqrbN4GmSG$2HhVk`ibt8B|MjDC(Y`tt_$Y9vT9;;pjJH-L&Z8 z_Z7*H3aTKknUy4Bb8}%Kp~$VypS<6TkL@7}9S?i&2dgHNyXZ}!u0X^K52lcxlM?q| z*O7R678rj{F$?7G=b`c8FJQTi)fr8`uc;5%(=#;?R*`C5VUNRmRrd8X?N+m~I>)iC zj}FbNw*}MMaO*%ZI z%Eqo|o@YB%6(m@_d}?O^DPb@Vd?mpsCN+Srqio)V-Og zJWn@Gcq-5*xle^&C%=leH|fBXUvYQ#aiFx8x;hc=42PhHV>V~Gax=o7Kc&!8o7YF_ zEB%(HK>qb_esIX`eLp!;K_)*79nR<`W%`*y?Wu?7TYmADgahM|5Xs&{SwI=>_kktH z1T#X1xumiEzJ{~4txBq8qfOamF;hIQiNL0EGTAPd2M|Acy_`*nD(ifXbO!c&jm-9_ zsnNk-kbP^cjnvjWn~4xWzLaapJ3BQBc{U8U_@G&MK6ghwIX#i8Z3E-`k=&b=0Mp!p z*tDDm4BMI{D^o;!-o&n`+F7+b+^^yK&fks-Z9B_a9lNm4R)X>6!)7yVvV zF$QjF3$N>+S{S27&E9aouof`vqPrLLwlXyI7_K*RJQt|k$2x(*PfR1Z&_y%_O-wQ? z?WI3;S~y2jkWwnTMiWi8dHU_haT+|*mlQgW8@iuMM){exe;c7A}5FhFo!5) z(-b79vUj(y@>e+@9_#^gPNd*&dO8knP0Gcb0meRM8?^v13$I0;reCi$8P?b9J3%(n zajKNsH_tA5>d_I{Cdw&2nhh%Z+E8*jnM~dyz<6`MZ4W_0agmxAbHjgCtqq+R3<~QE zzBoI56;C7^R-0wf)gp7k`~m~x>a zy=V01k%T;hWV0*NE diameter_) { + SetPropagateStaticness(); + } + diameter_ = diameter; + UpdateVolume(); + } + + void SetVolume(real_t volume) { + volume_ = volume; + UpdateDiameter(); + } + + void SetMass(real_t mass) { SetDensity(mass / volume_); } + + void SetDensity(real_t density) { + if (density > density_) { + SetPropagateStaticness(); + } + density_ = density; + } + + void SetPosition(const Real3& position) override { + position_ = position; + SetPropagateStaticness(); + } + + void ChangeVolume(real_t speed) { + // scaling for integration step + auto* param = Simulation::GetActive()->GetParam(); + real_t delta = speed * param->simulation_time_step; + volume_ += delta; + if (volume_ < real_t(5.2359877E-7)) { + volume_ = real_t(5.2359877E-7); + } + UpdateDiameter(); + } + + void UpdateDiameter() { + // V = (4/3)*pi*r^3 = (pi/6)*diameter^3 + real_t diameter = std::cbrt(volume_ * 6 / Math::kPi); + if (diameter > diameter_) { + Base::SetPropagateStaticness(); + } + diameter_ = diameter; + } + + void UpdateVolume() { + // V = (4/3)*pi*r^3 = (pi/6)*diameter^3 + volume_ = Math::kPi / real_t(6) * std::pow(diameter_, 3); + } + + void UpdatePosition(const Real3& delta) { + position_ += delta; + SetPropagateStaticness(); + } + + // empty unused functions + void ApplyDisplacement(const Real3& displacement) override{}; + + Real3 CalculateDisplacement(const InteractionForce* force, + real_t squared_radius, real_t dt) override { + return Real3{0, 0, 0}; + } + + private: + Real3 position_ = {0, 0, 0}; + real_t diameter_ = 0; + real_t volume_ = 0; + real_t density_ = 0; + Real3 velocity_ = {0, 0, 0}; +}; + +// Star subclass +class Star : public CelestialObject { + BDM_AGENT_HEADER(Star, CelestialObject, 1); + + public: + Star() : CelestialObject() {} + + explicit Star(real_t diameter) : CelestialObject(diameter) {} + + explicit Star(const Real3 position) : CelestialObject(position) {} + + ~Star() override = default; +}; + +// Planet subclass +class Planet : public CelestialObject { + BDM_AGENT_HEADER(Planet, CelestialObject, 1); + + public: + Planet() : CelestialObject() {} + + explicit Planet(real_t diameter) : CelestialObject(diameter) {} + + explicit Planet(const Real3 position) : CelestialObject(position) {} + + explicit Planet(const Star* main_star, real_t diameter, real_t distance, + real_t initial_speed) { + Real3 pos = main_star->GetPosition(); + pos[0] += distance; + Real3 initial_velocity = {0, initial_speed, 0}; + this->SetVelocity(main_star->GetVelocity() + initial_velocity); + this->SetPosition(pos); + this->SetDiameter(diameter); + } + + ~Planet() override = default; +}; + +// Satellite subclass +class Satellite : public CelestialObject { + BDM_AGENT_HEADER(Satellite, CelestialObject, 1); + + public: + Satellite() : CelestialObject() {} + + explicit Satellite(real_t diameter) : CelestialObject(diameter) {} + + explicit Satellite(const Real3 position) : CelestialObject(position) {} + + explicit Satellite(const Planet* main_planet, real_t diameter, + real_t distance, real_t initial_speed) { + Real3 pos = main_planet->GetPosition(); + pos[0] += distance; + Real3 initial_velocity = {0, initial_speed, 0}; + this->SetVelocity(main_planet->GetVelocity() + initial_velocity); + + this->SetPosition(pos); + this->SetDiameter(diameter); + } + + ~Satellite() override = default; +}; + +} // namespace astrophysics +} // namespace bdm + +#endif // CELESTIALOBJECT_H_ diff --git a/demo/solar_system/src/DataReader.h b/demo/solar_system/src/DataReader.h new file mode 100644 index 000000000..d26be2c06 --- /dev/null +++ b/demo/solar_system/src/DataReader.h @@ -0,0 +1,70 @@ +// ----------------------------------------------------------------------------- +// +// Copyright (C) 2021 CERN & University of Surrey for the benefit of the +// BioDynaMo collaboration. All Rights Reserved. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// +// See the LICENSE file distributed with this work for details. +// See the NOTICE file distributed with this work for additional information +// regarding copyright ownership. +// +// ----------------------------------------------------------------------------- +#ifndef DATAREADER_H_ +#define DATAREADER_H_ + +#include + +namespace bdm { +namespace astrophysics { + +// Define the structure to hold each object's data +struct ObjectData { + std::string name; + real_t mass; + real_t diameter; + real_t orbital_velocity; + real_t distance_from_parent; + real_t angle; + + void ScaleDistance(real_t scale_factor) { + diameter *= scale_factor; + distance_from_parent *= scale_factor; + orbital_velocity *= scale_factor; + } +}; + +inline void LoadDataFromFile(const std::string& filename, + std::vector& data) { + std::ifstream file(filename); + std::string line; + std::getline(file, line); // Skip the header line + + while (std::getline(file, line)) { + std::istringstream iss(line); + ObjectData obj; + std::string name; + double mass; + double diameter; + double orbital_velocity; + double distance_from_parent; + double angle; + + if (iss >> name >> mass >> diameter >> orbital_velocity >> + distance_from_parent >> angle) { + obj.name = name; + obj.mass = mass; + obj.diameter = diameter; + obj.orbital_velocity = orbital_velocity; + obj.distance_from_parent = distance_from_parent; + obj.angle = angle; + data.push_back(obj); + } + } +} + +} // namespace astrophysics +} // namespace bdm + +#endif // NASADATA_H_ diff --git a/demo/solar_system/src/Gravity.h b/demo/solar_system/src/Gravity.h new file mode 100644 index 000000000..70117b706 --- /dev/null +++ b/demo/solar_system/src/Gravity.h @@ -0,0 +1,116 @@ +// ----------------------------------------------------------------------------- +// +// Copyright (C) 2021 CERN & University of Surrey for the benefit of the +// BioDynaMo collaboration. All Rights Reserved. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// +// See the LICENSE file distributed with this work for details. +// See the NOTICE file distributed with this work for additional information +// regarding copyright ownership. +// +// ----------------------------------------------------------------------------- +#ifndef GRAVITY_H_ +#define GRAVITY_H_ + +#include "CelestialObjects.h" +#include "biodynamo.h" + +namespace bdm { +namespace astrophysics { + +// Gravitational constant in m^3/kg/s^2 +inline real_t gG = 6.674e-11; + +// set units of G constant +inline void SetUnits(char units = 'O', int times = 1) { + gG = 6.674e-11; + switch (units) { + case 'k': + gG *= pow(times * 1e-3, 3); + break; + case 'M': + gG *= pow(times * 1e-6, 3); + break; + case 'G': + gG *= pow(times * 1e-9, 3); + break; + default: + gG *= pow(times, 3); + break; + } +} + +// a behavior that applies Newtons Law for gravity and calculates the +// displacement for any celestial object +class Gravity : public Behavior { + Real3 acceleration_ = {0, 0, 0}; + real_t dt_ = Simulation::GetActive()->GetParam()->simulation_time_step; + + public: + Behavior* New() const override { return new Gravity(); } + + Behavior* NewCopy() const override { return new Gravity(*this); } + + void Run(Agent* agent) override { + CelestialObject* celestial = dynamic_cast(agent); + + Real3 celestial_pos = celestial->GetPosition(); + + auto* functor = new GravityFunctor(celestial, gG); + // for each agent calculate acceleration + Simulation::GetActive()->GetResourceManager()->ForEachAgentParallel( + *functor); + + dt_ = Simulation::GetActive()->GetParam()->simulation_time_step; + Real3 velocity = celestial->GetVelocity(); + Real3 acceleration = functor->acceleration; + Real3 position; + + // explicit Euler's method + velocity += acceleration * dt_; + position = celestial_pos + velocity * dt_; + + delete functor; + + celestial->SetVelocity(velocity); + celestial->SetPosition(position); + } + + class GravityFunctor : public Functor { + public: + CelestialObject* celestial; + real_t g; + Real3 acceleration; + + GravityFunctor(CelestialObject* celestial, real_t g) + : celestial(celestial), g(g), acceleration({0, 0, 0}) {} + + // override the function that overloads () operator to + // calculate acceleration based on given agent + void operator()(Agent* other_agent, AgentHandle) override { + CelestialObject* other_celestial = + dynamic_cast(other_agent); + + if (!other_celestial || other_celestial == celestial) { + return; + } + + Real3 dir = other_celestial->GetPosition() - celestial->GetPosition(); + real_t magnitude = dir.Norm(); + + dir.Normalize(); + real_t m1 = celestial->GetMass(); + real_t m2 = other_celestial->GetMass(); + real_t f = g * m1 * m2 / (magnitude * magnitude); + + acceleration += (f / m1) * dir; + }; + }; +}; +} // namespace astrophysics + +} // namespace bdm + +#endif // GRAVITY_H_ diff --git a/demo/solar_system/src/solar_system.cc b/demo/solar_system/src/solar_system.cc new file mode 100644 index 000000000..c68a45a39 --- /dev/null +++ b/demo/solar_system/src/solar_system.cc @@ -0,0 +1,51 @@ +// ----------------------------------------------------------------------------- +// +// Copyright (C) 2021 CERN & University of Surrey for the benefit of the +// BioDynaMo collaboration. All Rights Reserved. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// +// See the LICENSE file distributed with this work for details. +// See the NOTICE file distributed with this work for additional information +// regarding copyright ownership. +// +// ----------------------------------------------------------------------------- +#include "solar_system.h" + +int main(int argc, const char** argv) { + return bdm::astrophysics::Simulate(argc, argv); +} + +namespace bdm { +namespace astrophysics { + +Real3 RotateVector(const Real3& vec, real_t th, char axis) { + th *= M_PI / 180; + real_t costh = std::cos(th), sinth = std::sin(th); + + Real3 rot_vec; + switch (axis) { + case 'x': + rot_vec[0] = vec[0]; + rot_vec[1] = vec[1] * costh - vec[2] * sinth; + rot_vec[2] = vec[1] * sinth + vec[2] * costh; + break; + case 'y': + rot_vec[0] = vec[0] * costh + vec[2] * sinth; + rot_vec[1] = vec[1]; + rot_vec[2] = -vec[0] * sinth + vec[2] * costh; + break; + case 'z': + rot_vec[0] = vec[0] * costh - vec[1] * sinth; + rot_vec[1] = vec[0] * sinth + vec[1] * costh; + rot_vec[2] = vec[2]; + break; + default: + throw std::invalid_argument( + "Invalid rotation axis. Must be 'x', 'y', or 'z'."); + } + return rot_vec; +} +} // namespace astrophysics +} // namespace bdm diff --git a/demo/solar_system/src/solar_system.h b/demo/solar_system/src/solar_system.h new file mode 100644 index 000000000..f1ecf72cc --- /dev/null +++ b/demo/solar_system/src/solar_system.h @@ -0,0 +1,63 @@ +// ----------------------------------------------------------------------------- +// +// Copyright (C) 2021 CERN & University of Surrey for the benefit of the +// BioDynaMo collaboration. All Rights Reserved. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// +// See the LICENSE file distributed with this work for details. +// See the NOTICE file distributed with this work for additional information +// regarding copyright ownership. +// +// ----------------------------------------------------------------------------- +#ifndef SOLARSYSTEM_H_ +#define SOLARSYSTEM_H_ + +#include "CelestialObjects.h" +#include "DataReader.h" +#include "Gravity.h" +#include "biodynamo.h" + +namespace bdm { +namespace astrophysics { + +Real3 RotateVector(const Real3& vec, real_t th, char axis); + +inline int Simulate(int argc, const char** argv) { + Simulation simulation(argc, argv); + auto* rm = simulation.GetResourceManager(); + + std::vector data; + std::string filename = "planets.txt"; + + LoadDataFromFile(filename, data); + + constexpr real_t kDiameterScalingFactor = 10; + SetUnits('M'); + + // create sun + Star* sun = new Star(1392.7 * kDiameterScalingFactor); + sun->SetPosition({0, 0, 0}); + sun->SetMass(1988400e+24); + sun->AddBehavior(new Gravity); + rm->AddAgent(sun); + + // read data from .txt file and create the planets + for (auto object : data) { + auto* planet = new Planet(object.diameter * kDiameterScalingFactor); + planet->SetMass(object.mass * 1e+24); + planet->SetPosition({object.distance_from_parent, 0, 0}); + Real3 v = RotateVector({0, object.orbital_velocity, 0}, object.angle, 'x'); + planet->SetVelocity(v); + planet->AddBehavior(new Gravity); + rm->AddAgent(planet); + } + + simulation.GetScheduler()->Simulate(900); + return 0; +} + +} // namespace astrophysics +} // namespace bdm +#endif // SolarSystem.h From 72e12ab7d65cd1fb9ed8949a93fdc95d5249f199 Mon Sep 17 00:00:00 2001 From: christosneg Date: Mon, 10 Feb 2025 11:27:59 +0200 Subject: [PATCH 4/8] Created an algorithm in python that verifies results with real data for mars and solar_system demos --- demo/mars/check_results.py | 58 ++++++++++++++++++++++++++ demo/solar_system/check_results.py | 65 ++++++++++++++++++++++++++++++ 2 files changed, 123 insertions(+) create mode 100644 demo/mars/check_results.py create mode 100644 demo/solar_system/check_results.py diff --git a/demo/mars/check_results.py b/demo/mars/check_results.py new file mode 100644 index 000000000..3b693d086 --- /dev/null +++ b/demo/mars/check_results.py @@ -0,0 +1,58 @@ +import vtk +import os + +steps = 500 +ptuFiles = 8 +reader = vtk.vtkXMLUnstructuredGridReader() + +# first read all locations + +satellites = ["Phobos", "Deimos"] +locations = [[], []] +try: + # Read pvtu file for each time step and extract files + for i in range(steps): + k = 0 + for j in range(ptuFiles): + + # Find the output files + path = f"output/mars/Satellite-{i}_{j}.vtu" + if not os.path.exists(path): + raise Exception("Output file not found") + reader.SetFileName(path) + reader.Update() + output = reader.GetOutput() + + # Get all points from each file + for l in range(output.GetNumberOfPoints()): + locations[k].append(output.GetPoints().GetPoint(l)) + k = k + 1 + + # Scientific Data + # Since each step is 0.1 days + data = [76.6, 303.5] + stars = "**********************" + print(f'\n\n{stars}{stars}{stars}\n') + for satellite in satellites: + index = satellites.index(satellite) + # Find location after 1 full rotation + for step in range(1,steps): + + A = locations[index][step][1] > 0 + B = locations[index][step-1][1] < 0 + + if(A and B): + # Find the relative error + # Relative error = |x - x_0| / x + # where x is actual value + # and x is the measured value + relativeError = abs(step - data[index]) / data[index] + print(f" The relative error of {satellite}'s full rotation is {relativeError*100:.2f}%") + break + print(f'\n{stars}{stars}{stars}') + print(f"\n\nPrevious runs of this algorithm returned the following output:\n") + print(f" The relative error of Phobo's full rotation is 0.52%") + print(f" The relative error of Deimo's full rotation is 0.16%") + +except Exception as e: + print(f"error: {e}") diff --git a/demo/solar_system/check_results.py b/demo/solar_system/check_results.py new file mode 100644 index 000000000..443b44b30 --- /dev/null +++ b/demo/solar_system/check_results.py @@ -0,0 +1,65 @@ +import vtk +import os + +steps = 900 +ptuFiles = 8 +reader = vtk.vtkXMLUnstructuredGridReader() + +# first read all locations + +planets = ["Mercury", "Venus", "Earth", + "Mars", "Jupiter", "Saturn", "Uranus", + "Neptune", "Pluto"] +locations = [[], [], [], [], [], [], [], [], []] +try: + # Read pvtu file for each time step and extract files + for i in range(steps): + k = 0 + for j in range(ptuFiles): + + # Find the output files + path = f"output/solar_system/Planet-{i}_{j}.vtu" + if not os.path.exists(path): + raise Exception("Output file not found") + reader.SetFileName(path) + reader.Update() + output = reader.GetOutput() + + # Get all points from each file + for l in range(output.GetNumberOfPoints()): + locations[k].append(output.GetPoints().GetPoint(l)) + k = k + 1 + + # Scientific Data + # How many days/steps will it take for a full rotation + data = [88.0, 224.7, 365.24, 687.0, + 331.0, 10.747, 30.589, + 59800.0, 90560.0] + stars = "**********************" + print(f'\n\n{stars}{stars}{stars}\n') + for planet in planets: + index = planets.index(planet) + # Find location after 1 full rotation + for step in range(1,steps): + + A = locations[index][step][1] > 0 + B = locations[index][step-1][1] < 0 + + if(A and B): + # Find the relative error + # Relative error = |x - x_0| / x + # where x is actual value + # and x is the measured value + relativeError = abs(step - data[index]) / data[index] + print(f" The relative error of {planet}'s full rotation is {relativeError*100:.2f}%") + break + print(f'\n{stars}{stars}{stars}') + print(f"\n\nPrevious runs of this algorithm returned the following output:\n") + print(f" The relative error of Mercury's full rotation is 3.41%") + print(f" The relative error of Venus's full rotation is 0.31%") + print(f" The relative error of Earth's full rotation is 0.07%") + print(f" The relative error of Mars's full rotation is 0.44%") + print(f"\n\n ** The rest of the planets did not complete a full roation") + +except Exception as e: + print(f"error: {e}") \ No newline at end of file From 08653937072557377d5454bbf8dc1f31b68ae44a Mon Sep 17 00:00:00 2001 From: christosneg Date: Mon, 10 Feb 2025 11:52:28 +0200 Subject: [PATCH 5/8] Updated README.md files for mars and solar_system demos. --- demo/mars/README.md | 16 ++++++++++++++++ demo/solar_system/README.md | 17 +++++++++++++++++ 2 files changed, 33 insertions(+) diff --git a/demo/mars/README.md b/demo/mars/README.md index 4dc6f7781..a4a0ea1e0 100644 --- a/demo/mars/README.md +++ b/demo/mars/README.md @@ -55,3 +55,19 @@ bdm view ``` In order to clearly observe the satellites in relation to the planet, feel free to scale up the diameter or volume of the satellites in the simulation. +## Verify the Results + +By running the following commands: + +```bash +pip install vtk +``` + +```bash +cd [path_to_biodynamo]/demos/SolarSystem +python3 check_results.py +``` + +You can verify the accuracy of the simulation by comparing the time it takes for a celestial object to complete a full rotation. + +The data were acquired from the [NASA Mars Moons Facts](https://science.nasa.gov/mars/moons/facts/), which provides real-world orbital and rotational period measurements. diff --git a/demo/solar_system/README.md b/demo/solar_system/README.md index a97578d50..1ba788653 100644 --- a/demo/solar_system/README.md +++ b/demo/solar_system/README.md @@ -63,3 +63,20 @@ bdm run bdm view ``` In order to clearly observe the planets and the sun, feel free to scale up the diameter or volume of the satellites in the simulation. + +## Verify the Results + +By running the following commands: + +```bash +pip install vtk +``` + +```bash +cd [path_to_biodynamo]/demos/SolarSystem +python3 check_results.py +``` + +You can verify the accuracy of the simulation by comparing the time it takes for a celestial object to complete a full rotation. + +The data were acquired from the [NASA Planetary Fact Sheet](https://nssdc.gsfc.nasa.gov/planetary/factsheet/), which provides real-world orbital and rotational period measurements. From 7be10d4d9ae98914a13948a2cb861fa124b8aee2 Mon Sep 17 00:00:00 2001 From: christosneg Date: Wed, 5 Mar 2025 23:49:05 +0200 Subject: [PATCH 6/8] Fixed the scaling problem with solar_system demo by making the following changes 1) converted the units from Mega meters to Giga meters 2) The diameter of each agent was set relatively to earth instead of their actual size. This change will only affect the visualization of the simulation and not its accuracy. --- demo/solar_system/planets.txt | 20 ++++++++++---------- demo/solar_system/src/solar_system.h | 7 +++---- 2 files changed, 13 insertions(+), 14 deletions(-) diff --git a/demo/solar_system/planets.txt b/demo/solar_system/planets.txt index 91aa3b614..37787e8e4 100644 --- a/demo/solar_system/planets.txt +++ b/demo/solar_system/planets.txt @@ -1,10 +1,10 @@ -Planet Mass(10^24kg) Diameter(Mm) OrbitalVelocity(Mm/s) Distance_from_parent(Mm) Angle(deg) -MERCURY 0.330 4.879 0.0474 57900 7.0 -VENUS 4.87 12.104 0.0350 108200 3.4 -EARTH 5.97 12.756 0.0298 149600 0.0 -MARS 0.642 6.792 0.0241 227900 1.8 -JUPITER 1898 142.984 0.0131 778500 1.3 -SATURN 568 120.536 0.0097 1432000 2.5 -URANUS 86.8 51.118 0.0068 2871000 1.8 -NEPTUNE 102 49.528 0.0054 4495000 0.8 -PLUTO 0.0130 2.376 0.0047 5906000 17.2 +Planet Mass(10^24kg) Diameter(relative_to_earth) OrbitalVelocity(Tm/s) Distance_from_parent(Tm) Angle(deg) +MERCURY 0.330 0.383 0.0000474 57.90 7.0 +VENUS 4.87 0.949 0.0000350 108.2 3.4 +EARTH 5.97 1 0.0000298 149.6 0.0 +MARS 0.642 0.2724 0.0000241 227.9 1.8 +JUPITER 1898 11.21 0.0000131 778.5 1.3 +SATURN 568 9.45 0.0000097 1432.0 2.5 +URANUS 86.8 4.01 0.0000068 2871.0 1.8 +NEPTUNE 102 3.88 0.0000054 4495.0 0.8 +PLUTO 0.0130 0.187 0.0000047 5906.0 17.2 diff --git a/demo/solar_system/src/solar_system.h b/demo/solar_system/src/solar_system.h index f1ecf72cc..575d5330a 100644 --- a/demo/solar_system/src/solar_system.h +++ b/demo/solar_system/src/solar_system.h @@ -33,11 +33,10 @@ inline int Simulate(int argc, const char** argv) { LoadDataFromFile(filename, data); - constexpr real_t kDiameterScalingFactor = 10; - SetUnits('M'); + SetUnits('G'); // create sun - Star* sun = new Star(1392.7 * kDiameterScalingFactor); + Star* sun = new Star(109); sun->SetPosition({0, 0, 0}); sun->SetMass(1988400e+24); sun->AddBehavior(new Gravity); @@ -45,7 +44,7 @@ inline int Simulate(int argc, const char** argv) { // read data from .txt file and create the planets for (auto object : data) { - auto* planet = new Planet(object.diameter * kDiameterScalingFactor); + auto* planet = new Planet(object.diameter); planet->SetMass(object.mass * 1e+24); planet->SetPosition({object.distance_from_parent, 0, 0}); Real3 v = RotateVector({0, object.orbital_velocity, 0}, object.angle, 'x'); From fdb628808a3098da4baf3484e041206c97a1dfbb Mon Sep 17 00:00:00 2001 From: christosneg Date: Thu, 6 Mar 2025 00:05:21 +0200 Subject: [PATCH 7/8] Removed unused subclasses from CelestialObject.h. 1) removed Star subclass from mars demo 2) removed Satellite subclass from solar_system demo 3) removed all subclasses from newtons_law_test demo --- demo/mars/src/CelestialObjects.h | 24 -------- demo/newtons_law_test/src/CelestialObjects.h | 63 -------------------- demo/solar_system/src/CelestialObjects.h | 25 -------- 3 files changed, 112 deletions(-) diff --git a/demo/mars/src/CelestialObjects.h b/demo/mars/src/CelestialObjects.h index 10cddf869..f6c7f09a1 100644 --- a/demo/mars/src/CelestialObjects.h +++ b/demo/mars/src/CelestialObjects.h @@ -130,20 +130,6 @@ class CelestialObject : public Agent { Real3 velocity_ = {0, 0, 0}; }; -// Star subclass -class Star : public CelestialObject { - BDM_AGENT_HEADER(Star, CelestialObject, 1); - - public: - Star() : CelestialObject() {} - - explicit Star(real_t diameter) : CelestialObject(diameter) {} - - explicit Star(const Real3 position) : CelestialObject(position) {} - - ~Star() override = default; -}; - // Planet subclass class Planet : public CelestialObject { BDM_AGENT_HEADER(Planet, CelestialObject, 1); @@ -155,16 +141,6 @@ class Planet : public CelestialObject { explicit Planet(const Real3 position) : CelestialObject(position) {} - explicit Planet(const Star* main_star, real_t diameter, real_t distance, - real_t initial_speed) { - Real3 pos = main_star->GetPosition(); - pos[0] += distance; - Real3 initial_velocity = {0, initial_speed, 0}; - this->SetVelocity(main_star->GetVelocity() + initial_velocity); - this->SetPosition(pos); - this->SetDiameter(diameter); - } - ~Planet() override = default; }; diff --git a/demo/newtons_law_test/src/CelestialObjects.h b/demo/newtons_law_test/src/CelestialObjects.h index 10cddf869..0d275e2c1 100644 --- a/demo/newtons_law_test/src/CelestialObjects.h +++ b/demo/newtons_law_test/src/CelestialObjects.h @@ -130,69 +130,6 @@ class CelestialObject : public Agent { Real3 velocity_ = {0, 0, 0}; }; -// Star subclass -class Star : public CelestialObject { - BDM_AGENT_HEADER(Star, CelestialObject, 1); - - public: - Star() : CelestialObject() {} - - explicit Star(real_t diameter) : CelestialObject(diameter) {} - - explicit Star(const Real3 position) : CelestialObject(position) {} - - ~Star() override = default; -}; - -// Planet subclass -class Planet : public CelestialObject { - BDM_AGENT_HEADER(Planet, CelestialObject, 1); - - public: - Planet() : CelestialObject() {} - - explicit Planet(real_t diameter) : CelestialObject(diameter) {} - - explicit Planet(const Real3 position) : CelestialObject(position) {} - - explicit Planet(const Star* main_star, real_t diameter, real_t distance, - real_t initial_speed) { - Real3 pos = main_star->GetPosition(); - pos[0] += distance; - Real3 initial_velocity = {0, initial_speed, 0}; - this->SetVelocity(main_star->GetVelocity() + initial_velocity); - this->SetPosition(pos); - this->SetDiameter(diameter); - } - - ~Planet() override = default; -}; - -// Satellite subclass -class Satellite : public CelestialObject { - BDM_AGENT_HEADER(Satellite, CelestialObject, 1); - - public: - Satellite() : CelestialObject() {} - - explicit Satellite(real_t diameter) : CelestialObject(diameter) {} - - explicit Satellite(const Real3 position) : CelestialObject(position) {} - - explicit Satellite(const Planet* main_planet, real_t diameter, - real_t distance, real_t initial_speed) { - Real3 pos = main_planet->GetPosition(); - pos[0] += distance; - Real3 initial_velocity = {0, initial_speed, 0}; - this->SetVelocity(main_planet->GetVelocity() + initial_velocity); - - this->SetPosition(pos); - this->SetDiameter(diameter); - } - - ~Satellite() override = default; -}; - } // namespace astrophysics } // namespace bdm diff --git a/demo/solar_system/src/CelestialObjects.h b/demo/solar_system/src/CelestialObjects.h index 10cddf869..6348344d0 100644 --- a/demo/solar_system/src/CelestialObjects.h +++ b/demo/solar_system/src/CelestialObjects.h @@ -168,31 +168,6 @@ class Planet : public CelestialObject { ~Planet() override = default; }; -// Satellite subclass -class Satellite : public CelestialObject { - BDM_AGENT_HEADER(Satellite, CelestialObject, 1); - - public: - Satellite() : CelestialObject() {} - - explicit Satellite(real_t diameter) : CelestialObject(diameter) {} - - explicit Satellite(const Real3 position) : CelestialObject(position) {} - - explicit Satellite(const Planet* main_planet, real_t diameter, - real_t distance, real_t initial_speed) { - Real3 pos = main_planet->GetPosition(); - pos[0] += distance; - Real3 initial_velocity = {0, initial_speed, 0}; - this->SetVelocity(main_planet->GetVelocity() + initial_velocity); - - this->SetPosition(pos); - this->SetDiameter(diameter); - } - - ~Satellite() override = default; -}; - } // namespace astrophysics } // namespace bdm From 657e2c9aead150b568e695d73cce87791b739cb9 Mon Sep 17 00:00:00 2001 From: christosneg Date: Thu, 6 Mar 2025 00:32:45 +0200 Subject: [PATCH 8/8] The CelestialObject class in all demos now inherits from Cell class instead of Agent class. --- demo/mars/src/CelestialObjects.h | 95 ++------------------ demo/newtons_law_test/src/CelestialObjects.h | 95 ++------------------ demo/solar_system/src/CelestialObjects.h | 95 ++------------------ 3 files changed, 21 insertions(+), 264 deletions(-) diff --git a/demo/mars/src/CelestialObjects.h b/demo/mars/src/CelestialObjects.h index f6c7f09a1..1c8b951d4 100644 --- a/demo/mars/src/CelestialObjects.h +++ b/demo/mars/src/CelestialObjects.h @@ -21,20 +21,20 @@ namespace astrophysics { // creating the custom agent that // defines celestial objects -class CelestialObject : public Agent { - BDM_AGENT_HEADER(CelestialObject, Agent, 1); +class CelestialObject : public Cell { + BDM_AGENT_HEADER(CelestialObject, Cell, 1); public: // constructors - CelestialObject() : diameter_(1.0), density_(1.0) { UpdateVolume(); } + CelestialObject() { UpdateVolume(); } - explicit CelestialObject(real_t diameter) - : diameter_(diameter), density_(1.0) { + explicit CelestialObject(real_t diameter){ + this->SetDiameter(diameter); UpdateVolume(); } - explicit CelestialObject(const Real3& position) - : position_(position), diameter_(1.0), density_(1.0) { + explicit CelestialObject(const Real3& position){ + this->SetPosition(position); UpdateVolume(); } @@ -43,90 +43,9 @@ class CelestialObject : public Agent { Real3 GetVelocity() const { return velocity_; } - real_t GetDiameter() const override { return diameter_; } - - real_t GetMass() const { return density_ * volume_; } - - real_t GetDensity() const { return density_; } - - const Real3& GetPosition() const override { return position_; } - - real_t GetVolume() const { return volume_; } - - Shape GetShape() const override { return Shape::kSphere; } - void SetVelocity(const Real3& vel) { velocity_ = vel; } - void SetDiameter(real_t diameter) override { - if (diameter > diameter_) { - SetPropagateStaticness(); - } - diameter_ = diameter; - UpdateVolume(); - } - - void SetVolume(real_t volume) { - volume_ = volume; - UpdateDiameter(); - } - - void SetMass(real_t mass) { SetDensity(mass / volume_); } - - void SetDensity(real_t density) { - if (density > density_) { - SetPropagateStaticness(); - } - density_ = density; - } - - void SetPosition(const Real3& position) override { - position_ = position; - SetPropagateStaticness(); - } - - void ChangeVolume(real_t speed) { - // scaling for integration step - auto* param = Simulation::GetActive()->GetParam(); - real_t delta = speed * param->simulation_time_step; - volume_ += delta; - if (volume_ < real_t(5.2359877E-7)) { - volume_ = real_t(5.2359877E-7); - } - UpdateDiameter(); - } - - void UpdateDiameter() { - // V = (4/3)*pi*r^3 = (pi/6)*diameter^3 - real_t diameter = std::cbrt(volume_ * 6 / Math::kPi); - if (diameter > diameter_) { - Base::SetPropagateStaticness(); - } - diameter_ = diameter; - } - - void UpdateVolume() { - // V = (4/3)*pi*r^3 = (pi/6)*diameter^3 - volume_ = Math::kPi / real_t(6) * std::pow(diameter_, 3); - } - - void UpdatePosition(const Real3& delta) { - position_ += delta; - SetPropagateStaticness(); - } - - // empty unused functions - void ApplyDisplacement(const Real3& displacement) override{}; - - Real3 CalculateDisplacement(const InteractionForce* force, - real_t squared_radius, real_t dt) override { - return Real3{0, 0, 0}; - } - private: - Real3 position_ = {0, 0, 0}; - real_t diameter_ = 0; - real_t volume_ = 0; - real_t density_ = 0; Real3 velocity_ = {0, 0, 0}; }; diff --git a/demo/newtons_law_test/src/CelestialObjects.h b/demo/newtons_law_test/src/CelestialObjects.h index 0d275e2c1..d18727211 100644 --- a/demo/newtons_law_test/src/CelestialObjects.h +++ b/demo/newtons_law_test/src/CelestialObjects.h @@ -21,20 +21,20 @@ namespace astrophysics { // creating the custom agent that // defines celestial objects -class CelestialObject : public Agent { - BDM_AGENT_HEADER(CelestialObject, Agent, 1); +class CelestialObject : public Cell { + BDM_AGENT_HEADER(CelestialObject, Cell, 1); public: // constructors - CelestialObject() : diameter_(1.0), density_(1.0) { UpdateVolume(); } + CelestialObject() { UpdateVolume(); } - explicit CelestialObject(real_t diameter) - : diameter_(diameter), density_(1.0) { + explicit CelestialObject(real_t diameter){ + this->SetDiameter(diameter); UpdateVolume(); } - explicit CelestialObject(const Real3& position) - : position_(position), diameter_(1.0), density_(1.0) { + explicit CelestialObject(const Real3& position){ + this->SetPosition(position); UpdateVolume(); } @@ -43,90 +43,9 @@ class CelestialObject : public Agent { Real3 GetVelocity() const { return velocity_; } - real_t GetDiameter() const override { return diameter_; } - - real_t GetMass() const { return density_ * volume_; } - - real_t GetDensity() const { return density_; } - - const Real3& GetPosition() const override { return position_; } - - real_t GetVolume() const { return volume_; } - - Shape GetShape() const override { return Shape::kSphere; } - void SetVelocity(const Real3& vel) { velocity_ = vel; } - void SetDiameter(real_t diameter) override { - if (diameter > diameter_) { - SetPropagateStaticness(); - } - diameter_ = diameter; - UpdateVolume(); - } - - void SetVolume(real_t volume) { - volume_ = volume; - UpdateDiameter(); - } - - void SetMass(real_t mass) { SetDensity(mass / volume_); } - - void SetDensity(real_t density) { - if (density > density_) { - SetPropagateStaticness(); - } - density_ = density; - } - - void SetPosition(const Real3& position) override { - position_ = position; - SetPropagateStaticness(); - } - - void ChangeVolume(real_t speed) { - // scaling for integration step - auto* param = Simulation::GetActive()->GetParam(); - real_t delta = speed * param->simulation_time_step; - volume_ += delta; - if (volume_ < real_t(5.2359877E-7)) { - volume_ = real_t(5.2359877E-7); - } - UpdateDiameter(); - } - - void UpdateDiameter() { - // V = (4/3)*pi*r^3 = (pi/6)*diameter^3 - real_t diameter = std::cbrt(volume_ * 6 / Math::kPi); - if (diameter > diameter_) { - Base::SetPropagateStaticness(); - } - diameter_ = diameter; - } - - void UpdateVolume() { - // V = (4/3)*pi*r^3 = (pi/6)*diameter^3 - volume_ = Math::kPi / real_t(6) * std::pow(diameter_, 3); - } - - void UpdatePosition(const Real3& delta) { - position_ += delta; - SetPropagateStaticness(); - } - - // empty unused functions - void ApplyDisplacement(const Real3& displacement) override{}; - - Real3 CalculateDisplacement(const InteractionForce* force, - real_t squared_radius, real_t dt) override { - return Real3{0, 0, 0}; - } - private: - Real3 position_ = {0, 0, 0}; - real_t diameter_ = 0; - real_t volume_ = 0; - real_t density_ = 0; Real3 velocity_ = {0, 0, 0}; }; diff --git a/demo/solar_system/src/CelestialObjects.h b/demo/solar_system/src/CelestialObjects.h index 6348344d0..9e65161f0 100644 --- a/demo/solar_system/src/CelestialObjects.h +++ b/demo/solar_system/src/CelestialObjects.h @@ -21,20 +21,20 @@ namespace astrophysics { // creating the custom agent that // defines celestial objects -class CelestialObject : public Agent { - BDM_AGENT_HEADER(CelestialObject, Agent, 1); +class CelestialObject : public Cell { + BDM_AGENT_HEADER(CelestialObject, Cell, 1); public: // constructors - CelestialObject() : diameter_(1.0), density_(1.0) { UpdateVolume(); } + CelestialObject() { UpdateVolume(); } - explicit CelestialObject(real_t diameter) - : diameter_(diameter), density_(1.0) { + explicit CelestialObject(real_t diameter){ + this->SetDiameter(diameter); UpdateVolume(); } - explicit CelestialObject(const Real3& position) - : position_(position), diameter_(1.0), density_(1.0) { + explicit CelestialObject(const Real3& position){ + this->SetPosition(position); UpdateVolume(); } @@ -43,90 +43,9 @@ class CelestialObject : public Agent { Real3 GetVelocity() const { return velocity_; } - real_t GetDiameter() const override { return diameter_; } - - real_t GetMass() const { return density_ * volume_; } - - real_t GetDensity() const { return density_; } - - const Real3& GetPosition() const override { return position_; } - - real_t GetVolume() const { return volume_; } - - Shape GetShape() const override { return Shape::kSphere; } - void SetVelocity(const Real3& vel) { velocity_ = vel; } - void SetDiameter(real_t diameter) override { - if (diameter > diameter_) { - SetPropagateStaticness(); - } - diameter_ = diameter; - UpdateVolume(); - } - - void SetVolume(real_t volume) { - volume_ = volume; - UpdateDiameter(); - } - - void SetMass(real_t mass) { SetDensity(mass / volume_); } - - void SetDensity(real_t density) { - if (density > density_) { - SetPropagateStaticness(); - } - density_ = density; - } - - void SetPosition(const Real3& position) override { - position_ = position; - SetPropagateStaticness(); - } - - void ChangeVolume(real_t speed) { - // scaling for integration step - auto* param = Simulation::GetActive()->GetParam(); - real_t delta = speed * param->simulation_time_step; - volume_ += delta; - if (volume_ < real_t(5.2359877E-7)) { - volume_ = real_t(5.2359877E-7); - } - UpdateDiameter(); - } - - void UpdateDiameter() { - // V = (4/3)*pi*r^3 = (pi/6)*diameter^3 - real_t diameter = std::cbrt(volume_ * 6 / Math::kPi); - if (diameter > diameter_) { - Base::SetPropagateStaticness(); - } - diameter_ = diameter; - } - - void UpdateVolume() { - // V = (4/3)*pi*r^3 = (pi/6)*diameter^3 - volume_ = Math::kPi / real_t(6) * std::pow(diameter_, 3); - } - - void UpdatePosition(const Real3& delta) { - position_ += delta; - SetPropagateStaticness(); - } - - // empty unused functions - void ApplyDisplacement(const Real3& displacement) override{}; - - Real3 CalculateDisplacement(const InteractionForce* force, - real_t squared_radius, real_t dt) override { - return Real3{0, 0, 0}; - } - private: - Real3 position_ = {0, 0, 0}; - real_t diameter_ = 0; - real_t volume_ = 0; - real_t density_ = 0; Real3 velocity_ = {0, 0, 0}; };